Estimation of Nanodiamond Surface Charge Density from Zeta

Subscriber access provided by UIC Library

Article

Estimation of Nanodiamond Surface Charge Density from Zeta Potential and Molecular Dynamics Simulations Zhenpeng Ge, and Yi Wang J. Phys. Chem. B, Just Accepted Manuscript • DOI: 10.1021/acs.jpcb.6b08589 • Publication Date (Web): 24 Nov 2016 Downloaded from http://pubs.acs.org on November 25, 2016

Just Accepted “Just Accepted” manuscripts have been peer-reviewed and accepted for publication. They are posted online prior to technical editing, formatting for publication and author proofing. The American Chemical Society provides “Just Accepted” as a free service to the research community to expedite the dissemination of scientific material as soon as possible after acceptance. “Just Accepted” manuscripts appear in full in PDF format accompanied by an HTML abstract. “Just Accepted” manuscripts have been fully peer reviewed, but should not be considered the official version of record. They are accessible to all readers and citable by the Digital Object Identifier (DOI®). “Just Accepted” is an optional service offered to authors. Therefore, the “Just Accepted” Web site may not include all articles that will be published in the journal. After a manuscript is technically edited and formatted, it will be removed from the “Just Accepted” Web site and published as an ASAP article. Note that technical editing may introduce minor changes to the manuscript text and/or graphics which could affect content, and all legal disclaimers and ethical guidelines that apply to the journal pertain. ACS cannot be held responsible for errors or consequences arising from the use of information contained in these “Just Accepted” manuscripts.

The Journal of Physical Chemistry B is published by the American Chemical Society. 1155 Sixteenth Street N.W., Washington, DC 20036 Published by American Chemical Society. Copyright © American Chemical Society. However, no copyright claim is made to original U.S. Government works, or works produced by employees of any Commonwealth realm Crown government in the course of their duties.

Page 1 of 34

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

The Journal of Physical Chemistry

Estimation of Nanodiamond Surface Charge Density from Zeta Potential and Molecular Dynamics Simulations Zhenpeng Ge and Yi Wang∗ Department of Physics, Chinese University of Hong Kong, Shatin, N.T., Hong Kong E-mail: [email protected] Phone: +852 3943 6355

1 ACS Paragon Plus Environment


1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Abstract Molecular dynamics simulations of nanoparticles (NPs) are increasingly used to study their interactions with various biological macromolecules. Such simulations generally require detailed knowledge of the surface composition of the NP under investigation. Even for some well-characterized nanoparticles, however, this knowledge is not always available. An example is nanodiamond, a nanoscale diamond particle with surface dominated by oxygen-containing functional groups. In this work, we explore using the harmonic restraint method developed by Venable, et al., to estimate the surface charge density of nanodiamonds. Based on the Gouy-Chapman theory, we convert the experimentally determined zeta potential of a nanodiamond to an effective charge density (σef f ), and then use the latter to estimate σ via molecular dynamics simulations. Through scanning a series of nanodiamond models, we show that the above method provides a straightforward protocol to determine the surface charge density of relatively large (>∼100 nm) NPs. Overall, our results suggest that despite certain limitation, the above protocol can be readily employed to guide the model construction for MD simulations, which is particularly useful when only limited experimental information on the NP surface composition is available to a modeler.


Page 2 of 34

Page 3 of 34

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


Introduction The development of nanotechnology has led to a rapid increase in consumer products containing nanomaterials: from 2005 to 2014, the number of such products has grown by over 30 fold. 1 The increased exposure of our environment and biological systems to various nanoparticles (NPs) immediately raises concerns regarding their toxicity. 2–6 Understanding how NPs interact with a cell is therefore of fundamental importance to the development of nanoparticles that are both safer and more efficient in their designed functions. In recent years, computer simulation and modeling has become a standard technique in the study of NP interactions with biomolecules. 7–12 For instance, the cellular entry of nanoparticles has attracted considerable attention and both all-atom and coarse-grained molecular dynamics (MD) simulations, as well as dissipative particle dynamics simulations and continuum modeling, have been employed to study the interactions of NPs varying in size, shape and surface chemistry with a plasma membrane. 10 For MD simulations, a realistic model is often the prerequisite to obtain convincing results comparable with experiments. To this end, detailed knowledge of the surface composition of NPs, including the number of charged functional groups, is required by the modelers. Unfortunately, this information is not always available. For instance, in our own study 13 of nanodiamonds (NDs), which are nanoscale diamond particles with potential applications in biosensing, drug delivery and long-term tracking, 14–16 the exact composition of the ND surface is missing, although it is known from infrared spectra 17,18 that this surface is dominated by oxygen-containing functional groups, such as hydroxyl (OH) and carboxyl (COO− ). The surface charge density of nanodiamonds, which is largely determined by the density of carboxyl groups, cannot be measured from these experiments. For this reason, in our previous work, multiple ND models with different surface functionalization percentages were constructed and the range of their interaction strength with a model membrane was determined. 13 One metric routinely used to characterize the charged state of a nanoparticle is the zeta 3 ACS Paragon Plus Environment


1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 4 of 34

potential. 19–21 Zeta potential is the electrostatic potential on the plane of shear (or the slipping plane). In an ideal case, this plane forms a sheath enveloping a charged particle, and all the material inside the sheath, including the particle itself and the surrounding ions and liquid, moves as a single kinetic unit during an electrophoresis experiment. 19 In practice, the mobility of ions and water exhibits a gradational change over a certain region, 22 and the concept of a discrete slipping plane is at best an ‘effective’, or approximate one. From an electrophoresis experiment, zeta potential is normally obtained via the HelmholtzSmoluchowski equation: 19,23 ζ=

µη εr ε0

(1)

with µ ≡ ν/E, where ν is the drift velocity of the particle under an electric field E, and µ is its electrophoretic mobility, η is the viscosity of the aqueous solution with dielectric constant εr , and ε0 is the vacuum permittivity. Clearly, if zeta potential can be used to estimate the surface charge density of nanodiamonds, then an electrophoresis experiment can provide a valuable piece of information currently missing in their modeling. Since it is not the electrostatic potential on the particle surface, zeta potential cannot directly yield information on surface charge density. Thus, the key to estimating the latter from the former is to establish a link between the two quantities. In some recent studies, 24,25 such a link was established by first constructing a model of the system under investigation and then locating the position of the slipping plane from MD simulations. As shown by Heikkil¨ a et al., once the slipping plane is located, zeta potential can be computed from simulation trajectories and compared with the experimental measurement. 24 However, such an approach may suffer from the uncertainty of the estimated location of the slipping plane, the error in the electrostatic potential calculation, and, as noted in Venable, et al., 26 the lack of a consistent and general theoretical framework relating zeta potential with the Galvani potential obtained from MD simulations. Based on the Gouy-Chapman theory, 27 zeta potential can be directly related to an ‘effective charge density’, which is the sum of charge density on a flat surface and its surrounding 4 ACS Paragon Plus Environment

Page 5 of 34

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


ions enclosed by the slipping plane. Building on such a relation, Venable, et al. 26 developed the harmonic restraint method to measure the effective charge density of a lipid bilayer from MD simulations. The authors applied their method to revise the CHARMM non-bonded interaction parameters between Na+ and lipid oxygen atoms. Unlike the aforementioned approach, the harmonic restraint method does not involve the computation of electrostatic potential or the estimation of the location of the slipping plane. While it was first employed to correct the overbinding effect between ions and lipid headgroups, the method should be equally applicable on a nanoparticle. In this work, we explore using the experimentally measured zeta potential and the harmonic restraint method to estimate the surface charge density of nanodiamonds. Along this line, we constructed a series of ND models with an increasing number of carboxyl functional groups and computed their effective charges from all-atom MD simulations. The results are then compared with the effective charge of a representative nanodiamond determined from its experimentally measured zeta potential. Additionally, we varied the surface functionalization percentage of hydroxyls, in order to evaluate their role in ion binding to the ND surface. Our results indicate that, despite certain limitation, the above approach offers a robust and easy-to-use protocol to determine the surface charge density of a nanodiamond. Such a protocol is particularly useful given that only limited information on ND surface composition is available, in which case a relatively straightforward electrophoresis experiment can be conducted to provide additional guidance to the modelers.

Methods Theoretical Background The theory of electric double layer (EDL), which characterizes the distribution of ions and the electrostatic potential around a charged particle, was developed decades ago. 27 Here, ‘double layer’ refers to one layer of fixed charges on the particle surface and a second layer of excess counterions in the aqueous solution. Central to the EDL



1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 6 of 34

theory is the Poisson’s equation describing the electrostatic attraction (repulsion) of the counterions (coions) to the surface and the Boltzmann’s relation characterizing the diffusion of ions due to thermal motion. 28 The resulting non-linear Poisson-Boltzmann equation has an exact solution in the 1-dimensional case, which is known as the Gouy-Chapman equation for a flat charged surface:

σ=

p

8cN εr ε0 kB T sinh(

eψ0 ) 2kB T

(2)

where c stands for the ion concentration, N is the Avogadro constant, kB is Boltzmann’s constant, T is the temperature, and σ and ψ0 are the surface charge density and surface electrostatic potential, respectively. While the Gouy-Chapman equation describes the relation between ψ0 and σ, the very same relation exists between zeta potential (ζ) and the effective charge density (σef f ) described earlier:

σef f =

p eζ 8cN εr ε0 kB T sinh( ) 2kB T

(3)

The derivation of Eq 3 follows straightforwardly from that of Eq 2, 27 which is reviewed briefly in the supporting information (SI). As shown by the former equation, σef f can be determined solely from ζ measured at a given ion concentration. Based on this feature, we converted the experimentally determined zeta potential to an effective charge density and then used the latter quantity to estimate the surface charge density of a nanodiamond. Such an estimation was achieved with the harmonic restraint method developed by Venable, et al. 26 Specifically, we constructed a series of nanodiamond models at different surface charge densities, and simulated them in ionic solution with an applied electric field E along the x-axis (Fig 1 a). The nanodiamonds were constrained by a harmonic potential with a force constant k. Their displacement ∆d along the direction of the electric field was recorded


Page 7 of 34

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


throughout the simulations. Balancing the forces on a given ND particle yields:

EQef f = k∆d

(4)

where ∆d is the average of ∆d and Qef f is the effective charge of the nanodiamond. For the nanodiamond models used here, Qef f = σef f A, where A is the surface area of a ND slab (see the following section). The result of the above calculations is a curve of σef f vs. σ, from which the surface charge density corresponding to the target effective charge density can then be located. Model Construction In a recent experiment 29 studying the cellular entrance and intracellular translocation of nanodiamonds, the ND particles had a zeta potential of -28 mV and an average hydrodynamic radius of ∼115 nm. Generally, particles at this length scale cannot be modeled in their entirety with all-atom MD simulations. Hence, we constructed several flat ND ‘slabs’ to study a portion of the nanodiamond surface, which may be considered analogous to using a bilayer patch to study a portion of a large membrane vesicle. In both cases, the small average curvature on the particle surface renders it a reasonable approximation to adopt the Gouy-Chapman theory and to use a flat system in MD simulations. Altogether seven ND slabs with an increasing number of carboxyl groups were constructed, all of which have a size of ∼4.5 × 4.4 × 0.9 nm3 (Fig 1 a). Periodic bonds were used to link atoms across the periodic boundary in the x-y plane. The top and bottom surfaces of the ND slabs consisted of carbons with either 3 or 2 C-C bonds (at a ratio of 1:1). Due to steric conflict, the former type of surface carbons could only be functionalized with hydrogen; for the latter type of carbons, 0, 5, 7, 10, 12, 13 and 14 atoms on each surface were selected randomly and functionalized with carboxyl groups while the remaining atoms were functionalized with hydroxyls. The distance between any two carboxyls was at least 1 nm to ensure a relatively uniform charge distribution. Based on the number of carboxyls on each surface of the ND slabs, we will refer to them as 0c, 5c, 7c, 10c, 12c, 13c and 14c, 7 ACS Paragon Plus Environment


1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

respectively, where c stands for a carboxyl group. All ND slabs were then put into a 4.5 × 4.4 × 11.4 nm3 water box. Neutralizing Na+ ions, along with NaCl at a concentration of 0.17 M, were then introduced. The NaCl concentration was chosen to match the ionic strength of 1× phosphate-buffered saline (PBS) used in the electrophoresis experiment. 29 The surface charge density of a given ND slab was computed by dividing the number of carboxyls on a given surface by the area of that surface, i.e., σ = 5/(4.5 × 4.4) for the 5c ND slab. All analysis was performed by averaging results obtained from the two surfaces of a ND slab. Based on the 7c ND slab, we constructed six more systems (Fig S1) to examine the size effect of our calculation: In two systems, the thickness of the water layer was increased from 10 nm to 20 nm or 30 nm with the ND model unchanged; in another two systems, we increased the thickness of the ND slab from 0.9 nm to 1.8 nm or 2.7 nm while keeping the same number of carboxyls and hydroxyls; finally, we doubled or tripled the area of the ND slab, while keeping the density of carboxyls and hydroxyls unchanged. In all six systems, neutralizing Na+ ions were introduced, along with NaCl at a concentration of 0.17 M. We further constructed four more ND slabs to study the effect of hydroxyls on the effective charge of a nanodiamond. Based on the 5c ND slab, we reduced the number of hydroxyls to 75%, 50%, 25% and 0% of the original number, respectively, while keeping the number of carboxyls unchanged. Finally, in order to compare the result of the harmonic restraint method with other effective charge calculation methods, we also constructed two small, spherical nanodiamonds (SNDs) with a radius of 1 nm (SND1 ) and 1.5 nm (SND2 ), respectively (Fig 1 b). The surface of these spherical NDs was functionalized with 6 and 14 carboxyls, respectively. They were then placed in a 10 × 10 × 10 nm3 water box with neutralizing Na+ and NaCl at a concentration of 0.17 M.

Simulation Protocols Each of the ND slab systems mentioned above was first equilibrated for 30 ns without any electric field applied. Three sets of 50-ns simulations were then performed with different magnitudes of electric field and harmonic restraint force constant,


Page 8 of 34

Page 9 of 34

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


namely, E = 0.8 and k = 0.5, E = 0.8 and k = 1.0, as well as E = 1.6 and k = 1.0, where the units of E and k are kcal mol−1 ˚ A−1 e−1 and kcal mol−1 ˚ A−2 , respectively. A complete list of simulations conducted in this work is provided in Table S1. All simulations were performed with the 2.9 release of NAMD 30 and the CHARMM general force field 31 with TIP3P water. 32 Parametrization of the nanodiamond models follows our previous work. 13 A time step of 2 fs was adopted in all simulations with bonds involving hydrogen atoms constrained using RATTLE 33 and water geometries maintained using SETTLE. 34 The multiple-timestepping algorithm was used with short-range forces calculated every step and long range electrostatics calculated every two steps. The cutoff for short-range nonbonded interactions was set to 12 ˚ A, with a switching distance of 10 ˚ A. Assuming periodic boundary conditions, the Particle Mesh Ewald (PME) method 35 with a grid density of at least 1/˚ A3 was employed for computation of long-range electrostatic forces. Langevin dynamics with a damping coefficient of 1 ps−1 was used to keep the temperature constant at 310 K. Simulations of the ND slabs were performed under constant pressure with a fixed surface area and constant temperature conditions (NPAT), where a Nos´e-Hoover-Langevin piston 36 was used to maintain the pressure in the z-dimension at 1 atm. The parameters COMmotion and zeroMomentum were both set to their default values (no).

Results Validating the Harmonic Restraint Method on Nanodiamonds The harmonic restraint method was developed by Venable, et al. 26 to revise the Lennard-Jones interactions in the CHARMM force field between Na+ and lipid oxygen atoms of carboxyl, phosphate, and ester groups. Prior to applying it on a POPC bilayer, the authors validated their method by comparing the effective charge of sodium acetate measured from the harmonic restraint method with that computed via two other approaches. As nanoparticles are generally much larger and carry more charges than sodium acetate, the first task we undertook



1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

was to validate the harmonic restraint method on our nanodiamond models. To this end, we constructed two spherical NDs (SND1 and SND2 , see Fig 1 b) with a radius of 1 nm or 1.5 nm and a total surface charge of -6 or -14, respectively. We note that these two spherical NDs were only used to demonstrate that the harmonic restraint method can reliably produce Qef f (see below), whereas the Gouy-Chapman theory was not applied here. Although they are not perfectly spherical, the surface charge density of the two SNDs are approximately -0.478 e/nm2 and -0.495 e/nm2 , respectively, which fall within the range of the ND slabs’ σ studied in this work (Table S1). The effective charges of the above two SNDs were found to be -3.86 and -8.52, respectively, according to 50-ns harmonic restraint simulations (Table 1). Following Venable, et al., 26 we also performed the calculation with the drift velocity method and via the EinsteinSmoluchowski relation. In the former method, an electric field was applied and the drift velocity, ν, of a SND was determined from its displacement along the direction of E (Fig S2). As the electric force (EQef f ) is balanced by the friction force (νγ), the effective charge can be determined as Qef f = νγ/E = µγ, where γ is the (unknown) drag coefficient. By performing a series of simulations at different electric field (Fig 2), we determined the electrophoretic mobility µ via a least-square fit. Repeating the calculation for a spherical ND in pure water and assuming that γ remains the same under different ionic concentrations, we determined Qef f from the ratio of the measured µ values. The thus obtained Qef f values (-3.90 and -8.87, respectively) are found to be in good agreement with the results of the harmonic restraint method (Table 1). As an alternative to the harmonic restraint and drift velocity methods, once its electrophoretic mobility is known, the effective charge of a nanodiamond can also be obtained from Einstein-Smoluchowski relation, Qef f = µkB T /D0 , where D0 stands for the diffusion constant of the ND particle in the absence of an electric field. 37 With the calculation of D0 from a 30-ns simulation, Qef f was found to be -4.84 and -7.80 for SND1 and SND2 , respectively. Given that the computation of D0 introduces an additional source of error, we consider


Page 10 of 34

Page 11 of 34

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


the results obtained with the Einstein-Smoluchowski relation in reasonable agreement with those of the other two methods. The drift velocity method and the calculation based on the Einstein-Smoluchowski relation both rely on performing an electrophoresis experiment in silico. Therefore, these two approaches are only suitable for relatively small particles (radius

Estimation of Nanodiamond Surface Charge Density from Zeta

Recommend Documents