GPCR Oligomerisation Modulation by Conformational State and Lipid Interactions Revealed by MD Simulations and Markov Models

G protein coupled receptors (GPCRs) play key roles in cellular signalling pathways. GPCRs are suggested to form dimers and higher order oligomers in response to activation and other stimuli. However, a lack of structural details of the involvement of GPCR oligomerisation in receptor activation has held back our understanding of GPCR signalling in an in vivo context. Here, we combined extensive molecular dynamics simulations with Markov state modelling in order to understand the oligomerisation process of the Adenosine A2a receptor in response to receptor activation. The Markov state models reveal that activation leads to an increased propensity for oligomerisation, accompanied by a more interconnected oligomerisation network, when compared to the receptor in the inactive state. Coupling of the active state to a mini Gs protein shifted the oligomerisation pattern. Oligomerisation is also sensitive to changes in local membrane environment, such as the packing density of GPCRs in a membrane patch, which enhanced oligomerisation, and removal of PIP2 from the lipid bilayer membrane which disfavoured oligomerisation. These results support combinatory allosteric modulation of GPCR oligomerisation whereby the receptor may generate a specific oligomerisation profile in response to its local environment in order to favour the formation of a specific supra-molecular signalling complex.


Introduction
G Protein-Coupled Receptors (GPCRs) are the largest superfamily in the human genome. They are major drug targets due to involvement in many physiological processes. The canonical model of GPCR activation involves the formation of a ternary structure of a receptor, an agonist and a signalling partner 1 . However, this paradigm may require revision given increasing evidence of the existence of GPCR oligomerization and its biological importance 2 . For example, single-molecule imaging has revealed that a range of GPCRs (e.g. the dopamine D2 receptor and the M2 muscarinic receptor) can form metastable oligomers in living cells with lifetimes ranging from sub-milliseconds to seconds 3,4 . The degree of oligomerization varies among receptors. For example, at moderate expression levels, β2 adrenergic receptors form dimers, whilst cannabinoid CB1 receptor and opsin form tetramers 5 . The degree of oligomerization is positively correlated with receptor density and with the presence of an agonist for some GPCRs, including the serotonin 5-HT2c 6 , dopamine D2 7 and secretin receptors 8 . The biological importance of GPCR oligomerization may lie in the possibility of functional communication between adjacent protomers (i.e. individual GPCR molecules) within these assemblies. Such communication is suggested e.g. by bivalent drugs that simultaneously bind to two receptor molecules to generate diverse pharmacological profiles 9,10 , and by rescue of functionally deficient luteinizing hormone receptor via intramolecular interactions 11 .
More recently, super-resolution microscopy revealed that distinct clusters of GLP-1R (a Class B GPCR) form nanodomains to initiate endocytosis 12 , and that G protein coupling and β-arrestin mediated internalization are associated with the formation of CXCR4 oligomers 13 . A recent singlemolecule study links formation of µ-opioid receptor dimers with specific agonists and their downstream signals 14 .
The detection of G protein oligomers has pointed towards a scenario in which a signalling complex may consist of an oligomeric assembly of both GPCRs and G proteins. Single particle analyses showed that a heterotetramer of A1-A2a receptors exists as a dimer of dimers which simultaneously couples to both a Gi and a Gs protein 15 , and that addition of agonist resulted in a supramolecular complex of four receptors and four G proteins 16 . Similar to GPCR oligomerization, that of G proteins is dynamic and subject to modulation by such factors as palmitoylation 17 and the presence of ligands 18 .
cholesterol molecules and some fatty acid chains, are frequently seen at such dimer interfaces in e.g.
GPCR_oligo_main_text_v16 biorxiv.docx  P1Y12, A2aR and β2AR and others 26 . Indeed, cholesterol was proposed to have an impact on GPCR oligomerization either via alteration of bilayer properties, e.g. thickness, fluidity, acyl tail order, and/or by modulating the potential transmembrane interactions surface of the receptors by binding at specific lipid sites 27 .
Molecular dynamics (MD) simulations, in particular coarse-grained (CG) simulations, provide a useful tool to investigate protein:protein interfaces and the impact of membrane lipids on GPCR oligomerization. Rhodopsin was found to form interfaces with different association affinities corresponding to TM5/TM5 and TM1,2,8/TM1,2,8 interactions in CG simulations in POPC membranes 28 . Addition of cholesterol molecules to pure phospholipid membranes shifted the dominant dimer interface of CXCR4 from TM1/TM5-7 to TM3,4/TM3,4 in CG simulations 29 .
Increasing the concentration of cholesterol in POPC membranes correlates with enhanced plasticity and flexibility of dimeric association of the serotonin1A receptor, a modulatory effect exerted via both direct cholesterol binding and indirect membrane effects in CG simulations 30 . Kinetics of the homodimerization of opioid receptors have also been explored via CG simulations in a POPC/CHOL membrane and Markov state models, revealing that dimerization of µOR is dynamic with kinetically distinct dimeric assemblies 31,32 .
Despite these efforts, our current understanding of GPCR oligomerization still remains very limited.
In particular, we remain unclear how GPCR oligomerization responds to receptor activation, and how lipid interactions in more complex membranes containing multiple lipid species may modulate such responses. These missing pieces of the puzzle hold us back from a fuller understanding of GPCR signalling in an in vivo context, taking place within a crowded and complex membrane environment.
Here, we present a more complete view of GPCR activation in which receptor oligomerization, G protein coupling and lipid interactions are interconnected. We have employed Markov State Models (MSMs) to investigate the oligomerization behaviour of a Class A GPCR (the A2aR). By combining extensive ensembles of unbiased MD simulation data (over 2.6 ms of CG-MD in total) with MSMs we demonstrate that receptor activation results in an increased propensity to form oligomers when compared to the inactive state (in which the receptor remains predominantly in a monomeric form with occasional dimerization). Our results suggest potential combinatory allosteric modulation of GPCR signalling, in which the receptor may respond to the surrounding membrane environment to generate a unique oligomerization profile favouring the formation of specific supramolecular complexes. This helps us to understand the cell membrane-level context of the structural details of GPCR signalling complexity, presenting new possibilities to manipulate GPCR function.

MD simulations
We based initial simulations on large-scale membrane systems using a mixture of lipids to form an in vivo mimetic model of the plasma membrane environment within which multiple copies of the receptor can freely move, enabling both association and dissociation events to occur. We have focused on oligomerization of the Adenosine A2a receptor (A2aR) as this receptor has been demonstrated experimentally to form dimers and higher order oligomers 33,34 . Thus we placed 9 copies (positioned and oriented randomly in membrane plane relative to one another) of the A2aR in a 45 x 45 nm 2 membrane with 10 lipid species present (Fig 1A) to simulate the oligomerization process (i.e. the 9-copy systems). To study how the oligomerization may change in response to receptor activation, we generated 3 ensembles of simulations for the 9-copy systems, one each for the receptor in the inactive state (PDB id 3EML) , the active state (PDB id 5G53 receptor only) and the active state in complex with a mini Gs protein (PDB id 5G53 receptor and mini Gs). For each conformational state 10 repeat simulations of 50 µs duration were performed. We also performed simulations of the three conformational states with a higher protein density (~4% by area) by including of 16 copies of the receptor in a membrane of the same dimensions (termed 16-copy), as a number of studies (e.g. 5 ) have shown that altering receptor density can alter GPCR oligomerization.
Given recent studies, both computational and experimental [35][36][37] , demonstrating the interaction of phosphatidylinositol 4,5-bisphosphate (PIP2) with Class A GPCRs, we also performed simulations with setups identical to the 9-copy system but with a lipid bilayer devoid of PIP2 (9-copy NoPIP2) to see whether binding of PIP2 could modulate GPCRs via altering oligomerization.
Overall, 60 MD trajectories with a cumulative simulation time of 2.75 ms were generated ( Table 1).
The time evolution of the minimum distance between pairs of receptors illustrated that multiple association and dissociation events occur within the timescales simulated, and protein-protein interactions of different durations were sampled ( Fig. 2A and SI Fig. S1). On average, each trajectory sampled over 20 association or dissociation events in the 9-copy systems and over 60 in the 16-copy systems (Fig. 2B). The total number of association events in the 16-copy systems was approximately (16/9) 2 times higher than in the 9-copy system, indicating a second order reaction. Ensemble averages of the number of monomers, dimers, and higher order oligomers present suggested that all nine simulation systems approached an equilibrium towards the end of simulations, the characteristics of which were dependent on the membrane environment and the receptor conformational state ( Fig. 3 and SI Fig. S2). For example, comparing the 9-copy simulations of the Inactive vs. Active states, it can be seen that after 25 to 30 μs, the number of monomers drops to ~5 for the Inactive state compared to 2 or 3 for the Active state. The effect of omitting PIP2 from the simulation of the Active state of the receptor seems to decrease the receptor tendency to oligomerise such that the number of monomers remained at~5 (Fig. 3C). Overall, given these initial observations, we decided to analyse the kinetics of association and dissociation observed in the simulations in more detail using Markov State Models (MSMs) 38 .

Markov State Models
MSMs can be used to model oligomerisation in terms of a network of transitions between oligomeric states, based on information extracted from MD simulations 38 . The MD simulations were analysed in terms of oligomer formation (using a distance cut-off of 0.75 nm and hierarchical clustering)in order to build MSMs for each of the nine simulation ensembles ( Fig. 1B; see Methods for details).
In order to characterise the oligomeric states, clustering was performed of the oligomer structures occurring in the MD simulations; this was implemented separately within each oligomeric order if the order was lower than or equal to 5, whilst oligomers of order higher than 5 (i.e. hexamers and above) were distinguished only by their oligomeric orders. An MSM was estimated for each set of simulations using pyEMMA 39 , which reweights the transitions such that the equilibrium kinetics and stationary distributions can be recovered. After statistical validation, the MSMs were used to compute metastable conformations, equilibrium probabilities, kinetics, and oligomerization networks (see Methods for details). The following analyses (of equilibrium distributions, oligomerization kinetics and pathways) are based on the information obtained from the MSMs.

Receptor activation and packing density led to unique oligomerization fingerprints
The equilibrium distributions calculated from MSMs (see Table S1) support the tentative conclusions drawn from examination of the raw data (see above). Thus, it is clear that the active state favours oligomerization compared to the inactive state (e.g. in the 9-copy systems, for the active state, 25% of the population was in oligomeric form (including dimers and higher orders) compared to 7% for the inactive state). Furthermore, the active state is also associated with a shift in favour of higher order oligomers. Not surprisingly, comparison of the 9-copy and 16-copy simulations indicates that a higher receptor density in the membrane favours oligomerization. The presence of the mini-Gs protein bound to the receptor also noticeably shifted the oligomerization profiles in favour of higher order oligomers. Conversely, omission of PIP2 from the bilayer disfavours oligomerization, especially to form oligomers of order higher than 3 (albeit in shorter simulations). Taken together, the oligomerization fingerprints of equilibrium distribution suggest that A2aR activation may induce higher order oligomerization in a PIP2-dependent fashion. The positive correlation of oligomerization with both receptor activation and protein density is in agreement with observations on a number of GPCRs 6-8 .

Activation increased oligomeric plasticity
To better understand GPCR oligomerization in response to receptor activation, we characterized the oligomeric assemblies and compared the distributions of these assemblies between different conformational states. For the 9-copy systems, we calculated the distributions of oligomeric assemblies of order lower than pentamers (corresponding to 99.97%, 99.93% and 91.68% of the total populations of the inactive, active, and active + mini Gs states respectively). To describe these oligomers, we used the following metrics ( Fig. 4): (i) for dimers, the relative orientation angle pair (θ1, θ2) 40 , in which θ1 describes the orientation of the given monomer with respect to a reference monomer and θ2 describes the rotation of the second monomer around its z axis to match the reference monomer ( Fig 4A); (ii) for trimers, a bending angle φ defined by the centres of mass of the three monomers ( Fig 4B); (iii) for tetramers and pentamers, the projected lengths D1 and D2 of their 1 st and 2 nd principal axes onto the x and y axes ( Fig 4C). Thus, overall comparisons of the distributions of these three metrics for oligomeric assemblies from different simulations suggested that when in the active state the receptor showed higher oligomeric plasticity than when in the inactive state. This was the case whether or not the mini-Gs was bound to the active state receptor.
Comparable trends were observed in the 16-copy simulations, i.e. higher oligomerization plasticity was shown in the active state in all oligomeric orders (SI Fig. S3). This consistency suggests that the increased plasticity may be a consequence of the conformational changes during receptor activation, GPCR_oligo_main_text_v16 biorxiv.docx 24-Jun-20 8 which involve a prominent outward tilt of TM6 in the intracellular leaflet of the membrane. This would be expected to enlarge the potential oligomerisation surface of the receptor thus accommodating more extensive protein-protein associations in more diverse combinations. In the absence of PIP2 (see SI Fig. S4) the oligomeric distributions are altered relative to those from the simulations in the presence of PIP2, and in particular the absence of PIP2 seems to be associated with a degree of reduction in overall oligomeric plasticity.

The active state receptor exhibits a more diverse combination of association interfaces
Characterization of these oligomeric assemblies allowed us to determine the association interfaces and their relative frequencies from the MSMs. Eight association interfaces were identified in the 9copy systems (  Gs, whereas interfaces ICL3,TM5,TM6//H8,ICL1,TM1 and ICL3,TM5,TM6//ICL3,TM5,TM6 occurred more frequently due to the ample space provided for the two mini Gs by the tilting of helices TM6 and TM5.

Enhanced PIP2 interactions led to more stable protein-protein associations
Visual inspection of the MD simulation trajectories showed that PIP2 molecules were seen at the association interfaces mediating the protein-protein interactions ( Fig 6A). As it is difficult to incorporate lipid interactions in MSMs, we calculated the PIP2 interactions directly from simulation trajectories. Six PIP2 binding sites, all of which formed the core of the association interfaces, were observed ( Fig. 6B; We also examined the association interfaces in the 9-copy NoPIP2 systems from which the influence of PIP2 interactions was absent (SI Fig S6). Two prominent changes showed up when comparing the oligomerization interfaces in these two membrane environments: firstly the multiplicity of the interfaces decreased in the 9-copy NoPIP2 systems, indicating that the binding of PIP2 altered the association energy landscape to allow more patterns of association; secondly, the dominant interfaces shifted from those centring around the PIP2 binding site of ICL3_TM5_TM6 to those involving H8, TM1 and TM7 in NoPIP2 systems, indicating the lipid-mediated stabilization of this interaction.
Interestingly, the association interface H8,TM1,TM7//H8,TM1,TM7, which saw a significant increase of occurrence in the 9-copy NoPIP2 systems (Fig 6C), has been seen to bind to cholesterol both in MD simulations 29,35 and in crystal structures 19,41 . Thus the stability of this interface would be dependent on cholesterol-mediated interactions.

Active state showed a more interconnected oligomerization network
A simulated 10 ms Monte Carlo MSM trajectory was generated for the 9-copy ensemble (see e.g. Fig   7A and 7B). The resultant trajectories revealed a dynamic process of oligomerization in which A2aR switched between oligomeric states at sub-millisecond timescales and switched to more stable oligomeric assemblies when activated and when coupled to mini Gs. The trajectories showed a higher occurrence of dimers and higher order oligomers in the active state comparing to the inactive state.
We then calculated the lifetime and on-/off-rate based on the mean first passage time (MFPT). The computed lifetimes of monomer and dimer were 388 µs and 36 µs respectively for the inactive state compared with 359 µs and 380 µs respectively for the active state (SI Table S2), indicating a 10x stabilization of the dimer upon activation. The reaction rates (SI Tables S3 and S4) also revealed a more interconnected oligomerization network in the active state, with more frequent transitions seen between higher order oligomer and dimers, due to the higher stability of the active dimers (Fig. 7C).
In contrast, for the inactive state the high order oligomers tended to directly dissociate to monomers,.
Comparable trends were seen for the 16-copy simulation systems, i.e. a more interconnected oligomerization network for the active state with frequent transitions seen between different orders of oligomers and longer lifetimes for the dominant oligomeric states (SI Fig S7).

Discussion
By combining extensive MD simulations with Markov State models, we have explored the oligomerization process of the A2a receptor as a function of receptor activation state and density of receptors in a bilayer. This was made possible by the use of large scale CG simulation systems containing multiple copies of the receptors together with unbiased simulations on multi-microsecond timescales. One of the key findings from this study is that the oligomerization plasticity and dynamics are significantly increased when the A2aR is in the active state compared to when it is in the inactive state. Whilst the inactive receptor remained predominantly in the monomeric form whilst occasionally dimerising (with a large koff), the activated receptor showed increased population and longer lifetimes of dimers, alongside a more interconnected oligomerization network in which transitions between different orders of oligomers occurred more frequently. Comparing to the inactive state, the activated A2aR also exhibited a more diverse ensemble of oligomeric assemblies, resulting from an increased range of combinations of association interfaces which ultimately resulted from the increased intracellular surface around the outward tilted TM6. This was underscored by the shift towards TM6-associated interfaces in the active state oligomers. Given that this outward tilt of TM6 is a common feature in Class A GPCR activation, we suggest that the increase of oligomerization plasticity and dynamics might be an element of a shared mechanism for receptor activation. Our results also suggest an explanation of the difficulty in observing dimers or higher order oligomers in crystals.
The interactions of membrane lipids, e.g. cholesterols and polyunsaturated fatty acids, have been reported to affect GPCR oligomerization 42,43 , via altering the receptor surface and intercalating between protomers 27 . Here we report the effect of PIP2, which has previously been shown to interact with a number of Class A GPCRs 35-37 , on A2aR oligomerization. We saw from the simulations that PIP2 bound to the conserved PIP2 binding sites at the intracellular side of A2aR intercalating between protomers and enhancing protein-protein interactions. A direct correlation between the probability of association interfaces and the probability of PIP2 binding at these locations was revealed by comparing to the PIP2 interactions in our previous simulations of single copy of the receptor 35 . This was further confirmed by comparing to systems devoid of PIP2. PIP2 binding also contributed to the shift towards the TM6-associated interfaces of oligomerization in the active state. Given that our previous simulation studies have shown that the interactions of PIP2 are largely conserved across Class A GPCRs 36 , it is tempting to suggest that the effect of PIP2 on GPCR oligomerisation may also be shared. Such interactions could be important in determining GPCR oligomerisation e.g. in functionally important microdomains with a high local concentration of PIP2 44 . Interestingly, mutation of intracellular cationic residues, particularly those on ICL3, have been shown to play a critical role in formation of quaternary structures related to the functions of GPCR oligomers 45 .
Recent single-molecule imaging techniques have revealed that G proteins recruited by monomers and dimers of the β2AR elicited different signalling pathways 46 . The block of recruitment of G protein to β2AR dimers by the inverse agonist isoproterenol suggested that the dimers in this case were responsible for the basal activity of β2AR. The specific effect of FTY720-P, the active form of the pro-drug fingolimod, on S1PR1 was demonstrated to be dependent on receptor oligomerization causing a short-lived Gαi-mediated intracellular response that differed from that of the endogenous agonist S1P 47 . These studies suggest that the oligomerization of GPCRs may initiate signalling pathways different from those activated by the corresponding monomers. Signalling outcomes may also be dependent on the specific arrangement within oligomeric assemblies. Through a peptideinterfacing approach, A2aR-D2R heterotetramers have been shown to have different oligomeric interfaces when coupled to AC5 in the presence or absence of agonist 48 .
A number of simulation studies using the MARTINI force field have provided valuable insights into the association of membrane proteins 28,32,[49][50][51][52][53] and into protein-lipid interactions [54][55][56] . However, the current MARTINI model may also have some limitations 57,58 . In particular, the lack of size-specific ߝ for Lennard-Jones (LJ) potentials and the weak bond constants may sometimes lead to an overestimation of direct protein-protein interactions in the aqueous phase. Whilst the Inactive and Active simulation ensembles showed reasonable association-dissociation behaviour, in part due to the involvement of lipids in the protein-protein interface, the Active + mini Gs ensembles displayed a lower frequency of dissociation resulting from direct interactions of mini Gs in the aqueous phase.
This low sampling of dissociations in turn resulted in large uncertainties for properties calculated from MSMs of the Active + mini Gs state. Future studies of GPCR oligomerisation using MARTINI 3 (http://cgmartini.nl/index.php/martini3beta) which includes size-dependent LJ potentials may provide further information on e.g. the impact of G protein coupling on GPCR oligomerisation.
In the current study we have employed relatively simple models of the lipid mixture present in mammalian plasma membranes, as in our previous studies of GPCRs 35,59 . There have been recent advances in computational studies of more complex models of the lipid bilayer composition of cell GPCR_oligo_main_text_v16 biorxiv.docx 24-Jun-20 membranes 60,61 alongside experimental studies of the asymmetric lipid composition of plasma membranes 62 which will enable such models to be used in GPCR simulation studies in the future.

Conclusions
Reflecting upon our simulation results in the context of experimental studies [46][47][48] , we suggest that GPCR oligomerization may provide greater functional flexibility for the receptor signalling array via generating different supramolecular complexes that could initiate different signalling outcomes (Fig   8). From our MSM, we observed enhanced oligomerization plasticity facilitating the transition to oligomeric assemblies of the G protein-coupled state due to a larger sampling of oligomeric interfaces and faster transition dynamics between assemblies in the active state. In addition, distinctive oligomerization fingerprints, i.e. the level of oligomerization and the probabilities of various oligomeric assemblies, were elicited when the active A2aR was in different membrane conditions, which suggested that the oligomerization energy landscape was sensitive to the membrane environment. A number of studies have revealed that the coupling mechanism between GPCRs and G proteins is highly efficient, involving direct collisions with no intermediate 63 Table 1) using the insane.py script 67 . Electrolyte solution corresponding to~0.15 M NaCl was added and additional GPCR_oligo_main_text_v16 biorxiv.docx 24-Jun-20 13 sodium ions were added to neutralise the system. All the simulations were performed using GROMACS 5.1 68 . The CG simulations parameters were taken from ref 35 . A summary of simulations performed is provided in Table 1, amounting to a total of~2.6 ms of simulation data collected.
Oligomeric assembly characterisation. The minimum distance between each pair of proteins was monitored as a function of time (SI Fig. S1), from which the density of such distance distribution was estimated. The density estimation showed a first peak at~0.55 nm and a first trough at~0.7 nm. We found that a cut-off of 0.75 nm could best discriminate effective associations. Using this cut-off and hierarchical clustering, the number of oligomers of different orders were counted and each copy of the receptor was assigned to those identified oligomers separately for each frame. The identified oligomeric assemblies with the same oligomeric order were grouped together into an 'oligomer pool', and characterised, i.e. clustered, within each pool if the order is lower than 6. Since calculation of root mean square deviation (RMSD) based on structural coordinates is sensitive to the order of the receptor copies, we set out to assign an order of copies to the structures. We took a reference structures

MSM Analysis.
To assist understanding of the oligomerisation mechanism, the microstates resulted from Kmeans clustering were lumped into macrostates of different oligomeric orders SI where I = 1,2,..,9 in the 9-copy systems, and 1,2,…,16 in the 16-copy systems. The macrostate transition probabilities were obtained from the microstate probabilities with: is the macrostate equilibrium probability. The mean first passage time (mfpt) between the macrostates was calculated as described in ref 70 via the functionality in pyEMMA. The lifetime of a macrostate SI was calculated as mfpt(I,J), where ‫ܬ‬ ∈ Θ and ‫∉ܬ‬ ‫.ܫ‬ Lipid Interaction Analysis. Lipid binding sites, interaction counts and residence times were calculated using PyLipID (https://github.com/wlsong/PyLipID). A continuous lipid contact was defined using a dual cut-off scheme such that the interaction started when lipid headgroup moved to a given residue closer than 0.55 nm and ended when the headgroup moved farther than 1.0 nm. Lipid binding sites were derived from the lipid contact information using graph theory and community analysis. Graphs were constructed using residues as nodes and the frequency with which a given pair of residues interacted with the same lipid molecule as edges. The best partition function of the community library (https://github.com/taynaud/python-louvain) was used to detect community structures. Each calculated community was considered a binding site.
Lipid residence time of a binding site was calculated from the normalised survival time-correlation where ܶ is the total simulation time, Nj the total number of lipids of a given type, and ݊ (ߥ, ‫ݒ‬ + ‫)ݐ‬ is a function that takes the value 1 if a lipid occupies the given binding site for a duration of ‫ݐ‬after forming the contact at time ‫,ݒ‬ and takes the value 0 if otherwise.  Fig. 6B, to which the colour coding used above corresponds. * I = inactive; A = active; AmG = active + mini Gs.  Table 1. (B) Schematic workflow for constriction of the Markov State Model. The oligomers from an ensemble of simulation repeats were characterised and a transition matrix between the receptor monomers, dimers and these oligomers was calculated to generate a Markov state model using pyEMMA 39 . The number of association (defined as a pair of proteins coming closer than 0.75 nm) and dissociation (defined as a pair of proteins separating to further than 0.75 nm) events sampled in each trajectory for each of the 9-copy and 16-copy system ensembles. (C) The associations and dissociations have led to a dynamic equilibrium as illustrated by the time course of oligomer formation for the inactive state 9-copy, active state 9-copy, and active state 9-copy NoPIP2 simulations. Data from all the trajectories in the same simulation ensemble were plotted out to illustrate the ensemble trend.    (A) Snapshot from the 9-copy active simulations of the interaction interface ICL2,TM5//H8,TM1 with mediating PIP2 molecules shown in ball and stick. The receptor backbone beads were shown in cyan. (B) PIP2 binding sites on the active A2aR structure. The 6 binding sites were calculated using PyLipID (https://github.com/wlsong/PyLipID), and shown using different colours. The residues in each binding site are shown as atomic spheres with sphere scale proportional to their PIP2 residence times for the active state 9-copy simulations. (C) Comparison of the relative frequency of the H8,TM1,TM7//H8,TM1,TM7 interaction interface in the 9-copy systems and 9-copy NoPIP2 systems.  The activated GPCR generates an array of oligomeric assemblies which couples to signalling partners, generating an array of supramolecular signalling complexes capable of initiating different signalling pathways. This mechanism of combinatory modulation of GPCR is responsive to the lipid bilayer environment and the conformational state of the receptor. Receptor activation, PIP2 and/or mini-Gs binding, and an increase in receptor density in the bilayer all promote oligomerization.