Insights into Membrane Protein–Lipid Interactions from Free Energy Calculations

Integral membrane proteins are regulated by specific interactions with lipids from the surrounding bilayer. The structures of protein–lipid complexes can be determined through a combination of experimental and computational approaches, but the energetic basis of these interactions is difficult to resolve. Molecular dynamics simulations provide the primary computational technique to estimate the free energies of these interactions. We demonstrate that the energetics of protein–lipid interactions may be reliably and reproducibly calculated using three simulation-based approaches: potential of mean force calculations, alchemical free energy perturbation, and well-tempered metadynamics. We employ these techniques within the framework of a coarse-grained force field and apply them to both bacterial and mammalian membrane protein–lipid systems. We demonstrate good agreement between the different techniques, providing a robust framework for their automated implementation within a pipeline for annotation of newly determined membrane protein structures.


CGMD simulations
All systems were solvated with Martini 2.2 waters and Na + and Clions to a neutral charge and 0.15 M.
Additional elastic bonds of 1,000 kJ mol -1 nm -2 were applied between all backbone beads within 1 nm.
Electrostatics were described using the reaction field method, with a cut-off of 1.1 nm using a potential shift modifier, and van der Waals interactions were shifted from 0.9-1.2 nm. Bonds were constrained using the LINCS algorithm. For initial set up, self-assembly simulations were run to allow the membrane lipids to assemble into a bilayer around the protein. These were run over 100 ns in the NPT ensemble, with V-rescale temperature coupling at 323 K 1 and semi-isotropic Berendsen pressure coupling at 1 bar 2 . Following this, the systems are then simulated for an extended period (>3 µs) to permit the formation of specific proteinlipid interactions. These and all subsequent simulations used V-rescale temperature coupling at 323 K 1 and semi-isotropic Parrinello-Rahman pressure coupling 3 , unless specified otherwise. All simulations were run using Gromacs 2018 4,5 unless otherwise specified

Potential of mean force calculations
Steered MD simulations were run to allow construction of the 1D CV. For these, an umbrella pulling force of 1,000 kJ mol -1 nm -2 was applied along the desired CV (e.g. along the x-axis), with a rate of 0.01 nm ns -1 , defined in the mdp file as: and GL0, PO1, GL1, GL2, PO2, GL3 and GL4 beads of CDL. Umbrella potentials were applied to each independent simulation, using the code as described above but with pull_coord1_rate set to 0.
To prevent any rotational movement of the proteins in the PMF calculations, xy positional restraints of 100 kJ mol -1 nm -2 were applied to 3-4 backbone beads in each system: Kir2.

Free energy perturbation
For each FEP calculation, the target lipid was alchemically perturbed into a bulk lipid. For PIP2 to POPC, this involved a removal of beads PO1, PO2, RP1, RP2 and a conversion of bead RP3 from Martini type SP1 to Q0, including changing the charge from 0 to +1. The charge was also removed from 5 sodium beads to keep the system neutral at each λ state.
For CDL to POPC (mitochondrial) and POPE (bacterial), we initially tried splitting the CDL molecule into two separate residues, and growing in and converting the head group beads accordingly. This is the simplest path to follow, however removal of the bond between the residues proved difficult to remove without incurring very high energies, and the energies we obtained were highly variable. Instead, we chose to fully remove two of the acyl CDL chains (i.e. beads PO2, GL3, GL4, C1C, C2C C3C, C4C, C1D, C2D, D3D, C4D, C5D), and changed bead GL0 from type Nda to either Qd (for POPE) or Q0 (for POPC), switching from a 0 to +1 charge in both cases.
We used a dual topology approach, and FEP calculations were run using the following section added to the production mdp file: Absolute binding free energy calculations The details of the restraints used on the A2AR-cholesterol ABFE calculations are outlined below.
where RES_1 and RES_2 are reference positions on the receptor, Dum1 and Dum2 are the dummy particles and a_1 and a_8 are the first and last beads of the cholesterol molecule. Note that the dummy particles have the same initial coordinates as the a_1 and a_8 beads. As the receptor diffuses in the membrane, the dummy particles move with it without introducing additional energy into the system.
The RMSD restraints can then be established between the cholesterol particles and the dummy particles, using the following mdp setting: = 100 pull_coord8_rate = 0 ; repeat for groups [4][5][6] Restraint FEPs were then computed as described above.
ΔΔGbind values were calculated using the thermodynamic cycle in Supporting Figure