Advanced computational modeling for in vitro nanomaterial dosimetry
© DeLoid et al. 2015
Received: 28 July 2015
Accepted: 12 October 2015
Published: 24 October 2015
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.
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.
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.
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.
KeywordsEngineered nanomaterial Nanotoxicology in vitro dosimetry Nanosafety Hazard ranking Fate and transport modeling
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 . 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.
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 .
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)  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 .
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. . 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 . 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 .
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.
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.
PBS + 0.1 % BSA
0.175 ± 0.013
1.564 ± 0.009
RPMI + 10 % HS
0.755 ± 0.203
1.335 ± 0.000
RPMI + 0.5 % BSA
0.310 ± 0.091
1.420 ± 0.004
RPMI + 10 % HS
0.233 ± 0.018
1.251 ± 0.005
Alfa Aesar ZnO
RPMI + 10 % FBS
0.303 ± 0.122
1.650 ± 0.070
Simulation results for CFD and DG results for CeO2 and TiO2 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 CeO2 and 0.001 mg ml−1 for SiO2, resulting in a total of approximately 3.8 × 105 particles in each case. Even at these low values the CFD simulations for CeO2 and SiO2 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 SiO2. 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 Fe2O3, CeO2 and TiO2 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 Fe2O3 and CeO2, and 64 h for TiO2. 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 (cm2 ml−1 and cm2 cm−2). Calculation of number and surface area metrics from mass concentration and deposition metrics are described in methods.
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  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 , in which adsorption is characterized by the equilibrium dissociation constant, K D, as described in the methods section.
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 (SiO2, Fe2O3, TiO2 and CeO2), 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 CeO2 (〈d H〉v =982.1 nm, ρ EV =1.420 g cm−3), they provide slightly different predictions for TiO2(〈d H〉v =397.8 nm, ρ EV =1.251 g cm−3), and strongly different predictions for Fe2O3 (〈d H〉v =234.5 nm, ρ EV =1.335 g cm−3) and SiO2 (〈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 CeO2, 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 SiO2, which forms the smallest agglomerates.
Effect of agglomerate size and density on dosimetry
in which k B is the Boltzmann constant (kg m2 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.
Effect of Polydispersity
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)  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.
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 , and binding of wheat germ agglutinin to epithelial cell glycoproteins has been reported as 0.32 μM . 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  less than 0.1 % of the suspended material was either associated with or had transmigrated through cells after 24 h , 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 . 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 . 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. CeO2 at 24 h and TiO2 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 . 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 . 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 , empirical measurements of ENM dissolution by ICP-MS, such as those reported by Xia et al. for ZnO , 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.
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.
Nanomaterials and characterization.
ENMs investigated are listed in Table 1. SiO2, Fe2O3, and CeO2 ENM powders were generated in-house by flame spray pyrolysis using the Harvard Versatile Engineered Nanomaterial Generation System as previously described (VENGES) [41, 42]. TiO2 and ZnO ENM powders were obtained from EVONIK (Essen, Germany) and Alfa Aesar (Ward Hill, MA, USA) respectively.
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 . 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 . Briefly, sonication was performed in deionized water (DI H2O) using the critical dispersion sonication energy (DSE cr), which was determined as previously described for each ENM . 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 H2O 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 . 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 TiO2, 48 h for CeO2 and 72 h for Fe2O3). 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.
where C is the concentration (kg m−3), z is the position (m), and the proportionality constant D is the diffusion coefficient (m2 s−1), which is defined by the Stokes-Einstein equation (Eq. 2).
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.
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.
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.
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.
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
Calculation of particle number and surface concentrations
Implementation of the Langmuir isotherm adsorption
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
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 .
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.
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.
Open AccessThis 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.
- Nel A, Xia T, Mädler L, Li N. Toxic potential of materials at the nanolevel. Science. 2006;311:622–7.View ArticlePubMedGoogle Scholar
- Bakand S, Hayes A, Dechsakulthorn F. Nanoparticles: a review of particle toxicology following inhalation exposure. Inhal Toxicol. 2012;24:125–35.View ArticlePubMedGoogle Scholar
- 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 CeO2 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.PubMed CentralView ArticlePubMedGoogle Scholar
- 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.PubMed CentralView ArticlePubMedGoogle Scholar
- 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.View ArticlePubMedGoogle Scholar
- 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.PubMed CentralView ArticlePubMedGoogle Scholar
- 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.Google Scholar
- 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.View ArticlePubMedGoogle Scholar
- Cox DJ. Sedimentation of an initially skewed boundary. Science. 1966;152:359–61.View ArticlePubMedGoogle Scholar
- Cox DJ. Computer simulation of sedimentation in the ultracentrifuge. III. Concentration dependent sedimentation. Arch Biochem Biophys. 1967;119:230–9.View ArticlePubMedGoogle Scholar
- Cox DJ. Computer simulation of sedimentation in the ultracentrifuge. IV. Velocity Sedimentation of self-associating solutes. Arch Biochem Biophys. 1969;129:106–23.View ArticlePubMedGoogle Scholar
- Cox DJ. Computer simulation of sedimentation in the ultracentrifuge. V. Ideal and non-ideal monomer-trimer systems. Arch Biochem Biophys. 1971;142:514–26.View ArticlePubMedGoogle Scholar
- Cox DJ. Computer simulation of sedimentation in the ultracentrifuge. VI. Monomer-tetramer systems in rapid chemical equilibrium. Arch Biochem Biophys. 1971;146:181–95.View ArticlePubMedGoogle Scholar
- Cox DJ. Computer stimulation of sedimentation in the ultracentrifuge. VII. Solutes undergoing indefinite self-association. Arch Biochem Biophys. 1974;160:595–602.View ArticlePubMedGoogle Scholar
- 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.Google Scholar
- 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.View ArticlePubMedGoogle Scholar
- 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.View ArticlePubMedGoogle Scholar
- 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.PubMed CentralView ArticlePubMedGoogle Scholar
- 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.PubMed CentralPubMedGoogle Scholar
- 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.PubMedGoogle Scholar
- 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.Google Scholar
- van Holde KE, Johnson C, Ho PS. Principles of Physical Biochemistry. 2nd ed. Upper Saddle River, NJ, USA: Prentice Hall; 2005.Google Scholar
- Hinds WC. Aerosol Technology, Properties, Behavior, and Measurement of Airborne Particles. 2nd ed. New York, NY, USA: Wiley; 1999.Google Scholar
- 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.View ArticleGoogle Scholar
- 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.View ArticleGoogle Scholar
- Ounis H, Ahmadi G, McLaughlin JB. Brownian Diffusion of Submicrometer Particles in the Viscous Sublayer. J Colloid Interface Sci. 1991;143:66–277.View ArticleGoogle Scholar
- Tabakoff W, Wakeman T. Measured particle rebound characteristics useful for erosion prediction. J ASME. 1982;82:GT-170.Google Scholar