Dependence of interaction free energy between solutes on an external electrostatic field.

To explore the athermal effect of an external electrostatic field on the stabilities of protein conformations and the binding affinities of protein-protein/ligand interactions, the dependences of the polar and hydrophobic interactions on the external electrostatic field, -Eext, were studied using molecular dynamics (MD) simulations. By decomposing Eext into, along, and perpendicular to the direction formed by the two solutes, the effect of Eext on the interactions between these two solutes can be estimated based on the effects from these two components. Eext was applied along the direction of the electric dipole formed by two solutes with opposite charges. The attractive interaction free energy between these two solutes decreased for solutes treated as point charges. In contrast, the attractive interaction free energy between these two solutes increased, as observed by MD simulations, for Eext = 40 or 60 MV/cm. Eext was applied perpendicular to the direction of the electric dipole formed by these two solutes. The attractive interaction free energy was increased for Eext = 100 MV/cm as a result of dielectric saturation. The force on the solutes along the direction of Eext computed from MD simulations was greater than that estimated from a continuum solvent in which the solutes were treated as point charges. To explore the hydrophobic interactions, Eext was applied to a water cluster containing two neutral solutes. The repulsive force between these solutes was decreased/increased for Eext along/perpendicular to the direction of the electric dipole formed by these two solutes.

1. Introduction

Among the four fundamental interactions, the electrostatic force dominates non-covalent bond interactions between atoms in biomolecules such as proteins, DNA, and RNA. An external electrostatic field may alter the electrostatic interactions among atoms in proteins, and consequently change the stabilities of protein conformations or the binding affinities of protein-protein/ligand interactions. Changing the stability of the protein conformation may change the protein activity, and changing the binding affinities of the protein-protein/ligand interactions may change the regulation of the signal transduction network in cells or the expression of proteins. Although these effects may not cause diseases immediately, they may increase the possibility of diseases developing. Further, electromagnetic radiation is used to kill bacteria for food preservation. Potential problems have been studied [1,2] using protein experiments [3-7], cellular experiments [8-11], animal experiments [12], public health data [13], and computer simulations [14-16]. The results did not verify whether electromagnetic radiation is harmful to humans.

In most cases, the interactions in biomolecules or among biomolecules are regarded as the sum of the interactions among the atoms in the biomolecules. Understanding the interactions among atoms in biomolecules is therefore helpful in exploring problems such as the stability of protein conformations and the binding affinities of protein-protein/ligand interactions. In addition, understanding the effects of an external electrostatic field, −ext, on the interactions among atoms in biomolecules is helpful in exploring the effects of the external electric field on the stabilities of protein conformations or the binding affinities of protein-protein/ligand interactions. To explore the effects of ext on the mean force between two charged atoms in water, ext was decomposed into two components: one along the electric dipole formed by these two charged atoms, and the other perpendicular to this electric dipole. The effect of ext on the mean force between the two charged atoms in water was estimated by summing the effects on the mean force from these two components.

Most biomolecules exist in aqueous environments, and the interactions among atoms in water differ from those in vacuum. The solvent effect is significant for the stabilities of protein conformations [17-21] and the binding affinities of protein-protein/ligand interactions [22-25]. Numerous strategies have been developed for computing the solvation free energy [26-36]. By treating the charged atom as a point charge, the effect of ext on the interactions among charged atoms in water can be quickly estimated using Coulomb’s law. However, strategies using molecular dynamic (MD) simulations with explicit solvent models afford more accurate results.

For one or two charged atoms in a spherical water cluster, the dependence of the electric dipole of TIP3P water on the net electrostatic field at the oxygen atom of TIP3P water, and the dependence of the radial distribution function of TIP3P water on the mean force at the oxygen atom of TIP3P water, have been explored [33,37,38]. In this project, the source code of the Charmm package [39] was modified, and it was verified that it can be applied to an external electrostatic field. The external electrostatic field was applied to a simulation system containing one charged atom or two charged atoms in vacuum; the accelerations, velocities, and positions of atoms from MD simulations using the modified Charmm package were consistent with those obtained from analytical solutions (data not shown). The external electrostatic field along the x or y direction (EXext, EYext) was applied to a pure water cluster, and the relation between the dipole moment of TIP3P water and the net electrostatic field at the oxygen atom of TIP3P water was consistent with the results obtained in previous works [37,38]. The external electrostatic field along the x or y direction (EXext, EYext) was applied to a water cluster containing one or two charged atoms, and the dependences of the mean force and the potential of mean force (PMF) between the charged solutes on ext were studied using MD simulations [40]. The differences between the mean force estimated from a continuum solvent and that computed using MD simulations were discussed.

2. Results/Discussion

2.1. Dependence of p on Eext

On applying ext to a water cluster containing solutes, the water molecules were polarized by the ext. The electric force on the atoms of the solute from ext could be shielded by the polarized water molecules. To explore the effect of the electric dipole per water molecule, −, polarized by ext, ext was applied to the water cluster (Figure 1a), and and the net electric field on the water molecule, −net were computed from the trajectories of MD simulations using (1) and (3), respectively. The results showed that was proportional to ext in the region |Eext| < 50 MV/cm. The ratio of p to Eext was approximately 0.007 eÅ/(MV/cm), as pE = 0.007 eÅ/(MV/cm) * Eext (MV/cm). The p (black line) from MD simulations was compared with pE (gray line) (Figure 1b). With regard to the relationship between Enet and Eext, Enet was proportional to Eext. The proportionality constant was larger in the region |Eext| < 50 MV/cm than in the region |Eext| > 50 MV/cm (Figure 1c). Based on these results, the external electrostatic field along the x direction, −EXext = 40, 60, or 100 MV/cm, or the external electrostatic field along the y direction, −EYext = 50 or 100 MV/cm, was applied to a water cluster containing two solutes to compute the mean force and the PMF on the solute S2.

Figure 1

Dependences of p and Enet on Eext. (a) The external electrostatic field, −ext, was applied to a water cluster with a radius of 20 Å containing 1119 TIP3P [41] water molecules; (b) The time-averaged dipole moment per water molecule, − (black line), was computed from the trajectories of molecular dynamics (MD) simulations. The p (gray line) was plotted as 0.007 eÅ/(MV/cm) *Eext (MV/cm); (c) The time-averaged electrostatic field at the oxygen atom of TIP3P water, −net (black line), was computed from the trajectories of MD simulations.

The dependence of on the charged solute and net has been extensively studied [37,38]. A solute with a charge of −4.0 e, −3.0 e, −2.0 e, −1.0 e, −0.8 e, −0.6 e, −0.4 e, −0.2 e, +0.2 e, +0.4 e, +0.6 e, +0.8 e, +1.0 e, +2.0 e, +3.0 e, or +4.0 e was at the center of a spherical water cluster (Figure 2a), and and net were computed from the trajectories of MD simulations (Figure 2b). The solvent molecular polarizability ɛ0γmol was defined as dp/dEnet. Δp was proportional to ΔEnet in the region |Enet| < 25 kcal/(mol·eÅ) (1 MV/cm = 0.231 kcal/(mol·eÅ)), and ɛ0γmol = 0.0124 (mol·e2·Å2)/kcal (or γmol = 51.7 Å3[37]). For Enet > 50 kcal/(mol·eÅ), the net dipole of TIP3P is towards the direction of net, and the value of p approaches +0.49 eÅ. For Enet < −100 kcal/(mol eÅ), one of the hydrogen atoms was toward the anion, and the value of p approaches −0.35 eÅ. In the region −100 kcal/(mol·eÅ) < Enet < −50 kcal/(mol·eÅ), the net dipole of TIP3P or one of the hydrogen atoms could be toward the anion, so the value of p depends not only on Enet, but also depends on the solute charge and position in calculating net. The values of p therefore varied between −0.49 eÅ and −0.35 eÅ.

Figure 2

Dependence of p on Enet. (a) A solute with a charge of −4.0 e, −3.0 e, −2.0 e, −1.0 e, −0.8 e, −0.6 e, −0.4 e, −0.2 e, +0.2 e, +0.4 e, +0.6 e, +0.8 e, +1.0 e, +2.0 e, +3.0 e, or +4.0 e was at the center of a spherical water cluster of radius 20 Å containing 1118 TIP3P [41] water molecules. The van der Waals parameters of the solute assigned were the same as those for the oxygen atom of TIP3P water with ɛ = −0.1521 kcal/mol and Rmin/2 = 1.7682 Å; (b) The (Q2, r) and net(Q2, r) were computed from the trajectories of MD simulations. For |Q2| ≤ 4 e and r ≤ 10 Å, dependence of on net was shown (black line). Because the TIP3P water model is not a point dipole moment, does not only depend on net. The pE (gray line) was plotted as 0.0124 (mol·e2·Å2)/kcal * Enet [kcal/(mol·eÅ)].

The p from Eext was compared with that from the charged solute. The of water molecules polarized by ext (MV/cm) (Figure 1a) was compared with that polarized by a solute with charge Q(e) (Figure 2a). For example, when Eext = 50 MV/cm was applied to a pure water cluster, p = 0.35 eÅ (Figure 1b), which is 70% of the permanent electric dipole moment of TIP3P water. For the solute with a charge of +1.0 e in water, and the same van der Waals (vdW) parameters as the oxygen atom of TIP3P water, p at the first peak of the radial distribution function, solute distance 2.8 Å, was 0.35 eÅ [37,38]. The subconclusion is that p polarized by Eext = 50 MV/cm was similar to p at the first peak of the radial distribution function surrounding the solute with a charge of +1.0 e.

Eext shifts the equilibrium position of the p − Enet curve of the water molecule, and could reduce the electric shielding effect between charged atoms in macromolecules. Consider a macromolecule such as a protein in water solvent; the solvent molecules are polarized by the charged atoms in the macromolecules, and the polarized molecules shield the electrostatic interactions between the charged atoms in the macromolecules. The solvent molecular polarizability, γmol, describes the electric shielding effect between charged particles in dielectrics. In the case of no applied ext, the equilibrium position is at position A (Figure 2b). The Δ is proportional to Δnet in the region |Enet| < 25 kcal/(mol eÅ), and ɛ0γmol = 0.0124 (mol·e2·Å2)/kcal. Applying ext = 200 MV/cm to the water cluster leads to Enet = 66 kcal/(mol·eÅ). The equilibrium position is shifted to position B (Figure 2b). The ɛ0γmol for computing the electric shielding effect between charged atoms in macromolecules is 0.0006 (mol·e2·Å2)/kcal. Applying ext = 50 MV/cm to the water cluster leads to Enet = 29 kcal/(mol·eÅ). The equilibrium position is shifted to position C (Figure 2b). The γmol was 0.0007/0.0022 (mol·e2·Å2)/kcal in the direction of increasing/decreasing p. The subconclusion is that the electric shielding effect of water for computing the electric interactions between charged atoms in macromolecules could be reduced by Eext.

2.2. Dependence of Fnet(One_Atom) on Eext

To understand the effect of ext on the force on the charged solute, ext was applied to a water cluster of radius 20 Å, and the mean force on the solute S2 was computed using the trajectories of the MD simulations (Figure 3a). The net mean force, −net, was contributed by the external electrostatic field −E and the polarized water molecules −solv, using (9) and (10) [38,42]. The results showed that |Fnet| was small in the |Eext| < 50 MV/cm region because the E force was almost balanced by the solv force (Figure 3b). For |Eext| > 50 MV/cm, the dielectric water approached saturation, |Fsolv| increased slowly as Eext increaed, and d|Fnet|/dEext approached a constant (Figure 3b). To understand the dependences of the net and solv forces on the radius of the water cluster, ext was applied to a water cluster of radius 25 Å containing one solute with a charge of +1.0 e. The results showed that the net and solv computed from the water cluster of radius of 25 Å were similar to those obtained using a radius of 20 Å (Figure 3b).

Figure 3

Dependence of F(one_atom) on Eext. (a) The external electrostatic field ext was applied to a water cluster containing one charged solute S2. The solute S2 was at the center of a spherical water cluster of radius 20 or 25 Å containing 1118 or 2185 TIP3P [41] water molecules. To explore the general effect of ext on the atoms in biomolecules, the van der Waals parameters of the solute were assigned to be the same as the oxygen atom of TIP3P water with ɛ = −0.1521 kcal/mol and Rmin/2 = 1.7682 Å. For Q2 = +1 e (b) or −1 e (c), the ensemble average force on the charged solute contributed by Eext − FE (gray line), the polarized water molecules, −solv (dashed line), and the net force −Fnet (solid line) were computed from the trajectories of MD simulations with an amplitude of Eext ranging from 0.1 to 150 MV/cm. The radii of the water clusters were 20 Å (black line) and 25 Å (red line), respectively.

The TIP3P water molecule contained one oxygen and two hydrogen atoms. The vdW radius of the oxygen atom, Rmin/2 = 1.7682 Å, was larger than that of the hydrogen atom, Rmin/2 = 0.2245 Å. For the cation in water, the oxygen atom of water was closer to the cation. In contrast, for the anion in water, the hydrogen atom of water was closer to the anion. The radial distribution function of the oxygen or hydrogen atoms surrounding the cation therefore differed from that of those surrounding the anion [37]. However, the amplitude of net on the solute with a charge of +1 e was similar to that on the solute with a charge of −1 e (Figure 3c). This means that the net(one_atom; Eext) was independent of the sign of the charged solute.

2.3. Attractive Force between S1 and S2 Could Not Be Decreased by Applying an External Electrostatic Field along the Direction of the Electric Dipole Formed by S1 and S2

Polar interactions, such as those between hydrogen bond donors and acceptors, play a significant role in stabilizing protein conformations and protein-protein/ligand complex structures. On applying ext to macromolecules containing polar interactions, the electric dipole formed by the two atoms with opposite charges, S1 and the S2, prefers to align in the direction of ext, to reduce the potential energy, and ext pulls S1 and pushes S2 along the direction of ext (Figure 4a). The attractive force and the attractive interaction free energy between S1 and S2 decreased if S1 and S2 were treated as point charges in continuum dielectrics.

Figure 4

FX(two_atoms; EXext) − d curves. (a) An external electrostatic field along the x direction, EXext, was applied to a water cluster containing S1 and S2 solutes. The solute S1, with charges Q1 = −1 e positioned at (−d/2, 0, 0), and the solute S2, with charges Q2 = +1 e positioned at (+d/2, 0, 0), were in a spherical water cluster of radius 20 Å containing 1117 TIP3P [41] water molecules. The distance between the S1 and S2 solutes was from 2 to 8 Å. The vdW parameters of the S1 and S2 solutes were assigned to be the same as those of the oxygen atom of TIP3P water with ɛ = −0.1521 kcal/mol and Rmin/2 = 1.7682 Å; (b) The mean forces along the x direction, −FX, on the solute S2 were computed from the trajectories of MD simulations with external electrostatic fields EXext = 0 (solid black), 40 MV/cm (dashed gray), 60 MV/cm (dashed black), and 100 MV/cm (solid gray), respectively; (c) The potentials of mean force were computed from the mean forces in (b). FX(two_atoms; EXext) (black line) was compared with FXest(two_atoms; EXext) (gray line) for EXext = 40 MV/cm (d), 60 MV/cm (e), and 100 MV/cm (f); (g) EXext was applied to a water cluster containing the S2 solute. (h) The force on S2 was contributed by the polarized water molecules in the space occupied by S1. The dielectric polarization in the space occupied by S1 in (h) was the reverse of the dielectric polarization in the space occupied by S1 in (g).

For S1 and S2 exposed to water, the external electrostatic field along the x direction, −EXext, was applied to the water cluster containing S1 and S2 solutes (Figure 4a). The mean force on S2 along the x direction, −FX(two_atoms; EXext), was computed using the trajectories of the MD simulations. The results showed that FX(two_atoms; EXext = 0) was attractive (negative) in the 2.7 Å < d < 3.4 Å region, and the minimum value of FX(two_atoms; EXext = 0) was −3.7 kcal/(mol·Å) at the position d = 3.0 Å (Figure 4b). We also applied EXext = 40 or 60 MV/cm to a water cluster containing one solute with a charge of +1 e, and the mean force FX(one_atom; EXext) on the solute with a charge of +1 e was positive (Figure 3a). Was the attractive force on the solute S2 in Figure 4a decreased by application of the external electrostatic field? The results from the MD simulations showed that FX(two_atoms; EXext = 40 or 60 MV/cm) was more attractive than FX(two_atoms; EXext = 0) in the region 2.7 Å < d < 3.4 Å. The minimum values of FX(two_atoms; EXext) at position d = 3.0 Å were −4.4 and −5.0 kcal/(mol Å) for EXext = 40 and 60 MV/cm, respectively (Figure 4b).

The PMF, −PMF(two_atoms; EXext = 40 or 60 MV/cm), was calculated by integration of the mean force X from infinity. For comparison of the energies needed to escape the first well of PMF(two_atoms; EXext = 40 or 60 MV/cm), the second peak of PMF(two_atoms; EXext), was set at zero. The results showed that the depth of the first well of PMF(two_atoms; EXext = 40 or 60 MV/cm) was deeper than that of PMF(two_atoms; EXext = 0) (Figure 4c).

Treating S1 and S2 as point charges, FX(two_atoms; EXext) was estimated by summation of FX(two_atoms; EXext = 0) (Figure 4b) and FXnet(one_atom; EXext) (Figure 3b). FXest(two_atoms; EXext) was the sum of FX(two_atoms; EXext = 0) in Figure 4b and FXnet(one_atom; EXext) in Figure 3b. The results showed that FXest(two_atoms; EXext) was larger than FX(two_atoms; EXext), especially in the d < 3.4, 3.8, and 4.4 Å regions for EX = 40, 60, and 100 MV/cm, respectively (Figure 4d−f). This is because no water molecules can be polarized in the space occupied by S1. If S1 occupies the space, FX(two_atoms; EXext) should be estimated by summation of FX(two_atoms; EXext = 0) (Figure 4b), FXnet(one_atom; EXext) (Figure 4g), and FX(excluded_solvent; EXext) (Figure 4h), based on the superposition principle. FX(excluded_solvent; EXext) is the force on solute S2 contributed by the water in the space occupied by S1. The dielectric polarization in the space occupied by S1 in Figure 4h was the reverse of the dielectric polarization in the space occupied by S1 in Figure 4g. FX(excluded_solvent; EXext) was along the −x direction, therefore FXest(two_atoms; EXext) was larger than FX(two_atoms; EXext).

2.4. Attractive Force between S1 and S2 Was Unchanged and Increased by Applying EYext = 50 MV/cm and 100 MV/cm, Respectively

The external electrostatic field along the y direction, −EYext, was applied to a water cluster containing S1 and S2 solutes (Figure 5a). The mean force on S2 along the x direction, −FX(two_atoms; EYext), was computed using the trajectories of the MD simulations. The results showed that FX(two_atoms; EYext = 50 MV/cm) was similar to FX(two_atoms; Eext = 0), but FX(two_atoms; EYext = 100 MV/cm) was less than FX(two_atoms; Eext = 0) (Figure 5b). The difference between FX(two_atoms; EYext = 50 MV/cm) and FX(two_atoms; Eext = 0) at the position of the first minimum of FX(two_atoms; Eext = 0) was 0.3 kcal/(mol Å), and the difference between FX(two_atoms; EYext = 100 MV/cm) and FX(two_atoms; Eext = 0) at the position of the first minimum of FX(two_atoms; Eext = 0) was 2.2 kcal/(mol·Å) (Figure 5b).

Figure 5

FX(two_atoms; EYext) − d curves. (a) An external electrostatic field along the y direction, EYext, was applied to a water cluster containing S1 and S2 solutes. The others were the same as those in Figure 4a; (b) The mean forces along the x direction, −FX, on solute S2 were computed from the trajectories of MD simulations with external electrostatic fields EY = 0 (solid black), 50 MV/cm (dashed black), and 100 MV/cm (solid gray).

Eext was applied to the charged atom in the water cluster (Figure 2a), the force on the charged atom was along the direction of ext, and the force perpendicular to direction of ext was zero. S1 and S2 were treated as point charges. ext was applied perpendicular to the direction of the electric dipole formed by S1 and S2; the force on S1 and S2 from ext was along the direction of ext, and the attractive force and the interaction potential energy between S1 and S2 were unchanged. When Eext = 100 MV/cm was applied perpendicular to the direction of the electric dipole formed by these two atoms (Figure 5a), the attractive force between S1 and S2 increased (Figure 5b), and the interaction free energy also increased, as observed from MD simulations. This is because polarization of the water molecules was saturated on application of EYext = 100 MV/cm. The saturated water molecule is hard to polarize further by S1 and S2. The dielectric shielding effect between S1 and S2 was therefore reduced.

2.5. FY(Two_Atoms; EYext) Was Greater than FY(one_atom; EYext), Especially When d Was Small

EYext was applied to a water cluster containing two charged solutes (Figure 6a); the mean force on S2 along the y direction, −FY(two_atoms; EYext), was computed using the trajectories of the MD simulations. The results showed that FY(two_atoms; EYext = 50 or 100 MV/cm) was similar to FY(one_atom; EYext = 50 or 100 MV/cm), except in the region d < 3.0 Å (Figure 6b). The differences between FY(two_atoms; EYext) and FY(one_atom; EYext) at d = 6 Å were 0.2 and 0.3 kcal/(mol Å) for EYext = 50 and 100 MV/cm, respectively.

Figure 6

FY(two_atoms; EYext) − d curves. (a) The external electrostatic field along the y direction, −EYext, was applied to a water cluster containing two charged solutes. The other parameters were the same as those in Figure 4a; (b) The mean forces along the y direction, −FY, on solute S2 were computed from the trajectories of MD simulations with external electrostatic fields EY = 50 MV/cm (solid black) and 100 MV/cm (solid gray). FYnet(two_atoms; EYext) was compared to the force on solute S2 in Figure 3a with external electrostatic fields EY = 50 MV/cm (dashed black) and 100 MV/cm (dashed gray); (c) EYext was applied to a water cluster containing S2; (d) The force on S2 was contributed by the polarized water molecules in the space occupied by S1. The dielectric polarization in the space occupied by S1 in (d) was the reverse of the dielectric polarization in the space occupied by S1 in (c).

The force on S2 from S1 was along the −x direction. Treating S1 and S2 as point charges, FY(two_atoms; EYext) should be the same as FY(one_atom; EYext). However, the force on S2 in the water cluster containing two solutes was greater than the force on the charged solute S2 in the water cluster containing one solute (Figure 6b). This is because no water molecules in the space occupied by S1 can be polarized. FY(two_atoms; EYext) should be estimated by summation of FY(one_atom; EYext) (Figure 6c) and FY(excluded_solvent; EYext) (Figure 6d), based on the superposition principle. FY(excluded_solvent; EYext) is the force on solute S2 contributed by the water molecules in the region occupied by solute S1. The dielectric polarization in the space occupied by S1 in Figure 6d was the reverse of the dielectric polarization in the space occupied by S1 in Figure 6c. FY(excluded_solvent; EYext) was along the +y direction (Figure 6d), therefore FY(one_atom; EYext) was less than FY(two_atoms; EYext).

2.6. Dependence of FX(Two_Neutral_Atoms) on EXext and EYext

Hydrophobic interactions play a significant role in the stabilities of protein conformations and the binding affinities of protein-protein/ligand interactions. The effect of an external electrostatic field on the mean force between two neutral solutes was explored. Consider two neutral solutes in a water cluster (Figure 7a). The mean force on S2, −FX(two_neutral_atoms; Eext), was computed using the trajectories of the MD simulations. The results showed that FX(two_neutral_atoms; Eext = 0) was attractive in the region 3.2 Å < d < 5.0 Å, and repulsive in the region 5.0 Å < d < 6.0 Å (Figure 7a). The maximum attractive force was −0.8 kcal/(mol Å) at the d = 3.8 Å position, and the maximum repulsive force was 0.4 kcal/(mol Å) at the d = 5.6 Å position (Figure 7b).

Figure 7

FX(two_neutral_atoms; EXext/EYext)−d curves. (a) Two neutral solutes separated by a distance d along the x direction were in a TIP3P water cluster of radius 20 Å. An external electrostatic field along the x or y direction was applied to the water cluster. The other parameters were the same as those in Figure 4a; (b) The mean forces along the x direction, FX, on solute S2 were computed from the trajectories of MD simulations with Eext = 0 (solid black), EXext = 100 MV/cm (dashed black), and EYext = 100 MV/cm (solid gray); (c) The potentials of mean force were computed from the mean forces in (b); (d) For solutes S1 and S2 separated by a distance greater than 5.4 Å, the space between S1 and S2 can accommodate a water molecule, W1; (e) On applying EXext = 100 MV/cm to this water cluster, the water molecules will be polarized along the x direction. W1 will be pushed by the neighboring water molecules, W2 and W3; (f) On applying EYext = 100 MV/cm to this water cluster, the water molecules will be polarized along the y direction. W1 will be attracted by the neighboring water molecules, W2 and W3.

EXext = 100 MV/cm or EYext = 100 MV/cm was applied to a water cluster containing two neutral solutes (Figure 7a); the mean force on S2, −FX(two_neutral_atoms; Eext), was computed using the trajectories of the MD simulations. The results showed that FX(two_neutral_atoms; EXext = 100 MV/cm) and FX(two_neutral_atoms; EYext = 100 MV/cm) differed from FX(two_neutral_atoms; Eext = 0). FX(two_neutral_atoms; EXext = 100 MV/cm) was larger than FX(two_neutral_atoms; Eext = 0) in the region d < 4.4 Å, and less in the region 4.4 Å < d < 6.0 Å (Figure 7b). In contrast, FX(two_neutral_atoms; EYext = 100 MV/cm) was similar to FX(two_neutral_atoms; Eext = 0) in the region d < 4.4 Å, and FX(two_neutral_atoms; EYext = 100 MV/cm) was larger than FX(two_neutral_atoms; Eext = 0) in the region 4.4 Å < d < 6.0 Å (Figure 7b).

The PMF, −PMF(two_neutral_atoms; Eext), was calculated by integration of the mean force FX(two_neutral_atoms; Eext) from infinity. For comparison of the depths of the first wells of PMF(two_neutral_atoms; Eext) with EXext = 100 MV/cm or EYext = 100 MV/cm, the second peak of PMF was set to zero. The results showed that the depth of the first well of PMF(two_neutral_atoms; Eext = 0) was smaller than that of PMF(two_neutral_atoms; EXext = 100 MV/cm), and similar to that of PMF(two_neutral_atoms; EYext = 100 MV/cm) (Figure 7c).

For solutes S1 and S2 separated by a distance greater than 5.4 Å, the space between S1 and S2 can accommodate a water molecule (Figure 7d). When one of the water molecules, e.g., W1, stays at the position between solutes S1 and S2, W1 will push S2 along the +x direction. When EXext = 100 MV/cm is applied to this water cluster, the water molecules will be polarized along the x direction. W1 will be pushed by the neighboring water molecules, W2 and W3 (Figure 7e); the probability of one of the water molecules staying at the position of W1 when EXext = 100 MV/cm was applied was lower than that without application of an external electrostatic field. FX(two_neutral_atoms; EXext = 100 MV/cm) was therefore less than FX(two_neutral_atoms; Eext = 0) in the region 4.4 Å < d < 6.0 Å (Figure 7b). If EYext = 100 MV/cm is applied to this water cluster, the water molecules will be polarized along the y direction. W1 will be attracted by the neighboring water molecules, W2 and W3 (Figure 7f); the probability of one of the water molecules staying at the position of W1 under application of EYext = 100 MV/cm was larger than that without application of an external electrostatic field. FX(two_neutral_atoms; EYext = 100 MV/cm) was therefore larger than FX(two_neutral_atoms; Eext = 0) in the region 4.4 Å < d < 6.0 Å (Figure 7b).

3. Method

3.1. MD Simulations

The simulations were performed in an NVE ensemble using the CHARMM package [39] and spherical boundary conditions. The ion-water and water-water interaction energies were calculated by summation of the electrostatic and vdW pairwise energies with a non-bond cutoff of 99 Å. For TIP3P water, the charge states of the oxygen and hydrogen atom were −0.834 e and +0.417 e, respectively, and the vdW parameters of the hydrogen atom were ɛ = −0.046 kcal/mol and Rmin/2 = 0.2245 Å. The O–H bond length of TIP3P, 0.9572 Å, and the bond angle of H–O–H, 104.52°, were constrained during the simulations using the SHAKE algorithm [43]. The intrinsic electronic polarizability of the water molecule changed as a strong electric field [44] was not considered in this project. All atoms were propagated according to Newton’s equations using the leapfrog Verlet algorithm and a time-step of 2 fs at a mean temperature of 300 K. Each system was first minimized for 1000 steps, equilibrated for 200 ps, and subsequently subjected to 1 ns of production. The configurations were stored every 20 fs.

3.2. Application of Eext to Pure Water Cluster and Calculation of p and Enet from Trajectories of MD Simulations

Eext was applied to a pure water cluster (Figure 1a), and the electric dipole moment per water molecule, p, was calculated using the sum of the electric dipole moments of water molecules with an oxygen atom distance origin ≤ rcut over NC configurations/frames, divided by the number of water molecules N with an oxygen atom distance origin ≤ rcut over NC configurations/frames:

where qi is the charge on water atom i, ilm denotes the coordinates of atom i of water molecules m in configuration l, Olm denotes the coordinates of the oxygen atom of water molecule m in configuration l, n is the number of solvent molecules in the simulation system, Nc is the number of configurations/frames collected in equilibrium state in the MD simulations, and u(rcut − rOlm) is the Heaviside unit step function.

N was computed as the sum of the water molecules with oxygen atoms positioned at distance origin ≤ rcut:

where the first summation is over Nc configurations/frames, and the second summation is over the n solvent molecules in the simulation system.

The electrostatic field at the oxygen atom of TIP3P water, −net, was contributed by ext and water molecules with oxygen atom distance origins ≤ rcut over NC configurations/frames, divided by the number of water molecules N with oxygen atom distance origins ≤ rcut over NC configurations/frames as

where ɛ0 is the permittivity of free space, qi is the charge on water atom i, lmm′iO is the vector from atom i of water molecule m to the oxygen atom of water molecule m’ in configuration l.

3.3. For Charged Atom in Water Cluster, Calculation of p(r) and Enet(r) from Trajectories of MD Simulations

For one charged atom in a water cluster (Figure 2a), p(r) was calculated by summing the electric dipole moments of water molecules with oxygen atoms located between (r − Δr/2) and (r + Δr/2) over NC configurations/frames, divided by the number of water molecules N(r), as

N(r) in (4) was computed as the sum of the number of water molecules whose oxygen atoms were at a distance from the solute of between (r − Δr/2) and (r + Δr/2):

where the first summation is over Nc configurations/frames, the second summation is over the n solvent molecules in the simulation system, Olm denotes the coordinates of the oxygen atom of water molecule m in configuration l, and Δr is set at 0.1 Å.

Enet(r) was calculated by summing the electrostatic fields of water with its oxygen atom located between (r − Δr/2) and (r + Δr/2) over NC configurations/frames, divided by the number of water molecules N(r), as

where ɛ0 is the permittivity of free space, qi is the charge on water atom i, qj is the charge on the solute atom, Rlm′jO is the distance between the oxygen atom of water molecule m’ in configuration l and solute j, lmm′iO is the vector from atom i of water molecule m to the oxygen atom of water molecule m’ in configuration l.

3.4. Application of Eext to Water Cluster Containing One or Two Solutes and Calculation of Fsolv(r) and Fnet(r) from Trajectories of MD Simulations

The net mean force on solute S2 with charge Q2, −Fnet, was decomposed and contributed by the external electrostatic field −FE, solute S1 with charge Q1, −FS1, and the polarized solvent molecules, −Fsolv. FE was computed as FE = Q2E.

FS1 was the force acting on solute S2 at r2 because of solute S1 at r1, and can be computed as the sum of the electrostatic and vdW forces as

where R12 = r2 − r1, and the vdW parameters, ɛ12 and Rmin,12, were obtained using the standard combining rules.

Fsolv was the force acting on solute S2 at r2 because of the solvent molecules, and can be computed as the sum of the electrostatic and vdW forces as

where the first summation is over Nc configurations/frames, the second summation is over the n solvent molecules in the simulation system, qj is the charge on water atom j, Rlm2j is the vector from atom j of water molecule m to solute S2 at r2 in configuration l, and the vdW parameters, ɛ2j and Rmin,2j, were obtained using the standard combining rules.

4. Conclusions

To explore the athermal effect of ext on the stabilities of protein conformations and the binding affinities of protein-protein/ligand interactions, the dependence of the mean force between charged solutes or neutral solutes, S1 and S2, on ext was studied using MD simulations. The results showed that (1) ext shifts the equilibrium position of the p − Enet curve of the water molecule, and may reduce the dielectric shielding effect between charged atoms in macromolecules; (2) For ext along the direction of the electric dipole formed by S1 and S2, Eext = 40 or 60 MV/cm enhances the polar interactions between the two charged solutes; (3) For ext perpendicular to the direction of the electric dipole formed by S1 and S2, Eext = 100 MV/cm enhances the polar interactions between these two charged solutes; (4) The mean force and the PMF between two neutral solutes depend on ext.