Development of Artificial Neural Network Based Metamodels for Inactivation of Anthrax Spores in Ventilated Spaces Using Computational Fluid Dynamics

ABSTRACT

Linear, quadratic, and artificial neural network (ANN)-based metamodels were developed for predicting the extent of anthrax spore inactivation by chlorine dioxide in a ventilated three-dimensional space over time from computational fluid dynamics model (CFD) simulation data. Dimensionless groups were developed to define the design space of the problem scenario. The Hammersley sequence sampling (HSS) method was used to determine the sampling points for the numerical experiments within the design space. A CFD model, comprised of multiple submodels, was applied to conduct the numerical experiments. Large eddy simulation (LES) with the Smagorinsky subgrid-scale model was applied to compute the airflow. Anthrax spores were modeled as a dispersed solid phase using the Lagrangian treatment. The disinfectant transport was calculated by solving a mass transport equation. Kinetic decay constants were included for spontaneous decay of the disinfectant and for the reaction of the disinfectant with the surfaces of the three-dimensional space. To enhance the mixing of the disinfectant with the room air, a momentum source was included in the simulation. An inactivation rate equation accounted for the reaction between the spores and the disinfectant. The ANN-based metamodels were most successful in predicting the number of viable bioaerosols remaining in an arbitrary enclosed space. Sensitivity analysis showed that the mass fraction of the disinfectant, inactivation rate constant, and contact time had the most influence on the inactivation of the spores.

IMPLICATIONS

This investigation presents a framework for the development of user-friendly models; metamodels for the prediction of the number of viable spores remaining in an indoor room during disinfection from accurate but time-consuming CFD studies. During any decontamination event, to know when to stop pumping in the disinfectant and to know what level of log reduction of the spores have been achieved before even starting decontamination would provide valuable guidance. The neural network based metamodels can be applied to obtain quick and relatively accurate answers. This would be necessary when immediate information is required during emergencies.

INTRODUCTION

The deliberate release of anthrax spores, which occurred shortly after the events in September 11, 2001, resulted in 22 cases of anthrax infection and five deaths. This was a major factor in increasing national concerns for biological agent releases in the interior of building spaces.1 Biological agents have a high potential for causing massive civilian casualties.2 Most biological warfare agents, whether replicating (infectious) agents or derived toxins, have been weaponized for an aerosol assault.3 This mode is most likely to inflict widespread disease.3 Chlorine dioxide gas was used as the disinfectant during the decontamination of the Capitol Hill, Washington, DC.4 After the 2001 anthrax attacks, the decontamination effort costs were 100 million dollars.5 Decontamination of large indoor spaces and buildings following release of biological agents is challenging, as the response to the fall 2001 anthrax-release events indicate. The ability to efficiently and rapidly decontaminate rooms/buildings is limited by the lack of quantitative understanding of the behavior of biological agents and decontaminants in relation to building geometry, airflow pattern, and surface properties. Room decontamination also has significance in medical situations and laboratory design.

Numerous studies have been done regarding indoor airflow and particle transport in ventilated buildings. Lu et al.6 conducted experiments measuring particle concentrations in the two-zone room and validated their developed fluid dynamics model (CFD) model. Zhao et al.7 looked at the effect of displacement and mixing ventilation on particle concentration and deposition rate. CFD was used to compare traditional air displacement and under floor air distribution systems, which provided more detailed information at a lesser expense than experimental studies.8 Others have developed CFD models for larger buildings with multiple rooms9.10 and have done experimental studies looking at the effect of room furnishings on particle deposition.11 Historically, CFD modeling efforts have focused primarily on the simulation of airflow and dispersion of particles in apartment and buildings, not on disinfectant transport, interaction with bioaerosols, and decay in indoor air. Reshetin and Regens2 conducted simulation modeling of anthrax spore dispersion in a 50-story, high-rise building. The modeling results showed that a relatively small volume of aerosolized spores has the potential to disperse quickly throughout the high-rise building. Stuart and Wilkening5 investigated different mathematical models used to predict environmental degradation of the viability of potential bioweapons agents such as anthrax spores. The models behaved similarly at short times, <30 min; for longer-time phenomena, the different models gave vastly different predictions.5 Experimental studies on ClO₂ gas12 –15 have examined its efficiency in decontamination and have tried to understand its interaction with indoor surfaces.

Past studies provide detailed information and insight into the multiple factors influencing the dispersion of particles in a room, about spore dispersion and decay in buildings and about disinfectant transport and decay. The modeling results, CFD or otherwise, indicates how many particles remain in the room over a certain time period and do not incorporate disinfection. During any decontamination event, to know when to stop pumping in the disinfectant and to know what level of log reduction of the spores have been achieved before even starting decontamination would provide valuable guidance. But to obtain such accurate and quick information on the degradation of the spores would require extensive experimental work and the development of accurate but time consuming models such as a CFD model, specific to the scenario. This would be time consuming and require highly skilled personnel. This paper develops a framework for devising a user-friendly model that incorporates all the different aspects of the problem and is capable of predicting reasonably accurately and quickly how many viable spores remain during a decontamination event and in an emergency in an indoor space. Such a user-friendly model is known as a metamodel. “A metamodel serves as a surrogate or substitute for the more complex and computationally intensive simulation model” and was first proposed by Blanning.16 Metamodeling allows for wider exploration of the problem parameters, improves the understanding of the model to be generated, and enables further studies for solution optimization systems.17 ^– 19

In this study a comprehensive three-dimensional CFD model was developed and applied to numerically simulate the transport of the bioaerosols (e.g., spores of Bacillus anthracis) and their inactivation by a decontaminant (chlorine dioxide). Linear, quadratic, and artificial neural network (ANN)-based metamodels were developed and compared for predicting the number of viable bioaerosols remaining in an arbitrary enclosed space after sometime when a specific disinfectant has been applied. ANNs have increasingly been used instead of the conventional regression models to develop metamodels. This is because ANNs offer an alternative approach to conventional mathematical models for modeling complex systems. It supports multiple inputs and outputs and is capable of efficiently approximating any nonlinear function to an arbitrary degree with a single hidden layer. Recently ANN were trained and tested for optimizing ventilation systems in office environments.20 They were also applied for the prediction of particle number and dispersion in any arbitrary indoor environments and performed better than linear and quadratic model forms.21 They have been used to model the disinfection of waterborne bacteria,22 for predicting ground level ozone giving residents forecasts in advance23 –26 and for predicting the extent of groundwater contamination.17

METHODOLOGY

The development of the metamodel involved the following steps: (1) identifying the input variables and the objective(s) (problem description); (2) choosing an experimental design technique to generate scenarios for efficient sampling of the design space (design space of the input variables); (3) conducting numerical experiments using CFD (development of the CFD model); and finally (4) testing and choosing the form of the meta model (determination of the metamodel form).

Problem Description

The geometry of the three-dimensional (3D) room and the aerosol release are shown in Figure 1, which is representative of the problem scenario. The problem consists of three parts: airflow, aerosol transport and inactivation, and disinfectant transport and decay. To define the problem, 18 dimensional variables were selected. The details of the variables have been given in Table 1.27 The assumptions are that the room is bare of any furnishings or partitions and infiltration and exfiltration is negligible. Air flows in and out of the room through the inlet and outlet located at the opposite ends of the room. Spores are released from an arbitrary location in the room. For a specific period of time the spores are allowed to disperse in the room before a mixture of air and a disinfectant is released through the room inlet.

Table 1. List of variables

	Variable	Symbol	Unit
1	Initial number of the organisms released	N 0	spores/m^3
2	Number of viable organisms in the room after time	N	spores/m^3
3	Length of the room	L	m
4	Width of the room	W	m
5	Height of the room	H r	m
6	Distance between inlet and outlet	I	m
7	Cross-sectional area of the inlet and outlet	A	m^2
8	Velocity of mixture at the inlet	U mix	m/sec
9	Mass fraction of chlorine dioxide
10	Time given for the spores to initially disperse	τ	sec
11	Contact time of spores with chlorine dioxide	τc	sec
12	Inactivation rate of B. anthracis spores	k in	(L/mg)^0.5sec^−1
13	Surface decay rate of chlorine dioxide	k s	sec^−1
14	Bulk decay rate of chlorine dioxide	k b	sec^−1
15	Density of air	ρf	kg/m^3
16	Dynamic viscosity of air	μ	kg/m·sec
17	Density of spore	ρp	kg/m^3
18	Diameter of spore	d p	kg/m^3

Figure 1. Geometry of the room.

In a functional form the relationship between the variables in Table 1 can be written as follows where is the normalized number of viable spores remaining in the room,

(1)

By applying Buckingham Π theorem,28 14 dimensionless groups are obtained from eq 1:

(2)

Through dimensional analysis the dimensionality is reduced to a function that relates the 14 dimensionless groups. Rearranging the 13 dimensionless groups on the right hand side of the eq 2 leads to eq 3, where Q _mix (m³/sec) is the volumetric flow rate and θ (sec) is the residence time of the room calculated from , where V(m³) is the volume of the room. The dimensionless groups and are the same as the groups developed in previous studies by Hoque et al.21,27: the first four groups are the normalized length scales of the room; the last two normalizes the spore properties, particle diameter by the length scale of the room, and the spore density by the fluid density. is similar to the Reynolds number. is the ratio of the time the spores are in contact with only air and the residence time of the room. is the ratio of the time the spores are in contact with the air disinfectant mixture to the residence time of the room. The dimensionless groups account for the relationship between the decay rates of the disinfectant and the inactivation of the spores and the residence time of the room.

(3)

To determine the form of the function in eq 3, accurate numerical investigations are necessary.

Design Space of the Input Variables

The design space of the 13 input dimensionless groups identified above has to be explored to extract the maximum information about the relationship between the output, and the input groups. However due to the computational cost of running the CFD model, the experimental design must efficiently explore the design space. First the boundaries of the design space would have to be determined. To determine the boundaries of the design space, a maximum value and a minimum value have to be assigned to the individual dimensional variables. The focus of the current study was to determine the effects of the variables introduced due to inactivation and decay. To ensure this was done efficiently, other dimensionless groups that did not involve the variables , k _in, k _s, k _b, θ, and Q _mix were fixed at a constant value. The effect of the other dimensionless groups has been reported in a study by Hoque et al.21 The geometry of the room was fixed with dimensions of . The inlet and outlet area, A, was 0.5 m². The minimum mass fraction assigned to chlorine dioxide, , was 500 ppm and the maximum was 2000 ppm, whereas contact time of the spores, τ_c, with the disinfectant was ranged between 3 and 15 hrs. These estimates were based on the concentrations applied and the required length of time when disinfecting the Hart Senate Office Building, mail packages, and office building in Boca Raton, Florida.29 The range for k _s was determined from the estimated deposition velocities of chlorine dioxide on different surfaces.13 The reaction between decontaminants or pollutants and the indoor surfaces is characterized by mass transfer coefficients, which are referred to as deposition velocities.30,31 The range for k _b was estimated based on decay constants calculated from the regression analysis of the data obtained from experimental studies.1223514 The value for k _in was kept constant.12,14 Based on these values, the corresponding dimensionless groups in eq 3 were evaluated and the minimum and maximum levels chosen. Table 2 shows the minimum and maximum values of these dimensionless groups, labeled x ₁ to x₆, thus defining the boundaries of the design space. Sampling points at locations of this design space were next determined.

Table 2. Range for each of the dimensionless groups

		Groups	Maximum		Minimum
1.		L/W		1.67
2.		L/H r		2.08
3.		I/H r		0.6
4.		WH r/A		15
5.	x 1		2000		500
6	x 2	k sθ	3.64		0.00
7.	x 3	k bθ	0.43		0.01
8.	x 4	k inρ^1/2θ	1051		105
9.	x 5	Q mixρ/Lμ	7550		50
10.		τ/θ		1
11.	x 6	τc/θ	150		3
12.		L/d p		5.00 × 10^6
13.		ρp/ρf		861

Typically for physical experiments, classical experimental methods have been applied to determine the points or levels at which experiments would be conducted.18,32 For computer experiments, space-filling designs have been suggested as the better alternative.32 Four types of space-filling designs have been used more often in the literature: orthogonal arrays, Latin hypercube designs, Hammersley sequences, and uniform designs.18 Hammersley sequence sampling (HSS) has been found to provide better uniformity than Latin hypercube designs for a multi dimensional unit hypercube.18,33 It is an example of a low-discrepancy sequence. The concept of discrepancy is applied to choose the appropriate quasi-Monte Carlo sequence, of which the Hammersley sequence is an example. Twenty-five data points were generated for each dimensionless group through HSS.32,33 Table 3 gives details of the data. This was implemented in MATLAB.34 To obtain a more even distribution of the sampling points, two other conditions were imposed when sampling, one the correlation of the generated points for each input group less than 0.3 and the median of the generated points for each group was within 10% of the midpoint of the range of each of the dimensionless groups. The output, N/N ₀ was then determined by conducting numerical experiments at those locations.

Table 3. Design space of the chosen six dimensionless groups

Cases	x 1	k sθ x 2	k bθ x 3	k inρ^1/2θ x 4	Q mixρ/Lμ x 5	τc/θ x 6
1	530	0.73	0.15	191	693	24
2	590	0.91	0.09	376	882	62
3	650	0.40	0.26	815	435	66
4	710	2.91	0.06	526	538	87
5	770	2.27	0.34	143	936	70
6	830	1.37	0.24	332	378	129
7	890	2.32	0.20	933	756	6
8	950	2.33	0.04	260	951	134
9	1010	3.06	0.25	879	473	8
10	1070	0.29	0.17	965	635	49
11	1130	0.03	0.35	370	924	90
12	1190	0.68	0.07	199	670	111
13	1250	0.69	0.21	749	906	132
14	1310	1.02	0.19	876	1120	9
15	1370	1.35	0.06	992	420	30
16	1430	2.16	0.14	543	161	51
17	1490	1.93	0.40	597	658	84
18	1550	2.35	0.13	175	942	93
19	1610	3.35	0.32	820	588	118
20	1670	0.58	0.08	955	1060	112
21	1730	2.39	0.04	973	336	56
22	1790	1.48	0.25	298	714	5
23	1850	0.39	0.37	433	988	136
24	1910	3.49	0.13	569	351	17
25	1970	0.03	0.13	667	843	96

Development of the CFD Model

The geometry of the CFD model is the 3D room as shown in Figure 1. The dimensions assigned to the room were along the x-axis (longitudinal), y-axis (transverse), and z-axis (vertical), respectively. Air enters the room through the supply section having dimensions, located at a height higher than the exit section , which is located at the opposite wall of the room. The CFD model was utilized to predict flow fields, spore trajectories, disinfectant concentration, and decay and inactivation of the spores. The complete model comprised submodels for fluid flow and particle transport; disinfectant transport and spore inactivation. To ensure a more uniform gradient of the disinfectant concentration in the room, a momentum source was introduced for increasing turbulence and to enhance mixing.35 This was simulated by adding force terms, axial, normal, and circumferential to the momentum equation. The force was calculated from the equation, where is the mass displaced and is the velocity difference between the fluid velocity and the rotational velocity specified for the momentum source.35 The numerical simulations were conducted using CFD-ACE+,36 which allowed transient particle calculations.

Airflow Model

The 3D indoor turbulent flow was modeled here by Large Eddy Simulation (LES). LES is a technique of complexity between direct numerical simulation (DNS) and Reynolds averaged Navier-Stokes (RANS). DNS is the most accurate of the three methods but it requires the most computing power (high-speed and large-capacity computer) for airflow simulation in a building.37 RANS modeling applies different turbulence models to calculate the Reynolds stresses and close the system of Navier-Stokes equations. It generates the mean velocity field and not the instantaneous velocity field. The latter influences particle dispersion significantly. LES generates instantaneous flow information required for particle simulation.37 ,38 In LES, small-scale turbulence is filtered out from the Navier-Stokes equations, and a model is used to evaluate small scales. The filtered Navier-Stokes equations are solved for the large-scale motion, which is responsible for most of the momentum and energy transport. The filtered or resolved quantity is denoted by an overbar. The small scales are from subgrid-scale models, which in turn influence the large eddies. For incompressible flow of a Newtonian fluid, the filtered Navier-Stokes equations are

(4)

(5)

The effect of the small scales appears through the subgrid-scale (SGS) stress term, .39 The most commonly used subgrid-scale stress model is the Smagorinsky model, which is an eddy-viscosity model of the form, . The subgrid-scale stresses are related to the large-scale strain rate tensor by . The eddy viscosity, is obtained from the algebraic model, eq 6, where is the grid size and is the Smagorinsky constant. In the Smagorinsky model, the eddy-viscosity coefficient is maintained constant for the entire flow domain.

(6)

(7)

Particle Model

The particle trajectories were computed using the Lagrangian approach. The trajectory of the particle was determined by solving the following set of equations where is the velocity of the particle of mass m at location .

(8)

(9)

is the net gravitational force on a sphere and is represented by eq 10. The drag force, , is represented by the eq 11, where drag coefficient, velocity of airflow, and Cunningham correction factor.

(10)

(11)

The Cunningham correction factor is a function of the Knudsen number and is determined from , where λ is the mean free path. For particle motion with Reynolds' number less than 0.1, , when the Reynolds' number is from 1000 to , . In between, is calculated from: .40

Disinfectant Transport and Decay

Disinfection is a process designed for the deliberate reduction of the number of pathogenic microorganisms. In this study the transport and decay of ClO₂ disinfectant has been simulated when it is released into the inlet air stream of a room. Disinfectant mass fraction was calculated from the convective/diffusive mass transport equation:

(12)

where is concentration of the species, ClO₂, μ _t is turbulent viscosity, is the velocity vector, is the Cartesian coordinate vector with index , D is the molecular diffusivity of the species. Chlorine dioxide decay in the room space, , was modeled as a pseudo-first-order reaction, eq 13, where k _b (sec⁻¹) is the decay rate constant:

(13)

Chlorine dioxide reaction with surfaces was modeled as a first-order reaction, as shown in eq 14.

(14)

The rate of reaction based on unit surface area available for the current case was written as

(15)

where volume of the indoor space (m³), surface area (m²), over all gas transfer coefficient (sec⁻¹), and k = mass transfer coefficient or deposition velocity (m/sec).

Kinetics of Inactivation

The disinfection process involves the contact of a biological agent (for example, B. anthracis spores) with a dose of disinfectant (ClO₂) for a sufficient time period to achieve a desired level of inactivation. Throughout the disinfection contact period, both the disinfectant and viable spore concentrations would decrease. An inactivation rate equation for the reaction between viable anthrax spores and the disinfectant were derived based on the study of kinetic reactions involving chlorine dioxide at a fixed temperature and relative humidity.14 This study reported kinetics consistent with a simple Chick-Watson relationship.41

(16)

where r _in is the reaction rate defined above, T is the temperature, and RH is the relative humidity. The concentration C (decontamination agent) is a function of position within the contaminated zone and f is the viability of the anthrax spores. Thus the inactivation rate equation for the reaction between anthrax particles and the disinfectant has been generalized in this study as

(17)

where k _in is reaction constant, n is 0.5, C is the concentration of the disinfectant at the location of the spore, and f is the viability of the anthrax spores. Equation 17 was included into the simulation to track the viability of the spores.

Boundary Conditions and Numerical Methods

The control-volume method was employed to discretize the governing equations. A blended second-order upwind scheme was applied for the convective-diffusive terms and the implicit Euler scheme was used for the temporal terms in the governing equations. Dirichlet boundary conditions were specified at the inlet for the velocity components. The outlets were designated as pressure boundaries, with Dirichlet conditions specified for pressure and Neumann conditions (i.e., zero normal gradients) specified for all other dependent variables. The boundary conditions at the walls of the room represented no-slip conditions. ClO₂ enters the room mixed with the inlet air. Inlet mass fraction of ClO₂ was specified based on the case. The spores were initialized, slightly above the floor with a small upward velocity of 0.001 m/sec. When the spores reached the outlets, the calculation of the trajectories were terminated. When a spore strikes a surface, it was assumed that it lost half of its energy. This assumption accounted for the different types of surfaces, that is, walls or carpeted floor. Below is a summary of the assumptions:

1.	Incompressible and isothermal conditions
2.	No heat and mass transfer between the air and the particles
3.	The coefficient of restitution was set to 0.5
4.	Spherical, solid, and neutral particles
5.	Airborne particles have no influence on the surrounding airflow field, that is, one-way coupling
6.	Negligible particle-particle interactions
7.	Negligible resuspension.

Simulation Steps

Time-dependent flow fields and particle trajectories were simulated for the 25 cases generated by HSS. For each case, the inlet velocity, chlorine dioxide mass fraction, and kinetic rate parameters were first determined, and the air changes per hour were kept within a range of 1 to 10. This range encompasses the permitted air changes per hour for offices, residences, and hospitals. At the end of each simulation, the change of the viability of the spores was determined from the moment spores were released to the moment the simulation was terminated. The steps in the simulation were as follows:

1.	The flow field in the room was first established for one residence time.
2.	100,000 spores were released and tracked for one residence time.
3.	Chlorine dioxide was then mixed with the inlet air and released into the room and the momentum source turned on.
4.	The simulation was continued until at least a 5 log reduction of the viability of the B. anthracis spores was predicted.

For all cases, nonuniform grids were used with denser grids near the walls. The smallest mesh size was 0.01 m. The time step was 0.05 sec. The 25 cases needed 3 months of computational time to complete requiring on an average, 80 hrs for each case on a personal desktop computer with a 3-GB RAM and 3-GHz processor.

Determination of the Metamodel Form

Determining the architecture of the metamodel is an iterative process that involves choosing a form of the metamodel to represent the behavior of a given phenomenon. Then the parameters in that form are calibrated to the data to determine the best fit, and the generalizability (e.g., by cross-validation) assessed. Three metamodel forms were tested. They were the linear model, quadratic model, and neural networks. The sum of square of errors, , was calculated for each model form where e is the difference between the metamodel output and the CFD simulation output.

Linear and Quadratic Models

Equation 18 represents the multiple linear regression model where are the predictor variables, y is the response, 's are the regression coefficients, and is a random variable assumed to be independent and identically distributed as and is referred to as the random error term of the model.

(18)

The regression coefficient represents the expected change in response y per unit change in x when all the remaining independent variables are held constant. The linear model is linear because the expected value of y is a linear function of weighted sums of the 's. The linear regression model is attractive because the value of the 's can be estimated by analytical solution through the method of least squares. The quadratic form of the model as shown in eq 19 also includes the squared terms of the predictor variables.

(19)

Neural Networks

Neural networks grew out of research attempting to mimic the problem solving capacity of biological neural systems via learning and generalization. There are different model types of neural networks.42 The feed-forward multilayer perception (MLP) model type is the most widely used network and it consists of an input layer, one or more hidden layers, and an output layer. In each layer a specific number of nodes with an activation function are placed. Mathematically, a 3-layer ANN with n, m, and p as the number of input, hidden, output nodes, respectively, is based on the equation:

(20)

where are the output values and are the input values of the network, are the connection weights between the input and the hidden layers, are the connection weights between the hidden layer and the output layer, and f is the activation function. The activation functions could be linear, sigmoid, or hyperbolic tangent. The output of one node contributes to the input to the nodes in the next layer. The weights and biases are determined by optimization methods.17 ,42 The Levenberg-Marquardt numerical optimization technique was employed to update the weights and biases in the network.34,43

Model Assessment

To assess goodness of fit of a model, the following diagnostic plots were employed44: standardized residual versus metamodel-predicted values, y_k , and CFD results versus metamodel-predicted values, y_k . If the plot of model residuals versus model predictions shows a distinct nonrandom pattern then the equation or model is inappropriate. A plot of observed values (CFD results) versus metamodel predictions should be close to a 45° line if the model is accurate.44

For models with differing number of parameters, it has to be determined whether an increase in model parameters is justified by an improvement in fit. To determine this, an F test45 was performed on the basis of

(21)

In eq 21 df is the degrees of freedom, that is, number of data points minus number of parameters (subscript “1” refers to fit with fewer parameters). Critical values for the F distribution were obtained by consulting a standard table using and degrees of freedom. If , then it was concluded that increasing the number of parameters improved the network performance significantly.45

Cross-validation

One of the principal characteristics of a neural network is that if it is provided with a set of inputs and known outputs for a given system, the network organizes itself internally, that is, alters the various weights such that it predicts the expected output for any other given input. The internal process of self-organization of the system is training. A network is overtrained when the network focuses on the characteristics of individual data points rather than capturing the general patterns present in the entire training set. If a network has learned well, then it should classify an unseen pattern correctly and the net is said to generalize well. Cross-validation is a procedure applied to avoid overtraining and for selecting the optimal network. The cross-validation approach consists of a cyclic allocation of the data for training and testing and allows the use of all the data set for testing.46 ,47 The network is trained for a select set of data and tested against the remaining data set. The procedure is repeated until all the data are used. The SSE is computed for both training and testing for optimum convergence. Initially the SSE decreases with increasing number of hidden nodes for both training and testing. But when the network becomes overtrained, the SSE for the testing will start to increase, indicating that increasing the number of nodes is now reducing the generalizing ability of the network. At this point, the number of nodes in the hidden layer is not increased any further.22 ,48 The cross-validation method was applied to avoid overtraining and to select the final model form for the neural networks.

RESULTS AND ANALYSIS

Computational Results

The CFD model was validated against experimental data6 ,49 and was found to be in excellent agreement. Details of the validation studies of the CFD model have been given in Hoque et al.27 Grid sensitivity analysis was conducted to obtain grid-independent solutions.35 Table 3 lists the values of the six dimensionless groups for the 25 cases. For all cases, the rotational velocity of the momentum source was kept constant.35 For each case, the decrease of viability of the spores was tracked over time and recorded for every time step. In all 25 cases, the simulations were terminated well before the anticipated time determined from as greater than 5 log reduction was achieved before τ_c. For example, the theoretical mean residence time for case 14 was 728 sec and for case 22 it was 1269 sec. Spores were tracked for one residence time and then the disinfectant released. In Figures 2 and 3, the plot of the normalized number of viable particles remaining in the room over time for cases 14 and 22 has been shown.

Figure 2. Plot of normalized number of viable particles remaining in the room over time; ClO₂ was injected at 0 sec (case 14).

Figure 3. Plot of normalized number of viable particles remaining in the room over time, ClO₂ injected at 0 sec (case 22).

In the plots, 0 sec indicate when ClO₂ was released into the room. The plots have two lines, one showing the behavior of the spores when no inactivation is happening, that is, spore loss is only due to the airflow and the second line showing loss of spores due to both inactivation and convective flow. In both cases, the combined effect of inactivation and convective flow has the dominant effect, with a sharp decrease of the number of viable spores once disinfection has started. For case 14, nearly 5 log reduction was obtained within 60 sec of ClO₂ release, whereas for case 22, a 5 log reduction was obtained after 130 sec. However, was calculated to be 4 and 11 hrs for cases 14 and 22, respectively. Hence, instead of , for every case was recorded for every time step, where is the actual contact time between the disinfectant and the spores.

By recording 960 data points were obtained from these 25 cases. To map out the relationship and to develop predictive metamodels, the statistical and neural network toolbox of MATLAB34 was applied.

Linear and Quadratic Metamodel

Equation 22 is the linear metamodel. The input dimensionless groups were log transformed, as well as the output, the normalized number of viable particles in the room, . x ₆ represents . The SSE for the linear model was 76222 with an R ² of 0.55.

(22)

Figure 4a compares the predicted results from the linear model to the CFD-simulated results. Most of the data are not around the diagonal, indicating a lack of predictive ability by the linear model. Figure 4b is the plot of the standardized residuals against the predicted results from the model. The plot shows a trend and the data points are not evenly distributed about the horizontal axis at zero.

Figure 4. Linear model results. (a) CFD results against predicted results; (b) standardized residual plot.

Equation 23 is the quadratic metamodel. The R ² for the quadratic model was 0.7, an improvement over the linear model. Figure 5a and b show the residual plot and the plot of the predicted results from the quadratic model to the CFD results. The plots are similar to the plots obtained for the linear model, indicating the quadratic model is not suitable either.

Figure 5. Quadratic model results. (a) CFD results against predicted results; (b) standardized residual plot.

(23)

Neural Network

The neural networks developed have six nodes in the input layer corresponding to the six input dimensionless groups. One node is present in the output layer corresponding to the log normalized number of viable spores remaining in the room, . The chosen activation function was tangent-hyperbolic in the hidden layer(s). The purpose of training and testing the neural network is to identify the architecture that would be appropriate for use and to provide an assessment of the neural network performance. In this study two architectures were studied: one with single hidden layer and the second with two hidden layers. Theoretically the multiple hidden layer networks are not more powerful than single hidden layers but empirically networks with multiple hidden layers sometimes generalize better, with fewer parameters.50,51 For all cases, the neural network was trained until the sum of square of errors reached an assigned minimum value of 10⁻⁶ or remained constant at a low level for a number of iterations, 250,000. Increasing nodes leads to increasing number of adjustable parameters (weights and biases in the network). The maximum number of nodes included was such that the number of adjustable parameters did not exceed the number of data points in the training set.

Architecture with One Hidden Layer

The SSE for the neural network with one node in the hidden layer, ANN_6-1-1 was 11,107 with an R ² of 0.91. The number of adjustable parameters was nine in the neural network same as in the linear model. Figure 6a and b show the plot of the neural network predictions for with two nodes (R ² = 0.99) and three nodes (R ² = 0.995) in the hidden layer, respectively.

Figure 6. Training results (a) for ANN_6-2-1 and (b) for ANN_6-3-1.

Figure 7a and b show the corresponding residual plots. The performance of the neural network is significantly better than either the linear or the quadratic model for all the cases.

Figure 7. Standardized residual plot (a) for ANN_6-2-1 predictions and (b) for ANN_6-3-1 predictions.

To avoid overtraining the network and to determine when to stop increasing the number of nodes, the cross-validation method was applied. The steps applied to do the cross-validation were as follows: out of the 25 cases in the set, data belonging to 24 cases were cases were used to train the networks and were then validated against the remaining data of one case. The SSE was recorded. This was repeated for all 25 cases and the total sum of squares calculated by adding the individual sum of squares. Table 4 shows the results of the training and the cross-validation. In Figure 8, the SSE for both training and cross-validation as the number of nodes was increased from one to five in the hidden layer has been plotted. When the number of hidden nodes increased from three to four, the cross-validation sum of squares increased and continued to increase when the number of hidden nodes was five. Hence the optimum network for the architecture with one hidden layer was ANN_6-3-1 with 25 adjustable parameters.

Table 4. Change of SSE with increasing number of neurons for architecture with one hidden layer

Neurons	Adjustable Parameters	Training SSE	Cross-validation SSE
1	9	11107	27209
2	17	1583.5	25740
3	25	677.6	22989
4	33	275.2	121780
5	41	209.2	282390

Figure 8. Change of sum of square of errors with increasing number of neurons; architecture with one hidden layer.

Architecture with Two Hidden Layers

For the architecture with two hidden layers, the number of nodes in the hidden layers was adjusted such that the number of nodes in the second hidden layer was always less than or equal to the number of nodes in the first hidden layer. It has been shown that neural networks perform better when the second hidden layer does not have more nodes than the first hidden layer.21,51 The neural network, with one node in both hidden layers, ANN_6-1-1-1, has 11 adjustable parameters. The SSE obtained was 10,989 (R ² = 0.91). Figure 9a and b show the plots of the network predictions for ANN_6-2-2-1 and ANN_6-3-2-1. Figure 10a and b show the corresponding residual plots.

Figure 9. Training results (a) for ANN_6-2-2-1 and (b) for ANN_6-3-2-1.

Figure 10. Standardized residual plot (a) for ANN_6-2-2-1 predictions (b) for ANN_6-3-2-1 predictions.

The SSE decreased to 649 (R ² = 0.995) when the number of nodes in both hidden layers was increased to two. The adjustable parameters for this network, ANN_6-2-2-1, were 23 less than the number of parameters in the single hidden layer architecture, ANN_6-3-1, which has 25 but a higher sum of squares, 677.6. The architecture with two hidden layers appeared to be performing better than the architecture with one hidden layer. The number of nodes in the hidden layers was increased until there were four nodes in each hidden layer.

Cross-validation was then applied to determine the optimum network architecture. Table 5 lists the sum of squares for both training and cross-validation for all the two hidden layered architectures tried. Figure 11 plots these values against increasing number of nodes in the second hidden layer. The plot shows that increasing the number of nodes in the first hidden layer resulted in an overall decrease in the sum of squares irrespective of the number of adjustable parameters. The cross-validation results show that for networks with three and four nodes in the first hidden layer, overtraining occurred when the nodes in the second hidden layer were increased after a certain point. The cross-validation decreased from network ANN_6-2-1-1 to ANN_6-2-2-1. Based on the cross-validation results, the two networks that had comparatively lower sum of squares were ANN_6-4-2-1 and ANN_6-4-3-1. An F test showed that increasing the number of adjustable parameters from 41 to 47 when going from ANN_6-4-2-1 to ANN_6-4-3-1 resulted in a significant decrease to the sum of squares. Hence the optimum network chosen was ANN_6-4-3-1.

Table 5. Change of SSE with number of neurons (architecture: two hidden layers)

Neurons
First Hidden Layer	Second Hidden Layer	Adjustable Parameters	Training SSE	Cross-validation SSE
1	1	11	10989	65533
2	1	19	1578	58883
2	2	23	649	39165
3	1	27	445	23940
3	2	32	379	16299
3	3	37	308	3.82 × 10^10
4	1	35	287	61222
4	2	41	207	14059
4	3	47	189	15863
4	4	53	161	4.72 × 10^6

Figure 11. Change of sum of square of errors with increasing number of neurons in the second hidden layer; architecture with two hidden layers (number of neurons in the second hidden layer the number of neurons in the first hidden layer).

Between the final two architectures ANN_6-3-1 and ANN_6-4-3-1 the architecture with two hidden layers had a lower sum of squares for both training and cross-validation. Based on this and the very improved performance of the networks over the linear and quadratic metamodels, the final form of the metamodel chosen was the neural network ANN_6-4-3-1.

Sensitivity Analysis of the Input Dimensionless Groups

Figure 12 shows the sensitivity of the response, , when the magnitudes of each of the dimensionless groups were varied independently applying the neural network ANN_6-4-3-1. The magnitudes of the dimensionless groups were varied within the logarithm of the specified range (Table 3). The change of with increasing has been plotted in Figure 12a. As the concentration of the disinfectant was increased, the number of viable spores decreased. The plots in Figure 12b–f are for the other five dimensionless groups. For the dimensionless groups with the decay rate constants k _sθ and k _bθ, the number of spores increases with increasing decay. This occurs because more of the disinfectant is removed due to self-decay or due to surface interaction resulting in lower availability of disinfectant for the inactivation of spores. In case of , the number of viable spores decreased with increasing inactivation rate. Also, increasing contact time resulted in decreasing number of spores, as seen in Figure 12f. With increasing flow rate, decreases and then the plot levels off for the dimensionless group at the point where the number of viable spores in the room is approximately 20, indicating at some point the spores are not in the flow stream and are not being removed by increasing airflow.

Figure 12. Plots of sensitivity of the response of the output log(N/N ₀) applying ANN_6-4-3-1.

CONCLUSIONS

Numerical experiments were conducted using CFD for generating data to develop metamodels capable of predicting the number of viable spores remaining in a room after starting disinfection using chlorine dioxide. A CFD model was developed to simulate the inactivation of B. anthracis spores released in a room by chlorine dioxide including the variables mass fraction of chlorine dioxide, , surface and bulk decay rate of chlorine dioxide, , inactivation rate constant, , and contact time, , between spores and disinfectant. A momentum source was used in the room to enhance mixing of the disinfectant in the room. The CFD model was applied to simulate 25 cases generated applying HSS method.

The linear, quadratic, and neural-network-based metamodels were explored. The performance of the neural networks was significantly better than either linear or quadratic metamodel. Two types of architectures were explored for the neural networks: one with a single hidden layer, and the other with two hidden layers. Cross-validation was applied to determine the maximum number of nodes before overtraining occurred. ANN_6-3-1 was the best neural network for the single hidden layer architecture. For the two hidden layer architecture, the best neural network was ANN_6-4-3-1 based on the cross-validation sum of squares. The SSE obtained from the architecture with two hidden layers was lower than the SSE obtained from the architecture with one hidden layer. Hence the final form of the metamodel chosen was ANN_6-4-3-1. Sensitivity analysis studies were then conducted using this network. The analysis showed that the bulk and surface decay rate had similar sensitivity towards the log reduction of the viability of the particles. Mass fraction of chlorine dioxide, inactivation rate constant, and contact time had the most influence, whereas flow rate appeared to have the least impact on the inactivation of the spores.

The study outlines a systematic method for developing metamodels for predicting inactivation of spores from CFD simulations. The results show that this method can be applied to gain a better understanding into the disinfection kinetics in the indoor spaces and with spores. The metamodel developed can be used for providing information during emergencies. More information regarding the surface chemistry of chlorine dioxide would help in increasing model accuracy. The methodology can also be expanded by including more variables into the model that would account for additional complexities in the space such as partitions, multiple zones, etc. The presence and type of furniture can be incorporated into the model by including parameters such as area occupied by furniture, different surface roughness values, and possible rate constants, etc. The modeling could also include as a variable the range of impact of the momentum source on the mixing of the disinfectant and the influence of the impact for different room geometries.

ACKNOWLEDGMENTS

The work described in this paper was supported by the L. D. Betz Environmental Engineering Endowment at Drexel University.