Model of ligand-triggered information transmission in G-protein coupled receptor complexes

We present a model for the effects of ligands on information transmission in G-Protein Coupled Receptor (GPCR) complexes. The model is built ab initio entirely on principles of statistical mechanics and tenets of information transmission theory and was validated in part using agonist-induced effector activity and signaling bias for the angiotensin- and adrenergic-mediated signaling pathways, with in vitro observations of phosphorylation sites on the C tail of the GPCR complex, and single-cell information-transmission experiments. The model extends traditional kinetic models that form the basis for many existing models of GPCR signaling. It is based on maximizing the rates of entropy production and information transmission through the GPCR complex. The model predicts that (1) phosphatase-catalyzed reactions, as opposed to kinase-catalyzed reactions, on the C-tail and internal loops of the GPCR are responsible for controlling the signaling activity, (2) signaling favors the statistical balance of the number of switches in the ON state and the number in the OFF state, and (3) biased-signaling response depends discontinuously on ligand concentration.

Here, [K] and [P ] are kinase and phosphatase concentrations respectively. The forward and back reaction rates are k + = α + [K] and k − = β + [P ], respectively. In equilibrium, the forward and microscopically reversed rates are in detailed balance as are the back and microscopically reversed rates. This is not true in biologically relevant situations. In that case the forward and back rates are much larger than the respective microscopically reversed rates. The further into nonequilibrium the process is, the more efficient is the information transmission. To see this, note that the capacity , the maximum information flow, in a switch is given approximately by Qian and Roy (2012) The capacity C i in the switch is maximized when the reaction rates α + and β − in the clockwise direction is much greater than the microscopically reversed rates α − and β + in the counter-clockwise direction. This is just the condition for the process to be far from equilibrium. In equilibrium, the forward rates are equal to the reverse in detailed balance.
The further the rates are from detailed balance, the further the process is from equilibrium. The process is driven far from equilibrium by free energy forced into the system that is in excess of the thermal energy from the heat bath.
The free-energy surplus is a result of the concentration of ATP being held out of equilibrium with its reactants ADP and phosphate Qian and Reluga (2005). Under conditions in the cell, this energy is about 12 kcal/mol (Berg Jeremy et al., 2019, Sec. 15.2) much higher than the thermal energy of 0.6 kcal/mol ≈ 1/40 eV and much lower than the energy required to break covalent bonds, about 100 kcal/mol ≈ 4 eV . For these energies, the reaction rates obey the condition Qian (2007) α + β − α − β + ≈ 10 10 (S6) which, from Eq. S3, indicates that the capacity is quite close to its maximum loss-free value for biological processes. For biological systems driven by ATP, switches are nearly loss-free.

General BOIS Model
We assume one ACTIVE GTPC and N ACTIVE PdPCs such that M = N + 1. The steady-state entropy production rate σ for all the switches is Qian (2003) where T is the temperature of the heat bath in energy units, J i is the chemical flux of a single receptor passing through the receptor states of PdPC i, and µ i ≤ 0 is the free-energy cost for a receptor to pass through the PdPC, with similar expressions for the GTPC. Here, Pr( ON | i ) is the probability that ACTIVE switch i is ON, and Pr( OF F | i ) is the probability that the ACTIVE switch i is OFF. The probabilities for the INACTIVE states are all zero.
The entropy production, Equations S7, along with the normalization and free-energy constraints can be maximized with the method of Lagrange multipliers. Note that, in order to keep the value of the probabilities between zero and one, k ± i and k ± G cannot be smaller than J i and J G , respectively.
The total free-energy drop ∆G is finite. We quantify this by setting the average drop µ 0 over switches to be constant.
where µ i and µ G are the free-energy drops in the transit of PdPC and GTPC, respectively.
The entropy production rate with constraints in terms of Lagrange multipliers λ 1 and λ 2 is The maximum production is obtained by taking the partial derivatives of Eq. S11 with respect to µ i and J i and then setting the derivative to zero.

Specialized BOIS Model
The receptors in the OFF state are pooled in the specialized model. In that case, mass conservation in the general model, Eqs. S8 and S9, is replaced with where C is the receptor concentration of the C state, the part of the GPCR complex that excludes the ON states of the switches and the receptors unbound from ligands. The only difference between the specialized model and the general model is that Eq. S29 becomes with a comparable expression for G. The specialization only affects the value of the forward reaction rate k + i for ON and OFF states. The switch is OFF if

Switch Manipulation
Switch configurations can also change by interchange of switches. The first switch and the third switch in a configuration may swap locations on the C tail. If the one switch is ON and the other is OFF, then the configuration is altered by the exchange. Finally, the configuration can change by transmitting the information from one switch medium to another. For instance, a configuration may exist among the switches on the C tail of the GPCR complex. If the information is transmitted to the medium of effector coupling pathways, then the configuration will be altered. The same information will be transmitted, but the information will be described in a different way. The entire concept of switch location may be different in the two media.
The question remains as to what are the properties of the catalyst or catalysts that turn a switch ON and OFF. We have an example of a catalyst that does this. It has been observed that arrestin blockage of the G α pathway in a GPCR complex is a situation in which a PdPC catalyzes a GTPC switch to turn OFF ( (Violin et al., 2008, Eq. 7)). In the BOIS picture, this suggests that a phosphorylated site on the arrestin-bound C tail of the GPCR comes into contact with the GTPC, which is illustrated in Fig. S1B. The BOIS model predicts that switches must come in contact with each other in order to catalyze the switches to transition from ON to OFF and back. It is plausible that this can be achieved if the the C-tail/arrestin complex is flexible. It has been shown that arrestin-bound C tail is very flexible and can take on several conformations on and around the GPCR complex (Kahsai et al. (2018); Latorraca et al. (2018); Eichel et al. (2018)). In the GTPC case displayed in Fig. S1C, the ON state of a PdPC catalyzes the transition of the GTPC from ON to OFF by increasing the k − G . The GTPC OFF state inhibits k − 2 in a second PdPC, which turns that switch ON.
This suggests at least one mechanism for catalyzing the switches in general, illustrated in Fig. S1A. Here, a switch 1 in the ON state comes into contact with switch 3 in the ON state and catalyzes the transition of switch 3 from ON to OFF. The reaction takes k − 3 /J 0 from one to zero. The OFF state in switch 3 may also inhibit the reaction represented by k − 2 to turn switch 2 OFF. In other words, an ON state of one switch can catalyze another switch to turn OFF, while the OFF state can inhibit another switch to turn ON.
Figure S1. Speculative model: ON states on a switch catalyze another switch to go from OFF to ON, while OFF states of switch inhibit another switch's reaction rate k − i causing that switch to turn ON. A. Switch 1 comes within proximity to Switch 3. This is due to flexibility of C tail. The ON state of Switch 1 catalyzes Switch 3 to turn OFF (k − 3 becomes large). The OFF state of Switch 3 inhibits k − 2 and turns Switch 2 ON. B. Structure of the G α GT P complex (purple) and the arrestin-bound C tail of the GPCR complex during inactivation of the G α pathway. An initial configuration of phosphorylation switches represented by ON Recruiting Switch, recruits arrestin to the C tail of the GPCR complex. An ON phosphorylation site R * 1 on the arrestin/C-tail complex comes in contact with the G α subunit turning the GTPC OFF. The OFF state of the GTPC inhibits the reaction from ON to OFF of an OFF phosphorylation switch R 2 to turn it ON, thus maintaining the number of ON switches and OFF switches. C. Chemical schematic in the form of A. to describe the inactivation of the G α pathway.

Frontiers
These arguments have implications for transitions between and within the configuration groups illustrated in Fig. ?? for M = 4. If one switch catalyzes another switch to turn ON or OFF, then a transition is made between two groups in Fig. ??. If two switch states are changed as in Fig. S1A, then the switch configuration transitions between two configurations within a group. Changing the on-state of a single switch transitions between groups, while pairwise transition between switches within a group transition between two states within the group. Transitions within a group are particularly interesting. These transitions do not involve additional investment or loss of total information in the group as would occur if a sole switch were turned ON or OFF.
We propose here three mechanisms, C, P, and T for changing switch configurations while maintaining the number of ON switches. As mentioned, the first mechanism C changes of the switches by catalytically manipulating the return reaction rates k − i and k − G (Eqs. ?? and ??). We can represent this process mathematically as manipulation of a pair of ON/OFF states (0,1) and (1,0) and two pairs of states in which switches are both ON or both OFF (1,1) and (0,0).
Since the catalytic process of Fig. Fig. S1 only operates on two switches with opposite ON/OFF states, the C operator does nothing to the states (1,1) and (0,0), which implies where I is the identity operator. Two applications of C return you to your original state. In physics, this is known as charge symmetry (Lehnert (2016)).
The second mechanism P swaps the two switches from different sites. The transformed pair appears to be a mirror image of their pre-transformed configuration. This is a parity transformation. The parity operator P has the same set of transformations, Eqs. S31 and S35, as the catalytic operator C. This implies, not only is the operator PP invariant the two operators P and C are PC invariant (Lehnert (2016)).
The mechanisms C and P can both occur in the same population of switches, such as the switches on the C tail of the GPCR complex. The third mechanism T transitions between two distinct switch populations, such as the phosphorylation sites on the C tail and the conformational switches on the arrestin molecule (Latorraca et al. (2020)). In other words, the first transition in a cycle in Fig. ?? may occur on the C tail, while the next transition may occur as the information passes from the C tail to the arrestin molecule. The T transition can happen between any of the state pairs within a group (Fig. ??). A set of transitions is the identity operator if a cycle of transitions returns to its starting state.
where n is the total number of transitions T in the closed cycle. For n = 2, as would be the case if the chemical flow was from the C tail to the arrestin molecule and back, the cycle can be thought of as forward transition in time followed by a backward transition in time. This is known as time-reversal symmetry in physics (Lehnert (2016)). If we focus on the pair of switches, one with a location on the C tail and one with a location on the arrestin molecule that map into each other, then the mapping is equivalent to the C and P mappings: and T (1, 0) = (0, 1) if the switches change their ON/OFF state and the identity if they do not.
Therefore, any closed cycle of configurations within a group is invariant under some combination of the three operators C, P, and T This is known as CPT symmetry. What this means in non-mathematical language is that any state within a group can be accessed by a series of catalytic transitions between ON and OFF switch states, with an interchange of switch positions, and by transitioning from one switch medium to another (eg. C tail to arrestin). This transition can be of endogenous or pharmaceutical origins. Practically, this means that the three processes, turning switch pairs ON and OFF catalytically, swapping the locations of any two switches, and changing the locations of two switches on two separate switch media (C tail and arrestin) can change the code for effector coupling. Any state with the same number of ON states can be accessed by a combination of these processes. This increases the number of potential drug targets.

Simulation and Experiment
The MATLAB code that plots the figures is BOIS.mlx. Output inages are stored in the Graphs folder. The data files that the code reads are in the same folder. The Simbiology add-on must be installed to generate data files. The code for generating the BOIS data is in Simbiology project Central.sbproj. The code for generating Mass-Action data is in MassAction.sbproj.

Dose Response
We examine the relationship between forward and backward processes in the constrained BOIS model of Fig. S2B, which is used in the simulation. The flux relationship for ON switches, (R * i = 1), is Figure S2.
The total normalized receptor concentration R T is where R F is the normalized concentration of free receptors, those not bound to a ligand.
where θ is the Heaviside step function, L is the ligand concentration and K D is the ligand-receptor dissociation constant. The step function is required to keep the concentration R F equal to or greater than zero. The ligand-bound concentration C is A set of numerical experiments were performed to determine the threshold for activation of the switches. The architecture was composed of four switches, one GTPC and three PdPCs; of which one, two, three, or four were ON, while the remaining were OFF. All receptors were initially in state C. No receptors were initially in the ligand-free state or in the switch ON states. The total number of ON switches was M o = (1, 2, 3, 4). We found that switches did not activate (J 0 > 0) unless C > M o at the initial time t = 0. Once an ON (k − i = J 0 ) switch was activated, it remained activated irrespective of the concentration C. For cases in which the normalized concentration C did not initially exceed the number of ON switches M o , the simulator failed to evolve to a computable steady-state solution. This indicates that physically meaningful numerical solutions are not possible if the normalized receptor concentration in state C is less than M o during an initial increase in ligand concentration. In the initial sub-threshold situation, the concentration C may be nonphysically driven below zero, or the code is simply unable to find an integration time step for any finite value of J 0 . Once the ligand concentration exceeds the value in which C is greater than M o , then 2M o switches are activated with half the switches in the ON state and half in the OFF state. If the concentration C subsequently drops below M o , the ACTIVATED switches remain activated. No nonphysical values of concentration are generated and the integration runs smoothly.
In a second set of numerical experiments that mimicked assay experiments, receptors were initialized to lie in the ligand-free state R F instead of the ligand-bound OFF state C. The ligand concentration starts at zero and slowly increases to a value much larger than the dissociation constant K D . As the ligand concentration increases, receptors move from state R F to the GPCR-complex state C. At low values of C concentration, the number of receptors in the ON state of the switches is zero, no flux is permitted into the switch. When the complex receptor concentration C exceeds the ON switch number M o , the flux J 0 turns on for both ON and OFF switches. Once this event has occurred, the switches remains active even if the concentration C drops below the activation threshold. The activation threshold is given by We define EC 50 as the ligand concentration L at which the equality holds.
Runs with the ligand concentration normalized to EC 50 and for values of ON switch number M o = (1, 2) and for total normalized receptor concentration R T = (0.1, 1, 10, 100) are displayed in Fig. S3A. The figure illustrates dose-response curves for both PdPCs and GTPCs. Note that the curves are nearly universal when displayed as a function of the ratio of ligand concentration L to EC 50 . This implies that the dose-response curves are completely characterized by the quantity EC 50 defined in Eq. S56.

Differing Responses for Rigid and Flexible Proteins
Equation S56 can be interpreted in two ways depending on whether the receptor proteins are rigid or flexible. For rigid proteins, the number of ON switches is a constant and determined by the rigid bonds among protein elements. In that case, the ON number M o and the total receptor concentration R T determine the single common ligand concentration EC 50 to which all the dose-response curves are scaled. Equation  S56 indicates that for rigid proteins in which the number of ON switches is constant for the ensemble, the ON switches become ACTIVATED at a common value of L = EC 50 . In other words, all switches have the same value for EC 50 . This is not seen in experimental assays Rajagopal et al. (2011). Rather, different downstream pathways are associated with different values of EC 50 . This is further evidence that switches cannot be completely composed of rigid proteins.
The number of ON switches M o available for activation can be mutable in flexible proteins, however. The protein environment, including the ligand binding, may select the number of ON switches available to be activated at any given value of ligand concentration L. In that case, Eq. S56 can be interpreted as determining the number ON switches M o in the complex given the total receptor concentration R T and EC 50 . Moreover, it is well known Lefkowitz et al. (1993) that the activation of G proteins is a tertiary reaction among the ligand, receptor, and G protein. This implies that the downstream pathways affect the ligand dissociation constant K D . We allow for the possibility that the PdPC activation can also affect the ligand dissociation K D .
In the case of flexible protein conformations, the picture that emerges is: 1.At the initial time t = 0, the ligand concentration L is zero as are the GPCR -complex concentration C.
The number of ON switches is zero as is the number of OFF switches. 2.The ligand concentration L is increased until it reaches a point in which Eq. S56 for M o = 1 is satisfied.  1, 1, 10, 100). The dose response curves are nearly universal for this scaling. The curves can be characterized by the reference concentration R ref and the EC 50 . B. Curves of constant total receptor concentration. The total receptor concentration varies from R T = 1.5 to R T = 10.5 in steps of 1 in normalized concentration units. The curve for R T = 1.5 is the single point on the left. The curve for R T = 10.5 is the curve containing 10 points on the right. We see that increasing the number of active switches M = 2M o while holding the receptor concentration constant increases EC 50 for ligand-receptor complex.

3.At that point an ON (k
switch is ACTIVATED, thus satisfying the requirement that the number of ACTIVE ON switches is equal to the number of ACTIVE OFF switches.
4.The concentration of state C is drained into the new ACTIVE ON switch and the process repeats.

5.Switches are ACTIVATED as the ligand concentration is increased until the addition of a new ON switch
would require more receptors than are available at the given total receptor concentration R T .
We can make this argument quantitative. We define M max o as the maximum possible number of ON switches M o , both ACTIVE and INACTIVE.
where the excess receptor concentration ∆R T is Initial values are L = 0, M o = 0, R F = R T , and C = 0. The ligand concentration L increases until it reaches a value in which the threshold requirement, Eq. S43, is reached for M o = 1.
where the arguments of EC 50 and K D indicate that this is the first threshold that is met. When this ligand-concentration threshold is met, the receptors drain into the newly ACTIVATED ON state and the GPCR-complex concentration C becomes The first estimate for The ligand concentration L continues to increase until the second threshold is reached and another pair of switches become activated.
The second measurement of R T is The next estimate is In general, the n th estimate is The process continues until the final estimate of R T is or The flexible version of Eq. S43 is, for the n th vale of EC 50 /K D , . This is displayed as the black curves in Fig.S3B. The comparable expression for Rigid proteins is, from Eq. S56, This is displayed as red curves in Fig. S3B.

Observed Dissociation
Each pathway has its distinct ligand dissociation value EC 50 . We can also calculate a global dissociation K GD constant defined as the ligand value at which half of the receptors are free of ligands and half are bound with ligands. This is quite different from the EC 50 value calculated for a given downstream response to a given ligand. In that case, the dissociation EC 50 is the ligand concentration at which a switch becomes activated. Intuitively, it can be thought of as the ligand concentration at which a switch is halfway turned ON, at which the receptor concentration in the ON state is half the reference concentration.
This dissociation is also quite different from the dissociation constant K D determined by the reaction rates in ligand/receptor binding. Here, K D is associated with the reaction where the concentration of C, denoted by C, is related to the free-receptor concentration R F and the total receptor concentration R T by This is the case in which no switches are ACTIVATED. In that case, which is quite different from the the dissociation constants when switches are ACTIVATED. It seems to be very important to identify in experiments, which type of dissociation is being observed.
To calculate K GD , note that the concentration of ligand-free receptors is The observed value of the dissociation constant K obs D occurs when the ligand-free concentration is half the total concentration where M o (1/2)is the value of ON switch number that corresponds to half the total number of switches being ON. If M o is even, then can be written

Experimental Observations
Observations of EC 50 from Ref. Rajagopal et al. (2011) are given in Table S1.
Assay results on two different receptors, Adrenergic and Angiotensin II, are displayed. Two different downstream pathway concentrations were measured for each receptor, the G α pathway and one β-arr pathway. Other pathways may have been present, but they were not measured. The estimated reference receptor concentration R ref was taken to be the maximum Formoterol concentration for the Adrenergic receptor and Angiotensin II maximum concentration for the Angiotensin II receptor. It should be noted that the experimental observation used a common reference concentration for receptor switch states before its theoretical prediction by the BOIS model. The maximum normalized downstream concentrations R max for each pathway are given in columns marked G α and β arr. The EC 50 concentrations are measured in moles per liter. Several distinct ligands, listed in the table, were applied to each receptor. Table S1. Assay Data. Assay values for maximum stimulus for Adrenergic and Angiotensin receptors Rajagopal et al. (2011). The concentrations Gα and β arr are normalized to the reference value R ref . Ligand concentrations are molar. The ligand names correspond to the names used in Ref. Rajagopal et al. (2011). Balanced ligands for the Adrenergic receptor are Formoterol, Isoproterenol, Epinephrine, and Fenoterol. Balance ligands for the Angiotensin II receptor are Angiotensin II, TRV0120055, TRV0120056, A1, and S1C4. Adrenergic ligands that are Gα biased are Dobutamine, Norepinephrine, Clenbuterol, Salmeterol, and Salbutamol. No adrenergic ligands in the sample are β-arr biased. Angiotensin II ligands that are β-arr biased are TRV0120044, TRV0120045, TRV0120034, and S1G4G8. No angiotensin ligands in the sample are Gαbiased.