Conformational Transitions of Melittin between Aqueous and Lipid

Aug 16, 2018 - Stephen J. Fox† , Rajamani Lakshminarayanan∥⊥ , Roger W. ... Bioinformatics Institute (A*Star) , 30 Biopolis Street, #07-01 Matri...
0 downloads 0 Views 1MB Size
Subscriber access provided by Karolinska Institutet, University Library

B: Biophysics; Physical Chemistry of Biological Systems and Biomolecules

Conformational Transitions of Melittin between Aqueous and Lipid Phases: Comparison of Simulations with Experiments Stephen J. Fox, Rajamani Lakshminarayanan, Roger W. Beuerman, Jianguo Li, and Chandra S. Verma J. Phys. Chem. B, Just Accepted Manuscript • DOI: 10.1021/acs.jpcb.8b06781 • Publication Date (Web): 16 Aug 2018 Downloaded from http://pubs.acs.org on August 19, 2018

Just Accepted “Just Accepted” manuscripts have been peer-reviewed and accepted for publication. They are posted online prior to technical editing, formatting for publication and author proofing. The American Chemical Society provides “Just Accepted” as a service to the research community to expedite the dissemination of scientific material as soon as possible after acceptance. “Just Accepted” manuscripts appear in full in PDF format accompanied by an HTML abstract. “Just Accepted” manuscripts have been fully peer reviewed, but should not be considered the official version of record. They are citable by the Digital Object Identifier (DOI®). “Just Accepted” is an optional service offered to authors. Therefore, the “Just Accepted” Web site may not include all articles that will be published in the journal. After a manuscript is technically edited and formatted, it will be removed from the “Just Accepted” Web site and published as an ASAP article. Note that technical editing may introduce minor changes to the manuscript text and/or graphics which could affect content, and all legal disclaimers and ethical guidelines that apply to the journal pertain. ACS cannot be held responsible for errors or consequences arising from the use of information contained in these “Just Accepted” manuscripts.

is published by the American Chemical Society. 1155 Sixteenth Street N.W., Washington, DC 20036 Published by American Chemical Society. Copyright © American Chemical Society. However, no copyright claim is made to original U.S. Government works, or works produced by employees of any Commonwealth realm Crown government in the course of their duties.

Page 1 of 22 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

The Journal of Physical Chemistry

Title: Conformational Transitions of Melittin between Aqueous and Lipid Phases: Comparison of Simulations with Experiments Authors: Stephen J Fox , Rajamani Lakshminarayanan , Roger W Beuerman , Jianguo Li , Chandra S Verma #1

#*145

45

5

*123

Bioinformatics Institute (A*Star), 30 Biopolis Street, #07-01 Matrix, Singapore 138671 Department of Biological Sciences, National University of Singapore, 16 Science Drive 4, Singapore 117558 School of Biological Sciences, Nanyang Technological University, 60 Nanyang Drive, 4 Singapore 637551 Eye ACP, Duke-NUS Graduate Medical School, Singapore 169857 Anti-Infectives Research Group, Singapore Eye Research Institute, Singapore 168751 1

2

3

4 5

• •

# equal contributions * corresponding authors email: [email protected]; [email protected]

Abstract Peptides are promising drug candidates with advantageous therapeutic properties. However, their inherent flexibility makes the development of structure-activity relationships difficult. Molecular dynamics simulations have been widely used to study peptide conformations, but are limited by force field parameters. We explore the ability of nine combinations of commonly used protein, lipid and water force fields models (ff99/tip3p, ff14SB/tip3p, c22/tip3p, c22/tips3p, c36/tip3p, c36/tips3p, c36m/tip3p, c36m/tips3p, g53a6/spc) in capturing the conformational dynamics of the antimicrobial peptide melittin between the aqueous and model membrane environments. Circular Dichroism experiments of melittin displayed a structural transition from a random coil in aqueous solution to a helix in the presence of a model membrane. Of the protein/lipid/water models that we examined, c22 with the tips3p water model correctly recapitulated the experimentally observed disordered conformations in aqueous solution and helical conformations in the presence of the model membrane, followed by c36/tips3p. Hydration analysis revealed that the tips3p water model leads to stronger peptide-water interactions that in turn better describe solvation and its effects on conformational distributions in aqueous and membrane environments. The results of this study reveal the secondary structure preferences of various force fields and emphasize the role of hydration and micro-environment in modulating peptide conformations.

1 ACS Paragon Plus Environment

The Journal of Physical Chemistry 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Introduction Antimicrobial peptides (AMPs) are ubiquitous1,2. They are ~ 10-50 amino acids in length and are highly cationic, with charges usually in the range of +2 to +73. With resistance to small molecule antibiotics assuming alarming proportions world-wide, AMPs offer great potential as the next generation of antibiotics. They target bacterial membranes, which are harder for bacteria to remodel as a resistance mechanism, thus making them susceptible to this class of molecules4. Current AMPs have largely been discovered by chance or by mining databases for similarities in sequences and patterns such as cationicity. Success has been limited stemming from a poor understanding of the mode of action, and the lack of rigorous design principles. Biophysical experiments can provide some mechanistic information regarding the action of AMPs, but are limited in their ability to describe atomistic details. Atomistic molecular dynamics (MD) simulations have been used to great success in elucidating details of structural changes of molecules in space and time at different scales of resolution. However, their reliability strongly depends on the accuracy of the underlying force field parameters. This is particularly important since different force fields vary in their abilities to stabilize specific conformational states of peptides and proteins. These differences arise because of variations in the parameterization methods and solvent models used in the development of different force fields, and the differences in attempts to balance the non-bonded and bonded potentials. Despite these variations, current force fields are very good at reproducing crystal structures of globular proteins. However, this partly results from the short simulation times that have traditionally been accessible and partly because the force fields are optimized against the crystal structures. With increasing computational power, the time scales of simulations increases, and errors in the force fields becomes more noticeable. Hence, rigorous comparison of available force fields with experimentally available data is constantly required along with iterative refinement of the associated parameter sets5–9. There are several experimental methods that provide a variety of information regarding protein conformations and are used for verifying and improving the accuracy of force fields. Beauchamp et al9 used NMR chemical shifts to compare 11 different force fields. Their results showed that the best force fields had errors comparable to errors in experiments and suggest that further improvements in force fields will require improvements in NMR measurements 2 ACS Paragon Plus Environment

Page 2 of 22

Page 3 of 22 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

The Journal of Physical Chemistry

for proper benchmarking. Similarly, Lange et al8 examined the accuracy of 10 different force fields as well as the effects of the treatment of electrostatic interactions using experimental NMR chemical shifts for the proteins Ubiquitin and Protein G. Their results were highly dependent on the choice of the force field and the treatment of electrostatic interactions. They found that FF99SB in conjunction with PME gave the best results. However, they reported that longer time scales showed transitions to non-native conformational states, and/or persistence in high free energy states, skewing the obtained population frequencies. In 2015 Martin-Garcia et al7 reported progress when they found that long timescale simulations of these same two proteins with more modern force fields did reasonably well, although the balance of secondary structural elements (i.e helix and coil) was still unsatisfactory. Robertson et al10 evaluated several popular force fields in reproducing experimental NMR chemical shifts and experimental protein-ligand binding free energies of the flavodoxin/flavin mononucleotide system. Their results showed that all the force fields they tested performed well, although OPLS-AA/M was the best at computing the binding free energies using FEP calculations in close agreement with experimental data. Another parameter that needs to be explored is that of adequate sampling in MD simulations, especially when studying flexible proteins or small peptides. Marzinek et al11 compared the difference in the ability of various force fields to reproduce experimental structures, and additionally also explored the efficiency of several sampling methods using the small and highly flexible fusion peptide of a flavivirus envelope protein. Their results showed that temperature-based replica exchange molecular dynamics (TREMD) in conjunction with the amber ff99SB-ildn-Q force field produced the best results in reproducing the crystal structure. However, given that force fields are parametrised for globular proteins, the secondary structures of classical globular proteins render them less sensitive to errors inherent in force fields. Small peptides, and intrinsically disordered peptides, are however much more sensitive to the choice of the force field. There have been a number of attempts to improve the description of unfolded proteins. Yoo et al made several corrections to the non-bonded terms in amber and charmm force fields12. With their additions, they found that unfolded states were much better described, and with a smoothened free energy landscape, folding occurred much faster and the folded state was the free energy minimum. Additionally they also performed a re-parametrisation of pair-wise interactions to reduce artificial aggregation of water-soluble proteins13. The recent CHARMM36m force field has modifications to the CMAP potentials to enhance the conformational sampling of intrinsically disordered peptides14. We explore here the ability of force fields to capture the conformational dynamics of AMPs which are peptides known to largely adopt unstructured conformations in water and secondary/tertiary structures in the presence of 3 ACS Paragon Plus Environment

The Journal of Physical Chemistry 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

lipids (models of cell membranes). While studies have explored the abilities of force fields to recapitulate the known experimental structural information of AMPs, there are no reports on the conformational transitions they undergo between the aqueous and the lipid phases. In this study, we choose melittin as a model AMP peptide and examine the extent to which commonly used force fields can capture the conformational differences between the aqueous and membrane environments, using CD spectroscopic data as benchmark. Melittin is the principal toxic component in bee venom and displays antimicrobial properties by acting as a cell lytic membrane factor15,16. It is soluble in water and readily integrates into, and disrupts membranes. It is a small cationic peptide of 26 amino acids (sequence GIGAVLKVLTTGLPALISWIKRKRQQ). The experimental structure is mainly amphiphilic, with an asymmetry in polarity such that 13 of the first 20 amino acids are hydrophobic while 4 of the 6 C-terminal amino acids are charged (net charge is +5). Melittin is well studied, both experimentally and computationally, with the 3D structure being well characterised using x-ray crystallography and NMR17–19. Melittin exists as a tetramer at high concentrations or in the crystalline state. In this form, the hydrophobic residues line the inside of the cluster, whereas polar and charged groups are on the outside. The crystal structure of melittin (PDB id 2mlt, resolution 2.0 Å) shows a tetrameric assembly, with each monomer adopting a helical conformation, perturbed by a kink at proline14, dividing the overall structure into an N-terminal and a C-terminal helix. We first carry out CD experiments to study the differences in conformations adopted by melittin between the aqueous and membrane environments. Using the CD data as benchmark, we evaluate the ability of amber ff99SB, ff14SB, Gromos 53a6, charmm22, charmm36 and charmm36m force fields to recapitulate this data. It is known that melittin non-specifically targets various membranes, such as the mammalian membrane, the Gram positive bacterial membrane and the Gram negative bacterial membrane. As bacterial membranes are much more complicated, with various lipid compositions, we choose the simplest model of a membrane, that formed by POPC; this model has been used successfully in numerous studies20–22.

Methods Modelling and simulations: A structural model of melittin in a linear conformation was created using the tleap module in the Amber14 package23, with a positively charged N-terminus and a negatively charged C-terminus. For simulations in solvent, the linear peptide was placed in the centre of a rectangular box and solvated with TIP3P24 water molecules for the amber force fields, SPC25 for GROMOS53a6, and both TIP3P and TIPs3P for charmm force fields. The TIPs3P26 water model contains 4 ACS Paragon Plus Environment

Page 4 of 22

Page 5 of 22 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

The Journal of Physical Chemistry

additional VdW parametres on the hydrogens, and has been specifically parametrised for charmm force fields. Counterions were added to neutralize the system. For the simulations in membrane environment, a membrane patch consisting of 128 POPC lipids was used to represent the mammalian membrane. To construct the membrane, first we manually put the required number of lipid molecules on grids, resulting in a preassembled bilayer in the xy dimension. Then the preassembled bilayer was solvated with water and neutralized by adding counterions. Subsequently MD simulations were performed until convergence was seen for structural parameters such as area per lipid. Initially an unstructured conformation of melittin from the simulations in water was placed near the membrane surface and solvated the same way as above. Then the system was first subjected to 500 steps of energy minimization using the steep descent algorithm, followed by a 20 ps NVT simulation and subsequently the trajectory was propagated through production runs. Hamiltonian replica exchange was used as described in references27–31. We chose HREMD rather than TREMD because: (1) HREMD can efficiently sample the conformational space of peptides, but requires much fewer replicas than TREMD, making it computationally cheaper; (2) in the case of membrane simulations, TREMD fails because the membrane will likely collapse at high temperatures. In HREMD simulations, a number of replicas differing in their Hamiltonians were run in parallel. The interactions involving the peptide atoms (e.g., peptidepeptide and peptide-water interactions) in the Hamiltonians of the replicas were rescaled from 0% for the first replica (no rescaling) to 50% (the last replica). The first replica corresponds to the MD simulation of a normal Hamiltonian, while the last replica uses a virtual Hamiltonian with a much smoother free energy surface, allowing sampling of a much larger conformational space of the peptide. At certain time intervals, a Monte Carlo move was carried out to exchange the configurations of adjacent replicas, with the exchange acceptance rate determined by the Metropolis criterion32. If sufficient number of exchanges were accepted, the conformations sampled by the last replica can transfer to the first replica, resulting in significant enhancement in phase space sampling. In each HREMD simulation carried out in this study, 12 replicas were chosen, and the exchange was tried every 1 ps, which resulted in an acceptance ratio of 1020%, depending on the force field. In each HREMD simulation, the short-range electrostatic interactions and van der Waals interactions were truncated at 1.0 nm, while the long-range electrostatic interactions were calculated using PME33. The production simulation was carried out for 200 ns for each replica in the NPT ensemble, totalling 2.4 μs for each system. Temperature and pressure were maintained at 300K and 1 bar using the Nose-Hoover thermostat and ParrinelloRahman barostat, respectively. The secondary structures were determined by DSSP34. All simulations were carried out using the GROMACS4.6 package patched with PLUMED35. 5 ACS Paragon Plus Environment

The Journal of Physical Chemistry 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

The analysis of secondary structures was performed using the DSSP method. DSSP defines the protein secondary structure in terms of 3-turn helix, 5-turn helix, 5-turn helix, turn (hydrogen-bonded), bend (non-hydrogen bonded turn), residue isolated in beta-bridge, extended strand in parallel or anti-paralell betasheet, coil (none of the other classifications)34.

Circular Dichroism experiments: Far UV-CD spectra of melittin in PBS buffer (representing the aqueous environment) and in the presence of large unilamellar vesicles (LUVs) made up of a mixture of POPC+Cholesterol (in 3:1 ratio) in PBS to mimic the mammalian membrane environment were carried out on a Chirascan CD Spectropolarimeter using a 0.1 cm path length quartz cuvette at 20 °C. The Melittin:LUV molar ratio was maintained at 1:20. Spectra were recorded from 190 nm to 260 nm in 1.0 nm intervals and the final spectrum was taken to be the average of two scans. The final peptide concentration was maintained at 50 µg/ml.

Results and Discussion Experimental CD

Figure 1. CD of melittin in water and on the surface of a model POPC/Cholesterol membrane.

6 ACS Paragon Plus Environment

Page 6 of 22

Page 7 of 22 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

The Journal of Physical Chemistry

The CD spectrum was obtained for melittin in solution and in the presence of LUV (POPC/cholesterol=3:1). The data presented in Figure 1 shows that melittin largely adopts a random coil conformation in water (black) and a predominant α-helical conformation in the membrane-mimetic environment, as revealed by the two intense negative peaks around 208 and 222 nm, (blue). Estimation of the α-helical content by the method of Scholtz et al., indicated that in PBS it was 5.6%, and increased to 16.2% in the presence of the membrane-mimetic environment36. The data suggests a conformational transition to increasing helical conformations when melittin transfers from the aqueous phase to the membrane phase.

Secondary structure analysis Tables 1 and 2 present the secondary structures adopted by melittin during the REMD simulations in solution and on the membrane respectively, for the different force fields tested. The values are the percentages of residues found in a particular secondary structure, averaged over the final 100 ns of the REMD trajectories. As can be seen from Table 1, the secondary structures sampled in aqueous solvent favour a coil or bend conformation (both representative of unstructured peptides) in all the various force fields. Some force fields do populate beta sheet (ff99SB and g53a6) or alpha helical (ff14SB and c22/tip3p) conformations. Both C22 and c36(M) in combination with the charmm specific tips3p water model appear to perform the best in describing a random coil structure for melittin solvated in water, in good agreement with the CD data. Table 1. Average secondary structures sampled by melittin in water during the final 100ns of the REMD simulations. Standard deviations from block averaging (5 blocks from 100 ns) for helicity shown in brackets. ff99SB ff14SB g53a6 c22/tips3p c22/tip3p c36/tips3p c36/tip3p c36M/tips3p c36M/tip3p

Coil 46% 45% 41% 69% 37% 56% 68% 63% 53%

B-Sheet 13% 1% 21% 1% 0% 0% 3% 2% 6%

B-Bridge 3% 2% 2% 1% 1% 0% 2% 2% 4%

Bend 15% 15% 21% 22% 23% 37% 23% 26% 30%

Turn 16% 22% 14% 6% 19% 6% 4% 5% 6%

Helix 7%(4%) 16%(3%) 0%(0%) 1%(1%) 19%(2%) 0%(1%) 0%(0%) 0%(1%) 1%(1%)

Table 2. Average secondary structures sampled by melittin on the POPC membrane during the final 100 ns of the REMD simulations. Standard 7 ACS Paragon Plus Environment

The Journal of Physical Chemistry 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

deviations from block averaging (5 blocks from 100 ns) for helicity shown in brackets.

ff99SB ff14SB g53a6 c22/tips3p c22/tip3p c36/tips3p c36/tip3p c36M/tips3p c36M/tip3p

Coil 44% 47% 65% 35% 44% 61% 65% 70% 63%

B-Sheet 13% 1% 1% 0% 1% 1% 1% 1% 2%

B-Bridge 3% 3% 1% 0% 0% 1% 2% 2% 3%

Bend 21% 25% 20% 15% 15% 20% 19% 20% 21%

Turn 16% 19% 13% 9% 14% 10% 5% 5% 10%

Helix 4% (1%)) 5% (3%) 0% (0%) 41% (3%) 27% (6%) 7% (1%) 8% (2%) 2% (1%) 0% (1%)

The average secondary structure for melittin on the POPC membrane is presented in Table 2. The final column shows the average helicity seen during the simulations. As can be seen, the c22 force field predicts the highest helical populations, with 41% and 27% helicity for tips3p and tip3p water models, respectively, although the helices are not identical and span a conformational ensemble as is evident from Figure 3. The c36 force field, also predicts a certain amount of helicity (7% and 8% for tips3p and tip3p water models respectively). The percentage helicity sampled in the other force fields are all within 5%. In the HREMD simulations, the peptide encounters lower free energy barriers between metastable states due to rescaling of the interactions, resulting in much more efficient sampling of the phase space than in conventional MD simulations. Block averaging was utilised to check the convergence of secondary structures during the simulations. The final 100 ns of each simulation (that has been analysed) was further split into five 20 ns blocks, and the average helicity over each block was calculated. The standard deviations over the blocks are included in Tables 1 and 2. A change in the conformation adopted by a single amino acid being classified as ‘helix’ represents a 4% change. In addition, three simulations (c22/tips3p, c36/tips3p, and ff14SB) were extended for an additional 100 ns each, generating a total of 300 ns of production simulation for each. Figure 2 shows the averages of the block analysis from 0300 ns, with block size 20 ns. These structures are very dynamic and the amount of helicity can change by around 2-3 amino acids quite easily, changing the percentage by up to 8-12 %, resulting in high statistical fluctuations. For example, the percentage of helicity from 100-200 ns for ff14SB is ~4%. This increases to ~11% in the 200-300 ns block (a change of 2 amino acids). Likewise, c36/tips3p increases to ~13% (one additional helical amino acid). However, ff14SB and c36/tips3p predict much lower helicity than c22/tips3p does, consistent with the analysis of the 200 ns simulation data. 8 ACS Paragon Plus Environment

Page 8 of 22

Page 9 of 22 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

The Journal of Physical Chemistry

Although ff14SB does reproduce some helical structure on the membrane (511%), it is still too helical in solvent (compared to the CD data). In contrast, the charmm force fields predict reasonable structures in both aqueous and membrane phases. Although c22 is known to have a helical bias (which resulted in the refinements that were incorporated into c3637), and does indeed predict higher helicity on the membrane, it samples only 1% helix in bulk solvent with tips3p. Further analysis has been performed to elucidate the differences between the charmm force fields and water model combinations.

Figure 2 Block averaging over 300 ns HREMD simulations using c22/tips3p, c36/tips3p and ff14SB force fields.

9 ACS Paragon Plus Environment

The Journal of Physical Chemistry 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Figure 3. A. Crystal structure of melittin (PDB id: 2mlt). B. Superimposed snapshots of melittin (c22-cMAP/tips3p) on the membrane. In both, structures are shown as a cartoon depicting the backbone of the peptide. C. Histogram of c22/tips3p peptide backbone RMSD for the last 100 ns of REMD simulations. The crystal structure (PDB id: 2mlt) was used as the reference. Figure 3 shows representative helical structures of melittin: one is a linear helix, the crystal structure 2mlt.pdb, and the other is superimposed conformations sampled in the HREMD simulation of c22/tips3p force field. The various conformations sampled can also be seen from the wide range of RMSD in Figure 3.C (from 0.8 nm to 6.5 nm). The diverse helical conformations of melittin in the presence of POPC membrane suggests multiple modes of interactions in the early stages of melittin binding to POPC membranes. Overall the simulations reveal that melittin is highly flexible and very dynamic in the presence of POPC membrane.

10 ACS Paragon Plus Environment

Page 10 of 22

Page 11 of 22 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

The Journal of Physical Chemistry

The CHARMM family force fields use cMAP. cMAP is a grid correction energy map previously shown to improve agreement with experiment38. In order to check the importance of cMAP with respect to the melittin secondary structure, simulations with c22/tips3p were also performed without cMAP. The secondary structure results are presented in Table 3. As can be seen, in the simulation in solvent, several turns and a small amount of helix are sampled. In addition, a significant reduction in alpha-helicity and an increase in pi-helix is observed in the membrane simulations. These observations suggest the need for inclusion of cMAP corrections in the simulations with c22/tips3p. Table 3. Average secondary structures sampled by melittin in the final 100ns of the REMD simulations for c22/tips3p with no cMAP enabled. Coil B-Sheet B-Bridge Bend Turn Helix 5-Helix Solvent 47% 0% 2% 43% 3% 2% 3% POPC membrane 42% 1% 3% 28% 5% 5% 16%

Water interactions Since the c22 force field was the most successful in producing conformations of melittin in agreement with the CD experiments, we next decided to examine the influence on the sampling, of a few different water models.

Figure 4. Hydration number within 0.5 nm of melittin with the c22 and c36 force fields: analysis of HREMD simulation trajectories. The differences in hydration arise from: (1) differences in the interactions between melittin and the water model and (2) differences in the melittin conformations

11 ACS Paragon Plus Environment

The Journal of Physical Chemistry 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Different water models hydrate the peptide to different degrees. Waterbackbone interactions can significantly destabilise the helical conformation; indeed an updated tip3pm26 water model with an increased depth of the Leonard-Jones potential for the water hydrogens for improved protein-water interactions has been shown to yield better sampling of the conformations of intrinsically disordered peptides39. Figure 4 shows the number of water molecules within 0.5 nm of the peptide for both c22 and c36 force fields. In general, the hydration number depends on the peptide conformation, in addition to interactions with other non-solvent molecules. We observe that the helical conformations are the least solvated, favouring the formation of intra-peptide hydrogen bonds between backbone polar atoms. In contrast, the random coil is highly solvated. From figure 4, in the solvent environment c22/tip3p (helix 19%) displays a much smaller hydration number than c22/tips3p (helix 1%), c36/tip3p (helix 0%) and c36/tips3p (helix 0%). In the lipid environment, c22/tip3p (helix 27%) and c22/tips3p (helix 41%) also display smaller hydration numbers than c36/tip3p (helix 8%) and c36/tips3p (helix 7%). This is also evident in the differences in solvent-peptide interaction energies (Figure 5). There is a drop in interaction energy when the peptide is on the membrane, compared to in bulk solvent. The trend suggests that when the solvation interaction energies are weaker than a certain value, there is a propensity towards a helical structure, whereas larger interaction energies favour nonhelical states. This is in agreement with the above analysis indicating that if the hydration is too strong, it can result in the disruption of backbone interactions and disrupting the secondary structure.

12 ACS Paragon Plus Environment

Page 12 of 22

Page 13 of 22 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

The Journal of Physical Chemistry

Figure 5. Normal distribution of the solvent-peptide interaction energies for the different CHARMM22/water model combinations in solvent and on the membrane.

13 ACS Paragon Plus Environment

The Journal of Physical Chemistry 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Figure 6. 2D plots of the peptide distance from the membrane against solvation for c22 and c36 with tips3p or tip3p water models. Blue colour indicates low population and red colour indicates high population.

14 ACS Paragon Plus Environment

Page 14 of 22

Page 15 of 22 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

The Journal of Physical Chemistry

Figure 7. 2D plots of peptide solvation against the number of helical residues in the peptide for c22 and c36 with tips3p or tip3p water models. Blue colour indicates low population and red colour indicates high population.

To further understand the effect of solvation on melittin helicity in the membrane environment, we calculated the probability of the peptide helicity, solvation and distance from the membrane. From Figure 6, two things stand out: firstly, the systems with the tip3p water model have a higher propensity to remain adsorbed onto the surface of the membrane than when the tips3p water model is used, for both c22 and c36 force fields. This suggests that the tip3p water model leads to stronger peptide-membrane interactions. Secondly, melittin in tip3p is less solvated than when it is in tips3p, and this is consistent with the results of Figure 4. c22/tips3p results in the exploration of a larger volume of phase space, with the peptide spending significant amounts of time both on and near the membrane. This results in melittin adopting multiple structures in the membrane environment, which are consistent with the RMSD plots in Figure 3. When modelled using c36/tips3p, melittin does not spend much time adsorbed on to the membrane, with the highest population almost ~5 nm distant from the membrane surface. This suggests weaker peptidemembrane interactions in the case of c36. The increased hydration due to the tips3p water could arise from the inclusion of VdW’s on the hydrogen atoms. This is also observed in the interaction energies of melittin which are much 15 ACS Paragon Plus Environment

The Journal of Physical Chemistry 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 16 of 22

stronger with tips3p than with tip3p (figure 3). Proximity to the membrane can only partially explain the reduction in hydration that is seen in the simulations of c36 (figure 4). Plotting helicity against solvation (figure 7) shows a good correlation for c22 with both water models and in contrast, c36 does not show any trend with either water model, and also results in reduced helicity. These suggest that the c36 parameters disfavour helicity. Membrane interactions

For membranes, the most common properties are area per lipid and thickness, both of which can be compared to experimental values. Table 4. Membrane area per lipid (in nm2); experimental value is 0.643 nm2 11

tip3p

tips3p

C22 – membrane

0.605 ± 0.013

0.634 ± 0.014

C22 - membrane/peptide

0.562 ± 0.021

0.603 ± 0.014

Experimentally a POPC membrane has an area per lipid (apl) of 0.643 nm2 with membrane thickness of 3.91 nm at 30° . Membrane properties from the c22 simulations of the membrane, with and without the peptide, in the different solvation models are shown in table 4 (area per lipid) and Table 5 (membrane thickness). The area per lipid using tips3p is in agreement with the experimental value, while in the presence of the tip3p water model, the area per lipid is underestimated by ~5%. Moreover, in the simulations of only the membrane, tips3p produces the correct membrane thickness of 3.9nm, while tip3p results in a thicker membrane. Thus, tips3p is clearly more accurate in describing the biophysical properties of the POPC membrane. 11

Table 5. Membrane thickness (in nm); experimental value is 3.91nm 40

tip3p

tips3p

C22 - membrane

4.1 ± 0.3

3.9 ± 0.3

C22 - membrane/peptide

4.1 ± 0.4

4.1 ± 0.4

Conclusions We have examined the ability of a series of commonly used force fields and water models in predicting the conformational changes that characterize melittin 16 ACS Paragon Plus Environment

Page 17 of 22 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

The Journal of Physical Chemistry

between the aqueous and membrane environments, using replica exchange MD simulations and CD spectroscopy. CD data shows that in PBS melittin has virtually no helical content, which was observed to increase when in a membrane-mimetic environment. Among the different force fields tested, the charmm force fields with the tips3p water model were more consistent with the CD result, predicting a disordered conformation of melittin in aqueous solution. In the presence of a POPC membrane, c22 predicted the highest population of helical conformations, followed by c36. The helical conformation of CHARMM force fields largely arises from the cMAP corrections. For example, in the HREMD simulation of c22 force field, switching off the cMAP correction results in random conformations in the membrane environment. The water model was found to significantly affect the melittin conformation. The tip3p water model does not affect the peptide helicity in the membrane environment, but leads to overestimation of helicity in aqueous solution. Hydration analysis revealed that the tips3p water model leads to stronger peptide-water interactions, which results from the additional LJ contribution of hydrogen atoms in the tips3p water model. As the force field is a delicate balance of solvation and intramolecular interactions, the tips3p water model results in appropriate degrees of solvation, at least for a series of HREMD simulations of melittin in two different environments as explored here. C22 has higher helical propensity, generating structures somewhere between the highly helical x-ray and nmr structures and our CD values. The more recent version of the CHARMM force field, c36, predicts helicity slightly lower than the CD value, irrespective of the water model (tip3p and tips3p) used. In the current study, the other force fields fail to show an appropriate balance between simulations in water and in the presence of membrane. Amber99sb and gromos53a6 appear to be biased towards the beta-sheet conformation, and amber14sb overestimates helical structures in water, at least for melittin in our studies. Although most force fields are pretty accurate in simulating proteins, their accuracy in simulating small peptides is still challenging because of the sensitivity of their conformation to the microenvironment. This study emphasizes the dependence of the peptide conformation on the microenvironment and the importance of the choice of water model. It should be noted that the observations in this study are obtained using the antimicrobial peptide melittin; the conformations of other peptides with totally different structural properties (e.g., net charge, hydrophobicity, amphiphilicity etc.) may display different sensitivity to force fields. Nevertheless, the findings will hopefully stimulate further developments in the force fields.

Acknowledgements 17 ACS Paragon Plus Environment

The Journal of Physical Chemistry 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

SJF would like to acknowledge A*STAR for grants - JCO 1431AFG101 & IAF-PP (HBMS Domain) H17/01/a0/009; RL would like to acknowledge Co-operative Basic Research Grant from Singapore National Medical Research Council (grant number NMRC/CBRG/0048/2013). CV is the founder of sinopsee therapeutics, a biotech company developing molecules for therapeutic purposes; the current work has no conflict with the company. We than NSCC and A*STAR for computing resources

References (1)

Boman, H. G. Peptide Antibiotics and Their Role in Innate Immunity. Annu. Rev. Immunol. 1995, 13 (1), 61–92.

(2)

Ganz, T.; Weiss, J. Antimicrobial Peptides of Phagocytes and Epithelia. Semin Hematol 1997, 34 (4), 343–354.

(3)

Boman, H. G. Antibacterial Peptides: Basic Facts and Emerging Concepts. J. Intern. Med. 2003, 254 (3), 197–215.

(4)

Yeaman, M. R. Mechanisms of Antimicrobial Peptide Action and Resistance. Pharmacol. Rev. 2003, 55 (1), 27–55.

(5)

Sandoval-Perez, A.; Pluhackova, K.; Böckmann, R. A. Critical Comparison of Biomembrane Force Fields: Protein-Lipid Interactions at the Membrane Interface. J. Chem. Theory Comput. 2017, 13 (5), 2310–2321.

(6)

Rauscher, S.; Gapsys, V.; Gajda, M. J.; Zweckstetter, M.; De Groot, B. L.; Grubmüller, H. Structural Ensembles of Intrinsically Disordered Proteins Depend Strongly on Force Field: A Comparison to Experiment. J. Chem. Theory Comput. 2015, 11 (11), 5513– 5524.

(7)

Martin-Garcia, F.; Papaleo, E.; Gomez-Puertas, P.; Boomsma, W.; Lindorff-Larsen, K. Comparing Molecular Dynamics Force Fields in the Essential Subspace. PLoS One 2015, 10 (3).

(8)

Lange, O. F.; Van Der Spoel, D.; De Groot, B. L. Scrutinizing Molecular Mechanics Force Fields on the Submicrosecond Timescale with NMR Data. Biophys. J. 2010, 99 (2), 647–655.

(9)

Beauchamp, K. A.; Lin, Y. S.; Das, R.; Pande, V. S. Are Protein Force Fields Getting Better? A Systematic Benchmark on 524 Diverse NMR Measurements. J. Chem. Theory Comput. 2012, 8 (4), 1409–1414.

(10)

Robertson, M. J.; Tirado-Rives, J.; Jorgensen, W. L. Performance of Protein-Ligand Force Fields for the Flavodoxin-Flavin Mononucleotide System. J. Phys. Chem. Lett. 2016, 7 (15), 3032–3036.

(11)

Marzinek, J. K.; Lakshminarayanan, R.; Goh, E.; Huber, R. G.; Panzade, S.; Verma, C.; 18 ACS Paragon Plus Environment

Page 18 of 22

Page 19 of 22 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

The Journal of Physical Chemistry

Bond, P. J. Characterizing the Conformational Landscape of Flavivirus Fusion Peptides via Simulation and Experiment. Sci. Rep. 2016, 6. (12)

Yoo, J.; Aksimentiev, A. Refined Parameterization of Nonbonded Interactions Improves Conformational Sampling and Kinetics of Protein Folding Simulations. J. Phys. Chem. Lett. 2016, 7 (19), 3812–3818.

(13)

Yoo, J.; Aksimentiev, A. Improved Parameterization of Amine-Carboxylate and AminePhosphate Interactions for Molecular Dynamics Simulations Using the CHARMM and AMBER Force Fields. J. Chem. Theory Comput. 2016, 12 (1), 430–443.

(14)

Huang, J.; Rauscher, S.; Nawrocki, G.; Ran, T.; Feig, M.; De Groot, B. L.; Grubmüller, H.; MacKerell, A. D. CHARMM36m: An Improved Force Field for Folded and Intrinsically Disordered Proteins. Nat. Methods 2016, 14 (1), 71–73.

(15)

Lee, M.-T.; Sun, T.-L.; Hung, W.-C.; Huang, H. W. Process of Inducing Pores in Membranes by Melittin. Proc. Natl. Acad. Sci. 2013, 110 (35), 14243–14248.

(16)

Lee, M.-T.; Hung, W.-C.; Chen, F.-Y.; Huang, H. W. Mechanism and Kinetics of Pore Formation in Membranes by Water-Soluble Amphipathic Peptides. Proc. Natl. Acad. Sci. 2008, 105 (13), 5087–5092.

(17)

Lam, Y. H.; Wassall, S. R.; Morton, C. J.; Smith, R.; Separovic, F. Solid-State NMR Structure Determination of Melittin in a Lipid Environment. Biophys. J. 2001, 81 (5), 2752–2761.

(18)

Brown, L. R.; Braun, W.; Kumar, A.; Wüthrich, K. High Resolution Nuclear Magnetic Resonance Studies of the Conformation and Orientation of Melittin Bound to a LipidWater Interface. Biophys. J. 1982, 37 (1), 319–328.

(19)

Sharon, M.; Oren, Z.; Shai, Y.; Anglister, J. 2D-NMR and ATR-FTIR Study of the Structure of a Cell-Selective Diastereomer of Melittin and Its Orientation in Phospholipids. Biochemistry 1999, 38 (46), 15305–15316.

(20)

Chen, C. H.; Wiedman, G.; Khan, A.; Ulmschneider, M. B. Absorption and Folding of Melittin onto Lipid Bilayer Membranes via Unbiased Atomic Detail Microsecond Molecular Dynamics Simulation. Biochim. Biophys. Acta - Biomembr. 2014, 1838 (9), 2243–2249.

(21)

Irudayam, S. J.; Berkowitz, M. L. Binding and Reorientation of Melittin in a POPC Bilayer: Computer Simulations. Biochim. Biophys. Acta - Biomembr. 2012, 1818 (12), 2975–2981.

(22)

Benachir, T.; Lafleur, M. Study of Vesicle Leakage Induced by Melittin. BBA Biomembr. 1995, 1235 (2), 452–460.

(23)

Case, D. A.; Babin, V.; Berryman, J. T.; Betz, R. M.; Cai, Q.; Cerutti, D. S.; Cheatham, T. E.; Darden, T. A.; Duke, R. E.; Gohlke, H.; et al. {Amber 14} OR - University of California, San Francisco; 2014.

(24)

Jorgensen, W. L.; Chandrasekhar, J.; Madura, J. D.; Impey, R. W.; Klein, M. L. 19 ACS Paragon Plus Environment

The Journal of Physical Chemistry 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Comparison of Simple Potential Functions for Simulating Liquid Water. J. Chem. Phys. 1983, 79 (2), 926–935. (25)

Berendsen, H. J. C.; Postma, J. P. M.; Gunsteren, W. F. Van; Hermans, J. Interaction Models for Water in Relation to Protein Hydration. Intermol. Forces 1981, 331–342.

(26)

MacKerell, A. D.; Bashford, D.; Bellott, M.; Dunbrack, R. L.; Evanseck, J. D.; Field, M. J.; Fischer, S.; Gao, J.; Guo, H.; Ha, S.; et al. All-Atom Empirical Potential for Molecular Modeling and Dynamics Studies of Proteins †. J. Phys. Chem. B 1998, 102 (18), 3586– 3616.

(27)

Mu, Y. Dissociation Aided and Side Chain Sampling Enhanced Hamiltonian Replica Exchange. J. Chem. Phys. 2009, 130 (16), 164107.

(28)

Terakawa, T.; Kameda, T.; Takada, S. On Easy Implementation of a Variant of the Replica Exchange with Solute Tempering in GROMACS. J. Comput. Chem. 2011, 32 (7), 1228–1234.

(29)

Liu, P.; Kim, B.; Friesner, R. A.; Berne, B. J. Replica Exchange with Solute Tempering: A Method for Sampling Biological Systems in Explicit Water. Proc. Natl. Acad. Sci. 2005, 102 (39), 13749–13754.

(30)

Bussi, G. Hamiltonian Replica Exchange in GROMACS: A Flexible Implementation. In Molecular Physics; 2014; Vol. 112, pp 379–384.

(31)

Li, J.; Hu, Z.; Beuerman, R.; Verma, C. Molecular Environment Modulates Conformational Differences between Crystal and Solution States of Human βDefensin 2. J. Phys. Chem. B 2017, 121 (13), 2739–2747.

(32)

Yuji, S.; OKamoto, Y. Replica-Exchange Molecular Dynamics Method for Protein Folding. Chem. Phys. Lett. 314 1999 141–151 2007, 350 (November), 205–223.

(33)

Essmann, U.; Perera, L.; Berkowitz, M. L.; Darden, T.; Lee, H.; Pedersen, L. G. A Smooth Particle Mesh Ewald Method. J. Chem. Phys. 1995, 103 (19), 8577–8593.

(34)

Kabsch, W.; Sander, C. Dictionary of Protein Secondary Structure: Pattern Recognition of Hydrogen-bonded and Geometrical Features. Biopolymers 1983, 22 (12), 2577– 2637.

(35)

Hess, B.; Kutzner, C.; Van Der Spoel, D.; Lindahl, E. GRGMACS 4: Algorithms for Highly Efficient, Load-Balanced, and Scalable Molecular Simulation. J. Chem. Theory Comput. 2008, 4 (3), 435–447.

(36)

Scholtz, J. M.; Qian, H.; York, E. J.; Stewart, J. M.; Baldwin, R. L. Parameters of Helix– coil Transition Theory for Alanine-based Peptides of Varying Chain Lengths in Water. Biopolymers 1991, 31 (13), 1463–1470.

(37)

Best, R. B.; Zhu, X.; Shim, J.; Lopes, P. E. M.; Mittal, J.; Feig, M.; MacKerell, A. D. Optimization of the Additive CHARMM All-Atom Protein Force Field Targeting Improved Sampling of the Backbone φ, ψ and Side-Chain Χ1and Χ2Dihedral Angles. J. Chem. Theory Comput. 2012, 8 (9), 3257–3273. 20 ACS Paragon Plus Environment

Page 20 of 22

Page 21 of 22 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

The Journal of Physical Chemistry

(38)

Buck, M.; Bouguet-Bonnet, S.; Pastor, R. W.; MacKerell, A. D. Importance of the CMAP Correction to the CHARMM22 Protein Force Field: Dynamics of Hen Lysozyme. Biophys. J. 2006, 90 (4).

(39)

Boonstra, S.; Onck, P. R.; Van Der Giessen, E. CHARMM TIP3P Water Model Suppresses Peptide Folding by Solvating the Unfolded State. J. Phys. Chem. B 2016, 120 (15), 3692–3698.

(40)

Kučerka, N.; Nieh, M. P.; Katsaras, J. Fluid Phase Lipid Areas and Bilayer Thicknesses of Commonly Used Phosphatidylcholines as a Function of Temperature. Biochim. Biophys. Acta - Biomembr. 2011, 1808 (11), 2761–2771.

21 ACS Paragon Plus Environment

The Journal of Physical Chemistry 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

TOC graphic

22 ACS Paragon Plus Environment

Page 22 of 22