Balancing Force Field Protein-Lipid Interactions To Capture Transmembrane Helix-Helix Association.

Atomistic simulations have recently been shown to be sufficiently accurate to reversibly fold globular proteins and have provided insights into folding mechanisms. Gaining similar understanding from simulations of membrane protein folding and association would be of great medical interest. All-atom simulations of the folding and assembly of transmembrane protein domains are much more challenging, not least due to very slow diffusion within the lipid bilayer membrane. Here, we focus on a simple and well-characterized prototype of membrane protein folding and assembly, namely the dimerization of glycophorin A, a homodimer of single transmembrane helices. We have determined the free energy landscape for association of the dimer using the CHARMM36 force field. We find that the native structure is a metastable state, but not stable as expected from experimental estimates of the dissociation constant and numerous experimental structures obtained under a variety of conditions. We explore two straightforward approaches to address this problem and demonstrate that they result in stable dimers with dissociation constants consistent with experimental data.

Introduction

Although the folding of globular proteins is now relatively well-understood,[1] the mechanisms of membrane protein folding are less well characterized. Both from an experimental and from a computational perspective, membrane protein folding is much more challenging. Diffusion within a lipid bilayer is 2–3 orders of magnitude slower than in water, posing a problem for computer simulations, where the time scales accessible are limited.[2] Even obtaining reversible folding in experiments is more challenging than for globular proteins.[3] In addition, the mechanism by which membrane proteins are inserted and folded into membranes in cells involves numerous chaperones and membrane insertion machinery.[4] Despite these challenges to studying membrane protein folding, it is of considerable biomedical relevance: transmembrane domains play critical roles in many pharmaceutical applications, with G-protein coupled receptors being the most obvious example. Furthermore, there are many diseases associated specifically with misfolding of membrane proteins, with the cystic fibrosis transmembrane conductance regulator being a prime example.[5]

The prevalent model for spontaneous in vitro membrane protein refolding is the “two-stage model”, proposed by Popot and Engelman.[6−8] In this model, the helices first insert and fold in the membrane, followed by association and formation of tertiary packing interactions. It also appears that the mechanism by which helices are inserted into the inner membrane by the translocon has many similarities to this simple model.[4] The first of these processes, insertion, is accompanied by a favorable change in free energy and has been observed to occur spontaneously in simulations at high temperature.[9] The second process, helix association, is harder to study by simulation due to the high viscosity of the membrane, and to date has not been observed in unbiased all-atom simulations.

Perhaps the simplest and best characterized prototype for membrane protein folding in the assembly is the glycophorin A homodimer, which is formed by the association of two single-pass transmembrane helices. The small system size is amenable to all-atom simulation, whereas the wealth of available experimental data makes it suitable for assessing the accuracy of computational models for folding. Glycophorin A dimer structural models have been experimentally obtained by X-ray crystallography on samples derived from lipid cubic phase (LCP),[10] as well as by solid-state[11] and solution[12] NMR (Figure A). All experimental structures, obtained under a variety of conditions and with different membrane compositions, provide a highly consistent structural picture of a symmetric homodimer packing via a conserved GXXXG motif[13] (Figure A). Affinities estimated from analytical ultracentrifiguration[14] and FRET-based experiments[15−17] point to a stable dimer in a variety of lipid- and lipid-mimetic conditions, with the free energy of association spanning a rather wide range, from −16 to −51 kJ/mol[15] (using a reference concentration of 1/nm2), depending on the conditions.

Figure 1

(A) Experimental glycophorin structures, showing only residues Ser-69 to Arg-97 used in simulations. Structures are aligned on the PDB 5E4H crystal structure using backbone of residues 78–88: this ten-residue stretch was used to define the DRMS coordinate and crossing angle (the stretch is indicated by a straight black line next to the structure render and, repeated, under the sequence). For DRMS definiton, see eq in the Methods section. The solution NMR structure in a micelle has PDB code 1AFO, shown in blue; the lipid cubic phase crystal structure has PDB code 5EH4, shown in red; the solid state (ss) NMR structure by Smith et al. is shown in gray.[11] (B) Top- and side-view of the simulation box used. Water is shown as a continuous surface, lipids are not shown for clarity, and the protein is in cartoon representation. Note that the oblique viewing angle in the lower figure may make the protein appear to protrude into the solvent more than it does.

Given the abundant experimental data available for the glycophorin A dimer, a range of computational approaches has been taken to study it. Initially, implicit membrane models were used to characterize the dimerization energy landscape. For example, a promising study using an implicit membrane model by Sugita and Im successfully identified the native structure of GpA dimer via replica exchange simulations of the bound dimer.[18] Similarly, the implicit membrane model (IMM) from Lazaridis has been used to study GpA dimerization.[19] With an explicit representation of the membrane, a number of coarse-grained approaches have also been taken. Janosi et al. used a Monte Carlo-based approach with the coarse-grained MARTINI potential to study GpA dimerization.[20] Psachoulia used unbiased coarse-grained MD simulations to study the contacts mediating the GpA dimerization.[21] Sengupta and Marrink also computed a potential of mean force (PMF) of glycophorin A using umbrella sampling with MARTINI the force field.[22] Recently, a novel sampling approach was proposed by Hummer and co-workers to study the Mga2 helix-dimer behavior in the membrane[23] using MARTINI. Umbrella sampling at coarse-grained resolution can also be used to study dimerization of soluble protein dimers.[24] At a mixed coarse-grained and all-atom resolution, the PACE 1.5 force field combined a united-atom description of the protein with coarse-grained MARTINI lipids and solvent and was shown to result in contact formation between GpA helices in short, unbiased simulations,[25] although a free energy surface was not determined.

While the coarse-grained models have the advantage of being computationally inexpensive, an all-atom force field for both protein and membrane is the most detailed description and accurate model that is practical for running simulations and therefore expected to be most predictive for a range of systems. An all-atom description would be expected to capture better the native binding mediated by the GXXXG motif in GpA dimerization.

Toward this end, a landmark study using all-atom simulations by Chipot et al. employed a membrane-mimetic dodecane “slab” to speed up the “membrane” diffusion of glycophorin A in simulations, together with the adaptive biasing force enhanced sampling method to obtain a potential of mean force for helix association.[26] At the fully atomistic level, Pastor and co-workers used long, unbiased molecular dynamics simulations on the ANTON supercomputer to demonstrate the stability of the related ErbB1/B2 and EphA1 transmembrane helix dimers on a 1 μs time scale.[27] Most recently, Kuznetzov et al. used umbrella sampling to estimate the energy landscape of GpA dimerization in a POPC bilayer as a function of helix center-of-mass distance with the GROMOS 43a2 force field.[28]

Many of the above studies have determined free energies and PMFs for the association of the glycophorin and related TM helix dimers. Such calculations are challenging for several reasons, particularly when undertaken in membranes. The first is the choice of an appropriate reaction coordinate. PMFs determined by biasing along the interhelix center of mass distance, the most commonly used coordinate, typically show a deep minimum at short distances. However, it is known that distance alone is insufficient as a reaction coordinate for helix dimerization as it does not account for different helix–helix packing at close distances and does not capture the helix–helix crossing angle.[18] Recently, we have performed umbrella sampling simulations to determine the dimerization free energy for glycophorin A using the MARTINI 2.1 coarse-grained force field,[29] finding a deep minimum in the PMF for interhelix distance. However, as we have shown, most of this minimum was composed of non-native states in which the helices were approximately parallel rather than at the native crossing angle of −40° seen consistently across multiple experimental structures of GpA.[10] Steric trap experiments, used to determine the GpA association constant, cannot exclude the existence of non-native bound states. However, the fact that GpA dimerization-disrupting mutants cause a loss of binding in these experiments suggests that non-native bound state has a negligible population relative to native. We showed that an alternative one-dimensional reaction coordinate, the interhelical distance matrix RMSD, was able to distinguish a minor population of the correct native state in a coarse-grained simulation.[29]

A second critical challenge to determining PMFs is that of obtaining equilibrium sampling from simulations in a viscous membrane environment. In our earlier work, we used a stringent criterion for assessing sampling, namely, that the PMFs obtained from two different initial conditions (corresponding to natively bound and well-separated helices) should “converge” to the same result.[29] It turned out that very long umbrella sampling simulations (on a microsecond time scale or longer) were needed to satisfy this condition with the MARTINI model.

Here we investigate the free energy landscape for glycophorin A dimer formation with the widely used all-atom CHARMM36 force field.[30,31] We chose CHARMM36 because it is known to result in an excellent reproduction of membrane properties,[30,32] as well as being a high-quality force field for globular proteins.[31,33,34]

We use the interhelical DRMS to the native structure as a collective variable (see Methods section for the definition) for calculating helix–helix dimerization PMFs and we compare GpA dimerization PMFs to available experimental data. We obtain similar PMFs starting from both unbound and bound helices.

An important feature of all our PMFs is a substantial (10–20 kJ/mol) energy barrier for native dimerization, not seen in previous work. This barrier arises from the disruption of the tightly packed GXXXG motif and is distinct from the barrier for full dissociation. Perhaps surprisingly, we also find that the native dimer is unstable in phosphatidylcholine (POPC) membranes, in contrast to the available experimental Kd data.

We have tested several possible modifications of the force field to explore and address this discrepancy. The most promising approach is a marginal reduction of dispersion interactions between protein and lipids. This small correction is sufficient to render the native dimer the free energy minimum and results in reasonable agreement with available experimental dissociation constants. For reference, we compare the PMF derived from CHARMM36 to representative membrane force fields at the coarse-grained level of resolution (MARTINI 2.1). The all-atom model for protein and membrane (CHARMM36) results in a stable native dimer, which is also more stable than competing non-native dimers.

Methods

Molecular Simulation Methods

GROMACS version 4.6.7 (http://www.gromacs.org) was used for simulations,[35] with the PLUMED 2.2-hrex patch to enable the enhanced sampling functionality where needed.[36] Pressure was maintained at a reference of 1 bar via a Parrinello–Rahman barostat[37] with independent control in the X, Y, and Z directions, and with a coupling constant of 5 ps. Temperature was maintained via stochastic velocity rescaling[38] at 310 K with a relaxation time of 1 ps. Shifted Lennard-Jones interactions were cut off at 1.2 nm. Long-range Coulomb interactions were calculated with the particle mesh Ewald (PME) method,[39] using a grid spacing of 0.12 nm and a real-space cutoff of 1.2 nm.

NAMD simulations were performed with version 2.12,[40] in the NVT ensemble, using similar settings to GROMACS. Nonbonded interactions were truncated at 1.2 nm, with a force-based switching function applied between 1.0 and 1.2 nm, and the nonbonded pair-list was updated every 10 time steps, with a neighbor search cutoff of 1.6 nm. The PME method as used to compute electrostatic interactions, with a grid spacing of 0.1 nm and a sixth-order spline. Langevin dynamics with a friction of 1.0 ps–1 was used, with a time step of 2 fs.

Temperature replica-exchange simulations[41] were performed using the standard GROMACS implementation, spanning the range of temperatures shown in Table , and solute-tempering (REST2)[42] was performed using GROMACS 4.6.7 with the PLUMED 2.2-hrex patch. Protein and palmitoyl oleoylphosphatidylcholine (POPC) lipids were treated as the “hot” atoms in REST2 simulations. Exchanges between replicas were attempted every 1000 MD steps, the observed exchange acceptance probability was 0.05–0.20.

Table 1

Replica Exchange Simulation Setups Used in This Work

code	force field	method and ensemble	ladder spacing	components
Gromacs 4.6.7	CHARMM36	T-REMD NPT	45 windows, 310–450K	78 POPC
Gromacs 4.6.7	CHARMM36	T-REMD NPT	45 windows, 310–450K	78 POPC
				150 mM NaCl
Gromacs 4.6.7	CHARMM36m	T-REMD NPT	45 windows, 310–450K	78 POPC
Gromacs 4.6.7	CHARMM36	REST2 NPT	16 windows, 300–400K Teff	112 POPC
Gromacs 4.6.7	AMBER	REST2 NPT	45 windows, 300–450K Teff	112 POPC
NAMD 2.12	CHARMM36	T-REMD NVT	40 windows, 300–450K	78 POPC

The all-atom CHARMM36 protein and lipid force field was used in all simulations[43] together with TIP3P-CHARMM water model, except where otherwise noted.[44,45]

Glycophorin A All-Atom Dimer System

The initial wild-type structure of glycophorin was taken from PDB entry 1AFO.[12] Residues Ser-69 to Arg-97 were used in the simulations, with the flexible termini removed. The protein helix–helix dimer was inserted into a square, solvated POPC bilayer giving a final system with 112 lipid molecules (equal in each leaflet) and 3999 water molecules. The resulting system was approximately 6.4 nm in both the X and Y dimensions (Figure ).

For some simulations (where specified) a smaller, rectangular bilayer was used: the system then contained 78 POPC molecules, equally distributed between the two leaflets, and 3998 water molecules. This system was approximately 3.2 × 6.4 nm in the X and Y dimensions, respectively. Unless otherwise stated, the larger, square system with 112 POPC lipids was used in the simulations.

Potential of Mean Force Calculations

The glycophorin fragment mediating the tight helix–helix packing via the GXXXG motif was used to define the native DRMS collective variable. Heavy, nonsymmetric atoms of residues 78–88 were used to define the distance matrix. Only atom pairs with native distances between 0.1 and 0.6 nm were included in the definition of DRMS, defined aswhere d(x,x) denotes the distance between the coordinates x and x of atoms i and j in configuration X and summation runs over all atom pairs (i,j) specified in the in the native contact list, {contacts}.

An umbrella sampling simulation was setup with 40 replicas spaced between DRMS of 0 and 2.5 nm; the positions of umbrellas and spring constants of each umbrella are given in Supporting Information Table S1. All umbrella windows were run simultaneously with replica exchange moves between them to assist sampling of orthogonal coordinates. The crossing angle was defined as the angle between the vectors connecting the Cα atoms of residues 78 and 88 in each helix.

The WHAM method was used to perform unbiasing of the umbrella sampling trajectories,[46] using the Grossfield lab implementation, version 2.0.9 (http://membrane.urmc.rochester.edu/content/wham). The second half of each umbrella sampling simulation was used for unbiasing.

The dissociation constant was evaluated using the following equation:where Kd is in units of molecules·nm−2 and F(r) is the PMF for association on the radial center-of-mass distance coordinate r, with the entropy correction 2kBT ln(r) added. With that correction, there should be some distance b above which F(r) is constant, which defines the bound state. It is assumed that F(r) = 0 for large distances, so that the above integral converges, which can be ensured by adding a suitable constant to F(r).

Analysis and Visualization

Visualization was done using VMD,[47] and analysis was done using MDAnalysis.[48]

Results and Discussion

Sampling Dimer Formation

To determine the free energy surface for glycophorin dimerization, we perform umbrella sampling simulations of the dimer along the interhelical DRMS coordinate as described in our earlier work.[29] Two sets of simulations were performed: one started from the native, bound dimer (called “together”), and another set with all the replica windows initialized from the unbound configuration (called “separate”). Over the course of the simulations the helices remained stable, with the number of helical residues close to the experimentally determined structures (Supporting Information Figure S8). The resulting PMFs are shown in Figure A,B, projected onto the DRMS and interhelix distance coordinates, respectively. The PMFs started from the two initial conditions are challenging to converge, despite nearly 0.5 μs simulation per replica: the high viscosity of the lipid bilayer and the resultant slow diffusion in the plane of the membrane make this a very challenging sampling problem. However, it is clear that the two PMFs are converging, as indicated by monitoring the free energy difference between the bound and unbound states (Supporting Information Figure S1). Despite the incomplete sampling, the PMFs started from the “together” and “separate” initial conditions have consistent features.

Figure 2

Glycophorin A dimerization potentials of mean force along DRMS and center of mass (COM) distance collective variables (A, upper and lower panel, respectively). The thick curve is the result of pooling sampling from different initial conditions; thin curves demarcating shaded region are the obtained independently from each set of initial conditions (“together” and “separate”). Representative structures for three states (labeled 1–3) are in cyan, with the native structure in red and blue for reference (B). See also Supporting Information Figure S2 for a larger ensemble of representative structures, and Supporting Information Figure S4 for crossing angle distributions for each of the mimina.

There are two unexpected features of these PMFs. First, it appears that the most stable form of the glycophorin A in this force field is the dissociated state (Figure A), approximately 20–30 kJ/mol more stable than the native dimer state for the system size used, irrespective of the initial condition. Using eq , we obtain a dissociation constant of ∼10 molecules.nm–2, well outside the range of experimental estimates[15] of 10–3 to 10–9 molecules.nm–2.

Second, the native and non-native bound states (labeled “1” and “2” in Figure ) are separated by a substantial energy barrier, of approximately 10 kJ/mol on the DRMS coordinate. The barrier corresponds to the dissociation of the tightly packed GXXXG motif to form a non-native dimer state. Such a barrier was not observed in previous PMFs using the interhelical center of mass (COM) distance as a reaction coordinate[28] but may help to explain why a variety of transmembrane helix dimers are at least metastable in equilibrium simulations at 300 K.[27,49] Using one-dimensional Kramers theory[50] to estimate the rate k = Dωbω‡ exp[−Δ/kBT]/2π, with a barrier height Δ of 10 kJ/mol and curvatures of the bound state ωb and barrier ω‡ of ωb2 = ω‡2 = 200 kJ/mol/nm2 estimated from the PMF, and a helix lateral diffusion coefficient of 0.3 nm2/μs,[51] yields an approximate unbinding rate at 300 K of only ∼0.2 μs–1.

We observe a small barrier of ∼5 kJ/mol on the helix–helix COM distance coordinate in Figure B, relative to the ∼10 kJ/mol. The reason for the lower barrier along the COM distance is that it is degenerate relative to the interhelical DRMS coordinate. The apparent barrier at a helix–helix COM distance of around 0.7 nm in fact overlaps with a free energy minimum along DRMS (Figure A). That said, DRMS is also imperfect and appears to be degenerate with respect to the helix–helix crossing angle at a DRMS of around 0.25 nm (see minimum labeled with an asterisk in Figure B), suggesting that the true barrier to dissociation may be higher than 10 kJ mol–1. The smaller barrier is seen for the COM distance, which helps to explain why no barrier was seen in earlier work using that coordinate. Alternatively, the lack of a barrier in other studies may originate from the different force fields or membrane models used, or due to limited sampling. We note that this barrier arises from more than the entropic restriction arising from the requirement to orient the two helices correctly to bind in the native orientation. In Figure S7 we show the free energy surface projected in two dimensions onto the relative helix orientation and the DRMS: even in this projection, the barrier around the native minimum is still visible.

Figure 3

2D PMF projections along the helix–helix COM distance (A) and helix–helix crossing angle (B) against interhelical DRMS to native. A degenerate, non-native state has a center-of-mass distance of 0.7 nm and is labeled by “D” in (A). A non-native, positive crossing-angle state is labeled by an asterisk in (B). For crossing angle distributions of minima 1–3, see Supporting Information Figure S4. For additional projections, along helix–helix rotation angles, see Supporting Information Figure S7.

Because of the unexpected observation that the glycophorin A dimer was unstable in our simulations, we have performed additional simulations, using other enhanced sampling methods (rather than umbrella sampling) and also multiple simulation codes. To check the effect of the sampling method, we have also run glycophorin A simulations from the native state in the absence of an umbrella potential, in replica exchange simulations. The dimer embedded in a POPC bilayer appears unstable in the first replica, corresponding to physiological temperature, with dissociation occurring on a time scale of 50–300 ns (Figure ). This instability does not appear to depend on the particular method used: both standard temperature replica exchange (T-REMD) and solute-tempering replica exchange (REST2) result in the dimer dissociating on a 10–100 ns time scale (Figure ). Furthermore, the effect appears to be independent of system size: the larger square bilayer system and smaller rectangular systems both show similar instabilities (Figure ). Similar results are obtained with the NAMD simulation package, in which the CHARMM nonbonded cutoffs can be precisely replicated; in this case, the T-REMD runs were performed at constant volume owing to limitations of the standard REMD script provided with NAMD (pressure and volume changes between replicas are not included in calculating exchange probability), but the effect is reproduced nonetheless. Inclusion of 150 mM sodium chloride also does not change the result.

Figure 4

Glycophorin A dimer dissociates in simulations at physiological temperature (310 K). Fraction-of-dimer (“dimer” is defined as inter-helix DRMS < 1.0 nm) time series, with data smoothed by a running average over 5 ns (blue) or 10 ns (green) windows, from replica exchange simulations. Runs using the GROMACS simulation code with the CHARMM36 and CHARMM36m force fields in the NPT ensemble result in dimer dissociation with the T-REMD method (upper row). The same behavior is observed for CHARMM36 with NaCl ions at 150 mM concentration, or CHARMM36m[34] run with REST2 without ions (middle row, left and right column, respectively). Dissociation also occurns with the AMBER Slipids force field in GROMACS (lowest row, left column) and with the CHARMM36 force field in NVT REMD simulations using NAMD (lowest row, right column). Further details are available in Table .

Lastly, we considered whether the dimer instability was specific to the CHARMM36 force field or was also found in other force fields. A similar result is obtained with the recently published modification CHARMM36m.[34] Another widely used lipid force field is the AMBER Slipids model.[52,53] We find that combining Amber ff03w with Slipids seems slightly more stable but still results in dissociation after about 250 ns of REST2 simulation. The appearance of dissociated helix monomers at the physiological-temperature replica is inconsistent with available experimental data. Though a much longer simulation would be needed to get a converged estimate of the population of unbound dimer in the first replica, the available experimental data suggests that there should no detectable dissociated fraction observed for wild-type glycophorin A TM at physiological temperature.

It is known that bilayer properties obtained through simulations are somewhat sensitive to simulation parameters chosen, including both the long-range electrostatic interactions and the long-range contribution to the dispersion interactions. For example, Vattulainen et al.[54] observed freezing of DPPC bilayers when truncated electrostatics was used; using PME electrostatics allows the correct properties to be recovered. More recently, Mark at el. showed that using the new implementation of the GROMACS twin-range cutoff for van der Waals and Coulomb interactions results in DPPC forming a liquid-ordered phase above its transition temperature.[55] In these examples, it must be noted that DPPC is particularly sensitive to small changes in parameters as a consequence of being close to a phase transition at typical temperatures of interest, whereas the current simulations used POPC, well away from its phase transition temperature (Tc = 271 K).

Given the above sensitivity, a possible concern with the REST2/HREX methodology[42] is its effect on membrane properties, particularly in the replicas where the force field is most perturbed. The method applies scaling factors to the nonbonded (van der Waals, electrostatic) terms involving protein and lipids in the force field. Decreasing the scaling factor λ from a value of 1, which corresponds to the initial force field, is commonly thought as effectively making an affected group of atoms hotter, which can also be expressed in terms of an effective temperature, Teff (Supporting Information Table S2). One replica is always run with λ = 1.0, corresponding to the true force field at the temperature of interest (known as the neutral replica). To check that this approach reproduces equilibrium properties in the neutral replica and does not introduce artifacts in other replicas, we have performed long MD simulations, without replica exchange, of pure POPC bilayers with a range of scaling factors λ, with lipids as the hot group. We have checked a number of benchmark bilayer properties: area per lipid, lateral diffusion coefficient, and lipid order parameters, all of which appear to be within experimental error for the neutral replica (Supporting Information Figure S5A–C). As the scaling factor λ was decreased, the bilayer diffusion coefficient increased sharply (the reason REST2 is so effective), the disorder of the lipid tails increased, and the thickness and area per lipid decreased and increased, respectively. Importantly, however, the bilayer remained intact for all values of λ and no obvious artifacts were observed.

Force Field Adjustments

To address this force field limitation and propose a conservative correction, we consider how the current all-atom force fields are developed. In all-atom force fields, the protein and lipid parameters have typically been developed somewhat independently, through separate optimization to match experimental observables for lipid/water or protein/water systems. Provided that the overall framework of parametrization is similar for both the lipids and protein, the properties of mixed protein/lipid systems would be expected to be reasonable, but they are not used directly in parametrization.

The quality of a protein–lipid force field can be measured by comparing to experimental data.[56−58] However, this remains quite challenging due to the sparsity of available experimental data for mixed (i.e., protein/lipid) systems for which it is easy to obtain converged simulation results. Frequently used data reflecting partitioning of side-chain analogs or amino acids into the membrane does not directly test the protein–protein interactions within the bilayer and may also not be representative of the situation for larger peptides and proteins.[59−62] Thus, our aim in proposing a force field modification is to do the minimum needed to get a helix dimer whose stability is consistent with experiment.

Though it is possible that, for example, changing the lipid force field could lead to changes in structural properties of the membrane and indirectly stabilize the dimer, the lipid force field is already well established and appears to be in excellent agreement with available bilayer-related experimental data probing the bilayer structure.[32,63,64] We therefore do not think that the fault lies with the membrane properties and do not want our modification to affect the already well-optimized properties of lipid/water systems. For similar reasons, we do not want to modify the properties of protein/water systems.[31,33,43] This leaves only the possibility of changing parameters affecting protein–lipid interactions. Therefore, we hypothesize that there is a mismatch in protein–lipid nonbonded interaction strength, and we use a simple one-parameter scheme to adjust this energy scale by using the modified combination rulewhere ϵ is the Lennard-Jones contact energies for atom i, κ is a scaling factor, and ϵ is the off-diagonal Lennard-Jones contact energy for atoms i and j (i and j are protein and lipid atoms, respectively). Thus, the standard combination rule is adjusted by a factor κ, only for atom pairs in which one atom is a lipid type and one a protein type.

We then reweighted the potential of mean force for glycophorin dimerization, applying a scaling factor κ of 0.9 or 0.95 for the protein–lipid Lennard-Jones interactions (Figure A). We observed better agreement with experimental data and stabilization of the bound part of the PMF, including the native minimum. Because we sought only a minimal change, we did not pursue the use of a smaller scaling factor (e.g., κ = 0.8). Because the reweighting was quite noisy for κ = 0.9, we have confirmed that the correction has the desired effect by resampling with κ = 0.9, and we find the results are in agreement with the prediction from reweighting (Supporting Information Figure S6). With the altered protein–lipid interactions, the Kd of ∼10–3 molecules per nm2, just within the experimental range, reflects favorable association (Supporting Information Figure S3).

Figure 5

Adjusting the force field parameters to affect glycophorin A dimerization in simulation. Stabilizing the native dimer state by (A) scaling down protein–lipid van der Waals interaction (epsilon term) by factors of 1.0, 0.95, and 0.9. Scaling 0.9 was obtained from an independent simulations, scaling at 0.95 was obtained by reweighting the 1.0 simulation. The reweightings derived from the 1.0 and 0.9 PL-scaling runs are consistent (Supporting Information Figure S6). (B) Introducing a Cα–hydrogen bond-like term: reweighting simulations run at PL-1.0 (original force field, no scaling) with 0, 1, or 2 kJ/mol per formed Cα–hydrogen bond between the monomers. (C) Average number of protein–protein and protein–lipid contacts as a function of the collective variable (a “contact” was defined using heavy atoms only, on the interhelical distance pairs, within 5 Å).

We acknowledge the ad hoc nature of the applied force field modification. However, although changing the protein–lipid interaction strength may slightly affect other properties such as amino acid partitioning free energy into lipid bilayers, the changes we propose are still small relative to the large free energies typically associated with bilayer insertion. We have also considered whether more specific interactions may be responsible for the deviation from experiment. For example, Cα–hydrogen bonding has been suggested[65] to be important in stabilizing the GXXXG motif identified in transmembrane helix dimers[66] and other helix:helix interaction motifs in membranes.[67] These Cα–hydrogen bonds form between Cα–H hydrogen and acceptors with a π system, such backbone carbonyl groups.[68,69] In glycophorin A, such interactions can be identified in the GXXXG dimerization motif. To test whether including an additional energy term for such interactions could stabilize the native dimer, we have reweighted our potentials of mean force assuming a 1, 2, and 5 kJ/mol energetic gain per Cα–hydrogen bond formed. Hydrogen bond formation was defined as a protein Cα–hydrogen being within 0.32 nm of either a hydroxyl or carbonyl oxygen in the protein. This correction also stabilizes the native minimum, but in a qualitatively different way: whereas the scaling of protein–lipid interactions has a more global effect on the PMF, the hydrogen bond adjustment only affects the PMF in the immediate vicinity of the native state and has no effect on the rest of the PMF (Figure B). We cannot say only on the basis of glycophorin A TM region whether this type of correction will be applicable; that would require testing against a wider range of examples including those where Cα −H bonding does not play a role in stabilizing interfaces. Nonetheless, the magnitude of the required correction seems reasonable, because Cα −hydrogen bonds in proteins have been estimated to be in the range 7.9–10.6 kJ/mol[70] and some part of this must already be accounted for by the existing electrostatic terms in the force field.

Although the all-atom PMF obtained with CHARMM36 identified the native state as the most stable dimerized state, it nonetheless found that glycophorin A dimerization was overall unfavorable. In contrast, in an earlier study using the coarse-grained MARTINI 2.1 force field, we obtained the opposite result: the most stable state was dimeric, but most of the bound population was non-native and inconsistent with experimental evidence from mutagenesis.[29] Achieving a stable native dimer is relatively straightforward, as we have shown by some simple adjustments of the force field parameters; however, correct discrimination between native and non-native bound states appears to emerge more naturally in the all-atom model (Figure ), even though the native bound minimum is only slightly more stable than the non-native bound state in the modified force field.

Figure 6

Comparing the dimerization free energy surface for different force fields. 2D PMF projected along interhelical DRMS for (a) coarse-grained MARTINI force field, (b) all-atom CHARMM36 force field, and (c) all-atom CHARMM36 force field, with protein–lipid interaction scaled by 0.9.

Conclusion

For studying membrane protein folding in all-atom detail, it is essential to correctly capture the experimentally determined behavior of well-characterized reference systems. Membrane partitioning experiments for small molecules help to benchmark the free energy of membrane insertion of peptides and proteins.[59−62] Similarly, glycophorin A, as a minimal model of transmembrane helix association, is a prototype for the folding of helical membrane proteins, which are assembled from similar helix–helix interactions. It is also very challenging to obtain equilibrium properties for these systems, due to the bilayer viscosity. We have used a range of state of the art simulation methods to address the challenging problem of sampling glycophorin helix association.

We find that the CHARMM36 all-atom force field correctly preserves the helical stability of the dimer and can capture the native minimum, a feature that was not achieved with the coarse-grained models we considered. With the further introduction of a protein–lipid interaction correction, we are able to obtain the native state as the global free energy minimum, consistent with experimental data on wild-type glycophorin A dimer.

We were able to achieve a similar correction by introducing an energetic reward for formation of aliphatic hydrogen bonds. Although we favor the protein–lipid scaling, it is hard to discriminate between the two options based on currently available data. We anticipate that detailed information on binding kinetics should help to choose between these options.

A more extensive reparametrization of protein and lipid force fields would also be beyond the scope of the present work. Rather, our results suggest that glycophorin stability provides a new benchmark for assessing the accuracy of future force field developments.