Laboratory Aerosol Hazard Assessment Utilizing Material Balance Models

Abstract

Models that are simple to use and interpret and that predict aerosol concentration and integrated exposure dose are needed for local, timely, and cost-effective evaluation of risks associated with biocontainment operations. A material balance model is presented with explicit assumptions that provide a safety context for assessing potential laboratory operating and accident scenarios (1) to support determination of aerosol concentration when limited experimental data are available or (2) to assist in the interpretation of experimental results. The model is used to assess recent experimental data for large animal primary containment systems in animal biosafety level 4; it is also applied to estimate integrated dose from particle concentration for assessments of potential exposure of laboratory workers. The model is further shown to reproduce previously published laboratory accident results.

Aerosols are composed of particles that may include infectious or toxic material that are hazardous to laboratory workers. Models that are simple to use and interpret that predict aerosol concentration and integrated exposure dose are needed for local, timely, and cost-effective evaluation of aerosol risks associated with biocontainment operations. Material balance or mass balance models assume that the number of particles in a volume at some time in the future is the number in the volume now, plus the number that enters the volume, and minus the number that leaves the volume. These assumptions are a model of what happens to the particle concentration and not a simulation of the location of individual particles, the distribution of flows within a volume, particle-size dependencies, or deposition. Specifically, it is assumed in the material balance model that when particles enter or leave a volume, the concentration is instantaneously normalized across the entire volume; that is, the aerosol concentration is always uniform across the volume. Aerosol decay is not included in the model but could be approximated as a virtual volume where particles deposit, effectively removing deposited particles from the volume being modeled.

As an example, consider a worker in a laboratory when an aerosol release occurs and the concentration of aerosol in the room increases with time. With the material balance assumptions at a specific point in time, regardless of where the worker is located in the room, the concentration is the same. Provided that biocontainment laboratories have high air exchange rates, this assumption may be true over some time period due to mixing. However, if the worker is located relatively far from the aerosol source during the release, it can be assumed that the material balance assumptions provide an estimate of concentration that is higher than the actual concentration near the worker. If the material balance model–based concentration is used to assess integrated dose as a measure of occupational exposure, it would therefore be a conservative estimate; that is, the actual dose is likely lower than the estimate from the model. Similarly, if the worker is proximal to the release, we would expect the actual concentration to be higher than an estimate produced using the material balance assumption. It would not be satisfactory to utilize an assessment of potential occupational exposure based on this assumption for the worker in proximity to the source. One potential option to provide a more conservative estimate for the exposure of the worker in proximity to the source is to redefine the dimensions of the space being modeled to approximate the physical distance between the worker and the aerosol source—that is, the closer to the source the worker is, the smaller the volume used in the model. Due to the simplicity of the underlying assumptions of the material balance model, it is feasible to make these types of trade-offs.

A general analytic framework called the material balance aerosol hazard assessment (MBAHA) model was developed according to the described material balance assumptions for multiple connected volumes (eg, biocontainment rooms, biosafety cabinets, primary containment caging). Solutions for aerosol concentration and integrated dose through the MBAHA model involve well-understood mathematical expressions that can be evaluated analytically, albeit somewhat tediously, and numerically with a spreadsheet. The model analyses and spreadsheet implementation address several important biosafety applications in the risk management of aerosol hazards in biocontainment: concentration as a function of time across multiple interconnected volumes, integrated dose for potential occupational exposures, comparison of experimental data to MBAHA model predictions, and MBAHA model analyses of laboratory configurations of interest that have not been performed experimentally.

The description of the material balance model begins with a single volume but utilizes the notation needed for multiple volumes. To demonstrate its utility, the single-volume MBAHA model is applied to recent experimental aerosol data for an animal biosafety level 4 (ABSL-4) laboratory room presented by Landon et al ¹ to assess aerosol concentration and decay rate. The single-volume MBAHA model is also compared favorably to previously published models in the field, including estimation of integrated dose associated with laboratory accidents. A multivolume MBAHA model is described as well as a matrix-based framework for multiple interconnected volumes and aerosol sources. The equations are solved analytically for aerosol concentration with 2 volumes; the concentration is integrated analytically over time to estimate dose; and the MBAHA model is applied to experimental data on primary containment caging in ABSL-4 provided by Landon et al. Other conditions of interest that were not addressed experimentally by Landon et al were investigated, including aerosol flow out of the primary containment caging into the ABSL-4 room. Sufficient mathematical details are presented for implementation of the models in a spreadsheet. However, mathematical derivations are not essential to applying the model; that is, some may want to skip derivations, especially in a first reading.

Material Balance Method for a Single Volume

The aerosol concentration c₁(t) as a function of time t for a single room with volume v ₁, exiting airflow r ₁₀, and aerosol source s(t) with material balance assumptions is given by equation 1 (EQ-1):

$c_1\left(t\right)=c_1\left(0\right)e^{-\frac{r_{10}}{v_1}t}+\underset{0}{\overset{t}{\int }}{e}^{-\frac{r_{10}}{v_1}\left(t-\tau \right)}s\left(\tau \right)d\tau \equiv \vec{b}\left(t\right)+\vec{\sigma }\left(t\right)$

EQ-1

with the 2 exponential terms weighted by the initial concentration and by the source.

An aerosol source that instantaneously releases into the volume can be modeled mathematically as the generalized function, δ(t), which is zero everywhere except t = 0 and ${\displaystyle \int f\left(\tau \right)\delta \left(t-\tau \right)d\tau }=f\left(t\right)$ , and is referred to as the delta function. For $s^{\delta }\left(t\right)=s \delta \left(t-t_1\right)$ ,

${\sigma }^{\delta }\left(t\right)={\displaystyle \underset{0}{\overset{t}{\int }}{e}^{-\frac{r_{10}}{v_1}\left(t-\tau \right)}s\delta \left(\tau -t_1\right)d\tau }=\{\begin{array}{cc}0& t<t_1\\ s{e}^{-\frac{r_{10}}{v_1}\left(t-t_1\right)}& t_1\le t\end{array}$

EQ-2

which is of the same form as the c₁(0) term in EQ-1. So rather than solving for the c₁(0) term (or more generally, $\vec{b}\left(t\right)$ ), we solve for the $\vec{\sigma }\left(t\right)$ term and use a delta function source with s set to c₁(0). The delta source can be used to set initial conditions as well as represent release of an aerosol that happens over a short period (eg, dropping a culture plate and potentially releasing viable microorganisms as an aerosol).

For a rectangular source “[ ]” with $s^{[]}\left(t\right)=s$ for $t_1\le t\le t_2$ and 0 otherwise:

${\sigma }^{[]}\left(t\right)=\{\begin{array}{cc}0& t<t_1\\ \underset{t_1}{\overset{t}{\int }}{e}^{-\frac{r_{10}}{v_1}\left(t-\tau \right)}sd\tau =\frac{v_1}{r_{10}}s\left(1-e^{-\frac{r_{10}}{v_1}\left(t-t_1\right)}\right)& t_1\le t\le t_2\\ \underset{t_1}{\overset{t_2}{\int }}{e}^{-\frac{r_{10}}{v_1}\left(t-\tau \right)}sd\tau =\frac{v_1}{r_{10}}s{e}^{-\frac{r_{10}}{v_1}t}\left(e^{\frac{r_{10}}{v_1}{t}_2}-e^{\frac{r_{10}}{v_1}{t}_1}\right)& t_2\le t\end{array}$

EQ-3

Note that the second equation can be obtained by setting t₂ to t in the third equation. Once the source is done emitting, the exponential decay is at the same rate $r_{10}/v_1$ as previously determined. We use a rectangular source when there are aerosols released continuously or over a finite duration.

Validation with Recent Data from Soft-Walled Enclosures in ABSL-4

The exponential decay rate $0.410\pm 0.012 {\text{min}}^{-1}$ estimated from the experimental data for an ABSL-4 room in Landon et al ¹ is well represented by the model prediction (EQ-2):

$\frac{r_{10}}{v_1}=\frac{31.2{\text{ m}}^3/\text{min}}{74.98{\text{ m}}^3}=0.416\left[{\text{min}}^{-1}\right]$

EQ-4

where v₁ is the volume of the ABSL-4 room and r₁₀ is the high-efficiency particulate air (HEPA) filter flow rate leaving the room. When the soft-walled enclosure HEPA filter was running inside the ABSL-4 room, the rate estimated experimentally was $0.560\pm 0.032 {\text{min}}^{-1}$ . The overall rate can be approximated in the simple one-volume model by using the combined flow rates leaving the ABSL-4 room (31.2) and the soft-walled enclosure (8.5):

$\frac{r_{10}}{v_1}=\frac{31.2+8.5{\text{ m}}^3/\text{min}}{74.98{\text{ m}}^3}=0.529\left[{\text{min}}^{-1}\right]$

The approximation in the model is equivalent to even mixing between the enclosure and the room (ie, simply a larger single HEPA filter for the ABSL-4 room). To more appropriately model the concentration when the soft-walled enclosure is being used requires the 2-volume model, which is presented later.

Dose Calculations and Integration with Earlier Work

From a safety perspective, the cumulative exposure to the concentration (ie, dose) is a better measure of risk than concentration directly. We define dose d(t) as a function of concentration c(t), exposure duration, breathing rate b_r , and retention factor f_R :

$d\left(t\right)=\underset{0}{\overset{t}{\int }}c\left(\tau \right)b_r{f}_{R}d\tau$

EQ-5

Obviously, b_r and f_R could be more complex (eg, time or aerosol particle-size dependencies). Unless noted otherwise, we assign b_R = 15 [L/min] = 0.25 [L/s] and f_r = 1, as assumed in Bennett et al. ² A conservative estimate of dose with no air changes, r ₁₀ = 0, is also usually determined. For our simple 1-volume model of concentration, $r_{10}>0,$ from a delta source at time t ₁ (EQ-2),

$d^{\delta }\left(t\right)=\{\begin{array}{cc}0& t<t_1\\ b_r{f}_{R}s\frac{v_1}{r_{10}}\left(1-e^{-\frac{r_{10}}{v_1}\left(t-t_1\right)}\right)& t_1\le t\end{array}$

EQ-6

For a rectangular source from EQ-3, $d^{[]}\left(t\right)$ is

$d^{[]}\left(t\right)=\{\begin{array}{cc}0& t<t_1\\ b_r{f}_{R}s\frac{v_1}{r_{10}}\left(\left(t-t_1\right)+\frac{v_1}{r_{10}}\left(e^{-\frac{r_{10}}{v_1}\left(t-t_1\right)}-1\right)\right)& t_1\le t\le t_2\\ b_r{f}_{R}s\frac{v_1}{r_{10}}\left(\left(t_2-t_1\right)+\frac{v_1}{r_{10}}\left(e^{-\frac{r_{10}}{v_1}\left(t-t_1\right)}-e^{-\frac{r_{10}}{v_1}\left(t-t_2\right)}\right)\right)& t_2\le t\end{array}$

EQ-7

Table 1 provides a comparison of terminology from Dimmick et al ³ and Bennett et al ² with what is utilized in the MBAHA model. In some cases, the experiment-driven definitions have simplifying assumptions as compared with the MBAHA model.

Table 1.

Comparison of Model Definitions with Terminology Utilized in Other Publications.

Term	Other Pubs	Model	Description
Aerosol concentration	$AC=\frac{SS}{V}{T}_S$	ci(t)	Aerosol concentration (CFU/m3)
Source concentration	$SC=\frac{SS}{SF}$	$c_{0i}\left(t\right)$	Concentration in the source (agent/volume) where volume may be mL, mg, etc (CFU/mL)
Volume and source volume	V	vi	Volume surrounding the source (also referred to as “source volume”)
Spray factor 3	$SF=\frac{SS}{SC}$	$r_{0i}$	An experimentally driven collection of factors that provide rates converting source material to aerosol for different sources and scenarios (mL/s)
Spray factor 2	$SF=\frac{AC}{SC}$	$\frac{c_i\left(0\right)}{c_{0i}}=\frac{r_{0i}}{v_i}$	Removes the time-varying aspects of Dimmick et al 3 by assuming aerosols are instantly dispersed and includes room volume in the factor (mL/m3)
Source strength	${\text{S}}_{\text{m}}=SS=SF*SC$	$r_{0i}{c}_{0i}\left(t\right)=v_1{s}_1\left(t\right)$	Rate (CFU/min) that source material is converted to aerosol
Operation time	TS	t2 − t1	Time that the source is operating

Abbreviation: CFU, colony-forming units.

An assumption of prompt dispersal (eg, spray factor from Bennett et al ² ) across the entire working volume simplifies the models by eliminating time as a source parameter. This is treated as the delta source case in the MBAHA model and utilizes the more general definition by Dimmick et al ³ of source strength S_m for a single volume $r_{0i}{c}_{0i}\left(t\right)$ that includes time with a focus on rectangular sources in our analyses:

${\text{S}}_{\text{m}}\left(\text{t}\right)=SF\hspace{0.17em}\hspace{0.17em}{c}_{01}\left(t\right)=r_{01}{c}_{01}\left(t\right)\left[\frac{\text{CFU}}{\text{minute}}\right]=v_1{s}_1\left(t\right)$

This approach allows spray factor (SF) from Bennett et al ² to map simply to $SF=r_{01}$ .

The parameters k for room changes per hour [AC/H], S_m for source strength, and SF for spray factor found in Dimmick et al ³ map to the MBAHA model as

$r_{10}=k\left[\frac{\text{Rooms}}{\text{hour}}\right]v_1\left[\frac{{\text{m}}^3}{\text{Room}}\right]\frac{1}{60}\left[\frac{\text{hours}}{\text{min}}\right]$

$v_1{s}_1\left(t\right)=r_{01}{c}_{01}\left(t\right)={\text{S}}_{\text{m}}=SF\hspace{0.17em}\hspace{0.17em}{c}_{01}\left(t\right)$

and use the values S_m = 10 000 and v ₁ = 1000 m³. For a source that is continuously spraying (ie, a rectangular source with r ₁₀ > 0, 0 = t ₁ and t < t ₂ from EQ-3) the concentration in a single volume room is

$c_1\left(t\right)=\frac{v_1}{r_{10}}s\left(1-e^{-\frac{r_{10}}{v_1}\left(t-t_1\right)}\right)=\frac{60{\text{S}}_{\text{m}}}{k{v}_1}\left[1-e^{-\frac{k}{60}t}\right]$

EQ-8

The concentration from EQ-8 is plotted in Figure 1 and is identical to the graph of concentration found as Figure 2 in Dimmick et al^3(p260) (note the special case for k = 0, producing a concentration of $\frac{{\text{S}}_{\text{m}}}{v_1}t$ ).

Figure 1.

Utilizing equation 8 to reproduce the concentration results of Dimmick et al.^3(p260) AC/H, air changes per hour; CFU, colony-forming units.

Figure 2.

Comparison of dose models on 2- and 10-minute scales. The delta source is equivalent to the model used in Dimmick et al, ³ but the boxcar source implements the stated assumptions. On a 10-minute time scale, the results are essentially the same. PFU, plaque-forming units.

We also consider the dropped flask accident scenario assessed in Bennett et al.^2(p663) The original concentration is given by the product of the source concentration (5 × 10⁹ colony-forming units [CFU] /mL) and the spray factor (5.3 × 10^-7 mL/m³) or 2650 CFU/m³ = 2.65 CFU/L. The air change rate is 12 h^-1, breathing rate b_r = 15 L/min. Using this concentration as constant over exposure times of 30 seconds and 10 minutes, Bennett et al estimate dose as 19.9 = 2.65*15*30/60 and 398 = 2.65*15*10 CFU, respectively.

In our time-varying model of dose, the delta source matches the assumption that the aerosol instantaneously disperses across the volume of the room so

$c_1\left(t\right)=c_{01}\frac{r_{01}}{v_1}{e}^{-\frac{r_{01}}{v_1}t}=2.65e^{-\frac{12}{60}t}\left[\text{CFU}\right]$

and the associated dose is

$d\left(t\right)=\frac{15\text{*}2.65\text{*}60}{12}\left(1-e^{-\frac{12}{60}t}\right)\left[\text{CFU}\right]$

For 12 air changes per hour, we model that the dose at 30 seconds and 10 minutes is estimated as 11.6 and 171.9 [CFU]. As we show, at 10 minutes the 12 room air changes per hour have reduced the dose by more than a factor of 2.3 versus a room with no air flow (398 [CFU]). The ratio of reduced dose due to increased air flow is less for the 30-second dose (1.7 = 19.9/11.6).

Time to Reach a Specific Dose

Because we have an analytic expression for dose (EQ-6) and $s{v}_1=r_{01}{c}_{01} \left[m^{3}CFU\right]$ , we can solve for the time to reach a specific dose d_D for a delta source at time t = 0 and f_R = 1,

$t_D=-\frac{v_1}{r_{10}}\ln \left(1-\frac{r_{10}}{r_{01}{c}_{01}}\frac{d_D}{b_r}\right)$

EQ-9

or time to reach a specific concentration c_c from occurrence of a delta source (with EQ-2):

$t_C=-\frac{v_1}{r_1}\ln \left(\frac{c_c}{\text{s}}\right)=-\frac{v_1}{r_1}\ln \left(\frac{c_c{v}_1}{r_{01}{c}_{01}}\right)$

EQ-10

Therefore, to determine the time t_f in minutes to remove a fraction f (say 0.9) of the original concentration, use

$t_f=-\frac{v_1}{r_1}\ln \left(1-f\right)$

EQ-11

This can be used to produce a table of the fraction removed as a function of air changes per hour (note that the conversion to minutes requires a scaling of 60), shown in Table 2. In comparing this with Table 3.2 of the Health and Safety Executive Advisory Committee on Dangerous Pathogens,^4(p51) we observe a typo in the 0.9999 column where the values for 0.99999 appear; the correct values are in Table 2 here.

Table 2.

Time in Minutes to Remove Fraction of Initial Concentration with Equation 11 Scaled by 60 for AC/H Units.

	Fraction Removed
AC/H	0.9	0.99	0.999	0.9999	0.99999
6	23.03	46.05	69.08	92.10	115.13
7	19.74	39.47	59.21	78.95	98.68
8	17.27	34.54	51.81	69.08	86.35
9	15.35	30.70	46.05	61.40	76.75
10	13.82	27.63	41.45	55.26	69.08
11	12.56	25.12	37.68	50.24	62.80
12	11.51	23.03	34.54	46.05	57.56
14	9.87	19.74	29.60	39.47	49.34
16	8.63	17.27	25.90	34.54	43.17
18	7.68	15.35	23.03	30.70	38.38
20	6.91	13.82	20.72	27.63	34.54
25	5.53	11.05	16.58	22.10	27.63
30	4.61	9.21	13.82	18.42	23.03
40	3.45	6.91	10.36	13.82	17.27

Abbreviation: AC/H, air changes per hour.

a

Time to remove fraction of initial concentration: t = –60*ln(1 – f) / ACH [min].

Proximal and Distal Exposure

Dimmick et al^3(pp258-259) provide a sample source that generates aerosol for 0.5 minutes, 10⁹ plaque-forming units (PFU)/mL, laboratory room volume v ₁ = 2000 ft³, spray factor 5 × 10^-7 mL/min, b_R = 0.3 ft³/min breathing rate, and f_R = 0.3 retention factor. This provides a source strength of 500 PFU/min = (10⁹ PFU/mL) (5 × 10^-7 mL/min). Two scenarios were considered in the paper:

A worker remains proximal to the source (within 3 ft or a volume of 27 ft³). For a room with no air flow, this results in a local concentration of 9.26 PFU/ft³ = (500 PFU/min)(0.5 min)/(27 ft³). The dose accumulated in 2 minutes is estimated as 1.67 PFU = (9 PFU/ft³)(0.3 ft³/min)(0.3)(2 min).

A worker in the same laboratory volume but for 10 minutes and distal to the source (not within 3 ft) results in a diluted concentration of 0.125 PFU/ ft³ = (500 PFU/min)(0.5 min)/(2000 ft³). The dose accumulated in 10 minutes is estimated as 0.1125 PFU = (1.3 PFU/ft³)(0.3 ft³/min)(0.3)(10 min).

The calculations in Dimmick et al ³ begin once the aerosol source has completed generation and are equivalent to the MBAHA model of a delta source in a laboratory with no flow leaving the volume. Comparing the source output for the assumptions in Dimmick et al ³ with the rectangular source described that has duration of 0.5 minutes in a single volume shows the expected difference during the time that the rectangular source is generating particles (ie, concentration and dose are lower until the rectangular source catches up to the delta source). The delta source assumption therefore provides the more conservative dose estimate, as used in Dimmick et al ³ and shown in Figure 2 with 2- and 10-minute time scales.

To compare the value of being away from a source, the same delta source was utilized and the dose modeled for hemispherical volumes with radii ranging from 1 to 5 ft with a 1-ft increase in the radius of a hemisphere for each calculation (Figure 3). As expected, the closer to a source, the more reduction in dose there is for each foot moved away.

Figure 3.

Comparison of calculated dose received for each foot away from the source. Calculated by increasing the volume of an enclosing hemisphere.

Material Balance Method for Multiple Volumes

Consider a system with n volumes v_i [m³] (eg, laboratory rooms, biosafety cabinets, soft-walled enclosures), each with particle concentration c _i(t) [particles/m³], 1 ≤ i ≤ n. As shown in Figure 4, the volumes are interconnected with airflow r_ij [m³/min] leaving v _i and entering v _j, r_{i 0} [m³/min] leaving v_i through HEPA filtration (ie, a sink for the associated r_{i 0} c_i (t) particles leaving the system), and a source producing particles entering v_i with source concentration c _0i (t) [particles/m³] and source rate (also known as spray factor) r _0i [m³/min] (eg, virus shed by large animals), where t, m, min, and particles are time, meters, minutes, and CFU/PFU/other appropriate particulate measure, respectively. Depending on the application, volumes for the source concentration and source rate may be better expressed in other units (eg, mL) for liquid source material and spray devices. All flow rates in this model are assumed to move volumes of particles without changing pressures. It is also assumed that particles introduced into a volume are instantaneously dispersed across the entire volume—that is, concentration in each volume is uniform.

Figure 4.

Concentration is modeled as n interconnected volumes with particles exchanging pairwise among volumes and each volume having a HEPA filter (ie, sink) and a source. PFU, plaque-forming units.

The number of particles v _i c_i (t + dt) in volume v _i at time t + dt, where dt is small, can be expressed as the balance of material in the system, which must be conserved (particles staying in v _i, leaving to other volumes, entering from other volumes, leaving through a filter, and entering from a source):

$v_i{c}_i\left(t+dt\right)=v_i{c}_i\left(t\right)dt-\sum _{j=1}^n{r}_{ij}{c}_i\left(t\right)dt+\sum _{j=1}^n{r}_{ji}{c}_j\left(t\right)dt-r_{i0}{c}_i\left(t\right)dt+r_{0i}{c}_{0i}\left(t\right)dt$

Rearranging terms, taking the limit as dt goes to 0, and denoting the derivate with respect to time t as ${ċ}_i\left(t\right)$ ,

$v_i{ċ}_i\left(t\right)=-\sum _{j=1}^n{r}_{ij}{c}_i\left(t\right)+\sum _{j=1}^n{r}_{ji}{c}_j\left(t\right)-r_{i0}{c}_i\left(t\right)+r_{0i}{c}_{0i}\left(t\right),1\le i\le n$

This system of n equations can be represented as

$\overrightarrow{ċ}\left(t\right)=A\vec{c}\left(t\right)+\vec{s}\left(t\right)$

EQ-12

by the n × 1 column vectors $\vec{c}\left(t\right)={\left[\begin{array}{ccc}{c}_1\left(t\right)& \cdots & c_n\left(t\right)\end{array}\right]}^T,\overrightarrow{ċ}\left(t\right)={\left[\begin{array}{ccc}{ċ}_1\left(t\right)& \cdots & {ċ}_n\left(t\right)\end{array}\right]}^T$ , and $\vec{s}\left(t\right)={\left[\begin{array}{ccc}{r}_{01}{c}_{01}\left(t\right)/v_1& \cdots & r_{0n}{c}_{0n}\left(t\right)/v_n\end{array}\right]}^T,$ where T is the transpose function, and the n × n matrix A with elements a_ij and $1\le i,j\le n$ :

$a_{ij}=\{\begin{array}{cc}\frac{r_{ji}}{v_i},& i\ne j\\ \frac{1}{v_i}\left(-\sum _k{r}_{ik}+r_{ii}-r_{i0}\right),& i=j\end{array}$

EQ-13

Requiring t ≥ 0, the solution of this system of differential equations is given by

$\vec{c}\left(t\right)=e^{At}\vec{c}\left(0\right)+\underset{0}{\overset{t}{\int }}\hspace{0.17em}{e}^{A\left(t-\tau \right)}\vec{s}\left(\tau \right)d\tau \hspace{0.17em}$

EQ-14

There are established mathematical techniques for solving EQ-14, including eigenvalue approaches; for instance, see Dawkins ⁵ or Luenberger. ⁶ Utilizing eigenvalues λ⁽ⁱ⁾ and n × 1 eigenvectors ${\vec{k}}^{\left(i\right)}, 1\le i\le n,$ defined as $A{\vec{k}}^{\left(i\right)}={\lambda }^{\left(i\right)}{\vec{k}}^{\left(i\right)}$ where the eigenvalues solve the nth order polynomial $det\left(A-\lambda I\right)=0$ , where I is the n × n identity matrix.

$M=\left[\begin{array}{ccc}{\vec{k}}^{\left(1\right)}& \cdots & {\vec{k}}^{\left(n\right)}\end{array}\right]\text{ and }{e}^{\Lambda t}=\left[\begin{array}{ccc}{e}^{{\lambda }^{\left(1\right)}t}& & 0\\ & \ddots & \\ 0& & e^{{\lambda }^{\left(n\right)}t}\end{array}\right]\text{ then }{e}^{At}=M{e}^{\Lambda t}{M}^{-1}\text{ and}$

$\vec{c}\left(t\right)=M{e}^{\Lambda t}{M}^{-1}\vec{c}\left(0\right)+\underset{0}{\overset{t}{\int }}\hspace{0.17em}M{e}^{\Lambda \left(t-\tau \right)}{M}^{-1}\vec{s}\left(\tau \right)d\tau \equiv \vec{b}\left(t\right)+\vec{\sigma }\left(t\right)$

EQ-15

Fortunately, eigenvalue/eigenvector systems can be solved numerically (and are required for systems larger than n = 3). For 1- and 2-dimensional systems, the math to solve is straightforward, though sometimes tedious, and provides insight into the contributions of individual parameters as well as allowing additional analyses like estimation of dose. The results have already been presented for the 1-dimensional MBAHA model which is the multivolume model with n = 1, $a_{ij}=a_{11}=-\frac{r_{10}}{v_1}={\lambda }^{\left(i\right)}$ , and ${\vec{k}}^{\left(1\right)}=\left[1\right]=M=M^{-1}$ . The 2-dimensional (ie, n = 2) MBAHA model is presented for determining concentration. Dose calculations are not presented, but note that dose can be estimated by numerical integration of concentration or analytically by using the definition of dose in EQ-5 applied with matrix calculations to EQ-15.

Two-Volume Material Balance Model

For a 2 × 2 matrix A, $0=det\left(A-\lambda I\right)$ can be solved for eigenvalues and eigenvectors analytically:

${\lambda }^{\left(1,2\right)}=\frac{a_{11}+a_{22}}{2}\pm \frac{1}{2}\sqrt{{\left(a_{11}+a_{22}\right)}^2-4\left(a_{11}{a}_{22}-a_{12}{a}_{21}\right)}$

EQ-16

To determine the eigenvectors,

$\vec{0}=\left(A-\lambda I\right)\vec{k}=\left[\begin{array}{c}\left(a_{11}-\lambda \right)k_1+a_{12}{k}_2\\ a_{21}{k}_1+\left(a_{22}-\lambda \right)k_2\end{array}\right]$

The M matrix of eigenvectors and its inverse are

$M=\left[\begin{array}{cc}{a}_{12}& {\lambda }^{\left(2\right)}-a_{22}\\ {\lambda }^{\left(1\right)}-a_{11}& a_{21}\end{array}\right];$

$M^{-1}=\frac{1}{a_{12}{a}_{21}-\left({\lambda }^{\left(1\right)}-a_{11}\right)\left({\lambda }^{\left(2\right)}-a_{22}\right)}\left[\begin{array}{cc}{a}_{21}& -{\lambda }^{\left(2\right)}+a_{22}\\ -{\lambda }^{\left(1\right)}+a_{11}& a_{12}\end{array}\right]$

There are special cases to maintain numerical stability. The case where when no mixing occurs across volumes (a ₁₂ = a ₂₁ = 0) is not addressed, because this is simply 2 one-volume models. Three cases are addressed here: (1) $a_{12}\ne 0, a_{21}\ne 0;$ (2) $a_{12}=0, a_{21}\ne 0$ ; and (3) $a_{12}=0,a_{21}\ne 0, a=a_{11}=a_{22}$ .

1) For $a_{12}\ne 0,\hspace{0.17em}\hspace{0.17em}{a}_{21}\ne 0$ :

$M{e}^{\Lambda t}{M}^{-1}=\left[\begin{array}{cc}{a}_{12}{e}^{{\lambda }^{\left(1\right)}t}& \left({\lambda }^{\left(2\right)}-a_{22}\right)e^{{\lambda }^{\left(2\right)}t}\\ \left({\lambda }^{\left(1\right)}-a_{11}\right)e^{{\lambda }^{\left(1\right)}t}& a_{21}{e}^{{\lambda }^{\left(2\right)}t}\end{array}\right]M^{-1}$

$=\frac{\left[\begin{array}{cc}{a}_{12}{a}_{21}{e}^{{\lambda }^{\left(1\right)}t}+\left(-{\lambda }^{\left(1\right)}+a_{11}\right)\left({\lambda }^{\left(2\right)}-a_{22}\right)e^{{\lambda }^{\left(2\right)}t}& a_{12}\left(-{\lambda }^{\left(2\right)}+a_{22}\right)e^{{\lambda }^{\left(1\right)}t}+a_{12}\left({\lambda }^{\left(2\right)}-a_{22}\right)e^{{\lambda }^{\left(2\right)}t}\\ \left({\lambda }^{\left(1\right)}-a_{11}\right)a_{21}{e}^{{\lambda }^{\left(1\right)}t}+\left(-{\lambda }^{\left(1\right)}+a_{11}\right)a_{21}{e}^{{\lambda }^{\left(2\right)}t}& \left({\lambda }^{\left(1\right)}-a_{11}\right)\left(-{\lambda }^{\left(2\right)}+a_{22}\right)e^{{\lambda }^{\left(1\right)}t}+a_{12}{a}_{21}{e}^{{\lambda }^{\left(2\right)}t}\end{array}\right]}{a_{12}{a}_{21}-\left({\lambda }^{\left(1\right)}-a_{11}\right)\left({\lambda }^{\left(2\right)}-a_{22}\right)}$

$\begin{array}{l}M{e}^{\Lambda t}{M}^{-1}\vec{s}\left(t\right)\\ =\frac{\left[\begin{array}{c}{a}_{12}{a}_{21}{e}^{{\lambda }^{\left(1\right)}t}{s}_1\left(t\right)+\left(-{\lambda }^{\left(1\right)}+a_{11}\right)\left({\lambda }^{\left(2\right)}-a_{22}\right)e^{{\lambda }^{\left(2\right)}t}{s}_1\left(t\right)+a_{12}\left(-{\lambda }^{\left(2\right)}+a_{22}\right)e^{{\lambda }^{\left(1\right)}t}{s}_2\left(t\right)+a_{12}\left({\lambda }^{\left(2\right)}-a_{22}\right)e^{{\lambda }^{\left(2\right)}t}{s}_2\left(t\right)\\ \left({\lambda }^{\left(1\right)}-a_{11}\right)a_{21}{e}^{{\lambda }^{\left(1\right)}t}{s}_1\left(t\right)+\left(-{\lambda }^{\left(1\right)}+a_{11}\right)a_{21}{e}^{{\lambda }^{\left(2\right)}t}{s}_1\left(t\right)+\left({\lambda }^{\left(1\right)}-a_{11}\right)\left(-{\lambda }^{\left(2\right)}+a_{22}\right)e^{{\lambda }^{\left(1\right)}t}{s}_2\left(t\right)+a_{12}{a}_{21}{e}^{{\lambda }^{\left(2\right)}t}{s}_2\left(t\right)\end{array}\right]}{a_{12}{a}_{21}-\left({\lambda }^{\left(1\right)}-a_{11}\right)\left({\lambda }^{\left(2\right)}-a_{22}\right).}\end{array}$

With a delta source in volume 1, ${\vec{s}}^{{\delta }_1}\left(t\right)={\left[\begin{array}{cc}{s}_1\delta \left(t-t_1\right)& 0\end{array}\right]}^T$ , $\vec{\sigma }\left(t\right)$ is 0 for t < t₁ and for t₁ ≤ t:

${\vec{\sigma }}^{{\delta }_1}\left(t\right)=\underset{0}{\overset{t}{\int }}\hspace{0.17em}M{e}^{\Lambda \left(t-\tau \right)}{M}^{-1}{\vec{s}}^{{\delta }_1}\left(\tau \right)d\tau =\frac{\left[\begin{array}{c}{a}_{12}{a}_{21}{e}^{{\lambda }^{\left(1\right)}\left(t-t_1\right)}{s}_1+\left(-{\lambda }^{\left(1\right)}+a_{11}\right)\left({\lambda }^{\left(2\right)}-a_{22}\right)e^{{\lambda }^{\left(2\right)}\left(t-t_1\right)}{s}_1\\ \left({\lambda }^{\left(1\right)}-a_{11}\right)a_{21}{e}^{{\lambda }^{\left(1\right)}\left(t-t_1\right)}{s}_1+\left(-{\lambda }^{\left(1\right)}+a_{11}\right)a_{21}{e}^{{\lambda }^{\left(2\right)}\left(t-t_1\right)}{s}_1\end{array}\right]}{a_{12}{a}_{21}-\left({\lambda }^{\left(1\right)}-a_{11}\right)\left({\lambda }^{\left(2\right)}-a_{22}\right)}$

For ${\vec{s}}^{{\delta }_2}\left(t\right)={\left[\begin{array}{cc}0& s_2\delta \left(t-t_2\right)\end{array}\right]}^T$ , $\vec{\sigma }\left(t\right)$ is 0 for t < t₂ and for t₂ ≤ t it can be shown that

${\vec{\sigma }}^{{\delta }_2}\left(t\right)=\underset{0}{\overset{t}{\int }}\hspace{0.17em}M{e}^{\Lambda \left(t-\tau \right)}{M}^{-1}{\vec{s}}^{{\delta }_2}\left(\tau \right)d\tau =\frac{\left[\begin{array}{c}{a}_{12}\left(-{\lambda }^{\left(2\right)}+a_{22}\right)e^{{\lambda }^{\left(1\right)}\left(t-t_2\right)}{s}_2+a_{12}\left({\lambda }^{\left(2\right)}-a_{22}\right)e^{{\lambda }^{\left(2\right)}\left(t-t_2\right)}{s}_2\\ \left({\lambda }^{\left(1\right)}-a_{11}\right)\left(-{\lambda }^{\left(2\right)}+a_{22}\right)e^{{\lambda }^{\left(1\right)}\left(t-t_2\right)}{s}_2+a_{12}{a}_{21}{e}^{{\lambda }^{\left(2\right)}\left(t-t_2\right)}{s}_2\end{array}\right]}{a_{12}{a}_{21}-\left({\lambda }^{\left(1\right)}-a_{11}\right)\left({\lambda }^{\left(2\right)}-a_{22}\right)}$

With a delta function at 0 =t₁ = t₂, s₁ = c₁(0), s₂ = c₂(0), and EQ-15:

${\vec{\sigma }}^{{\delta }_1}\left(t\right)+{\vec{\sigma }}^{{\delta }_2}\left(t\right)=\underset{0}{\overset{t}{\int }}\hspace{0.17em}M{e}^{\Lambda \left(t-\tau \right)}{M}^{-1}\left({\vec{s}}^{{\delta }_1}\left(\tau \right)+{\vec{s}}^{{\delta }_2}\left(\tau \right)\right) d\tau =M{e}^{\Lambda t}{M}^{-1}\vec{c}\left(0\right)=\vec{b}\left(t\right)$

which is the same as setting initial concentrations. Focus is then placed on $\vec{\sigma }\left(t\right)$for solution.

A similar approach can be utilized for a rectangular source in volume 1 that has

${\vec{s}}^{\left[1\right]}\left(t\right)=\{\begin{array}{cc}{\left[\begin{array}{cc}0& 0\end{array}\right]}^T& ,t\le t_1\\ {\left[\begin{array}{cc}{s}_1& 0\end{array}\right]}^T& ,t_1\le t\le t_2\\ {\left[\begin{array}{cc}0& 0\end{array}\right]}^T& \text{, }{t}_2<t\end{array}$

And with $\underset{t_1}{\overset{t_2}{\int }}\hspace{0.17em}{e}^{\lambda \left(t-\tau \right)}d\tau =\left(e^{\lambda \left(t-t_1\right)}-e^{\lambda \left(t-t_2\right)}\right)/\lambda$ , ${\vec{\sigma }}^{\left[1\right]}\left(t\right)=\underset{0}{\overset{t}{\int }}\hspace{0.17em}M{e}^{\Lambda \left(t-\tau \right)}{M}^{-1}{\vec{s}}^{\left[1\right]}\left(\tau \right)d\tau$ is

$\begin{array}{c}\begin{array}{cc}\left[\begin{array}{c}0\\ 0\end{array}\right]& ,t\le t_1,\end{array}\\ \begin{array}{cc}\frac{\left[\begin{array}{c}{a}_{12}{a}_{21}\left(e^{{\lambda }^{\left(1\right)}\left(t-t_1\right)}-1\right)\frac{s_1}{{\lambda }^{\left(1\right)}}+\left(-{\lambda }^{\left(1\right)}+a_{11}\right)\left({\lambda }^{\left(2\right)}-a_{22}\right)\left(e^{{\lambda }^{\left(2\right)}\left(t-t_1\right)}-1\right)\frac{s_1}{{\lambda }^{\left(2\right)}}\\ \left({\lambda }^{\left(1\right)}-a_{11}\right)a_{21}\left(e^{{\lambda }^{\left(1\right)}\left(t-t_1\right)}-1\right)\frac{s_1}{{\lambda }^{\left(1\right)}}+\left(-{\lambda }^{\left(1\right)}+a_{11}\right)a_{21}\left(e^{{\lambda }^{\left(2\right)}\left(t-t_1\right)}-1\right)\frac{s_1}{{\lambda }^{\left(2\right)}}\end{array}\right]}{a_{12}{a}_{21}-\left({\lambda }^{\left(1\right)}-a_{11}\right)\left({\lambda }^{\left(2\right)}-a_{22}\right)}& ,t_1\le t\le t_2,\hspace{0.17em}\text{and}\end{array}\\ \begin{array}{cc}\frac{\left[\begin{array}{c}{a}_{12}{a}_{21}\left(e^{{\lambda }^{\left(1\right)}\left(t-t_1\right)}-e^{{\lambda }^{\left(1\right)}\left(t-t_2\right)}\right)\frac{s_1}{{\lambda }^{\left(1\right)}}+\left(-{\lambda }^{\left(1\right)}+a_{11}\right)\left({\lambda }^{\left(2\right)}-a_{22}\right)\left(e^{{\lambda }^{\left(2\right)}\left(t-t_1\right)}-e^{{\lambda }^{\left(2\right)}\left(t-t_2\right)}\right)\frac{s_1}{{\lambda }^{\left(2\right)}}\\ \left({\lambda }^{\left(1\right)}-a_{11}\right)a_{21}\left(e^{{\lambda }^{\left(1\right)}\left(t-t_1\right)}-e^{{\lambda }^{\left(1\right)}\left(t-t_2\right)}\right)\frac{s_1}{{\lambda }^{\left(1\right)}}+\left(-{\lambda }^{\left(1\right)}+a_{11}\right)a_{21}\left(e^{{\lambda }^{\left(2\right)}\left(t-t_1\right)}-e^{{\lambda }^{\left(2\right)}\left(t-t_2\right)}\right)\frac{s_1}{{\lambda }^{\left(2\right)}}\end{array}\right]}{a_{12}{a}_{21}-\left({\lambda }^{\left(1\right)}-a_{11}\right)\left({\lambda }^{\left(2\right)}-a_{22}\right)}& \text{, }{t}_2<t.\end{array}\end{array}$

For a rectangular source s ₂ during $t_1\le t\le t_2$ in volume 2, ${\vec{\sigma }}^{\left[2\right]}\left(t\right)$ for t₂ ≤ t is

$\frac{\left[\begin{array}{c}{a}_{12}\left(-{\lambda }^{\left(2\right)}+a_{22}\right)\left(e^{{\lambda }^{\left(1\right)}\left(t-t_1\right)}-e^{{\lambda }^{\left(1\right)}\left(t-t_2\right)}\right)\frac{s_2}{{\lambda }^{\left(1\right)}}+a_{12}\left({\lambda }^{\left(2\right)}-a_{22}\right)\left(e^{{\lambda }^{\left(2\right)}\left(t-t_1\right)}-e^{{\lambda }^{\left(2\right)}\left(t-t_2\right)}\right)\frac{s_2}{{\lambda }^{\left(2\right)}}\\ \left({\lambda }^{\left(1\right)}-a_{11}\right)\left(-{\lambda }^{\left(2\right)}+a_{22}\right)\left(e^{{\lambda }^{\left(1\right)}\left(t-t_1\right)}-e^{{\lambda }^{\left(1\right)}\left(t-t_2\right)}\right)\frac{s_2}{{\lambda }^{\left(1\right)}}+a_{12}{a}_{21}\left(e^{{\lambda }^{\left(2\right)}\left(t-t_1\right)}-e^{{\lambda }^{\left(2\right)}\left(t-t_2\right)}\right)\frac{s_2}{{\lambda }^{\left(2\right)}}\end{array}\right]}{a_{12}{a}_{21}-\left({\lambda }^{\left(1\right)}-a_{11}\right)\left({\lambda }^{\left(2\right)}-a_{22}\right)}\hspace{0.17em}\hspace{0.17em}\text{.}\hspace{0.17em}$

EQ-17

With ${\vec{\sigma }}^{\left[2\right]}\left(t\right)=0$ for t ≤ t₁, and the t₁ ≤ t ≤ t₂ result obtained by setting t₂ = t in the equation above. The equation from the rectangular functions during t₁ ≤ t ≤ t₂ can also be used to represent a constantly emitting source.

2) For $a_{12}=0,\hspace{0.17em}\hspace{0.17em}{a}_{21}\ne 0$ :

${\lambda }^{\left(1,2\right)}=\frac{a_{11}+a_{22}}{2}\pm \frac{a_{11}-a_{22}}{2}=\{\begin{array}{c}{a}_{11}\\ a_{22}\end{array}$

EQ-18

To determine the eigenvectors

$\vec{0}=\left(A-\lambda I\right)\vec{k}=\left[\begin{array}{c}\left(a_{11}-\lambda \right)k_1\\ a_{21}{k}_1+\left(a_{22}-\lambda \right)k_2\end{array}\right]=\{\begin{array}{c}\left[\begin{array}{c}0\\ a_{21}{k}_1+\left(a_{22}-a_{11}\right)k_2\end{array}\right],{\lambda }^{\left(1\right)}\\ \left[\begin{array}{c}\left(a_{11}-a_{22}\right)k_1\\ a_{21}{k}_1\end{array}\right],{\lambda }^{\left(2\right)}\end{array}$

The M matrix of eigenvectors and its inverse are

$M=\left[\begin{array}{cc}{a}_{22}-a_{11}& 0\\ -a_{21}& 1\end{array}\right];M^{-1}=\frac{1}{a_{22}-a_{11}}\left[\begin{array}{cc}1& 0\\ a_{21}& a_{22}-a_{11}\end{array}\right]$

$M{e}^{\Lambda t}{M}^{-1}=\left[\begin{array}{cc}\left(a_{22}-a_{11}\right)e^{{\lambda }^{\left(1\right)}t}& 0\\ -a_{21}{e}^{{\lambda }^{\left(1\right)}t}& e^{{\lambda }^{\left(2\right)}t}\end{array}\right]M^{-1}=\left[\begin{array}{cc}{e}^{{\lambda }^{\left(1\right)}t}& 0\\ \frac{-a_{21}{e}^{{\lambda }^{\left(1\right)}t}+a_{21}{e}^{{\lambda }^{\left(2\right)}t}}{a_{22}-a_{11}}& e^{{\lambda }^{\left(2\right)}t}\end{array}\right]$

For the delta function source s: ${\vec{\sigma }}^{{\delta }_1}\left(t\right)$ is 0 for t < t₁ and for t₁ ≤ t and ${\vec{s}}^{{\delta }_2}\left(t\right)$ with t₂:

${\vec{\sigma }}^{{\delta }_1}\left(t\right)=\underset{0}{\overset{t}{\int }}\hspace{0.17em}M{e}^{\Lambda \left(t-\tau \right)}{M}^{-1}{\vec{s}}^{{\delta }_1}\left(\tau \right)d\tau =\left[\begin{array}{c}{s}_1{e}^{{\lambda }^{\left(1\right)}\left(t-t_1\right)}\\ \left(-e^{{\lambda }^{\left(1\right)}\left(t-t_1\right)}+e^{{\lambda }^{\left(2\right)}\left(t-t_1\right)}\right)\frac{a_{21}{s}_1}{a_{22}-a_{11}}\end{array}\right]$

EQ-19

${\vec{\sigma }}^{{\delta }_2}\left(t\right)=\underset{0}{\overset{t}{\int }}M{e}^{\Lambda \left(t-\tau \right)}{M}^{-1}{\vec{s}}^{{\delta }_2}\left(\tau \right)d\tau =\left[\begin{array}{c}0\\ s_2{e}^{{\lambda }^{\left(2\right)}\left(t-t_2\right)}\end{array}\right]$

For rectangular sources, first in volume 1 alone, ${\vec{\sigma }}^{\left[1\right]}\left(t\right)=\underset{0}{\overset{t}{\int }}\hspace{0.17em}M{e}^{\Lambda \left(t-\tau \right)}{M}^{-1}{\vec{s}}^{\left[1\right]}\left(\tau \right)d\tau$ is

$\{\begin{array}{cc}\left[\begin{array}{c}0\\ 0\end{array}\right]& ,t\le t_1,\\ \left[\begin{array}{c}\left(e^{{\lambda }^{\left(1\right)}\left(t-t_1\right)}-1\right)\frac{s_1}{{\lambda }^{\left(1\right)}}\\ \frac{-a_{21}{s}_1\left(e^{{\lambda }^{\left(1\right)}\left(t-t_1\right)}-1\right)}{\left(a_{22}-a_{11}\right){\lambda }^{\left(1\right)}}+\frac{a_{21}{s}_1\left(e^{{\lambda }^{\left(2\right)}\left(t-t_1\right)}-1\right)}{\left(a_{22}-a_{11}\right){\lambda }^{\left(2\right)}}\end{array}\right]& ,t_1\le t\le t_2,\\ \begin{array}{c}\left(e^{{\lambda }^{\left(1\right)}\left(t-t_1\right)}-e^{{\lambda }^{\left(1\right)}\left(t-t_2\right)}\right)\frac{s_1}{{\lambda }^{\left(1\right)}}\\ \frac{-a_{21}{s}_1\left(e^{{\lambda }^{\left(1\right)}\left(t-t_1\right)}-e^{{\lambda }^{\left(1\right)}\left(t-t_2\right)}\right)}{\left(a_{22}-a_{11}\right){\lambda }^{\left(1\right)}}+\frac{a_{21}{s}_1\left(e^{{\lambda }^{\left(2\right)}\left(t-t_1\right)}-e^{{\lambda }^{\left(2\right)}\left(t-t_2\right)}\right)}{\left(a_{22}-a_{11}\right){\lambda }^{\left(2\right)}}\end{array}& \text{, }{t}_2<t.\end{array}$

EQ-20

For a rectangular source s ₂ located in volume 2, for t₂ ≤ t (note that 0 = a₁₂ implies 0 = r₂₁ and volume 2 is isolated from volume 1 since the source is in volume 2):

${\vec{\sigma }}^{\left[2\right]}\left(t\right)=\{\begin{array}{cc}\left[\begin{array}{c}0\\ 0\end{array}\right]& ,t\le t_1\\ \left[\begin{array}{c}0\\ \left(e^{{\lambda }^{\left(2\right)}\left(t-t_1\right)}-1\right)\frac{s_2}{{\lambda }^{\left(2\right)}}\end{array}\right]& ,t_1\le t\le t_2\\ \left[\begin{array}{c}0\\ \left(e^{{\lambda }^{\left(2\right)}\left(t-t_1\right)}-e^{{\lambda }^{\left(2\right)}\left(t-t_2\right)}\right)\frac{s_2}{{\lambda }^{\left(2\right)}}\end{array}\right]& \text{, }{t}_2<t\end{array}$

EQ-21

3) For $a_{12}=0,\hspace{0.17em}\hspace{0.17em}{a}_{21}\ne 0,\hspace{0.17em}\hspace{0.17em}a=\hspace{0.17em}{a}_{11}=\hspace{0.17em}{a}_{22}$ :

${\lambda }^{\left(1,2\right)}=\frac{a_{11}+a_{22}}{2}=a$

The single eigenvector is

$\vec{0}=\left(A-\lambda I\right)\vec{k}=\left[\begin{array}{c}\left(a-\lambda \right)k_1\\ a_{21}{k}_1+\left(a-\lambda \right)k_2\end{array}\right]=\{\begin{array}{c}\left[\begin{array}{c}0\\ a_{21}{k}_1\end{array}\right],{\lambda }^{\left(1\right)}\\ \left[\begin{array}{c}0\\ a_{21}{k}_1\end{array}\right],{\lambda }^{\left(2\right)}\end{array},\vec{k}=\left[\begin{array}{c}0\\ 1\end{array}\right]$

For this matrix A, the solution to $\overrightarrow{ċ}\left(t\right)=A\vec{c}\left(t\right)$ is given as follows (β₁ and β₂ constants from initial conditions):

$\vec{c}\left(t\right)={\beta }_1{e}^{at}\vec{k}+{\beta }_2\left(t{e}^{at}\vec{k}+e^{at}\vec{u}\right)$

where $\vec{u}$ is a solution to $\left(A-\lambda I\right)\vec{u}=\vec{k}$ . For this example,

$\vec{k}=\left[\begin{array}{c}0\\ 1\end{array}\right]=\left(A-\lambda I\right)\vec{u}=\left[\begin{array}{c}0\\ a_{21}{u}_1\end{array}\right]\therefore \vec{u}=\left[\begin{array}{c}\frac{1}{a_{21}}\\ 0\end{array}\right]$

$\vec{c}\left(t\right)=c_2\left(0\right)e^{at}\left[\begin{array}{c}0\\ 1\end{array}\right]+a_{21}{c}_1\left(0\right)\left(t{e}^{at}\left[\begin{array}{c}0\\ 1\end{array}\right]+e^{at}\left[\begin{array}{c}\frac{1}{a_{21}}\\ 0\end{array}\right]\right)$

$\vec{c}\left(t\right)=e^{at}\left(\left[\begin{array}{c}{c}_1\left(0\right)\\ c_2\left(0\right)\end{array}\right]+t\left[\begin{array}{c}0\\ a_{21}{c}_1\left(0\right)\end{array}\right]\right)$

Model with No Particles Flowing from Soft-Walled Enclosure Into ABSL-4 Room

The 2 × 2 model provides an analytic framework that represents 2 volumes (eg, BSL-4 room and soft-walled enclosure), removal of particles at a specified rate for each volume (eg, HEPA filtration), arbitrary initial concentrations in each volume, instantaneous sources of particles, and continuous sources of particles with specified start and stop spray duration. Since the math is available in close form, it can be used for additional analyses, including dose estimation. The model consistency is demonstrated with the experimental measurements.

For the 2-volume, n = 2, experiment with a soft-walled enclosure in a room, use EQ-13 and EQ-16:

$A=\left[\begin{array}{cc}\frac{-r_{12}-r_{10}}{v_1}& \frac{r_{21}}{v_1}\\ \frac{r_{12}}{v_2}& \frac{-r_{21}-r_{20}}{v_2}\end{array}\right]$

${\lambda }^{\left(1,2\right)}=-\frac{1}{2}\left(\frac{r_{12}+r_{10}}{v_1}+\frac{r_{21}+r_{20}}{v_2}\right)\pm \frac{1}{2}\sqrt{{\left(\frac{-r_{12}-r_{10}}{v_1}+\frac{r_{21}+r_{20}}{v_2}\right)}^2+4\frac{r_{12}}{v_1}\frac{r_{21}}{v_2}}$

based on the experimental design and the observation that particles do not appear to flow from the enclosure into the room. The variables in the model are shown in Table 3.

Table 3.

Variables for the Model When Particles Do Not Flow from the Soft-Walled Enclosure Into the Room.

Variable	Value (Units)	Description
v 1	67.32 (m3)	Volume 1: the ABSL-4 room (74.98 [m3]) – enclosure (7.66 [m3])
v 2	7.66 (m3)	Volume 2: the soft-walled primary containment enclosure
r 10	31.2 (m3/min)	Rate of flow leaving ABSL-4 room (ie, HEPA filtered)
r 20	8.5 (m3/min)	Rate of flow leaving the soft-walled enclosure (ie, HEPA filtered)
r 12	8.5 (m3/min)	Rate of flow from volume 1 to 2 (ie, pulled into the soft-walled enclosure)
r 21	0 (m3/min)	Rate of flow from volume 2 to 1 (ie, leaking from the soft-walled enclosure). When r21 = 0, then a12 = 0.

Abbreviations: ABSL, animal biosafety level; HEPA, high-efficiency particulate air.

The resulting matrix A and eigenvalues (EQ-16) are

$A=\left[\begin{array}{cc}\frac{-r_{12}-r_{10}}{v_1}& 0\\ \frac{r_{12}}{v_2}& \frac{-r_{20}}{v_2}\end{array}\right]=\left[\begin{array}{cc}-0.5897& 0\\ 1.11& -1.1097\end{array}\right]$

${\lambda }^{\left(1\right)}=a_{11}=\frac{-r_{12}-r_{10}}{v_1}=-0.5897;{\lambda }^{\left(2\right)}=a_{22}=\frac{-r_{20}}{v_2}=-1.1097$

As shown in EQ-18 and EQ-19, a source $\vec{s}={\left[\begin{array}{cc}{s}_1\delta \left(t\right)& s_2\end{array}\delta \left(t\right)\right]}^T$ has concentrations

$\left[\begin{array}{c}{s}_1{e}^{\left(\frac{-r_{12}-r_{10}}{v_1}\right)t}\\ s_1\frac{\frac{r_{12}}{v_2}}{\frac{-r_{20}}{v_2}+\frac{r_{12}+r_{10}}{v_1}}\left(-e^{\left(\frac{-r_{12}-r_{10}}{v_1}\right)t}+e^{\left(\frac{-r_{20}}{v_2}\right)t}\right)+s_2{e}^{\left(\frac{-r_{20}}{v_2}\right)t}\end{array}\right]$

verifying the rate in volume 1 that was approximated in EQ-10. The concentration and rate for volume 2 are a little more complicated with the difference of 2 exponentials. Since ${\lambda }^{\left(1\right)}>{\lambda }^{\left(2\right)}$ , the exponential term with λ⁽¹⁾ will dominate as t increases and the soft-walled enclosure is removing particles at approximately the same rate as the ABSL-4 room, which is why the decay appears the same for the first 30 to 60 minutes as the room recovers from hallway particles as seen in the experimental data plots by Landon et al. ¹ Encoding the relevant equations for delta and rectangular sources in Visual Basic in Excel and using parameters consistent with the experiments, a model calculation for the experimental data is presented in Figure 5.

Figure 5.

Two-volume model with no particle flow from the soft-walled enclosure into the room, background of 1000 particles/m³, delta function sources at t = 0 of strength 3 000 000 in both volumes, and rectangular sources of strengths 200 000 000 and 400 000 000 at times 70-71 and 190-192, respectively, released in the enclosure. The model results compare favorably with the experimental data presented in Landon et al. ¹ ABSL-4, animal biosafety level 4.

In a similar manner, Figure 6 presents model results of the concentration in the ABSL-4 room for (1) a 1-minute source in the ABSL-4 room with no soft-walled enclosure operating and (2) the same source inside the operating soft-walled enclosure. Also of significance in this model is the faster return to background concentration from the introduction of hallway particles for the ABSL-4 room when the enclosure HEPA filter is operating versus when the soft-walled enclosure is absent.

Figure 6.

Modeled aerosol concentration in the animal biosafety level 4 (ABSL-4) room with “hallway” contamination for <30 minutes and a source from 45 to 46 minutes both with no soft-walled enclosure and with the source inside an operating enclosure. The shape and rates compare favorably with experimental data in Landon et al. ¹

Model with Particles Flowing from the Soft-Walled Enclosure Into ABSL-4 Room

To investigate the potential impact of particles escaping from the soft-walled enclosure into the ABSL-4 room, 3 sources were modeled: (1) an initial concentration (delta source at time t = 0), (2) a 1-minute rectangular source inside the soft-walled enclosure, (3) a 1-minute source outside the soft-walled enclosure and a delta source outside the soft-walled enclosure. The delta source at t = 0 and 1-minute sources are scaled to match the experimental data presented by Landon et al ¹ for the hallway contamination and first challenge. Two different rates of particles flowing from the soft-walled enclosure into the room (ie, r ₂₁) are used: one set to produce a noticeable variation above background but within the maximum variation (r ₂₁ = 0.00085, approximately 4000 particles/m³ above background) observed experimentally and with a rate equal to one-tenth of the soft-walled enclosure HEPA flow rate returned to the ABSL-4 room (r ₂₁ = 0.85, which is also 1000 times larger than the other model calculation). The results of the model calculation are presented in Figure 7. To help explain this, EQ-12 is examined for a rectangular source in volume 2, and notice that the resulting concentration in volume 1 is weighted by $a_{12}=r_{21}/v_1$ . So as r₂₁ was changed, c ₁ changed roughly proportionally.

Figure 7.

Model calculation for non-zero particle escape (r ₂₁) from the soft-walled enclosure into the animal biosafety level 4 (ABSL-4) room for 2 different values of r ₂₁. A 10-minute source is used in the soft-walled enclosure (60-70 minutes) to help display the different responses with the peak “leakage” concentration approximately scaling with r ₂₁. For a 1-minute source in the ABSL-4 room, there are no significant differences in performance.

Conclusions

New experimental data on primary containment enclosures in ABSL-4 and data available from the published work of other groups were utilized to demonstrate a general material balance biocontainment model. A biocontainment room by itself (n = 1) and a biocontainment room holding a primary containment enclosure (n = 2) can both be represented by the MBAHA model to calculate concentration and integrated exposure dose. The analytic 2-dimensional model applies to all 2-volume interactions, including biosafety cabinets within laboratory rooms. The general model can be applied for an arbitrary number of volumes, n, by numerical methods on a computer that can solve matrix equations with eigenvalue decompositions.

By using simple assumptions that are amenable to spreadsheet calculation and demonstrating correlation with experimental data, a general model has been provided for biocontainment professionals to estimate particle concentration, dose, and exposure time to reach a specified dose for a wide range of scenarios.