Conformation-dependent dynamics of macromolecular ions in the gas phase under an electrostatic field: A molecular dynamics simulation

Abstract

Highly charged macromolecular ions exhibit various conformations in gas. Intramolecular charge-to-charge interaction induces a transformation from a globular structure into a stretched one. The change in the molecular conformation brings a complex dynamic behavior of ions under an electrostatic field. In the present study, we visualized the movement of a monovalent and multiply charged straight chain macromolecules, polyethylene glycol (PEG), by a molecular dynamics (MD) simulation. The simulation showed that a singly charged PEG ion (899 < M_W < 4,643) takes a globular conformation. The electrical mobilities of these ions determined from the migration distance under an electric field were compared with the experimental data and those determined by the classical Mason–Schamp theory under an assumption of spherical shape. As a result, we obtained a good agreement between the MD, theoretical, and experimental data for the monovalent ions. We also found that the MD simulation successfully predicts the electrical mobility of the multiply charged stretched PEG ions, but the classical theory fails. We were able to visualize the periodic bending and stretching motion by the MD simulation. This unique motion results from the localization of charges on the PEG molecule and may have a significant effect on the dynamic behavior of macromolecule ions in gas.

Copyright © 2019 American Association for Aerosol Research

1. Introduction

When airborne ions are placed in an electrostatic field, they migrate in the direction of electric field while experiencing drag by the collisions of gas molecules. The migration velocity under unit electrostatic field, that is, the electrical mobility, is a basic transport property of ions that reflects their size, charge state, and molecular structure. Ion mobility spectrometry (IMS) has been widely used for analyzing the molecular structures of bio-molecules by coupling it with electrospray ionization (ESI) (Ude, Fernandez de la Mora, and Thomson 2004; Ude and Fernandez de la Mora 2005; Trimpin et al. 2007; Brunelli, Flagan, and Giapis 2009; Larriba and Fernandez de la Mora 2012; Devine et al. 2017). The molecular structure analysis by the ESI/IMS is based on the hypothesis that the electrical mobility of airborne molecules is closely related to the molecular structure in a solution (Loo et al. 1997; Greig et al. 1995). However, previous studies (Ude, Fernandez de la Mora, and Thomson 2004; Trimpin et al. 2007; Larriba and Fernandez de la Mora 2012) suggested that the electrosprayed highly charged polymer chain takes various molecular conformations in gas as a result of charge-to-charge interactions in the molecule. According to the molecular dynamics (MD) simulation (in vacuum) of the molecular structure of polyethylene glycol (PEG) (Ude, Fernandez de la Mora, and Thomson 2004; Larriba and Fernandez de la Mora 2012), the molecular conformation of PEG (11 kDa) changed from a globular to a stretched conformation as the number of charge increased. Such transformation of molecular structure results in a complex dynamic behavior, causing a difficulty in the molecular structure analysis by the measurement of electrical mobility. The complexity in dynamic behavior may also have introduced a significant uncertainty in the filtration efficiency measurement using highly-charged macromolecular ions as test particles (Omori et al. 2017).

Molecular dynamics simulations are a powerful tool to estimate not only the molecular structure but also the transport properties of molecules (Liang et al. 2004; Ren, Shi, and Wang 2014). Our group has reported the applicability of an MD simulation for calculating the molecular diffusion coefficient in an ionic liquid system (Higashi et al. 2017) and the evaporation rate of electrosprayed nanodroplets (Higashi et al. 2015). In the present study, we applied the MD simulation to predict the electrical mobility of highly charged macromolecular ions accounting for the drag caused by the collisions of surrounding gas molecules (nitrogen).

2. Method

The simulations were carried out using a commercial MD software (Material Explorer 5.0 [FUJITSU Co., Ltd]). Several kinds of positively charged tetra-alkyl ammonium (TAA) ions (M_W = 74–579) and PEG ions (M_W = 899–4,643) were chosen as test molecules. The Dreiding force field (Mayo, Olafson, and Goddard 1990) was adopted to describe inter-and intra-molecular interaction (see supplementary information).

Figure 1 shows an example of a simulation cell. 200 N₂ molecules were randomly located in the cubic cell (L = 20 nm) which corresponds to the atmospheric pressure condition. Table 1 summarizes the molecular ions (TAA and PEG) used in the simulation. Molecular weights (M_w) of TAA ions were varied by the alkyl chain length from tetra methyl ammonium (M_w = 74) to tetra decyl ammonium (M_w = 579). PEG ions with 20, 30, 40, 60, 80, and 105 degrees of polymerization were also used (M_w = 899–4,643). The number of charges of the TAA ions were set to unity and the partial charge was calculated by the three-dimensional charge equilibration method (Pople and Gordon 1967). The PEG ions were multiply charged by attaching the ammonium ions (NH₄⁺) to the PEG molecule. The value of partial charges of the neutral PEG molecules and singly charged ammonium ion were taken from the reference (Consta and Chung 2011) and from NIH database values (https://pubchem.ncbi.nlm.nih.gov/compound/223#section=Top), respectively, as shown in Table 2. The molecular image in Table 1 shows an equilibrium structure of the singly charged PEG899 after a certain time of MD calculations (>0.5 ns). As shown in this image, NH₄⁺ ion was electrically fixed by surrounding oxygen atoms with partial negative charge. Similar electrical fixing of the NH₄⁺ ions was also observed for the multiply charged PEG ion systems depending on the PEG M_w. For example, four ions (z = 4) could be mounted on PEG899, while z = 7 on PEG4,643. Thus, PEG + zNH₄⁺ complex could be regarded as a single molecular ion. Consta and Chung also showed similar equilibrium conformation for sodiated (Na⁺-bonded) PEG molecules after certain time of pre-calculation (Consta and Chung 2011). Based on their report, we also conducted the calculation without electric field to obtain equilibrium molecular ion structures (typically for >0.5 ns). This pre-calculation was carried out at least for three times for each condition.

Figure 1. An example of the molecular dynamics simulation, 200 nitrogen gas molecules and an ion (TAA or PEG). The arrows are electric field direction of each component (E_x, E_y, and E_z) and total one (E).

Table 1. Calculation conditions for the TAA and PEG molecules.

	Tetra-alkyl ammonium (TAA)	Polyethylene glycol (PEG)
Molecular weight [−]	74–579	899–4,643
Charge number [−]	1	1–7
Chemical formula
Image of molecular model

Table 2. Partial charges on the PEG and NH4 molecules.

Molecule	Atom	Charge (e)
PEG (Consta and Chung 2011)	H at chain end	0.435
	O at chain end	−0.7
	CH2 at chain end	0.265
	CH2	0.25
	O	−0.5
NH4 (https://pubchem.ncbi.nlm.nih.gov/compound/223#section=Top)	H	−0.8
	N	0.45

To calculate the electrical mobility, an electrostatic field, E_x = E_y = E_z = 10 kV/cm (E in Figure 1), was applied as an external field by introducing the external electric field acting on all the partial charges across the entire system in the field’s direction. Calculation with the electrostatic field to x, y, and z directions has an advantage to compare the calculated velocities to three directions. The electric potential V is equal to product of electric field strength, E, and length of simulation cell, L (V = EL). The electrical mobility calculation was typically carried out for 5 ns. The number of calculation steps were set to 10⁷ (time step of calculation was 0.5 fs). In this simulation, a constant particle number, constant temperature (T = 298 K), and constant volume (Canonical ensemble) were applied in a simulation cell with periodic boundary conditions. In order to control the temperature of the system, velocity rescaling method (Woodcock 1971) was used. The rescaling was performed at each step. Two common integrators using the velocity Velret method and the Gear-predictor-corrector method are installed in the present MD simulator (Material Explorer 5.0). These methods have been widely used in molecular dynamics simulation because of the low calculation cost. In the present simulation, the gear-predictor-corrector method was chosen because of the higher accuracy. The accuracy of the calculation was insured by the gear 5th order predictor corrector method. The calculation was carried out five times for each condition.

3. Results and discussion

3.1. Singly charged ion

Figure 2(A), and Movie 1 in supplemental material show the time variation of the snapshot and a visualization of the calculation result, respectively. As shown in Figure 2(A), the singly charged PEG3,538 ion migrated in the same direction as the electric field, collided with nitrogen molecules, and retained its globular conformation. Such time variation of migration velocity was plotted in Figure 3. The open round symbol was the velocity averaged over five times calculation. As a result of acceleration of the ion in the small calculation cell, the nitrogen molecules also accelerated along the field direction, which causes an overall translation of the system. In order to cancel this effect, the relative velocity with respect to the center of gravity of nitrogen molecules was taken as the ion migration velocity. The solid line in Figure 3 represents theoretical velocity of the ions when the external forces acting on the molecular ion are electric force and drag force, which is expressed by equation S-6 in supplementary information. As the line in Figure 3 shows, the velocity rapidly increases for periods of up to 2 ns because the electrostatic force accelerating the ions is larger than the fluid resistance. Simultaneously, the fluid resistance increases as the velocity increases. After 2 ns, the mobility is presumably asymptotic to the constant value because the electric force and drag force are balanced. In general, the relaxation time (∼2 ns) is short enough to measure the experimental electrical mobility as a constant value.

Figure 2. Position of z charged PEG3,538 ion and nitrogen gas molecules against time in the simulation cell, (A) z = 1 (B) z = 5, (C) z = 7.

Figure 3. Time variation of electrical velocity of singly charged PEG3,538. Open circle and solid line represents the simulation result and the theoretical velocity.

As shown in Figure 3, the steady-state simulated migration velocity almost agreed with the theoretical calculation, though the random collisions with the N₂ molecules induced a slight jittering up and down. To compare this simulated velocity with experimental data, the electrical mobility was calculated by the following equation: (1) $$Z_{\mathrm{p}}=\frac{v_{\text{relaxation}}}{E},$$(1)

where, v_relaxation is the steady-state velocity averaged over the period from 3 to 5 ns.

Figure 4 shows the electrical mobility of singly charged TAA and PEG ions as a function of the M_w. The simulated results are compared with the experimental results taken from the reference (Ude and Fernandez de la Mora 2005; Trimpin et al. 2007). The simulated electrical mobility and experimental electrical mobility are represented by open and closed symbols, respectively. As shown in Figure 4, the simulated electrical mobility of TAA and PEG are in good agreement with the mobility determined experimentally. The line in Figure 4 represents the theoretical calculation result for a singly charged spherical PEG ion by assuming the bulk density of PEG, ρ = 1,260 kg/m³, using the following Mason–Scamp equation (Mason and McDaniel 1988): (2) $$Z_{\mathrm{p}}=\frac{3ze}{16N_0}\sqrt{\frac{2\pi }{\mu kT}}\frac{1}{\Omega },$$(2)

(3) $$\Omega =\xi \frac{\pi {\left(d+d_{\mathrm{g}}\right)}^2}{4},$$(3)

where N₀ is the number of gas molecules [m⁻³], μ is the reduced mass of the ion and gas molecules [kg], k is the Boltzmann constant [J/K], T is temperature [K], Ω is the collision cross section, d is the ion diameter [m], and d_g is the diameter of the gas molecules [m]. The values in Figure 4 were set as d_g = 0.30 nm and ξ = 1.36, obtained from the measurements of Millikan (Millikan 1923; Davies 1945). The simulated results also reasonably follow the theoretical line with the bulk density because the singly charged PEG in this simulation took an almost spherical conformation as predicted by MD simulation. For this reason, the MD simulation appears to serve capably for predicting the mobility of the molecular ion.

Figure 4. Molecular weight dependence of electrical mobility: experimental, simulated, and theoretical.

3.2. Multiply charged ion

Figure 2(B) and (C) show the simulation results for multiply charged PEG ions as a comparison with the singly charged ions (Figure 2A). The multiply charged ions (z = +5 and z = +7) also migrated in the same direction as the electric field, but the repulsive force between the charges on the ions caused them to take a straight-chain conformation rather than the spherical conformation of the singly charged ions. This transformation depending on the number of charges.

Figure 5 shows the time variation of the migration velocity of the singly (z = +1, closed circles) and multiply charged (z = +7, open circles) PEG ions (M_W = 3,538 Da). Each plot represents the averaged value of the calculation results (5 times) at each calculation time step. The solid line in Figure 5 shows the steady state value of electrical velocity (per charge) for the singly charged spherical ion expressed as v_sph = Z_p,sphE, where Z_p,sph is the experimentally measured mobility for singly charged PEG ions (Ude, Fernandez de la Mora, and Thomson 2004). The vertical axis of the 7-charged was set to be 1/7 of the simulated migration velocity (velocity per charge) for the comparison. As shown in Figure 5, the average value of the migration velocity of z = +1 was close to the theoretical value of the spherical ions, while that of multiple charged ion was much smaller than the spherical assumption. Such variation of the electrical velocity may be caused by the change in the molecular conformation from globular (z = +1) to stretched (z = +7). It is worth noting that the singly charged ion shows higher fluctuation of the velocity than that of the multiply charged ion. These fluctuations are triggered by the random collisions with nitrogen molecules. On the other hand, the normalized velocity of 7-charged ion was relatively stable because seven times higher electrical force was applied to the 7-charged ion. To compare to experimental results and investigate the effect of the structural transformation on the electrical velocity, we calculated the electrical mobility from the migration velocity. The open and closed symbols in Figure 6, respectively, represent the simulated and previous experimental results (Ude, Fernandez de la Mora, and Thomson 2004) for each number of charge. The solid line in Figure 6 represents the theoretical values for the spherical ion (Equation 2(2) $$Z_{\mathrm{p}}=\frac{3ze}{16N_0}\sqrt{\frac{2\pi }{\mu kT}}\frac{1}{\Omega },$$(2) ) with the bulk density, ρ = 1,260 kg/m³. The broken line in Figure 6 is the theoretical line including the influence of polarization (Kumar et al. 2014; Larriba and Hogan 2013; Higashi et al. 2015). To take into account for ion-induced polarization of gas molecules, the collision cross section was represented as follow: (4) $$\Omega =L\xi \frac{\pi {\left(d+d_{\mathrm{g}}\right)}^2}{4},$$(4)

where, L is the correction factor for polarization expressed as follow: (5) $$\begin{array}{c}L=1+{\Psi }_{\text{pol}}\left(\frac{1}{3.1}+\frac{1}{\xi }\left(\frac{1}{16}+\frac{4}{33}{\Psi }_{\text{pol}}\right)\right),\\ \text{for}{\Psi }_{\text{pol}}<1,\end{array}$$(5)

where, Ψ_pol is the polarization energy to thermal energy ratio expressed as follow: (6) $${\Psi }_{\text{pol}}=\frac{2{\alpha }_{\text{pol}}{z}^2{e}^2}{\pi \left(d+d_{\mathrm{g}}\right)2kT},$$(6)

where, α_pol is polarizability of gas molecule approximated as 1.7 × 10⁻³⁰ m³ for N₂ and O₂. Figure 6 also presents a snapshot of molecular image during the transition of the conformation (M_W = 3,538). As shown in the molecular model image of z = +7, the molecule took the form of a fibrous structure consisting of two components: a ball-shaped structure surrounding the ammonium ion and a fibrous structure stretching from ball to ball. The length of the fibrous structure gradually shortened as the charge number decreased, until the molecule finally assumed the globular shape shown in the images for z = +1, and +2. The tendency shown in the conformation was reflected in the difference between simulated and theoretical mobilities. The simulated mobility of the z = +1 PEG ion agreed with the theoretical line of the spherical ion. The mobility gradually diverged from the theoretical line as the number of charges increased. More importantly, the simulated mobility values of the multiply charged molecular ions agreed well with the experimental results. On this basis, we can assert that our simulation can also reproduce the conformation and dynamics of the multiply charged stretched ions under an electrostatic field.

Figure 5. Migration velocity per charge number of PEG3,538 ion against time. Open/Closed symbols and solid line were represented singly/septuple charged ion and the experimental migration velocity (Ude, Fernandez de la Mora, and Thomson 2004), respectively. The vertical axis of the 7-charged ion was set to be 1/7 of the simulated migration velocity (velocity per charge) for comparison.

Figure 6. Relationship between the electrical mobility and charge number (PEG3,538).

Next, we investigated the dynamics of the straight-chain (z = +7, M_W = 3,538 Da) PEG ions in an electrostatic field in detail. The orientation θ was defined as the angle made by the molecular stretching direction and the electric field. In the calculation of θ, the direction of the stretched PEG ion was determined by the line connecting the two ammonium ions at the both ends of PEG as shown by the lines in Figure 7(A–D). The time variation of this orientation is plotted in the upper part of Figure 7. As Figure 7 demonstrates, none of the PEG ions were orientated in a parallel direction, though theta fluctuated between 20° and 90°. More specifically, the two waves (i) and (ii) seen at 3.5–5 ns showed almost the same wavelengths and heights. Hence, the orientation was determined by another factor in addition to the random collision with gas molecules. As seen in the snapshots in Figure 7, the highly charged PEG ion was stretched perpendicularly to the electric field because the charges on the ion received the force uniformly. The charges were in fact unevenly distributed on the PEG ion, imposing a torque. The PEG ion rotated in an arbitrary direction (forward rotation) as a result, as shown in positions A to B. The macro molecular ion was flexible, however, resulting in a torque in the reverse direction that strained the molecule, as shown in B. The strain imposed during this rotation gradually decreased the torque in both the forward direction and reverse directions, although the decrease in the reverse-direction torque was less pronounced. As a result, the molecule rotated in the reverse direction and its orientation returned to the original position A. Such torque for rotation of forward and back direction may be changed by the non-uniformly distributed ions on the PEG molecule. The interaction between the charge and the electric field causes dipole moment to bend the PEG molecule. Therefore, the degree and frequency of the pendulum rotation may be dependent on the field strength and dipole moment. We regret that this hypothesis could not be reproduced in the simulation, however, as the dynamics was also affected by the collisions with the N₂ molecules. Our findings might indicate that the straight-chain molecular ion was transported with the periodic bending and stretching motion.

Figure 7. Time dependence of the orientation of the 7-charged PEG3,538 ion (upper subfigure) and the conformations and the orientation of PEG ion at the specific calculation time A, B, C and D (lower subfigure).

4. Conclusion

An MD simulation of singly and multiply charged molecular ions (TAA and PEG) in N₂ gas was performed. The conformation of the molecular ions was visualized, and the electrical mobility was calculated. The migrating singly charged TAA and PEG ions took almost spherical shapes, and the simulated mobility was in good agreement with previous experimental data as well as theory. The multiply charged PEG ions took various conformations depending on number of charges and M_w. A reasonable agreement between calculated electrical mobility and previous experimental data was obtained. Consequently, our MD simulation was confirmed to capably reproduce the dynamics and conformation of the multiply charged ions under transportation in an electric field. We also investigated the periodicity of the orientation. We thus, conclude that the highly charged PEG ion vibrated as it migrated under an electrostatic field in the gas phase, and that the vibration might influence the transport properties.

Additional information

Funding

This work was supported by JSPS KAKENHI Grant Number 15H04171, and JST CREST Grant Number JPMJCR18H4, Japan