7548
J. Phys. Chem. B 2005, 109, 7548-7556
Molecular Dynamics Simulations of Helix-Forming, Amine-Functionalized m-Poly(phenyleneethynylene)s Bamidele Adisa and David A. Bruce* Department of Chemical Engineering, Clemson UniVersity, Clemson, South Carolina 29634-0909 ReceiVed: NoVember 11, 2004; In Final Form: February 8, 2005
We present here the results of all-atom and united-atom molecular dynamics (MD) simulations that were used to examine the folding behavior of an amine-functionalized m-poly(phenyleneethynylene) (m-PPE) oligomer in aqueous environment. The parallelized GROMACS MD simulation code and OPLS force field were used for multiple MD simulations of m-PPE oligomers containing 24 phenyl rings in extended, coiled and helix conformations separately in water to determine the minimum energy conformation of the oligomer in aqueous solvent and what interactions are most important in determining this structure. Simulation results showed that the helix is the preferred minimum energy conformation of a single oligomer in water and that Lennard-Jones interactions are the dominant forces for the stabilization of the helix. In addition, these solvophobic interactions are strong enough to maintain the helix conformation at temperatures up to 523 K.
Introduction Numerous heteroatom-containing polymers, especially biobased polymers and proteins, can exhibit a wide range of physical properties and chemical reactivity depending on the conformation adopted.1,2 The formation of a specific polymer conformer is governed by the system temperature and the presence of numerous intermolecular interactions of varying types (e.g., hydrogen bonding, van der Waals, Coulombic, etc.).3 Studies have also shown that the hydrophobicity of side groups greatly influences the nature and extent of polymer folding.2-5 To control or predict the properties and reactivity of these systems, one must first be able to predict the conformation that will be adopted for a given sequence of monomers. In an effort to better understand the nature of helix formation, which is one of several conformations that are commonly adopted by bio-based polymers and protein segments, Moore and others2,6 examined a simpler system of polymers, m-poly(phenyleneethynylene)s or m-PPEs (example shown in Figure 1), where a range of side group hydrophobicities were examined. These oligomers were observed to undergo a solvent-driven transition from a random conformation in chloroform to a helical conformation in acetonitrile. On the basis of their results, they concluded that solvophobic interactions are the dominant force in helix stabilization. Additional m-PPE oligomers with varying functional groups were synthesized by Arnt and Tew.7 They investigated the ability of these m-PPEs to adopt amphiphilic conformations at an air-water interface. Although, they did not perform detailed studies on the folding behavior of these materials in solution, they did not rule out the possibility of these materials adopting helical conformations.7 However, the observed emission spectra differed from those reported by Prince et al.6 for helical m-PPEs in solution. Arnt and Tew also performed experiments to investigate the conformational changes of amphiphilic m-PPEs in aqueous solution.8 Their results indicated that m-PPEs without alkyl substituents adopted helical conformations in aqueous solution, whereas m-PPEs having long alkyl substituents were more likely to take on nonhelical conformations. * Corresponding author contact information. E-mail:
[email protected]. Phone: 864-656-5425. Fax: 864-656-0784.
Figure 1. Monomeric unit of the m-PPE oligomer. C5 is connected to C21 of the previous unit, and C21 is connected to C5 of the next unit. R is n-C4H9.
Further studies examining the folding behavior of m-PPEs functionalized with highly polar or ionic groups have indicated that the m-PPE oligomers adopt a helical conformation in highly polar solvents (e.g., water) and a more open or random-coil structure in less polar solvents (e.g., alcohols).9,10 Specifically, Tan et al.9 showed via spectroscopic measurements that sulfite functionalized m-PPE oligomers (m-PPE-SO3), which are conjugated polyelectrolytes, adopt a helical conformation in water and a random-coil conformation in methanol. In addition, Inuoye et al.10 have recently reported that poly- and oligo(methynylpyridine)s are driven into helical conformations by hydrogen-bonding interactions with saccharides, with the helical structures having the same chirality as the bound saccharides. Furthermore, experimental studies have also shown that arylene ethynylene strands self-assembled into double-helical structures when pyridine was added.11 In recent years, attention has also been drawn to ortho-linked PPEs (o-PPEs). There have been computational12 and experimental13-15 studies performed on o-PPEs and the results of these studies suggest that the ortho-linked polymers are also capable
10.1021/jp0407122 CCC: $30.25 © 2005 American Chemical Society Published on Web 03/22/2005
MD Simulations of m-Poly(phenyleneethynylene)s of adopting helical conformations in solvents of appropriate polarity. However, these systems are more likely to experience ring strain and increased steric interactions upon folding into helical structures. Such interactions tend to reduce the free energy of folding, thus making helical conformations less favored. Though the experimental work of Moore, Tew and others2,6-8 clearly demonstrated that solvophobic interactions influence secondary structure formation in m-PPEs, they did not provide a detailed atomic-level description of the intramolecular repulsions or intermolecular attractions that drive polymer systems to take on particular conformations. In general, the elucidation of the significance of each of these phenomena experimentally is difficult, hence, the need for atomic-scale computational studies of these folding phenomena. In an effort to better understand the interatomic forces influencing helix formation, this paper focuses on the results of explicit atom (including solvent) and united atom MD simulations that examined the preferred conformations of m-PPE oligomers having alkyl ether and amine functional groups attached to the oligomer backbone. Monte Carlo (MC) and molecular dynamics (MD) techniques are powerful tools for describing complex chemical systems at the atomic scale. These force field-based techniques provide an effective means for understanding and predicting macroscopic behavior that develops at time scales (1-100 ns) beyond those easily sampled by quantum computational methods, including density functional theory (DFT) methods.16 Additionally, coarse grained (e.g., bead and spring or lattice) models,16,17 which are capable of predicting even longer time scale events (1-10 µs), do not provide explicit representations of all atoms in the system; hence, they are unable to describe the atomic-level interactions that determine the potentially chiral 3-D conformation of a complicated polymer system. A limited number of computational studies have examined the folding behavior of m-PPE oligomers. Nelson et al.2 conducted computational studies on a series of m-PPEs, carboxylic acid-functionalized oligomers in water and methyl benzoate analogues in chloroform, using MD simulations that incorporated a modified version of the Optimized Potentials for Liquid Simulation (OPLS).18 Their results, which were corroborated by experimental studies, showed that helix formation occurred with chain lengths greater than 7 monomer units. Similarly, Elmer and Pande19 performed a series of MD simulations on a 12-mer m-PPE oligomer using the NAMD220 MD code, the CHARMM1921 parameter set, and an implicit solvent model. Their results demonstrated that the m-PPE folds into a helix via a series of intermediate partially coiled states. Sen22 also conducted MD simulation studies of helix formation of m-PPE oligomers having ester functional groups using the CHARMM-2621,23 simulation package. His models illustrated that helical conformations were favored over random-coil structures when the oligomer was surrounded by TIP3P24 water molecules. More recently, Lee and Saven25 have reported results of simulations conducted using NAMD220 with the CHARMM-1921 force field on a helical m-PPE oligomer in aqueous solution. The m-PPE oligomer was found to maintain its helical structure throughout the simulation. All of the aforementioned MD simulations showed that m-PPE oligomers with carboxylic acid and ester side groups form stable helical structures in aqueous solvents, but no studies have examined m-PPE oligomers with amine side groups. Further, some of the previous work involved the use of implicit solvents; thus, it was impossible to gain insights into the ordering of solvent molecules around the polymer and the influence that this might have on helix stabilization. Because proteins are known to contain amine functional groups, the amine-terminated
J. Phys. Chem. B, Vol. 109, No. 15, 2005 7549 oligomers modeled in this study, which are identical to those studied by Arnt and Tew,7,8 provide a simple system for further investigating the various interactions that influence the formation of extended, stabilized secondary structures in aqueous media. The chemical structure of the monomeric unit of the oligomer used in this study is shown in Figure 1. The MD studies described in this work used the Groningen Machine for Chemical Simulation (GROMACS, version 3.2) software package.26-28 The GROMACS software employs a Message Passing Interface (MPI) parallelization method, which enables it to be used on cluster-type supercomputers (e.g., Beouwulf clusters). This algorithm couples a simple onedimensional domain decomposition method with a cyclical calculation procedure.27 Although it was developed for biological molecules (e.g., proteins and lipids) with complex bonded interactions, GROMACS is still suitable for MD simulations of nonbiological polymer systems because of its efficient calculation of nonbonded interactions. Its advanced algorithm and optimized coding also make it the fastest MD simulation package available.28 Because the number of nonbonded interactions scales with the square of the number of interaction sites present,16 considerable computational time can be saved if the number of interaction sites can be reduced. Hence, we examined a unitedatom model, OPLS-UA,18 as well as an all-atom model, OPLSAA,29 and compared the results so as to determine if the united atom model can adequately describe the interatomic interactions affecting structure formation. At present, computational speeds do not easily allow for the performance of MD simulations longer than tens of nanoseconds, which is significantly shorter than the time scale (microseconds) needed for longer chain polymers to fold into ordered secondary structures, such as the helix structure.22 To overcome this shortfall, a procedure similar to that implemented by Sen22 was employed. With this method, multiple MD simulations were performed using m-PPE oligomers containing 24 phenyl rings (12 residues, where a residue includes 2 phenyl rings) in extended, coiled and helix conformations separately in water at 300 K, and the results from each of these simulations were compared. The results were analyzed to determine what conformation the m-PPE oligomer prefers in aqueous solvent and what interactions were most important in the stabilization of the helix conformation. In addition, MD simulations were performed on the m-PPE oligomers in the helix conformation at different temperatures to determine the thermal stability of the helix structure. Computational Details Force Field and Parametrization. The system potential energy was calculated using the GROMACS implementation of the OPLS-AA force field.18,26-31 With this force field, eq 1, the system potential energy is a classical function of the internal degrees of freedom and the nonbonded distances between atomcentered sites.
EPE )
kr(r - r0)2 + ∑ kθ(θ - θ0)2 + ∑ bonds angles 5
∑ kξ(ξ - ξ0)2 + n)0 ∑Cn(cos(ψ))n + impropers fij ∑i ∑ j>i
{
qiqje2 rij
+ 4ij
[( ) ( ) ]} σij rij
12
-
σij rij
6
(1)
In the first three bonded energy terms, the force constants for intramolecular deformations (kr, kθ, kξ) define the magnitude
7550 J. Phys. Chem. B, Vol. 109, No. 15, 2005 of the energy required to move the internal coordinates (r, θ, ξ) away from their unstrained default values (r0, θ0, ξ0). These nontorsional bonded interactions are modeled by harmonic terms for bond stretching, angle bending, and out of plane deformations for planar groups. The fourth term defines the proper torsions in terms of the specific dihedral angle (ψ) and RyckaertBellemans potential parameter Cn, where n ) 0, 1, ..., 5.32 The final term in eq 1 quantifies the nonbonded energy associated with all intermolecular pairs of sites, i and j, separated by a distance rij. These nonbonded interactions are modeled by Coulombic and 6-12 Lennard-Jones terms, where q is the partial atomic charge, and σ and are the Lennard-Jones parameters. Standard combination rules are used such that σij ) (σiiσjj)1/2 and ij ) (iijj)1/2. The nonbonded energy term is evaluated for all intermolecular interactions as well as for all intramolecular interactions between sites separated by three or more bonds. The scaling factor fij is 1.0 for all nonbonded interactions except for the 1,4-intramolecular interactions, where fij is 0.5 to permit the use of the same parameters for both the intermolecular and intramolecular interactions.29 Atom types were assigned using the default OPLS-AA parameter set for both the united-atom and all-atom models, except for the aliphatic CH2 and CH3 groups, where implicit hydrogens were assumed for the united-atom model and parameters were obtained from the OPLS-UA parameter set.18 Partial atomic charges were initially assigned using the default OPLS18,29 charges for appropriate atom types. However, these charges were modified to provide closer agreement with those obtained from quantum mechanical (QM) energy minimized structures. The QM minimization was performed using Jaguar 4.0 software,33 run on an SGI Onyx2 InfiniteReality2 Dual Rack, and employed density functional theory (DFT) using the hybrid B3LYP34,35 density functional and the 6-31G** 36 split valence basis set. Partial atomic charges were obtained by fitting the electrostatic potential (ESP) to each atom center and constraining it to replicate the molecular dipole. Generation of Starting Structures. The m-PPE oligomers containing 24 phenyl rings were manually generated in the extended, coiled and helix conformations using Materials Studio 2.1,37 and the resulting output coordinate files were modified to make them compatible with GROMACS. The helical structure was generated with 6 phenyl rings per turn of the helix. The resulting structures were energy-minimized using steepest descents method with a maximum step size of 0.01 nm and tolerance of 10 kJ‚mol-1‚nm-1. This was done to remove strain in the polymer backbone and alleviate high energy close contacts. The energy-minimized structures were then solvated with SPC38 water to a density of 1 kg‚L-1. Solvation39 of m-PPE oligomers was achieved by placing each oligomer at the center of a rectangular periodic box of SPC38 water molecules where the distance between the solute and the edge of the box was at least 1.5 nm. The solvated system was initially energyminimized using steepest descents method, and then minimized further using the L-BFGS method of Nocedal.40 This minimization protocol helped to remove any close contacts between the oligomer and the surrounding water molecules. The resulting periodic systems were the starting configurations for all MD simulations. MD Simulation. All MD simulations were done using the GROMACS 3.2 simulation package. In the simulations, the leapfrog algorithm41,42 was used to integrate Newton’s equations of motion with a time step of 2 fs. Periodic boundary conditions were applied and nonbonded force calculations employed a grid system for neighbor searching. In this system, only the atoms
Adisa and Bruce TABLE 1: Radius of Gyration (Rg) for the Helix, Extended and Coiled Conformations of the m-PPE Oligomer, Averaged over the Last 1 ns of MD Trajectory Rg (nm) structure coileda extended helixa
all-atom
united-atom
1.00 ( 0.06 1.09 ( 0.20b 0.83 ( 0.01
1.00 ( 0.03 1.00 ( 0.15a 0.83 ( 0.01
a Standard deviation obtained from 1 simulation. b Standard deviation average for 16 simulations.
in the neighboring grid cells are considered when building a new neighbor list. Neighbor list generation was performed after every 5 steps. A twin-range27 cutoff was used for both LennardJones and Coulombic calculations. A cutoff radius of 0.9 nm was used for short-range forces, which were calculated every simulation step, and a cutoff radius of 1.4 nm was employed for long-range forces, which were calculated during neighbor list generation. In each simulation, the temperature was controlled by employing a weak coupling to a heat bath, with a coupling constant of 0.1 ps.43 The m-PPE oligomer and water were independently coupled to the heat bath. Initial velocities were randomly assigned from a Maxwell distribution at the selected simulation temperature. The LINCS algorithm44 was used to constrain all bonds in nonwater molecules, and the SETTLE algorithm45 was used to constrain all water bond lengths and angles to their equilibrium values. All MD simulations were continued for 5 ns and atomic trajectories were saved every 5000 steps (10 ps intervals) for analysis. To have a statistical sampling of the behavior of an extended polymer structure in water, four different all-atom extended structures were simulated using different initial velocity vectors of the atoms to give a total of 16 simulations. Most simulations were conducted at 300 K; however, dynamics of the all-atom helical oligomer system, using a time step of 1 fs, was also studied at 348, 373, 423, 473, and 523 K. Results and Discussion The default OPLS-AA/UA18,29 bond stretching and angle bending parameters were used as they were in close agreement with those of other popular force fields21,23-28,46,47 and predicted structures having geometries similar to those observed by X-ray diffraction. Similarly, the default values for the torsional energy and nonbonded interaction terms were used because they had been optimized for the OPLS18,29 force field. The OPLS-AA/ UA torsional parameters were translated into Ryckaert-Bellemans32 parameters to make them compatible with GROMACS. For the partial atomic charges, the default OPLS-AA/UA charges were used for most of the atoms; however, some charges were adjusted to have neutral subunits or charge groups and also have them be in close agreement with the charges obtained from the quantum energy minimization (see Table 1S in Supporting Information). These adjustments were in the range of (0.05 e. It is important to state that these adjustments in charges had to be kept small else reparametrization of all nonbonded and torsional energy terms would have been required because these terms were parametrized simultaneously.29 Non-OPLS dihedral terms were also added to account for conjugation along the polymer backbone. Experimental and theoretical studies48-51 have shown that the barrier for rotation in aryl groups linearly conjugated by alkyne linkages is very low; hence, the orientation of the aromatic rings relative to one another varies between coplanarity and orthogonality. In this
MD Simulations of m-Poly(phenyleneethynylene)s
J. Phys. Chem. B, Vol. 109, No. 15, 2005 7551
Figure 2. Evolution of the radius of gyration (Rg) for the all-atom helix, coiled, and highest (L11 and L12) and lowest (L5 and L16) energy extended structures of the m-PPE oligomer.
kJ‚mol-1,
study, we used a torsional barrier of 2.3924 which was obtained experimentally for tolane by Okuyama et al.,51 for the C12-C16-C19-C20 (ψ1) and C23-C21-C7′-C8′ (ψ2) dihedrals. We also viewed the trajectories of the helical m-PPEs using VMD52 and observed some flexibility in these dihedrals while the global helical structure was still maintained. This is in agreement with the results observed by Lee and Saven.25 The dihedral angles ψ1 and ψ2 for the all-atom helical structure averaged over the last 1 ns of trajectory were 6.76 ( 3.53° and 6.92 ( 3.65°, respectively, and ψ1 and ψ2 for the united-atom helical structure were 4.07 ( 3.35° and 10.39 ( 3.48°, respectively. We observed that the flexing of these dihedrals caused the height of the helix to vary slightly throughout the simulation whereas the diameter of the helix remained relatively constant. The results of m-PPE MD simulations were analyzed to determine the mechanism for secondary (helical) structure formation as well as to quantify the energy barriers and time scale for structural transformations. Several methodologies for analyzing atomic trajectories provide insights into structure formation. For example, the radius of gyration and solvent accessible surface area provide a measure of the degree of compactness of the polymer, whereas the radial distribution function for the oxygen atoms provides information about the presence of helical structure units. All of these results indicate that the amine-functionalized oligomer prefers to adopt compact helical conformations that reduce the interaction of the alkyl side chains with the polar solvent (water), while maximizing the favorable interactions between the primary amine groups and water. The radius of gyration, Rg, which provides a measure of oligomer compactness, can be calculated as follows:
Rg )
(∑∑ ) 2 i||ri|| mi imi
Figure 3. Evolution of the root-mean-square deviation (RMSD) for the all-atom helix, coiled, and highest (L11 and L12) and lowest (L5 and L16) energy extended structures of the m-PPE oligomer with reference to their initial energy-minimized structures.
(L5 and L16) energy extended structures. It can be seen that the Rg for the helix structure remains relatively constant over the course of the MD simulation, indicating that it did not change significantly from its initial energy-minimized structure. However, the Rg for the coiled structure shows a steady decrease over the first 2.5 ns of trajectory until the structure becomes trapped in a local energy minima, yielding a steady value for Rg that is above that of the helix. Similarly, the Rg for the L16 structure shows a steady decrease over the first 3.5 ns of trajectory before reaching an average steady value. The Rg profiles for the L5, L11 and L12 structures show several peaks and valleys before becoming relatively steady after 4 ns. The fluctuations in Rg for these systems are indicative of the end segments of the oligomer having to become extended so that it can recollapse into a lower energy structure that exhibits greater interaction between phenyl rings (i.e., a more helical structure). It is important to note that the average Rg value for the helix is lower than those for the coiled and extended structures, see Table 1. The relaxation times for the m-PPE oligomer in the helix, coiled and extended conformations were compared by performing root-mean-squared deviation (RMSD) on the coordinates of the final structures at each time frame of the trajectories with respect to the initial energy-minimized structures. The equation used to calculate the RMSD can be expressed as follows:
RMSD(t,t0) ) N
1/2
(2)
where mi and ri are, respectively, the mass of atom i and the position of atom i with respect to the center of mass of the molecule. Figure 2 shows the evolution of the Rg over time for the all-atom helix, coiled, and highest (L11 and L12) and lowest
[ ∑ || 1
N
M i)1
mi ri(t) - ri(t0)
|| ]
1/2
2
(3)
where M ) ∑ mi and ri(t) and ri(t0) are the positions of atom i)1 i at time t and at the start of the simulation (t0), respectively. Figure 3 shows the evolution of the RMSD over time for the all-atom helix, coiled, and highest (L11 and L12) and lowest (L5 and L16) energy extended structures. It indicates that the RMSD remained constant over the course of the trajectory for the helix structure, whereas the RMSD for the coiled structure reached a relatively constant value after 2.5 ns. In addition, the
7552 J. Phys. Chem. B, Vol. 109, No. 15, 2005
Adisa and Bruce
TABLE 2: Root Mean Square Deviation (RMSD) for the Helix, Extended and Coiled Conformations of the m-PPE Oligomer, Averaged over the Last 1 ns of MD Trajectory RMSD (nm) structure a
coiled extended helixa
all-atom
united-atom
1.10 ( 0.11 1.28 ( 0.17b 0.40 ( 0.04
0.79 ( 0.07 1.05 ( 0.13a 0.39 ( 0.04
a Standard deviation obtained from 1 simulation. b Standard deviation average for 16 simulations.
L16 extended structure showed lower RMSD variability after about 3.5 ns, whereas the RMSD profiles for the L5, L11 and L12 extended structures show several peaks and valleys that correspond to those observed in the Rg profiles before becoming relatively constant after 4 ns. It is also important to note that the average RMSD value for the helix is considerably lower than those for the coiled and extended conformations, see Table 2. Furthermore, the low variability of the RMSD for the helix throughout the trajectory shows that it did not change considerably from its initial energy-minimized structure, thus indicating that it represents a minimum energy conformation of the m-PPE oligomer in aqueous solution. Similarly, Figure 4 shows the evolution of the potential energies over time for the all-atom helix, coiled, and highest (L11 and L12) and lowest (L5 and L16) energy extended structures. It shows the gradual lowering of the potential energies of the coiled and extended structures over the first 2.5 and 4.0 ns of trajectory, respectively, after which steady-state values are attained, indicating that most of the conformational irregularities present in the initial energy-minimized structures have been removed. Similar trends were observed with simulations of the united-atom structures. For these calculations, the potential energy is defined as the intramolecular potential energy of the m-PPE oligomer and it includes all bonded (angle, proper and improper dihedrals) and nonbonded (Lennard-Jones and Coulombic) interactions between all the atoms of the molecule.
Figure 4. Evolution of the potential energy for the all-atom helix, coiled, and highest (L11 and L12) and lowest (L5 and L16) energy extended structures of the m-PPE oligomer.
Figure 5. Top (a) and side (b) views for the energy-minimized helix. Top (c) and side (d) views for the post-simulation all-atom m-PPE helix. Gray, white and black represent carbon, oxygen and nitrogen atoms, respectively. Hydrogen atoms are omitted for clarity.
Figure 6. (a) Representative picture of the m-PPE oligomer in the coiled conformation after energy minimization. (b) Post-simulation picture of the m-PPE oligomer in the coiled conformation. Only the oligomer backbones are shown for clarity.
The bond energies are not included because constraints were applied to all bonds. We can see from Figure 4 that the potential energy of the helix structure is lower than the potential energies of the coiled and extended structures. This suggests that the helix structure is the lowest energy conformation for the m-PPE oligomer in aqueous solution, whereas the less compact structures resulting from extended MD simulations of the coiled and extended structures are the result of these oligomers becoming trapped in local energy minima. However, it is important to state that these results show that a helical conformation is favored when a single m-PPE is placed in water, which is analogous to an infinitely dilute solution. It is entirely possible that at higher concentrations of the m-PPE in aqueous solution, different conformations and extended assemblies of the oligomers are observed due to the steric interactions between the alkyl side groups in the helix interior. Hence, future simulation studies focused on the prediction of minimum energy conformations for dissolved polymers should examine whether the presence of other oligomers influences the conformation of each individual oligomer and consequently, the assemblies of these structures. The initial energy-minimized and final structures of the helix and coiled oligomers are shown in Figures 5 and 6, respectively. It can be seen in the post-simulation helical structure that the phenyl rings exhibited a staggered vertical alignment, thus deviating from the ideal model in which the phenyl rings of adjacent turns were perfectly aligned. Similarly, the extended structures formed compact structures that were similar to the post-simulation helix structure, in which the nonpolar alkyl
MD Simulations of m-Poly(phenyleneethynylene)s
J. Phys. Chem. B, Vol. 109, No. 15, 2005 7553
TABLE 3: Components of the Potential Energy of the m-PPE Oligomer in the Helix, Extended and Coiled Conformations, Averaged over the Last 1 ns of MD Trajectory structure
anglec
dihedralc
coiled extendedb helixa
887.9 ( 40.6 889.2 ( 47.9 896.6 ( 45.7
1189.9 ( 15.9 1181.2 ( 15.8 1180.1 ( 21.3
coileda extendeda helixa
530.8 ( 38.8 534.9 ( 34.0 533.2 ( 37.5
329.3 ( 11.4 331.4 ( 9.9 325.0 ( 10.6
a
improper dihedralsc
Lennard-Jonesc
Coulombicc
potential energyc
All-Atom 232.6 ( 79.3 228.3 ( 89.7 230.5 ( 22.0
-167.9 ( 102.9 -157.1 ( 114.9 -410.3 ( 31.7
-576.6 ( 11.4 -583.3 ( 11.4 -584.2 ( 10.4
1566.4 ( 147.3 1558.3 ( 174.3 1312.7 ( 55.5
United-Atom 240.2 ( 64.4 232.8 ( 54.4 237.4 ( 20.7
-95.8 ( 68.6 -117.3 ( 98.6 -353.6 ( 35.1
-392.6 ( 9.1 -373.5 ( 9.7 -348.0 ( 10.1
611.8 ( 122.9 608.3 ( 125.9 394.0 ( 51.3
Standard deviation obtained from 1 simulation. b Standard deviation average for 16 simulations. c Units are kJ‚mol-1.
groups were in the inner surface and the polar amine groups were exposed on the outer surface. From Figure 7, it can be seen that extended structure L5 formed three complete loops of the helix, whereas structures L12 and L16 each formed 2 complete loops of the helix after 5 ns of simulation. These results suggest that the average time scale for an extended structure to form a well ordered helix is greater than 5 ns. The length of time needed to form an ordered helix is primarily governed by the probability of the oligomer having sufficient energy to uncoil from a lower energy, compact structure into a higher energy extended structure that could then reorganize into a preferred helical conformation. From a simulation standpoint, this energy impetus may be provided by increasing and then decreasing the temperature of the simulation (simulated annealing) or by a replica exchange process (e.g., replica exchange MD53) that involves similar systems of atoms being modeled at a range of temperatures. To determine which type of interactions affected helix stabilization, a comparative study of the energetics of the m-PPE structures was done. Table 3 shows the different components of the average potential energies of the oligomer in the helix, coiled and extended conformations. It can be seen that the bonded interaction terms in all three conformations are quite similar. This indicates that in the transition of an extended or coiled structure to the helix conformation, the bond angles and dihedrals in the oligomer undergo minimal deformation. Similarly, the Coulombic interactions are similar for the three conformations in the all-atom model, whereas the helix structure has a slightly higher Coulombic energy than those for the coiled and extended structures in the united-atom model. However, the Lennard-Jones interactions lead to the helix conformation
Figure 7. (a) Representative picture of the m-PPE oligomer in the extended conformation before energy minimization. (b) Post-simulation pictures of the highest (L11 and L12) and lowest (L5 and L16) energy extended structures. Only the oligomer backbones are shown for clarity.
being strongly favored over the coiled and extended conformations in both the all-atom and united-atom models; hence, it can be seen that optimization of the Lennard-Jones interactions plays a dominant role in the stabilization of the helix, which is in agreement with published results.2,22 The importance of the Lennard-Jones interactions is also seen in Figure 8, where the bonded and Coulombic interactions remain virtually unchanged throughout the trajectory, whereas the Lennard-Jones interactions show a steady decrease in energy due to the reduction in unfavorable solvent-polymer interactions. This trend was also observed in the coiled structures. Because solvophobic interactions, which appear as Lennard-Jones interactions, are the primary driving force for helix formation, it is likely that changes in polymer side group polarity or solvent polarity could lead to entirely different polymer conformations being the preferred minimum energy structure. Another way of monitoring the reduction of unfavorable solvent-polymer interactions is to examine the average solvent accessible surface area (SASA) of the oligomer using the calculation procedure developed by Eisenhaber et al.54 The SASA is equivalent to the molecular surface drawn by a probe (or solvent) sphere, when the probe is rolled over the van der Waals surface of the molecule.55-60 This corresponds to a van der Waals surface in which the atomic radii are increased by the radius of the probe. The probe used in this study was a water molecule, assumed to be a sphere of radius 0.14 nm. The total SASA is the sum of the hydrophobic SASA and hydrophilic SASA. An atom is referred to as hydrophobic if it has an absolute charge value equal to or less than 0.2 e, (|e| e 0.2). Figure 9 shows the evolution of the total SASA over time for the all-atom helix, coiled, and highest (L11 and L12) and lowest (L5 and L16) energy extended structures of the m-PPE oligomer. It can be seen that the total SASA for the helix structure remains relatively constant over the course of the trajectory, whereas the coiled and extended structures show a gradual decrease in
Figure 8. Evolution of the components of potential energy of the lowest energy (L5) all-atom m-PPE extended structure.
7554 J. Phys. Chem. B, Vol. 109, No. 15, 2005
Adisa and Bruce
TABLE 4: Solvent Accessible Surface Areas (SASA) for the Helix, Extended and Coiled Conformations of the m-PPE Oligomer All-Atom SASA (nm2) (hydrophobic atoms) structure a
coiled extendedb helixa
SASA (nm2) (hydrophilic atoms)
(nm2)Total SASA
t0
t4-5c
t0
t4-5c
t0
t4-5c
33.82 33.64 19.58
24.47 ( 3.57 26.30 ( 3.95 20.19 ( 0.58
3.92 3.86 2.95
4.12 ( 0.32 4.15 ( 0.33 3.83 ( 0.27
37.74 37.51 22.53
28.59 ( 3.74 30.45 ( 4.15 24.14 ( 0.60
United-Atom 2
SASA (nm ) (hydrophobic atoms) structure coileda extendeda helixa a
SASA (nm2) (hydrophilic atoms)
(nm2)Total SASA
t0
t4-5c
t0
t4-5c
t0
t4-5c
30.12 38.48 16.31
21.91 ( 2.29 21.71 ( 3.17 16.94 ( 0.65
8.60 8.72 5.87
7.73 ( 0.43 7.59 ( 0.46 7.37 ( 0.27
38.72 47.20 22.18
29.64 ( 2.52 29.30 ( 3.54 24.31 ( 0.68
Standard deviation obtained from 1 simulation. b Standard deviation average for 16 simulations. c Averaged over the last 1 ns of trajectory.
total SASA over the first 2.5 ns and 4 ns of trajectory, respectively, before reaching a steady value. The hydrophobic, hydrophilic and total SASA of the m-PPE oligomer in the helix, coiled and extended conformations at the beginning and end (average of last 1 ns) of the simulation are listed in Table 4. It can be seen that the coiled and extended structures experience significant reductions in total SASA, whereas the helix shows a slight increase. The large reduction in the total SASA for both the coiled and extended structures can be attributed to these structures becoming more compact during the course of the MD simulation, whereas the slight increase in the total SASA for the helix can be attributed to a slight conformational change from the complete overlap of the phenyl rings of adjacent turns to a staggered structure in which there is partial overlap. This deviation causes a slight increase in the surface area accessible to the solvent and the internal volume of the helix, so as to better accommodate the bulky hydrophobic alkyl side groups. Figure 10 shows plots of the Rg and potential energy versus the total SASA obtained from the 16 all-atom extended structure simulations. The data indicate that the Rg and potential energy increase with the total SASA. This is expected because increases
in SASA require that more of the hydrophobic polymer be exposed to unfavorable interactions with the polar solvent. We also observed that the initial structure of the extended oligomer had more effect on final conformation than the initial velocity vectors of the atoms; hence, simulations of polymer folding should examine a variety of unrelated starting conformations so as to provide a better sampling of phase space. The effect of temperature on helix stability was studied via analysis of all-atom m-PPE oligomer trajectories from several canonical simulations at different temperatures. Though several techniques have been used to monitor the extent of helix formation, for example Elmer and Pande19 monitored the dihedral angles, we chose instead a more straightforward approach, which involves a comparison of the radial distribution functions between the oxygen atoms (O34) of the m-PPE oligomer in the helix conformation. The radial distribution function (rdf), gAB(r) between particles of type A and B is defined as
Figure 9. Evolution of the total solvent accessible area (SASA) for the all-atom helix, coiled, and highest (L11 and L12) and lowest (L5 and L16) energy extended structures of the m-PPE oligomer.
Figure 10. Plots of Rg (a) and potential energy (b) versus the total SASA obtained from the 16 all-atom m-PPE extended structure simulations. Data averaged over the last 1 ns of MD trajectory.
gAB(r) )
〈FB(r)〉 〈FB〉local
)
1
1
〈FB〉local NA
NA NB δ(r
∑ ∑ i∈A j∈B
ij
- r)
4πr2
(4)
where 〈FB(r)〉 is the particle density of type B at a distance r
MD Simulations of m-Poly(phenyleneethynylene)s
Figure 11. Computed radial distribution functions g(r) between the oxygen atoms (O34) of the m-PPE oligomer in the helix conformation at different temperatures and corresponding distances between oxygen atoms (inset).
around particles A, and 〈FB〉local is particle density of type B averaged over all spheres around particles A with radius rmax. Figure 11 shows the radial distribution functions between the oxygen atoms (O34) of the helix structure from MD simulations carried out at different temperatures. It can be seen that the shape of the rdfs between the O34 atoms of the helix at the various temperatures are quite similar. The peaks a, b and c correspond to the distances shown in the inset in Figure 11. Peak a, centered at 0.45 nm, represents the average distance between oxygen atoms of adjacent helix turns, b, centered at 0.79 nm, corresponds to the average distance between the oxygen atoms within the same helix turn, and c, at 0.92 nm, represents the average distance between oxygen atoms that are two helix turns apart. Furthermore, the m-PPE oligomer retains its helix conformation at the end of all simulations done at the temperatures between 300 and 523 K, thus indicating that the helix is stable within the temperature range investigated in this study. Conclusions This molecular modeling study has shown that a helical structure is the preferred conformation of an amine-functionalized m-PPE oligomer in aqueous environment at infinite dilution conditions. In addition, we observed that Lennard-Jones interactions are the main driving forces for helix formation and that these solvophobic forces are sufficient to keep the oligomer in a helical conformation at temperatures as high as 523 K. We also observed that both the united-atom and all-atom models predicted similar behavior for the m-PPE oligomer. Thus, this study clearly shows that MD simulations provide an effective means for predicting the preferred conformations of phenyleneethynylene oligomers in aqueous solutions. Acknowledgment. We acknowledge the financial support of the National Science Foundation (NSF) (CAREER 9985022) and the ERC Program of NSF under Award Number EEC9731680. Supporting Information Available: OPLS and quantum charges for amine-functionalized m-PPEs (Table 1S) and results of m-PPE extended structure simulations (Table 2S). This material is available free of charge via the Internet at http:// pubs.acs.org. References and Notes (1) Anfinsen, C. B. Science 1973, 181, 223-230. (2) Nelson, J. C.; Moore, J. S.; Wolynes, P. G.; Saven, J. G. Science 1997, 277, 1793-1796.
J. Phys. Chem. B, Vol. 109, No. 15, 2005 7555 (3) Dill, K. A. Biochemistry 1990, 29, 7133-7155. (4) Friedman, R. A.; Honig, B. Biophys. J. 1995, 69, 1528-1535. (5) Kauzmann, W. AdV. Protein. Chem. 1959, 14, 1-63. (6) Prince, R. B.; Saven, J. G.; Wolynes, P. G.; Moore, J. S. J. Am. Chem. Soc. 1999, 121, 3114-3121. (7) Arnt, L.; Tew, G. N. J. Am. Chem. Soc. 2002, 124, 7664-7665. (8) Arnt, L.; Tew, G. N. Macromolecules 2004, 37, 1283-1288. (9) Tan, C. Y.; Pinto, M. R.; Kose M. E.; Ghiviriga, I.; Schanze, K. S. AdV. Mater. 2004, 16, 1208-1212. (10) Inouye, M.; Waki, M.; Abe, H. J. Am. Chem. Soc. 2004, 126, 20222027. (11) Orita, A.; Nakano, T.; An, D. L.; Tanikawa, K.; Wakamatsu, K.; Otera, J. J. Am. Chem. Soc. 2004, 126, 10389-10396. (12) Blatchly, R. A.; Tew, G. N. J. Org. Chem. 2003, 68, 8780-8785. (13) Jones, T. V.; Blatchly, R. A.; Tew, G. N. Org. Lett. 2003, 5, 32973299. (14) Shotwell, S.; Windschief, P. M.; Smith, M. D.; Bunz, U. H. F. Org. Lett. 2004, 6, 4151-4154. (15) Grubbs, R. H. and Kratz, D. Chem. Ber. Recl. 1993, 126, 149157. (16) Leach, A. R. Molecular Modelling: Principles and Applications, 2nd ed.; Prentice Hall: Essex, U.K., 2001. (17) Binder, K. Monte Carlo and Molecular Dynamics Simulations in Polymer Science; Oxford University Press: New York, 1995. (18) Jorgensen, W. L.; Tirado-Rives, J. J. Am. Chem. Soc. 1988, 110, 1657-1666. (19) Elmer, S.; Pande, S. V. J. Phys. Chem. B 2001, 105, 482-485. (20) Kal, L.; Skeel, R.; Bhandarkar, M.; Brunner, R.; Gursoy, A.; Krawetz, N.; Phillips, J.; Shinozaki, A.; Varadarajan, K.; Schulten, K. J. Comput. Phys. 1999, 151, 283-312. (21) Brooks, B. R.; Bruccoleri, R. E.; Olafson, B. D.; States, D. J.; Swaminathan, S.; Karplus, M. J. Comput. Chem. 1983, 4, 187-217. (22) Sen, S. J. Phys. Chem. B. 2002, 106, 11343-11350. (23) MacKerell, A. D., Jr.; Wiorkiewicz-Kuczera, J.; Karplus, M. J. Am. Chem. Soc. 1995, 117, 11946-11975. (24) Jorgensen, W. L.; Chandrasekhar, J.; Madura, J. D.; Impey, R. W.; Klein, M. L. J. Chem. Phys. 1983, 79, 926-935. (25) Lee, O.; Saven, G. S. J. Phys. Chem. B 2004, 108, 11988-11994. (26) Berendsen, H. J. C.; van der Spoel, D.; van Drunen, R. Comput. Phys. Comm. 1995, 91, 43-56. (27) Berendsen, H. J. C.; van der Spoel, D.; van Drunen, R.; van Buuren, A. R.; Apol, E.; Meulenhoff, P. J.; Tieleman, D. P.; Sijbers, A. L. T. M.; Feenstra, K. A.; Lindahl, E.; Hess, B. Gromacs User Manual Version 3.1.1; Nijenborgh 4, 9747 AG Groningen: The Netherlands. Internet: www.gromacs.org, 2002. (28) Lindahl, E.; Hess, B.; van der Spoel, D. J. Mol. Mod. 2001, 7, 306-317. (29) Jorgensen, W. L.; Maxwell, D. S.; Tirado-Rives, J. J. Am. Chem. Soc. 1996, 118, 11225-11236. (30) Rizzo, R. C.; Jorgensen, W. L. J. Am. Chem. Soc. 1999, 121, 48274836. (31) Price, M. L. P.; Ostrovsky, D.; Jorgensen, W. L. J. Comput. Chem. 2001, 22, 1340-1352. (32) Ryckaert, J. P.; Bellemans, A. Faraday Discuss. Chem. Soc. 1978, 66, 95-99. (33) Jaguar 4.0; Schrodinger, Inc.: Portland, Oregon, 2000. (34) Becke, A. D. J. Chem. Phys. 1993, 98, 5648-5652. (35) Lee, C.; Yang, W.; Parr, R. G. Phys. ReV. B 1988, 37, 785-789. (36) Pople, J. A.; Krishnan, R.; Binkley, J. S.; Seeger, R. J. Chem. Phys. 1980, 72, 650-654. (37) Materials Studio version 2.1; Accelrys Inc., 2001. (38) Berendsen, H. J. C.; Postma, J. P. M.; von Gunsteren, W. F.; Hermans, J. In Intermolecular Forces; Pullman, B., Eds.; Reidel: DordrechtHolland, 1981. (39) Eisenberg, D.; McLachlan, A. D. Nature 1986, 319, 199-203. (40) Nocedal, J. Math. Comput. 1980, 35, 773-782. (41) Hockney, R. W. Methods Phys. 1970, 9, 136-211. (42) Potter, D. Computational Physics; Wiley and Sons: London, 1973. (43) Berendsen, H. J. C.; Postma, J. P. M.; DiNola, A.; Haak, J. R. J. Chem. Phys. 1984, 81, 3684-3690. (44) Hess, B.; Bekker, H.; Berendsen, H. J. C.; Fraaije, J. G. E. M. J. Comput. Chem. 1997, 18, 1463-1472. (45) Miyamoto, S.; Kollman, P. A. J. Comput. Chem. 1992, 13, 952962. (46) Weiner, S. J.; Kollman, P. A.; Case, D. A.; Singh, U. C.; Ghio, C.; Alagona, G.; Profeta, S., Jr.; Weiner, P. J. Am. Chem. Soc. 1984, 106, 765784. (47) Weiner, S. J.; Kollman, P. A.; Nguyen, D. T.; Case, D. A. J. Comput. Chem. 1986, 7, 230-252. (48) Schmieder, K.; Levitus, M.; Dang, H.; Garcia-Garibay, M. A. J. Phys. Chem. A. 2002. 106, 1551-1556. (49) Sipachev, V. A.; Khaikin, L. S.; Grikina, O. E.; Nikitin, V. S.; Trætteberg, M. J. Mol. Struct. 2000. 523, 1-22.
7556 J. Phys. Chem. B, Vol. 109, No. 15, 2005 (50) Liberles, A.; Matlosz, B. J. Org. Chem. 1971. 36, 2710-2713. (51) Okuyama, K.; Hasegawa, T.; Ito, M.; Mikami, N. J. Phys. Chem. 1984. 88, 1711-1716. (52) Humphrey, W.; Dalke, A.; Schulten, K. J. Mol. Graph. 1996, 14, 33-38. (53) Sugita, Y.; Okamoto, Y. Chem. Phys. Lett. 1999, 314, 141-151. (54) Eisenhaber, F.; Lijnzaad, P.; Argos, P.; Sander, C.; Scharf, M. J. Comput. Chem. 1995, 16, 273-284. (55) Richards, F. M. Annu. ReV. Biophys. Bioeng. 1977, 6, 151-176.
Adisa and Bruce (56) Richmond, T. J.; Richards, F. M. J. Mol. Biol. 1978.119, 537555. (57) Richmond, T. J. J. Mol. Biol. 1984.178, 63-89. (58) Kaliannan, P.; Gromiha, M. M.; Elanthiraiyan, M. Int. J. Biol. Macromol. 2001, 28, 135-141. (59) Hermann, R. B. J. Phys. Chem. 1972, 76, 2754-2759. (60) Deanda, F.; Pearlman, P. S. J. Mol. Graph. Model. 2002, 20, 415425.
