A Collision Coupling Model Governs the Activation of Neuronal GIRK1/2 Channels by Muscarinic-2 Receptors

The G protein-activated Inwardly Rectifying K+-channel (GIRK) modulates heart rate and neuronal excitability. Following G-Protein Coupled Receptor (GPCR)-mediated activation of heterotrimeric G proteins (Gαβγ), opening of the channel is obtained by direct binding of Gβγ subunits. Interestingly, GIRKs are solely activated by Gβγ subunits released from Gαi/o-coupled GPCRs, despite the fact that all receptor types, for instance Gαq-coupled, are also able to provide Gβγ subunits. It is proposed that this specificity and fast kinetics of activation stem from pre-coupling (or pre-assembly) of proteins within this signaling cascade. However, many studies, including our own, point towards a diffusion-limited mechanism, namely collision coupling. Here, we set out to address this long-standing question by combining electrophysiology, imaging, and mathematical modeling. Muscarinic-2 receptors (M2R) and neuronal GIRK1/2 channels were coexpressed in Xenopus laevis oocytes, where we monitored protein surface expression, current amplitude, and activation kinetics. Densities of expressed M2R were assessed using a fluorescently labeled GIRK channel as a molecular ruler. We then incorporated our results, along with available kinetic data reported for the G-protein cycle and for GIRK1/2 activation, to generate a comprehensive mathematical model for the M2R-G-protein-GIRK1/2 signaling cascade. We find that, without assuming any irreversible interactions, our collision coupling kinetic model faithfully reproduces the rate of channel activation, the changes in agonist-evoked currents and the acceleration of channel activation by increased receptor densities.

INTRODUCTION GIRK (G protein-activated Inwardly Rectifying K + -channel) channels play fundamental physiological roles, such as control of heart rate and neuronal excitability, and are also highly implicated in disease such as addiction, depression, bipolar disorder, and cardiac arrhythmias (Dascal, 1997;Hibino et al., 2010;Luscher and Slesinger, 2010;Voigt et al., 2014;Mayfield et al., 2015;Rifkin et al., 2017). Opening of the channel is achieved by a highly recognized signaling cascade involving agonist binding to a G-protein Coupled Receptor (GPCR), which in turn activates heterotrimeric G-proteins (of the Ga i/o -type), to promote 'release' of Gbg subunits. Then, direct interaction of Gbg subunits with GIRK leads to channel opening and the appearance of the agonist-evoked current, I evoked [(Logothetis et al., 1987;Reuveny et al., 1994;Lim et al., 1995;Slesinger et al., 1995); reviewed in (Clapham and Neer, 1997;Dascal and Kahanovitch, 2015)].
Despite more than four decades of studies, the details behind this prototypical activation scheme remain highly debated. For instance, it remains unclear how the different signaling proteins are arranged at the membrane to bring about robust and efficient channel opening, fast activation kinetics and, importantly, signaling specificity (the strong preference for GIRK activation by Gbg derived from G i/o rather than G s or G q ). Two contrasting mechanisms have been proposed [reviewed in (Hein and Bunemann, 2009)]. The first, denoted collision coupling (Tolkovsky and Levitzki, 1981;Tolkovsky et al., 1982;Shea et al., 2000), assumes unrestricted diffusion of GPCRs, G-proteins, and effectors in the plasma membrane. After GPCR activation, GIRK activation occurs through random collisions with proteins of this cascade (Vorobiov et al., 2000;Hein et al., 2005;Touhara and MacKinnon, 2018). According to this model, a single receptor may activate several G-proteins (as in visual system; (Arshavsky et al., 2002)) and, therefore, an increase in the number of receptors is expected to accelerate activation kinetics of I evoked . Indeed, this was observed for several GPCR-GIRK cascades (Vorobiov et al., 2000;Wellner-Kienitz et al., 2000;Hein et al., 2005;Kahanovitch et al., 2017). The second mechanism posits long-lived "preformed" complexes of GPCRs, G-proteins, regulatory proteins (e.g., Regulators of G Protein Signaling; RGS) and GIRKs in various combinations (Huang et al., 1995;Lavine et al., 2002;Fowler et al., 2006;Jaen and Doupnik, 2006;Dupre et al., 2007;Rusinova et al., 2007;Dupre et al., 2009;Nagi and Pineyro, 2014;Tateyama and Kubo, 2018). In support, several lines of evidence indicate that Gbg and Ga i/o can associate with GIRKs, as early as in the endoplasmic reticulum ( (Rebois et al., 2006;Robitaille et al., 2009); reviewed in (Zylbergold et al., 2010)), recruit G-proteins to the plasma membrane Kahanovitch et al., 2014), and possibly remain associated at the plasma membrane (Fowler et al., 2006;Rebois et al., 2006;Riven et al., 2006;Berlin et al., 2011;Kano et al., 2019). Further, in neurons, cardiomyocytes and heterologous expression models, some GIRKs exhibit an agonistindependent (basal) current (I basal ) that is Gbg-dependent, suggesting some kind of pre-coupling of GIRK with the G protein before the receptor has been engaged (reviewed in (Dascal and Kahanovitch, 2015)). We have also shown that GIRK1/2, but not the GIRK2 homotetramers, recruit Gproteins to the plasma membrane, favoring Gbg over Ga Kahanovitch et al., 2014). The preferential association with Gbg explains the high, Gbg-dependent I basal of GIRK1/2, contrasting the small and Gbg-independent I basal of GIRK2 homotetramers (Rubinstein et al., 2009;Yakubovich et al., 2015). Together, these findings support the existence of dynamic G protein-GIRK complexes; however, whether they require permanent association is debated (Yakubovich et al., 2015). It is conceivable that different signaling cascades may proceed at different modes, namely collision or preformed modes, or a mixture of the two (e.g., only G protein and effector are pre-coupled).
The preformed complex model can seamlessly account for specificity (i.e., preferential activation of the channel by Gbg released from a particular type of Ga-subunit) as well as for the high speed of GPCR-induced activation of GIRKs (Hille, 1992;Huang et al., 1995), limited mainly by the kinetics of G protein cycle (Vorobiov et al., 2000;Lohse et al., 2008;Hein and Bunemann, 2009). However, specificity can also be quantitatively described in purely kinetic terms, i.e. collision coupling (Touhara and MacKinnon, 2018). For instance, if a particular Ga-type, namely Ga i/o , is quicker to provide Gbg to the channel than other Ga's, it would appear as though the channel solely responds to Gbg derived from that specific Ga (Touhara and MacKinnon, 2018). Indeed, heterologous overexpression of proteins of the b-adrenergic-Ga s cascade can lead to activation of GIRK via Ga s -derived Gbg (Lim et al., 1995;Bender et al., 1998;Touhara and MacKinnon, 2018). These results show that the system can indeed proceed, at least in some instances, via a collision coupling mechanism. Lastly, at high levels of expression of the reactants (proteins participating in the cascade), especially the GPCR, the kinetics and magnitude of effector activation via a collision coupling cascade would be indistinguishable from those attained by a preformed complex (Lauffenburger and Linderman, 1996). Therefore, the nature and concentrations of the reactants strongly affect the speed and specificity of the responses, thereby making it hard to distinguish between different modes of activation. It is therefore critical to study these mechanisms by systematic "titration" of the interactors.
In the current work, we set out to understand the mode of coupling in the classical M2R-Ga i/o -GIRK cascade, by combining electrophysiological, fluorescence, and biochemical measurements in Xenopus oocytes with kinetic modeling. Specifically, we studied this cascade by systematically varying and quantifying surface densities of proteins involved in it, and monitored outcomes on GIRK activation. Next, we combine a Thomsen-Neubig-like mathematical model of GPCR activation (Thomsen and Neubig, 1989;Zhong et al., 2003) with extant models of GIRK activation by Gbg to quantitatively describe GIRK activation in detail. We find that a collision coupling model faithfully reproduces both the fast activation kinetics of agonist-induced GIRK responses, and their dependence on GPCR surface density.

Ethical Approval and Animals
Experiments were approved by Tel Aviv University Institutional Animal Care and Use Committee (permits M-08-081 and M-13-002). All experiments were performed in accordance with relevant guidelines and regulations. Xenopus laevis female frogs were maintained and operated as described (Dascal and Lotan, 1992). Briefly, frogs were kept in dechlorinated water tanks at 20 ± 2°C on 10 h light/14 h dark cycle, anesthetized in a 0.17% solution of procaine methanesulphonate (MS222), and portions of ovary were removed through a small incision in the abdomen. The incision was sutured, and the animal was held in a separate tank until it had fully recovered from the anesthesia and then returned to post-operational animals' tank. The animals did not show any signs of post-operational distress and were allowed to recover for at least 3 months until the next surgery. Following the final collection of oocytes, after 4 surgeries at most, anesthetized frogs were killed by decapitation and double pithing.

DNA Constructs and mRNA Injection
cDNA constructs of GIRK subunits, Gb 1 , Gg 2 , M2R, YFP-GIRK1 were described in detail in previous publications (see (Tabak et al., 2019) for a detailed list). Fluorescent proteins (CFP and YFP) contained the A207K mutation to prevent dimerization (Zacharias et al., 2002). DNAs of M2R-CFP and M2R-YFP were produced by inserting the PCR product of the human M2R (Lechleiter et al., 1990), flanked by EcoRI on both sides, in pGEM-HJ vector containing CFP (cerulean) or YFP flanked by EcoRI and HindIII, yielding M2R-C/YFP CT . The M2R-Ga i3-C351G tandem cDNA was created by ligating the M2R cDNA sequence in frame with the Ga i3-C351G cDNA, connected via a 6 nucleotide sequence GAATTC (EcoRI restriction site). Thus, the full primary sequences of M2R and Ga i3-C351G are connected by a 2-amino acid linker, Glu-Phe. The DNA of GluR1 L507Y (Stern-Bach et al., 1998) was generously provided by Y. Stern-Bach (Hebrew University). All DNAs were cloned into pGEM-HE, pGEM-HJ or pBS-MXT vectors, which are high expression oocyte vectors containing 5' and 3' untranslated sequences of Xenopus b-globin (Liman et al., 1992), as previously described Berlin et al., 2011). mRNA was transcribed in vitro as described in (Dascal and Lotan, 1992) and precipitated overnight at -20°C with 4 M LiCl. mRNAs were divided into 1 to 2 ml aliquots and stored at -80°C. 50 nl of mRNA were injected into the equatorial part of oocytes, two to three days before the experiments.

Electrophysiology
Whole-cell GIRK currents were measured using two-electrode voltage clamp (TEVC) with Geneclamp 500 (Molecular Devices, Sunnyvale, CA, USA), sampled at 1 kHz and filtered at 200 Hz, at room temperature (21-23°C), as previously described (Rubinstein et al., 2009), using agarose cushion electrodes (Schreibmayer et al., 1994) filled with 3M KCl, with resistances of 0.1-0.6 MW. Data acquisition and analysis were done using pCLAMP software (Molecular Devices). For recording, oocytes were placed in a fast-perfusion chamber (see Figure S1A and below). Holding potential was set to -80 mV. Basal GIRK currents (I basal ) were measured by switching from physiological solution (low K + , ND96: 96 mM NaCl, 2 mM KCl, 1 mM CaCl 2 , 1 mM MgCl 2 , 5 mM Hepes, pH 7.5) to a high K + solution (high K + , in mM: 24 KCl, 2 NaCl, 1 CaCl 2 , 1 MgCl 2 , HEPES, pH 7.5). For recording evoked currents (I evoked ), solution was then switched to high K + solution containing 10 mM acetylcholine (ACh). Addition of 5 mM Ba 2+ (blocker of GIRK) was typically applied at the end of each recording to isolate non-GIRK currents and to calculate net GIRK currents. Total current (I total ) was assessed by summing basal and evoked currents (I total = I basal + I evoked ).
In the perfusion system employed in these experiments, each perfusion tube (inlet) is directly incorporated into the bath chamber (shaped like a thin elongated bar), rather than via a manifold, ending~3 mm from the oocyte in large diameter openings (see Figure S1A). The suction has also been incorporated directly into the bath, 1 mm above the level of the oocyte, to reduce the bath solution volume and to allow fast exchange of the solution. In order to access the solution exchange time, we employed an AMPA receptor (AMPAR) mutant; GluR1 L507Y , which lacks fast desensitization (Stern-Bach et al., 1998). The activation time constant (t act ) of AMPAR activation is below 1 ms (Grosskreutz et al., 2003) and can therefore be considered essentially instantaneous compared to the slower kinetics of GIRK activation. 1 ng RNA of GluR1 L507Y was injected into Xenopus oocytes. 50 nl of 20 mM solution of EGTA was injected 2 hours before experiment, to prevent the appearance of Ca 2+ -dependent Clcurrents (Dascal, 1987). AMPAR was activated by applying saturating glutamate concentration (1 mM) to the bathing solution. A representative recording of AMPAR current is shown in Figure S1B. The rising phase of the response to glutamate, fitted to a single exponential function, was 88.6 ± 14.5 ms (n=9), indicating that the average solution exchange rate time constant of our perfusion system was about 90 ms.

Whole Oocyte Radioactive Quinuclidinyl Benzilate (QNB*) Labeling
Whole oocyte binding experiments were performed as described in (Ben-Chaim et al., 2003). Briefly, three days following mRNAs injection, oocytes were dropped into 200 ml of 0.67 nM QNB* (Halvorsen and Nathanson, 1981). After 1 min of incubation, the oocyte was taken to a washing chamber that contained 4000 ml of ligand-free solution (washing stage) and rapidly (after 1-2 s) removed to the scintillation liquid by use of a custom device. The device is composed of plastic holder that enables insertion of a pipetor with a standard pipette tip (volumes up to 200 ml) trimmed 4 mm from the edge. This ensures extraction of the oocyte with minimal amount of liquid. Then, individual oocytes were placed in vials to which 4 ml of scintillation fluid was added and counted with Packard 2100TR TriCarb Liquid Scintillation Analyzer. Specific binding was determined by subtracting the binding from native, uninjected oocytes.

Fluorescence Imaging
Imaging of fluorescence in the plasma membrane (PM) of whole oocytes was performed as previously described (Berlin et al., 2010;Berlin et al., 2011). Briefly, whole oocytes were placed in a glass-bottom dish, and all images were obtained from the animal hemisphere close to oocyte's equator (see Berlin et al, 2010; Figure 2-micrographs showing homogenous fluorescence in the animal pole). Imaging experiments were performed on a confocal laser scanning microscope (LSM 510 Meta, Zeiss, Germany) with 20x or 5x objective lenses, digital zoom = 2, pinhole 3 Airy units, equipped with a HFT 405/514/633 beam splitter. CFP was excited by 405 nm and emission was collected in the wavelength interval of 449-500 nm, peak emission 481 nm; YFP was excited by 514 nm and emission collected in the interval of 524-609 nm, peak emission: 534 nm. Analysis was done using Zeiss LSM software. The net intensity of fluorescence in the PM was measured by averaging the signal obtained from three standard regions of interest along the membrane with background subtraction (Berlin et al., 2010). In each experiment, uninjected oocytes were tested for intrinsic fluorescence with the use of either lasers: 405 and 514 nm excitation. In all confocal imaging procedures, care was taken to completely avoid saturation of the signal. In each experiment, all oocytes from the different groups were studied using constant LSM settings.

Statistical Analysis
Statistical analysis was performed using SigmaPlot software (Systat Software, San Jose, California, USA). Data are presented as mean ± S.E.M. Two group comparisons were done using two-tailed student's t-test. Multiple group comparison was done using oneway analysis of variance (ANOVA), with post hoc Tukey tests. Activation kinetics (t act ) were obtained by fitting evoked currents by a mono-exponential fit.

Modeling and Simulation
Steady state calculations for estimation of initial conditions were done with Matlab for Windows (Mathworks Inc., Natick, Massachusetts). System of algebraic equations is shown in Supplemental Material 2 and was solved numerically. Timecourse simulation was done utilizing Berkeley Madonna 8.3.23.0 (R. Macey and G. Oster, University of California, Berkeley) for Windows. System of ordinary differential equations was generated based on schemes shown in Figures 4A, B and Figure 5B, and solved numerically by the 4 th order Runge-Kutta method.

Estimation of GIRK-Gbg Interaction Parameters
Reported GIRK-Gbg affinity values span several orders of magnitude depending on estimation method (nM -mM range, (Dascal, 1997)). Considering this we utilized available crystal structure of GIRK2-Gbg (Protein Data Bank number 4KFM, (Whorton and MacKinnon, 2013)) and two structures of GIRK1-Gbg complex generated by homology modeling and proteinprotein docking procedure (Mahajan et al., 2013) for structurebased prediction of protein-protein interaction free energy and affinity. In particular, the above mentioned structures were submitted to PRODIGY server (Xue et al., 2016) that calculated both free energy of interaction and K D values utilizing algorithm which is based on correlation between number of interfacial contacts at the interface of a protein-protein complex and its experimental binding affinity together with properties of the noninteracting surface as described by Kastritis et al. (2011). The k on (association rate constant values) were predicted utilizing Transcomp software (Qin et al., 2011). This software utilizes transient-complex theory developed by Alsallaq and Zhou (2008) for predicting protein-protein association rate constants. Coordinates supplied by structural data supplied in PDB format are used to generate the transient complex and rate constant is calculated based on the electrostatic interaction energy in the transient complex.

Collision Coupling Between M2R and G i/o in GIRK Cascade in Xenopus Oocytes
In previous publications, we presented evidence for catalytic collision coupling between M2R and Ga z (Vorobiov et al., 2000) and GABA B receptors with endogenous or coexpressed Ga i/o (Kahanovitch et al., 2017) in the GPCR-G protein-GIRK cascade reconstituted in Xenopus oocytes. To examine whether this is also the case for the M2R-G i/o -GIRK cascade, and for the following quantitative description and kinetic modeling of the cascade, we first characterized the mode of M2R-G i/o coupling using our previously developed strategy (Vorobiov et al., 2000). Specifically, we initially assessed how changes in GPCR (M2R) A B FIGURE 1 | Increasing expression levels of M2R accelerates the activation of GIRK1/2. (A)a representative GIRK1/2 activation experiment. Oocytes were injected with the following RNAs: GIRK1 and GIRK2 (2 ng/oocyte each) and 500 pg/oocyte M2R. I evoked was elicited by 10 µM ACh. Inset-zoom-in on the activation phase of I evoked (black plot) and a mono-exponential fit (red). (B)kinetics of GIRK1/2 activation. Oocytes expressed a constant amount of GIRK1/2 (2 ng RNA/oocyte), with increasing levels of M2R, in the range 10-1000 pg/oocyte, and t act was determined by monoexponential fitting as shown in A (N=2-7 experiments, n= 13-25 cells). concentration impacts the kinetics of activation (t act ) of I evoked of heterotetrameric GIRK1/2 channels. Here, we used a saturating concentration of acetylcholine (ACh; 10 µM) to activate different densities of ectopic M2R, whereas the signal transduction from GPCR to the channel relied on endogenous G i/o proteins. We employed a fast perfusion system with a time constant of solution exchange below 100 ms (see Methods and Figure S1). Importantly, in neurons and cardiomyocytes the t act of I evoked is in the range 200-700 ms (Pott, 1979;Sodickson and Bean, 1996). Thus, our measurements of t act in the oocyte introduces an overestimation of GIRK1/2 activation kinetics. However, this overestimation is relatively minor, especially at low densities of M2R (see below). Our results show that the increase in the amount of mRNA of M2R per oocyte (i.e., increase in surface density) speeds up the activation of GIRK1/2 ( Figure 1A, red plot-mono-exponential fit from which we extract t act ), with a corresponding decrease in the time constant of activation ( Figure 1B). These results are consistent with those obtained for M2R-G z -GIRK and GABA B R-G i/o -GIRK cascades in this heterologous model (Vorobiov et al., 2000;Kahanovitch et al., 2017).

Estimation of Membrane Protein Density
For detailed kinetic analysis of GIRK1/2 activation, we sought to estimate the densities of proteins involved in our system, explicitly GIRK1/2, M2R and G-protein subunits-Ga and Gbg, by a quantitative approach previously developed in our lab (Yakubovich et al., 2015). Briefly, channel density is typically calculated based on the maximal GIRK1/2 current (i.e., I bg ) measured in oocytes that coexpress Gbg at saturating concentration Channel density was calculated from I total as explained in the text. The correlation between fluorescence and number of YFP-GIRK1 molecules is shown with superimposed linear regression line, extended to origin of coordinates. The regression equation was y=4x, i.e. one channel/µm 2 corresponds to fluorescence intensity of 4 AU. Note that, since each channel has two YFP molecules, the calibration factor in this experiment is: 1 YFP molecule/µm 2 = 2 AU. (E) -estimating the surface density of M2R-YFP, for RNA concentrations of 1, 2 and 5 ng/oocyte. YFP fluorescence, in AU, is shown on the left Y-axis. M2R-YFP surface density (right axis) was calculated using the calibration factor derived from YFP-GIRK1/GIRK2 measurements. The relationship between M2R-YFP RNA dose and M2R-YFP surface density was fitted with linear regression, extended to the origin of coordinates, in the form y=62.5x. (usually 5 ng/oocyte of Gb RNA and 1-2 ng Gg RNA). Under these conditions, the channel's open probability (P o ) is~0.105. Channel density in M2R-expressing oocytes can also be calculated from the total current obtained upon activation with saturating 10 µM dose of ACh (I total ). I total is the sum of agonist-independent GIRK1/2 current, I basal , and the ACh-elicited I evoked . We found that, for GIRK1/2, there is a stable relationship between I total and I bg over a wide range of channel densities such that, on average, I bg = 2I total (Yakubovich et al., 2015). If I bg or I total are known, GIRK1/2 density could be estimated using a modification of the classical equation (Hille, 2002): where i single is the single channel current and N is the number of functional channels in the PM. The channel's surface density is defined as N/S, where S is the surface area of the cell (210 7 µm 2 , deduced from oocyte's capacitance of~200 nF (Dascal, 1987)). Based on data and calculations from (Yakubovich et al., 2015), under the conditions used here (24 mM K + external solution), the surface density for GIRK1/2 or YFP-GIRK1/GIRK2 can be estimated using a simple translation factor: Eq:2 density(channels=mm 2 ) = 0:79 I bg (mA) = 1:58 I total (mA) In most experiments reported here, oocytes were injected with 1-2 ng or GIRK1 and GIRK2 m-RNA each. This generally elicited strong channel activity that corresponds to a "high density" group with an average surface density of~21 GIRK1/ 2 channels/µm 2 (Yakubovich et al., 2015). There is a possibility of formation of functional GIRK2 homotetramers under these experimental conditions. However, the basal current we measured ranged between 3 to 5 µA (see Figure 2A). In oocytes, injection of 1 ng mRNA of GIRK2 gives rise to basal currents of 0.05 -0.2 µA (Rubinstein et al., 2009). Therefore, under our experimental conditions, it is most likely that the predominant channels recorded are indeed GIRK1/2 heterotetramers. Moreover, the preferred stoichiometry of the homologous GIRK1/4 channel is two subunits of GIRK1 and two subunits of GIRK4, rather than GIRK4 homotetramers (Silverman et al., 1996).
Next, we employed YFP-GIRK1/GIRK2 as a "molecular ruler" to translate surface densities of the channel, obtained from current, to fluorescence measurements. Here, YFP density was assumed to be twice that of the channel, since each GIRK1/2 heterotetramer is believed to contain two GIRK1 subunits, by analogy with GIRK1/4 (Silverman et al., 1996;Corey et al., 1998). First, we determined the conditions for optimal channel expression for the calibration procedure. We injected increasing amounts of YFP-GIRK1, GIRK2 and Gbg RNAs and observed a linear relationship between I bg and YFP-GIRK1 fluorescence over the range of channel RNA doses of 0.2-1 ng/ oocyte, in line with the assumption that fluorescence corresponds to functional channels ( Figure S2). Linearity was lost at high RNA doses, suggesting that at high expression levels, some channels are non-functional (possibly not at membrane).
In the experiment shown in Figure 2, we expressed a range of doses of YFP-GIRK1/2 with 1 ng M2R RNA and measured surface levels of channel fluorescence and total GIRK currents in response to ACh (Figures 2A, C, D). We also expressed a range of doses of M2R-YFP and monitored YFP surface levels ( Figures  2B, E). A linear relationship between YFP-GIRK surface density and I total of YFP-GIRK1/GIRK was observed at low doses of RNA (below 1 ng) and this range was used for the estimation of YFP-GIRK1 density for calibration purposes ( Figure 2D). In the same experiment and with identical imaging settings, we measured the PM expression of M2R-YFP at different RNA doses and, using the calibration factor from Figure 2D, calculated the PM density of M2R-YFP. Figure 2E shows that the relationship between the RNA dose and M2R-YFP density was linear in the range 1 -5 ng RNA, yielding receptor densities of~20 to~300 M2R molecules/µm 2 .
There are reports showing that the density of ion channels can be higher at the animal hemisphere or enriched around the injection site, which may effectively increase the density of PM proteins (M2R) in these areas (Robinson, 1979;Lopatin et al., 1998;Machaca and Hartzell, 1998). The assumption of homogeneity of M2R and GIRK distribution is therefore an approximation, which gives very good agreement with experiment. Note that, even if all the receptors and channels are located exclusively in the animal hemisphere, the surface density will only be changed two-fold.

Quantifying the Relationship Between M2R-YFP Surface Density and GIRK1/2 Activation Parameters
To compare the PM expression of different M2R constructs used in this study, we injected a range of RNA doses of wild-type (wt) M2R, M2R-YFP and M2R-CFP and measured the number of QNB binding sites in the PM of intact oocytes utilizing the methodology developed by Ben-Chaim et al. (Ben-Chaim et al., 2003). All three M2R constructs rendered similar number of QNB binding sites ( Figure 3A), showing that they express at similar levels. We could therefore extend the M2R-YFP RNAdensity relationship shown in Figure 2E towards M2Rwt and M2R-CFP.
We next studied the impact of M2R-CFP surface density on activation parameters of YFP-GIRK1/GIRK2. We expressed a range of M2R-CFP receptor densities (with a constant amount of channel's RNA), by injecting 1-15 ng RNA/oocyte, and monitored the PM level of M2R-CFP ( Figure 3B) and YFP-GIRK1 ( Figure 3C) along with currents amplitude and t act of I evoked (Figures 3D-F). The level of YFP-GIRK1/GIRK2 remained unchanged at all doses of M2R RNA except when the receptor was injected at excessively high doses, 25 ng. This yielded a decrease in the PM level of the channel ( Figure 3C), likely due to non-specific competition of RNAs for the same pool of ribosomes (Richter and Smith, 1981) or a trafficking interference. We have, therefore, adjusted the amplitude of I evoked for the change (even if slight) in channel's PM level ( Figure 3F). Together, we find that increase in M2R density (validated by fluorescence) is associated with a sharp rise in both the speed and amplitude of channel activation ( Figures 3D-F). Notably, maximal amplitude of I evoked is obtained at lower PM densities of the receptor than those required to obtain the fastest activation (lowest t act ; compare Figures 3E, F). These observations are in-line with the predictions of the catalytic collision coupling model. Unexpectedly, we note that, though M2R-CFP expresses at equivalent levels as M2Rwt ( Figure 3A), it exhibits slower kinetics at all RNA doses (compare Figures 1D,  3D-F). Therefore, M2R-CFP was solely used to assess surface density of the receptor, but not for quantitative description of the native cascade, where we use M2Rwt.

Modeling: G Protein Cycle Model
For quantifying the M2R-G i/o -GIRK1/2 cascade, we elaborated the Thomsen et al. (1989) model of G-protein cycle by combining it with the ternary complex model developed by De Lean et al. (1980). Of note, a similar approach was used by Falkenburger et al. for the description of another muscarinic receptor and cascade, namely the M1R-G q -phosphoinositide signaling mechanism and regulation of the KCNQ channels (Falkenburger et al., 2010;Hille et al., 2014). A schematic representation of the G-protein cycle model is shown in Figure 4A. List of reactions with corresponding rate constants  is shown in Table 1. To maintain microscopic reversibility, we incorporated into the G-protein cycle model a G-protein dissociation step (i.e. RGa GTP Gbg dissociation, reaction 7) and also a GPCR-independent dissociation step (reaction 10) (Sarvazyan et al., 1998;Yakubovich et al., 2005). It must be emphasized that the rate constants for the latter reaction have been incorporated in the model of Touhara and MacKinnon (2018). Ga GTP Gbg dissociation rate has been reported by Hollins et al. (2009) and association rate could be estimated based on the microscopic reversibility assumption (Table 1, reaction 7). Furthermore, GDP/GTP exchange is split into two reversible reactions. Rate constants of GDP and GTP binding have been determined experimentally (Higashijima et al., 1987;Zhong et al., 2003). We also incorporated GaGbg nucleotide free state in the process which leads from Ga GDP bound state to GTP bound state as proposed by Ross (Ross, 2008). In order to estimate kinetic parameters of reactions of the model that have not been determined experimentally, we used the open-source software tools PRODIGY for the estimation of K D values (Xue et al., 2016) and TRANSCOMP for estimating association rate constants (Qin et al., 2011). Lastly, we scrutinized the crystal structure of M2R with heterotrimeric G-protein (PDB: 6OIK) (Maeda et al., 2019) for deriving the kinetic parameters related to M2R-G protein coupling (Table 1).

Modeling: Channel Activation Model
We made use of our previously described GIRK1/2 gating model (Yakubovich et al., 2015), denoted "graded contribution model" ( Figure 4B). This model is based on the assumption that a channel that is occupied by 1 to 4 Gbg molecules can reach the open conformation, but the contribution of each Gbg-occupied state is different; the higher Gbg occupancy, the higher the contribution to open probability. This model is based on studies of Sadja et al. and Ivanova-Nikolova et al. (Ivanova-Nikolova et al., 1998;Sadja et al., 2002) on a highly homologous channel; GIRK1/4. To estimate the parameters of interaction between GIRK1/2 and Gbg for the graded contribution model, we proceeded in a similar approach as done for the unknown parameters of M2R-G-protein interaction, namely we analyzed two structural models of GIRK-Gbg complex. The first is the crystal structure of GIRK2 in complex with Gbg (Whorton and MacKinnon, 2013) (PDB: 4KFM) and the second is a docking model for the GIRK1-Gbg complex (Mahajan et al., 2013) (termed "best scoring model", or BS). Both models were  Table 1. (C) -comparison of the experimentally observed and predicted AChevoked currents, with GIRK-Gbg interaction parameters from two structural models, 4KFM and BS. (D) -representative analysis of the time course of GIRK1/2 activation. The experimental result (dotted black line) is from an oocyte injected with 0.5 ng M2R RNA. Superimposed are simulated currents according to graded contribution model of GIRK1/2 activation. The experimental parameters in this cell were: I basal = 15.5 µA, I total = 22.8 µA. Estimated channel density was 36 channels/ µm 2 . Initial concentrations of Ga and Gbg available to the channel are: for the case of 4KFM model 3.23 molecules Ga i /channel and 0.45 molecules Gbg/channel; for the case of BS model, 3.24 molecules Ga i /channel and 0.46 molecules Gbg/channel. Each plot represents the recorded or simulated current normalized to its maximum. (E)mean t act values from all experiments with wild-type M2R, superimposed on data obtained from fitting of simulated time-courses according to different models. Each experimental point shows mean value of t act ± SEM from one experiment (n=3-7 oocytes). Simulated time courses were generated for the case described by Yakubovich et al. as "high density group", i.e. I basal =13.36 µA, I total =17.2 µA, n=21 channels/µm 2 . Amounts of available Ga and Gbg molecules per one channel, that are required to obtain the observed I evoked , were calculated using the graded contribution model: with the BS structure-based parameters, 3.65 Gbg and 0.39 Ga molecules/channel; with the for 4kfm structure-based parameters, 3.62 Gbg and 0.38 Ga molecules/channel. subjected to analysis by PRODIGY and TRANSCOMP software from which we obtained K D and k on values of GIRK-Gbg interaction. For calculating initial values for all channel and Gproteins states (i.e. before agonist application), a system of algebraic equations was numerically solved assuming that the reaction between GIRK-G protein is in steady-state (Supplemental Material 1). Time-course of activation of GIRK1/2 was simulated as a solution of system of ordinary differential equations (Supplemental Material, Eq. 19-34).

Simulation of GIRK1/2 Activation Time-Course and Amplitude
In order to validate our model, we simulated the experimentally observed GIRK1/2 activation by a step application of 10 µM ACh. We chose a representative recording in which 0.5 ng/oocyte M2R mRNA was injected. This corresponds to a PM expression of about 31 receptors/µm 2 , according to our calibrations (see Figure  2E). We used available values of Ga and Gbg for the activation of GIRK from its I basal and I total as previously described (Yakubovich et al., 2015) (see Supplemental Material 1). We then simulated the evoked current of GIRK1/2 and compared with experimental values (Figures 4C, D). We find that our developed model satisfactorily reproduces the fast kinetics and the amplitude of I evoked , with kinetic parameters of GIRK-Gbg interaction obtained from both structural models tested (4KFM and BS). Subsequently, we ran time-course simulations over a wide range of M2R receptor densities, 0.1 -100 receptors/µm 2 . The initial parameters of I basal and I total for these simulations were adopted from the "high density group" of GIRK1/2 channel expression obtained by injection of 1-2 ng (21.7 channels/µm 2 , I basal = 13.36 µA and I total = 17.2 µA; (Yakubovich et al., 2015)). In each simulation, t act was extracted from mono-exponential fit of the activation phase of I evoked , and the calculated values of t act versus receptor density were superimposed on the experimentally measured t act obtained from a large number of experiments with M2Rwt ( Figure 4E). We find that that our model satisfactorily predicts the acceleration of activation rate with the increase in GPCR density. It also faithfully reproduces the real kinetics of the receptor, namely it reproduces the fastest kinetics of activation obtained when using high receptor densities.
(*) k on and k off estimated by PRODIGY and TRANSCOMP. (**) k on adjusted, see text.
red plot), whereas the amplitude of I evoked persistently increased with higher RNA dose of M2R-Ga i3-G351G ( Figure 5A, black plot). These results confirm that increase in RNA dose of the tandem is accompanied by an increase in its surface density and, more importantly, indicate that low doses of the tethered receptor-Ga pair cannot provide enough Gbg to activate the large number of channels. We then proceeded to develop a scheme of GPCR-Ga i3 tandem-mediated activation of GIRK1/2 ( Figure 5B) and subsequently simulated time-course of GIRK1/2 activation by a range of GPCR-Ga i3 tandem densities. Three possible scenarios were simulated: 1) M2R-Ga concatemer was assumed to have the same affinity to Gbg as Ga, and no change in Gbg concentration was assumed with concatemer expression; 2) M2R-Ga concatemer was assumed to have same affinity to Gbg as Ga, and 1:1 increase in Gbg concentration was assumed with concatemer expression; 3) M2R-Ga concatemer was assumed to have 10 fold lower affinity to Gbg than Ga and no change in Gbg concentration was assumed with concatemer expression. Simulated evoked currents and t act values show that, whereas I evoked values increase with the increase in GPCR-Ga i3 density, t act values remain constant ( Figures 5C, D).
These simulations thereby fully recapitulate the outcomes on channel activation via a preformed complex between the GPCR and the G protein and provide a unique and contrasting picture than what we obtain when all components are free (i.e., non-fused).
Thus, these observations further support the collision-coupling nature of M2R to G protein signaling in the M2R-GIRK cascade reconstituted in Xenopus oocytes.

Comparison to Cooperative Gating Model
We next tested another detailed model of the GPCR-GIRK cascade based on collision-coupling published by Touhara and MacKinnon (2018), termed here cooperative gating model. This model incorporates the Thomsen et al. model of the G-protein cycle (Thomsen et al., 1988;Thomsen and Neubig, 1989), and the receptor independent G-protein heterotrimer dissociation (Sarvazyan et al., 1998;Sarvazyan et al., 2002). This model assumes that channel activation is cooperative, i.e., each Gbg binds to the channel with stronger affinity than the previous one. It also assumes that GIRK can open only when all four Gbgbinding sites have been occupied. In the model, we applied affinities and rate constants from (Touhara and MacKinnon, 2018),

A B D C
FIGURE 5 | Physical tethering of M2R and Gai3 converts a collision coupling mechanism to a preformed-complex: experiment and simulation. (A) -Incremental expression of fused M2R-Ga i3 C351G (PTX-insensitive) in the presence of coexpressed A-protomer of pertussis toxin (PTX; 0.2 ng RNA/oocyte) shows increase in I evoked with growing amounts of injected RNA (black plot), but kinetics of activation remain unchanged (red plot). Right-Histogram of evoked currents (black) and t act (red) of GIRK1/2 coexpressed with M2R-wt. Result shown are from one experiment; number of cells (n) tested in each group are shown above experimental points or in the bar. No significant change in t act was found (one way ANOVA, P = 0.154). Spearman correlation coefficient calculated for analysis of evoked currents~1, P = 0.0167. (B)scheme of G-protein utilized to simulate GIRK1/2 activation by M2R-Ga i3 C351G. Blue arrows and numbering denote reactions that are shared with M2R wt activation pathway, as described in Figure 4. Red arrows denote reactions present only in the current scheme. The numbering of reactions and the rate constants are the same as in Table 1. (C)simulated I eoked values obtained assuming a range of expression level of fused M2R-Ga i3 C351G. (D)summary of t act obtained from fitting time-course of activation of GIRK1/2 by range of M2R-Ga i3 C351G densities with mono-exponential function. Three possible scenarios were simulated for analysis of M2R-Ga i3C351G experiments. Black bars; M2R-Ga concatemer is assumed to have same affinity to Gbg as Ga, and no change in Gbg concentration is assumed with concatemer expression. Red bars; M2R-Ga concatemer is assumed to have same affinity to Gbg as Ga, and 1:1 increase in Gbg concentration is assumed with concatemer expression. Green bars; M2R-Ga concatemer is assumed to have 10-fold lower affinity to Gbg than Ga and no change in Gbg concentration is assumed with concatemer expression. For simulation it was assumed that GIRK1/2 is expressed at levels similar to "intermediate density group" described in Yakubovich et al. (2015) i.e. under pre-expression conditions there are~9.7 channels/µm 2 and 3.5 Gbg molecules/channel. It is also assumed that under PTX expression conditions most endogenous Ga i3 is ADP-ribosylated and subsequently degraded. including K D of binding of the first Gbg to GIRK of 60 µM, with a cooperativity factor µ=0.3 for the binding of each next Gbg, to generate a system of differential equations analogous to Eq. 19-34 for the simulation of GIRK1/2 activation. The densities of Gbg and Ga were as for standard conditions, specifically~21 channels/µm 2 and 31 M2R receptors/µm 2 (as for simulation of "high density group" (Yakubovich et al., 2015); 0.5 ng/oocyte RNA of M2R corresponds to~31 receptors/µm 2 ; see Figure 2). Since GIRK1/2 has a substantial GPCR-independent but Gbgdependent I basal [see (Rubinstein et al., 2009) for details], a certain excess of free Gbg over Ga in the vicinity of the channel must be assumed (Yakubovich et al., 2015), whatever the mathematical approach used to describe channel's behavior.
Using the cooperative gating model equations (Eq.14-18) and parameters of GIRK-Gbg interaction, as described by (Touhara and MacKinnon, 2018), we calculate that~7 free Gbg molecules are needed to account for the I basal measured at "high density". We next tested various pairs of Gbg and Ga densities and selected those which most closely recapitulated experimental measures for I basal and I total ( Figure S3). Furthermore, for analysis, we selected the minimal number of Gbg (out of tested pairs), with the corresponding amount of Ga molecules per channel that best reproduced the experimental data ( Figure S3, see legend).
Using the abovementioned initial conditions, we simulated the time-course of GIRK1/2 activation for a range of M2R densities similar to as shown in Figure 4E. Similar to data shown in Figure 3 (for M2R-CFP), Figures 1B and 4E, the amplitude of the simulated I evoked increased with receptor density,whereas t act showed persistent decrease with the growing densities. However, the t act obtained was much faster than that observed in our experiments or predicted by our model, and varied less (i.e., of narrow range; Figure 6B). In order to determine the nature of this discrepancy, we compared similar kinetic steps between our model and that of cooperative gating model. We noticed that the GTP hydrolysis rate used by Touhara et al. is 100 times faster that used in other reports (Zhong et al., 2003). We also note that this fast rate of GTP hydrolysis is characteristic for conditions in which RGS is present, such as observed in cardiac myocytes (Doupnik, 2015). In our case, Xenopus oocytes lack RGS and thus slower GTP hydrolysis rates are expected (Keren-Raifman et al., 2001). Indeed, substituting slower hydrolysis rate constant into cooperative gating model restored channel activation kinetics to as obtained in our experiments ( Figures 6C, D).

Summary of Findings
The main goal of this work was to develop a detailed kinetic model for the GPCR-G-protein-effector cascade based on experimental data obtained in a prototypical expression system, the Xenopus oocyte. As many key features for this cascade are still missing, we deemed it quintessential to measure reactant densities in the plasma membrane along a detailed description of current features (amplitudes and kinetics). We used fluorescently labeled GIRK1/2 as a 'molecular ruler" (Yakubovich et al., 2015) and conducted a controlled and quantitatively monitored expression of the component proteins of the signaling cascade, the GPCR (M2R) (B) -t act of mono-exponential fit of simulated I evoked obtained for a range of M2R densities. For A and B, the rate constants were taken from Touhara et al. (1) and it was assumed that Gbg =16 molecules/channel and Ga = 9 molecules/channel (red circles). (C) -Simulated I evoked obtained from simulation with k hydrolysis = 0.02 s -1 and a range of M2R densities. (D) -t act of mono-exponential fit of time-course from simulation of I evoked with k hydrolysis = 0.02 s -1 and a range of M2R densities. For calculations done in (C, D), k hydrolysis was assumed to be 0.02 s -1 and Gbg = 9 molecules/channel and Ga = 2 molecules/channel (blue circles). and the effector (GIRK1/2). G-protein concentrations were measured in a previous work using two independent methods, the fluorescence measurements with YFP-GIRK1/2 as caliper, and by quantitative Western blots (Yakubovich et al., 2015). Using a titrated expression approach (as previously developed by us; (Vorobiov et al., 2000)), we demonstrate that incremental expression of both wild-type and CFP-tagged M2R, M2R-CFP, is accompanied by acceleration (decrease in t act ) of channel activation, i.e., I evoked . In contrast, activation kinetics of GIRK responses elicited by an M2R-Ga i3 tandem, which necessarily mimics the case of a "preformed GPCR-G protein complex", were fast and invariable over a wide range of expression levels of the fusion protein. We combined these observations and data to develop a comprehensive mathematical model for the M2R-G-protein-GIRK1/2 signaling cascade with free M2R based on free diffusion of all components, and a modification of this model for the case of M2R-Ga i "preformed complex". Our model faithfully recapitulates and predicts all the quantitative aspects of GIRK1/2 activation explored here: the acceleration of activation with increasing densities of M2R, the t act and the amplitude of I evoked , and the dependence of the amplitude, but not the kinetics, of GIRK currents elicited via the activation of the M2R-Ga i3 fusion protein. Our results strongly support the collision coupling mode of signaling between M2R and the G i/o protein in the GPCR-G protein-GIRK cascade reconstituted in Xenopus oocytes. More broadly, our model, and a similar collision coupling type of model of Touhara et al. demonstrate that a purely diffusion-limited coupling mechanism can fully account for the fast kinetics of GIRK responses in excitable cells, without the need to assume a preformed complex. Importantly, our results emphasize the utility of our approach of controlled incremental GPCR expression to distinguish between different coupling modes in G protein-mediated signaling cascades.

Model of G Protein Cycle
The mathematical approach used to describe the GPCR-Gprotein-effector cascade is well elaborated (Lamb and Pugh, 1992;Turcotte et al., 2008;Falkenburger et al., 2010). Several models have been developed specifically for the analysis of GIRK activation. The first model of GIRK activation that also incorporated a G protein cycle of GIRK activation was published in 1988 (Breitwieser and Szabo, 1988). Subsequently, updated models have been developed by groups of Kurachi and Mackinnon Murakami et al., 2013;Touhara and MacKinnon, 2018). Though all these models implement Thomsen-Neubig style G-protein activation model (Thomsen and Neubig, 1989;Zhong et al., 2003), they differ in certain key details. In particular, the model proposed by Murakami et al. (2013) assumes two affinity states of GIRK to Gbg and two affinity states of GPCR to agonist; however, the GPCR-agonist affinity states are not kinetically interconnected and unrelated to coupling to G-protein. These are incompatible with the well-established dependency of agonist affinity upon the G protein-GPCR association (Maguire et al., 1976;De Lean et al., 1980;Haga et al., 1986;Weis and Kobilka, 2018). The model proposed by Touhara and MacKinnon (2018) does not include an explicit bimolecular reaction of agonist binding to GPCR. Both models describe Ga GTP Gbg heterotrimeric complex dissociation as an irreversible reaction, which precludes the implementation of the microscopic reversibility principle (Colquhoun et al., 2004).
The G-protein cycle model presented here ( Figure 4A) is a logical development and, in a way, a synthesis of previously proposed models with certain improvements. In particular, receptor-agonist-G-protein interaction is formulated as a complete ternary complex model (De Lean et al., 1980). Both Ga GTP -Gbg and Ga GDP -Gbg interactions are described as reversible reactions. To enable the implementation of microscopic reversibility in the cycle, we have excluded the obligatorily irreversible GTP hydrolysis ( Figure 4A, reaction 9) from the main cycle. Including only reversible reactions in the circular parts of G-protein activation model makes the model more thermodynamically plausible. Moreover, we do not assume that the proteins in the cascade are in unlimited supply, and the equations of the model completely take into account "depletion" of proteins by excess of their binding counterparts. The only exception is GTP and GDP that are assumed to be in unlimited supply ("free ligand approximation", e.g. as is done for agonist in standard descriptions of agonist-receptor interactions).
Combining G-protein cycle model with GIRK1/2 gating model successfully reproduced the experimental observations. In particular, simulation of M2Rwt expression experiment demonstrated similar decremental t act dependence on receptor density (Figure 4), a feature that was subsequently nearly abolished when modeling the M2R-Ga tandem activation process ( Figure 5). These findings strengthen the notion that, in Xenopus leaves oocytes, M2R and G-proteins are not in "preformed complex", rather interact reversibly.

Channel Activation Models
The process that leads to GIRK opening following the binding of Gbg is still unclear (Glaaser and Slesinger, 2017). Currently there are three detailed GIRK gating models developed chronologically by Kurachi group (Hosoya et al., 1996;Murakami et al., 2010;Murakami et al., 2013), our group Yakubovich et al., 2015) and MacKinnon's group Wang et al., 2016;Touhara and MacKinnon, 2018). The model proposed by Kurachi's group is based on Monod-Wyman-Changeux allosteric model of GIRK activation and formulates two binding states of Gbg for each channel subunit. To complete the model, the authors used sub-nanomolar affinity for the GIRK-Gbg interaction. The model developed by MacKinnon's group includes a detailed binding reaction of Gbg to GIRK2 and GIRK1/4 and is based on elegant experiments with purified GIRK2 and Gbg incorporated into bilayer membranes. In their work, they suggest that Gbg binding is cooperative and that only channels occupied by four Gbg undergo the closed-open transition . The model presented in the current work is based on sequential binding of Gbg molecules to GIRK1/2 with graded contribution of each Gbg-occupied state to open probability. Notably, this is based on ample experiments using the homologous heterotetrameric GIRK1/4 (Ivanova-Nikolova and Breitwieser, 1997;Ivanova-Nikolova et al., 1998;Sadja et al., 2002). Considering the different affinity values, reported in the literature, for GIRK-Gbg interaction (see below), we derived the K D values of this interaction from analyzing the available crystal structure of GIRK2-Gbg (Whorton and MacKinnon, 2013) and the docking model of GIRK1-Gbg (Mahajan et al., 2013). Notably, both estimates rendered K D values (100-200 nM) on par with those measured in biochemical experiments, 100-800 nM (Krapivinsky et al., 1995;Doupnik et al., 1996). Furthermore, our model of channel activation, combined with the G protein cycle model, reliably reproduced the amplitudes of GIRK1/2 activations when we used the surface densities of the channel, GPCR and G-protein directly measured under physiological conditions. We note that, with an appropriate adjustment for the oocyte expression system, the cooperative gating model of Touhara and MacKinnon (2018) also reproduced the experimental data reported here, in terms of GIRK1/2 currents, t act and the dependence of t act on M2R density ( Figure 6). Interestingly, simulations done with the original parameters of the G-protein cycle model of Touhara et al. initially yielded a relatively shallow dependence of t act on M2R density, despite the fact that it is a collision coupling-type model. Analysis of this discrepancy lead us to the finding that t act dependence on receptor density even in "collision-coupling model" can be masked by a rapid rate of GTP hydrolysis, as can be obtained when RGS proteins are present ( Figure 6).
A more consistent discrepancy of this model's estimate with our estimates is related to estimations of the number of Gbg and Ga molecules required for channel activation. The endogenous levels of Gbg in oocyte's plasma membrane is in the range of 20-40 molecules/µm 2 , and is further increased to~80 molecules/µm 2 (Yakubovich et al., 2015) upon overexpression of GIRK1/2 that recruits additional Gbg to the plasma membrane Kahanovitch et al., 2014). Calculations with the Touhara et al. model showed that a minimum of 7 free Gbg molecules need to be available at any time to account for the basal, GPCR-independent GIRK1/2 activity; at least 9 Gbg and 2 Ga are needed to account for the observed I total (I basal +I evoked ) ( Figure S3). For 21 GIRK1/2 channel/µm 2 (as observed after injection of 1 ng RNA of each subunit; (Yakubovich et al., 2015)),~75 Gbg molecules per m 2 are required to attain the total GIRK1/2 current according to our model, but~180 Gbg molecules per m 2 are needed with the cooperative gating model. The differences in Gbg estimates stem mainly from the distinct GIRK-Gbg affinity estimates used: 60 µM for the first Gbg bound, with progressively improved affinity for each following Gbg in Touhara and MacKinnon (2018), vs.~0.15 µM in our model ( Table 1). An even greater discrepancy may be expected if one uses the K D of~300 µM for GIRK2-first Gbg interaction, as estimated experimentally in lipid bilayer in the presence of Na + (the natural condition in a living cell's cytosol) . Touhara and MacKinnon (2018) also note that, with such low affinity of interaction between GIRK and Gbg, Gbg surface densities needed for GIRK activation should be much higher than the physiological range. This leads them to propose that the GPCRs, G-proteins and GIRK channels interact in hot spots (Sungkaworn et al., 2017), where all components of the cascade are highly concentrated. Another possible explanation is that the actual affinity of GIRK-Gbg interaction is higher than 200-300 µM K D despite the measurements in lipid bilayer experiments  or in solution by nuclear magnetic resonance (Yokogawa et al., 2011).
Notably, these studies employed a mutant Gg lacking lipid modification (geranylgeranylation, in the case of Gg 2 ) and it is established that lipid modification of Gg is an important determinant of high-affinity interaction between Gbg and its binding partners such as Ga and phosducin (Myung et al., 1999;Lukov et al., 2004). Further study is needed to better determine the affinity of GIRK-Gbg interactions in living cells.

Collision-Coupling Versus Preformed Complex
There is a longstanding debate regarding the existence, or lack, of diffusion-dependent steps in the GPCR-G-protein cycle. A large body of data suggests pre-formed (or dynamic) complexes between these proteins. For instance, measurement of diffusion coefficients of GPCRs demonstrate non-homogeneity, pointing to partial restriction of diffusion and possible organization of GPCRs in "islands" (Daumas et al., 2003;Suzuki et al., 2005;Baker et al., 2007). Similar restriction in lateral mobility is noted for G protein (Kwon et al., 1994). These observations are supported by studies showing that immobilized GPCRs can activate G-protein molecules with the same rate as mobile GPCRs (Lober et al., 2006). Going downstream in the cascade, the dissociation of Ga from Gbg during the activation process also remains questioned. Whereas the dogma states full dissociation and diffusion of the latter, it has been shown that some Ga and Gbg may not undergo dissociation, rather undergo spatial rearrangement after activation (Bunemann et al., 2003;Lambert, 2008). Onwards, complexes between G-protein and effectors, such as GIRK, have been noted Zylbergold et al., 2010) as well as even larger supramolecular complexes consisting of GPCRs, G-proteins and modulating molecules, such as RGS, have been demonstrated (Lavine et al., 2002;Doupnik, 2008).Despite this body of work, pure collision coupling has also been demonstrated in many other cases (see introduction). Notably, the distinction between the two different modes is not trivial. In the current study, we elaborate our protocol using "titration" of proteins densities at the membrane, specifically those of the GPCR (Vorobiov et al., 2000) and show that it allows to quantitatively distinguish between the modes of GPCR-G protein coupling.
We compared the kinetic properties of two different settings of M2R-G-protein interaction. The first consists of M2R-Ga fusion protein; a one-to-one relationship between the GPCR and Ga is enforced, giving rise to a bona fide "preformed complex". Notably, in the scenario, we also assume that the complex necessarily includes Ga GDP -Gbg (in view of their very high nanomolar affinity), resulting in a full GPCR-G protein "preformed" complex. The second setting involves independent, untethered proteins. These two scenarios reveal that, whereas the kinetics of activation of GIRK by M2Rwt are highly dependent on the receptor density, those of M2R-Ga fusion are not. This clear distinction strengthens the idea that, at least in Xenopus oocytes, M2Rwt can indeed diffuse in the plasma membrane to activate several G-protein molecules (and GIRK subsequently). This conclusion is supported by studies, specifically conducted in Xenopus oocytes, demonstrating that, M2Rs and G-proteins are not permanently co-localized and diffuse unrestrictedly in the plasma membrane (Hein and Bunemann, 2009).
In summary, the current concept of GPCR-G-protein effector signaling may be schematically presented as three possible, and not mutually exclusive, paradigms. The first, the preformed complex model, is expected to follow first order kinetics, in which the rate of activation is concentration-independent. In the second, the catalytic collision coupling model, the rate of activation is anticipated to be highly concentration-dependent, at least for diffusion-limited cases. Of note, the dependence of activation kinetics on receptor density might also be influenced by G-protein inactivation rate. Indeed, we demonstrate this by employing the simulation of Touhara et al. with the slower GDP hydrolysis rate ( Figure 6). The third, is a mixture of the two. Importantly, the relationship of activation rate and concentration of the reactants is not trivial and is formulated, in most cases, as differential rate laws rather than as integrated rate law utilized for first order kinetic processes (useful for preformed complexes). It was previously shown that coupling reaction rate constant under diffusion limited conditions is equal to diffusion transport constant (Lauffenburger and Linderman, 1996) and is therefore expected to be dependent on receptor density (Mahama and Linderman, 1994;. Together, the dependence of GIRK's activation rate on the density of M2R most likely originates from this phenomenon as a result of collision coupling of GPCR and G-protein. The recently proposed "hot spot" interaction model of GPCR-G-protein activation (Suzuki et al., 2005;Sungkaworn et al., 2017) represents a particular case of collision coupling model (but not preformed complexes) and elegantly describes cases in which there is relatively low affinity between the reactants, because restrictions of molecules within a tight hot spot is expected to robustly increase their effective concentration. Further studies based on stochastic analysis of GPCR-G-protein-GIRK system and measurement of GIRK-Gbg affinity will deepen our understanding of the above described phenomena.

DATA AVAILABILITY STATEMENT
The raw data supporting the conclusions of this article will be made available by the authors, without undue reservation.

ETHICS STATEMENT
The animal study was reviewed and approved by Tel Aviv University Institutional Animal Care and Use Committee (permits M-08-081 and M-13-002).

AUTHOR CONTRIBUTIONS
SB, EA, RH-J, UK-electrophysiological and imaging experiments. DY, HP-generation of the model. ND, SB, DYplanning of the study. DY, ND, SB-writing the first draft of the manuscript. All authors contributed to the article and approved the submitted version. SUPPLEMENTAL FIGURE S2 | Linear relationship between Gbg-activated YFP-GIRK1/GIRK2 current and total surface density is observed in the range 0-1 ng of YFP-GIRK1 RNA. (A)representative confocal images of oocytes expressing YFP-GIRK1/GIRK2. Channels were expressed by injecting the indicated doses of RNA of YFP-GIRK1 and GIRK2 (1:1) and activated by coexpression of 5 ng Gb and 1 ng Gg RNA. (B)correlation between total GIRK surface density and the Gbg-dependent GIRK current, I bg . I bg is the total agonist-independent current in Gbg-expressing oocytes. It was measured in 24 mM K + solution by subtracting the non-GIRK currents remaining after inhibition of >95% GIRK activity by 5 mM Ba 2+ (Rubinstein et al., 2007). Total surface density is reflected in YFP fluorescence levels, functional channel density is proportional to I bg . Correlation between I bg and fluorescence is linear for RNA doses of YFP-GIRK1/GIRK2 up to 1 ng/oocyte of each subunit. n = 6-12 for each point.  Figures 6C, D. All calculations except the 9:2 Gbg/Ga pair (blue bar) have been made for the condition of fast hydrolysis of GTP, k hydrolysis =2 s -1 . The 16:9 Gbg/Ga pair was the lowest dose of G protein subunits per channel that reproduced I total under this condition. Lower doses of Gbg/Ga produced lower I total . The 9:2 Gbg/Ga pair was the lowest G protein subunits combination that reproduced I total under the condition of low GTP hydrolysis rate (k hydrolysis =0.02 s -1 ).