Interplay of G Protein-Coupled Receptors with the Membrane: Insights from Supra-Atomic Coarse Grain Molecular Dynamics Simulations

G protein-coupled receptors (GPCRs) are central to many fundamental cellular signaling pathways. They transduce signals from the outside to the inside of cells in physiological processes ranging from vision to immune response. It is extremely challenging to look at them individually using conventional experimental techniques. Recently, a pseudo atomistic molecular model has emerged as a valuable tool to access information on GPCRs, more specifically on their interactions with their environment in their native cell membrane and the consequences on their supramolecular organization. This approach uses the Martini coarse grain (CG) model to describe the receptors, lipids, and solvent in molecular dynamics (MD) simulations and in enough detail to allow conserving the chemical specificity of the different molecules. The elimination of unnecessary degrees of freedom has opened up large-scale simulations of the lipidmediated supramolecular organization of GPCRs. Here, after introducing the Martini CGMD method, we review these studies carried out on various members of the GPCR family, including rhodopsin (visual receptor), opioid receptors, adrenergic receptors, adenosine receptors, dopamine receptor, and sphingosine 1-phosphate receptor. These studies have brought to light an interesting set of novel biophysical principles. The insights range from revealing localized and heterogeneous deformations of the membrane bilayer at the surface of the protein, specific interactions of lipid molecules with individual GPCRs, to the effect of the membrane matrix on global GPCR self-assembly. The review ends with an overview of the lessons learned from the use of the CGMD method, the biophysical−chemical findings on lipid−protein interplay.


INTRODUCTION
Transmembrane signaling through G protein-coupled receptors, GPCRs, is involved in many fundamental biological processes. GPCR signaling is initiated by the activation of the transmembrane receptor (R) by an extracellular stimulus, followed by the binding of its cognate G protein on the intracellular side triggering the exchange of GDP for GTP within the G protein, which in turn leads to a multitude of intracellular chain reactions. 1 Understanding GPCR signaling is an inherent multiscale problem. On one end of the scale, structural details of the activation mechanism of both the receptor (R) and the G protein (G) are important as emphasized by numerous recent studies using extensive dedicated experimental and computational resources on the conformational states of R, 2−6 G, 7,8 and R−G, 9,10 picturing a dynamic process. 5,11−14 In addition, highresolution spectroscopy and computational studies have looked into the details of ligand binding in relation to activation. 12,15−18 On the other end of the scale, receptor signaling may be described on a cellular level using a system biology approach based on analytical models. 19−22 These models have recently been dealing with the incorporation of GPCR oligomerization. 19 The challenge of connecting the atomically detailed scale to the system's level remains. Additional methods have emerged to probe this intermediate range of length and time scales. In the field of conventional experiments, labeling and singlemolecule methods have shown a lot of progress; for a recent review see ref 23. In the field of computational modeling, coarse grain molecular dynamics (CGMD) has become an important tool to study biomolecular processes, 24−26 including the lateral organization of membrane proteins. 27 In this review, we focus on CGMD studies of GPCRs based on the Martini model, 28 which, as we will show, has proven to be a very powerful tool to study GPCR signaling. The nature of the model is specifically pertinent to length and time scales where protein−lipid interplay is important and therefore makes it unique for this type of studies. Although other CG approaches exist none have been used in the study of GPCRs in lipid membrane with such resolution. The applications using the model have concentrated on the effect of lipid−protein interactions on GPCR behavior in biological membranes. GPCR oligomerization state is at the heart of most of the studies discussed in this review since it has been shown to affect their function quite dramatically, although the relevance in vivo is still a matter of debate. 29−32 We first describe the foundations of the Martini CG in MD simulations (referred to as Martini CGMD method, including the associated MD setups typical of such simulations) and associated techniques, underlining the potential of this approach and the main limitations. Then we describe in detail a few series of articles, which were grouped by theme but often coincide with the laboratory of origin and a few more isolated studies. We follow by discussing the lessons collectively learned on both the methodological and the membrane protein biophysics principle perspectives. We end the review with short concluding remarks and perspectives in the field.
Other recent reviews in this field include reviews on the progress of molecular dynamics simulations for the study of GPCRs, 33−35 a review of more general computational methods for predicting the structure of GPCRs and their interactions with ligands with an emphasis on allosteric issues, biased signaling and oligomerization, 36 reviews covering the different length and time scales relevant to GPCR signaling, 37,38 and some reviews on GPCRs related to their interactions with the membranes. 39,40 Others have focused on the simulation of biomembranes at both atomistic and CG resolutions. 27 Martini coarse grain (CG) force field is an extremely versatile model currently most adapted to biological systems. 42,43 It has rapidly gained popularity since its first use on simple lipidic systems. 42,44 The success of the model rests on its grounding on simple and intuitive physicochemical principles, which make it easily transferable to a wide range of systems and fields from biochemistry to materials sciences. 28 The model has recently been applied to complex systems such as mimics of realistic biological membranes with up to 62 different lipid types and proteins. 45,46 As a CG model, Martini reduces the resolution of the representation of a system of interest by discarding degrees of freedom (DoF). It is understood that the questions under investigation with such model must not strongly depend on those DoF. In a similar manner, as atomistic models neglect the electronic DoF, the Martini model averages atomic properties to chemical entities and neglects individual atoms. Its resolution is supra-atomic, with on average four heavy atoms (nonhydrogen) and associated hydrogens grouped together into superatoms or CG beads.
The simplicity of the Martini model relies on the definition of a limited set of building blocks: 19 different types of CG beads, covering chemical group properties from superpolar to very hydrophobic. 43 The transferability relies on the parametrization of those building blocks using thermodynamic data, in particular, the partitioning free energy of the CG beads between various media with different polarity or hydro-phobicity. The propensity of chemical groups (building blocks for larger biological molecules) to partition in different environments, in combination with more specific interactions, is at the heart of many biophysical processes. 47 In the Martini model the (nonbonded) interactions use simple Lennard− Jones potentials and are defined by what we call the "interaction matrix". 43 For molecules for which the mapping requires more than a single bead (lipids, ligands, amino acids and proteins, nucleotides and DNA, etc.), effective bonded potentials are built such that the geometrical features and the conformational flexibility of that particular molecule, as observed in simulation at the atomistic resolution or experiments, are reproduced. 48 Bead-type assignment (typography) is then determined and refined to reflect partitioning data and energy profiles if available. 48 Topologies for the Martini model are now available for a large variety of lipids, 43,49 amino acids, 50 sugars, 51 nucleotides, 52 fullerene, 53 and polymers, 54,55 and tools have been developed to build topologies for lipids, proteins, and DNA and construct systems where they are mixed. 49 These scripts are accessible on the Martini Web site (cgmartini.nl). More details about the Martini model can be found in our original papers 42,43 and reviews. 24,28,48 Examples of the Martini topology for a typical lipid, peptide, and a GPCR are shown in Figure 1.
2.1.2. ElNeDyn Approach for Proteins. A current limitation of the Martini CG force field is the lack of a systematic representation of the directionality of interactions between polar groups. This aspect of molecular interactions is important for small local dipole−dipole interactions such as hydrogen bonds. 56 H-bonding directionality is of particular importance for the protein backbone stability. 57−59 The secondary structures elements such as α-helices and β-sheets entirely depend on them. Implementations of local dipoles in Martini have been used to mimic the electrostatic shielding of water 60 and to improve the behavior of polar amino acid side chains, 61 but a satisfactory approach for the protein backbone has not yet emerged. Thus, the use of additional constraints is mandatory to maintain the secondary and in most cases the ternary structure of proteins. Bond and Sansom first introduced Mapping examples for a few biological molecules: 1,2-dipalmitoyl-sn-glycero-3phosphocholine, DPPC; water; a protein helical fragment in atomistic (AA) and CG resolutions. In the three cases, the AA resolution is depicted as ball and stick and the CG resolution by the underlying bond structure in black sticks with the beads or super atoms as transparent van der Waals spheres. In the CG helical fragment, the backbone beads are highlighted by red full spheres. The type of bead is indicated in gray on the side of the beads. Adapted from ref 48 with permission. Copyright 2013 Springer. (B) CG model of a GPCR, rhodopsin. The atomistic structure is shown as a cartoon with the seven helical segments colored magenta, the sheets on the extracellular side (top) of the protein in yellow, and the loops in gray. A transparent surface is used to indicate the van der Waals space occupied by the protein. The coarse grain model is shown using van der Waals spheres for all beads. The color code indicates the polarity of the beads as defined in the model. 50 The elastic network used in the ElNeDyn approach 65 is shown together with the trace of the protein (black sticks). The graph reports the root-mean-square-deviation (RMSD) of the rhodopsin relative to its starting structure using either the standard Martini model version 2.1 50 or the ElNeDyn approach. 65 the idea of using an elastic network (EN, a set of strings connecting the protein backbone interaction sites) to maintain the fold of a protein 62,63 in their pioneering work on using the Martini CG to study protein/lipid interplay. 64 Systematic parametrization of the Martini EN was subsequently performed in our group and coined ElNeDin for Elastic Network en Dinamicá(Spanish for in dynamics). 65 ElNeDyn (English version) is now a very standard protocol for proteins consisting of more than a few interacting helices. See Figure 1B for an illustration on a GPCR. Typically β-sheets are extremely deformable in standard Martini, while α-helices require strong bonded terms. 50 More sophisticated ENs have been developed to better capture the complexity of some Views from the extracellular, cytoplasmic, and membrane side are shown. This system is typically used to study the details of the protein/membrane interactions, such as membrane deformation and lipid binding sites. Lipids are shown in gray, white, and blue stick for the tails, glycerol backbone, and head, respectively. Water phase is omitted in the extracellular and cytoplasmic views and shown as blue spheres in the side view. (B) Sixteen receptors embedded in a membrane bilayer. Starting structure (t = 0 μs) with random orientations of the receptors is depicted together with the progression of the systems in time depicted by snapshots at 4, 12, and 20 μs*. (C) Sixty-four receptors embedded in a membrane bilayer after 100 μs* of selfassembly. (D) Two receptors embedded in a lipid bilayer. This system is commonly used in the determination of PMF using the umbrella sampling (US) approach with the weighted histogram analysis method (WHAM). In the US the distance between the proteins, d, and their relative orientations, ϕ1 and ϕ3, is controlled using harmonic potentials in order to probe defined interfaces. This is done by the virtual bond algorithm (VBA), which defines a distance and two dihedral angles, ϕ1 and ϕ3, or alternatively two angles. See the Methods and Systems section for details. Dimensions of the systems are indicated. Systems were constructed using rhodopsin as receptor. In A−C, these proteins are depicted by orange tubes positioned at the location of the 7 transmembrane helices and the helix 8 is shown in black. In A, the lipids are shown in gray. In B and C, periodic boundary conditions are shown in gray with the proteins in green. Protein:lipid ratios were 1/466, 1/328, and 1/100 in the monomeric (A), dimeric (D), and 16-(B) and 64- protein flexibility, e.g., a virus capsid, 66 membrane proteins, 67 and more recently GPCRs, μOR, 68 and CXCR4. 69 2.1.3. Note on the Martini CG Time Scale. In the parametrization of Martini molecules in general and for the EN in particular, it is important to match the time scales of the (atomistic) simulationson which the EN parameters are built uponwith the CG simulation. The time scales in CG models, their dynamics, are increased in most cases relative to atomistic models. Notably, the removal of atomistic frictions will smoothen the free energy landscape. A factor of 4 was derived from the increase of water diffusion and alkyl chains conformational dynamics. 42,43,70−72 Accordingly, a similar scaling factor was used to collect fluctuations in atomistic simulations and compare them to CG simulations to building an EN4 times longer AA than CG simulations are required. Using a different value for this factor will affect the ratio of conformational flexibility and dynamics of the CG protein vs the lipids chains and, in turn, might affect the heat transfer between them. This change might be most directly pictured in the loop regions of proteins. However, as an elastic object the protein might also react differently in response to pressure from a membrane bilayer or other external forces. These effects have not been systematically quantified and might be minor. In this review, time using the factor 4 are indicated by an asterisk (*).

Methods and Systems
In this section, we quickly describe the techniques that have been used, in conjunction with the Martini CG model, in the articles discussed in this review.
2.2.1. Molecular Dynamics Simulations. Molecular dynamics (MD) simulations in general but of biological systems in particular have tremendously evolved since their early start. 73,74 Developments of both computer technologies and algorithms efficiency have contributed to this progress. The Martini CG model has been originally developed for the GROMACS software, 75 a versatile user-friendly package dedicated to high-throughput MD simulations and analysis. Implementations are also available for NAMD, 76 GROMOS, 72 and Desmond. 77 A standard MD simulation is literally and simply the propagation of a system in time. This numerical performance is achieved by iteratively integrating the equations of motion Newton's second law. The particles of the system move according to forces acting on them at a given timedescribed by a force field and thus resulting from the interactions between them. The trajectory of a system constitutes the raw data, and only its analysis reveals the system's behavior. The quality and sophistication of the analysis of the trajectory are therefore determinant.
The MD simulation of a system is defined by its configuration and the thermodynamic state it is performed in (temperature, pressure, etc.), but it also depends on the starting conformation of the system as it will often only explore a restricted region of the phase space. CG models and Martini, in particular, are attempts to alleviate this limitation by reducing the number of DoF and allow the system to explore more systematically equilibrium conditions. This has been illustrated by keystone simulations such as the formation of complex lipid phase in early Martini studies 42,44 and more recently on a multitude of systems. 28 In the review, we use the Martini CGMD simulation/method to refer to the use of the Martini CG model in MD simulations, which includes the typical setups of such simulations. A CGMD simulation, however, is not a different technique than regular MD simulations.

Systems.
Most of the studies discussed in this review have used relatively simple systems, which are summarized in Figure 2. The simpler system consists in a receptor embedded in a membrane bilayer. The bilayer composition may range from a single lipid component to a more sophisticated mixture aiming at mimicking native membrane compositions. Details are given for each study when appropriate. A single-receptor system is typically used to characterize the interplay between specific lipid molecules, the bilayer matrix, and the receptor. Lipid binding sites, membrane adaptation to the presence of a receptor, lipid, and receptor lateral mobilities may be determined. This simple system is often used as a unit cell to build more complex or crowded ones.
Khelashvili and co-workers used a system notably different from the common lipid bilayer. They built a lipidic cubic phase (LCP) with monoacylglycerol molecules and characterized its geometrical properties for three different chain lengths, composition, and temperature ( Figure 2F). 78,79 2.2.3. GPCR Self-Assembly Simulations. Self-assembly simulations have been part of the Martini CGMD studies of GPCRs from the start 80 and have been regularly used since then. 81−86 This type of simulation aims at studying a system in which multiple receptors coexist, as in membranes. They consist in embedding a large number of receptors (2−144) in a preformed lipid bilayer and regularly spaced to maximize their dispersion in the membrane ( Figure 2B), often built from the repeat of a unit cell containing a single receptor. From that dispersed configuration one may follow the assembly of the receptors. Unfortunately, on the time scale accessible to the Martini CGMD simulations, it is mainly the association of the receptors that is observed. Only limited events of dissociation have been reported. The supramolecular organization of the receptors described by Martini CGMD simulations on time scales up to 100 μs is most likely not representative of biological time scales and their complex dynamics in physiological conditions. A more realistic lipid composition would also be beneficial. These self-assembly simulations are typically out-of-equilibrium simulations. However, they provide extremely valuable information on the receptor behavior in a lipid bilayer, e.g., their propensity to oligomerize, the interfaces that are accessible, and the kinetics of the contacts. The comparison of simulations performed with different conditions, e.g., lipid composition (thickness, viscosity) and receptor, allows one to characterize the oligomerization pattern pertaining to lipid/receptor interplay. 46,80−86 In the search for relevant protein−protein contacts, Wassenaar et al. developed a protocol, "docking assay for transmembrane components" (DAFT), 87 to predict effectively the most accessible interfaces of transmembrane (TM) helices and proteins. DAFT is a clever way to setup systematic selfassembly simulations of proteins/peptides. It consists of automatically setting up a large number of simulations with the molecules randomly oriented and embedded in a lipid bilayer. The protocol was tested with the Martini force field and uses MD simulations with the associated setup. It can extend to systems with up to nine different protein monomers assuring the best-starting structures with optimized periodic boundary conditions to sampling partners' interfaces. The protocol was validated on simple transmembrane helices and the GPCR rhodopsin. The results showed that for a simple system such as a transmembrane helix (glycophorin A) with a unique dimeric interface the method is able to reliably predict the complex and the effect of debilitating mutant. It was also able to predict the trimeric arrangement of the GCN4-derived peptide MS1, a transmembrane helix. For more complex systems such as GPCR oligomerization where multiple interfaces may form, DAFT, as regular self-assembly simulations, has difficulties to generate interfaces with populations reflecting their relative thermodynamic stability. The protocol relies on a large number of short (∼1 μs) simulations that can only capture events accessible on that time scale.
The lack or limited amount of binding/unbinding events in self-assembly simulations in general but of GPCR in particular, by definition, precludes the access to thermodynamic quantities such as the relative binding free energies of specific interfaces. 81 Another point to consider when interpreting the populations of interfaces obtained from self-assembly simulations is a possibly strong kinetic contribution to the receptor approach. This point was recently clearly demonstrated by Provasi et al., 85 see below.
Nevertheless, self-assembly simulations have revealed important lipid/protein interplays and will continue to be of great use. The combination of self-assembly simulations of multiple receptors revealing cooperative features and supramolecular organizations [80][81][82]85,86 and DAFT generating more systematic information on accessible interfaces will be of interest to GPCRs 41 but also for other transmembrane proteins. 88, 89 2.2.4. Biased Simulations: US/WHAM and Metadynamics. An area of great interest for GPCRs is the search for structural keys to the supramolecular organization of rhodopsin or the oligomerization pattern of GPCRs. The interfaces involved in contacts between receptors and their associated strength (relating to their binding free energy) are central in that matter. In principle, the relative strength of protein interfaces can be extracted from the radial distribution function (RDF) or the potential of mean force (PMF; RDF(r) = exp[-PMF(r)/k B T]) computed from an unbiased simulation. However, the simulations must have sampled the conformational space accessible to the system at equilibrium. In the present case, that means binding and unbinding events of all possible protein interfaces. This degree of convergence is not yet accessible within the current time scale of Martini CGMD simulations. The self-assembly described above provides great information about the interfaces but not their relative strength.
Alternatively, one can determine the potential of mean forces (PMF) between two receptors using the umbrella sampling (US) technique. 90 It discretizes the reaction coordinate, the distance between receptors, into bins that are explored independently using a biasing umbrella potential to maintain the distance close to the value at each bin (umbrella window) and assuring overlap between consecutive bins. The deviation from the ideal distribution of the distance reflects the (un)stability of the system in that particular window, and the knowledge of the biasing potential allows one to generate unbiased distributions corresponding to each umbrella window. These are combined afterward using the weighted histogram analysis method (WHAM). 91−93 Filizola and co-workers and Periole and co-workers used this approach on opioid receptors and on rhodopsin, respectively (see below), with similar strategies to explore particular interfaces. The relative orientations of the proteins were "restrained" using additional harmonic potentials but with slightly different restraining potentials on the protein orientations. 37 The dihedral angle restraints used by Periole and co-workers (part of the virtual bond analysis, VBA, 94 Figure  2D) present the advantage to cover the full range of protein orientations and avoid redundancy ( Figure 4). Periole and coworkers used a six-dimensional WHAM (although only 3 were actively used) to be able to correct for these additional biases in the final PMF. 81 Filizola and co-workers added a layer of sophistication to the characterization of the GPCR interfaces. They used welltempered metadynamics 95 to enhance sampling when generating starting configurations and free energy maps of the relative orientation of the receptors and also combined with the US technique at each umbrella simulation. Well-tempered metadynamics adds Gaussian biases onto defined collective variables, here angles maintaining the receptor orientations to remove energy barriers along them. The biases might be collected and used to reconstruct the free energy profile along the collective variables. Filizola and co-workers also corrected for the use of metadynamics when combined with US. 96

MARTINI CGMD SIMULATION STUDIES ON GPCRS
In this section, we review the studies that have used the Martini CGMD simulation approach and added techniques described above to study GPCRs. These studies have mostly been performed in a handful of laboratories keeping the technique quite close to the original one so that it is easier to relate them. Toward the end of this section we review a few more recent and isolated studies that further illustrate the exciting potential of the CGMD approach for GPCRs.

Rhodopsin Supramolecular Organization
The Martini CGMD simulation venture into the exploration of GPCR signaling was pioneered by the study of the supramolecular organization of rhodopsin, 80 a light-sensitive receptor involved in the visual phototransduction. This study was motivated by the report of an intriguing highly ordered organization of the photoreceptor in rod outer segments (ROS) disc membrane from AFM images. 97 Furthermore, our colleagues from Sakmar's laboratory at Rockefeller University/ NYC and Brown's laboratory from the University of Arizona at Tucson were able to strictly correlate the spatial distribution of rhodopsins in model membranes with their photoactivation. 98 Of particular interest was the similar response of rhodopsin's spatial organization (FRET efficiency) and photoactivation to the change of membrane thickness. The results strongly suggested that the hydrophobic mismatch between the receptors and the lipid bilayer governs their degree of association and thereby their activation.
A first set of Martini CGMD simulations 80 described the selfassembly of 16 rhodopsins embedded in a lipid bilayer with a range of thicknesses matching the experimental data. 98 The receptors were observed to spontaneously self-assemble and doing so to a degree depending on the membrane thickness, in perfect agreement with the experimental data from Botelho et al. A maximal dispersion of the receptor was observed for an intermediate lipid thickness matching the hydrophobic thickness of rhodopsin, while thinner and thicker bilayer increased the receptor's propensity to form contacts. The correlation between the hydrophobic mismatch and the assembly propensity of the receptor was not a complete surprise to either experimental or computational approaches as hydrophobic mismatch had been pointed out as a driving force for membrane protein association in general 99 and for rhodopsin in particular. 98 The Martini CGMD simulations matched quite well the FRET signal reported earlier.
An exiting novel finding of the CGMD study was the demonstration of nonhomogeneous membrane deformation around the protein ( Figure 3A and 3B), i.e., the membrane bilayer adapted to the protein/membrane interface variably at different regions of the protein surface. In addition, the regions of the protein surface where the membrane deformed varied with the thickness of the membrane; in other words, the protein/membrane interface heterogeneity varied with the hydrophobic thickness of the membrane. Of particular note was that the regions on the protein surface where the membrane deformed the most in response to hydrophobic mismatch strongly correlated with the location where the protein was seen making protein−protein contacts upon self-assembly. At that point the details of the protein−protein contacts were not analyzed, as the statistical relevance was considered not sufficient to be conclusive.
On the basis of the assumption that protein/membrane hydrophobic mismatch drives integral protein assembly by reducing the membrane deformation, 99 these results suggested that protein contacts would form favorably at specific locations of the protein surface and that these locations might vary with the membrane thickness. While the former interpretation has been confirmed by several approaches (see below), the latter has not yet been observed. Inhomogeneous membrane deformation around rhodopsin was confirmed by atomistic resolution MD simulations ( Figure  3C). 100 On the basis of the analysis of these higher resolution simulations the authors described the free energy associated with the presence of an hydrophobic mismatch between the surface of the protein and the lipid bilayer thickness as the sum of two main components: the membrane deformation to match the hydrophobic surface of the protein to the best of its abilities and the residual hydrophobic mismatch (RHM) occurring when the protein remains exposed to an unfavorable environment. The residual hydrophobic mismatch was found to be an important contribution to the system free energy. See below for more details on the studies from Weinstein, Khelashvili, and coworkers on the membrane/protein interplay. 39,40 In a follow-up study, Periole et al. 81 probed more specifically the relative strength of different interfaces in rhodopsin dimers. This work was intended at building a model of rhodopsin organization in its native environment, the ROS disc membrane, depicted by AFM images 97 as a highly ordered set of rhodopsin row-of-dimer. The details of the interfaces could not be resolved from the AFM data, and the nature of the dimer interface remained highly debated. A row-of-dimers model was built based on the images that showed a symmetric TM4/5 rhodopsin dimer. 101 This apparent static arrangement of rhodopsins was quite controversial as apparently contradicting earlier biophysical experiments. 102 However, it was extremely exciting to have access to structural data relating to the supramolecular organization of rhodopsin in native conditions and build models to reveal their interactions at a pseudo atomistic resolution.
Periole et al. 81 probed rhodopsin interfaces using two complementary approaches: self-assembly simulations, similar to the earlier study, 80 and calculation of the potential of mean force (PMF) as a function of the receptor separation ( Figure  4). While the first approach allows observing the formation and characterization of the interfaces most accessible to the receptors as they assemble on the time scale simulated, the second approach allows a precise quantification of the relative strength of the interfaces. Ideally, one would use the first type of simulations (multiple replicas) from which potential interfaces may be identified, and subsequently, the second method may be used to compare them.
The self-assembly simulations were carried out using conditions where the receptor interactions would be the least affected by the membrane deformation due to the presence of a hydrophobic mismatch. 80 The simulations confirmed the preferential linear arrangement of rhodopsins when embedded in a membrane bilayer and revealed a few preferential interfaces ( Figure 2B and 2C). The interfaces involving TM1 and H8 (simultaneously) or TM5 combined symmetrically and asymmetrically were predominantly formed. The interfaces involving the other sides of rhodopsin, around TM4 and TM6, did not form in a significant amount. The PMFs rationalized these observations. The limited formation of interfaces centered on TM4 and TM6 was due to the presence of an energy barrier to their formation ( Figure 4E). This energy barrier resulted from the trapping of lipids at the interface and the stabilization of a metastable state where lipids lubricate the interface. The interfaces forming in the self-assembly simulations (TM1/H8 and TM5) did not show an energy barrier to their formation ( Figure 4E). The PMFs led to the definition of two types of interfaces: some weak with an energy barrier to their spontaneous formation and others strong with a deep minimum at the interface and no energy barrier to complexation ( Figure  4E).
Interestingly, the PMFs demonstrated a striking stability of the symmetric TM1/H8 dimer interface compared to the other interfaces probed, an interface long discarded for the small size of protein burial associated with it. 103 It was however not a complete surprise as this interface had previously been observed in two-dimensional 104 (2D) and three-dimensional 105 (3D) electron microscopy (EM) and X-ray crystallography of opsin, rhodopsin, and metarhodopsin I and II. 106−108 This data motivated an earlier study in which we were able to show that the symmetric TM1/H8 interface existed in the native membrane based on a combination of chemical cross-linking experiments, partial proteolysis, and high-resolution liquid chromatography-mass spectroscopy. 109 Many other GPCRs have also been crystallized with the TM1/H8 interface. 110−112 On the basis of the TM1/H8 rhodopsin dimer interface it was possible to construct a row-of-dimer model that fulfilled the structural feature extracted from the AFM images. In addition to the TM1/H8 dimer interface, this model also complied with the existence of lubricated interfaces, centered on TM4 and TM6, that would stabilize interfaces between two dimers in a row. We further discuss this organization later on.
On a biophysical perspective on protein/membrane interplay, the stability of the TM1/H8 interface in the rhodopsin dimer is of primary importance, because it has a smaller protein burial than the others. It is actually quite systematic that within the interfaces probed the strong ones have smaller burial than the weak ones ( Figure 4 A and 4E and ref 81). This observation strongly challenges the use of buried accessible surface area as a predictor of the strength of membrane-embedded protein− protein interfaces. Other forces must apply. We found a few amino acid interactions forming networks at the interface. One study recently suggested that the residual hydrophobic mismatch might be a significant contribution to GPCRs interface stabilization (see section 3.3 for more details).

Rationalization of Opioid Receptors Oligomerization
Filizola and co-workers, inspired by the Martini CGMD simulation technique 80 described above and its combination with metadynamics 95 to enhance sampling, published a nice series of studies probing GPCR interfaces. They emphasized on δ-, μ-, and κ-opioid receptors (OR). These studies were a logical continuation of their previous interest in the determining GPCR interfaces. 113 They started by characterizing the binding profile of the δopioid receptor (δOR) at the TM4 (also referred to as TM4/3) interface calculating the first PMF of a GPCR interface. 114 They subsequently compared it to the PMF of the TM4/5 interface, allowing them to rationalize cross-link experiments. 115 The receptors were embedded in a POPC lipid bilayer with 10% cholesterol, and for the model of δOR they used a similar approach to that developed for rhodopsin. 65,80 Only the EN defining the framework of the protein was modified to give more conformational flexibility to the loop regions. 114 Although the sampling might have been slightly limited due to short simulation times accessible at the time (as shown by the roughness of the PMF curves, Figure 4C), these PMFs provided a fascinating view of the interfaces profile and their relative strength. The simulations gave the TM4 (or TM4/3) interface of δOR as favorite compared to the TM4/5 interface. However, both interfaces had no energy barrier to binding ( Figure 4C), and the free energy stabilization in the dimer configuration suggested a half-life time on the order of seconds, 0.2 and 4.4 s for TM4/5 and TM4 (or TM4/3), respectively.
In a later study, Johnston and Filizola extended their work on OR using PMFs to look at interfaces of μOR and κOR as found in the crystallographic structures. 116 These include the TM1/ H8 for both receptors and the TM5/6 for μOR. The TM1/H8 interface uses contacts of TM1 on the side of TM2, similar to the case of rhodopsin. The results show that the symmetric interface involving TM5/6 is more stable than the one involving TM1/H8 and by a significant amount when compared to the other interfaces and other receptors ( Figure  4C). The TM1/H8 interface of μOR has a similar strength as the TM1/H8 interface of other GPCRs (β 1 AR and β 2 AR, described below and using the TM7 side of TM1 for contacts; Figure 4B). Striking is the increase of the TM1/H8 interface stability for κOR by a factor close to 5 in terms of depth of the free energy well compared to μOR ( Figure 4C). The authors attributed this increase to a different conformation of the helices 1, 2, and 8 in the two receptors. TM1 is pointing away from the helical bundle in κOR, while it is close to it in μOR as in most receptors. This observation is a clear illustration that the local structure of the receptor determines the interaction between them, adding to the importance of the details of an interface in determining its strength.
In their most recent work on OR, Provasi et al. 85 used the self-assembly approach to looking at the interfaces preferentially formed by three receptors from the opioid subfamily they have studied previously using PMFs: δOR, μOR, and κOR. They investigated homomeric complexes for the three receptors and heteromeric between the pairs δOR/μOR and δOR/κOR. As noted by the authors, although this approach does not allow the calculation of free energies and thereby prevents the comparison of the relative strength of the interfaces formed (based on their populations) as with PMFs, 81 it reveals the interfaces accessible on the time scale simulated. In all five systems the opioid receptors formed filiform (elongated) structures using the small sides of the receptors: centered either on TM1/2 or TM5 alone and combined with TM4 or TM6 ( Figure 4B). Significant populations of symmetric interfaces were observed only in the cases of TM1/2/H8 and TM5. The frequency of appearance of the interfaces depends on the receptors simulated, but their confidence intervals were broad, and the differences might not be significant. This lack of significance of the different populations is demonstrated by the fact that they do not reflect the results from the PMFs determined previously. 116 For instance, the symmetric TM5/6 interface for μOR is given to be more stable than the TM1/2/H8 one ( Figure 4C), but it is not significantly more populated. According to the PMFs, an even larger population would be expected for the TM1/2/H8 interface of kOR. Also of note is the absence of the TM4 (or TM4/3) interface for the δOR, while the TM4/5 is present. The PMFs indicated that they were not so dissimilar ( Figure 4C). 114,115 Finally, the TM3 and TM7 were not involved in interfaces.
An interesting aspect of this last study from Provasi et al. 85 is the attempt to correlate the varying dynamics of the lipids around the receptors with the interface-specific dimerization rates (k on ) in all systems. Their analysis suggests that the heterogeneity of the lipid dynamics around the OR defines different zones of viscosity that in turn modulate the kinetics of association of the receptors. They notably observed regions on the surface of the receptors associated with "jammed" lipids having extended persistent times, thereby preventing the receptors from approaching each other. The authors argue that this behavior may rationalize the absence of the TM5/6 symmetric interface in the case of μOR that we noted above. This argument would actually be more general and reduce both TM4/5 and TM5/6 interfaces, while TM1/2/H8 would be favored kinetically. "Jammed" lipids were also observed in the case of rhodopsin. 81 In that case, they stabilized a metastable state of the TM4/4 interface in which the interface was lubricated by a layer of lipids. 81 These kinetically trapped lipids resulted in an energy barrier in the PMFs ( Figure 4E) and prevented the spontaneous association of rhodopsins.
In further studies it would be of great interest to investigate the correlation between the heterogeneity of the lipid diffusion (and the associated membrane viscosity) around the receptors reported by Provasi et al., which suggests stronger lipid− protein interaction, and the membrane deformation resulting from the attempt of the membrane to compensate for the hydrophobic mismatch between the lipid bilayer thickness and the hydrophobicity of the protein surface. Periole et al. reported varying heterogeneous membrane deformation around rhodopsin dependent on the membrane bilayer thickness 80 and noted the correlation between the locations of these deformations with the protein contacts made during self-assembly simulation. Such deformations were also reported for β 1 AR and β 2 AR and their differences used to rationalize the oligomerization patterns of these receptors based on the existence of a residual hydrophobic mismatch, 82 as discussed in the following section.

GPCRs Stability in LCP To Rationalize the "in Meso" Crystallization
Khelashvili and co-workers used the Martini CGMD approach 80 to study the details of the interactions between GPCRs and their environment in order to explain their behavior in various lipid phases 78,79 and identify forces driving their assembly. 82 Before engaging in actual Martini CGMD simulations, Mondal et al. developed a method to characterize and quantify the energy associated with membrane deformations around multihelical proteins using MD simulations at atomistic resolution. 100 The authors confirmed the existence of membrane thickness-dependent heterogeneous deformation around the rhodopsin as reported earlier using Martini CGMD simulations ( Figure 3). 80 They further used these atomistic simulations to calibrate and validate a threedimensional continuum (CT) model referred to as 3D-CTMD for its combination with MD. This model integrates two key contributions to the energy cost of hydrophobic mismatch between the protein and the membrane bilayer: 117 the membrane deformation aiming at minimizing the mismatch and the residual hydrophobic mismatch (RHM)the contribution from an incomplete hydrophobic matching resulting in energetically unfavorable exposure of the protein at the protein−lipid interface. The energy term associated with RHM is derived from the exposed surface area (extracted from the MD simulation) and transfer energies between hydrophobic and polar environments of protein residues. The authors identified the free energy associated with RHM as a contribution at least as important as the one from the membrane deformation itself; both contribute to a few k B Ts. Furthermore, they showed that by averaging out the deformations, the consideration of the radial heterogeneity of the membrane deformation leads to an underestimation of the actual contribution. They also show that the RHM contribution varies with the helical segment of rhodopsin and suggest that their relative contributions on TM1, TM2, TM4, and TM5 in different bilayer thicknesses would explain the contacts between receptors observed in self-assembly simulations. 80 Taking advantage of the increased sampling and system size offered by the Martini CGMD approach and its proven ability to reproduce the system's behavior observed at atomistic resolution, 100 Khelashvili et al. 78 applied the 3D-CTMD approach to rationalize the mechanism of "in meso" crystallization, a method that recently became popular for GPCRs. The challenge was to be able to simulate a receptor in a lipidic cubic phase (LCP, Figure 2F) and lamellar phase long enough to extract statistically significant quantification of the hydrophobic/hydrophilic exposure of residues at the protein/ lipid interfacedata mandatory for the 3D-CTMD approach. 100 Both systems could be simulated to convergence using rhodopsin and monoolein as GPCR and lipid models, respectively. The results revealed that the highly curved LCP was actually able to accommodate better to the rhodopsin radially inhomogeneous hydrophobic surface than the lamellar phase, and thus, LCP reduced the residual hydrophobic mismatch energetic cost. Rhodopsin is thus well accommodated in an LCP environment. The authors further argued that rhodopsin's high cost for mobility in LCP will not be counterbalanced by the low RHM and thus will favor a monomeric state of the protein. In contrast, in the lamellar phase the lower cost for mobility of the receptor and the higher RHM would drive the system toward receptor oligomerization. Although the mechanism of transition of the protein from the LCP to the lamellar phase is not yet clear, the oligomerization of the receptor in lamellar stacks may lead to bulk crystals as well as other steps involved in the crystallization process, e.g., addition of precipitant. 78 Of particular interest in this study is the identification of TM1 and TM5 as the regions on the rhodopsin surface with higher residual exposure in the lamellar phase. Assuming that hydrophobic mismatch is a driving force of protein oligomerization in the membrane, this observation suggests that these two regions would be primarily involved in protein contact in the lamellar phase and thus in the crystals. This hypothesis was supported by the analysis of GPCR crystal structures issued from the in meso technology. Out of 14 protein−crystal contacts (from 12 structure files 78 ) 10 use TM1 and/or TM5 in "canonical" interfaces in which the receptors are parallel. It remains, however, to quantify to what extent the different receptors experience similar residual exposure on TM1 and TM5 as rhodopsin does in monoolein. Additional support for this hypothesis comes from the self-assembly simulations of rhodopsin in regular lipid bilayer, in which TM1 and TM5 were found involved in the most prominent protein contacts. 81 This might also be the case for the δ-, μ-, and κ-opioid receptors as discussed above and reported by Provasi et al. 85 82 In a subsequent study, Khelashvili and co-workers further investigated the reasons for the success of the in meso crystallization method. In this very well conducted work Johner et al. 79 first characterized the geometrical features of the Pn3m LCPs, such as mean curvature, surface area, and Gaussian curvature of the midplane and water−lipid interface for different values of the lattice constant (dimension of the lattice). They performed this characterization for three monoacylglycerol (MAG) molecules with different tail lengths: 14, 16, and 18 carbon atoms, 7.7, 7.9, and 9.9 MAG, respectively. Besides validating their CG model of the lipid, this analysis revealed the inhomogeneity of the lipid bilayer properties in the LCP. Relevant here is the variation of the membrane thickness reaching its thinnest level at the saddle point of the LCP, where the proteins are thought to reside. The authors went on by investigating the behavior of the adenosine A2A receptor, A 2A R, into a LCP built with 9.9 MAG. They specifically looked at how a different lattice constant could accommodate A 2A R and A 2A R engineered with a thermostabilized apocytochrome b562 from E. coli (M7W, TM102I, K106L, referred to as BRIL 118 Figure  5G). The contact made was similar to the one found in the crystal ( Figure 5G). 118 In the context of the in meso crystallization method, Johner et al. interpreted these results as (i) the addition of precipitant (modifying the lattice constant) would destabilize the GPCR in the LCP and promote its migration to an adjacent lamellar phase and (ii) the formation of contact between closely spaced copies of the proteins in LCP (as seen for A 2A R-BRIL) would possibly promote the stacking of the protein before their migration to the stacked lamellar phase. The results also rationalize the dependence of the in meso success on the engineering (size, shape, and orientation) of the GPCR.
3.4. β 1 and β 2 Adrenergic Receptors Oligomerization Pattern β 1 and β 2 adrenergic receptors (AR) are of particular interest because these two receptors are highly similar (67% sequence identity) but do experience distinct association patterns in the membrane. 119,120 Johnston et al. 96 compared their interfaces centered on TM4/3 (similar to the one named TM4 in their previous study on δOR, see above and Figure 4B) and on TM1/H8. The PMFs were this time obtained using Martini CGMD simulations coupled to well-tempered metadynamics 95 to augment the sampling (see Methods and Systems section). The comparison of the PMFs of both receptors and for both interfaces demonstrated the two receptors to behave identically within the error bars, at odds with FRAP experiments giving a more stable dimer for β 2 AR. 119 The PMFs also showed that for both receptors the interface TM1/H8 is stronger than the TM4/3. The authors interpreted the results as due to the fact that the receptors might be diffusing as dimers associated through the TM1/H8 interface and dimers of dimers would weakly interact using other interfaces, e.g., TM4/3, thereby reflecting the FRAP data. 119 The similarity of the PMFs for the two receptors also contrasts with a later work on these receptors using the same CGMD technique and to which the authors contributed. In that study significant differences were observed between β 1 AR and β 2 AR interactions with the bilayer environment, consonant with experimental observations of oligomerization patterns of β 1 AR and β 2 AR 82 (see the next paragraph for more details). A possible explanation for this discrepancy is the calculation of the PMFs for the TM1/H8 interface of β 1 AR and β 2 AR in which the contacts made by TM1 face the side of TM7 instead of TM2 (TM1/H8* in Figure 4B).
Mondal et al. applied their Martini CGMD/3D-CTMD combination to the case of β 1 AR and β 2 AR. 82 In the case of β 2 AR the authors show that the association between the receptors reduces the energy penalty from residual hydrophobic mismatch (RHM, Figure 6A). They further show that the magnitude of the RHM reduction is a function of the part of the protein involved in the contact, with TM1, TM4, and TM5 most affected ( Figure 6C). The analysis actually points to only a few pertinent residues. The correspondence of these regions with the most observed contacts formed in a self-assembly simulation suggested once more that the RHM is a major driving force toward the determination of protein contact upon oligomerization. The authors then made the comparison with the RHM found for β 1 AR, mainly pointing to TM1 as a hot spot while strongly reduced on the TM4/5 interface compared to β 2 AR ( Figure 6B). This difference provides a potential explanation for the oligomerization patterns of the two receptors. 119,120 β 1 AR would form mainly dimers at the TM1 interface, and β 2 AR higher order oligomers would form also involving TM4 and TM5 ( Figure 6). 82 A detailed analysis pointed to only a few residues to explain the difference in RHM.
In further studies it would be interesting to confirm this hypothesis by looking at the contact zones of β 1 AR in a selfassembly simulation. It would also be interesting to performed self-assembly simulations of β 2 AR mutants lacking the residues stabilizing the interfaces. The relevance of the few residues pointed out by this study could also be tested by mutagenesis experiments.
Ghosh et al. performed a self-assembly simulation of 16 copies of β 2 AR in a DSPC lipid bilayer. 84 Analysis of the simulation showed an overall similar behavior as previously described for other receptors. [80][81][82]85 Notably, the receptors form a string-like cluster (slightly more branched than for other receptors, but it could result from a longer simulation time), and the contact between receptors involves TM1, H8, TM5, and TM6. There is one significant difference in the simulations of Ghosh et al. They did not use any restraints to maintain the ternary structure of the helical bundle leading to a rmsd of β 2 AR from 5 to 8 Å. These deformations are huge, and although the helical bundle might still look like the original receptor, it has certainly adopted a structure with significant differences. One needs to keep in mind that the active state of β 2 AR is only ∼2.9 Å away from its inactive structure. (This value was obtained by comparing the active and inactive states of β 2 AR (restricted to residues 32−175, 179−230, and 265− 341) as found in 3D4S and 3SN6 PDB entries. A similar value, 2.85 Å, has been reported for rhodopsin. 108 ) Figures 9 and S8 in the manuscript by Ghosh et al. illustrate these structural differences. The figures depict a β 2 AR dimer formed during a simulation and its comparison with the dimer found in a crystal structure (PDB ID 3PXO 108 ). In the CG dimer, one of the interacting monomers has H8 perpendicular to the membrane plane. The authors interpreted the stabilization by ∼25 kcal/ mol of the TM1/H8 interface obtained in the Martini CGMD simulation compared to the experimental model by the requirement for "subtle rearrangement of the TM helices in order to form proper oligomeric assemblies, which was achieved through the CGMD simulation". These structural changes of H8 orientation are likely to be too drastic to be realistic.

Cholesterol Involvement in GPCRs Assembly
Sengupta and co-workers used the Martini CGMD simulation approach to study the interaction of cholesterol (chol) with GPCRs and rationalize its effect on their assembly. These studies follow the extended work of Chattopadhyay on the subject. Chol, a significant component of biological membranes, has been shown to play a critical role on membrane proteins function 121 and GPCRs supramolecular organization and function, 122−124 but its mechanism of action is still a matter of debate. Chol would act directly on GPCRs through specific interaction with the receptors 125−128 or indirectly through modification of the membrane bilayer mechanical properties. 123,129−131 In a first study, Sengupta and Chattopadhyay 132 characterized the interaction pattern of chol with the serotonin 1A receptor (5-HT 1A R). They ran Martini CGMD simulations of a single receptor embedded into a 1-palmitoyl-2-oleoyl-snglycero-3-phosphocholine (POPC) membrane bilayer containing 0%, 9%, 30%, or 50% chol. Overall, chol molecules interact with, or explore, most of the receptor surface during the simulations but spend the most time at, or bind to, specific locations. Notably, one of these locations corresponds to one of the three highly conserved chol recognition amino acid consensus (CRAC) motifs. 127 It is located on TM5. The other two CRAC motifs were not observed arguably due to the exposure of these sites to the aqueous phase. Chol is found to interact with 5-HT 1A R on time scales from nanoseconds to microseconds in the case of preferred occupancy sites. This pattern of interactions is in line with the concept of "nonannular" binding of chol to 5-HT 1A R. 125 Nonannular binding refers to lipids that do not frequently exchange with other lipids in the bulk membrane.
In a subsequent work, Prasanna et al. 83 studied the interaction of chol with β 2 AR, a receptor also known to have interactions with chol. 110,111 Here, the authors used selfassembly Martini CGMD simulation of a pair of receptors embedded in a POPC lipid bilayer with 0%, 9%, 30%, or 50% of chol. Analysis of these simulations led the authors to conclude that the presence and increase of the chol content in the membrane systematically affects the interface the receptors use to assemble. In the absence of cholesterol, β 2 ARs almost exclusively interact in a symmetric manner though TM4/5. With 9% and 30% chol, β 2 AR assembles less exclusively using TM4/5 as TM1/2 gets involved. At 50% chol, β 2 AR switches to an almost exclusive mode of assembly using TM1/2 in a symmetric interface. The reported effect of chol on the interface of assembly is significant, and the statistics are rather convincing. To rationalize this effect of chol, the authors analyzed the chol occupancy at the seven helices of the receptor. On the basis of their analysis the authors propose that the presence of chol in the membrane prevents the use of TM4/5 as an interface because it occupies an interaction site located at the center of TM4. The occupancy site corresponds to the position of a cholesterol molecule found in a crystal structure. 110,111 An interesting aspect of Prasanna et al.'s results is the apparent contrast of the proposed mechanism of cholesterol stabilization of β 2 AR dimer at the TM1/H8 interface by the simulations with the one that could be suggested by inspection of the β 2 AR structures. 110,111 From the crystal structure one could speculate that chol stabilizes the TM1/H8 interface by its location at the interface itself and thus acts as glue ( Figure 7A). This behavior has been suggested for other systems. 133 In contrast, the simulations point at an alternative mechanism by which chol would favor the TM1/H8 interface by destabilizing the TM4/5 interface. The latter is predominant in the absence of chol ( Figure 7B). The two mechanisms are, however, not mutually exclusive. From the data presented, a stabilization of the TM1 interface by direct interaction with chol cannot be excluded. In fact, the occupancy of TM1 and TM2 at 50% chol increases significantly from the unbound to the bound receptor situation (Figure 4 in ref 83). Other recent studies found chol at the TM1/H8 interface of GPCR dimers. 46,112 This effect of cholesterol on GPCRs assembly is of great biological and biophysical importance, and it might well be the case that chol modulates the interfaces of β 2 AR by either or both of the mechanisms described above (Figure 7). It is informative to compare with the study from Mondal et al. 82 based on a similar Martini CGMD approach and described above. First, there is an apparent agreement. Specifically, both studies give a mixture of TM1/2 and TM4/5 as the interface of β 2 AR with 9% or 10% cholesterol for Prasanna et al. 83  There are, however, a few points in the study of Prasanna et al. 83 that make the interpretation of the results not straightforward. The lack of palmitoyl chains attached to the receptor in Prasanna et al. might have an impact on the interaction of cholesterol with the receptor and of the receptors together. 110,111 The pattern of increase of chol occupancy described by the authors is not straightforward in the graph. Most differences are within the error bars. Notable is the mention of an increase of chol contacts with TM4 between 30% and 50% of chol. Other helices experience similar changes. Moreover, the main difference in chol contacts is observed between 9% and 30% chol, while the main difference of interacting helices is observed between 30% and 50% chol (Figures 3 and 4 in ref 83).
Of particular concern in the studies carried out by Sengupta and Chattopadhyay is the lack of use of a protocol to maintain the ternary structure of the GPCR. Most proteins deform significantly as the Martini CG force field is not capable of keeping the secondary and ternary structure without the use of an elastic network. We illustrate this deformation on the example of rhodopsin in Figure 1, which is in line with the work of Horn et al. 134 Such deformation was also reported for another study of β 2 AR using the regular Martini force field without ElNeDyn (described above). 84 Typically consecutive helices reorient under the pressure of the environment, leading to atomistic root-mean-square-deviations (rmsd) reaching ∼6− 7 Å. Such deformations should be considered as significant at the protein level. As a comparison, the rmsd between the inactive and active receptors, which include an important helical motion, is only ∼2.9 Å (see section 3.4 for more details). Visual inspection of the graphical representations of β 2 AR in Prasanna et al. (Figure 2 in ref 83) suggests significant helix reorientations. The deformations of the helical bundle might well explain differences of interfaces observed with different amounts of chol in the membrane and may affect interactions with cholesterol as shown for rhodopsin. 134 Despite these issues, simulating the effect of cholesterol on GPCRs assembly is of great biological and biophysical importance and led Sengupta and Chattopadhyay to propose a new paradigm 135 in which the cholesterol interactions with GPCRs would be described as a set of hot spots instead of actual binding sites in an "undulating energy landscape characterized by shallow minima with small energy minima". These "high occupancy sites" instead of strong interactions would contribute to moderating GPCRs cross-talk and drug efficiency.
During the reviewing process of this review, Prasanna, Sengupta, and Chattopadhyay published their latest work on the effect of cholesterol on GPCRs oligomerization. 136 In this new study, Prasanna et al. follow up on their earlier work 132 on 5-HT 1A R discussed above. They applied their receptor selfassembly strategy in POPC model bilayers with four chol concentrations (0%, 9%, 30%, and 50%) described above in the case of β 2 AR. The two papers are technically very similar, and the overall conclusions for 5-HT A1 R are very close to those obtained for β 2 AR, namely, the presence of chol affects the mode of interaction of the receptors (their interface), and this effect depends on the concentration of chol. A few differences are also noticeable. First, as described by the authors, in the case of 5-HT A1 R, the presence of chol seems not only to destabilize the dimer interface (TM1) observed in pure POPC to favor different interfaces, TM4/5 and TM5/6, but also to weaken the interactions or their specificity quite significantly. This point is clearly illustrated by the absence of a predominant interface in the presence of chol in the case of 5-HT 1A R but not for β 2 AR (compare Figure 3a . The authors use this lack of preferred interface for 5-HT 1A R in the presence of chol to argue for a functional modulation of the receptor interface plasticity by chol. Second, and unfortunately not discussed by Prasanna et al., the prominent interface observed in the 5-HT A1 R dimer in the absence of cholesterol, TM1/2, and destabilized by chol is the same interface that is absent in the presence of chol and favored by chol in the case of β 2 AR. Chol has thus an opposite effect on the two receptors. If this prediction were to be correct, it would imply that cholesterol not only has the ability to affect the propensity and interface of GPCR oligomerization but also to do so on a receptor-specific manner. This ability might well constitute a natural strategy to selectively activate receptors based on their lipid environment, which may vary dynamically as a function of the cell membrane compartment they are embedded in. Although this latest work of Prasanna et al. is technically sound it is also tainted by the lack of any structural restrains to maintain the ternary structure of the receptors. The chol−receptors interactions might be impaired by the receptor deformation.
In an even more recent study, Prasanna et al. investigated the interactions of glycosphingolipids, the ganglioside GM1, with 5-HT A1 R. 137 GM1 is found to strongly interact with the protein surface, showing preferential sites. Notably, GM1 is found to bind at a consensus "sphingolipid binding domain" previously described. 138 This work is the first computational study to identify GM-GPCR interactions, and it illustrates the challenges of acquiring proper statistics for complex membrane where lipids stick to the protein more than POPC or cholesterol do. Techniques allowing more efficient sampling will be necessary for such systems.

Receptors Distribution in a Plasma Membrane-Like Bilayer
Sansom and co-workers contributed significant studies to the development and popularization of the Martini CGMD approach for membrane proteins [62][63][64]139 and their interactions with the matrix. 140−145 Notably, although they now routinely use the Martini force field and the ElNeDyn approach they developed their own version of the protein force field independently and keep deriving new topologies 146 to innovate in modeling using Martini CGMD simulations. 147 Recently, Koldso and Sansom studied the organization of lipids and proteins in a mimic of the plasma membrane, PM, in a couple of papers. 46,148 In their most recent work, 46 they addressed the organization and the dynamics of a plasma membrane model in three conditions: free of protein, PM system, containing multiple copies of a single transmembrane helix, gp130 system, the transmembrane domain of gp130, or of a GPCR, the sphingosine 1-phosphate receptor 1, S1P 1 system. Up to now this set of systems is by far the more complex simulated using a (pseudo) atomistic resolution. This would not be possible with any other model than the Martini CG model combined with MD simulations. The plasma membrane model used has an asymmetric lipid composition, PC:PE:Sph:GM3:chol and PC:PE:PS:PIP2:chol in the outer and inner leaflets, respectively, with the ratio 40:10:15:10:25 and 10:40:15:10:25. In the system containing gp130 and S1P 1 , 576 and 144 copies were inserted, respectively. The three systems were simulated for 10 μs, leading to system sizes after relaxation ranging from 104 to 125 nm. 46 Although the simulation length is arguably short (illustrated by the unsteady evolution of the oligomer state of the proteins at 10 μs), their analysis led to extremely valuable information and provided an unprecedented view of such biological system.
Koldso and Samson first noted the reduction of the membrane large-scale undulations observed in the PM system to smaller magnitude ones in the gp130 system and mostly absent in the S1P 1 system. The possibility of a kinetic effect (due to the slowdown of the lipid lateral diffusion observed in the presence of the proteins) was discarded since the diffusion was affected by a factor of ∼1.5 while the undulations by a factor larger than 5 (estimated from the similarity of magnitude at 2 and 10 μs in the PM and S1P 1 systems, respectively). Second, the lateral diffusion of the lipids in the simulations was found reduced by the presence of the proteins and more within the S1P 1 system than in the gp130 system. The slowing down of lipid dynamics in a crowded environment was predicted by both experimental 149,150 and computational 85,144,151−153 studies, which suggested that the effect is most pertinent for these lipids interacting specifically with the proteins and thus moving in tandem with them, leading us to the third observation in the Koldso and Samsom work relevant to us here: the predominant interaction of chol and PIP2 lipid molecules with S1P 1 . This was first notable by the larger reduction of their lateral diffusion. The authors quantified these interactions by a contact analysis revealing high contact frequency of PIP2 with basic residues on the intracellular side and of chol with most hydrophobic residues in the transmembrane span of S1P 1 ( Figure 8D). Both chol and PIP2 formed annular interactions with S1P 1 . 46 The head of chol molecules formed more specific contacts with the intracellular ends of TM1−3 and TM5−6 ( Figure 8D). Interestingly, the authors reported the presence of a cholesterol at the interface of S1P 1 dimers when using TM1 as an interface ( Figure 7A), and they discussed this observation in the context of the previously suggested relevance of chol for the stabilization of this interface in other GPCR, namely, β 2 AR 111 and mGlu1. 112 See below for more discussion.
Finally, the authors described the oligomerization pattern of S1P 1 in the PM. On a 10 μs time scale, they found S1P 1 largely as a monomer with transient formation of dimers, trimers, and higher oligomeric states. The author noted the contrast of this behavior with the one reported for rhodopsin, 80,81 β 2 AR, 82 opioid, 85 and A 2A and D 2 86 receptors where the receptors were predominantly forming linear arrays or filiform structures (Figures 2 and 9). At this point it is not clear if this difference results from (i) more dynamic or weaker receptor interactions for S1P 1 , (ii) their slower lateral diffusion and yet unsteady oligomeric distribution at 10 μs, or (iii) a competition between protein−protein and protein−lipid interactions. As suggested,  chol interactions might stabilize the S1P 1 interactions but also generate a kinetic barrier as proposed for opioid receptors. 85 It is also not clear how PIP2 is behaving. It might prevent the receptors from approaching each other as Sengupta and coworkers found chol behaving in the case of β 2 AR. 83

DHA Effect on GPCR Oligomerization
Guixa-Gonzaĺez et al. 86 recently used the Martini CGMD approach to investigate the mechanism by which docosahexaenoic acid (DHA) may affect the oligomerization of adenosine A 2A and dopamine D 2 GPCRs. DHA is an omega-3 polyunsaturated fatty acid (PUFA) 22 carbons long with 6 double bonds. DHA has been shown to be essential for proper brain function, and a low level of DHA in the brain was linked to patients with mental 154 and neurological disorders. 155,156 Similarly, earlier studies reported the importance of DHA for vision by affecting the function of the visual photoreceptor rhodopsin, a specialized GPCR. 157−163 These two biological processes have in common that they take place in cell membranes with an extremely high DHA content: 40 and 60% of PUFA in cerebral gray matter 164,165 and rod outer segments (ROS) membrane, respectively. 166 experiments to investigate if and by which mechanism DHA might affect the assembly of A 2A and D 2 receptors. The study consists mainly in the comparison of two conditions: low-and high-DHA-content membrane models corresponding to a disease-like and healthy-like patient, respectively. BRET experiments did not detect a difference in receptors interactions in low-and high-DHA-content membranes. In contrast, receptors self-assembly simulations showed a clear increase (∼20%) of receptor contacts in high-compared to low-DHAcontent membranes. The difference in time resolution of the two approaches, microseconds for Martini CGMD and milliseconds for BRET, 86 led the author to propose that DHA has a kinetic effect on the receptor assembly. CGMD simulations would capture this feature, while BRET would be blind to it. They further analyzed their simulations to characterize the mechanism by which DHA operates this kinetic effect.
The analysis of the simulations led the authors to the following observations. First, the receptors arrange in a linear array forming one or two contacts as reported for other receptors in similar time scales. [80][81][82]85 The author tested the stability of this filiform arrangement of the receptors on a longer time scale, possibly the longest to date: 260 μs. The receptors formed filiform structures extending toward shortbranched arrangements ( Figure 9A). Second, DHA avidly surrounds the receptor surfaces ( Figure 9A and 9B), which the authors confirmed by MD simulations at an atomistic resolution and in line with results on rhodopsin. 134,160,161,175 Third, in monomeric systems (only one protein) DHA increases the lateral and rotational diffusion of the protein but not of the lipids. Lipid diffusions, however, were lowered in the crowded protein environment (18 proteins) with the maximal effect observed for SDPCthe most proteininteracting lipid of the set present. These observations suggest that the lipids are following the protein dynamics and most in a crowded environment. Fourth, in a system with a higher DHA concentration (31%, high-DHA), associated with a 3 times increase of lateral diffusion, the author did not observe a faster protein oligomerization rate. The proteins actually took longer to associate. The authors interpreted these observations as the ability of a healthy-like membrane to form local or partial phase separation but not the high-DHA system. This local phase separation (DHA prefers DHA lipids) would be a possible mechanism for DHA to favor protein contacts in healthy membranes. DHA-coated proteins would come together easier in a locally/partially phase-separated system. Note, the authors used the terminology "local phase separation" to describe a phenomenon that would better be referred to as the formation of local domains enriched in particular lipid types.
This mechanism of DHA action on the A 2A R and D 2 R would be different than for rhodopsin where DHA was hypothesized to shift the MI−MII equilibrium by making more free volume available in the membrane and thereby the conformational change associated with the formation of MII less costly energetically. 157,158,176 In summary, the increase of receptor oligomerization observed in Martini CGMD simulations of healthy-like membranes results from a combination of an increase of the receptor translational and rotational diffusion, a decrease of the membrane surface exploration of the receptor due to partial or local phase separation (limited to DHA-rich regions), and an increase of the effective receptor-interacting diameter due to DHA coating and thus able to sense each other at longer distances. 86

Lipids Interactions with Rhodopsin
Grossfield and co-workers used the Martini CGMD simulation approach to looking at lipid distributions and interactions around rhodopsin. 134 This study followed a couple of publications on the subject using an atomistic resolution [160][161][162]175,177,178 and inspired by a large set of experimental data. 98,122,124,131,157−159,176,179−191 These works characterized the effect of the various lipid components specific to the rod outer segment (ROS) membrane and known to affect the function of rhodopsin. These effects are often measured by following the meta I−meta II (MI−MII) equilibrium, the last two photointermediates of rhodopsin. Only MII is able to bind its cognate G protein transducin. The MI−MII equilibrium is sensitive to ROS membrane properties such as its lipid headgroup composition (PC, PE, and PS), PUFA (DHA in particular), and cholesterol contents. These effects are nicely summarized in a paper by Horn et al. 134 and references therein. In short, MII (active rhodopsin) is favored by an increase in the negative curvature of the membrane (provided by PE and DHA) and an acidic membrane surface (provided by PS). Chol favors MI, but its content decreases as the disks age (become functional), and it stabilizes the disk membrane by compensating PE and DHA negative curvature. Horn et al. 134 built a Martini CG model similar the one they used previously in atomistic MD simulations: a single receptor embedded into a 2:2:1 molecular ratio of SDPC:SDPE:chol. [160][161][162]175 This particular composition aims at mimicking the ROS lipid composition. 179,192−196 Two systems were simulated: one with rhodopsin (based on the PDB ID 1U19 197 ) and one with opsin, representing the activated receptor (based on the PDB ID 3CAP 106 ). In the activated receptor, a conformational change involves the movements of the cytoplasmic side of TM5 (inward) and TM6 (outward) toward the helical bundle ( Figures 8A, 10A, and 10B). They performed 16 independent simulations of both systems each 1.6 μs, which represents an increase of system size by a factor of ∼3 and of simulation length by about 10 when compared to the atomistic data. In their analysis, Horn et al. demonstrated on many occasions the high degree of convergence of the sampling performed, illustrating one of the powers of Martini CGMD simulations: statistical significance of complex systems.

Chemical Reviews
Horn et al. performed detailed analysis of the simulations and presented it in a very clear and convincing manner. 134 The Martini CGMD simulations basically confirmed the trends observed in the atomistic simulations 161 but provided much better statistics leading to conclusive observations. The increased sampling also revealed chol interaction sites. First, DHA chains are found with higher density at the protein surface forming a ring. Second, stearoyl chains are excluded from the first shell into a second, most likely resulting from the attachment of both DHA and stearoyl chains to a unique headgroup. Third, a slight but statistically significant preference of PE headgroup was found for the protein surface with a possible preferred region of interaction for rhodopsin close to TM3/4/5 but less marked for opsin. Fourth, a cholesterol binding site was observed at the TM1/H8 interface behind the palmitoyl chains attached at Cys322 and 323 (Figure 8) in both rhodopsin and opsin. It interacts vigorously with the palmitoyl chains. This site corresponds to the location of cholesterol in crystal structures of other GPCRs. 110,111 Another site, found in both systems, is located at the cytoplasmic side of TM3/4/5. It was observed in atomistic simulations of adenosine A 2A receptor. 198 A third site was found in opsin but not in rhodopsin at the cytoplasmic end of TM5 and TM6. This region is marked by a large conformational change of TM5/6 and ICL3. This emphasizes the sensitivity of the cholesterol interaction with the receptor to the protein conformational details. Overall, the clusters and residues involved in cholesterol binding closely resemble the groups of residues identified in the previous atomistic simulations. 161 It should be noted that a cholesterol interaction site at the cytoplasmic side of TM1/2/3, previously identified by long atomistic simulations, 177 was not observed. Also, a few sites observed in β 2 AR and A 2A R in atomistic simulations were not found. It is actually not clear if they should.

Stability of Experimental and Theoretical GPCR Interfaces
During the reviewing process of this review, Baltoumas et al. 199 published a study based on Martini CGMD simulations that probe consistently GPCRs dimer interfaces derived from experimental or theoretical approaches. The interfaces of visual receptor rhodopsin, opioid, adrenergic and metabotropic glutamate receptors, and CXCR4 receptors from different species were thus studied in a systematic manner, resulting in a total of 21 systems simulated starting from known structures. Overall, the authors describe that the structural diversity of receptor interfaces found in experimental and theoretical models reduces in the simulations to converge toward consensus interfaces. The consensus conformations strongly resemble the ones previously described in the literature for Martini CGMD simulations, further emphasizing the relevance of TM1/2/H8 and TM5/6 interfaces, and TM4/5 to a lesser extent for GPCRs oligomerization. Most notably, the interfaces "loosely packed" or containing cholesterol molecules rear- Figure 10. Conformational heterogeneity in the TM5, ICL3, and TM6 region. Comparison of (A) two inactive rhodopsin (rho) structures, 1U19 and 2I35, and (B) inactive and active (rho*) rhodopsin, 2I35 and 2X72, respectively. Views from the membrane on TM6 and from the cytoplasmic side are shown. (C) Steric clashes at the TM5/6 dimer interface with the inactive structure of rhodopsin (1U19, center) and active rhodopsin (2X72, right). Dimer conformation built with the inactive rhodopsin structure 2I35 is shown for comparison (left). On this arrangement multiple salt bridges are formed and shown by yellow links between negatively (green) and positively (blue) charged side chains. Zoom in of the region where the interactions or the clash occurs is given for the three dimer's arrangements. ranged significantly, increasing protein contacts and minimizing the involvement of chol while keeping similar receptors orientation and improved protein packing. This study also revealed the presence of hydrogen bonds and aromatic residues (π−π) interactions at the interface of the receptors and close to the membrane/water interface. Finally, correlations between the dimer interfaces and regions of the receptors associated with its function were found, suggesting a possible regulatory mechanism of the function of GPCRs by their oligomerization state.

Methodological Insights: Pros and Cons of the Martini CGMD Approach
In the field of transmembrane proteins and most particularly for GPCR, Martini CGMD simulations have revealed a multitude of exciting findings in the extended range of size and time scales inaccessible to traditional all-atom approaches. Not surprisingly, these findings gravitated around the interplay between GPCRs and the membrane matrix (length scale in which Martini CGMD simulations excel) in relation to their oligomerization and supramolecular organization and the involvement of specific lipids (determinant for GPCRs).
The Martini CGMD simulations discussed in this review illustrate the main advantage of the method: providing an unprecedented view of the GPCRs embedded in membrane bilayer in different phases, enabling the rationalization of experimental data and revealing new features. PMFs may reveal the relative strength of interfaces between receptors, while selfassembly simulations indicate the most accessible. The complex composition of lipid bilayer allows unraveling preferences, binding sites, and protein/membrane biophysical and biochemical interplays.
There are, however, a few aspects of the method that are still limiting and some that need to be considered carefully. We discuss here the issues pertaining specifically to GPCR. General issues of the Martini CGMD approach have been described previously. 24,28,48,200 4.1.1. Conformational Restriction: Needed To Avoid the Collapse of the Receptor but Prevents Finding the Most Stable Interface in a Particular Condition. The conformational restriction imposed by the use of the ElNeDyn approach 65 is of great importance. We have shown in the example of rhodopsin (Figure 1) that the regular Martini model leads to a ∼6.5 Å rmsd of rhodopsin from the starting conformation. Others have shown similar deformations for β 2 AR 84 and rhodopsin. 134 This magnitude of deviation will be observed with any receptor and potentially any protein, globular or integral. ElNeDyn keeps it within ∼2−2.5 Å rmsd. Only a few exceptions have been reported, such as mechano-sensitive channels. 201 The drawback of the use of ElNeDyn is the actual impossibility of modeling significant conformational changes of the protein. ElNeDyn does, however, include fluctuations relative to small changes and reproduces the plasticity of the protein as an object. 65 However, conformational changes such as a receptor activation or even a loop conformational change cannot be observed in the simulations. Filizola and co-workers used a modified version of ElNeDyn where the force constant of the EN applied to the loop regions was reduced to reproduce the expected increased mobility. This approach, although appealing, has to be used with precaution. Reducing the force constant of the EN bonds as low as 250 kJ mol −1 nm −2 is not appropriate in the Martini force field. The bond may get so weak that it gets overpowered by the local nonbonded forces (LJ interactions), leading to the overlapping of consecutive beads. In addition, the behavior of unstructured protein segments in the Martini force field has not been investigated in terms of conformational flexibility. The loop conformational changes have thus to be considered with care. Alternatively, one could use a recent approach proposing to derive dihedral restrains from AA MD simulations to guide the conformational search of small amino acid sequences. 202 Work is in progress to develop a Martini model with flexible loop regions and more generally a polarized flexible backbone. 203 Including conformational details of loop regions in the search of GPCR interfaces is, however, important. One example lies in the work of Johnston and Filizola on the opioid receptors discussed in this review. 116 They indeed observed a significant increase in the stabilization of the TM1/2/H8 symmetric interface for kOR compared to any other GPCR reported ( Figure 4C and 4D). The outward pointing of TM1 from the helical bundle was suggested to be the cause of this increase. Another example of the sensitivity of the system to the protein structure resides in the interactions of chol with rhodopsin and opsin reported by Horn et al. 134 The different conformation of the TM5/ICL3/TM6 region in the two structures allowed chol to bind only to opsin ( Figure 8A).
We illustrate this point further with the PMFs of an alternative rhodopsin interface performed recently. We have been curious about why the symmetric interface using TM5/ 6found in the more recent structures of rhodopsin, opsin, μOR, and κOR and found more stable than TM1/H8 in μOR 116 was literally absent from our self-assembly calculations. 81 The comparison of the structure used in self-assembly simulations and the structure of opsin ( Figure 10A) reveals the presence of a steric clash that does not allow formation of the interface observed in the latter structure. The clash occurs between the intracellular loops 3 (ICL3) between TM5 and TM6 ( Figure 10C). The PMF calculated using a CG model based on this alternative ICL3 conformation indicates that the TM5/6 interface is more stable than the TM1/H8 one ( Figure  4E).
Of potential interest in the search of alternative approaches to ElNeDyn is the one coined Dom-ElNeDyn. 204 In this approach the EN to maintain the protein fold is decomposed in domains. These domains do not share elastic bonds and are thus free to move one relative to another. This could be used to allow a section of the receptor to move, i.e., TM5 and TM6. The techniques allowing heterogeneous ENs 66 or the combination of multiple ENs 205 could also be advantageous.
4.1.2. Still Some Limitation on the System Size and Time Simulated. Some of the simulations discussed in this study have pushed the limits of computational approaches in terms of both methodology and computer time dedicated to them. However, it is amazing how the field is evolving and that the first Martini CGMD simulations of GPCRs that contained 16 protomers simulated for 8 μs* 80 can be already regarded as relatively easy experiments. At the time it was performed it took a couple of months on a supercomputer to be completed and thereby represented a major achievement. Today one can run this type of simulation without major investment. Nowadays, major achievements consist of running systems containing from 4 81 to 9 46 times more proteins for about 10 81 to 30 86 times longer and in more realistic membrane systems containing (many) different lipids and protein types. 45,46,206,207 Although the systems are getting bigger and more complex, it is important to keep in mind that the biological processes that the Martini CGMD approach allows one to tackle often happen on longer time scales than actually currently simulated. A typical example is the lack of bind/unbinding events in the selfassembly simulation discussed in this review. One cannot interpret the interface populations in terms of relative stability. One has to perform PMF calculations in which the interfaces are predetermined by restraining the relative orientation of the proteins (see above). These simulations are quite costly and tricky to perform, as one has to check the behavior of the system in most windows within the range of receptor distances explored with the umbrella sampling. The amount of sampling (computer time) needed for reaching convergence in some windows is significant as the system might get kinetically trapped, e.g., due to trapping of lipids. This approach has also been suggested to overestimate the energetics 208−210 when sampling is not sufficient.
More realistic systems are also more computationally demanding and therefore limited to a few tens of microseconds. 45,46 4.1.3. Providing More Resolution Using Multiscale Approaches. The coarse grain resolution of the Martini model poses a limitation since the details of the atomistic interactions are lacking. When needed, atomistic resolution can be added using so-called backmapping methods. 211−214 The use of such technique has been applied to refine and/or validate protein/ lipid interfaces found by Martini CGMD simulations. 96 The early studies on rhodopsin 78,80 and later on β 2 AR 82 have shown heterogeneous deformation of the membrane around the receptors. These deformations are in perfect agreement with simulations using atomistic resolution. 100,217 Notably, Mondal et al. showed that the assumption of a radially symmetric deformation of the membrane at the protein surface in their 3D-CTMD method would severely underestimate the free energy cost of membrane deformation. They estimated the effect of heterogeneous deformation to almost 3 k B T. 100 The heterogeneity of the membrane deformations will depend on the membrane properties itself. On the basis of Martini CGMD simulations, we reported that the deformations of the membrane bilayer at the surface of rhodopsin vary with the membrane thickness ( Figure 3A and 3B). 80 Simulations using an atomistic resolution confirmed the delocalization of the deformations with the membrane thickness ( Figure 3C). 100 In the Martini CGMD simulations the regions at the protein surface where the lipid's interactions with the protein would deform the membrane the most correlate with the preferred regions of protein−protein contact, 80 in line with the concept that hydrophobic mismatch contributes to the oligomerization of membrane protein by reducing the energy cost due to membrane deformation. Weinstein, Khelashvili and co-workers extended this view by suggesting that in addition to the energy cost due to membrane deformation, which actually seems minimal for bilayers of physiological thickness, 100 the presence of residual hydrophobic mismatch (RHM, remaining hydrophobic mismatch due to incomplete matching by membrane deformation) would contribute significantly to the free energy of the system. 78,79,82,100 RHM results from the incomplete hydrophobic matching between the membrane and the protein when a too high cost for membrane deformation occurs or when two consecutive helices have drastically different hydrophobic properties. This would leave regions of the protein exposed to an unfavorable environment. They proposed that in the case of β 1 AR and β 2 AR the difference of RHM would explain their different oligomerization pattern observed experimentally. 119,120 In the case of rhodopsin the difference in RHM in lipid cubic and lamellar phases would rationalize the recent success of the in meso crystallization method. 78,79 LCP reduces the RHM on rhodopsin, thereby providing a stable environment. A subsequent destabilization of the LCP would drive the receptors to a lamellar phase for crystallization.

Lipid Trapping Results in Kinetic Barriers
. Two distinct types of interfaces were characterized in the case of rhodopsin. 81 One type of interface has a deep well at the distance of contact and no free energy barrier to assembly, e.g. TM1/H8, TM4/5, and TM5/6. The second type of interface does not experience a significant stabilization when formed and shows a clear energy barrier to their formation. The analysis of the simulations revealed the trapping of lipids at the interface stabilizing a metastable state where the interface is lubricated by a layer of lipids. 81 In the case of opioid receptors, Provasi et al. reported potentially related phenomena. 85 They observed heterogeneous lipid dynamics, and thus membrane viscosity, around the receptors. They could correlate the zones of slow lipids at the protein surface with low propensities (k on ) of the receptor to use that interface to assemble. The authors interpreted this correlation as the presence of a kinetic barrier to the formation of complexes with zones associated with jammed lipids. 4.2.3. Evidence of Specific Lipid Binding Sites Mediating Dimerization. Direct, specific, and localized interactions of cholesterol molecules at the surface of receptors have been reported in Martini CGMD studies for rhodopsin, 134 5-HT 1A , 132 and β 2 AR 83 ( Figure 8). These interactions are in line with reports using atomistic MD simulations, 161,177,198,218−220 receptor crystal structures, 110−112 and experimental data. In their study on β 2 AR, Sengupta and coworkers 83 suggested that cholesterol, through occupying certain locations on the protein surface with higher frequency, blocks interfaces from engaging in protein−protein contacts. Given the ubiquitous presence of cholesterol in biomembranes, this blocking behavior of cholesterol is of general importance for membrane protein complex formation. It would, therefore, be important to consolidate this observation by finding other occurrences of this behavior, maybe for other GPCRs. It is notable that this "blocking" behavior is complementary to the "gluing behavior" reported recently in the case of cardiolipin acting as a glue between the proteins constituting the respiratory chain complexes. 133,207 A gluing behavior of chol would also be consistent with the finding of its tight binding at the interface of GPCR dimer structures. 110−112

Several Factors
Contribute to GPCRs Self-Assembly into Linear Aggregates. All the contributions to receptor/membrane bilayer interplay described in the  101 and a potential alternative interface (Figure 4), respectively. (C and D) Binding mode of transducin (cognate G protein of rhodopsin) to the three row-of-dimer models shown in B using (C) the canonical orientation (can, build from the β 2 AR-Gαs complex structure 232 ) and (D) an alternative orientation (alt 233 ). See text for more details on the models. (E) Illustration of the pecking motion experienced by G protein (Gt) that it might use during its search for activated rhodopsin to overcome physical barriers. Lipid bilayer is show in gray (aliphatic chains) and blue (head groups). Gtαβγ trimer subunits are colored as α in green, β in red, and γ in yellow; full trimer is shown only every 8 μs*. α-Helical C-terminus of Gtα (Gα/αCT) is colored in orange and depicted so that the magnitude (∼2 nm) of its pecking motion is visible. previous subsections result in forces that will affect if not govern the propensity for the receptors to remain as monomers or to associate into dimers and higher ordered structures. Thus far, all studies of GPCR using the Martini CGMD approach show the receptors forming dimers and higher ordered structures with the exception of the PM-like membrane composition but again on a short time scale. 46 The receptors actually assemble with a predominance of linear arrays (filiform) structures 80−82,84−86 with the appearance of some small branched structures on a very long time scale. 86 The most straightforward explanation is the preference of the receptors to interact through the "small" sides of the receptors: centered either on TM1 or on TM5 (Figure 4). This behavior has been most clearly demonstrated in the case of rhodopsin. 81 It was shown that the interfaces involving TM1/H8 and TM5 (TM5, TM4/5) were found to be highly involved in the receptor interfaces formed in self-assembly simulations. The PMFs as a function of the receptor distance proved them to be much more stable than the other interfaces. 81 The interfaces involving TM4 (noted TM4/3 in other studies) and TM6 (TM6/7 might be a better representation) were not observed upon self-assembly. The PMFs of these interfaces revealed an energy barrier to their formation and a metastable state in which lipids lubricate the interface TM4. Linear aggregates are formed, most likely because the short sides can form direct contacts whereas the larger sides remain lubricated.
In their study of β 1 AR and β 2 AR, Mondal et al. rationalized the filiform organization of β 2 AR observed in self-assembly simulations on the basis of the presence of residual hydrophobic mismatch (RHM) on the "small" sides of the receptors, TM1/H8, and TM5 and not on the other larger sides. 82 It is not completely clear why the TM4/3 (also called TM4, see Figure 4) interface was not found with a notable RHM and more populated in the simulations since it was found to be relatively stable and barrier free by Johnston et al. 96 Furthermore, the RHM patterns on β 1 AR and β 2 AR allowed Mondal et al. 82 to rationalize the different behavior of these receptors in terms of oligomerization: while β 1 AR would mainly form dimers, β 2 AR would engage in large and more dynamic oligomers. 119,120 In the case of opioid receptors, Provasi et al. 85 also observed filiform structures mainly involving the small sides of the receptors: centered on TM1 (TM1/2/H8) and TM5 (TM4/5 and TM5/6). Unfortunately, the authors could not conclude on the relative strength of the interfaces formed. However, they did describe a striking anticorrelation between the involvement of an interface and the lipid dynamics (the membrane viscosity) at that interface. In other words, they propose the existence of a strong kinetics component in the determination of formation of the receptor contacts in self-assembly simulations. The fluidity of the lipids or viscosity of the membrane bilayer is shown to vary with the surface of the receptor and correlate with the kinetic of formation of interfaces (k on ). The more fluid the membrane, the more prompt the interface is to interact, favoring small interfaces. 4.2.5. Protein Burial Is Not Appropriate for Measuring Protein Interface Strength. The use of the protein burial associated with a protein interface in the membrane environment as commonly done for globular complexes appeared not to be a reliable tool to predict interface strengths. In the studies discussed in this review on GPCRs, the stronger interfaces have systematically a lower protein burial than weaker ones ( Figure  4). 81,96,116 This is at odds with common practice for soluble proteins. 221 It is likely that other forces are at play in membrane proteins. Many are discussed in this review.

Toward the Role of a Row-of-Dimers Organization for Rhodopsin Signaling
This highly ordered view of the rhodopsin organization is quite compelling ( Figure 11A). It is difficult to imagine rhodopsin freely diffusing as suggested by early biophysical experiments, but this mobility has been the object of debate that is out of the scope of this review. 102 It is also quite challenging to imagine transducin, Gt, rhodopsin cognate G protein, searching its way to the receptors, finding the activated one, and thereafter following the cascade of intracellular biochemical processes. These aspects have been discussed previously. 101,222−228 Here we discuss novel aspects provided by Martini CGMD simulations.
4.3.1. Row-of-Dimers Leads to a Unique Side of Rhodopsin Exposed to Bulk Membrane. One direct consequence of the highly symmetric supramolecular organization of rhodopsin that is not often considered in the functional models is that it leaves only one side of rhodopsin exposed to the membrane bulk, discarding the ends of the rows. If one wants to consider an alternative model as the sliding of Gt along the row-of-dimers, 226 the side of rhodopsin exposed to the membrane bulk is of great interest. It would be the only side viewed by Gt approaching. This exposed face is entirely determined by the conformation of the intra or functional dimer ( Figure 11B), emphasizing the importance of the search for the main interface for rhodopsin and potentially other GPCRs. In such organization, the search of Gt for an activated receptor, rho*, is radically simplified. It would ensure that Gt approaches the receptor by the same side. Note that the exposed side is not necessarily the one to bind, but it might contain the recognition side for Gt to engage in the binding mechanism.
In regard to the Gt search of an activated receptor, it is interesting to note the pecking motion observed in Martini CGMD simulation ( Figure 11E). This motion allows the Cterminus of its α-subunit, which binds the receptor within the groove open upon activation, to go over potential physical barriers on the order of 2−3 nm from the membrane surface ( Figure 11E). In the mean time, Gt is anchored to the membrane by its post-translational lipid modification of its α and γ subunits ( Figure 11). 4.3.2. Row-of-Dimer with the TM1/H8 Dimer Structure. On the basis of the Martini CGMD simulation work on rhodopsin we built a structural model following the row-ofdimer arrangement reported from AFM images and satisfying the cell dimensions determined from these images of rhodopsin in disk membrane prepared from mouse retina (Figure 11). 97 The duality of the interfaces (strong vs weak and lubricated, Figure 4E) fitted perfectly a model of rhodopsin organized in rows-of-dimers ( Figure 11A). The stronger interface, TM1/H8, was used as the main dimer interface (intra) and the weak ones, TM4 and TM6, were used in between the dimers (inter) in a row ( Figure 11B). Using the TM1/H8 interface as the main dimer interface was supported by its strength 81 and its presence in the ROS membrane. 109 This model satisfies the structural restrains extracted from the AFM images ( Figure 11A and  11B). 81 This model differs from the model proposed earlier by the authors of the AFM images, which utilizes the TM4/5 interface for the intradimer contact ( Figure 4) and TM2 and TM6 for the interdimer contacts. 101,226 Functional models have been derived from this particular arrangement. 101,226,227 4.3.3. TM5/6 Interface a Possibility for Rhodopsin? In searching for potentially stable interfaces for rhodopsin, 81 the TM5/6 did not appear as a potential candidate. It was only recently observed in an experimental structure of μOR 229 and was shown to be more stable than the TM1/H8 interface for μOR based on a PMF study using Martini CGMD simulations. 116 It was also reported to spontaneously form in self-assembly Martini CGMD simulations of homo-and heterodimeric interfaces of opioid receptors. 85 Also of potential interest is its involvement in a rearrangement of the interface of a family C GPCR. 230 The question is why did not the TM5/6 interface form in our extensive self-assembly simulations of rhodopsin? 81 As quickly discussed above, the construction of a model of a rhodopsin dimer reflecting the TM5/6 interface found for μOR (PDB ID 4DKL 229 ) revealed a severe clash between the ICL3, the intracellular loop between TM5 and TM6, of the rhodopsin model structure used in our Martini CGMD simulations (based on PDB ID 1U19 197 ). Among the rhodopsin structures deposited in the PDB, the most recent structures have an alternative conformation of TM5/ICL3 with an additional helical turn on the cytoplasmic side of TM5 and a different conformation of ICL3 (PDB ID 1GZM 231 and 2I35 107 ). Using such conformation for TM5/ICL3 (based on PDB ID 2I35) in a new CG model of rhodopsin in the dimer arrangement reflecting the TM5/6 orientation indicates the disappearance of the steric clash between ICL3 ( Figure 10C). The determination of the PMF of that clash-free interface demonstrated an intriguing stability of the TM5/6 interface ( Figure 4E). The interface seems stabilized by a double salt bridge on the cytoplasmic side ( Figure 10C), leading to a more stable interface than the TM1/H8 one, also shown for μOR. 116 At this point it is not clear how the dimer based on the TM5/6 interface could accommodate the conformational change associated with rhodopsin activation. There would be a steric clash between the two monomers upon activation of one of them ( Figure 10B and 10C). To avoid the unfavorable overlap, the movement of TM5 and TM6 upon rhodopsin activation would have to trigger either the dissociation of the dimer or its reorganization. Such reorganization is possible and has been recently proposed to explain FRET data for the metabotropic glutamate receptor. 230 4.3.4. Gt Binding to Rhodopsin in the Row-of-Dimers with Different Dimer Interface. Here, we search for structural clues from looking at how Gt fits on the rows-ofdimer models built with different rhodopsin interfaces. For this exercise we use the three most popular rhodopsin interfaces, TM1/H8 (stable in Martini CGMD simulations, seen in crystal structures and observed in ROS membrane), TM4/5 (built from AFM images and cross-link experiments), and TM5/6 (stable in Martini CGMD simulations and observed in crystal structures), and two activated rhodopsin/transducin complexes (rho*/Gtαβγ): one built according to the β 2 AR/Gs structure, 232 the canonical model (rho*/Gt can ), and an alternative orientation, 233 rho*/Gt alt . The later model was built prior the publication of the β 2 AR/Gs complex, and we use it here as an alternative model for a conceptual purpose.
The Mukhopadhyay et al. 233 motivation for building this alterative model, rho*/Gt alt , originates in that the structural model of opsin in complex with a variant of the C-terminus of alpha-subunit of Gt suggested a strong clash between Gtαβγ· GDP and the membrane bilayer, which was then interpreted by a large packing change of the alpha5 helix of Gt upon signal transduction from rhodopsin. 234 Mukhopadhyay et al. built the rho*/Gtαβγ·GDP from bits and pieces using available highresolution structures (see Mukhopadhyay et al. 233 for details). The main difference between the rho*/Gt can and the rho*/Gt alt is a rotation of the whole Gtαβγ by about 120°relative to rho while keeping the C-terminus of Gtα bound to rhodopsin. This rotation displaces the N-terminus helix of Gtα from the groove in between TM2 and TM4 (ICL1 and ICL2) in rho*/Gt can to the groove between TM4 and TM5 (ICL2 and ICL3) in rho*/ Gt alt .
The three row-of-dimer models are depicted together with the original AFM images 97 and a model resulting from a Martini CGMD study 81 (Figure 11A and 11B). The relative orientation of the receptors can be appreciated by H8 highlighted in orange. The comparison of the three models free of Gt leads to interesting observations. The first concerns the distribution of rhodopsin protrusions (TM6) marked by yellow transparent spheres in Figure 11B. While TM1/H8 (3.87 and 4.76 nm) and, to a lesser extent, TM4/5 (3.75 and 5.15 nm) models are close to the experimental distances 97 (3.8 and 4.2 nm) distribution ( Figure 11A), the model TM5/6 (3.91 and 1.97 nm) is not. The protrusions are too close. Second, the relative orientations of the pairs of dots in a dimer along rows in TM1/H8 and TM4/5 models indicate that while TM1/H8 respects the alignment found in the AFM images (γ = 85), TM4/5 does not. TM4/5 would have a left shift (γ = 95). Thus far the TM1/H8 model seems thus to be the best candidate.
The accommodation of rho*/Gt can and rho*/Gt alt by the three models also provides interesting observations ( Figure  11C and 11D). We will work here with two assumptions: one is that the lipid anchors of Gt (post-translational modification on the N-terminus of Gtα and C-terminus of Gtγ) should be embedded in the membrane bulk or close to it; the second is that the functional complex involves possibly a dimer of rhodopsin but not exclusively. 235−237 In the case of Gt can , only the TM5/6 model puts Gt lipid anchors in the bulk membrane. In the same TM5/6 model, some of Gtα would interact with the second monomer of the intradimer. In the TM1/H8 and TM4/5 models, the Gt can places its lipid anchors at least partially on top of the neighboring rhodopsin in the next dimer ( Figure 11C) as most of Gtα and Gtβ. Gt alt is intriguing because although it has a ∼120°rotation of Gt compared to the canonical binding, this rotation makes the TM1/H8 model a good candidate. Notably, it places Gt lipid anchors into the membrane bulk and removes clashes of Gtβ with the neighboring rhodopsins dimer. In the cases of TM4/5 and TM5/6 models, the lipid anchors would again be at least partially on top of the neighboring dimers. In the case of the TM4/5 model, the Gt alt binding mode places Gtαβγ in interaction with the rhodopsin dimer in an apparent alignment.
As a note, it is interesting to mention that if we had considered only a single rhodopsin dimer all interfaces and Gt binding modes would have accommodated Gt lipid anchors in the membrane bulk and only Gt alt bound to TM4/5 would have had important contact (potential clash) with the partner rhodopsin. The others would possibly accommodate Gt without major contact with the second receptor.
In summary, the analysis of Gt binding using two modes onto three models of rhodopsin in a row-of-dimer arrangement does not lead to convincing observations pointing to a most likely complex. All have at least a major improbable feature. This leaves us with the possibility that the rhodopsin might change its interface between active and inactive states, in a similar fashion as recently reported for another receptor. 230 The activated metarhodopsin could break the row-of-dimers and dissociate from it as a monomer that binds and activates transducin. However, no evidence has pointed to such a mechanism for rhodopsin.

Insights into the Oligomerization State of non-Rhodopsin GPCRs
The set of studies reviewed in this manuscript has unfortunately not brought many conclusive data on the actual oligomerization or supramolecular organization of GPCRs. It results from the relatively limited time scale of the simulations (maximum to date, 260 μs*) during which mostly assembly of receptors has been reported leading to linear or short-branched structures. Only studies using more complex membranes to mimic native lipid composition seem to picture GPCRs (S1P 1 ) as more dynamic but so far only from simulations with short time scale. 46 In the case of rhodopsin, expected to be highly ordered, the Martini CGMD simulations have not reported the spontaneous formation of such organization in the time scale accessible. The studies of non-rhodopsin GPCRs have however been able to provide clear preferences for certain interfaces leading to rationalization of experimental observations. An interesting aspect of the range of interface strengths among the GPCRs tested ( Figure 4) is the relative similitude of the interfaces probed. Only a few exceptions are markedly different and possibly not realistic. This is in contrast to rhodopsin, which seems to have a more discrete and defined set of interfaces (Figure 4). This difference might reflect the more general idea that GPCR oligomers are more dynamic while rhodopsin is more ordered.

CONCLUDING REMARKS
Overall, the studies described in the review represent an extremely valuable set of simulations of GPCRs. They reveal ways of GPCRs to interact with a membrane lipid bilayer, forces that affect their behavior as monomers, and trigger them to form oligomers and higher order structures, leading to new hypotheses for better functional models and inspiring the design of new experiments to probe them. We can only wish to get similar data for more GPCRs and extracted more systematically in order for them to be more conclusive when combined. These studies also highlighted generic protein/lipid interplays that should extend to membrane protein biophysics in general.
It makes no doubt that the membrane environment is a very complex media 238 that goes beyond the original fluid mosaic model 239 and that Martini CGMD simulations will be an essential tool in uncovering the principles governing this fundamental challenge in biological chemistry. 240 The nature of the data produced by this method will increase in complexity to become closer and closer to realistic membrane compositions with time scales reaching experimental observables. It is incredibly exciting to envision the data that will be accessible in 5−10 years.

Notes
The author declares no competing financial interest.

ACKNOWLEDGMENTS
This work was supported by The Netherlands Organization for Scientific Research (NWO). I thank Thomas Huber, Siewert Jan Marrink, and Tom Sakmar for comments on the manuscript. I also thank S.J.M. for his unconditional support and scientific inspiration and T.S. for his constant encouragement and belief in the power of the CGMD method, and I would like to particularly acknowledge the invaluable contributions of T.H. in the development of the Martini CGMD simulations for rhodopsin and the incredible satisfaction and pleasure it is to share the venture and credits of this research with him.