Parameterization of Palmitoylated Cysteine ... - ACS Publications

Nov 16, 2017 - Hess, B.; Lindah, E. Gromacs: High Performance Molecular. Simulations through Multi-Level Parallelism from Laptops to Super- computers...
1 downloads 16 Views 1MB Size
Subscriber access provided by READING UNIV

Article

Parameterization of Palmitoylated Cysteine, Farnesylated Cysteine, Geranylgeranylated Cysteine and Myristoylated Glycine for the Martini Force Field Yoav Atsmon-Raz, and D. Peter Tieleman J. Phys. Chem. B, Just Accepted Manuscript • DOI: 10.1021/acs.jpcb.7b10175 • Publication Date (Web): 16 Nov 2017 Downloaded from http://pubs.acs.org on November 20, 2017

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 free 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 accessible to all readers and 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.

The Journal of Physical Chemistry B 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 34

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

Parameterization of Palmitoylated Cysteine, Farnesylated Cysteine, Geranylgeranylated Cysteine and Myristoylated Glycine for the Martini Force Field Yoav Atsmon-Raz1 and D. Peter Tieleman*1 1

University of Calgary, Department of Biological Sciences, Centre for Molecular Simulation,

2500 University Drive NW, Calgary, Alberta, Canada, T2N 1N4, *Corresponding author: [email protected], Phone: (403) 220-2966

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

ABSTRACT Peripheral membrane proteins go through various post-translational modifications that covalently bind fatty acid tails to specific amino acids. These post translational modifications significantly alter the lipophilicity of the modified proteins and allow them to anchor to biological membranes. Over 1000 different proteins have been identified to date that require such membrane-protein interactions in order to carry out their biological functions including members of the Src and Ras superfamilies that play key roles in cell signaling and carcinogenesis. We have used all-atom simulations with the CHARMM36 force field to parameterize four of the most common post-translational modifications for the Martini 2.2 force field – palmitoylated cysteine, farnesylated cysteine, geranylgeranylated cysteine and myristoylated glycine. The parameters reproduce the key features of clusters of configurations of the different anchors in lipid membranes as well as the water-octanol partitioning free energies of the anchors, which are crucial for the correct reproduction of the expected biophysical behavior of peripheral membrane proteins at the membrane-water interface. Implementation in existing Martini set up tools facilitates the use of the new parameters.

2 ACS Paragon Plus Environment

Page 2 of 34

Page 3 of 34

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

INTRODUCTION Cell membranes form complicated and heterogeneous mixtures of membrane-associated proteins and lipids which interact with each other. The structural interface through which the proteins interact with the membrane divides them into two groups - integral membrane proteins and peripheral membrane proteins (PMPs). Both groups differ from one another in their function, structural properties and their biochemical interactions with biological membranes - integral membrane proteins use entire domains to immerse themselves within the bilayer while many PMPs use post-translationally lipidated amino acids to anchor themselves to a bilayer or to electrostatically interact with it via a polybasic domain while keeping their surfaces exposed to the intra/extra-cellular environment.1 In most cases, such membrane anchors (MAs) are added through the S-acetylation of a cysteine sidechain by a thioester bond that binds palmitic acid to the cysteine residue through a process known as palmitoylation or through the formation of a thioether bond with an isoprenoid lipid through a process known as prenylation, which involves the addition of a farnesyl group or a geranylgeranyl group to the CAAX motif.2 The specificity of which of these moieties is added to a protein is dependent on the residue in the X position – if X is a serine, methionine, alanine or glutamine it typically becomes farnesylated and if it is a leucine, it typically becomes geranylgeranylated. Similarly, N-termini positioned glycines, may undergo an analogousprocess known as N-myristoylation through the addition of a myristic acid to its amine group. Most importantly, each of the aforementioned MAs plays a critical role in various biological processes such as cell cycle signalling,3 nuclear envelope architecture, heterochromatin organization4–6, cell differentiation, apoptosis, 7,8 lipid trafficking, lateral domain formation in membranes and cell adhesion.9 Mutations of these particular residues have been related to the development of lung, colon and pancreatic cancers.10–19 Furthermore, they have also been found to be of potential therapeutic importance as has been previously noted for myristoylation as a possible target for antifungal drugs17, HIV-1 infection via membrane binding of the retroviral gag protein,20 and as potential targets for the treatment of cancer.21 Similarly, palmitoylation has been implicated in neurological pathologies such as Huntington’s,22,23Alzheimer’s disease 24–26, bipolar disorder 27,28 and X-linked mental retardation29,30 and defects in the CAAX motif of prenylation prone proteins can lead to tumorigenesis as well as to the development of mutations 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

in A-type lamins.4,31,32 Currently, there are over 500 S-acylated proteins that have been experimentally identified in humans with over 100 of this group undergoing prenylation and 100 other proteins that have been found to undergo myrisotylation.33 Therefore, significant effort has been put into the development of predictive lipidation algorithms and structural databases such as such as CSS-palm 4.034(http://csspalm.biocuckoo.org/), Myristoylator35 (http://web.expasy.org/myristoylator/) and PRENbase36 (http://mendel.imp.univie.ac.at/sat/PrePS/PRENbase) which are used to predict additional protein lipidation targets.2 However, structure determination of peripheral membrane proteins is challenging due to the flexible, dynamic nature of their lipid environment and the propensity towards structural disorder of the anchored regions of these proteins.21,37–41 Consequently, this leads to considerable difficulty in elucidating the underlying details behind the function and structure of the different MA insertion mechanisms which characterizes PMPs such as electrostatic interactions between a protein’s polybasic domain with the membrane’s lipid headgroups,42–4445,46 or a process that involves multiple consecutive anchor insertions into a bilayer.18 Both of these mechanisms are currently poorly understood, and requires alternative methods that will allow to investigate the dynamics and conformational changes of the PMP and its interacting environment which may be provided by the use of computational approaches such as all-atom (AA) molecular dynamics simulations (MD) as well as coarse-grained (CG) MD simulations. Currently, CG-MD simulations have reached time and length scales in which PMPs can be simulated in complex lipid mixtures while recent developments in X-ray crystallography, neutron diffraction and quantitative mass spectrometry have gained spatio-temporal resolution that has provided significant overlap between computational and experimental biophysics. CGMD simulations require parameters that accurately reproduce atomistic dynamics and need to be systemically parameterized in a self-consistent manner in order to allow a reliable comparison between different systems. Previous studies have reported parameterizations of myristoylated glycine for the anchoring mechanism of HIV-1 MA protein,47 palmitoylated and farnesylated cysteine48 for the structural 4 ACS Paragon Plus Environment

Page 4 of 34

Page 5 of 34

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

effect that NRas has on bilayer phase separation49 while geranylgeranylated cysteine was parameterized for a study on the yeast exocyst complex.50 In all of these cases, the applied parameters were taken from previously existing topologies in the Martini force field that were not explicitly designed for these post-translational modifications or alternatively, assigned from isolated molecular structures such as for farnesyl methyl ether.49 Consequently, it is challenging to compare between these anchors in terms of their partitioning preferences or reliably utilize them in tandem as is relevant for proteins that undergo multiple lipidations such as N-Ras, HRas6, SNAP25,51Fyn, Lck and Hck52 and membrane dissociation as is the case for the reversible cycle of cysteine palmitoyl53,54 and for glycine myristoyl.17 Specifically, the partitioning between an aqueous and a fatty environment (or lipophilicity) of each of these modified residues is critical to the function of the PMPs in which their lipophilicity plays a key role in their localization to target compartments within the cell and onto lipid rafts in the plasma membrane.55,56 Furthermore, their partitioning between the aqueous and fatty environments are also an important physicochemical criteria for any potential drug candidate and is one of the guidelines for Lipinski’s Rule of 5.57,58 In this study, we report a self-consistent set of parameters for the Martini force field59,60 which we have modeled from AA MD simulations that were produced using the CHARMM36 force field. 61–63 Our parameters provide a basis for future studies that will focus on lipid-protein interactions and the relationship between membrane composition and the dynamics and function of PMPs.

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

Figure 1. Coarse grained representations of A) Cysteine palmitoyl, B) Cysteine farnesyl, C) Cysteine geranylgeranyl and D) Glycine myrsitoyl in the Martini force field. The P5, Nda, C5, C3 and C1 type beads are respectively colored in red, magenta, blue, lime and cyan. E) The applied mapping scheme we used to represent in coarse-grained resolution each of the four lipidated amino acids.

METHODS All atom simulation protocol In order to obtain the necessary distributions to parameterize the bonded parameters of each of the four lipidated MAs we prepared a 10 x 10nm2 POPC bilayer with 148 lipids per leaflet with the CHARMM-GUI membrane builder.64–67 Structures of MFCIH and GMFCIH peptides were prepared with the Accelrys Discovery Studio Visualizer 4.168 and inserted into each of the four lipid tails in separated setups by removing a POPC lipid from the center of the upper leaflet and inserting the anchor within the formed cavity. The resulting structures were minimized in 50,000 steps using the steepest descent algorithm, equilibrated in 1ns NPT simulation and then simulated for 100 ns production in the NPT ensemble with the CHARMM36 force field 61–63 and the TIP3P69 water model with the GROMACS 5.0.x70–74 software package. A time step of 2 fs was used. Temperature was controlled by coupling to a Nose-Hoover thermostat75,76 set at 310K with a coupling constant of τT = 1.0ps and the pressure was coupled via a semi-isotropic 6 ACS Paragon Plus Environment

Page 6 of 34

Page 7 of 34

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

Parrinello-Rahman barostat77 with a coupling constant of τp = 5.0 ps and a reference value set at 1 bar with a compressibility value of 4.5 × 10−5 bar−1. For neighbor searching, we used the Verlet group scheme which was updated every 20th step, with short-range LJ and Coulomb interactions switched off between 1.0 and 1.2 nm and long-range electrostatics were handled with the particle mesh Ewald summation scheme with a force-switch Van-der-Waals modifier applied.78,79 The force field parameters and topologies for the peptides, the POPC lipids and the modified amino acids were obtained directly from the February 2016 GROMACS implementation of the CHARMM36 force field61–63 and are available from the Mackerell group website (http://mackerell.umaryland.edu/charmm_ff.shtml). We note that the original parameters for the lipidated amino acids were published explicitly in the CHARMM format and were imported and integrated into the GROMACS format with the cgenff_charmm2gmx.py script which is also available on the MacKerell group website. All distributions for the bond and angle parameters were respectively acquired with the gmx distance and gmx gangle GROMACS tools and were calculated from the centers of masses of the relevant atom groups to allow a direct comparison with the Martini distributions.

Coarse grained mapping schemes Our chosen mapping schemes (Figure 1) were determined in accordance with the chemical moieties which were originally described for the Martini force field60 and its extension for proteins.80 Specifically, for the side chain of the modified cysteines, we mapped the thio-ether and thio-ester bridges to a C5 bead which has been previously used to represent 1-propanethiol60 instead of an Na bead that was used in the work of de Jong et al.48 For the prenylated chains of farnesylated cysteine and geranylgeranylated cysteine, we respectively assigned four and five C3 beads while for the aliphatic chains of palmitoylated cysteine and myristoylated glycine we assigned four C1 beads. For the amino moiety of myristoylated glycine we have assigned an Nda bead instead of the P5/P4 beads that is more commonly used for amino acids in the Martini force field in order to account for the N-acetylation of the N terminal amine group while the C termini were set a neutral charge.80 For the three lipidated cysteines, both C and N termini were set with a charge of 0. 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

Coarse grained simulations protocol We have produced CG topologies and initial structures for the MFCIH and GMFCIH peptides for the Martini force field59,60 by applying the martinize script74 to our initially constructed AA structures and prepared a CG POPC bilayer with INSANE81 from which we removed a single POPC in accordance with the AA bilayer. We then minimized the system for 10000 steps with the steepest descent algorithm, followed by an NPT equilibration period of 100 ns and a production trajectory of 2 µs. The temperature was coupled to a v-rescale thermostat82 which was set at a reference value of 310K with a coupling constant of τT = 1.0ps. Pressure was coupled via a semi-isotropic Parrinello-Rahman (PR) barostat77 with a coupling constant of τp = 12.0 ps, a reference value set at 1 bar with a compressibility value of 3.0 × 10−4 bar−1 with no dispersion correction being applied, as recently proposed by de Jong et al. as optimal choices for CG-MD simulations when using the Martini force field with GROMACS 5.70–74 Structural analysis – clustering and density profiles To obtain the conformational clusters of each of our trajectories, we computed the RMSD matrix of the modified amino acids followed by clustering. We used the GROMOS clustering method with a cutoff of 0.25 nm for both CG and AA simulations.83 Further structural analysis was carried out by calculating the density profiles for each of the heavy atoms in each one of the four residues’ fatty chains. Free Energy Calculations In the Martini force field, non-bonded interactions are optimized to reproduce the free energies of partitioning between water and octanol, which is experimentally determined as a partitioning coefficient. and describes its lipophilicity .60,84 The value of the partitioning coefficient can be obtained from equation (1): ∆∆ீ

ೀೈ − ଶ.ଷ଴ଷோ் = ݈‫ܲ݃݋‬ைௐ (1)

The difference between the free energies of solvation of water and octanol is described by

∆∆‫ܩ‬ைௐ = ∆‫ܩ‬ை − ∆‫ܩ‬ௐ . Both values can be computed by following an alchemical reaction that 8 ACS Paragon Plus Environment

Page 8 of 34

Page 9 of 34

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

describes the decoupling of each solute from its environment using the coupling parameter λ to obtain the free energy of solvation in each solvent in terms of the solute’s decoupling by integrating the value of ർ

డ௎ሺఒሻ డఒ

඀ over the range of all the values of λ.85 Consequently, we may ఒ

compute the free energy of partitioning between the two solvents by finding the value of the difference between the solvation energies as described in equation (2) : ఒୀଵ

∆∆‫ܩ‬ைௐ = ∆‫ܩ‬ை − ∆‫ܩ‬ௐ = ‫׬‬ఒୀ଴ ቆർ Here ർ

డ௎ሺఒሻ డఒ

డ௎ሺఒሻ డఒ



ఒ,௢௖௧௔௡௢௟

−ർ

డ௎ሺఒሻ డఒ



ఒ,௪௔௧௘௥

ቇ ݀ߣ (2)

඀ is the time averaged first-order derivative with respect to the coupling parameter λ ఒ

of the potential energy function describing the total solvent-solute interaction as a function of the coupling parameter λ which ranges between the values of 0 to 1, with 0 expressing a completely decoupled state of the solute and 1 describing a fully interacting state. Separate simulations of each λ value were performed and during each simulation the derivatives of the system’s energy were calculated using its neighboring λ states. In order to avoid any singularities in the potentials for the interactions,86 the van der Waals (VDW) interactions were treated with a soft-core potential with GROMACS soft core parameters sc_alpha = 0.5, sc_sigma=0.3 and sc_power = 1, while all bonded interactions were interpolated linearly. This procedure combines the pathways of each solvent into a single thermodynamic cycle (Figure 2) that yields the free energy of partitioning between the water-saturated octanol and pure water.

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 2. The thermodynamic cycle to calculate the free energy of partitioning between water saturated 1-octanol and water Since all four modified residues have no electrostatic charges in their martini representations we did not calculate the electrostatic component of their CG structures. For the calculation of the free energy of partitioning, we prepared a 5x5x5 nm3 cubic box of water composed of 1700 regular Martini water beads with 10% of them randomly reassigned as Martini “anti-freeze” beads. For each of the water-saturated octanol boxes we prepared a 5x5x5 nm3 cubic box with octanol which was assembled with PACKMOL87 in order to form a randomized binary mixture with a 1-octanol:water molar fraction ratio of ~0.74:0.2688,89, resulting in 35 Martini water beads, as each CG bead represents 4 real water molecules,90 and 382 Martini octanol CG molecules. The CG topologies and parameters of the water and octanol are available from the Martini website (http://cgmartini.nl). This procedure was then repeated for each of the beads in the anchors (C1/C3/C5/Nda/P5) and the results were compared to previously reported results by Bereau and Kremer on individual beads (Figure S1).91 We then prepared 21 simulations with uniformly distributed values of λ for each water-saturated octanol box and 21 simulations for each pure water box. For each value of λ we ran a simulation of a pre-equilibrated solvent box, with the lipidated residue inserted into its center so it would fully interact with the solvent with a Berendsen pressure coupling scheme and a coupling constant of 5 ps.92 Each simulation was then energy minimized (steepest descent, 10000 steps) and equilibrated in NVT for 100 ns and in NPT for another 100 ns with the appropriate λ parameter and finalized by a production run of 1 µs per window using Langevin Dynamics.93 10 ACS Paragon Plus Environment

Page 10 of 34

Page 11 of 34

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

All eight peptide-solvent (four 1-octanol boxes and four water boxes for each modified residue) systems were then copied and modified by removing their amino backbone bead as well as the first sidechain bead for all of the cysteines to probe possible discrepancies in free energy values that may be related to the amino acids themselves. In order to allow a comparison with AA simulations, we also computed the free energies of partitioning for water saturated 1-octanol and pure water with the CHARMM3661–63 force field under GROMACS 5.0.7.70–74 For these calculations, we used the same lipidated amino acids parameters which were applied in the AA anchored peptide simulations as described in the previous section. The initial structures for the lipidated amino acids were obtained by removing all of the amino acids from the MFCIH peptide besides the cysteine residue while leaving its charged amine and carboxyl termini, while for the myristoylated glycine we removed all of the residues from the GMFCIH peptide except glycine. We note that the sequence selection was randomly chosen in advance as to avoid biasing the obtained parameters in favor of any particular PMP and were validated via BLAST to ensure that they indeed do not represent any specific case.94 However, we are aware that previous studies have demonstrated that the position and the composition of the peptide may affect the partitioning between the liquid-ordered and liquid-disordered domains in a membrane,95 but are expected to be of minor importance in comparison to the significantly increased contribution to the hydrophobic propensity by the lipidated residues.96 Both the topology and parameters for 1-octanol were obtained from the CGenFF forcefield.97 AA solvent boxes were prepared with PACKMOL87 for each of the three modified cysteines with identical sizes to those that were previously used for the coarse-grained simulations, with 140 TIP3P water molecules and 382 octanol molecules. For the AA simulations, we designed an alchemical reaction coordinate path that initially decoupled the partial electrostatic charges from every atom via 10 windows, followed by a full decoupling of the VDW interactions through an additional 21 simulations for the AA watersaturated octanol and pure water boxes. The initial structures were then energy minimized (steepest descent, 10000 steps), equilibrated in the NVT ensemble for 100 ps, followed by 100 ps NPT equilibration and a final production run of 9 ns for each window with the same free energy and soft core parameters as we have previously used for the CG free energy simulations. 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

Pressure was controlled with isotropic PR coupling with a coupling constant of 3 ps with all other simulation parameters for the AA solvent boxes the same as for the POPC-anchored peptide. Similar systems were prepared for the fatty tails systems, which we copied the previous AA peptide-solvent boxes and removed their amino backbone atoms as well as the Cβ and Sγ atoms of the modified cysteine residues. In the CG membrane anchor boxes this modification was repeated by removing their backbone beads as well as the first sidechain beads for the lipidated cysteines. We have detailed the list of our studied systems in table 1. Each completed set of simulations was analyzed and validated with the alchemical Gromacs toolset to calculate free energies of partitioning and associated errors using MBAR. 98,99 Table 1. A detailed list of all simulated systems. System

AA/CG

MF-Cpalmitoyl-IH peptide embedded in POPC bilayer MF-Cfarnesyl-IH peptide embedded in POPC bilayer MF-Cgeranylgeranyl-IH peptide embedded in POPC bilayer Gmyrsitoyl-MFCIH peptide embedded in POPC bilayer MF-Cpalmitoyl-IH peptide embedded in POPC bilayer MF-Cfarnesyl-IH peptide embedded in POPC bilayer MF-Cgeranylgeranyl-IH peptide embedded in POPC bilayer Gmyrsitoyl-MFCIH peptide embedded in POPC bilayer Cpalmitoyl in water box Cfarnesyl in water box Cgeranylgeranyl in water box Gmyrsitoyl in water box Cpalmitoyl in water box Cfarnesyl in water box Cgeranylgeranyl in water box Gmyrsitoyl in water box Cpalmitoyl in water octanol box Cfarnesyl in water octanol box Cgeranylgeranyl in water octanol box Gmyrsitoyl in water octanol box Cpalmitoyl in water octanol box

AA

Production time (ns) 100

AA

100

AA

100

AA

100

CG

2000

CG

2000

CG

2000

CG

2000

AA AA AA AA CG CG CG CG AA AA AA

9 9 9 9 1000 1000 1000 1000 9 9 9

AA CG

9 1000

12 ACS Paragon Plus Environment

Page 12 of 34

Page 13 of 34

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

Cfarnesyl in water octanol box Cgeranylgeranyl in water octanol box Gmyrsitoyl in water octanol box palmitoyl in water box farnesyl in water box geranylgeranyl in water box myrsitoyl in water box palmitoyl in water box farnesyl in water box geranylgeranyl in water box myrsitoyl in water box palmitoyl in water octanol box farnesyl in water octanol box geranylgeranyl in water octanol box myrsitoyl in water octanol box palmitoyl in water octanol box farnesyl in water octanol box geranylgeranyl in water octanol box myrsitoyl in water octanol box

CG CG

1000 1000

CG AA AA AA AA CG CG CG CG AA

1000 9 9 9 9 1000 1000 1000 1000 9

AA AA

9 9

AA

9

CG

1000

CG CG

1000 1000

CG

1000

RESULTS Our initial approximation for the bonded parameters of the modified cysteine residues and the myristoylated glycine was made with reference to the atomistic simulations of our membraneembedded MFCIH and GMFCIH peptides. Initial choices for the mapping schemes were determined from the chemical moieties which were used to parameterize the CG beads of Martini60 and were followed by an iterative process of validation and optimization of each mapping scheme and their respective bonded parameters. All of our finalized parameters have been provided in the supporting information as well as in the latest version of the martinize script (https://github.com/raziel81/martinize.py/). Angle and distance distribution plots derived from the AA simulations and their overlap with the CG simulations are available in the supporting information as well. All distributions obtained from the AA simulations approximately yield unimodal distributions that were recapitulated in the Martini forcefield59 and were fitted accordingly. We note that due 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

to Martini’s smoother potential, we could not reproduce bimodal distributions such as the BBC1-C2 angle of cysteine palmitoyl, cysteine geranylgeranyl and cysteine farnesyl in which we aimed to fit our distributions to overlap as much as was possible with the AA distribution. Further examination of the conformational landscape of our studied MAs has been carried out by calculating the RMSD matrix of each of the modified residues with respect to their initial minimized structures throughout their AA and CG trajectories and comparing the averaged structures of their computed clusters. From the AA simulations, we identified the primary structure of the three modified cysteines to be of a “kinked rod” with the position of the kink being at the thio-ester/ether bond. For the myristoylated glycine its primary structure favors a curved bow-like orientation which it sustained throughout its trajectory and was accordingly reproduced with martini (Figure 3, Tables 3-4). Interestingly, such an observation has been previously reported for DOPC acyl-chain methyl groups from neutron diffraction studies as well. 100

Furthermore, a second structure in the AA simulations was identified in each of the four MAs in which the fatty chain has curled back towards the peptide’s adjacent hydrophobic sidechains. The population of this structure was only 10% of cysteine geranylgeranyl’s conformational ensemble and 6% for cysteine palmitoyl (Figure 4), while for both cysteine farnesyl and glycine myristoyl it has a population of less than 0.5% of the trajectory’s conformational landscape. However, in the CG simulations we have only obtained a significant population of the aforementioned secondary clusters for cysteine geranygeranyl with a population of ~ 2% of its conformational ensemble. Density profiles for the fatty chains of each of the membrane anchors are shown in Figure 5. The most peripheral atoms of cysteine palmitoyl, cysteine geranylgeranyl and glycine myristoyl have a propensity for curling towards the upper leaflet’s interface with the water, especially in cysteine geranylgeranyl. However, in contrast to cysteine geranylgeranyl and cysteine palmitoyl, cysteine farnesyl’s heavy atoms produced narrower distributions. Previously, Singh and Tieleman have observed that that cysteine favors partitioning into the interface region between water and a POPC bilayer as well as between water and cyclohexane when its positioned in the center of the Wimley-White pentepeptide.101 Yet in our current work, we observed that once cysteine is lipidated, it strongly favors the hydrophobic interface instead. 14 ACS Paragon Plus Environment

Page 14 of 34

Page 15 of 34

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

In another work, de Jong et al. have run CG MD simulations on models of H-Ras and N-Ras in which H-Ras was palmitoylated at residues C181 and C184 and farnesylated at C186 while NRas was only palmitoylated at C184 palmitoylated and farnesylated at C186.48 Their results have demonstrated similar partitioning behavior for cysteine palmitoyl as in our work, in which cysteine palmitoyl favors the hydrophobic region of the bilayer. However, their study has also suggested that cysteine farnesyl partitions into the water phase while in our current study, we have shown that it partitions deeper into the hydrophobic region of the bilayer, at least within the context of a small model peptide. Furthermore, we have analyzed the penetration depths of each of the MAs in our CG simulations, with respect to the POPC bilayer center of mass, and found that cysteine farnesyl penetrates less deeply into the bilayer then the cysteine palmitoyl (Table 2) which is in qualitative agreement with the results of de Jong et al. as well.48

Table 2. Coarse-grained simulations analysis of center of mass penetration depth of the MAs with respect to the POPC bilayer COM.

Mean Standard deviation

Cysteine FarnesylBilayer COM ∆z [nm] 1.22 0.23

0.94

Cysteine PalmitoylBilayer COM ∆z [nm] 0.92

Glycine MyristoylBilayer COM ∆z [nm] 1.20

0.26

0.22

0.25

Cysteine GeranygeranylBilayer COM ∆z [nm]

Table 3. Cluster sizes of each of the POPC-embedded membrane anchors trajectories from the AA simulations.

1

Cysteine Farnesyl 99.5

Cysteine Geranygeranyl 86.0

Cysteine Palmitoyl 93.2

Glycine Myristoyl 99.9

2

0.5

10.0

5.9

0.1

3

-

2.0

0.8

-

4

-

1.1

0.1

-

5

-

0.8

-

-

Cluster #

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 34

Table 4. Cluster sizes of each of the POPC-embedded membrane anchors trajectories from the CG simulations.

1

Cysteine Farnesyl 100

2

-

1.9

0.2

-

3

-

0.1

0.1

-

4

-

-

-

-

5

-

-

-

-

Cluster #

Cysteine Geranygeranyl 98.0

Cysteine Palmitoyl 99.7

Glycine Myristoyl 100.0

Figure 3. Illustration of the central structure of each cluster as was calculated from the all-atom (A) and coarse-grained simulations (B) 70,72

16 ACS Paragon Plus Environment

Page 17 of 34

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 4. Snapshots in which representations the central structures of the two major clusters of cysteine geranylgeranyl are illustrated from all-atom simulations (A+B) and coarse-grained simulations (C+D)

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

Figure 5. All-atom density profiles of A) Cysteine farnesyl, B) Cysteine palmitoyl, C) Cysteine geranylgeranyl and D) Glycine myristoyl. The color gradient is coded so that with respect to the peptide’s backbone, the most peripheral heavy atoms are red and the heavy atoms which are the most proximal to the backbone are coded in blue.

The free energy of partitioning of specific beads between different solvents is a primary thermodynamic property to reproduce in Martini.60 Specifically, the partitioning of a solute between water and 1-octanol is considered a primary descriptor in drug development57 and has been utilized as a key property in the development of its more recent extensions for carbohydrates102 and DNA103 which have utilized an “alchemical cycle” through thermodynamic integration85 as have similarly done in our current study. Previously, Bereau and Kremer have observed a discrepancy in the free energies of water-octanol partitioning for each of the originally defined CG beads in Martini and the values which they have obtained from a similar cycle such as was previously used for Martini carbohydrates and DNA.91 In our current study, 18 ACS Paragon Plus Environment

Page 18 of 34

Page 19 of 34

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

we have repeated these calculations with our own simulation protocol and reproduced them for the individual beads that we used for all four MAs (Figure S1). All four MAs as well as their cysteine/glycine-separated fatty tails, were simulated in water saturated 1-octanol and pure water solvation boxes with Martini and were compared to their allatom structures via CHARMM36 forcefield. For the full structures of the lipidated amino acids the free energies of octanol-water partitioning with Martini deviated within the range of ~10-25 kj/mol while for the cysteine/glycine-separated tails the deviations were far smaller within the range of