Article pubs.acs.org/JPCC
Polarization Effect and Electric Potential Changes in the StimuliResponsive Molecular Monolayers Under an External Electric Field Jun Zhao,†,‡,◊ Xingyong Wang,†,◊ Nan Jiang,§ Tianying Yan,∥ Zigui Kan,† Paula M. Mendes,⊥ and Jing Ma*,† †
School of Chemistry and Chemical Engineering, Institute of Theoretical and Computational Chemistry, Key Laboratory of Mesoscopic Chemistry of MOE, Nanjing University, Nanjing 210093, P. R. China ‡ School of Physics and Optoelectronic Engineering, Yangtze University, Jingzhou, Hubei 434023, P. R. China § School of Pharmacy, Nanjing Medical University, Nanjing, Jiangsu 210029, P. R. China ∥ Institute of New Energy Material Chemistry, Collaborative Innovation Center of Chemical Science and Engineering (Tianjin), Nankai University, Tianjin 300071, P. R. China ⊥ School of Chemical Engineering, University of Birmingham, Edgbaston, Birmingham B15 2TT, U.K. S Supporting Information *
ABSTRACT: Under the exerted external electric field (Eex), significant polarization effects are demonstrated in the switching process of positively charged oligopeptide chains on the gold surface by using the polarizable force-field-based molecular dynamics (MD) simulations. The changes of electrostatic environment during the conformational switching are described by calling density functional theory (DFT) calculations of atomic charges or fragment-centered dipole moments. The charge-variable polarizable force field (Q-POL) model shows a good agreement with the experimental observations of both the open-circuit state (OC) without Eex and extended/bent conformation (ON/OFF) with upward/ downward Eex, while the nonpolarizable force field (Non-POL) with the fixed atomic partial charges fails to reproduce experiments for the ON state. The charged oligopeptide chains (with each residue furnished with a positive charge) challenge the application of the coarse-graining polarizable force field based on the fragment dipoles (Dfrag-POL). Through tracing the dynamics of charge centers of both counterions and oligopeptide chains in the switching process, a qualitative picture is depicted for understanding the polarization phenomena in the dielectric monolayer in OC/ON/OFF states. With the applied electric field, the rearrangement of ions leads to the induced internal E-field, Ein, in the direction opposite to Eex. The values of total electric field, Etotal, and the interfacial potential Φx were evaluated. The relationship between the applied electric potential and the interfacial potential is built for the switching OC/ON/OFF states. The analysis of the induced internal E-field and interfacial potential may shed light on the further modification of theoretical models to better understand the electrical-induced switching mechanism.
1. INTRODUCTION Electroswitchable surfaces are now attracting much research interest due to their potential applications in building smart materials for drug delivery, bioseparation, data storage, and microelectronic device fabrication.1−4 Under the electric field, a switchable surface was first built by Lahann and Langer1 to trigger the surface switching from hydrophilic to hydrophobic. Since then, an unprecedented ability to manipulate the interactions of peptides, DNA, proteins, and cells with surfaces has been achieved by employing electrical stimuli on selfassembled monolayers (SAMs) with a number of different electroresponsive groups.2,5−9 Among those switchable surfaces, a kind of SAM, consisting of the mixed long and short chains together on the substrate, is frequently used in the design of switchable surfaces,2,3,10−12 as shown in Figure 1(a). The long chains have electroswitchable units, either on the head or in the © 2015 American Chemical Society
side chain, which will respond to external electric perturbations. The short chains can not only keep the ordered SAM structure but also provide sufficient spatial freedom for conformational changes of the long chains. These advances in fabrication of smart surfaces have stimulated the interests of theoretical simulations and designs. By introducing the electric field effect into molecular dynamics (MD) simulations, Pei and Ma3 predicted electric-field-induced switching behaviors of a series of the mixed long (ω-carboxyalkyl) and short (alkyl) chains covered H−Si(111). The picture of constructing reversible switching SAMs with the mixed long/short chains was extended to biointerfaces to modulate protein binding to surfaces.2,12 The Received: May 20, 2015 Revised: September 15, 2015 Published: September 17, 2015 22866
DOI: 10.1021/acs.jpcc.5b04805 J. Phys. Chem. C 2015, 119, 22866−22881
Article
The Journal of Physical Chemistry C
Figure 1. (a) Schematic illustration of the switching of SAMs between extended and bent conformations, denoted as ON and OFF states here. (b) The surface models used in the MD simulations. The lengths of the long and short chain are obtained with the Non-POL model. The long purple chain represents the biotin-4KC chain, and the short red chains denote TEGT.
performed to investigate the conformational changes under the instantaneously switching surface potential.18 A simplified analytical model was also proposed to describe the dynamic behavior of oligonudeotide rods in DC and AC electric fields.19 In addition to the E-field-induced conformational changes, a large number of anions and cations moved toward two different surfaces at strong electric fields, and such nonequilibrium ion distributions induced an internal electric field, which was addressed in the study of polyelectrolyte brushes in the presence of external electric fields.20 The other interfacial electric properties, such as solvation effects on electronic polarization at the air−water interface,21−23 potential of zero charge of the electrode,24 electric double layers of ionic liquids,25 induced dipoles, and interfacial potential26 of the water/metal or ion/ solid interface under the applied electric field, have also attracted extensive attention. In comparison with the real experimental conditions, however, the electric field model used in the theoretical simulation was often simplified by a constant electric potential (using a homogeneous electric field)3,12,14,17,27,28 or a constant surface charge (modifying the surface with either negatively or positively charged groups).29−33 In fact, the electric field exerted on the smart surfaces is not homogeneous any more, as a result of the induced internal electric field, which will be analyzed in details in this work. In addition, the interfacial environment of the SAMs is often far more complicated than what has been analyzed for
electrical-induced conformational changes of the oligopeptide:ethylene glycol-terminated thiol mixed SAMs were investigated by us using the polarizable force field based MD simulations12 (Figure 1(b)). In addition, some other electric field responsive systems aroused intensive attention of theoretical chemists. Several different methods have been adopted to study the conformational changes in response to the applied external electric field. Using CHARMM force field based MD simulations and surfaceenhanced Raman spectroscopy, Liu and co-workers13 studied the external electric-field-induced conformational change of dodecapeptide probes tethered to a nanostructure metallic surface. The effects of applied electric field on the conformational transition of poly(ethylene glycol)-terminated alkanethiol SAM on gold were investigated by using a parallel molecular dynamics method.14 The combination of first-principle quantum mechanics (QM) and the mesoscale lattice Monte Carlo method was used to analyze the possibility of electric-field-based conformational changes and simulate the switching as a function of voltage and voltage rate.15 The surface-supported SAMs with azobenzene derivatives have also been demonstrated to switch between cis and trans conformations under applied external electric fields using the quantum chemical calculations16 or reactive molecular dynamics simulation.17 By treating the DNA as a chain of elastically connected beads with full hydrodynamic interactions between monomers, Brownian dynamic simulations were 22867
DOI: 10.1021/acs.jpcc.5b04805 J. Phys. Chem. C 2015, 119, 22866−22881
Article
The Journal of Physical Chemistry C Table 1. Optimized Geometry and Intramolecular Hydrogen Bonds
HB1 conformers a b c a
rO···H/Å 1.58 (1.57) 1.60 (1.60) 1.87 (1.82)
a
end-to-end distance
torsion angle
∠O−H−N/°
rO···H/Å
HB2 ∠O−H−N/°
dO23−S8/nm
φ/°
163.9 (159.6) 173.6 (169.8) 158.3 (154.7)
1.58 (1.56) 1.91 (1.96) 1.65 (1.78)
176.4 (169.4) 159.3 (152.2) 172.8 (176.2)
2.38(2.35) 2.23(2.21) 1.62(1.59)
61.1 (55.6) 60.2 (42.1) 76.5 (83.9)
Values outside and inside the parentheses are calculated at B3LYP and M06-2X functionals, respectively.
model is employed as in our previous work (ref 12). To model the gold substrates, five layers of gold atoms were cut from the Au(111) surface, and they were fixed during the simulations. The charge for the gold atoms was set to zero, and charge transfer from the SAMs to the surface was neglected. Two-dimensional rhombic periodic boundary conditions and a 34.60 × 34.60 × 77.42 Å3 slab model were applied throughout our simulations. To reproduce the experiment conditions, 2115 water molecules and 4 chloride ions were adopted to simulate the phosphate-buffered saline solution. All MD simulations, on the basis of both nonpolarizable and polarizable force fields (FFs), were performed in the canonical (NVT) ensemble. The temperature was controlled at 298 K by the Andersen thermostat.37 The equations of motion were integrated using the velocity Verlet algorithm with the time step of 1 fs. The Discover module in the Materials Studio package38 and the Amber39 were employed to run the MD simulations. 2.2. Nonpolarizable Model. First of all, the performance of MD simulations depends mainly on the FF selected, so we tested four kinds of FFs, cvff (consistent-valence FF), compass (condensed-phase optimized molecular potentials for atomistic simulation studies), pcff (polymer consistent FF), and Amber/ FF03 (the third-generation Amber force field) with the default parameters. During the MD runs, the atomic charges were fixed no matter how the conformation changed under the external electric field, thus these classical FF models with default parameters are called the nonpolarizable (Non-POL) models. It has been previously shown that the cvff performs best, according to the energy scan of the switching barrier of the biotin-4KC chain.12 To further test the performance of various FFs, we calculated the relative energy for three local minimum conformations of biotin-4KC, which were obtained through DFT optimizations. In order to test the performance of the density functionals, two kinds of DFT functionals, B3LYP and M06-2X, were selected to optimize the geometries, respectively. As shown in Table 1, with these two functionals, the superimposition of the optimized
water/metal or ion/solid systems. It is hence interesting to investigate the interfacial potential and the induced electric field of SAMs in response to the applied E-field. Herein, we present a systematic theoretical exploration of the switching mechanism of a recently reported electroswitchable oligopeptide surface12 by molecular dynamics (MD) simulations with polarization models. DFT calculations are adopted to obtain the atomic partial charges or fragment centered dipole moments at regular intervals during the simulation, and this can be regarded as a quasi on-the-fly polarization model. In comparison with the fragment-based polarizable FF (Dfrag-POL) models,34−36 our Q-POL method is also shown to be more suitable for treating systems with strong polarizations. This is of vital importance in the study of complex biological systems, where polarizations caused by either the aqueous environment or charged chemical groups are ubiquitous. Moreover, the internal induced E-field, which is often ignored during MD simulations, has also been investigated. The interfacial potential distributions have been plotted to display the difference between the experimental external electrical potential and the simulated external electrical potential. The relationship between the applied electrical potential and interfacial potential is employed to reach a reasonable understanding of the studied electricalinduced switching behavior of positively charged oligopeptide chains on a gold surface.
2. COMPUTATIONAL METHODS 2.1. Surface Model and MD Simulation Details. The nomenclature for the components of the studied system is shown in Figure 1b. Four positively charged lysine (K) residues worked as the switching unit, and a terminal cysteine (C) was used for attachment to the gold surface. They were represented by “4KC”. The short chain was tri(ethylene glycol)-terminated thiols, denoted as “TEGT”. The biotin was attached on the head of the 4KC long chain, with a purpose to bind with the neutravidin protein. The conformational switching of biotin-4KC was simulated by using molecular dynamics. The same slab surface 22868
DOI: 10.1021/acs.jpcc.5b04805 J. Phys. Chem. C 2015, 119, 22866−22881
Article
The Journal of Physical Chemistry C
polarization models have been developed. By treating the local change in charge density around an atom with induced dipoles or charged particles, the induced dipole models44−52 and the Drude oscillator models53−69 have been proposed for the study of a wide range of problems. The other kind of model is the fluctuating charge (FQ) method,70−97 in which the values of atomic charges are treated as dynamical variables in response to the polarization, which are derived from the principle of electronegativity equalization.98−106 The transferability of the force field parameters needs to be tested when those polarization models are used to describe the switching behaviors of the charged oligopeptide surfaces under an external electric field because most models do not consider the impacts of the external electric field when constructing force field parameters. Another kind of polarization model is based on semiempirical or ab initio quantum mechanical (QM) calculations.34−36,107−116 Taking advantage of linear-scaling fragment-based QM calculations and the multilayer coarse-graining model, several polarization models have been proposed to treat electrostatic interactions of solvated α-conotoxin peptides.34,36 In spite of these successes, few studies have been devoted to examining the polarization effects under electric fields. Here, we adopted two kinds of polarizable (POL) models, as schematically shown in Scheme 1. One polarizable model, QPOL, is based on the variable partial charges. The electrostatics energy, Eelectrostatics, is calculated according to eq 1
structures indicates both functionals give similar predictions on the geometry. It is noteworthy to point out that the intramolecular hydrogen bonds (labeled as HB1 and HB2 in Table 1) exist in three local minima conformers a−c. Obviously, the folding of the structures in conformers a−c may be ascribed to the intramolecular hydrogen bonding interactions between two O atoms in biotin and N−H bonds of NH3+ groups in the third lysine (lys3) and the forth lysine (lys4), respectively. The geometry parameters of the hydrogen bonds HB1 and HB2 are listed in Table 1. To better illustrate the conformational changes between the folded (OFF) and fully stretched (ON) states, the end-to-end distances, d, between atom O23 (of biotin) and S8 (of thiol) and the torsion angles φ of C48−C57−C59−N61 (around the linkage of lys3 and lys4 residues) are also given in Table 1. Among the three conformers, the bent conformer c has the shortest end-to-end distance, d = 1.62 nm, far smaller than that (4.1 nm) for the fully stretched state. As shown in Table 2, the relative energies predicted by the two DFT functionals for the three local minima a−c are also in the Table 2. Relative Energies for Three Conformations a−c Calculated with Different Methodsa method E kcal mol−1 DFT B3LYP/6-31G(d) M06-2X/6-31G(d) Non-POL-FF cvff pcff compass Amber/FF03 POL-FF Q-POL-cvff(NBO)b Q-POL-cvff(Mulliken)b Q-POL-cvff(MK-ESP)b Q-POL-cvff(CHelpG-ESP)b Q-POL-Amber/FF03c Dfrag-POL-compass (Frag I)d Dfrag-POL-compass (Frag II)d Dfrag-POL-Amber/FF03 (Frag I)c Dfrag-POL-Amber/FF03 (Frag II)c
conformers a
b
c
0 0
4.56 2.76
7.56 10.38
0 0 0 0
8.04 3.91 −0.68 19.93
14.51 15.31 14.60 33.00
0 0 0 0 0 0 0 0 0
6.02 8.64 9.78 12.51 18. 10 6.88 12.04 16.78 16.57
12.79 51.13 49.18 23.05 54.30 25.68 23.93 51.09 51.01
Eelectrostatics =
∑ i,j
qi q
j
4πε0rij
(1)
where rij is the interatomic distance and qi, qj are variable charges of the ith and jth atoms, respectively, which are calculated from DFT on the dynamically charged conformations sampled from MD simulations. In the present Q-POL polarizable model, atomic charges for the positively charged biotin-4KC molecule were obtained by DFT calculations at the M06-2X/6-31G(d,p) level of theory. All DFT calculations were carried out with the Gaussian 09 program package.117 Different kinds of population analysis, such as Mulliken, NBO (natural bond orbital), MK-ESP (MK-electrostatic-potential), and CHelpG-ESP, were employed to estimate the atomic charges for the biotin-4KC molecular chain. Inspecting from Table 1, we can find that the relative energies of conformers a−c obtained with NBO charges are in the best agreement with the DFT results. Thus, in our following MD simulations, the Q-POL-NBO method will be utilized to investigate the switching mechanism and the polarization effects. All the atomic charges calculated from DFT are listed in Table S1 in the Supporting Information. As shown in Scheme 1(a), our strategy involves several steps. In the first step, geometry optimization was performed for a single biotin-4KC molecule at the M06-2X/6-31G(d,p) level of theory. Then the NBO charge for each atom was calculated at the optimized structure by M06-2X/6-31G(d,p). It should be mentioned that when biotin-4KC was tethered to the gold surface the hydrogen atom in the thiol end would be removed. Herein we tackle this problem by dispersing the charge of this hydrogen atom (q(Hend)) among all the other atoms in biotin4KC. Assuming the initially calculated charge of the ith atom in biotin-4KC except Hend is qinit(i), then the corresponding dispersed charge will be qinit(i) + q(Hend)/N, where N is the number of atoms in biotin-4KC except Hend. Then the atomic partial charges for biotin-4KC in MD simulation were replaced
a
The structures of a−c in the calculations by FF methods were from geometry optimizations at M06-2X/6-31G(d). bThe natural bond orbital (NBO), Mulliken, and electrostatic potential (ESP) atomic charges calculated at M06-2X/6-31G(d) were adopted in the Q-POL calculations. cThe scale factor of the electrostatic parameters was set to 0.10. dThe scale factor of the electrostatic parameters was set to 0.50.
same order. Both B3LYP and M062X functionals show the minimum a is more stable than the other two in vacuum. Almost all the selected force fields give the same order of the relative energies as DFT. 2.3. Polarizable Model. Classical force field (FF) based molecular simulation has been widely used in the modeling of both biomolecules and nanomaterials and can provide atomic details that are often inaccessible in experiments. In conventional FF methods, polarization effects cannot be considered well since the electrostatic potential is simply expressed as a sum of pairwise Coulombic interactions between the atom-centered point charges with fixed parameters.40−43 However, polarization effects cannot be neglected in some cases. Until now, a great variety of 22869
DOI: 10.1021/acs.jpcc.5b04805 J. Phys. Chem. C 2015, 119, 22866−22881
Article
The Journal of Physical Chemistry C Scheme 1. Flowchart of the Procedure of the (a) Q-POL and (b) Dfrag-POL Polarization Modelsa
a
The magnitude of the dipole moments for conformers a−c are listed with the green arrows denoting the directions of the dipole moment. HB1 and HB2 refer to two intramolecular hydrogen bonds, respectively.
by these NBO charges. To validate this strategy, we adopted a simple gold cluster model having three gold atoms connected to biotin-4KC (Figure S1 in the Supporting Information) and calculated the NBO charges for it. It was shown that the three gold atoms had a total charge of +0.32 e, which was insignificant compared to the total charge of the system (+4 e). We also compared the partial charges for every atom obtained by both methods, as shown in Table S2. The largest difference was less than 0.08 e. Therefore, considering the computational cost, we adopted this strategy. The same charge generation procedure was done for TEGT, and the atomic charges for TEGT were fixed to the initial values throughout the MD simulation since TEGT was neutral and the atomic charges did not change significantly when applying the electric fields with different strength and directions. For the charged oligopeptide monolayer on the gold surface, some independent Q-POL MD runs have been performed, with the partial charges renewed every 10, 50, and 100 ps, respectively. No distinct differences are observed about conformational changes (end-to-end distance: Figure 2a; torsion angle: Figure S2a). In subsection 3.1.2, it will also be shown that the different
charge updating periods lead to little changes in the charge fluctuations along time evolution. In Figure 2b and Figure S2b, the ultrafine charge-updating interval of 0.2 ps was adopted for an isolated biotin-4KC chain. For the studied complex SAM system with thousands of atoms, a compromise has been made between computational costs and accuracy: the atomic partial charges will be updated every 100 ps in the following MD simulations for QPOL models. Starting from the initial conformations as shown in Table 1, energy minimization was carried out, followed by a long enough MD simulation after which the system reached equilibrium. We have carried out nearly 1 ns presimulations for Non-POL and QPOL models. Subsequent 400 ps simulations without and with upward/downward external electric field have been then carried out to collect data of OC/ON/OFF states for Non-POL and QPOL models, respectively. The OC condition was done by prolonging a 400 ps simulation with no electric field applied after equilibrium. Next, the electric fields in both +z and −z directions were applied to model the positive and negative surface potentials in the experiment. During all the simulations, the 22870
DOI: 10.1021/acs.jpcc.5b04805 J. Phys. Chem. C 2015, 119, 22866−22881
Article
The Journal of Physical Chemistry C
directions (designated as μX, μY, and μZ, respectively) are represented by
μX =
∑ qiXi i
μY =
∑ qiYi i
μZ =
∑ qiZi i
atomic charges for biotin-4KC were updated by DFT calculations every 100 ps. This was achieved by a Fortran program that interfaced Gaussian and Discover in this procedure: (i) extract the atomic coordinates of biotin-4KC from the structural file after 100 ps simulation; (ii) add a hydrogen atom with the coordinate of the Au atom that tethers the thiol of botin4KC and generate a Gaussian gif file; (iii) carry out DFT calculations and obtain NBO charges; (iv) disperse the charge of the added hydrogen atom among all the other atoms in biotin4KC; (v) insert the updated atomic charges into the Discover input files; (vi) perform for the next 100 ps MD simulation. These six steps were then repeated until the simulation was completed (Scheme 1(a)). The detailed information on MD simulation process for Non-POL and Q-POL models was illustrated in Figure S3(a) of the Supporting Information. It is also interesting to make a comparison with the recently developed fragment-based method (Dfrag-POL), which uses the fragment-centered dipole moments to calculate dipole−dipole interactions.34−36 Dfrag-POL was originally proposed to treat the large-sized systems by cutting the whole systems into several fragments and describing the electrostatical energy in terms of 118−120 the fragment-centered dipole moments μFrag The diI . μμ pole−dipole electrostatic energy Eelectrostatics is calculated according to eq 2
∑k I
