- Research
- Open access
- Published:

# Advanced computational modeling for *in vitro* nanomaterial dosimetry

*Particle and Fibre Toxicology*
**volume 12**, Article number: 32 (2015)

## Abstract

### Background

Accurate and meaningful dose metrics are a basic requirement for *in vitro* screening to assess potential health risks of engineered nanomaterials (ENMs). Correctly and consistently quantifying what cells “see,” during an *in vitro* exposure requires standardized preparation of stable ENM suspensions, accurate characterizatoin of agglomerate sizes and effective densities, and predictive modeling of mass transport. Earlier transport models provided a marked improvement over administered concentration or total mass, but included assumptions that could produce sizable inaccuracies, most notably that all particles at the bottom of the well are adsorbed or taken up by cells, which would drive transport downward, resulting in overestimation of deposition.

### Methods

Here we present development, validation and results of two robust computational transport models. Both three-dimensional computational fluid dynamics (CFD) and a newly-developed one-dimensional Distorted Grid (DG) model were used to estimate delivered dose metrics for industry-relevant metal oxide ENMs suspended in culture media. Both models allow simultaneous modeling of full size distributions for polydisperse ENM suspensions, and provide deposition metrics as well as concentration metrics over the extent of the well. The DG model also emulates the biokinetics at the particle-cell interface using a Langmuir isotherm, governed by a user-defined dissociation constant, *K*
_{D}, and allows modeling of ENM dissolution over time.

### Results

Dose metrics predicted by the two models were in remarkably close agreement. The DG model was also validated by quantitative analysis of flash-frozen, cryosectioned columns of ENM suspensions. Results of simulations based on agglomerate size distributions differed substantially from those obtained using mean sizes. The effect of cellular adsorption on delivered dose was negligible for *K*
_{D} values consistent with non-specific binding (> 1 nM), whereas smaller values (≤ 1 nM) typical of specific high-affinity binding resulted in faster and eventual complete deposition of material.

### Conclusions

The advanced models presented provide practical and robust tools for obtaining accurate dose metrics and concentration profiles across the well, for high-throughput screening of ENMs. The DG model allows rapid modeling that accommodates polydispersity, dissolution, and adsorption. Result of adsorption studies suggest that a reflective lower boundary condition is appropriate for modeling most *in vitro* ENM exposures.

## Background

Assessing the potential health risks associated with exposures to the vast number and variety of engineered nanomaterials (ENMs) entering manufacturing workplaces and now present in myriad consumer products is a daunting task that requires fast, inexpensive and reliable screening strategies [1–5]. Although animal testing may provide the most relevant and reliable assessment of a nanomaterial’s biological activity, the associated cost and throughput constraints, as well as humane concerns, render this approach impractical [1, 4, 6]. A screening approach using *in vitro* cell-based assays provides a logical and efficient alternative [7]. However, *in vitro* testing has often produced results that are inconsistent with those of corresponding *in vivo* studies and even of other *in vitro* studies of the same ENMs [4–8]. Although the many differences between monocultured cell and whole animal experimental systems may account for much of this disparity, it is likely that failure to match *in vitro* and *in vivo* doses, as a result of inadequate characterization of ENM powders and suspensions, and more importantly, failure to account for transport of ENM particles or agglomerates in suspension over time, is in part responsible for these disparities [9–13].

Until recently most *in vitro* studies have reported dose in terms of either an initial administered mass concentration or of a total administered mass [11, 14, 15]. The former assumes that sedimentation either does not occur, or is negligible, and the latter assumes that it is complete, with all of the suspended material instantly transported to the cells at the bottom of the cell culture well. The reality lies between these extremes, and depends upon the intrinsic physicochemical properties of the suspended material, the extrinsic properties of suspending media, and the time course of the exposure.

Of late, greater emphasis has been placed upon achieving a better understanding of the exposures experienced by cells *in vitro*, and of the measures necessary for accurate and meaningful dosimetry [12, 16–20] (+1 for > =20). These include several standardized methodologies, from generation of stable suspensions, to accurate physical characterization of formed agglomerates in suspension, to appropriate modeling of the transport of particles or agglomerates during exposure. An integrated approach for *in vitro* nanodosimetry based on these methodologies is depicted in Fig. 1.

Most industrially-relevant ENMs form agglomerates when dispersed in liquid. The size, or size distribution, and effective density of these agglomerates determine their transport properties (sedimentation and diffusion coefficients), and thus control the rate at which the material settles, which in turn determines the ENM mass, surface and particle number delivered to cells as a function of time [12, 16, 17]. ENM agglomeration state can also influence the nature and extent of biological effects [21, 22]. As we have recently reported, the agglomerate size distribution, and more importantly the agglomerate stability (re-agglomeration rate), are dependent upon the dispersion protocol and properties of the dispersant [11, 17]. Initial sonication of ENM powder in deionized water at or above an ENM-dependent critical delivered sonication energy (*DSE*
_{cr}), followed by dilution into culture media, produces suspensions with remarkably stable size distributions, whereas the size distributions of suspensions created at less than *DSE*
_{cr} can change substantially (suggesting re-agglomeration) over time [11].

Once a stable suspension is created, particle sizes are typically determined by either dynamic light scattering (DLS), disk centrifugation, or tunable resistive pulse sensing technology (TRPS) [13] and are routinely reported in *in vitro* studies. Effective density of agglomerates can differ from bulk ENM density by many fold due to the formation of protein coronae and intra-agglomerate trapping of media [12, 17]. Despite its important influence on transport and thereby on delivered cellular dose, it is rarely measured and reported in published *in vitro* studies. This is primarily due to methodological limitations that until recently required expensive instrumentation such as an Analytical Ultracentrifuge (AUC). The volumetric centrifugation method (VCM), recently developed by the authors, provides a fast and practical alternative for determination of effective density of ENMs in suspension [12].

The first computational fate and transport model developed for ENM dosimetry was the In vitro Sedimentation, Diffusion and Dosimetry model (ISDD) reported by Hinderliter et al. [16]. The ISDD model utilizes an analytical solution of the Mason Weaver equation to determine, in a suspension column of a given height, and a given concentration of particles or agglomerates of a given diameter and density, the per area mass, surface area, and number of particles, as well as the fraction of total suspended material deposited as a function of time. This provided a ground-breaking improvement in dosimetry accuracy and enabled meaningful comparisons and safety hazard rankings among ENMs of different physiochemical properties. The ISDD initially estimated effective density using the Sterling equation, based on a theoretical fractal model and estimated fractal dimensions for ENM agglomerates [23]. More recently the ISDD has been modified to utilize the empirically-based and more accurate VCM to measure effective density, with the integrated approach now referred to as VCM-ISDD [17].

Although the VCM-ISDD model importantly allows assessment of relative deposition over time, it employs a perfectly adsorptive (sticky) lower boundary condition, the result of which is the prediction of complete deposition for all materials, within a few hours to several days, depending on particle size and density. This is inconsistent with our observations of nanoparticle suspensions over extended periods of time (months). With the exception of cases in which particles are sufficiently large or dense that the role of diffusion relative to sedimentation is negligible, or cases in which strong adsorption or uptake by cells creates a perpetual sink at the bottom of the well, the concentration gradient created by sedimentation would drive net diffusion upward, such that an equilibrium is eventually obtained, and complete deposition would not occur. The VCM-ISDD model is also limited to modeling transport of a single agglomerate size at a time. Although transport of a polydisperse suspension can be modeled by summing results of individual simulations for each size category, this is unwieldy for high-resolution size distributions, which from DLS, for example, may contain up to 300 size species, and moreover it precludes the possibility of adjusting for dynamic cross-species interactions and effects on transport parameters. In addition, output from the VCM-ISDD is limited to the absolute or fractional deposition over time, and does not provide concentration metrics, which more closely represent what cells “see.” Finally, the VCM-ISDD does not provide mechanisms for modeling of ENMs that undergo dissolution with time, which can not only reduce deposition of solid ENM and add a dissolved dose component, but may also dynamically change relevant transport parameters of particles. Thus although the VCM-ISDD provides useful measures and relative comparisons of transport and effective dose among different materials, a number of refinements are needed.

Here we present results from two advanced numerical transport models that address these issues: 1) a new one-dimensional iterative distorted grid (DG) model, and 2) the gold standard three-dimensional computational fluid dynamics (CFD) approach. Each of these models provide deposition as well as concentration of ENMs as a function of time, both at the bottom of a culture well (the cell microenvironment) and as a function of vertical position within suspension column. Both models accommodate simultaneous simulation of polydisperse suspensions. The DG model further provides modeling of particles that undergo dissolution over time, as well as a variable ‘stickiness’ boundary condition at the bottom, implemented using a Langmuir isotherm process. As we demonstrate below, the results of the distorted grid model, which was implemented in MATLAB, and typically runs in a few minutes on standard desktop computer, are in close agreement with those from CFD. Finally, we validated the DG model by experimentally measuring particle deposition profiles in thinly-sectioned frozen columns of industry-relevant ENM suspensions.

## Results

### Comparison of Distorted Grid and CFD model simulations

The Distorted Grid (DG) model is based on a model previously developed for analysis of protein systems via their behavior in an ultracentrifuge [24–30], which was adapted for particle transport and implemented in MATLAB. In this model the column of suspended particles is divided into *n* compartments separated by *n* + 1 boundaries, and successive brief rounds of simulated sedimentation and diffusion are performed, with transfer of suspended material between compartments. A detailed mathematical description of the DG model design is provided in the methods section.

In the Computational Fluid Dynamics (CFD) model, particles are initially assigned to compartments of a 3-dimensional grid representation of the suspension column, and a solution of the Navier–Stokes equation is used to calculate the movement of individual particles between compartments. A complete description of the CFD model used in these studies is presented in methods.

Concentration profiles (ENM mass concentrations as a function of distance from the bottom of the well) at 12 h, and fractional deposition (fraction of mass present in the bottom 10 μm of the well) as a function of time, as computed by the DG and CFD models, for CeO_{2}(in RPMI + 0.5 % BSA) and SiO_{2} (in PBS + 0.1 % BSA), are compared in Fig 2a. A thickness of 10 μm was somewhat arbitrarily chosen for the disk at the bottom of the well representing the cell microenvironment. We based this choice on thickness typically observed for typical adherent cells in culture. However both the DG and CFD models allow this thickness value to be specified by the user, and smaller or larger thicknesses may be appropriate depending on the type of cell or system being studied. Properties of ENM powders and of their forms in liquid suspension, including effective densities (*ρ*
_{EV}), are given in Table 1. The volume-weighted size (*d*
_{H}) distributions (calculated from particle number distributions from DLS as described in methods) of liquid suspensions used in these simulations are shown in Additional file 1: Figure S3.

Simulation results for CFD and DG results for CeO_{2} and TiO_{2} are nearly identical (within <5 %) (Fig. 2). Because the CFD model tracks movement of each individual particle it is highly computationally-intensive, which places practical constraints on the number of particles, and thus the model geometry and concentration of ENMs, that can be simulated in a reasonable time with available computational resources. Thus, in our CFD simulations we employed a column height of 1 mm, and concentrations of 0.1 mg ml^{−1} for CeO_{2} and 0.001 mg ml^{−1} for SiO_{2}, resulting in a total of approximately 3.8 × 10^{5} particles in each case. Even at these low values the CFD simulations for CeO_{2} and SiO_{2} required, respectively, 12 and 48 h of compute time on a multi-processor supercomputer. By contrast, the corresponding simulations using the DG model were each completed in less than a minute on a standard single-processor PC.

The results of both the DG model and CFD can be affected by selection of compartment (DG) or grid (CFD) size. To ensure grid size-independence, simulations were performed over a range of compartment sizes for SiO_{2}. Whereas the final predicted deposition or bottom of column concentration increased substantially as compartment and grid sizes were decreased from 50.0 to 5.0 μm, results were independent of compartment/grid sizes for values less than 5 μm (data not shown). Likewise, simulations were run over a range of iteration time intervals to verify that results were independent of time interval between 0.1 to 1.0 seconds.

### Experimental Validation of Distorted Grid model by frozen section analysis

Concentration profiles predicted by the DG model were compared with those obtained from frozen sections of corresponding columns of nanoparticle suspensions. Columns of Fe_{2}O_{3}, CeO_{2} and TiO_{2} suspensions (powder and liquid suspension properties shown in Table 1) were flash frozen and cryosectioned following incubations (to allow transport) at room temperature of 24 h for Fe_{2}O_{3} and CeO_{2}, and 64 h for TiO_{2}. Material concentrations within individual sections were determined by spectrophotometry using standard curves prepared from frozen-and-thawed suspensions of known mass concentrations. The process for creating suspension columns and for obtaining cryosectioned slice concentrations is described in the methods section and illustrated in Additional file 1: Figure S2. For all three ENMs tested the concentration profile predicted by the DG model was in close agreement with that obtained from frozen sections (Fig. 2b).

### DG model-derived mass, particle number and surface dose metrics

Although ENM mass concentration and deposition per unit area are the most commonly-used *in vitro* dose metrics, it is also important to consider exposures in terms of particle number and surface area [9, 10, 15]. Not only is it possible that biological effects may in some cases track more closely with these values, but because of differences in agglomerate formation the relative magnitudes of number and surface area metrics for suspensions of different materials undergoing transport may differ substantially from the corresponding mass metrics. The DG model provides both concentration and deposition outputs in terms of number of particles (N ml^{−1} and N cm^{−2}) and ENM surface area (cm^{2} ml^{−1} and cm^{2} cm^{−2}). Calculation of number and surface area metrics from mass concentration and deposition metrics are described in methods.

DG model-derived mass, particle number and ENM surface area dose metrics as a function of time (from 0 to 120 h) for four ENMs (SiO_{2}, Fe_{2}O_{3}, TiO_{2} and CeO_{2}), with initial mass concentration, *C*
_{0} = 0.1 mg ml^{−1}, suspension column height of 3 mm, and volume averaged particle diameter (〈*d*
_{H}〉_{v} x, Table 1) are shown in Fig. 3. Because we define deposition (not equivalent to binding, which is discussed below) as the total mass, particle number or surface area in the bottom 10 μm of the column, it is calculated as the product of concentration and the volume of that disk. The shapes of the concentration curves are therefore identical to those of the corresponding deposition curves, differing only in scale and units, and both can be represented in the same graph using appropriately scaled axes (Fig. 3). Interestingly, the relative magnitudes and ranking of the four materials by concentration (or deposition) differs for the three different classes of metrics. For example, mass concentration and deposition are at all time points greatest for CeO_{2} and smallest for SiO_{2}, whereas particle number and surface area concentration and deposition of CeO_{2} are much less than those of SiO_{2}. These differences are not a result of transport differences *per se*, but are rather due to differences in agglomeration state. Specifically, since the agglomerates of SiO_{2} (〈*d*
_{H}〉_{v} =149.9 nm) are much smaller than those of CeO_{2} (〈*d*
_{H}〉_{v} = 982.1 nm), the total number and surface area of SiO_{2} for a given mass concentration are considerably greater than the corresponding values for CeO_{2}.

### Effect of particle-cell binding/adsorption – well bottom boundary condition

The data presented to this point, including those from validation experiments, were obtained using a reflective boundary condition at the bottom of the well: particles arriving at the bottom were not bound or removed from the system, but remained present, undergoing Brownian motion, and contributing to the concentration gradient driving diffusion. The opposite extreme boundary condition is one in which all particles that reach the bottom of the well adhere to the cells there, effectively removing them from the concentration gradient. The real world condition likely lies somewhere between these two extremes, and is ENM-, media-, and cell type-dependent. Some of the ENM agglomerates that reach the cell microenvironment in an *in vitro* exposure system may adhere to cells, or in the case of phagocytic cells may be internalized, which in either case would eliminate their contribution to the diffusion gradient opposing sedimentation, and thereby increase the net rate of transport toward the bottom of the well. Quantification of particle-cell binding and uptake is an emerging research area [22] and a complex topic that will be addressed in future work, but because binding can, as we will demonstrate, have a large effect on net transport, it is important to incorporate this variable in fate and transport models. Initially we can employ reasonable estimates of binding, based on our understanding of the particle- and protein-protein interactions that are likely to be involved, which can be later refined as specific data for particle-cell binding kinetics become available. Accordingly, an adsorption strength (stickiness) parameter was incorporated in the DG model to simulate boundary conditions with variable levels of adsorption at the bottom of the well. This is implemented as a Langmuir isotherm process [31], in which adsorption is characterized by the equilibrium dissociation constant, *K*
_{D}, as described in the methods section.

The effect of cellular adsorption on total (free + cell bound) ENM mass concentration in the lower 10 μm exposure volume, as well as the cell-bound ENM mass per unit well bottom area (mg cm^{−2}) for SiO_{2}, Fe_{2}O_{3} , and TiO_{2} (all at *C*
_{0} = 0.1 mg ml^{−1}, for a 3 mm column height, using volume averaged particle diameter (〈*d*
_{H}〉_{v}, Table 1) ), with various values of *K*
_{D}, are shown in Fig. 4. It is clear from these results that very small values of *K*
_{D}, in the nanomolar range, produce sizeable changes in exposure concentrations (at the bottom of the well). Further decreases in *K*
_{D} did not increase predicted well bottom concentrations beyond those for the smallest *K*
_{D} values shown. Such small values correspond to high-affinity binding typical of specific protein interactions, and are several orders of magnitude smaller than the millimolar or micromolar values typical of non-specific interactions [32–34]. Interactions between ENM agglomerates and cell surface biomolecules are most likely of the latter, weak non-specific type. We therefore assume that in most cases the effect of particle-cell binding on transport is negligible, and the reflective boundary condition was therefore used in all other simulations reported here.

Because the VCM-ISDD model imposes a perfectly sticky (sink) lower boundary condition, we would expect it to predict faster and more complete deposition of material than the DG model with a reflective boundary, particularly for smaller particles, where diffusion plays a relatively larger role. To test this we compared dose metrics predicted by the DG and VCM-ISDD models for four materials (SiO_{2}, Fe_{2}O_{3}, TiO_{2} and CeO_{2}), at initial mass concentration, *C*
_{0} = 0.1 mg ml^{−1}, suspension column height = 3 mm, and volume averaged particle diameter (〈*d*
_{H}〉_{v}, Table 1). The results of these simulations are shown in Additional file 1: Figure S4. Whereas both models are in almost perfect agreement regarding the transport of CeO_{2} (〈*d*
_{H}〉_{v} =982.1 nm, *ρ*
_{EV} =1.420 g cm^{−3}), they provide slightly different predictions for TiO_{2}(〈*d*
_{H}〉_{v} =397.8 nm, *ρ*
_{EV} =1.251 g cm^{−3}), and strongly different predictions for Fe_{2}O_{3} (〈*d*
_{H}〉_{v} =234.5 nm, *ρ*
_{EV} =1.335 g cm^{−3}) and SiO_{2} (〈*d*
_{H}〉_{v} =149.9 nm, *ρ*
_{EV} =1.564 g cm^{−3}). Specifically, the VCM-ISDD model predicts deposition fraction approaching 1.0 for all materials, whereas the DG model predicts complete deposition only in the case of CeO_{2}, which forms the largest agglomerates. For the three other materials the DG model predicts that the deposited fraction approaches an equilibrium value well below 1.0, and less than 0.2 for SiO_{2}, which forms the smallest agglomerates.

### Effect of agglomerate size and density on dosimetry

The sedimentation velocity, *v*
_{s}, of a particle of hydrodynamic diameter *d*
_{H} (m) under gravity is defined as

where *g* is the acceleration due to gravity, *g* (m s^{−2}), *ρ*
_{EV} is the effective density of the particle (kg m^{−3}), *ρ*
_{media} is the media density (kg m^{−3}), and *η* is the media dynamic viscosity (kg m^{−1} s^{−1}). The rate of mass transport by diffusion is proportional to the diffusion coefficient *D* (m^{2} s^{−1}), which is defined by the Stokes-Einstein equation as:

in which *k*
_{B} is the Boltzmann constant (kg m^{2} s^{−2} K^{−1}), *T* is the absolute temperature (°K), and *η* is the media dynamic viscosity (kg m^{−1} s^{−1}). Sedimentation velocity is thus proportional to the difference between the densities of the sedimenting particle and the suspending media, and to the square of particle diameter (Eq. 1), whereas diffusion is independent of density and linearly related to particle diameter (Eq. 2). Transport modeling should therefore predict that more dense and particularly larger particles concentrate more rapidly at the bottom of a column of suspension than less dense and smaller particles.

We simulated transport using the DG model for hypothetical particles of constant size (*d*
_{H} = 150 nm) ranging in density from *ρ*
_{EV} =1.1 to 2.5 g cm^{−3}, and for particles of constant density (*ρ*
_{EV} =1.5 g cm^{−3}) ranging in size from *d*
_{H} = 150 to 500 nm. Results of these simulations are shown in Fig. 5. As expected, for particles of constant size, predicted transport was more rapid for particles of higher density, and likewise for particles of constant density, predicted transport proceeded more rapidly for larger particles.

### Effect of Polydispersity

Since most ENM suspensions are polydisperse, it is important to account for the effect of polydispersity when modeling transport. The DG model allows input of a particle size distribution consisting of paired arrays of particle hydrodynamic diameters and fractional contributions to the total agglomerate mass (e.g., from DLS). We compared concentration profiles and dose metrics obtained from the DG model using volume-weighted polydisperse size/fraction arrays with those obtained using volume-weighted average sizes (*d*
_{Hv}), for SiO_{2}, Fe_{2}O_{3}, TiO_{2} and CeO_{2} suspensions (24 h simulations, *C*
_{0} = 0.1 mg ml^{−1}, column height = 3 mm, volume averaged particle diameters shown in Table 1, volume-weighted distributions shown in Additional file 1: Figure S3). The results, shown in Fig. 6, reveal substantial differences between modeling with mean size and modeling using size/fraction arrays. In the case of CeO_{2}, for example, well bottom concentration and fraction deposited reached maximum plateaus at 4 h when transport was modeled using an average diameter, but had only reached approximately two thirds maximum at the same time when modeled using the full size distribution. In the Case of TiO_{2}, average size simulations resulted in higher predicted concentration and deposition at early time points, but lower predicted values at later times, relative to modeling with polydispersity. For both Fe_{2}O_{3} and SiO_{2,} average size simulations substantially under-predicted concentration and deposition at all times relative to the polydisperse simulations.

### Effect of solubilization

For some nanomaterials it is important to account for solubilization of the material over time during *in vitro* exposure. Xia et al. (2008) [35] found that for 10 μg ml^{−1} (122.8 μM) suspensions of ZnO in DMEM with 10 % FBS, up to 74 % of the ZnO was solubilized within 10 minutes (corresponding to a dissolved concentration of 90 μM), and that in water initial dissolution was less extensive (20 μM), but continued linearly over 3 h to a maximum of 60 μM. Solubilization of ENM during an *in vitro* exposure not only reduces the concentration of solid ENM to which cells are exposed, and adds an exposure of solubilized or ionized material, but may also change the transport properties of agglomerate species. Although more complex relationships between solubilization and agglomerate state and agglomerate transport properties are possible, as a first approximation we modeled solubilization as a reduction of agglomerate size in proportion to the extent of dissolution.

Solid ENM mass exposure concentration (in bottom 10 μm) and dissolved ENM concentrations predicted by DG under various solubilization scenarios for ZnO (*C*
_{0} = 0.01 mg ml^{−1}, column height = 3 mm, *d*
_{H} = volume-averaged mean (Table 1)) for 24 h transport duration are shown in Fig. 7. ZnO solid mass exposure concentration (mg ml^{−1} in bottom 10 μm of well) and solubilized ZnO concentration (mg ml^{−1}) for initial solubilization to 0, 30, 60 and 90 μM (0.0, 2.442, 7.327 and 4.885 μg ml^{−1}, dissolved fraction = 0.0, 0.2442, 0.4485, and 0.7327) with no further dissolution during transport are shown in Fig. 7a and b. Solid mass and dissolved exposure concentrations with initial solubility of 20 μM (1.63 μg/ml, dissolved fraction = 0.163) and continuous constant dissolution of 0, 2, 4 and 6 μM h^{−1} (0.0, 0.16, 0.32 and 0.48 μg ml^{−1} h^{−1}, fraction dissolved = 0.0, 0.016, 0.032 and 0.048 h^{−1}) are shown in Fig. 7c and d. Solid and dissolved exposure mass concentrations with initial solubility of 20 μM and further linear increase in dissolution up to 12 h to maximum total dissolution of 20, 30, 60 and 90 μM (1.63, 2.442, 7.327 and 4.885 μg ml^{−1}, dissolved fraction = 0.163, 0.2442, 0.4485, 0.7327) at 12 h, and remaining constant thereafter, are shown in Fig. 7e and f.

## Discussion

We have experimentally validated the predictions of the DG model by analysis of frozen and sectioned columns of test materials, and shown that the results of the DG model compare closely with the CFD model. Although CFD may be considered the gold standard for transport modeling, it is computationally-intensive, requiring expensive software and hardware, and long compute times, whereas the DG model requires only MATLAB, and a typical personal computer, and runs in a few minutes or less.

We have demonstrated here that the suspension-cell boundary condition at the bottom of the well is a critical determinant of predicted transport and dose. With a reflective boundary, sedimentation generates a gradient that drives diffusion upward, which reducing net downward transport and leading to an equilibrium gradient beyond which no further net transport occurs. Conversely, when the boundary is modeled as perfectly adherent, particles are sequestered to a greater or lesser extent depending on the value of *K*
_{D}, reducing, and in the extreme case inverting the diffusion-driving gradient, and thus increases net transport toward the bottom. Our results show that values of *K*
_{D} in the nanomolar range or lower are required to appreciably affect delivered dose metrics (Fig. 4). *K*
_{D} values of this order are typical of specific protein-protein interactions. For example, *K*
_{D} for mouse IgG anti-biotin antibody binding to biotin-BSA targets was found to be 2.3 nM [36], and binding of wheat germ agglutinin to epithelial cell glycoproteins has been reported as 0.32 μM [37]. By contrast, non-specific protein-protein interactions are considerably weaker, with *K*
_{D} values in the millimolar range [32–34]. Interactions between nanoparticle corona proteins and cell surface biomolecules in cell culture would likely fall predominantly into the latter non-specific category, and cell surface adsorption of particles would therefore not be expected to significantly alter transport. Moreover, the same proteins constituting the ENM protein corona and participating in interactions with cell surface biomolecules, would also be highly abundant in free form within the media. To the extent that significant affinity existed between these proteins and cell surface biomolecules, the relatively abundant free protein molecules would occupy most of the available binding sites. In addition, we previously showed that when epithelial cells were incubated with metal oxide and silica-coated metal oxide ENMs [38] less than 0.1 % of the suspended material was either associated with or had transmigrated through cells after 24 h [18], consistent with non-specific binding and negligible effects of adsorption on transport. These data and considerations taken together suggest that a reflective boundary condition is most likely for suspensions of most industry relevant metal and metal oxide ENMs suspended in culture media. Nevertheless, further work is needed to fully assess the role of adsorption and uptake under various conditions and with various materials and cells. We have previously shown, for example, that in the absence of a protein corona (i.e. in protein-free media), metal oxide ENM agglomerates adhere more tightly to cell surfaces than they do when protein is present [21]. Such interactions between naked metal or metal oxide particles and cell surface proteins could conceivably obtain *K*
_{D} values sufficiently small as to affect transport. Likewise, conformational changes in serum proteins caused by ENM adsorption could potentially increase their affinity for cell surface proteins, leading to greater binding and uptake of ENMs [39]. Further investigation to assess *K*
_{D} for various ENM/cell systems are underway in our lab.

The DG model provides accurate modeling of transport for polydisperse suspensions by treating particle size as a pair of arrays (diameter and corresponding fraction of total mass), calculating transport of each size species separately, and summing the resulting mass concentrations over all species. As might be expected, and shown in Fig. 6, results obtained using polydisperse size distribution were in most cases substantially different from those obtained using a volume-weighted mean size. These differences were considerable for all materials examined, in spite of the relatively small polydispersity indices (PdI ≤ 0.3) measured by DLS for three of the four ENMs (Table 1). Although the differences for some materials at certain time points were negligible (e.g. CeO_{2} at 24 h and TiO_{2} at 12–15 h), it is clear that taking polydispersity into consideration improves accuracy of predicted dose metrics.

We should note that the effective density obtained by VCM and used in these simulations is an average density across all agglomerates, and it is likely that a distribution of effective densities is in fact superimposed upon the agglomerate size distribution. In addition, we have assumed that the composition and effective density of agglomerates are constant over the course of the simulated incubation. Although the close match between simulation and validation measurements presented here suggest that these approximations are reasonable, further studies to provide more detailed effective density characterization as a function of time and size are warranted.

In a recently-proposed dosimetry model, dynamic agglomeration of silver nanoparticles was estimated using stochastic modeling, based upon theoretical collision rates and estimated attractive forces [40]. Though promising, the detailed physicochemical characterization required for these estimates are not readily obtainable for the vast numbers of ENMs that require testing. Furthermore, we have previously shown that when ENM suspensions are prepared by sonication in DI water above critical delivered sonication energy (*DSE*
_{cr}), followed by dilution into culture media, the size distributions of the resulting suspensions are stable over time [11]. However, in cases where substantial concentration- or time-dependent changes in agglomeration state occur, corresponding adjustments to the model could easily be made, incorporating empirical measurements of agglomeration state over time and/or concentration, rather than relying upon theoretical estimates. Because of the iterative and compartmental design of the DG model, it is well suited to such modifications. Indeed, the associating protein system-transport modeling for which its antecedent was designed includes re-calculation of molecular species concentrations (based on stoichiometry, association constants and the law of mass action) at the end of each iteration [24–30]. An analogous approach could be employed to accommodate concentration or time-dependent changes in ENM agglomeration state.

The majority of ENMs are insoluble in aqueous solution, but when dissolution does occur, as it does with ZnO, this must be accounted for in dosimetry and transport modeling. Although it may be possible to theoretically estimate rates of dissolution over time in some specific cases [40], empirical measurements of ENM dissolution by ICP-MS, such as those reported by Xia et al. for ZnO [35], can readily be used to obtain dynamic dissolution data that can then be incorporated in to the transport model. In addition to decreasing solid ENM concentrations, solubilization may affect the size or distribution of sizes of agglomerates. The DG model allows specification of initial dissolution as well as either a constant dissolution rate or dissolution over a specified time to maximum fraction (Fig. 7). Beyond initial dissolution, assumed to have occurred prior to characterization of the suspension, further dissolution was modeled as a decrease in size of agglomerates corresponding to the decrease in mass. Additional studies to characterize size distribution over time for soluble materials will be helpful in refining the model for transport of such materials.

## Conclusions

In summary, both the CFD and the newly-developed DG models provide versatile and robust tools for accurately determining *in vitro* ENM concentration and deposition metrics over time. The DG model builds upon important earlier contributions in this area, allowing nanotoxicologists to account for polydispersity and solubilization of ENMs in suspension, as well as the effect of particle-cell binding at bottom of the well. Such advanced fate and transport models will enable nanotoxicologists to develop integrated and standardized dosimetric approaches, and will help to improve the accuracy and reliability of *in vitro* toxicological assays for engineered nanomaterials.

## Methods

### Nanomaterials and characterization.

ENMs investigated are listed in Table 1. SiO_{2}, Fe_{2}O_{3}, and CeO_{2} ENM powders were generated in-house by flame spray pyrolysis using the Harvard Versatile Engineered Nanomaterial Generation System as previously described (VENGES) [41, 42]. TiO_{2} and ZnO ENM powders were obtained from EVONIK (Essen, Germany) and Alfa Aesar (Ward Hill, MA, USA) respectively.

Specific surface area, *SSA*, was determined by the nitrogen adsorption/Brunauer-Emmett-Teller (BET) method using a Micrometrics Tristar 3000 (Micrometrics, Inc., Norcross, GA, USA). Equivalent primary particle diameter, *d*
_{BET}, was calculated, assuming spherical particles, as

where *ρ*
_{p} is the particle density, which was obtained for each particle from the densities of component materials, at 20 °C, reported in the CRC handbook of Chemistry and Physics [43]. Particle crystal size and diameter was also determined by X-ray diffraction using a Scintag XDS2000 powder diffractometer (Scintag Inc., Cupertino, CA, USA), reported here as *d*
_{XRD}.

### ENM dispersal and characterization in suspension

Dispersions were prepared based on a protocol recently developed by the authors [11]. Briefly, sonication was performed in deionized water (DI H_{2}O) using the critical dispersion sonication energy (*DSE*
_{cr}), which was determined as previously described for each ENM [11]. ENMs were dispersed at 5 mg cm^{−3} in 3 ml of solute in 15 ml conical polyethylene tubes, by sonication with a cup horn Branson Sonifier S450-A (Branson Ultrasonics Corporation, Danbury, CT) (maximum power output 400 W at 60 Hz, continuous mode, output level 3, power delivered to sample: 1.25 W). Stock DI H_{2}O suspensions were then diluted to final concentrations in either RPMI + 10 % FBS or PBS + 0.5 % BSA and vortexed for 30 seconds. Dispersions were analyzed by DLS for determination of hydrodynamic diameters (*d*
_{H}) and polydispersity indices (PdI) using a Zetasizer Nano-ZS (Malvern Instruments, Worcestershire, UK).

### Effective density by volumetric centrifugation method (VCM)

Effective densities of ENMs were determined as previously described [12]. Briefly, 1.0 ml samples of 100 μg cm^{−3} suspensions of ENMs were dispensed into packed cell volume (PCV) tubes (Techno Plastic Products, Trasadingen, Switzerland) and centrifuged at 2,000 × *g* for one hour to pellet the suspended material in the bottom (capillary) of the PCV tube.

Agglomerate pellet volumes, *V*
_{pellet}, were measured using the *easy-measure* device from the PCV tube manufacturer. This device roughly resembles a thick steel ruler. The front face is etched along the top edge with graduations at 0.025 μl intervals from 0 to 5 μl. On the back face is a shelf that declines from near the top edge of the ruler at the left end, to near the bottom edge at the right end. The PCV tube rests with its bottom (capillary) end on the ramp, while being held in a sliding holder that wraps around the ruler. The front of the holder contains a lens that magnifies the view of the PCV tube capillary and the ruler graduation marks. The tube in its holder is slid along the ramp until the top edge of the pellet in the capillary is aligned with the top edge of the ruler, and the volume is read from the graduation mark coincident with that horizontal position.

Media densities were calculated from the mass of a 50 ml sample by subtracting the weight of a 50 ml volumetric flask from the weight of the same flask containing 50 ml of media. Effective agglomerate densities were calculated from *V*
_{pellet} values of triplicate samples for each ENM, assuming a stacking factor, SF, of 0.634.

### Validation by frozen sections

Cylindrical receptacles for suspensions were created as follows. Polypropylene 15 ml conical tubes were cut down to approximately 5 cm in height and lined with polyethylene film. Approximately 3 ml of liquid paraffin wax was dispensed into each tube. A smooth 7 mm diameter hardwood dowel was inserted vertically into the paraffin in each tube, and maintained perpendicular to the bench surface until cooled, after which the dowel was removed, leaving a smooth cylindrical cavity of the same dimensions (7 mm in diameter, ~4 cm in height).

Nanomaterial suspensions were dispensed into the paraffin receptacles to create a column 10 mm in height (184 μl for the cylinder 7 mm in diameter), tubes were sealed with plastic wrap and incubated at 22 º C (the temperature employed for all simulations) for the designated times (24 h for TiO_{2}, 48 h for CeO_{2} and 72 h for Fe_{2}O_{3}). Plastic-wrapped paraffin molds containing transported suspensions were carefully removed from tubes, flash frozen by immersing in liquid nitrogen for one minute, and moved to a pre-cooled (−30 °C) cryostat (CM-3000, Leica Biosystems, Buffalo Grove, IL, USA). Plastic wrap was removed and frozen paraffin broken away to extract the solid pellets of suspension, which were then mounted top side down in OCT fluid on specimen chucks. Pellets were sectioned at 25 μm and ten sections were combined per sample, representing a total thickness of 250 μm each. Samples were allowed to thaw and their ENM concentrations determined from optical absorbance and standard curves obtained with a NanoDrop 2000 UV–vis spectrophotometer (Thermo Scientific, Wilmington, DE, USA).

### Distorted Grid model

#### Input and Initialization

Input data required by the DG model include the initial mass concentration of the ENM, density of the bulk ENM material, effective density of agglomerates (e.g. from the VCM), height of the suspension column, the size or size distribution of the suspended particles (i.e. diameter and corresponding fraction of total solute for each particle species *j*), media density and viscosity, and temperature.

The vertical cylindrical column representing the cell culture well is divided by *n* + 1 horizontal boundaries into *n* compartments of equal height. The top and bottom of the well are the ‘impermeable’ boundaries 1 and *n* + 1, respectively; thus each compartment *i* is enclosed between boundaries *i* and *i* + 1. The initial distribution of each particle species *j* is described by assigning a value *C*
_{
i,j
} to the mean concentration in each compartment. Since the initial suspension is homogeneous throughout the well, this is equal to the initial concentration *C*
_{0,j
} of species *j* in the well for all cells. Thus,

#### Modeling diffusion

According to Fick’s first law a substance flows from a region of higher concentration to a region of lower concentration at a rate (per unit cross-sectional area) proportional to the magnitude of the concentration gradient at the boundary of the two regions:

where *C* is the concentration (kg m^{−3}), *z* is the position (m), and the proportionality constant *D* is the diffusion coefficient (m^{2} s^{−1}), which is defined by the Stokes-Einstein equation (Eq. 2).

During a short time interval, Δ*t*, the mass of solute species j moving upward out of compartment *i*, and across boundary *i*, *M*
_{out}, can be represented as

where *A* is the cross-sectional area of the column, *D*
_{
i,j
} is the diffusion coefficient of particle species *j* at boundary *i*, *C*
_{
i − 1,j
} and *C*
_{
i,j
} are the concentrations of particle species *j* in the two compartments, and Δ*z* is the distance between the centers of the two compartments. Likewise, the mass of particle species *j* moving upward into compartment *i*, across boundary *i* + 1, *M*
_{in}, is given as

These mass movements of material into and out of a compartment in single round of diffusion are depicted in Additional file 1: Figure S1 b. The net change in concentration of species *j* in compartment i, Δ*C*
_{
i,j
} can be defined as

Substituting Eqs. 6 and 7 into Eq. 8, and simplifying gives

Note that the cross-sectional area *A* cancels and is thus not needed for this calculation. In each iteration of the model, Eq. 9 is applied to each compartment, and the new concentration of species *j* in compartment *i* after time Δ*t* is calculated as:

Although the diffusion coefficients for a given species *j* at each boundary *D*
_{
i,j
} are considered identical for the systems modeled here, this may not always be the case. Here we have assumed that diffusion coefficients are independent of concentration and that possible cross-diffusion coefficients are negligible. However, the model has been implemented with an array of diffusion coefficients (one for each boundary) for each agglomerate species *j*, allowing future accommodation of more complex systems. The implementation and also includes an option to add a non-linear concentration dependence factor for the diffusion coefficient, represented as

where *k* is the concentration dependence factor. For all simulations reported here the value of *k* was set to its default value of 0. We also assume for polydisperse suspensions that the concentrations of different species do not affect diffusion coefficients of one another, and that such systems can therefore be modeled by simultaneously but independently simulating transport of a each species. We further assume that particles or agglomerates do not dissociate and re-associate during transport, and that surface charge effects (attraction and repulsion) are negligible. However, the segmented, iterative design of the model would allow such considerations to be accommodated relatively easily, by adjusting concentrations and diffusion coefficients at the end of each short round of simulated transport.

#### Modeling sedimentation

The sedimentation velocity, *v*
_{s}, of a particle under gravity is defined by Eq. 1. The sedimentation coefficient of a particle, *S*, is defined as the ratio of a particle’s sedimentation velocity to the acceleration applied to it. Thus, for a particle sedimenting under gravity,

where *g* is the acceleration due to gravity. The downward vertical displacement of a particle with sedimentation coefficient *S* in Δ*t* seconds is given by

In the previously-described method for simulation of velocity sedimentation that is the basis for this model [26–31], sedimentation was modeled by displacing compartment boundaries in the direction of sedimentation, and subsequently re-calculating solute concentrations in the bounded compartments. Because of this boundary migration the model did not adequately handle the situation at the bottom of the cell (bottom of the well in our model), since over the course of the simulation a number of boundaries would accumulate at the bottom of the cell and collapse associated compartments. In our model, since the concentration in the compartments at the bottom of the cell are of the most interest, it was necessary to revise the sedimentation component of the model so that solute (particles) could be moved between compartments without permanently moving compartment boundaries. This is accomplished by calculating, in each round of simulated sedimentation, the distance by which, for each particle species *j*, each compartment *i* would be displaced during the simulation time interval Δ*t*:

where Δ*z*
_{
i,j
} is the downward displacement of compartment *i*, and *S*
_{
i,j
} is the sedimentation coefficient of particle species *j* at the center (along the *z* axis) of compartment *i*. By the sedimentation of particle species *j* into compartment *i* from compartment *i* − 1 above, the concentration in compartment *i* is increased by the product of the concentration in compartment *i* – 1 and the fraction of the compartment height that is displaced into compartment *i*:

where *h* is the height of a compartment. Likewise, by sedimentation of particle species *j* out of compartment *i*, the concentration in compartment *i* is decreased by

The new concentration is thus calculated for compartment *i* by adding Δ*C*
_{in} to, and subtracting Δ*C*
_{out} from its previous concentration, which yields

A single round of sedimentation for one compartment is depicted in Additional file 1: Figure S1b. In order to prevent compartment displacement in one iteration from completely overtaking the next compartment for the user-selected Δ*t*, the displacement that would be experienced by the particle species with the largest sedimentation coefficient in that time is calculated at the outset of the simulation. If that displacement is greater than *h*/2, then the largest time step allowable (producing displacement *h*/2) is calculated and used in place of the specified time step.

As with diffusion, although in the simulations described here it was assumed that the sedimentation coefficients at each boundary are identical, and that there is no dependence of sedimentation coefficients on concentration, the model is implemented with an array of sedimentation coefficients (one per boundary per species) and allows the introduction of concentration dependence. The sedimentation coefficient frequently varies in a non-linear way with concentration [44]. This relationship is often described by:

where *k* is a solute-dependent constant. For all simulations reported here the value of *k* was set to its default value of 0.

#### Iteration and output

The simulation proceeds by alternating rounds of diffusion (Eqs. 9 and 10) and sedimentation (Eq. 17) of duration Δ*t* until the selected cumulative time is obtained. At user-selected intervals during the simulation, as well as the end of simulation, the concentrations in each compartment, as well as derived mass, particle number and surface area dose metrics are exported to an excel file.

### Corrections to frictional coefficients

The DG model includes corrections to particle frictional coefficients in the form of ratios of corrected to ideal Stokes frictional coefficient, *f*/*f*
^{0}. Dividing the diffusion coefficients and sedimentation velocities by these ratios provides the associated correction for transport calculations.

As particle size approaches the mean free path of the media, *λ*, (≈ 2.5 × 10^{−10} m for water), slipping of solvent molecules at the particle surface (where ideally solvent has zero velocity) requires a correction referred to as the “slip”, or Cunningham correction factor, *C*
_{
c
} [45, 46], such that:

where

The magnitude of this correction for typical ENM particles and agglomerates is almost negligible (e.g. *C*
_{
c
} ≈ 1.006 for *d* = 100 nm), and only becomes significant for particles with diameters <5 nm (*C*
_{
c
} ≈ 1.12 for *d* = 5 nm, and *C*
_{
c
} ≈ 1.30 for *d* = 2 nm). In the DG model this correction is applied by default unless otherwise specified by the user.

Additional corrections for non-spherical shape, solvation, and surface roughness can also be included. Whereas diameters obtained by DLS, being calculated from measured diffusion coefficients using the Stokes-Einstein equation (Eq. 2), are by definition hydrodynamic diameters (*d*
_{H}), which account for these effects, diameters measured by other methods (TEM, TRPS) may more accurately be considered diameters of idealized (unsolvated, smooth) spheres of equivalent volume (*d*
_{E}). A correction of the form *f*/*f*
^{0} for non-spherical shape, commonly referred to as a dynamic shape factor, *χ*, can be specified in the model input either directly, or by indicating one of four specific shapes for which *χ* is known or can be calculated from known equations. These shapes include a cube, for which *χ* = 1.08 [45], and a prolate ellipsoid, oblate ellipsoid and circular cylinder with a ratio of major to minor axis length *P* = *a*/*b*, for which *χ* can be calculated as follows [44]:

In most cases these corrections are relatively small. For example, for a prolate ellipsoid, oblate ellipsoid and cylinder with *P* = 3, *χ* = 1.11, 1.10 and 1.27, respectively.

For other shapes and solvation and roughness factors, estimation of *f*/*f*
^{0} is a complex topic beyond the scope of this paper. The DG model nevertheless provides options for user-specified shape, solvation and roughness correction factors.

### Calculation of volume-weighted hydrodynamic diameter distribution and average

The volume-weighted fraction, *f*
_{V,i
}, for each diameter, *d*
_{H,i
}, corresponding to the number fraction, *f*
_{N,i
}, obtained from DLS were calculated as

and the volume-weighted average hydrodynamic diameter, 〈*d*
_{H}〉_{V}, as

### Calculation of particle number and surface concentrations

The effective density of an ENM agglomerate, *ρ*
_{EV}, which we can measure by the VCM method, can be expressed as a weighted average of the media and raw nanomaterial densities:

where *f*
_{V,NM} is the volume fraction of the agglomerate consisting of the nanomaterial, *ρ*
_{NM} is the density of raw nanomaterial, and *ρ*
_{m} is the density of the media. Solving for *f*
_{V,NM} yields

The volume concentration of raw nanomaterial, *C*
_{V,NM}, can be expressed as the product of *f*
_{V,NM} and the volume concentration of the agglomerate, *C*
_{V,agg}, which is also equivalent to the mass concentration of raw nanomaterial, *C*
_{NM}, divided by its density:

Dividing through by *f*
_{V,NM} yields

Dividing by the volume of a single agglomerate, *πd*
^{3}_{H}
/6 (assuming spherical shape) gives the agglomerate particle number concentration as

where *d*
_{H} is the hydrodynamic diameter of the agglomerate. Multiplying this by the surface area of a spherical agglomerate, *πd*
^{2}_{H}
, yields the agglomerate surface area concentration,

Finally, assuming that agglomerates are homogeneous (i.e. that primary nanomaterial particles are evenly distributed throughout the agglomerate), the fraction of agglomerate surface made up by nanomaterial surface is identical to the volume fraction of nanomaterial in the agglomerate, namely *f*
_{V,NM}, and the nanomaterial surface concentration, *C*
_{S,NM} is the product of *C*
_{S,agg} and *f*
_{V,NM}. Thus,

### Implementation of the Langmuir isotherm adsorption

Binding of agglomerates implemented as a Langmuir isotherm adsorption process, in which the equilibrium dissociation constant, *K*
_{D}, is defined as

where *k*
_{d} and *k*
_{a} are the rate constants for desorption and adsorption, *θ* is the surface coverage, defined as the fraction of surface sites occupied by adsorbed particles, and [*P*] is the molar concentration of particles [31]. Rearranging Eq. 33, the fraction of surface sites occupied, *θ*, can be expressed as

To implement Langmuir isotherm adsorption in the DG model, at the end of each round of simulated diffusion and sedimentation, the particle (agglomerate) molar concentration [*P*] (mol L^{−1}) in the bottom compartment is calculated from the particle (agglomerate) mass concentration *C*
_{
P
} as

where *N*
_{
A
} is Avagadro’s number, *ρ*
_{EV} is the effective density (agglomerate density) and *r* is the particle radius. The fraction of surface sites occupied is then calculated from Eq. 34. The factor of 10^{−3} in the numerator provides conversion of units from moles m^{−3} to moles L^{−3}. Because the nature and distribution of the surface adsorption sites are hypothetical in our model, we assume that the fraction of sites occupied is equivalent to the fraction of surface area occupied. Particles available for adsorption are limited to those present in the bottom simulation compartment with a height of Δ*z*. The fraction of bottom surface that can be covered by the particles within this compartment can thus be estimated as

where *A*
_{m} is the cross-sectional area per mass of particle:

The fraction of particles in the bottom compartment that are bound is then calculated as

The concentration (mg cm^{−3}) of particle bound within the bottom compartment is then calculated as

and the concentration of free particle is

Since bound particles do not contribute to the concentration gradient driving diffusion, the concentration of free particles is used for calculating diffusion in the bottom compartments.

### Computational Fluid Dynamics model

For CFD mass transport predictions, the Eulerian–Lagrangian solution of the Navier–Stokes equation was used to model the fluid phase as a continuum, while treating the dispersed particles as a discrete phase, and tracking each particle in the Lagrangian coordinate system within the calculated flow field. Particle-to-particle interactions were neglected and it was assumed that the dispersed phase occupied a much smaller volume fraction than the fluid phase. The Eulerian–Lagrangian approach calculates the particle trajectories by integrating the forces acting on each particle in a Lagrangian reference frame [47]. In the corresponding force balance equation, the particle inertia is equal to the sum of other forces acting on the particle:

Here, \( \overrightarrow{F} \) is an additional acceleration term, \( \overrightarrow{g} \) is gravitational acceleration, *ρ* is the fluid density, *ρ*
_{p} is the particle density, \( \overrightarrow{u} \) is the fluid phase velocity, \( {\overrightarrow{u}}_{\mathrm{p}} \) is the particle velocity, and \( {F}_{\mathrm{D}}\left(\overrightarrow{u}-{\overrightarrow{u}}_{\mathrm{p}}\right) \) is the drag force per unit particle mass, the general expression for which is

Here, *μ* is the fluid viscosity, *d*
_{p} is the particle diameter, *C*
_{
D
} is a drag correction factor, and *Re* is the relative Reynolds number, which is defined as:

In the model design employed in this work, laminar diffusion, Brownian motion and gravitational force were considered, while other forces, such as thermophoretic Force and Saffman’s lift Force were neglected. Because *d*
_{p} < 500 nm, Brownian motion was included in the additional force term (\( \overrightarrow{F} \)) [48, 49], and the Stokes-Cunningham Drag Law, which accounts for non-zero relative velocity of fluid at the particle surface (“slip”), was chosen [50]. In this specific drag law, *F*
_{D} is defined as:

where the factor *C*
_{
c
} is the Cunningham correction factor, as defined above (Eq. 20). Note that *C*
_{
c
} in the denominator here replaces the *C*
_{D}
*Re*/24 term in the general drag force expression (Eq. 42).

In this study, the simulation domain encompassed as a cylinder space with constant dimensions of 1.5 mm diameter and 1.0 mm height. To reduce the computational resources required, transport was modeled within a *π*/4 radian sector of the circular cylinder representing the cell culture well. A structured grid with a total of 3,586,000 elements was employed to represent the computational spatial domain. Grid independency analysis was carried out by varying grid numbers in all three dimensions. Performance at different grid aspect ratios, ranging from 1.0 to 3.0 was analyzed to minimize computational resources. For the final simulations, grid sizes of 2.5 μm and 5.0 μm were selected for *z* direction and *x,y* directions of Cartesian coordinates, respectively, resulting in the uniform grid aspect ratio of 2.0.

In the sector-shaped simulation domain, a symmetry boundary condition was employed for the cutting planes (flat vertical sector surfaces). All other surfaces, including the top, bottom and outer curved cylinder surfaces, were modeled as no-slip wall boundaries with zero roughness and a constant temperature of 22 °C. A reflective boundary condition was imposed at all walls for the discrete phase mode [47, 50], with a default constant value of 1.0 for both normal and tangential coefficient of restitution [47].

A three-dimensional structured mesh with the finite volume method and double precision option were utilized to discretize the computational domain, and to describe the mass and momentum transport for each cell. The SIMPLE algorithm was used to solve the pressure and velocity components. Spatial discretizations for momentum and energy equations employed the second-order upwind scheme. Time-dependent terms were interpolated using the implicit second-order interpolation scheme. All particles were injected at *t* = 0 and assigned an initial velocity magnitude of 0.0 m s^{−1}.

## References

Nel A, Xia T, Mädler L, Li N. Toxic potential of materials at the nanolevel. Science. 2006;311:622–7.

Bakand S, Hayes A, Dechsakulthorn F. Nanoparticles: a review of particle toxicology following inhalation exposure. Inhal Toxicol. 2012;24:125–35.

Pirela S, Molina R, Watson C, Cohen JM, Bello D, Demokritou P, et al. Effects of copy center particles on the lungs: A toxicological characterization using a Balb/c mouse model. Inhal Toxicol. 2013;25:498–508.

Krewski D, Acosta D Jr, Andersen M, Andersen H, Bailar JC 3rd, Boekelheide K, et al. Toxicity testing in the 21st century: a vision and a strategy. J Toxicol Environ Health B Crit Rev. 2010;13:51–138.

Pyrgiotakis G, McDevitt J, Gao Y, Branco A, Eleftheriadou M, Lemos B, et al. Mycobateria inactivation using Engineered Water Nanostructures (EWNS). Nanomedicine. 2014;10:1175–83.

George S, Xia T, Rallo R, Zhao Y, Ji Z, Lin S, et al. Use of a high-throughput screening approach coupled with in vivo zebrafish embryo screening to develop hazard ranking for engineered nanomaterials. ACS Nano. 2011;5:1805–17.

National Science and Technology Council Committee on Technology, Subcommittee on Nanoscale Science, Engineering and Technology, NNI (National Nanotechnology Initiative). Environmental, Health, and Safety (EHS) Research Strategy 2011. Available at http://www.nano.gov/sites/default/files/pub_resource/nni_2011_ehs_research_strategy.pdf.

Demokritou P, Gass S, Pyrgiotakis G, Cohen JM, Goldsmith W, McKinney W, et al. An in vivo and in vitro toxicological characterization of realistic nanoscale CeO

_{2}inhalation exposures. Nanotoxicology. 2013;7:1338–50.Gangwal S, Brown JS, Wang A, Houck KA, Dix DJ, Kavlock RJ, et al. Informing selection of nanomaterial concentrations for ToxCast in vitro testing based on occupational exposure potential. Environ Health Perspect. 2011;119:1539–46.

Oberdörster G. Nanotoxicology: in vitro-in vivo dosimetry. Environ Health Perspect. 2012;120:A13.

Cohen J, DeLoid G, Pyrgiotakis G, Demokritou P. Interactions of engineered nanomaterials in physiological media and implications for in vitro dosimetry. Nanotoxicology. 2012;7:417–31.

Deloid G, Cohen JM, Darrah T, Derk r, Rojanasakul l, Pyrgiotakis G, Wohlleben W, Demokritou P. Estimating the effective density of engineered nanomaterials for in vitro dosimetry. Nat Commun. 2014, doi:10.1038/ncomms4514.

Pal AK, Aalaei I, Gadde S, Gaines P, Schmidt D, Demokritou P, et al. High resolution characterization of engineered nanomaterial dispersions in complex media using tunable resistive pulse sensing technology. ACS Nano. 2014;8:9003–15.

Teeguarden JG, Hinderliter PM, Orr G, Thrall BD, Pounds JG. Particokinetics in vitro: dosimetry considerations for in vitro nanoparticle toxicity assessments. Toxicol Sci. 2007;95:300–12.

Rushton EK, Jiang J, Leonard SS, Eberly S, Castranova V, Biswas P, et al. Concept of assessing nanoparticle hazards considering nanoparticle dosemetric and chemical biological response metrics. J Toxicol Environ Health. 2010;73:445–61.

Hinderliter PM, Minard KR, Orr G, Chrisler WB, Thrall BD, Pounds JG, et al. ISDD: A computational model of particle sedimentation, diffusion and target cell dosimetry for in vitro toxicity studies. Part. Fibre Toxicol. 2010;7, doi:10.1186/1743-8977-7-36

Cohen JM, Teeguarden JG, Demokritou P. An integrated approach for the in vitro dosimetry of engineered nanomaterials. Part Fibre Toxicol 2014;11, doi:10.1186/1743-8977-11-20

Cohen JM, Derk R, Wang L, Godleski J, Kobzik L, Brain J, Demokritou P. Tracking translocation of industrially relevant engineered nanomaterials (ENMs) across alveolar epithelial monolayers in vitro. Nanotoxicology. 2014;8, doi:10.3109/17435390.2013.879612

Pal AK, Bello D, Cohen J, Demokritou P. Implications of in vitro dosimetry on toxicological ranking of low aspect ratio engineered nanomaterials. Nanotoxicology 2015;12, doi:10.3109/17435390.2014.986670

Cohen, JM, DeLoid, GM, Demokritou, P. A critical review of in vitro dosimetry for engineered nanomaterials. Nanomedicine (Lond). 2015; Epub ahead of print.

Pyrgiotakis G, Blattmann CO, Pratsinis S, Demokritou P. Nanoparticle-nanoparticle interactions in biological media by atomic force microscopy. Langmuir. 2013;29:11385–95.

Pyrgiotakis G, Blattman CO, Demokritou P. Real-time nanoparticle-cell interactions in physiological media by atomic force microscopy. ACS Sustain Chem Eng. 2014;2:1682–90.

Sterling Jr MC, Bonner JS, Ernest AN, Page CA, Autenrieth RL. Application of fractal flocculation and vertical transport model to aquatic sol-sediment systems. Water Res. 2005;39:1818–30.

Cox DJ. Sedimentation of an initially skewed boundary. Science. 1966;152:359–61.

Cox DJ. Computer simulation of sedimentation in the ultracentrifuge. III. Concentration dependent sedimentation. Arch Biochem Biophys. 1967;119:230–9.

Cox DJ. Computer simulation of sedimentation in the ultracentrifuge. IV. Velocity Sedimentation of self-associating solutes. Arch Biochem Biophys. 1969;129:106–23.

Cox DJ. Computer simulation of sedimentation in the ultracentrifuge. V. Ideal and non-ideal monomer-trimer systems. Arch Biochem Biophys. 1971;142:514–26.

Cox DJ. Computer simulation of sedimentation in the ultracentrifuge. VI. Monomer-tetramer systems in rapid chemical equilibrium. Arch Biochem Biophys. 1971;146:181–95.

Cox DJ. Computer stimulation of sedimentation in the ultracentrifuge. VII. Solutes undergoing indefinite self-association. Arch Biochem Biophys. 1974;160:595–602.

DeLoid G, Computer Simulation of Velocity Sedimentation for Mixed Associating Protein Systems, 1984, Thesis, Kansas. State U, https://ia800300.us.archive.org/19/items/simulationstudie00delo/simulationstudie00delo.pdf

Czepirski L, Balys MR, Komorowska-Czepirska E. Some generalization of Langmuir adsorption isotherm. Internet J Chem. 2000;3:14.

Johansson H, Jensen MR, Gesmar H, Meier S, Vinther JM, Keeler C, et al. Specific and nonspecific interactions in ultraweak protein-protein associations revealed by solvent paramagnetic relaxation enhancements. J Am Chem Soc. 2014;136:10277–86.

Bashir Q, Scanu S, Ubbink M. Dynamics in electron transfer protein complexes. FEBS J. 2011;278:1391–400.

Tang C, Ghirlando R, Clore GM. Visualization of transient ultra-weak protein self-association in solution using paramagnetic relaxation enhancement. J Am Chem Soc. 2008;130:4048–56.

Xia T, Kovochich M, Liong M, Mädler L, Gilbert B, Shi H, et al. Comparison of the mechanism of toxicity of zinc oxide and cerium oxide nanoparticles based on dissolution and oxidative stress properties. ACS Nano. 2008;2:2121–34.

Landry JP, Fei Y, Zhu X. Simultaneous measurement of 10,000 protein-ligand affinity constants using microarray-based kinetic constant assays. Assay Drug Dev Technol. 2012;10:250–9.

Wang W, Yang Y, Wang S, Nagaraj VJ, Liu Q, Wu J, et al. Label-free measuring and mapping of binding kinetics of membrane proteins in single living cells. Nat Chem. 2012;4:846–53.

Gass S, Cohen JM, Pyrgiotakis G, Sotiriou GA, Pratsinis SE, Demokritou P. A safer formulation concept for flame-generated engineered nanomaterials. ACS Sustain Chem Eng. 2013;1:843–57.

Schnitzer JE, Bravo J. High affinity binding, endocytosis, and degradation of conformationally modified albumins. Potential role of gp30 and gp18 as novel scavenger receptors. J Biol Chem. 1993;268:7562–70.

Mukherjee D, Leo BF, Royce SG, Porter AE, Ryan MP, Schwander S, et al. Modeling physicochemical interactions affecting in vitro cellular dosimetry of engineered nanomaterials: application to nanosilver. J Nanopart Res. 2014;16:2616.

Demokritou P, Büchel R, Molina RM, DeLoid GM, Brain JD, Pratsinis SE. Development and characterization of a Versatile Engineered Nanomaterial Generation System (VENGES) suitable for toxicological studies. Inhal Toxicol. 2010;22:2107–16.

Sotiriou GA, Diaz E, Long MS, Godleski J, Brain J, Pratsinis SE, et al. A novel platform for pulmonary and cardiovascular toxicological characterization of inhaled engineered nanomaterials. Nanotoxicology. 2012;6:680–90.

Haynes WM. CRC Handbook of Chemistry and Physics. 92nd ed. Boca Raton, FL, USA: CRC Press/Taylor and Francis; 2011.

van Holde KE, Johnson C, Ho PS. Principles of Physical Biochemistry. 2nd ed. Upper Saddle River, NJ, USA: Prentice Hall; 2005.

Hinds WC. Aerosol Technology, Properties, Behavior, and Measurement of Airborne Particles. 2nd ed. New York, NY, USA: Wiley; 1999.

Allen MD, Raabe OG. Re-evaluation of Millikan’s oil drop data for the motion of small particles in air. J Aerosol Sci. 1982;6:537–47.

ANSYS Inc. ANSYS-FLUENT Theory Guide. Release 14.0. ANSYS, Inc. Canonsburg, PA 2011.

Mendes PJ, Pinto JF, Sousa JM. Numerical Simulation of In-Vitro Dispersion and Deposition of Nanoparticles in Dry-Powder-Inhaler Aerosols. J Nanosci Nanotechno. 2010;10:2791–7.

Ounis H, Ahmadi G, McLaughlin JB. Brownian Diffusion of Submicrometer Particles in the Viscous Sublayer. J Colloid Interface Sci. 1991;143:66–277.

Tabakoff W, Wakeman T. Measured particle rebound characteristics useful for erosion prediction. J ASME. 1982;82:GT-170.

## Acknowledgements

The authors thank Georgios Sotiriou (Particle Technology Laboratory, ETH Zurich) for the electron microscopy imaging, BET and XRD analysis of the ENM samples, and J. Teeguarden (Pacific Northwest National Laboratory, Richland, WA, US) for providing the MATLAB software implementation of the ISDD dose delivery model used in this work. This research project was supported by NIEHS grant (ES-0000002), NSF grant (ID 1235806) and the Center for Nanotechnology and Nanotoxicology at The Harvard School of Public Health. This work was performed in part at the Harvard Center for Nanoscale Systems (CNS), a member of the National Nanotechnology Infrastructure Network (NNIN), which is supported by the National Science Foundation under NSF award no. ECS-0335765.

## Author information

### Authors and Affiliations

### Corresponding authors

## Additional information

### Competing interests

None of the authors have any competing financial interests in the work described in this manuscript.

### Authors’ contributions

GD designed and implemented DG model, performed simulations and cryosection validation experiments, analyzed data, and co-wrote manuscript. JC contributed to experimental design and data analysis, and co-wrote manuscript. GP contributed to experimental design and data analysis, prepared nanoparticles, and co-wrote manuscript. SP: contributed to experimental design, performed particle characterization, and co-wrote manuscript. AP contributed to experimental design, performed particle characterization, and co-wrote manuscript. JL designed and conducted CFD analyses, and co-wrote manuscript. JS contributed to CFD experimental design, supervised CFD analyses, and co-wrote manuscript. PD supervised study, contributed to experimental design and data analysis, and co-wrote manuscript. All authors read and approved the final manuscript.

## Additional file

### Additional file 1: Figure S1.

Distorted Grid model. Figure S2 Cryosection validation method. Figure S3 Volume size distributions from DLS. Figure S4 DG vs VCM-ISDD. (DOCX 2479 kb)

## Rights and permissions

**Open Access** This article is distributed under the terms of the Creative Commons Attribution 4.0 International License (http://creativecommons.org/licenses/by/4.0/), which permits unrestricted use, distribution, and reproduction in any medium, provided you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license, and indicate if changes were made. The Creative Commons Public Domain Dedication waiver (http://creativecommons.org/publicdomain/zero/1.0/) applies to the data made available in this article, unless otherwise stated.

## About this article

### Cite this article

DeLoid, G.M., Cohen, J.M., Pyrgiotakis, G. *et al.* Advanced computational modeling for *in vitro* nanomaterial dosimetry.
*Part Fibre Toxicol* **12**, 32 (2015). https://doi.org/10.1186/s12989-015-0109-1

Received:

Accepted:

Published:

DOI: https://doi.org/10.1186/s12989-015-0109-1