Characterization of âStarâ Droplet Morphologies ... - ACS Publications

Article pubs.acs.org/JPCA

Characterization of “Star” Droplet Morphologies Induced by Charged Macromolecules Mahmoud Sharawy† and Styliani Consta*,†,‡ †

Department of Chemistry, The University of Western Ontario, London, Ontario Canada N6A 5B7 Department of Chemistry, University of Cambridge, Lensfield Road, Cambridge, CB2 1EW, United Kingdom

‡

S Supporting Information *

ABSTRACT: “Star” morphologies of charged liquid droplets are distinct droplet conformations that, for a certain charge squared to volume ratio, have lower energy than their spherically shaped analogues. For these shapes to appear, the charge should be carried by a single ionic species. A typical example of a charge carrier that we employ in this study is a fully charged double-stranded oligodeoxynucleotide (dsDNA) in an aqueous and an acetonitrile droplet. We characterize the structure and dynamics of the star-shaped droplets. We find that by increasing the charge squared to volume ratio, the droplet evolves from spherical to “spiky” shapes, by first passing from droplet sizes that undergo enhanced shape fluctuations relative to those of the larger spherical droplets. These fluctuations mark the onset of the instability. We also find that in the spiky droplet, the orientation of the solvent molecules in the first shell about the dsDNA is very close to that in the bulk solution. However, this orientation is substantially different farther away from the dsDNA. With regards to dynamics, the motion of the spikes is reflected in the autocorrelation functions of rotationally invariant order parameters that show a damped oscillator form of decay, indicative of the elastic motion of the spikes. We compare the formation of spikes with that of the ferrofluids and the dielectric materials in an electric field, and we conclude that they represent a different entity that deserves its own characterization. The study provides insight into the manner in which the charge distribution may give rise to well-controlled droplet morphologies and calls for experiments in this direction.

■

INTRODUCTION The instabilities of droplets induced by an electric field have received continuous attention starting from the seminal works of Rayleigh,1 followed by those of Taylor.2,3 The stability of charged droplets plays significant role in many practical applications such as electrospray ionization (ESI) coupled to mass spectrometry (MS),4−6 atmospheric aerosols,7 liquid helium droplets charged with surface-bound electrons,8,9 and jets.10,11 The simplest model that describes the charge induced instabilities was developed by Rayleigh in 1882.1 Rayleigh’s model considers a spherical droplet that is uniformly charged and may undergo shape fluctuations. The model is based on the assumptions that the droplet is a conductor and that its volume is constant (the surface area may change). In this model, the energy of the fluctuating droplet is written as the sum of the electrostatic energy and the surface energy. Linear stability analysis of the energy12,13 leads to the Rayleigh criterion for stability, which reads as Q R 2 = 64π 2ε0γR 0 3

spherical shape. eq 1 defines the Rayleigh limit of a charged droplet. When the surface interactions overcome the electrostatic repulsion between the charges, then a droplet is found “below the Rayleigh limit”. When the strength of the interactions is reversed then a droplet is found “above the Rayleigh limit”. Rayleigh’s model has been used in the stability studies of atmospheric aerosols and electrospray produced droplets among other systems (charged colloids,14,15 atomic nuclei). In previous research,16−21 we have established by molecular simulations and analytical theory that different charge distributions lead to a variety of droplet morphologies above or below the Rayleigh limit. When separated ions are present in the droplet, the instability manifests by fragmentation of the droplet and the release of the excess charge by small clusters of solvated ions.16,17,22 In the presence of macroions with high bound charge such as peptides or proteins,18 nucleic acids,23 carbon buckyballs,24 and charged poly(ethylene glycol),19,25 the charge distribution in the droplet is different from that of relatively fast-diffusing single ions. So far we have classified three general droplet morphologies,18,20 which are, the “popping” model of a contiguous extrusion of the macroion

(1)

where QR denotes the charge of the droplet at the Rayleigh limit, ε0 the vacuum permittivity coefficient, γ the surface tension, and R0 the droplet radius. eq 1 comes from the l = 2 spherical harmonic mode, which is the first mode that may become unstable in the energy perturbation to that of the © 2016 American Chemical Society

Received: August 22, 2016 Revised: October 10, 2016 Published: October 19, 2016 8871

DOI: 10.1021/acs.jpca.6b08486 J. Phys. Chem. A 2016, 120, 8871−8880

Article

The Journal of Physical Chemistry A

Figure 1. Typical snapshots from molecular simulations18,23 of droplets containing a macroion and complexes of macroions. (a)Completely protonated polyhistidine (blue macroion) in an aqueous droplet.18 The water molecules are depicted by the red (oxygen) and white (hydrogen) spheres. The size of the polyhistidine has been enlarged for visualization purposes. (b) [dsDNA]20− in an aqueous droplet of approximately 1700 H2O molecules.23 The droplet is found above the Rayleigh limit. The droplets are characterized by the formation of solvent spikes on the surface. (c) [dsDNA]20− surrounded by Na+ ions (blue spheres) and Cl−ions (green spheres) in an aqueous droplet of approximately 1500 H2O molecules.23 The droplet is found above the Rayleigh limit. Solvent spikes appear on the surface of the droplet. (d−f) Demonstration that the formation of spikes on the surface of a droplet is a general behavior that appears when a nonfissile charged species is solvated in the interior of a droplet. The droplets are composed of 2148 TIP3P modeled water and a spherical fictitious ion positively charged at the center of the droplet. The charge of the central ion is (d) 18 e (where e is the elementary proton charge), (e) 27 e, and (f) 36 e.18,21

from a droplet,18−20 the “pearl” model of subdroplets along a macroion,18,20 each with charge as predicted by the Rayleigh limit, and the “spiky” droplet morphology.18 We have found that the contiguous extrusion of the macroion from a droplet takes place below the Rayleigh limit.19,21 The spiky droplet morphology appears above the Rayleigh limit when a macroion carries a nonfissile charge. In Figure 1 we show examples of droplets that attain the star morphologies for a certain charge squared to volume ratio. Parts a−c of Figure 1 depict realistically (all-atom) modeled peptides18 and nucleic acids23 in aqueous droplets. In parts d−f of Figure 1, the charge carrier is a spherical ion found at the center of an aqueous droplet.18,21 Parts d−f of Figure 1 demonstrate that the formation of spikes is a general phenomenon that has its origin in the nonfissility of the charge. We note here that ordered star structures have been found in different types of systems varying in size from the atomistic scale18 to self-assembled structures as described in ref 26. Research is being conducted for finding the underlying principles for these star formations.26 Even though we have found the “spiky” droplet morphologies in previous research18 a detailed characterization of their structure and dynamics is still missing from the literature. Since the “star” droplet structures are stable they can be created experimentally therefore, knowledge of their structure is important for experimental realization. In our studies we have chosen a fully charged doublestranded oligodeoxynucleotide (dsDNA) in aqueous and acetonitrile droplets as a representative realistic system to demonstrate the structural and dynamic features of the “spiky” droplet morphologies. The reasons that we choose to study the DNA complexes have been presented in ref 23. Both water and

acetonitrile solvents are commonly used in experiments. Even though both solvents are polar they have considerably different physical properties that are summarized in Table 1. Table 1. Parameters of the Solventsa solvent

d (g/cm3)

dN

γ

ϵr

μ (D)

BP (K)

H2O CH3CN

1.00 0.786

33.46 11.53

73 29

80 36

1.9 3.9

373 350

γ denotes the experimental value of the surface tension (mN/m) at room temperature; dN, the number density (number of molecules/ nm3); μ, the dipole moment (gas molecule); and ϵr, the dielectric constant (relative permittivity). BP is the boiling point. a

We find that the solvent structure around a macroion in the spiky droplets is distinct relative to that found in the bulk solution27 or in the bulk interfaces.28,29 Our findings regarding the structure of the charged droplets are general and hold for other possible candidates of high charge as well such as peptides and carbon buckyballs.24

■

MODEL AND COMPUTATIONAL METHODS Model. We simulated droplets that were composed of a double-stranded DNA in water and in acetonitrile. The dsDNA that was used in this study was d(A11)2 composed only of adenine (A) and thymine (T), with a net charge of −20 e ([dsDNA]20−). The 5′ ends lack a phosphodiester linkage at the 5′ positions and are terminated by a hydroxyl group. The helical structure was initially in B-tract conformation and made one turn. 8872


Article


nanodrops containing the d(A11 )2 . In these tests, we thermalized one from each system duplicate to T = 350 K for 2 ns, decoupled it from the thermostat, and then we continued the simulation without temperature control. The dynamics was compared with the second duplicate that was simulated at constant temperature until the d(A11)2 was in the gas phase. The one difference in the dynamics between the two set-ups, without and with temperature control was the cooling of the droplet without the temperature control. This had no profound effect on the results we report here. We decided that the best method to adopt in this investigation was to use the Nosé−Hoover thermostat with all runs. To test the sensitivity of the results on the MD simulation package, a simulation was also performed using NAMD37 version 2.9 molecular dynamics program. The program was used to simulate a 4000 H2O nanodrop that contained the d(A11)2. This simulation was under similar conditions to those simulated with GROMACS, that is to say constant temperature at T = 350 K, AMBER force field (ff12SB parameters), and nonperiodic boundary conditions. The results were typical of those found in the droplets that were simulated with GROMACS using AMBER99SB-ILDN. Moreover, a single-strand G-quadruplex d(A3G3T2AG3T2AG3T2AG3A2) (PDB ID 2HY9) was simulated in aqueous TIP3P droplets. Hereafter, we use GQX to denote the G-quadruplex. The GQX has a 25 phosphodiester linkages that give a total charge of −25 e. The simulation protocol for GQX was the same as that of the [dsDNA]20−. For comparisons in the structure of the solvent, bulk simulations of [dsDNA]20− were performed in 11139 TIP3P modeled water molecules at T = 350 K. The electorstatic interactions were treated with particle mesh Ewald with Coulomb cutoff at 4 nm. We added a counter charge of Na+ ions to neutralize the net charge on the dsDNA. The Na+ ions were restrained beyond the cutoff distance: at positions away from the macroion. The equations of motion were integrated with the velocity-Verlet the integrator and the temperature was controlled with the Nosé−Hoover thermostat. The simulation box dimensions were ≈7 × 7 × 7 nm. The system was equilibrated for 0.5 ns and the production run was for 10 ns.

The nucleic acids were solvated by water molecules modeled by the transferable intermolecular potential function (TIP3P)30 with flexible bonds and angle bending. The acetonitrile parameters are not readily available for AMBER force field. Thus, we implemented the parameters31 of a six-site acetonitrile into the AMBER port in GROMACS. To validate the proper implementation of these parameters comparisons of the radial distribution function (RDF) in the bulk solution are presented in Figure S1 in the Supporting Information. We note here that the parametrization simulations in ref 31 were performed at temperature 298 K and pressure of 1 bar with isothermal compressibility taken from the experiments. In this model, the partial charges were generated by quantum chemistry calculations and not by trying to fit certain preselected solution properties. Once the partial charges were determined, the force field was tested for its quality by reproducing certain bulk solution properties. It was found that their force field can reproduce reasonably well the experimental density, very well the heat of vaporization and well the solution RDFs and the dipole moment of the acetonitrile molecules. It has been shown in our previous articles18 and it is confirmed in the present study that the formation of the spikes is a general behavior that does not depend on the fine details of the force field. For instance, both TIP3P or SPC/E water models will show the formation of the spikes. We use the acetonitrile model to point out the differences from water due to the difference in the dielectric constant between the two liquids. Computational Methods. Molecular dynamics (MD) simulations were carried out with the GROMACS32−34 package version 4.6.5 using AMBER99SB-ILDN35 force field port. This force field includes optimized potential terms for bondstretching and angle bending vibrations, which are described by an ideal harmonic oscillator, as well as torsional strains, Coulombic electrostatic interactions, and Lennard-Jones (LJ) interactions. The cross interactions between two sites i and j are estimated by the combining rules ϵij = (ϵiϵj)1/2 and σij = (σiσj)1/2, where ϵi and σi denote the depth of the LJ potential and the diameter of the atomic species i, respectively. The following protocol was applied for droplets: the d(A11)2 initially in vacuo, was solvated with TIP3P H2O or CH3CN shells; each solvation step was repeated with a different seed number so as to generate a duplicate of different solvent positions. The solvation was followed by energy minimization; and 200 ps of position restraint MD, which allowed the macroions to migrate to the most energetically favorable positions. After these preparation steps, the production simulations were carried out in a series of 0.5 ns constant temperature runs. Because the evaporated solvent molecules that were far from the droplet had negligible interactions with the solute, the evaporated solvent molecules were removed every 0.5 ns when they were located outside a radius of 70 nm from the solute. The Nosé−Hoover thermostat36 was applied at a temperature of 350 and 320 K to the H2O and CH3CN droplets, respectively. No bond constraints were applied. The time step used in all the simulations was 0.5 fs, which is 13 times smaller than the fastest, O−H, vibration in the system. The electrostatic interactions were calculated without cutoff in nonperiodic boundary conditions (in vacuo). The simulations were terminated when the two strands of the d(A11)2 were separated. In the most of the realizations the total time length of the runs ranged between 5 to 18 ns. The effect of the Nosé−Hoover algorithm on the dynamics of the simulation was tested on 6500 water molecules

■

RESULTS AND DISCUSSION Structure of Charged Droplets. The droplet evolution during the desolvation of the dsDNA has several distinct features: The first prominent feature is the change in the droplet morphology from spherical, when the system is below the Rayleigh limit, to spiky when the system is above the Rayleigh limit. The second prominent feature is the strong orientation of the solvent molecules dipole moments around the macroion. These two features are discussed in the following sections. Description of the Droplet Morphologies. Typical snapshots of the changes in the shape of a droplet containing a [dsDNA]20− ion in the course of evaporation are shown in Figure 2, parts a and b, for water and acetonitrile, respectively. These changes are accompanied by the unwinding and the separation of the [dsDNA]20−. A movie that shows the “on-thefly” development of the spiky droplet morphologies and the motion of the spikes is included in the Supporting Information (Movie S1). The changes in the droplet shape in relation to the conformational changes of [dsDNA]20− have been presented in our previous article.23 In this study, we focus on quantifying the “spiky” droplet morphologies and explaining their physical 8873


Article


below the Rayleigh limit. This critical charge was estimated by eq 1 with the surface tension value of 0.0432 N/m at T = 350 K for TIP3P water model.38,39 In comparison to the aqueous droplet, a droplet composed of 11500 CH3CN molecules and a [dsDNA]20− ion (Figure 2b(i) has a radius of ≈6.1 nm. Using the Rayleigh criterion eq 1 with a surface tension of 0.029 N/m (experimental value at T = 320 K), the critical charge is estimated to be −12 e; hence, the CH3CN droplet is also found initially below the Rayleigh limit. At a certain charge squared to volume ratio the droplets enter the instability regime. In our simulation setup the charge to volume ratio is reached by gradual evaporation of individual water molecules. The onset of the instability is marked by enhanced fluctuations that are large but short-lived relative to the spikes shown in Figure 2. The transition to the “spiky” shapes is not sharp and we attribute this behavior to the finite size of the droplet. It is noted that it is surprising and remarkable that a droplet composed of a dielectric with no free charges shows the instability at a size close to the theoretical prediction. For γ = 0.0432 N/m at T = 350 K and for a charge of −20 e, the theoretical prediction for the onset of instability is at 6000 H2O molecules. The simulation finding is for 5500 H2O molecules. To this number we add 400 H2O molecules equivalent to the size of the dsDNA. Thus, we observed the onset of instability at 5900 H2O molecules. To test the Rayleigh prediction we have also studied different systems such as that of an aqueous droplet solvating a G-Quadruplex. The quadruplex that we have chosen to simulate has a net charge of −25 e, with estimated Rayleigh limit at 9300 H2O molecules. The enhanced fluctuations appear for droplets with approximately 9500−9300 H2O molecules. Figure 2a,b(ii−iv) for both H2O and CH3CN droplets show the droplet morphologies above the Rayleigh limit, respectively. In this regime, the droplet shape deforms. In aqueous droplets, the deformation starts by a forming a droplet with a single spike at approximately 4820 H2O molecules (if we consider that the volume of the dsDNA is approximately that of 400 H2O molecules, the first spike appears at approximately 5200 H2O molecules). The shape deformations as the droplet shrinks appear with a gradual increase in the number of spikes as shown in Figure 2a(ii−iv). In contrast, Figure 2b shows that in the CH3CN droplets the deformations are not as pronounced as that of water in a shrinking droplet. The electric field orients the CH3CN dipoles in a head-to-tail arrangement. Acetonitrile has a lower dielectric constant than that of water; therefore, it is expected that the polarization of the solvent, due to orientation, is less pronounced than that of water. The spiky droplet morphologies are stable shapes with lower energy than that of the corresponding spherical droplet. This is evidenced by the fact that a particular shape is maintained for a range of droplet sizes. For instance, in the dsDNA−H2O droplet 1-, 2-, 3-, 4-, 5-, and 6-spike morphologies are maintained in a range of several hundreds of water molecules at approximately 4215, 3794, 2661, 1809, 1278, and 700 H2O molecules, respectively. The spikes are dynamic, they coalesce and divide to form larger or smaller extrusions, respectively. When the spikes divide, they separate in a zipping-down manner starting from the tip. We note that we find instances where the droplet morphology oscillates between two different number of spikes, three-four or four-five until it settles to a permanent shape. Regarding the number of the spikes we compare the data reported in ref.21 for an aqueous droplet that has a spherical

Figure 2. Typical snapshots in the course of evaporation of (a) an aqueous droplet with an embedded [dsDNA]20− at T = 350 K and (b) an CH3CN droplet with a [dsDNA]20− at T = 320 K. The nucleotides are colored in red for A (adenine) and blue for T (thymine).

origin. At first, we summarize the evolution of the droplet shape with respect to the number of molecules in the system (the charge remains constant). The initial state of the aqueous nanodrop (Figure 2a(i), which has a radius ≈3.6 nm, shows a clear stability demonstrated by its spherical shape. Rayleigh’s critical charge for this droplet size is −22 e, therefore the nanodrop is found 8874


Article


spikes and 12 spikes in droplets with a spherical charge in the middle (Figure 1d−f). Typical geometries are that four spikes form a tetrahedron, five a tetragonal pyramid and trigonal bipyramid, six an octahedron, and 12 an icosahedron. These structures show striking similarities to those that electrons form in the Thomson problem. In the Thomson problem, M electrons constrained on a unit sphere attain a certain geometry by minimizing their electrostatic energy. The Thomson problem and the formation of spikes are fundamentally different problems, but we attribute the similarities to the fact that the spikes arrange themselves so as their distance between one another is maximized. For instance in aqueous droplets the spikes at the termini of the dsDNA are found in staggered conformations (Figure S2 in the Supporting Information) in order to maximize their distance. We note here the differences between the two problems. In the MD simulations of the spiky droplets there is temperature, thus the number of spikes and their length are variables. In the course of the MD runs we find events of birth-death of spikes as well as spikes of different length. Even though there are fluctuations in the number of spikes and their length, we can still clearly identify the major structures and their similarity with those in Thomson problem. For certain number of spikes MD may detect more than one structures. The different structures appear because of the fluctuations in the positions of the spikes. An example is shown in Figure 4 where the droplets with five rays can form a

positive charge at the center to that of the current calculations for the aqueous [dsDNA]20− droplets. This comparison is presented in Figure 3. The simulations discussed in ref 21 were

Figure 3. Number of spikes vs √X. The open circles correspond to aqueous droplets with a positive spherical ion at the center21 and the red squares to aqueous droplets with [dsDNA]20−.

performed for droplets composed of 2148 TIP3P modeled water molecules at T = 300 K and surface tension approximately equal to 0.0523 N/m. In those simulations the number of H2O molecules remained constant and the charge varied. The droplet was found in a spherical cavity so that an effective vapor pressure was established. In the simulations presented here, the charge is the droplet remains constant and the number of H2O molecules varies. The current simulations are performed at T = 350 K (γ ≈ 0.0432 N/m). To compare the number of spikes, we plot this number as a function of the square root of X, where X is defined as X=

Q2 64π 2ε0γR 0 3

(2)

Evidently, at the Rayleigh limit, X = 1. Alternatively, X = (Q/ QR)2 when the number of molecules remains the same in the droplets as in the case of the single ion−water simulations in a spherical cavity discussed above. Moreover, X = NR/N, (where N denotes the number of solvents molecules and NR those at the Rayleigh limit) when the charge remains the same in the droplet and the number of solvent molecules changes as in the case of the solvent evaporation in the [dsDNA]20−-water droplets presented here. In writing X = (Q/QR)2 or X = NR/N, we assume that the surface tension remains the same in the droplets with various charge squared to volume ratios. Figure 3 reveals that the [dsDNA]20− droplets have the same number of spikes as the positively charged droplet for the same value of X1/2. This is an impressive finding because it shows that even though the structure of the spikes in the negatively and positively charged systems is different, their number is determined by the amount of charge only. Even more, the complicated structure of the [dsDNA]20− vs that of a simple spherical ion also does not affect the number of spikes. Therefore, the formation of spikes is a general finding that does not depend on the molecular details of the macroion. As the dsDNA containing droplet evaporates, the clear structure of the thicker spikes is lost because of the many string-like formations of the solvent. For this reason for the dsDNA-H2O droplet we compared the data in Figure 3 up to a maximum of six spikes. We note here that the spikes attain certain geometries. By visual inspection we could distinguish the geometries up to nine

Figure 4. Typical snapshots of droplets containing a central ion with charge Q = 24 e and 2148 TIP3P modeled water molecules21 at T = 300 K. The rays can form a tetragonal pyramid (top panel) and a trigonal bipyramid (bottom panel).

tetragonal pyramid and a trigonal bipyramid. We cannot minimize the energy of the spiky droplet structures because the lowering of the temperature will increase the surface tension and thus, the dominant interactions in the system will be modified. 8875


Article


two different droplet shapes (with the same number of solvent molecules). Molecular Structure of the “Star” Morphologies. A closer look on the structure of water is presented in Figure 7.

The global shape of the droplet can be systematically analyzed by using rotational invariant descriptors. A rotational invariant descriptor is given by ⎛ 4π Q L = ⎜⎜ ⎝ 2L + 1

⎞ |Q LM|2 ⎟⎟ ⎠ M =−L M=L

∑

(3)

where N

Q LM =

∑a = 1 YLM(rcm ̂ (a)) N

(4)

where YLM (r̂cm (a)) are the spherical harmonics and r̂cm is the unit vector defined from the center of mass of the droplet to atom a. Figure 5 shows how the various QL distinguish the spiky shapes. As shown, the largest differentiation appears by using

Figure 5. QL vs L for various droplet morphologies.

the Q3 and Q6 descriptors. The combination of two order parameters (Figure 6) can better distinguish the droplet morphologies. An example of this combination is shown in Figure 6 for Q3 and Q6. Other combinations may also be used such as Q2 and Q6. A combination of order parameters can be used in finding the free energy barrier of transitions between

Figure 7. Orientational probability distribution F(θ, r) of solvent dipole moments about the center of mass of DNA. The profiles corresponds to (a) H2O about d(A11)2 in a droplet of 1680 H2O molecules (66 nm3) as shown in Figure 2a, (b) H2O about d(A11)2 in bulk, and (c) CH3CN about d(A11)2 in a droplet of 1200 CH3CN (137 nm3). The calculations were averaged over 500 ps.

Figure 7 shows the probability density of the cosine of the angle (θ) defined between the dipole moment vector of H2O (or CH3CN molecule) (denoted hereafter by μ⃗s) and the vector pointing from the center of mass (COM) of dsDNA to the oxygen site of a water molecule (or the methyl group of CH3CN). In these calculations the solvent surrounding the nucleic acid was divided into 0.5 nm wide shells. Figure 7a shows that at a distance greater than 3 nm from the COM, the μ⃗s are colinear with the axis from the COM of the dsDNA to the oxygen site of a water molecule. The H2O molecules in the shells between 3 to 6 nm are found in the spikes as depicted in Figure 8. Figure 8 shows only the μs orientation of the 180°. Thus, the molecules in the spikes “see” the dsDNA like a point

Figure 6. Distinction of the droplet morphologies using a combination of two order parameters. 8876


Article


Frequently, the water molecules at the spikes of the charged droplets easily form ring structures40 or other organized shapes by local motion. Typical shapes of clusters at the tips of the spikes are shown in Figure 9. These clusters are reminiscent of

Figure 8. Radial distribution of the 180° dipole moment orientation of H2O about the [dsDNA]20− in a typical “spiky” droplet. The droplet comprised 1680 H2O molecules and d(A11)2. The boundary of the spiky morphology is shown by schematic. The distribution was calculations by averaging the H2O molecules orientation for 500 ps. The white region at the center corresponds to the space occupied by the [dsDNA]20−.

charge. In the range of 1.03.0 nm, the distribution of θ is broad and peaks at 150−180° (see also Figure 8). At this distance, the H2O molecules are found in the closest solvation shells to the [dsDNA]20−, thus they sense the electric field by the local charge of dsDNA. Therefore, water dipoles are strongly oriented toward specific charged regions of phosphates groups as they surround the dsDNA. In summary, the water close to the dsDNA “sees” the duplex as polyions, while the water at far distance from the dsDNA “sees” the DNA as a point charge. Even though the H2O molecules are strongly oriented in the spikes, they have freedom in their motion. For comparison purposes, we calculated the distribution of the angle θ around the [dsDNA]20− in bulk aqueous solution (Figure 7b). In distance greater than 3 nm the H2O dipole moments attain all possible orientations almost uniformly. Here, there is no surface. The presence of additional water layers, quenches the strong orientation of the 180° degrees that appears in the droplets. In the distance less than 3.0 nm, the angle distribution in bulk shows a similar frequency (1.4) to that in the spiky droplet (Figure 7a). However, in the spiky droplet H2O molecules near the backbone of dsDNA are on the surface. This configuration results in overlapping of highly oriented H2O dipoles toward the COM with H2O dipoles in near shells that are oriented toward the local charges. We compare the structure of H2O with that of CH3CN at the same ratio of Q2/γR03 ≈ 1600, that corresponds approximately to 2.7 times of the Rayleigh limit ratio for both solvents. In CH3CN droplets, the orientation of their dipole moments around the dsDNA is not as pronounced as that of the H2O molecules. At a distance greater than 3.5 nm the 180° H2O dipole moment orientation still prevails. At shorter distances the angle distribution is almost uniform. The lack of strong orientation of CH3CN in the electric field of the dsDNA reflects the lower dielectric constant of CH3CN than that of H2O. The orientation of CH3CN will depend on the temperature. We expect that at temperature below 320 K where the thermal motion is quenched, or at a higher charge, the electric field of the dsDNA will have a stronger effect in orienting the CH3CN molecules.

Figure 9. Typical cluster formations at the tips of the spikes. In the snapshots, the lengths of the H-bonds are indicated.

the ring, book, and prism structures that have been found for small water clusters in the gaseous state.27,41−44 Short strings of water with approximately 3−4 H2O molecules and rings were also identified on the surface of pristine aqueous droplets. These clusters are more subtle in comparison to those that appear in a charged droplet. In the pristine aqueous droplets the strings are shorter than those of the charged droplets and the rings are short-lived. In the aqueous droplet we compared the length of the hydrogen bonds and the oxygen−oxygen radial distribution functions in droplets of different sizes as well as in various regions of droplets (spikes vs central part of the droplet). We also compared these quantities with that of a pure, zero net charge, aqueous droplet. The comparisons did not show any difference of these quantities in the various systems for the TIP3P model of water. This could be a result of the limitations of the TIP3P model in the presence of highly charged ions. Dynamics. Figure 10 shows the dipole moment (μ⃗s) autocorrelation function (ACF) of the water molecules surrounding the dsDNA. In general the smaller the droplet the slower the ACF decays. The dipole ACFs for a spherical, 3spike and 4-spike droplet is fitted by a sum of three exponential functions (the parameters of the fitting are presented in the Supporting Information). As it has been already discussed in the previous sections (Figure 7 and Figure 8), the H2O molecules are strongly oriented close to the dsDNA and in the spike region. Therefore, the H2O dipole moment decorrelates by diffusion. There are several mechanisms for decorrelation. The fastest relaxation time varies between 1.8 ps (spherical droplet) to 1.2 ps for the three-spike shape. We attribute the short time in the local motion of the water dipoles. This time is similar for the different shapes. There is an intermediate time scale in the range of 5 ps (three-spike droplet) to 25 ps (four8877


Article


Figure 10. Dipole moment autocorrelation function (ACF) of H2O (C(t)) in a droplet of 6500 H2O molecules, which is spherical (solid line), 2827 H2O molecules, which attains a three-spike shape (dotted line), and a 1678 H2O molecules, which attains a four-spike shape as shown in Figure 8 (dashed line).

spike droplet). We attribute the intermediate time to the decorrelation of the H2O dipole moment by the collective motions in the spikes and the fast diffusion of the molecules within the spikes. Finally, the longest time, which is a few hundreds of picoseconds is attributed to several motions. The rotation of the droplet is one source of the slow relaxation time. In our systems the occasional removal of the evaporated water molecules masks the real rotational time. Monitoring of the motions of the molecules in the spikes vs the core, showed that the molecules exchange within a few hundreds of picoseconds (Movie S2 in the Supporting Information) between the two regions. Indicative trajectories of the diffusion of a single molecule is shown in Figure S3 and Movie S3 of the Supporting Information. The diffusion of the molecules will be affected by the motions of the spikes. These motions include floppy movements of the spikes, which are collective motions of the molecules, as well as the quenching of the spikes and the birth of new ones. We expect that the diffusion time will vary considerably for the various morphologies. For instance for certain ratios of droplet charge squared to volume a droplet may appear with a definite number of spikes. These are the systems of which we show the ACFs in Figure 10. However, there are going to be droplet charge and volume parameters where the droplet will not settle to a definite number of spikes (the number of spikes may vary; for the thicker spikes between 1 and 6 the number of spikes may vary at most by one). In Figure 11, the time ACF for the Q3 descriptor is shown with a fitting to the function

Figure 11. Time autocorrelation of Q3 for droplets with one and two spikes (from top to bottom), respectively. The ACF for 3 spikes is presented in Figure S4 of the Supporting Information. Details are discussed in the text.

position. The ACF shows a decaying oscillatory behavior, which is indicative of the elastic motion of the spikes.

■

CONCLUSIONS The formation of spikes in droplets due to the presence of a nonfissile charged species such as a protein or nucleic acid has been unexpected. These structures appear beyond a charge squared to volume ratio of the droplet that marks the onset of instability. The spiky structures are dynamic and persist over long times limited only by the rate of solvent evaporation. We characterized their shapes by rotationally invariant order parameters and studied their dynamics. Analysis of the orientation of the water molecules around a fully charged dsDNA shows that close to the DNA the orientation of the water molecules is determined by the local structure of the DNA. Beyond a certain distance, the H2O molecules orient around the central macroion as they would orient around a point charge. Autocorrelation functions of the order parameters reveal the elasticity of the system through the decaying oscillations of the structural quantifiers of the spikes. Destabilization of liquid interfaces in the presence of magnetic or electric fields has some commonalities. In every case the surface tension acts as a stabilizing factor, while electroor magneto-static forces tend to disrupt the interface. Spherical or planar interfaces are, in the absence of the fluctuations, the extrema of the system energy, which is written as a sum of electrostatic, surface energy and gravitational energies. The Rayleigh criterion is derived from an analysis of the interfaces stability due to small perturbations, leading to eq 1 applicable to a variety of models. As the field increases, a planar or a spherical interface must attain a new configuration that minimizes the system energy.

y(x) = exp( −x /a0)(cos(x /a1) − (a1/a0) sin(x /a1))

with parameters of (a0, a1) = (181.177, 36.3714), (124.901, 37.834) in units of ps for one and two spikes, respectively. For three spikes, the ACF is shown in Figure S4 of the Supporting Information. The ACF shows behavior similar to that of a damped harmonic oscillator. This is supported by visual inspection of the trajectories of droplets with 1, 2, and 3 spikes. In the simulations we observe persistent structures with fixed number of spikes moving about an apparent equilibrium 8878



■

To our knowledge, so far there are no analytical results describing the new minimal energy surface in the general case. Taylor’s results demonstrate that the only self-similar solution is a unique cone with an angle of ≈100 deg for an an aqueous droplet at an electric field.3 The conical solution also predetermines the field around the cone tip that may of may not be compatible with that of an original problem of minimization of the energy of an interface. The morphology of liquid surfaces in the presence of electric field beyond stability transition includes jets, droplet spines and patterns in ferrofluids.45−47 Review of the existing results leads the authors to conclude that there is no compact surface that minimizes the free energy of a surface that includes only the surface tension and electro- or magneto-static energies. The main evidence for this comes from the analysis of the numerical simulation results as analytical results are not available. In all the cases mentioned above additional factors are required to stabilize the interface. In the case of jets, the interface is not stationary, but it is a steady dynamic state stabilized by viscous flow in the jet. Patterns in ferrofluids liquids are determined by the complex nonlinear magnetization of the liquid in the applied magnetic field. It is interesting to note that ferrofluids do not admit Taylor conical solution. At the tip of a cone, the balance of forces requires that the magnetic force on the interface is to be compensated by the surface tension force. On the other hand, the surface tension force diverges as one approaches the tip while the magnetostatic force reaches its maximum value due to saturation of the magnetization. Thus, the interfaces in ferrofluids must be always smooth surfaces. This in turn validates the use of γ × area surface energy term. The patterns observed in ferrofluids should be contrasted with the shapes observed in charged droplets. The presence of sharp tips and ordered structure of solvent is evident. Because of the differences of these structures from the ferrofluids and other known sharp protrusion (jets) we think that a description of these minimum energy shapes by the analysis of only the electrostatic and surface terms would be incomplete and likely to fail. The shapes of the “star” morphologies are well controlled by the amount of charge in the central ion. The presented study provides the background analysis of their structure and dynamics that it is useful for their experimental realization. Potential applications of the “star” droplet morphologies may be envisioned in material science and catalysis.

■

Article

AUTHOR INFORMATION

Corresponding Author

*(S.C.) E-mail: [email protected]. Notes

The authors declare no competing financial interest.

■

ACKNOWLEDGMENTS

■

REFERENCES

S.C. acknowledges an International Incoming Marie Curie Fellowship of the European Commission held at the Department of Chemistry, University of Cambridge. Sci-Net, SHARCNET and Compute Canada are acknowledged for providing the computing facilities. S.C. thanks Prof. D. Frenkel, Department of of Chemistry, University of Cambridge, for his infinite curiosity about physics, the inspiring discussions on the charged systems, and comments and Dr. S. S. Xantheas, Pacific Northwest National Laboratory, USA, for his comments and encouragement. Discussions with Dr. Anatoly Malevanets are very greatly appreciated.

(1) Rayleigh, Lord XX. On the equilibrium of liquid conducting masses charged with electricity. Philos. Mag. Ser. 5 1882, 14, 184−186. (2) Taylor, G. Disintegration of water drops in an electric field. Proc. R. Soc. London, Ser. A 1964, 280, 383−397. (3) Taylor, G.; McEwan, A. The stability of a horizontal fluid interface in a vertical electric field. J. Fluid Mech. 1965, 22, 1−15. (4) Dole, M.; Mack, L.; Hines, R.; Mobley, R.; Ferguson, L.; Alice, M. Molecular meams of macroions. J. Chem. Phys. 1968, 49, 2240−2249. (5) Iribarne, J.; Thomson, B. Evaporation of small ions from charged droplets. J. Chem. Phys. 1976, 64, 2287−2294. (6) Thomson, B.; Iribarne, J. Field-induced ion evaporation from liquid surfaces at atmospheric-pressure. J. Chem. Phys. 1979, 71, 4451− 4463. (7) Davis, E.; Bridges, M. The Rayleigh limit of charge revisited: light scattering from exploding droplets. J. Aerosol Sci. 1994, 25, 1179− 1199. (8) Leiderer, P. Nonlinear Effects - The Multielectron Dimple; Springer: 1997. (9) Leiderer, P.; Scheer, E.; Kono, K.; Lin, J.-J.; Rees, D. Stability of Surface State Electrons on Helium Films. J. Low Temp. Phys. 2016, 183, 258. (10) Eggers, J.; Villermaux, E. Physics of liquid jets. Rep. Prog. Phys. 2008, 71, 036601. (11) Guerrero, J.; Rivero, J.; Gundabala, V. R.; Perez-Saborid, M.; Fernandez-Nieves, A. Whipping of electrified liquid jets. Proc. Natl. Acad. Sci. U. S. A. 2014, 111, 13763−13767. (12) Hendricks, C.; Schneider, J. Stability of a conducting droplet under the influence of surface tension and electrostatic forces. Am. J. Phys. 1963, 31, 450−453. (13) Consta, S.; Malevanets, A. Disintegration mechanisms of charged nanodroplets: novel systems for applying methods of activated processes. Mol. Simul. 2015, 41, 73−85. (14) Deserno, M. Rayleigh instability of charged droplets in the presence of counterions. Eur. Phys. J. E: Soft Matter Biol. Phys. 2001, 6, 163−168. (15) Thaokar, R.; Deshmukh, S. Rayleigh instability of charged drops and vesicles in the presence of counterions. Phys. Fluids 2010, 22, 034107. (16) Consta, S. J. Mol. Struct.: THEOCHEM 2002, 591, 131−140. (17) Consta, S.; Mainer, K. R.; Novak, W. Fragmentation mechanisms of aqueous clusters charged with ions. J. Chem. Phys. 2003, 119, 10125−10132. (18) Consta, S. Manifestation of Rayleigh instability in droplets containing multiply charged macroions. J. Phys. Chem. B 2010, 114, 5263−5268.

ASSOCIATED CONTENT

S Supporting Information *

The Supporting Information is available free of charge on the ACS Publications website at DOI: 10.1021/acs.jpca.6b08486. (a) Comparison of radial distribution functions to validate the implementation of the acetonitrile parameters in GROMACS, (b) examples of the geometries of aqueous droplets containing a [dsDNA]20−, and (c) fitting parameters of the order parameters ACFs (PDF) Movie S1, formation of spikes during the desolvation of the 11-mer all A·T dsDNA in a −20 e charged aqueous droplet (MPG) Movie S2, diffusion of water molecules between the inner and outer (spike) solvation shells of dsDNA (MPG) Movie S3, instance of a diffusion path of a water molecule MPG 8879


Article

The Journal of Physical Chemistry A (19) Consta, S.; Malevanets, A. Manifestations of charge induced instability in droplets effected by charged macromolecules. Phys. Rev. Lett. 2012, 109, 148301. (20) Consta, S.; Malevanets, A. Classification of the ejection mechanisms of charged macromolecules from liquid droplets. J. Chem. Phys. 2013, 138, 044314. (21) Consta, S.; Oh, M. I.; Malevanets, A. New mechanisms of macroion-induced disintegration of charged droplets. Chem. Phys. Lett. 2016, 663, 1. (22) Ichiki, K.; Consta, S. Disintegration mechanisms of charged aqueous nanodroplets studied by simulations and analytical models. J. Phys. Chem. B 2006, 110, 19168−19175. (23) Sharawy, M.; Consta, S. How do non-covalent complexes dissociate in droplets? A case study of the desolvation of dsDNA from a charged aqueous nanodrop. Phys. Chem. Chem. Phys. 2015, 17, 25550−25562. (24) Paliy, M.; Consta, S.; Yang, J. Interactions between Carbon Nanoparticles in a Droplet of Organic Solvent. J. Phys. Chem. C 2014, 118, 16074−16086. (25) Consta, S.; Chung, J. K. Charge-Induced Conformational Changes of PEG-(Na+)n in a Vacuum and Aqueous Nanodroplets. J. Phys. Chem. B 2011, 115, 10447−10455. (26) Glaser, M.; Schnauß, J.; Tschirner, T.; Schmidt, B. S.; MoebiusWinkler, M.; Käs, J. A.; Smith, D. M. Self-assembly of hierarchically ordered structures in DNA nanotube systems. New J. Phys. 2016, 18, 055001. (27) Ludwig, R. Water: from clusters to the bulk. Angew. Chem., Int. Ed. 2001, 40, 1808−1827. (28) Kuhne, T. D.; Pascal, T. A.; Kaxiras, E.; Jung, Y. New insights into the structure of the vapor/water interface from large-scale firstprinciples simulations. J. Phys. Chem. Lett. 2011, 2, 105−113. (29) Chen, X.; Yang, T.; Kataoka, S.; Cremer, P. S. Specific ion effects on interfacial water structure near macromolecules. J. Am. Chem. Soc. 2007, 129, 12272−12279. (30) Jorgensen, W. L.; Chandrasekhar, J.; Madura, J. D.; Impey, R. W.; Klein, M. L. Comparison of simple potential functions for simulating liquid water. J. Chem. Phys. 1983, 79, 926. (31) Grabuleda, X.; Jaime, C.; Kollman, P. A. Molecular dynamics simulation studies of liquid acetonitrile: New six-site model. J. Comput. Chem. 2000, 21, 901−908. (32) Lindahl, E.; Hess, B.; van der Spoel, D. GROMACS 3.0: a package for molecular simulation and trajectory analysis. J. Mol. Model. 2001, 7, 306−317. (33) Hess, B.; Kutzner, C.; van der Spoel, D.; Lindahl, E. GROMACS 4: Algorithms for Highly Efficient, Load-Balanced, and Scalable Molecular Simulation. J. Chem. Theory Comput. 2008, 4, 435−447. (34) Pronk, S.; Páll, S.; Schulz, R.; Larsson, P.; Bjelkmar, P.; Apostolov, R.; Shirts, M. R.; Smith, J. C.; Kasson, P. M.; van der Spoel, D.; et al. GROMACS 4.5: a high-throughput and highly parallel open source molecular simulation toolkit. Bioinformatics 2013, 29, 845. (35) Lindorff-Larsen, K.; Piana, S.; Palmo, K.; Maragakis, P.; Klepeis, J. L.; Dror, R. O.; Shaw, D. E. Improved side-chain torsion potentials for the Amber ff99SB protein force field. Proteins: Struct., Funct., Bioinform. 2010, 78, 1950−1958. (36) Martyna, G. J.; Klein, M. L.; Tuckerman, M. Nosé-Hoover chains: The canonical ensemble via continuous dynamics. J. Chem. Phys. 1992, 97, 2635−2643. (37) Phillips, J. C.; Braun, R.; Wang, W.; Gumbart, J.; Tajkhorshid, E.; Villa, E.; Chipot, C.; Skeel, R. D.; Kalé, L.; Schulten, K. Scalable molecular dynamics with NAMD. J. Comput. Chem. 2005, 26, 1781− 1802. (38) Vega, C.; de Miguel, E. Surface tension of the most popular models of water by using the test-area simulation method. J. Chem. Phys. 2007, 126, 154707. (39) Sharawy, M.; Consta, S. Effect of counterions on the charging mechanisms of a macromolecule in aqueous nanodrops. J. Chem. Phys. 2014, 141, 104321.

(40) Sheriff, F. A.; Consta, S. Charge-induced instabilities of droplets containing macromolecular complexes. Can. J. Chem. 2015, 93, 173− 180. (41) Liu, K.; Brown, M.; Carter, C.; Saykally, R.; Gregory, J.; Clary, D. Characterization of a cage form of the water hexamer. Nature 1996, 381, 501−503. (42) Keutsch, F. N.; Saykally, R. J. Water clusters: untangling the mysteries of the liquid, one molecule at a time. Proc. Natl. Acad. Sci. U. S. A. 2001, 98, 10533−10540. (43) Pérez, C.; Muckle, M. T.; Zaleski, D. P.; Seifert, N. A.; Temelso, B.; Shields, G. C.; Kisiel, Z.; Pate, B. H. Structures of cage, prism, and book isomers of water hexamer from broadband rotational spectroscopy. Science 2012, 336, 897−901. (44) Wang, Y.; Babin, V.; Bowman, J. M.; Paesani, F. The water hexamer: cage, prism, or both. Full dimensional quantum simulations say both. J. Am. Chem. Soc. 2012, 134, 11116−11119. (45) Stephen, P. S. Low viscosity magnetic fluid obtained by the colloidal suspension of magnetic particles. US Patent 3,215,572, 1965. (46) Cowley, M.; Rosensweig, R. E. The interfacial stability of a ferromagnetic fluid. J. Fluid Mech. 1967, 30, 671−688. (47) Richter, R. Mag (net) ic Liquid Mountains. Europhys. News 2011, 42, 17−19.

8880


Characterization of âStarâ Droplet Morphologies ... - ACS Publications

Recommend Documents

Characterization of âStarâ Droplet Morphologies ... - ACS Publications