Abstract
A mathematical model for the binding and function of metabotropic glutamate receptors was developed, with the aim to gain new insights into the functioning of these complex receptors. These receptors are homodimers, and each subunit is composed of a ligand binding [Venus flytrap (VFT)] domain and a heptahelical domain (HD) responsible for Gprotein activation. Our mechanistic model integrates all structural information available so far: the various states of the VFT dimer (openopen, closedopen, and closedclosed), as well as the fact that a single HD is active at a time. To provide the model with parameters with biological meaning, two published experimental studies were reanalyzed. The first one reports a negative cooperativity in agonist binding (J Biol Chem 279:35526–35534, 2004), whereas the other indicates a positive cooperativity in agonistmediated response (Nat Struct Mol Biol 11:706–713, 2004). The former study allowed us to explain the mechanistic features associated with VFT recognition by agonists and antagonists integrating a negative allosteric interaction for agonist binding. The second study helped us to quantitatively describe the functional dynamics of transduction of the VFT occupation into functional response, confirming a putative positive cooperativity at the level of receptor coupling efficacy. This model will help both to better understand the functioning of these receptors and to characterize the mechanism of action of various types of allosteric modulators. Moreover, this model may be of general utility for oligomeric systems in which the ligand binding and effector domains correspond to distinct structural domains.
The two main neurotransmitters, glutamate and GABA, activate not only ionotropic receptors but also Gproteincoupled receptors (GPCRs), called the metabotropic glutamate (mGlu) and GABA_{B} receptors, respectively. These receptors belong to the class C of the large GPCR family (Pin et al., 2003) and play essential roles in the central nervous system by regulating fast excitatory and inhibitory transmission. As such, these receptors are the subject of intense research for the development of new drugs targeting the central nervous system.
In addition to their sequence divergence with the other GPCRs, mGlu and GABA_{B} receptors have peculiar structural characteristics. These receptors form constitutive dimers that are stabilized by a disulfide bridge in the case of mGlu receptors. Moreover, each protomer of a mGlu dimer is composed of three main structural domains: an extracellular Venus flytrap (VFT) domain where agonists bind, a transmembrane heptahelical domain (HD) responsible for Gprotein activation, and a cysteinerich domain (CRD) that interconnects the VFT and the HD both structurally and functionally (Rondard et al., 2006). These structural features make these receptors complex proteins and raise several issues regarding how agonist binding in the VFT leads to Gprotein activation by the HD.
Important information on the functioning of these proteins has been obtained from mutagenesis and structural studies (Bessis et al., 2000, 2002; Galvez et al., 2000; Kunishima et al., 2000; Tsuchiya et al., 2002). It is now recognized that agonist binding in the VFT stabilizes the closed state of the VFT, which stabilizes in turn a new position of the VFTs relative to one another in the dimer. Such a relative movement of the VFTs has been proposed to favor a new relative position of the HDs leading to the activation of one of them (Hlavackova et al., 2005).
Mathematical models, empirical or mechanistic, have been widely used to study complex systems in biology. The former models pose the main advantage of being simple, offering, therefore, the easy determination of crucial parameters, which, although lack physical meaning, serve to characterize the general functioning properties of these systems. In contrast, mechanistic models aimed at mimicking mathematically the expected functioning of the protein complex can help to validate details of the proposed mechanism and can provide important information to better understand the mechanism of action of modulatory compounds. Indeed, not only the mGlu receptor structure is complex, but also its functional properties. For example, whereas positive cooperativity functional effect of the agonist was reported (Kniazeff et al., 2004) by analyzing the response of these dimeric receptors, a negative cooperativity between agonists binding sites was found (Suzuki et al., 2004) using binding experiments. Moreover, in addition to agonists and antagonists, a number of allosteric modulators with either positive (PAM) or negative (negative allosteric modulator) effects have been identified, the PAMs enhancing either agonist affinity or agonist potency or both (Goudet et al., 2004).
In this study, we present a mechanistic model for the ligand binding and the functioning of the dimeric mGlu receptors. This model accommodates very well both binding and functional data published on these dimeric receptors and rationalizes the different kind of cooperativity reported for agonist binding and agonistmediated functional effects.
Materials and Methods
Modeling ConcentrationSignal Curves. Empirical and mechanistic models are commonly fitted to curve data points to provide parameter estimates for comparison between either ligands or receptors. In the empirical case, the parameters lack physical meaning, whereas in the mechanistic, the parameters are the chemical constants of the process. Concentrationsignal curves (where signal stands for bound ligand or functional response) are normally depicted using the logarithm function for the ligand concentration, leading to sigmoid curves. In this article, a mechanistic model was developed, which, when convenient, was handled empirically by grouping chemical constants into empirical parameters.
Curve Shape Analysis. Quantitative curve shape analysis is necessary for accurate comparisons between concentrationsignal curves. The following pharmacologic descriptors can be used for the shape analysis of f(x) curves, where f stands for the signal and x = log[A], being [A] ligand concentration (Giraldo et al., 2002; Giraldo, 2003):

The right asymptote of f as x increases: .

The left asymptote of f as x decreases: .

The location of the curve or midpoint (x_{50}): x for f = Left + 1/2 (Right – Left). Note that we use left and right instead of the more common terms bottom and top because for inverse agonists in functional studies, right is lower than left, and bottom and top would not be appropriate.

The slope of the curve at the midpoint (x_{50}). That is, the value of the first derivative of the function at x_{50}: (df/dx)_{x}_{50}. This descriptor is directly related to the Hill coefficient at x_{50}: n_{H50} = 4/[(Right – Left) · ln 10] × (df/dx)_{x}_{50}.

The index of asymmetry of the curve is the distance between x_{50} and the point of inflection (x_{I}). The point of inflection is a point on a curve at which the curvature changes from convex to concave or vice versa. At x_{I}, the first derivative of the f(x) function is a maximum or minimum, and, as a consequence, the second derivative is equal to zero. The location of the point of inflection serves for the assessment of the symmetry of the curve. An f(x) curve is symmetric if the point of inflection matches the midpoint, x_{I} = x_{50}, and asymmetric if it does not, x_{I} ≠ x_{50}.

The number of inflection points to distinguish between monophasic curves (one inflection point) and biphasic curves (three inflection points).

The grade of positive or negative cooperativity of one bound ligand onto the binding of the next for a dimeric receptor: the difference between the slope at the midpoint for the concentrationsignal function (f) under study and the monophasic f = 1/(1 + 10^{x50–x}), where the latter presents a slope parameter of one.
The Problem of Parameter Estimation in Overparameterized Models: Evolutionary Algorithms versus Classic Gradient Fittings. Classic fittings by gradient nonlinear procedures pose the drawback of their dependence on initial parameter values. These procedures are not appropriate for models where many parameters are included, and many local minima are supposed to exist because many solutions, including the global minimum, may be ignored. Evolutionary algorithms (EAs) (Eiben and Smith, 2003) can be used as a viable alternative to these problems. EA explore the complete parameter space by including, in a computational program, the mechanisms of reproduction, mutation, and the Darwinian principles of natural selection. Mimicking biological evolution, these programs iteratively generate better and better solutions by creating new generations of parameter estimates. Because EAs follow stochastic rather than gradient methods, the possibility for solutions to become trapped in local minima are lower (Moles et al., 2003). Because of the considerable number of parameters involved in our mechanistic model, an inhouse evolutionary algorithm (Roche et al., 2006) was used. We expect that a beta version of the program will be accessible to researchers in the near future (in the web page of some of the authors; http://servet.uab.es/biomathematics).
Statistics. Statistical comparisons between two groups (WT and mutated receptors) were performed by unpaired Student's t test, with the inclusion of Sidàk correction for multiple comparison tests. A value of p < 0.05 was considered statistically significant.
Results
Mathematical Modeling of mGlu Molecular Mechanism. Structural studies identified three states for the dimeric VFT domain: openopen (OO), closedopen (CO), and closedclosed (CC), which, if connected to the HD, give no, partial, and full activity, respectively (for review, see Pin et al., 2004). These different states were included in our modeling schema (Fig. 1). The distribution of these various states is governed by equilibrium constants according to a conformational induction approach (for discussion on conformational and selection approaches, see Supplemental Fig. 1). Thus, K_{1} and K_{2} are equilibrium dissociation constants for ligand binding to open states, whereas X_{i} and Y_{i} are induction constants for the closure of open states from either free or occupied receptors, respectively. As illustrated in Fig. 1, the equilibrium between inactive and active HD dimers is governed by a constant (L) depending on the state (OO, CO, or CC) of the VFT dimer.
After a parsimony procedure, we developed first the mathematical equations for ligand binding to the VFT domain; then, we included the connection to the HD to achieve fully a mechanistic description of receptor function.
Modeling the Binding of Ligands to the VFT Domain. The binding of a ligand A to the VFT receptor model depicted in Fig. 1 is described by eq. 1: where x = log[A], and x_{50} (the concentration of ligand for halfoccupation of the receptor binding sites) and x_{d}[determines the cooperativity and the shape (number of phases) of the curves] are empirical parameters depending on the equilibrium constants included in the mechanism (see eq. 7 in the Appendix).
To analyze ligand binding, we used eq. 1 following both empirical and mechanistic approaches. Three conditions can be considered for the x_{d} parameter, x_{d} = 0, x_{d} < 0, and x_{d} > 0. If x_{d} = 0, then f = 1/(1 + 10^{x}^{50–x}), and the ligand binds the dimeric receptor as it would do for a monomeric receptor (absence of cooperativity and monophasic curve). The other two conditions are related to the presence of cooperativity and then the possible appearance of biphasic curves (see below).
To illustrate the meaning of cooperativity in terms of eq. 1, Fig. 2 depicts the fractional binding for three ligands with a common x_{50} (–9) value and different x_{d} (0, –3, +3) values. Relative to the reference x_{d} = 0 (solid line), for x_{d} =–3, a steeper curve (long dashed line) in the midpoint is obtained, indicating positive cooperativity, whereas for x_{d} = +3, a flatter curve in the midpoint is obtained, indicating negative cooperativity (short dashed line). To compare the steepness of the f(x) curves, the Hill coefficient at the midpoint (n_{H50}) can be used (see Materials and Methods). In the simulations (see Fig. 3), and using eq. 8 in the Appendix (Giraldo, 2003), values of 1, 1.998, and 0.002 were determined for n_{H50} with x_{d} = 0, –3, and +3, respectively.
The Sign and Extent of Cooperativity and the Number of Phases of Saturation Binding Curves. The value of n_{H50} is commonly used for cooperativity characterization, and values of 1, greater than 1, and lower than 1 are attributed to absent, positively present, and negatively present cooperativity, respectively; this is consistent with our model (Fig. 4; paragraph above). It is worth noting that an n_{H50} greater than one is an indication of receptor oligomerization, but values equal or lower than one are also compatible with multiple receptor binding sites. The relationship between oligomerization and cooperativity has been discussed previously [for instance, see Park et al. (2002) for M_{2} muscarinic receptors; Strange (2005) for D_{2} dopamine receptors; Urizar et al. (2005) for glycoprotein hormone receptors; and Vilardaga et al. (2008) for the crosstalk between α_{2A}adrenergic and μopioid receptors (for review, see Milligan and Smith, 2007; Springael et al., 2007)]. In particular, it is worth mentioning the proposal that binding experiments can be used to prove the existence of receptor dimers, suggested on the basis of experiments on vasopressin and oxytocin receptors showing positive or negative cooperativity (Albizu et al., 2006).
The steepness, measured by n_{H50}, and the number of phases of saturation binding curves are not independent properties. As illustrated in Fig. 2, for ligands showing either a neutral or a positive cooperativity, a monophasic curve is obtained, whereas ligands having a negative cooperativity give a biphasic curve. This observation is further confirmed by the number of inflection points on these curves (Fig. 3). Although a single point of inflection at x =–9 (monophasic curves) for the neutral (Fig. 3A) and positive cooperative (Fig. 3B) ligands is observed, two additional points, at x =–5.7 and x =–12.3 (see Appendix eq. 9), are obtained for the negative cooperativity (Fig. 3C) ligand.
However, as shown in Fig. 4, the number of phases for ligands with negative cooperative binding can only be detected if x_{d} > log2 (Fig. 4, red line). The relationship between x_{d} and n_{H50} (eq. 8) reveals that for negative cooperativity (n_{H50} < 1), two phases can be detected (three inflection points) when n_{H50} < 0.67 but not for 0.67 ≤ n_{H50} < 1.
Accurate characterization of curve shapes is fundamental for the correct analysis of binding experiments. Equation 1, in its empirical form (two parameters, x_{50} and x_{d}), has proved instrumental for the detection of the sign and extension of cooperativity and the identification of the possible existence of two phases, in the case of negative cooperativity. It is worth noting that negative cooperativity can be accounted by other mathematical models as, for instance, the twoindependent sites model or a monovalent receptor interacting with a Gprotein. However, positive cooperativity cannot be explained without considering the receptor as multivalent (Mattera et al., 1985; Christopoulos and Kenakin, 2002; Albizu et al., 2006; Franco et al., 2006).
Cooperativity Binding Effects and the mGlu Molecular Mechanism. After exploring in detail using an empirical approach the effects of cooperativity on the saturation binding curves, the question arises about what mechanistic determinants control this property according to the VFT binding model represented in Fig. 1. As shown in the Appendix (eq. 7), the cooperativity parameter x_{d} was defined as x_{d} = log c_{1}/(2√_{c2),} where c_{1} and c_{2} depend on the equilibrium constants included in the binding mechanism. In fact, each of the cooperativity conditions can correspond to several sets of binding constants; accordingly, there is not a unique solution to describe positive, negative, or absence of cooperativity. Nevertheless, some particular cases are worthy of discussion:

K_{1} = K_{2}, Y_{5} = X_{1}, Y_{6} = X_{2}, Y_{1} = Y_{2} = Y_{5}, Y_{3} = Y_{4} = Y_{6}, x_{d} = 0, and the ligands show no cooperativity.

K_{1} > K_{2}, Y_{5} = X_{1}, Y_{6} = X_{2}, Y_{1} = Y_{2} < Y_{5}, Y_{3} = Y_{4} < Y_{6}, x_{d} < 0, and the ligands show positive cooperativity.

K_{1} < K_{2}, Y_{5} = X_{1}, Y_{6} = X_{2}, Y_{1} = Y_{2} > Y_{5}, Y_{3} = Y_{4} > Y_{6}, x_{d} > 0, and the ligands show negative cooperativity.
We see that absence of cooperativity can be obtained if the three horizontal branches of the VFT part in Fig. 1 have the same weight; positive cooperativity, if the lowest branch (doubly occupied receptor) has higher weight than the middle branch (singly occupied receptor); and negative cooperativity, if the lowest branch has lower weight than the middle branch. The above equilibrium constants relationships are qualitatively consistent with the cooperativity concept, but, to obtain a quantitative explanation of the mechanism and provide biologically meaningful parameters, an application to experimental data is needed.
Reanalyzing Experimental Binding Data with the VFT Mathematical Model. One study (Suzuki et al., 2004) examined experimentally the cooperativity between VFT binding sites upon ligand binding. The authors used the purified soluble VFT dimer of mGlu1 as a model, and four ligands [two agonists, glutamate and quisqualate, and two antagonists, (S)MCPG and LY367385] were investigated. The Hill analysis of the titration curves showed cooperativity only for glutamate (negative cooperativity). It is interesting to note that in the presence of calcium ions, the values of x_{50} decreased for agonists but not for antagonists, indicating that the positive effects of calcium ions on receptor affinity are specific for agonist binding and, therefore, are associated with the stabilization of the active (closed) conformation of the VFT protomers (both CO and CC states). In addition, the Hill coefficient of glutamate binding at the midpoint changed from 0.55 to 0.70 upon addition of calcium ions, revealing an effect on the cooperativity of glutamate binding. For quisqualate, Hill coefficients of 1.04 and 0.92 were obtained in the absence and in the presence of calcium ions, respectively, and a parallel left shift of the curve upon addition of calcium was observed (Table 1, experimental x_{50} and n_{H50} values in italics).
These results can be well accommodated within our model. Table 1 shows some combinations of parameter values compatible with the above experimental findings. Although the values are arbitrary in absolute terms, the relation between them follows a plausible pharmacologic criterion, thereby allowing for a quantitative exploration of mechanistic hypotheses. Glutamate in the absence of calcium ions was taken as a reference, and the following conditions were stated: 1) closure of VFT is rare in the absence of ligand (X_{1} = X_{2} = 10^{–4}); 2) binding to OO states involves negative cooperativity (K_{2} > K_{1}); 3) glutamate in an open VFT does not affect the closure of the unoccupied associated partner (Y_{1} = X_{1}); 4) as a positive agonist, glutamate in an open VFT domain induces its closure (Y_{2} = 5); 5) positive cooperativity was assigned to the closure induction of an occupied VFT from the previously closed neighbor (Y_{3} = 10^{3} > Y_{2}); 6) because of the same arguments, Y_{5} = 5 and Y_{6} = 10^{3}; and 7) the four constants Y_{1} to Y_{4} are linked, and Y_{4} = Y_{1}Y_{3}/Y_{2} = 2 × 10^{–2}.
It is worth noting that we distinguish between two cooperativity concepts: binding cooperativity and induction cooperativity. The former involves the binding to inactive (open) states, and the latter points to the induction of active (closed) states. In the context of the VFT system, we hypothesize that only a closed state can induce the closure of the associated VFT and the latter must be occupied to facilitate the process; such an hypothesis is supported experimentally (Kniazeff et al., 2004). In addition, negative binding cooperativity and positive induction cooperativity were proposed for glutamate. With the above constants, a Hill coefficient of 0.54 (using eq. 8) and an x_{50} of –5.35 were obtained (for these and ensuing results, see Table 1). These values are in agreement with the above experimental results (Suzuki et al., 2004). In the presence of calcium ions, closed states are stabilized. Thus, Y_{2}, Y_{3}, Y_{5}, and Y_{6} were augmented. The calculated Hill coefficient increases to 0.68 (eq. 8), and the curve shifts to the left (x_{50} =–6.05).
For quisqualate, first in the absence of calcium ions, we followed a similar rationale. We maintained the constants for the binding to open states (K_{1} and K_{2}) as in glutamate, but, due to the greater agonist capability of quisqualate, we increased the constants for induction of closed states (Y_{2}, Y_{3}, Y_{5}, and Y_{6}). The consequences on the pharmacologic curveshaped descriptors relative to glutamate were an increase of the Hill coefficient from 0.54 to 1 and a left shift of x_{50} from –5.35 to –6.1. A further increase of Y_{2}, Y_{3}, Y_{5}, and Y_{6} constants to account for the presence of calcium ions led to a further left displacement of the curve (x_{50} =–6.44) and a lowering of the Hill coefficient (n_{H50} = 0.90). We see that a concomitant increase of the induction constants for the lower (Y_{5}, Y_{6}) and middle (Y_{2}, Y_{3}) rows of the VFT model (Fig. 1) can yield different results on the Hill coefficient, depending on the weight of one row relative to the other.
Finally, the behavior of an antagonist (LY367385) has been simulated by giving the same values to the binding constants for the open states (K_{1} = K_{2} = 10^{–6}, absence of binding cooperativity) and low values for the constants for induction of closed states (Y_{2}, Y_{3}, Y_{5}, and Y_{6}). The parameter values yield a Hill coefficient of one and an x_{50} equal to –6. For this antagonist, increasing 2 orders of magnitude the values of Y_{2}, Y_{3}, Y_{5}, and Y_{6} constants, to account for the presence of calcium ions, had no effects neither in the Hill coefficient nor in the location of the curve along the xaxis.
These analyses clearly illustrate that our VFT model allowed a correct description of several experimental saturation binding curves. Such mechanistic analyses revealed important information regarding the functioning of these complex dimeric receptors. Essentially, agonists and antagonists were appropriately differentiated by their different propensity to induce the closure of the VFT binding sites relative to the basal state. A quantitative description of cooperativity was provided, which allowed, in addition, an analysis in terms of binding and induction cooperativity concepts. This formalism may explain some striking results as, for instance, the observed negative cooperativity of some agonists, which may be interpreted as the sum of negative binding (to open states) cooperativity and positive induction (of closed states) cooperativity, the latter being the characteristic that identifies a full agonist. However, to properly discuss agonist behavior, the transduction of binding into functional response is needed.
Modeling the mGlu Function. Transduction of VFT binding into functional response involves the activation of the HD. Both experimental data (Hlavackova et al., 2005; Damian et al., 2006) and theoretical calculations (Filizola et al., 2006) suggest that only one HD per dimer is in the active conformation at a time. This is also supported by the recent results by Bayburt et al. (2007) for the rhodopsin dimer. These authors show that indeed a single rhodopsin in a dimer can reach the active MII state at a time. Likewise, White et al. (2007) also provide interesting data showing the asymmetrical functioning of the purified NTS1 receptor dimer. The simplest model accomplishing this finding is an asymmetric twostate dimer model consisting of two identical protomers (RR) for the inactive and two different protomers (RR*) for the active dimer state, whose relative populations are governed by an equilibrium constant, say L (eq. 11).
In mGlu receptors the Gproteinactivating domain, HD, and the ligandbinding VFT are linked by the CRD, allowing the functional coupling between them (Rondard et al., 2006). Although agonistinduced closure of one of the VFTs is required to activate the HD, the closure of the two VFTs is necessary for full activation. In our model, the coupling between VFT and HD dimers was defined by assuming that the equilibrium constant L depends on the state of the VFT dimer, with L_{1}, L_{2}, and L_{3} for OO, OC, and CC, respectively, and where the L_{1} < L_{2} < L_{3} relationship is expected (see Fig. 1 and eq. 11). Thus, in our model, the active RR* HD is present in each of the VFT states but with a different propensity of formation. Furthermore, we did not explicitly assign the R or R* to a specific HD subunit related to the state (O or C) of its associated VFT. This is in agreement with recent findings (Brock et al., 2007) showing that VFT agonist stimulation involves an intersubunit rearrangement resulting in the activation of either HD (cisor trans) with the same efficiency.
The fractional functional response (f_{R*}) is defined as the fraction of receptor concentration in the active form (meaningbearing R* state). where [RR_{xyT}] and [RR*_{xyT}] stand for total receptor (VFTCRDHD) concentrations in inactive and active HD dimer states, respectively, and xy denotes either OO, CO, or CC VFT dimer states. The mechanistic definition of the a_{i} parameters can be found in Appendix eq. 12.
Assessing the Shape of ConcentrationEffect Curves. Quantitative characterization of the shape of the fractional response given by eq. 2 may provide useful information about the ligandreceptor interaction. Thus, using the transformation x = log[A], theoretical basal and maximal or minimal responses can be calculated as the left and right asymptotes of f_{R*}, respectively (eqs. 3 and 4).
The expressions for basal and maximal or minimal responses make sense because none of the liganddependent equilibrium constants appear in the basal response, whereas only constants for induction of active closed states for the doubly occupied receptor (Y_{5}, Y_{6}) appear in the maximal or minimal responses. It is worth mentioning that none of the constants involving the singly occupied receptor dimer appears in the expression for the right asymptote in agreement with an occupation of all receptor sites because [A] increases infinitely. Equations 3 and 4 measure the efficacy of the system in the absence and in the presence of the ligand, respectively (note the formal similarity between both equations). If we assume that the coupling constants L_{i} are not liganddependent, Y_{5} and Y_{6} are the only constants responsible for the intrinsic efficacy of a ligand. A definition of neutral antagonism, positive agonism, and inverse agonism can be obtained by a comparison between basal and maximal (or minimal) responses: Right = Left, Right = maximum > basal, and Right = minimum < basal, respectively. Note that for an inverse agonist, the term maximum changes to minimum. A neutral antagonist results by making Y_{5} = X_{1} and Y_{6} = X_{2}; a positive agonist, by making Y_{5} > X_{1} and/or Y_{6} > X_{2}; and an inverse agonist, by making Y_{5} < X_{1} and/or Y_{6} < X_{2}. For illustration, Table 2 shows some combinations of parameters yielding to either positive or inverse agonism.
The middle row of the VFT binding part in Fig. 1 (the induction of active closed states from single receptor occupation), which is not present in the definition of efficacy (maximal or minimal response – basal response), affects the potency (A_{50}) of the agonist (Appendix eq. 15). Agonist potency can be investigated for any ligand different from a neutral antagonist, where the left and right asymptotes take the same value (the denominator of eq. 15 is 0), and A_{50} is indeterminate. The ± sign in eq. 15 results for the possibility of A being either a positive or an inverse agonist. A systematic variation on the equilibrium constants included in the fractional response (eq. 2) can be found in the Supplemental Fig. 2). The collection of curves on display offers a mechanistic explanation for a broad range of concentrationeffect profiles, including full and partial agonism, inverse agonism, monophasic, biphasic, and bellshaped curves.
Reanalyzing Experimental Functional Data with the mGlu Mathematical Model. In the binding section of this study, experimental data were reanalyzed to encompass the theoretical pharmacological space within realistic limits. The same rationale is followed here, and a recent experimental work (Kniazeff et al., 2004) involving functional studies in a fulllength (mG5C1mG5C2) receptor using quisqualate as agonist was selected. In this study, the authors observed that mutating the VFT in the binding site (the socalled YADA mutant) led to a loss of agonistinduced activity of the receptor. However, when a single VFT was mutated in the dimer (the socalled mG5C1mG5C2YADA mutant), agonist binding in the WT VFT allowed agonist interaction in the mutated VFT to further increase receptor activity. The resulting doseresponse curve was biphasic, and two x_{50} values (x_{501} and x_{502} for the first and second phase, respectively) were identified. The x_{501} of quisqualate for the YADA heterodimer was close to that measured on the wildtype receptor, whereas the second response yielded a maximum of ∼80% of that measured with the control receptor. Similar biphasic shapes were obtained when performing some other mutations (the socalled YATA or SATA mutants). It is interesting to note that the x_{501} of quisqualate did not depend on the mutants used, whereas in contrast, the x_{502} largely varied among the mutants. Altogether, these results suggested that mutations damage the mechanism of agonistinduced closure of the mutated protomer and that the mechanism can be, albeit only in part, recovered by the activation of the associated WT subunit, with the first phase of the response curve resulting from agonist binding in the wildtype subunit only (Kniazeff et al., 2004).
Figure 5 shows the experimental data for WT (Fig. 5A) and singlemutated VFT (Fig. 5C) as solid circles. The Empirical Models 5 (Hill equation) and 6 (sumoftwo fractional Hill equations) were fitted to the former (Fig. 5A, solid line) and the latter (Fig. 5C, solid line) receptor systems, respectively.
For the WT VFT (Fig. 5A), the parameter estimates under eq. 5 were as follows: Basal = 0.46, Maximum = 0.99, x_{50} = –7.24, and n_{H} = 1.26 (for these and ensuing results, see Tables 3 and 4). It is worth mentioning the unusual very high value found for the basal response in this experiment, being the typical ones around 20 to 25% of the maximal response. The value, greater than one, for the Hill coefficient (n_{H}) suggests that positive cooperativity is originally present in the WT receptor if quisqualate is used as agonist. For the singlemutated VFT, the parameter estimates under eq. 6 were as follows: Basal = 0.41, Maximum = 0.89, p = 0.30, x_{501} =–6.93, n_{H1} = 1.33, x_{502} =–3.56, and n_{H2} = 1.09. After mutation, a decrease in the maximal response and a split of the curve into two phases are observed (Fig. 5C). As shown earlier (Kniazeff et al., 2004), the value for x_{501} is close to that measured on the WT receptor. It is interesting to note that this similarity is also observed for n_{H1}. Both properties indicate that the first response comes from quisqualate binding in the wildtype subunit only and that the responsegenerating machinery for the WT protomer is not significantly affected by mutation. As mentioned above, a demonstration for this comes from the different x_{502} values and a common x_{501} value obtained with different mutants (see Fig. 5d of Kniazeff et al., 2004). Finally, the decrease of both the maximal response and the second Hill coefficient suggests that the functional activity provided by the mutated protomer is partially reduced.
To provide a mechanistic interpretation of the experimental data, the WT curve was examined using the model developed for mGlu function (eq. 2). The high number of parameters included in the equation and the correlation possibly existing between them preclude the use of classic gradient nonlinear fittings. Accordingly, a stochastic evolutionary algorithm was used (see Materials and Methods). To do this, we chose first a reference state by taking the values listed in Table 1 for the quisqualateVFT interaction in the absence of calcium together with plausible values for the VFTHD allosteric constants (basically, the L_{1} < L_{2} < L_{3} relationship). To assure sufficient sampling of parameter population, a ±3 interval was chosen for each optimized parameter, and 100 independent runs were performed. Table 3 shows the mean and S.D. of each of the parameters. Essentially, we found negative cooperativity for the binding to the OO sites (K_{1} < K_{2}), positive cooperativity for the induction of closed states both in singly (Y_{2} < Y_{3}) and in doubly (Y_{5} < Y_{6}) occupied VFT dimers, and the expected L_{1} < L_{2} < L_{3} relationship for the transduction of occupation into response. Figure 5B shows the theoretical curve produced by that run of 100 whose parameters are closest (Euclidean distance) to the mean values. Visual comparison with the curve produced by the empirical onesite model (Fig. 5A) indicates a similar fitting, which is confirmed by the x_{50} and n_{H50} curveshaped descriptors (Table 4).
For the mutated VFT, we used the values inferred for the WT receptor as a starting point, and our evolutionary algorithm was employed to obtain the optimized parameters. For simplicity, we considered first that only some of the constants, those directly associated with the mutated protomer, could change after mutation; that is, K_{2}, Y_{3}, and Y_{6}. Furthermore, because the mutation of the VFT of one protomer may alter the ability of the dimer to transmit the signal, L_{1}, L_{2}, and L_{3} were allowed to change as well. A systematic analysis was performed to identify the equilibrium constants mainly affected by the mutation; successive independent fittings were carried out ranging from a single optimized parameter to all possible combinations of parameters (that is, from two to six optimized parameters). In all cases, a ±3 interval was allowed for each of the parameters varied. It is interesting to note that when only one optimized parameter was contemplated in the fitting, three of them clearly differentiated from the others in the capacity of reducing the sum of squares of the error: L_{3}, Y_{6}, and K_{2}, in decreasing order. Moreover, the combinations of two parameters leading to best fittings were (L_{3}, K_{2}), (L_{3}, Y_{6}), and (K_{2}, Y_{6}). These results indicate that the main effects of VFT mutation on functional response are as follows: 1) an impairment of the allosteric VFTHD interaction associated to the CC state (L_{3}), 2) a diminution of the capacity of the first closed protomer (supposedly the WT) to induce the closing of the second protomer (supposedly the mutated) in the doubly occupied VFT (Y_{6}), and 3) a decrease in the affinity of the ligand for the second binding site (K_{2}).
To test statistically the above hypotheses on parameter modification after receptor mutation, 100 independent runs were performed for the mutated receptor, with all fit parameters free and a ±3 interval for each of the parameters (Table 3). Statistical comparison between WT and mutated receptor by Student's t test, including Sidàk correction for multiple comparisons, confirmed the mechanistic proposals suggested above. As expected, changes on K_{2}, Y_{6}, and L_{3} reached statistical significance. The other parameters that significantly changed on receptor mutation were L_{1}, Y_{2}, X_{1}, and X_{2}. As can be seen from eq. 3, L_{1} ≈ Basal/(2(1 – Basal)) if X_{1} and X_{2} are much lower than 1. Then, the change on L_{1} reflects the observed change on the basal response. The model predicts a half lowering of Y_{2} after mutation, which suggests that the mutant protomer hampers the closing of the occupied WT. This indicates that a mutation on the recognition site of one protomer affects the binding and activation capacity of the mutated unit but also the intrinsic efficacy of the WT partner. The lower values for X_{1} and X_{2} on the mutated relative to the WT receptor are consistent with above findings. Finally, to illustrate the goodness of fitting of our mechanistic approach, Fig. 5D shows the theoretical curve produced by that run of 100 for which the parameters are the closest (Euclidean distance) to the mean values. The apparent resemblance with the fitting produced by the empirical twosite model (Fig. 5C) is quantitatively confirmed by the location and Hill coefficients parameters (Table 4). However, only the mechanistic model provides a detailed view of the important steps of the activating process affected by the mutations.
Functional Dynamics of WT and Mutated Receptors: The Distribution of VFT States upon Ligand Binding.Figure 6 shows the relative distribution of VFT states for both WT and mutated receptors (Kniazeff et al., 2004) using the same mechanistic constants values as in Fig. 5, B and D. We see that WT and mutated receptor profiles show similarities and differences. The maximal asymptote as [A] decreases (on the left) corresponds to the free open (OO) state (solid red line). The fractional concentration of this state ranges between 1 and 0, reaching the asymptotic minimal value for lower [A] values in the case of the WT. Two bell curves appear in both cases, which correspond to singly occupied O_{A}O (long dashed red line) and C_{A}O (longdashed blue line) VFT states. These curves are broader for the mutated receptor, indicating that the importance of these states spans over a longer [A] range for this receptor genotype. It is interesting to note that the O_{A}O state, which is slightly present in the WT, contributes significantly to the VFT distribution in the mutated receptor. The maximal asymptotes as [A] increases (on the right) correspond to doubly occupied states. However, although in the case of the WT all the doubly occupied receptors are in the fully active C_{A}C_{A} form (short dashed green line), a distribution of states is obtained for the mutated receptor, namely, fractional 0.75, 0.20, and 0.05 values for the fully active C_{A}C_{A} (short dashed green line), partially active O_{A}C_{A} (short dashed blue line), and slightly active O_{A}O_{A} (short dashed red line), respectively, are found. In addition, the significant presence of C_{A}C_{A} starts at lower concentrations for the WT. Comparison between Figs. 6B and 5D shows that the functional intermediate plateau of the mutated receptor is produced mainly by the C_{A}O state because it is both more abundant and more efficacious (L_{2} > L_{1}) than the O_{A}O state. Figure 6B helps understand the functional dynamics of mutated receptor. Basal response is due to the OO state; accumulative addition of A leads first to single occupation, which can be in either O_{A}O and C_{A}O forms, where mainly the WT protomer is occupied in the singly occupied mutated receptor; the relevant importance of the O_{A}O species compared with the WT receptor is an indication that mutation of a protomer affects also the propensity of closure of the WT partner, and increasing further [A] leads to doubly occupied receptors in the order [C_{A}C_{A}] > [O_{A}C_{A}] > [O_{A}O_{A}], the two latter states being concentrations nonnegligible in contrast to the WT receptor.
The Cooperativity Issue and the Ligand Recognition Mechanism: The Application of the Model to Prediction of Ligand Properties. As quantitatively shown above, the combination of two states (open and closed) and domain dimerization provides the receptor with a flexible mechanism for ligand recognition and signal transduction. In the present article, two conceptually different cooperativity classes for the mGlu receptor have been proposed, binding and induction, the former related to binding to inactive VFT open states, and the latter related to induction of active VFT closed states. The model presented herein suggests that, if we assume three categories, absent (0), positive (+), and negative (–), for both binding and induction cooperativities, nine types of ligands could, in principle, exist from the {binding, induction} combination of cooperativities. One of the combinations {–, +} was identified for glutamate and quisqualate full agonists (Table 1). The sign of the cooperativities indicates that although the binding of the first ligand to the OO state diminishes the binding of the second ligand to the same state (negative binding cooperativity), the closure of a VFT protomer facilitates the closure of the second (positive induction cooperativity). It is interesting to note that other combinations, as for instance {+, –}, are theoretically possible and would yield to a potent partial agonist, i.e., a ligand with lower maximal response but a significantly leftshifted functional curve relative to a more efficacious full agonist. Figure 7 displays the concentrationeffect curves of two agonists with opposite cooperativity effects. One agonist (solid curve) poses negative binding cooperativity (K_{1} = 10^{–5}, K_{2} = 10^{–3}) and positive induction cooperativity (Y_{2} = Y_{5} = 10, Y_{3} = Y_{6} = 10^{2}), the other agonist (dashed curve), on the contrary, bears positive binding cooperativity (K_{1} = 10^{–5}, K_{2} = 10^{–7}) and negative induction cooperativity (Y_{2} = Y_{5} = 10, Y_{3} = Y_{6} = 10^{–2}). Comparison of the curves shows a slight leftshift displacement including an increase of the slope and a notable decrease of the maximal response of the latter curve relative to the former one, converting a full agonist into a partial but not less potent (in the common response range) agonist. Structureactivity studies identified the socalled APTC, a new family of mGlu orthosteric ligands (Schann et al., 2006). Among the series of synthesized ligands, one (FP0429) was shown to be a full mGlu4 agonist and a partial mGlu8 agonist. Because these receptors are highly homologous, the question on the structural features responsible for ligand specific properties arose. Sitedirected mutagenesis of two residues within the binding site that differ between mGlu4 and mGlu8 switched the agonist definition of FP0429 from full to partial agonist and vice versa (Frauli et al., 2007). It is interesting to note that docking analysis in mGlu VFT model attributed a role in agonist binding to one of the residues, whereas the other participated in the closure, thus the functional activity (Frauli et al., 2007). This explanation is consistent with our conceptual distinction between binding and functional (induction) cooperativities, suggesting that approaches based on the cooperativity issue may be of interest in the drug discovery process. To this end, and with the aim of helping chemists in their structureactivity analysis, Table 5 shows the association of the mechanistic parameters of the model (the equilibrium constants of Fig. 1) with the typical pharmacological properties basal response, efficacy, potency, and slope (the latter property with reference to biphasic and bellshaped curves).
Discussion
In the present study, a mathematical modeling of mGlu binding and function was made by constructing two models, one for the VFT and the other for the HD dimer domains. For simplicity, the analysis of ligand binding to VFT orthosteric sites was conducted without the inclusion of HD. A detailed exam of the shape of saturation binding curves in terms of the structure of the VFT lobes was performed, allowing the differentiation between positive and negative cooperativity and the description of biphasic curves. For the analysis of mGlu function, a coupling constant L, which embodies the probability of yielding an active HD dimer, was included. This constant modulates the population of active HDs upon selective binding to VFT states. According to experimental data, only one HD is activated at a time within a dimer. Thus, R*R* was precluded, and an asymmetric twostate dimer model (RR,RR*) was proposed.
The model allowed for the quantification of agonist efficacy and potency and the interpretation of pharmacologic curve profiles in mechanistic terms. In addition to theoretical simulations (supplemental data), two published experiments, one involving binding (Suzuki et al., 2004) and the other functional (Kniazeff et al., 2004), were satisfactorily reanalyzed. An important outcome of the analysis was the mechanistic distinction between binding and induction cooperativities, the former related to the affinity of the second ligand to the OO state with respect to the binding of the first molecule and the latter to the induction of closure of the second protomer after closing of the first one. The reanalysis of the binding study (Suzuki et al., 2004) provided evidence that full agonists are characterized by positive induction (of closed states) cooperativity, which, depending on the sign and magnitude of the binding (to open states) cooperativity, may lead to observed negative, null, or positive cooperativity. The reanalysis of the functional study (Kniazeff et al., 2004) identified the binding/transduction parameters that were mainly affected by receptor mutation and the different functional dynamics (VFT states distribution) of WT and mutated receptors upon agonist concentration. Moreover, based on the conceptual distinction between binding and induction cooperativities, the functional profile of a theoretically potent partial agonist was obtained, which suggests a possible application of the model in structureactivity studies.
The model also quantitatively illustrates some of the advantages for a receptor of being a dimer. A receptor failure that would make a receptor fully inactive in the case of a monomer receptor can be partially compensated by an associate protomer in the case of a dimer receptor (Kniazeff et al., 2004). In addition, the interdependence between binding sites makes a dimer receptor more efficient than the sum of two single monomers. First, positive induction cooperativity facilitates the closure of the second site, increasing the efficiency of the system. Second, negative binding cooperativity for slightly active open sites (OO) biases the receptor sites distribution toward partially active (CO) and fully active (CC) receptor sites.
It has been suggested that, in addition to the conformation of each HD within the dimer, the relative positioning between the heptahelical protomers plays an important role in signal transduction. Thus, at least two conformations for the active RR* dimeric state, (RR*)_{a} and (RR*)_{b}, with a probability of occurrence depending on the VFT state, are conceivable. This variety of active conformations could explain the multiplicity of pathways associated to VFT activation. It has been found that only CC leads to full activation of G_{q}, whereas the CO state leads both to G_{s} coupling and to partial G_{q} activation (Kniazeff et al., 2004; Tateyama and Kubo, 2006). This level of detail, which might be useful in biochemical experiments involving more than one Gprotein, has not been considered necessary for the purposes of the present study, where only one conformation for the RR* active state was included.
It is remarkable that by taking an agonist as a reference, a partial but potent virtual agonist was devised by our model by inverting the binding and induction cooperativities of the reference ligand. In this regard, we expect that the mechanistic components included in the model may help chemists in their structureactivity studies to quantify the effects of drugs. A paradigmatic related example can be found in the discovery of burimamide (Black et al., 1972), the first H_{2} receptor antagonist, which was obtained by taking the structure of the endogenous agonist histamine as the chemical starting point and progressively removing its agonist properties along a structureactivity pathway including partial agonists as signposts.
The inclusion in the model of VFT and HD domains leaves the model ready for the analysis of allosteric compounds, which, by binding to the HD, modulate the binding and/or function of orthosteric compounds acting on the VFT domain. In this regard, we expect the model to be able to account for the effects of PAMs and negative allosteric modulators on constitutive activity, agonist efficacy, agonist affinity, etc. Moreover, a quantitative explanation of the functional differences between Ca^{2+} and Gd^{3+} on these receptors represents another challenge to the model. The analysis of experiments involving these and other molecular interventions as, for instance, the molecular design of antagonists and inverse agonists by altering the cooperativity properties as suggested above, will allow not only check the validity of the model but open new possibilities to tune its parameters. This will be the subject of further work. Finally, the results suggest that our model, initially conceived for mGlu receptors, may apply to any other receptor system composed of an extracellular agonist binding domain and a transmembrane functional domain.
Appendix
The Fractional Binding Function for the VFT Domain
where the total concentration of receptors, [R_{t}], and the concentration of bound ligand, [A_{bound}], are defined as follows: and
Because the f binding function was defined relative to the total number of binding sites (2[R_{t}]), f ranges between 0 and 1. It is worth noting that eq. 7, expressed as a function of c_{1} and c_{2} parameters, is the same as that obtained for the socalled twostate dimer model (Franco et al., 2006), and, for this reason, the x_{d} parameter is related to a cooperativity index as empirically defined in the previous model (Casadó et al., 2007). However, the definition of the parameters is different here because, in the present model, the states of the protomers (either closed or open) within the dimer molecule can be distinguished.
The Hill Coefficient at the Midpoint
The Hill coefficient at the midpoint, n_{H50}, is related to the first derivative at the midpoint, (df/dx)_{x}_{50}, by (see Materialsand Methods):
Determining the Number of Phases of Sigmoid Curves
The number of phases of a sigmoid curve depends on the number of points of inflection. Monophasic curves show one point of inflection whereas biphasic curves present three. Mathematical analysis of eq. 7 shows that, in general, the number of points of inflection depend on the value of x_{d} relative to log 2. For x_{d} ≤ log 2, there is one point of inflection at x = x_{50}; for x_{d} > log 2, there are three points of inflection, one at x = x_{50} and the other two at
The Asymmetric HD Activation Model
Equation 10 shows the equilibrium between inactive (RR) and active (RR*) HDs. where R and R* stand for the inactive and active HDs within the dimer, respectively, and 2L = ([RR*]/[RR]) is the macroscopic equilibrium constant for the equilibrium between HD dimer states.
In our model, it was assumed that the equilibrium constant L depends on the state of the VFT domain. Thus, three apparent constants for the equilibrium between inactive (RR) and active (RR*) HDs are defined: where {, , and } and {[RR_{ooT}], [RR_{ocT}], and [RR_{ccT}]} stand for total active and inactive HD connected to OO, OC, and CC VFT, respectively.
The Fractional Functional Response
The fractional functional response (f_{R*}) is defined as the fraction of receptor concentration in the active form: where and a_{i} = c_{i}/c_{6} for i = 1 to 6 using the equilibrium constants depicted in Fig. 1 and eqs. 10 and 11.
Pharmacological Descriptors of Functional Response Curves
Using the transformation x = log[A], the theoretical basal response is calculated as the left asymptote of f_{R*}, and the theoretical maximal or minimal response is obtained as the right asymptote of f_{R*}.
The potency (A_{50}) of the agonist is calculated as [A] for Response = Left + [(Right – Left)/2]. where a = a_{1} – a_{3}a_{4}; b = a_{3}a_{4}a_{5} –2a_{2}a_{4} + a_{1}a_{5}; and c =–a_{4}(a_{1} – a_{3}a_{4}); and the ± sign in eq. 15 results for the possibility of A being either a positive or an inverse agonist.
The sensitivity of the receptor to an increment in the agonist concentration is measured by the first derivative of the receptor function (eq. 16).
The Hill coefficient at the midpoint can be calculated from eq. 16 (Giraldo et al., 2002; Giraldo, 2003).
Acknowledgments
We thank Anna Castellano for technical assistance.
Footnotes

This study was supported in part by the Ministerio de Educación y Ciencia, Spain (Research Grant SAF200406134).

Article, publication date, and citation information can be found at http://jpet.aspetjournals.org.

doi:10.1124/jpet.107.133967.

ABBREVIATIONS: GPCR, Gproteincoupled receptor; mGlu, metabotropic glutamate; VFT, Venus flytrap; HD, heptahelical domain; CRD, cysteinerich domain; PAM, positive allosteric modulator; EA, evolutionary algorithm; WT, wild type; OO, openopen; CO, closedopen; CC, closedclosed; (S)MCPG, (S)αmethyl4carboxyphenylglycine; LY367385, (+)2methyl4carboxyphenylglycine; OC, openclosed.

↵ The online version of this article (available at http://jpet.aspetjournals.org) contains supplemental material.
 Received November 7, 2007.
 Accepted February 19, 2008.
 The American Society for Pharmacology and Experimental Therapeutics