Preferential Coupling of Dopamine D2S and D2L Receptor Isoforms with Gi1 and Gi2 Proteins—In Silico Study

The dopamine D2 receptor belongs to rhodopsin-like G protein-coupled receptors (GPCRs) and it is an important molecular target for the treatment of many disorders, including schizophrenia and Parkinson’s disease. Here, computational methods were used to construct the full models of the dopamine D2 receptor short (D2S) and long (D2L) isoforms (differing with 29 amino acids insertion in the third intracellular loop, ICL3) and to study their coupling with Gi1 and Gi2 proteins. It was found that the D2L isoform preferentially couples with the Gi2 protein and D2S isoform with the Gi1 protein, which is in accordance with experimental data. Our findings give mechanistic insight into the interplay between isoforms of dopamine D2 receptors and Gi proteins subtypes, which is important to understand signaling by these receptors and their mediation by pharmaceuticals, in particular psychotic and antipsychotic agents.


Introduction
Dopamine receptors belong to rhodopsin-like G protein-coupled receptors (GPCRs) and share the molecular architecture typical for this family of proteins. They consist of seven transmembrane (TM) helices, which are connected by three extracellular and three intracellular loops (ECL1, ECL2, ECL3, and ICL1, ICL2, ICL3, respectively). Dopamine receptors are divided into D 1 -like (D 1 and D 5 receptors) and D 2 -like (D 2 , D 3 , and D 4 receptors) subfamilies based on their activation or inhibition of adenylyl cyclase. These aminergic GPCRs are involved in regulation of many physiological conditions, ranging from voluntary movement and reward to hormonal regulation and hypertension [1], and thus they are important molecular targets for the treatment of a number of diseases: Parkinson's disease, restless legs syndrome, sexual dysfunction, dementia, depression, bipolar disorder, Huntington's disease, and schizophrenia [2]. Among dopamine receptor subtypes, D 2 is the most widely explored in medicinal chemistry, being a target for antipsychotics, drugs against Parkinson's disease, and many other drugs.
The dopamine D 2 receptor inhibits adenylyl cyclase as it is coupled to G i/o protein. Two alternatively spliced variants are produced from the D 2 receptor gene and code for the D 2L and D 2S isoforms, which are 444 and 415 amino acids in length, respectively [3]. These isoforms share comparable pharmacological features and are expressed in the same cell types, with a ratio that usually favors the expression of the longer isoform [3]. The D 2S , however, is dominant in the cell bodies and projection axons of the dopaminergic cell groups of the mesencephalon and hypothalamus, while the D 2L is more strongly expressed by neurons in the striatum and nucleus accumbens, brain structures targeted by dopaminergic fibers [4]. The D 2L isoform differs from D 2S by the insertion of 29 amino acids in the ICL3 of the receptor. This loop is engaged in the coupling of the receptor to various intracellular partners. The available data show that the D 2 isoforms have differentiated affinities to different G proteins, implying that these receptors might have different roles in vivo [3]. For instance, it was found that dopamine D 2S and D 2L receptors may differentially contribute to the actions of antipsychotic and psychotic agents in mice [5].
In 1993, Montmayeur et al. showed that the 29 amino acids insertion present in D 2L allows it to interact specifically with G i2 [6]. Moreover, in vivo studies indicate that the D 2L isoform is preferentially coupled to G i2 and D 2S prefers G i1 over G i2 [7]. In 1994, Sengoles performed a study with the mutant G iα proteins and found that the D 2L isoform signaled through the G i3 protein and D 2S through G i2 [8]. In 2004, Sengoles et al. created two-point mutations of the D 2S receptor in ICL3 by random mutagenesis (R233G and A234T) [9]. The mutant receptors exhibited a change in the G i subtype coupling preference in comparison to a native D 2S receptor, suggesting the importance of the ICL3 sequence for G protein interaction. Therefore, while available data gives some premises, details of preferential coupling of the isoforms with G i protein subtypes remain unclear.
The interaction of D 2S or D 2L isoforms with an agonist results in the process of receptor activation and subsequent binding with G protein. The outward movement of TM6 is the most pronounced conformational change on the cytoplasmic side of the receptor occurring during activation as reported for β 2 adrenergic receptor [10]. This change is assisted by the outward movement of TM5 and a slight inward adjustment in the position of TM3 and TM7 to accommodate space for the binding of G protein [10]. The mechanisms of activation are, however, common in the whole family A of GPCRs.
The process of signal transmission from the agonist binding site on the extracellular part of the receptor to its intracellular part is enabled by molecular switches [11]. There are four widely approved molecular switches believed to take part in activation of class A GPCRs: two switches involving movements of specific side chains including (1) the W6.48 tryptophan toggle switch (the CWxP motif in TM6) accompanied by a transmission switch involving neighboring amino acids and (2) the Y7.53 tyrosine toggle switch (the NPxxY motif in TM7); as well as two other switches operating by breaking of bonds/interactions linking TMs: (3) the ionic lock involving helices TM3 and TM6, and (4) the 3-7 lock linking helices TM3 and TM7 (residue numbers according to Ballesteros-Weinstein nomenclature [12]) [11,13].
In spite of about 350 GPCR crystal structures including 85 active structures [14] in the Protein Data Bank, the full experimentally solved structures of receptors with long and flexible ICL3 are not yet available. Due to limitations of experimental methods, at present, such structures can be only obtained using molecular modeling approaches as it has been recently done, e.g., for the human serotonin 5-HT 2A receptor [15]. In this work, we aimed to construct complete models of D 2S and D 2L isoforms in complex with dopamine and both G i1 and G i2 proteins and, subsequently, to examine possible G protein coupling preferences of these isoforms and the receptor activation/deactivation hallmarks upon agonist binding and G protein coupling with molecular dynamics (MD) techniques.

Construction of the Models of D 2L and D 2S Receptors in Complex with G i1 and G i2 Proteins
First, we constructed models of the studied receptors in active conformation: D 2L receptor in complex with G i1 protein (L1), D 2L receptor in complex with G i2 protein (L2), D 2S receptor in complex with G i1 protein (S1), and D 2S receptor in complex with G i2 protein (S2) using homology modeling with Modeller v.9.19 and Yasara tool for loop modeling. The attempts to model the ILC3 loop with Robetta [16] or I-Tasser [17] were not successful. The crystal structure of human β 2 -adrenergic receptor in complex with a heterotrimeric G s protein (PDB ID: 3SN6 [18]) was used as a template for the helical bundle and G protein while the crystal structures of human dopamine D 2 , D 3 , and D 4 receptors (PDB IDs: 6CM4 [19], 3PBL [20], and 5WIU [21], respectively) were used as templates for the extracellular loops. The sequence identity between the dopamine D 2 receptor and β 2 adrenergic receptor is 37% and the sequence similarity 59%. The final models were selected from the generated populations of 100 models each (four final models from each population, 16 models in total). Figure 1 shows best-scored models of D 2S and D 2L receptors with the ICL3 loop. The transmembrane region of the receptors displays the features of the active conformation. The ICL3 loop connects the intracellular termini of TM5 and TM6. As predicted by the Predict Protein [22] online server, both the short and long variant of the loop are constituted by a few α-helical regions connected by disordered regions. The best-scored receptor-G protein complexes were validated using ProCheck [23], Verify3D [24], and ERRAT [25], as well as by molecular docking. The results of validation of the 16 best-scored models are presented in Table 1. It can be seen that, according to the used tools, the obtained models are characterized by good quality. According to ProCheck, over 90% of residues in all the models fall into the favored region and only about 1% of amino acids are located in the outlier region. As visualized by the Ramachandran plots, no outliers were found in the transmembrane regions of the receptor models. Most outliers were glycine and proline residues located in loop regions of G proteins. Moreover, according to Verify3D, the final models had over 80%of residues that had an average 3D-1D scores ≥0.2 and, according to ERRAT, the overall quality factor in all chosen models exceeded 90.

Ligand Docking
A natural agonist-dopamine-was docked with Molegro software to sixteen best-scored receptor-G protein complexes models. The final docking poses of the ligand to four models ( Figure 2) were selected based on the scoring functions values, visual inspection and literature data. As typical for the orthosteric ligands of aminergic GPCRs, the protonatable nitrogen atom of dopamine forms a salt bridge with the carboxylic group of the D3.32 in TM3. Moreover, meta-and para hydroxyl groups of the catechol ring form hydrogen bonds with the side chains of S5.43 and S5.46 in TM5, respectively. In the case of L1 complex, an additional hydrogen bond interaction of the para-hydroxy group of dopamine with the main chain of V3.33 was found. In the case of L1 and S1 complexes, the hydrogen bond was also observed between the protonatable nitrogen atom of the ligand and the side chain of T7.38 or Y7.34, respectively.

Molecular Dynamics Simulations
For the preliminary MD simulations, 16 best-scored models (four from each population) in complex with dopamine were selected. Each ligand-receptor complex was immersed in a native-like membrane and subjected to minimization, equilibration, and production MD simulations of 200 ns. RMSD (root-mean-square deviation), RMSF (root-mean-square-fluctuation), and the protein-ligand short-range Lennard-Jones and Coulombic interaction energy analyses were performed. Figure 3 shows plots of four selected systems that were subjected to further processing. RMSD and RMSF values for the C α atoms were calculated using Gromacs tools for preliminary 200 ns simulations to check for the stability of the models. The average RMSD values were: L1-4.76 Å, L2-4.09 Å, S1-4.97 Å, S2-4.56 Å, and RMSF: 1.95 Å, 1.80 Å, 1.79 Å, 2.00 Å, respectively. As expected, the ICL3 regions are characterized by the highest RMSF values compared to all other regions and they are the most flexible parts of the systems. The selected complexes were subjected to the production phase of molecular dynamics (1 µs) in three replicas each. Moreover, we simulated four corresponding systems without dopamine for comparison.  Figure 4 shows plots of the center-of-mass (COM) distances between the C α atoms of the receptor residues: Y7.53 in highly conserved NPxxY motif, L2.46, I3.46, T5.54, V6.40, and the C α atoms of the α5-G α (C-terminus) G protein residues: D351, C352, G353, L354, F355 (in both G i protein subtypes) during 1 µs MD simulations. In order to illustrate the tendency to decrease or increase the distance, the best fit line of all points was drawn. Figure 4 clearly shows that the distance between the studied residues decreases for the L2 and S1complexes. In contrast, for the L1 and S2 complexes, the studied distance increases. As the validation of these results, the same distance was measured for the systems without dopamine ( Figure 5) and it increased in all cases. These results suggest the preferential coupling of the D 2L isoform with G i2 protein and D 2S with G i1.   Figure 6 presents plots for α5-G α protein fluctuations from the entire simulation. The average values for S2 indicate greater fluctuations in this region compared to other complexes. In the case of L2, the average fluctuation values show that the studied region is more stable during 1 µs molecular dynamics simulations. To approximate the interaction between G protein and D 2 receptor, the energy of interaction was calculated with the Gromacs tools. Table 2 shows the values of average short-range Lennard-Jones potential energy for all replicas and receptors without dopamine.  Furthermore, in order to investigate the interactions between G proteins and dopamine D 2 receptor isoforms, the distance matrices consisting of the distances between residue pairs were calculated with gromacs tools ( Supplementary Information, Figure S1). The last 200 ns of simulations were taken into account. On the distance map, points corresponding to the length from 0-3 Å are marked as black dots. Interactions between α5-G α protein and D 2 receptor residues are surrounded by a circle. The role of M6.36 is worth noting as it appears on all the maps except for simulations without agonist and all S2 simulations. In the case of S2 F7.56 appears in all simulations with the agonist. R6.58 appears in all simulations of L2 and S1. Moreover, the role of P139 of ICL2 is also apparent as it does not appear only for S1. In this case, R3.50 of DRY motif seems to play a role. It can be concluded, that the performed analysis did not allow to indicate the residues from dopamine D 2 receptor isoforms which are responsible for the selectivity of interactions with respective G i protein subtypes as all the indicated residues are shared by both isoforms. Probably, the positions and/or conformations of identified residues is different in four investigated complexes which can explain their differentiated importance for interactions with G i protein subtypes. In the case of L2 in the second replica with agonist, L2 without agonist and S1 without agonist, significant flexibility can be observed. Interestingly, for the χ1 of H6.55 remarkable differences were found. For L1 and S2 in all simulations with and without agonist, the dihedral shows large oscillations for most of the simulation time. In simulations of L2 and S1, this conformational change does not occur so often, although those regions are flexible, and the dihedral angle oscillated at around 180 • . One of the hallmarks of GPCR active conformation is described by the transmission switch (former Trp rotamer toggle switch) [11,13]. In short, in all GPCR crystal structures with agonists, movements of TM5 and TM6 can be observed, including relocation of conserved residues Trp6.48 and Phe6.44 toward Pro5.50. It results in the bending of TM6 at the CWxP motif in active structures in contrast to straightened TM6 in inactive structures. Figure 12 shows TM6 conformational changes for the receptors with and without dopamine after 1 µs molecular dynamics simulations. L1 with an agonist showed slight differences between the starting conformation and the conformation after MD simulations. The difference appeared below the rotamer toggle switch, W6.48, where a slight bending in the two out of three replicas was observed. In the second replica and for the receptor without agonist, a tendency to straighten TM6 is evident. An interesting observation was found for L2 and S2 complexes. In the first case, significant helical bending can be seen below the W6.48 for all simulations with dopamine. For S2, all simulations showed a tendency to straighten TM6 and shift M6.36 into the protein interior. For S1, only in the one replica and the receptor without agonist TM6 clearly changed the spatial conformation. MD simulations also showed that, for receptors without dopamine, the tendency to straighten TM6 is notable. These observations further confirm the preferential coupling of the D 2S isoform with G i1 protein while D 2L isoform can be coupled mainly to the G i2 protein and also to the G i1 protein.     The distance between TM5 and TM7 at Cα atoms of Y5.58 and Y7.53 was also examined ( Figure 13). Table 3 shows the distances for receptors after 1 µs MD simulations. It is clearly seen that the distances between Y5.58 and Y7.53 for L1 and S2 were significantly higher than for L2 and S1. In comparison to the inactive state D 2 crystal structure (PDB: 6CM4 [19]: 19.4 Å), the Y5.58-Y7.53 distance in the S2 structure showed the highest values, indicating a structural similarity with the inactive state of the receptor, while L2 showed the smallest distance value in all simulations further confirming preferential G i subtype protein coupling.

Structural Dynamics and Conformational Changes of Helices and Switches
It is also well-known that activation of GPCRs correlates with the formation of a continuous internal water pathway as a hydrophobic layer of amino acid residues next to the characteristic NPxxY motif forms a gate that opens to form a continuous water channel only upon receptor activation [26]. Rearrangement of TM3, TM6, and TM7 affects the formation of the water channel visible in Figures 14 and 15. The water channel remained open for all replicas of L2 with dopamine, one replica of L1 with dopamine, and two replicas of S1 with dopamine. For S2, in all replicas, water molecules were blocked at the hydrophobic barrier consisting of five residues: L2.42, I2.43, L3.43, I3.46, and M6.36 (orange in Figure 7). This situation also happened in the case of L1 in two replicas where TM5, TM6, and TM7 conformation changes were not clearly visible. This data is consistent with the results described above. Surprisingly, for L2 and S1 without dopamine, the water channel remained open during the MD simulations as well as for L2 with dopamine for all replicas and S1 in two replicas, which indicates a slower deactivation compared to L1 and S2. In simulations without dopamine, sodium ion (purple in Figure 15) diffused from the solvent layer into the deep orthosteric pocket of the receptor, which is consistent with the results described above.   Figure 14. Internal water channel with marked residues (orange) forming a hydrophobic barrier for receptors with dopamine after 1 µs MD simulations. The helices are marked as pink springs and the sodium ions as a purple sphere.
Another important feature that distinguished receptor conformational changes after 1 µs MD simulations was the distance between the R3.50 of the DRY motif and E6.30, forming the interhelical ionic lock involving TM3 and TM6. Because of the fact that receptor models were constructed on the active 3SN6 template, the distance values, given in Table 4, are much higher than in the inactive crystal structure of the dopamine D 2 receptor (3.0 Å) showed in Figure 16. However, differences in distances are notable. The smallest distance appeared in all S2 replicas. All receptors without dopamine showed lower values than their counterparts with dopamine. The highest values(>20 Å) appeared for all replicas of L2 and S1 supporting our previous observations.  Referring to the literature data indicating that R233 from the ICL3 of D 2S isoform is important for interaction with G protein [8], we also investigated the interactions between the α subunit of G protein and R233 in ICL3 of both isoforms. Table 5 presents the average distances between the side chain COM of D309, L310, and K312 in G iα1 or K307, D310, and L311 in G iα2 and side chain COM R233 in ICL3 during the last 100 ns of simulations. In the case of L2 and S2, the smallest distance is observed. Moreover, in two replicas of L2 and one replica of S2, the hydrogen bond in the last step of the simulation is formed between D310 from the G protein and R233 from the ICL3 of the receptor.

Principal Component Analysis
For additional support, MD simulations were analyzed by principal component analysis (PCA)for all the trajectories concatenated in four populations separately for L1, L2, S1, and S2 with their counterparts without dopamine. The conformational changes of TM6 and TM5 in combination with TM7 were considered by non-mass-weightened analysis. The projections were for the last 100 ns of the MD simulations. Covariance matrices of individual subspaces were compared, using the Gromacs tools, and are given in Table 6. The overlap is 1 if matrices are identical. It is 0 when the sampled subspaces are totally orthogonal [27]. Our result allows the comparison of subspaces with a slight deviation only.  Figure 17 shows the conformational space of TM6 explored in MD simulations, described by PCA projection of the motion of the TM6 (fragments covering the last 100 ns of simulations are presented to improve clarity). Clusters for receptors without agonists are negative along PC1, which corresponds to straightening the TM6 at W6.48. For L1 and S2, there are significant deviations from the start conformation. L2 shows a shift toward positive values for PC1 (bending TM6 at W6.48) while S2 takes negative values. We can see that TM6 of S2 adopts a similar conformation state as TM6 for the receptor without an agonist. The changes for L1 and S1 show slight shifts compared to the starting structure, which is also shown in Figure 12 with the final TM6 conformations. Figure 18 shows the motion of the protein in phase space along PC1 and PC2 values for TM5 in combination with TM7 at Y5.58 and Y7.53. Clusters for receptors without an agonist for L1, L2, and S2 are negative along PC1. It is significant that clusters of S2 with agonist have mostly negative values, which means adopting conformation similar to the receptor without dopamine. The changes for S1 show slight shifts compared to the starting structure. Thus, PCA results support our conclusions described above.

Discussion
The aim of our work was to construct full models of the dopamine D 2 receptor D 2S and D 2L isoforms in complex with a natural agonist, dopamine, and to study the coupling of these isoforms with G i1 and G i2 proteins, as the experimental data about D 2S and D 2L isoforms and G i protein subtype preference remains unclear [6,8]. Although models of the dopamine D 2 receptor in active conformation with or without the respective G protein are already available in the literature [28][29][30][31], this is, to our best knowledge, the first time full D 2S and D 2L isoforms, including ICL3 loop, have been modeled. Modeling of ICL3 in the case of GPCRs with a long ICL3 is a challenge as there are no templates for this highly flexible protein fragment. However, intracellular loops, in particular ICL3, play an important role in receptor activation due to their interaction with G proteins and the influence on G protein preferential coupling [32][33][34]. Recently, a full model of human serotonin 5-HT 2A receptor has been reported as an example of this challenging loop modeling [15].
Before any further investigations could be performed, an in silico model needs to be properly validated. The protocol for the construction of a homology model of the dopamine D 2 receptor in active conformation in complex with G protein using multiple templates is based on our earlier experience with the modeling of the µ opioid receptor in active conformation [35][36][37][38]. This protocol turned out to be successful as the comparison of the µ opioid receptor model with the crystal structure revealed the C α RMSD of whole structures of 2.60 Å, while removal of the most disordered regions (N-terminus, C-terminus, ICL3) decreased RMSD to 1.91 Å, below the crystal resolution (2.10 Å) [37]. The short and long ICL3 loop were modeled based on the predicted secondary structure and the interactions of these loops with respective regions of G proteins are in accordance with experimental data [9]. The final models were selected from populations of 100 models each and assessed using internal Modeller scoring functions as well as validated applying widely used tools for this purpose. According to all the studied factors, the models were characterized by good quality. The docking of dopamine to the models, in agreement with available literature data, further confirms the correctness of the models. Next, we performed preliminary molecular dynamics studies of 16 models (four models of each type) to select, after 200 ns of simulations, the best-scored four systems for further investigations. Finally, our simulations were performed in a native-like membrane environment, which allows us to reproduce correctly molecular events concerning GPCR functioning and to follow subtle aspects of these receptors' early activation at the molecular level [39,40].
As there are no reports about structural aspects of the full-length dopamine D 2 receptor concerning its interactions with respective G proteins, we investigated how the length of the ICL3 loop affects the interactions with the G i1 and G i2 proteins and receptor activation processes. It turned out that differences in loop length contribute to the different behavior of the receptor when it binds to a specific G protein. The increasing distance between the intracellular parts of dopamine receptor D 2L isoform transmembrane helices, and α5-G i1 , the binding region of G protein, indicates that G i1 protein is sliding out of the receptor-binding surface and the signal transduction started by agonist binding is stopped. The corresponding distance for the D 2L and G i2 complexes decreases slightly, which may indicate that the active state is maintained. Simulations without dopamine showed that the described protein regions move away from each other, which may also indicate the deactivation process and which constitutes an additional validation of our simulations. In the simulations for D 2S isoform, it can be seen that the D 2S receptor moves away from the α5-G i2 protein. The decrease in the distance between α5-G α helix and the intracellular parts of receptor transmembrane helices for L2 and S1 complexes in all replicas suggest preferential coupling of the D 2L isoform with G i2 protein and D 2S isoform with G i1 protein, which is in agreement with experimental data obtained by Montmayeur et al. [6] and Grünewald et al. [7].Thus, the length of the ICL3 loop of dopamine D 2 receptor isoforms governs coupling with the respective G i protein subtypes, which can be of importance to explain different in vivo roles of these isoforms. Experimental studies confirmed the importance of ICL3 loop N-and C-termini for G protein coupling [41,42].
Many studies have shown that GPCRs share a set of residues called molecular switches that are involved in receptor activation and signal transduction [26,43]. The rearrangement of the hydrophobic residues in TM2, TM3, and TM6 helices are involved in signal propagation [44]. In the inactive state of dopamine, D 2 receptor residues L2.42, I2.43, L3.43, I3.46, and M6.36 (orange on Figure 7) constitute a hydrophobic pocket which prevents the creation of hydrogen bond network involved in the interaction with G protein. Our molecular dynamics simulations show that in the case of systems where the distance to respective G proteins is reduced, the hydrophobic barrier is broken and the flow of water molecules in the water channel is enabled, which is characteristic for the active state of GPCRs [45,46]. Interestingly, some systems without dopamine also maintained an open hydrophobic barrier, which may indicate a slow receptor deactivation, and further support the assumption that breakdown of the intra-receptor water chain in many dopamine-bound L1 and S2 complexes was caused rather by non-compatible allosteric signals from binding partners than by chance. Moreover, PCA results show that one of those systems-S1-maintains many similarities to agonist-bound systems, which may indicate that in this particular combination the allosteric effect of G protein coupling on the receptor conformation is considerable. We also examined the conformational changes of TM6 and found that its bending is correlated with the distance between TM5 and TM7 in the place of two highly conserved tyrosine residues Y5.58 and Y7.53.Systems with the tendency to straighten TM6 were also characterized by the largest distances between TM5 and TM7.The stronger bending leads to closer interaction between Y5.58 and Y7.53, which was suggested to stabilize the active state of the receptor via the water-mediated hydrogen bond [47,48].Our results regarding the conformation of TM6 further support preferential coupling of dopamine receptor isoforms with respective G i protein subtypes.
We also examined the behavior of several microswitches, which are important for the GPCR activation process. We analyzed changes in the dihedral angle of the conserved W6.48 [43,[49][50][51], F6.44 (called transmission switch [52][53][54]), Y7.53 of NPxxY motif [32,55], and H6.55 (a crucial residue for dopamine D 2 receptor activation [56,57]). It was found that H6.55 is characterized by a high degree of conformational flexibility. We also investigated the behavior of the residues forming ionic lock between R3.50 and E6.30 and noticed the correlation between the ionic lock distance and other hallmarks of receptor active or inactive state. The lowest values were observed for complexes without dopamine and for the D 2S G i2 , so we can assume that these models tend to adopt inactive conformation of the receptor.

Receptor Model Construction
The sequences of the human D 2L and D 2S receptors, G αi1 and G αi2 proteins were obtained in FASTA format [58] from the UniProt database (https://www.uniprot.org/).The crystal structure of the human β 2 adrenergic receptor complexed with a heterotrimeric G s protein (PDB ID: 3SN6 [18]) was used as a template for the helix bundle for homology modelingof D 2L and D 2S receptor isoforms in the active conformation as well as the template for G protein. In addition, the crystal structures of human dopamine D 2 , D 3 , and D 4 receptors (PDB IDs: 6CM4 [19], 3PBL [20], and 5WIU [21], respectively) were used as templates for the extracellular loops. Multiple sequence alignment was carried out with MUSCLE (Multiple Sequence Comparison by Long-Expectation) [59].
The homology models of D 2L and D 2S receptors in active conformation in complex with the respective G proteins were built using Modeller v. 9.19 (AndrejŜali, San Francisco, USA) [60]. The models of the D 2L and D 2S ICL3 loops were generated with Yasara software with restrictions imposed on secondary structure predicted by PredictProtein [22] online server (https://www.predictprotein.org/). Ten models were created in each of the four model populations (L1, L2, S1, S2) differing in the loop conformation. The most probable loop models for the long and short loops were selected on the basis of their potential interactions with G proteins and favorable orientation in the intracellular area (with no overlap with the membrane).

Molecular Docking
The structure of orthosteric ligand, dopamine, was modeled using the Hartree-Fock approach and 6-31G* basis set of Spartan v. 10 VI.0.1(Wavefunction, Inc., Irvine, California, USA) [61]. MolegroVirtual Docker 6.0 software (Molexus IVS, Odder, Denmark) [62] was used for docking simulations of flexible ligand dopamine into the rigid receptor models. The actual docking simulations were performed using the following settings: number of runs = 100; maximal number of iterations = 10,000; maximal number of poses = 50; and the poses representing the lowest value of the scoring function (MolDockScore) were further analyzed as previously reported [63,64]. The most probable docking pose of dopamine was selected from the poses where a protonatable nitrogen atom of dopamine formed an electrostatic interaction with the conserved aspartate from the third transmembrane helix, D3.32, taking into account available literature data on other interactions between dopamine and D 2 receptor [65].

Molecular Dynamics
Sixteen final ligand-receptor complexes were subjected to molecular dynamics with Gromacs version 2018.4 [66] in native-like conditions. The membrane environment for the complexes was prepared using the Charmm-GUI Membrane Builder server [67]. The complexes were immersed in an asymmetric membrane consisting of nine types of lipids in the proportions appropriate for membrane rafts [68] containing cholesterol, sphingomyelin, DOPE, DOPC, DOPS, PLPC, POPC, POPE, POPG, and aqueous phase: TIP3P water molecules with 0.15 M NaCl. A 3SN6 crystal structure was used as the receptor orientation template in a membrane. An Amber03 force field [69] was used for receptors, Slipids (Stockholm lipids) [70] for the membrane, and General Amber Force Field (GAFF) [71] for ligands. EPS charges were obtained by RESP ESP charge Derive Server [72] and processed with the ACPYPE server [73] to gain ligand topologies. The properly protonated receptor structures were obtained from the H++ server [74]. A template receptor (β 2 adrenergic receptor in complex with G s protein) in a membrane was first minimized using 500 steps. Then, it was equilibrated in 1 ns NPT simulations using the Berendsen barostat to control volume fluctuations, followed by 5 ns NVT simulations and 10 ns NPT simulations using the Parrinello-Rahman barostat. The template receptor in the equilibrated system was changed to 16 dopamine-receptor (4 × L1, 4 × L2, 4 × S1, and 4 × S2). Each system was again minimized, and equilibrated under 1 ns NVT and 5 ns NPT as above with protein and ligand position restrains (force constants of 10,000 kJ mol −1 nm −2 ) on the heavy atoms. The molecular dynamics simulations of 16 systems were performed for 200 ns using a time-step of 2 fs. Four systems (one of each type: L1, L2, S1, and S2) were selected based on the comparison of RMSD, RMSF, and energy of dopamine-receptor interactions (the protein-ligand short-range Lennard-Jones and Coulombic interaction energy). These systems were subjected to the production phase of molecular dynamics (1 µs) in three replicas each. As a reference, these systems without dopamine were also simulated. Standard Gromacs tools were used for the analysis of the results.

Conclusions
In summary, we used for the first time in silico approaches to construct models of full D 2S and D 2L dopamine D 2 receptor isoforms and studied their coupling with G i protein subtypes. Our results indicate the preferential coupling of the D 2L isoform with G i2 protein and D 2S isoform with G i1 protein, which is in accordance with experimental data. The G i protein subtype preference is further supported by different hallmarks of receptor active state, including conformation of microswitches, the conformation of TM6, and the formation of the water channel. The results in this study give mechanistic insight in the interplay between isoforms of dopamine D 2 receptors and G i proteins subtypes, which is important to understand signaling by these receptors and their mediation by pharmaceuticals-in particular, psychotic and antipsychotic agents.

Conflicts of Interest:
The authors declare no conflicts of interest. The funders had no role in the design of the study; in the collection, analyses, or interpretation of data; in the writing of the manuscript, or in the decision to publish the results.