Identifiability of equilibrium constants for receptors with two to five binding sites

The function of receptor proteins embedded in cell membranes is controlled by the binding of ligands to highly specific binding sites. Following ligand binding, the receptors exert conformational changes, highly specific as well, that evoke secondary responses in either the membrane or the cell interior. The most simple-structured ligands are ions, e.g., Ca²⁺ ions, but most ligands are organic molecules, including highly diverse neurotransmitters, or peptide- and proteohormones.

In ionotropic receptors, the binding of ligands controls the open probability of a pore for ions. To realize a pore structure, many of these receptors are oligomers composed of between two and five subunits. These subunits can be either identical or different, resulting in homomeric or heteromeric receptors, respectively. The oligomeric nature of the receptors is what causes the number of binding sites to be predominantly larger than one. The consequence is that the activation of the receptors is controlled by the binding of several ligands and the resulting conformational changes by each binding step influence each other, resulting in the phenomenon of cooperative receptor activation. Moreover, the cooperative conformational changes of the subunits feed back to the binding steps themselves in a reciprocal fashion (Colquhoun, 1998; Kusch et al., 2010). This complexity of processes severely limits quantitative studies on receptor activation.

For >50 yr, the allosteric Monod–Wyman–Changeux (MWC) model (Monod et al., 1965) has been used to interpret the activation of symmetric proteins composed of identical subunits. The applicability of this model is based on the simple elegance of the assumption that the subunits perform a joint “allosteric” conformational change of all subunits and that the equilibrium of this allosteric step is shifted in proportion to the number of ligands bound. To keep the models manageable, often-defined stoichiometric factors have been used. This leads to an extremely low number of required parameters for describing the activation process, one equilibrium association constant, one equilibrium constant for the allosteric conformational change, and one allosteric factor. The MWC model has been widely applied, including various ionotropic membrane receptors (e.g., GABA_A receptors [Steinbach and Akk, 2019], nicotinic acetylcholine receptors [nAchRs; Lee and Sine, 2005], and also cyclic nucleotide-gated [CNG] channels [Goulding et al., 1994]) that are gated from the intracellular side. Moreover, the MWC model was extensively used for quantifying the oxygen binding in the tetrameric hemoglobin molecule (Eaton et al., 1999), though it is composed of different types of subunits, two α-chains and two β-chains in the case of adult hemoglobin (Fermi et al., 1984). However, the simplifying assumptions of using only a few parameters and fixed stoichiometries often severely limited detailed conclusions on the interaction of the binding sites and subunits in the framework of the MWC model, either for current measurements (Colquhoun, 1998) or binding measurements (Middendorf and Aldrich, 2017a, 2017b), even when applying state-of-the-art Bayesian strategies (Hines et al., 2014). For heteromeric receptors, the identifiability becomes even worse because any analysis would require more than one equilibrium association constant. Some remedy to overcome these problems in ligand-gated ion channels (LGICs) was to study subunit cooperativity by sophisticated single-channel analyses of burst activity, e.g., in homopentameric (Beato et al., 2007) or heteropentameric glycine receptors (GlyRs; Burzomato et al., 2004). An elegant alternative single-channel approach was to coexpress subunits with impeded binding sites and a second mutation in the pore to identify the subunit composition. This approach has been applied to homotetrameric CNG (Liu et al., 1998) and BK channels (Niu and Magleby, 2002), activated by intracellular cGMP and Ca²⁺, respectively.

Recently, we quantified activation gating in heterotetrameric CNG channels composed of two CNGA2 subunits, one CNGA4 and one CNGB1b subunit, by systematically impeding the four binding sites and subjecting the resulting 16 concentration–activation relationships (CARs) to a global fit analysis of 16 complex respective kinetic schemes (Schirmeyer et al., 2021). This approach allowed us to describe receptor activation by 32 equilibrium association constants for the 4 binding sites and 4 closed–open isomerization constants. Consistency of the subunit composition was reached by concatenating the subunits. The high degree of complexity could be managed by transposing the models to 4-D hypercubes. Due to multiple microscopic reversibilities and by equating the constants of the two A2 subunits, this fit required only 17 free parameters, which is, however, still a formidable number. Nevertheless, the precision of our fit was surprisingly high.

We therefore studied on theoretical grounds how this unexpectedly good identifiability of parameters could be reached (Benndorf et al., 2022). We showed that with this strategy indeed multiple functional parameters of the individual subunits can be determined because the constraints in this approach are sufficiently high. Furthermore, we employed a stochastic fit approach using scaled unitary (SU) start vectors to avoid any bias arising from a specific start vector. We showed that 23 parameters could be determined reliably, including 15 virtual equilibrium-association constants, 4 equilibrium constants for the closed–open isomerizations, and 4 disabling factors for the mutations. The 15 virtual equilibrium-association constants allowed us to determine the 32 equilibrium-association constants specific to the subunits of a tetrameric channel at each degree of ligand binding.

It was therefore intriguing to generalize our strategy to receptors with lower and higher numbers of binding sites. Biologically, the incidence of receptors with four subunits and four binding sites seems to dominate over the other constellations. For ionotropic receptors, typical tetrameric examples with four binding sites are AMPA receptors (Sobolevsky et al., 2009) or CNG channels (Li et al., 2017). Representative trimeric receptors with three binding sites are purinergic P2X receptors (Kawate et al., 2009) or acid-sensing ion channels (ASICs; Jasti et al., 2007). Pentameric ionotropic receptors belong mainly to the class of cys-loop receptors, including the families of nAchRs, GlyRs, and serotonin receptors (5HTRs). Many of these receptors are formed by two α-subunits and three other homolog subunits. For example, the classical nAChR of both the electric eel and the electric fish (Changeux et al., 1970; Miledi and Potter, 1971) is composed of two α-subunits, one β-subunit, one γ-subunit, and one δ-subunit (Weill et al., 1974; daCosta and Baenziger, 2013), forming only two binding sites between the adjacent subunits α and γ as well as α and δ. On the other hand, there are also homopentameric relatives with five binding sites, e.g., the α7 nAchR (McKay et al., 2007) or the α1 GlyR (Zhu and Gouaux, 2021). The key determinant for our analysis is not the number of subunits but the number of functional binding sites.

Herein, we elaborate for channels of known composition with two to five binding sites how promising it is to quantify receptor activation in great detail by our strategy, to impede the affinity of defined binding sites, and to globally fit all possible CAR combinations with either wild type or impeded binding sites by complex systems of coupled kinetic schemes. In other words, we empirically study the identifiability of multiple equilibrium constants for synthetic parameter sets by globally refitting multiple simulated CARs. We demonstrate that complex models with more binding sites and possible CARs can provide more accurate and precise results than simple models with less binding sites and possible CARs.

To introduce the terminology used herein, we first consider the heteromeric allosteric (HA) model for the tetrameric channel as described previously (Schirmeyer et al., 2021) and which is also shown in its “systematic representation” in Fig. S1 A with respect to the number of binding sites. The model contains 32 different equilibrium association constants for the closed states, K_x, describing the binding of the ligand to the four binding sites for each occupation pattern of the other binding sites, and the 4 closed–open isomerization constants, E_x, for the transitions of the single to the quadruple liganded channel, thereby assuming that E_x is equal at an equal degree of liganding. Notably, due to microscopic reversibility, all equilibrium association constants for the open channel are implicitly given by the respective K_x of the closed channel and the two related E_x. At zero ligand, P_o was set to the negligible value 10⁻⁵. A tetrameric HA model is termed in the following as the “HA⁴” model.

For the previous global fit analysis, mutational factors, fd_x, were introduced for each mutated binding site, describing the reduction of the affinity. With all combinations of wild type and mutated binding sites, this led to maximally 16 CARs and corresponding 16 specific HA⁴ models. These 16 HA⁴ models build together the complex heterotetrameric allosteric CHA⁴ model which was finally subjected to a global fit. This terminology is used here throughout. The superscript indicates the number of binding sites.

Notably, the multifaceted network of binding steps in Fig. S1 A could be elegantly analyzed by transforming it into a 4-D hypercube (Fig. S1 B) and introducing 15 independent virtual equilibrium association constants, Z_x (Fig. S1 C), determining the 32 K_x by simple ratios (Schirmeyer et al., 2021). This representation of a model will be termed herein as “cube representation” and will also be used for channels with two to five dimensions.

Because the number of binding sites can deviate from the number of subunits and to generalize our approach, we categorize the LGICs herein by the actual number of binding sites and not by the number of actual subunits. Hence, the superscript at the end indicates the number of actual binding sites. For example, a tetrameric CNG channel with four binding sites would be described by a CHA⁴ model and a pentameric nAchR with two binding sites by a CHA² model.

For the LGICs with two to five binding sites, four types of respective HA models are considered herein. To keep track, we use the nomenclature “HA^x” model, where x denotes the number of ligands. To include the effects of the mutations, each mutation generates its own HA^x model for each CAR. The number of HA^x models included in a complex CHA^x model is given by n_CAR,max (Eq. 21) and is found to be 4, 8, 16, and 32 for x = 2, 3, 4, and 5, respectively.

For all computations, we assumed that microscopic reversibility (Colquhoun and Hawkes, 1995) is obeyed, resulting in the ligand binding to a defined number of independent equilibrium constants on which the others depend (Colquhoun et al., 2004).

In our previous analysis on heteroterameric channels, we eased the considerations by interpreting the 16 closed states of our HA⁴ models as corners of a 4-D hypercube in which the transitions, specified by the 32 K_x, form the 32 edges (Fig. S1 B; Benndorf et al., 2022).

Herein, we essentially follow this logic: given state 0, the 4-D hypercube defines 15 virtual equilibrium association constants (violet lines Z₁–Z₁₅ in Fig. S1 C), which allow the computation of the 32 K_x by respective ratios of two Z_x (see below). Microscopic reversibility is directly obeyed. For disabling the binding sites, we assume that the affinity of a binding site is reduced by a factor fd_x < 1 (x = 1…4). As in our previous approach, the fd_x factors were considered to affect solely the affinity at a binding site for the ligand. For the equilibrium constants of the closed–open isomerizations, E_x (x = 1…4; Fig. S1 A, bottom), it is assumed that they are equal at an equal degree of liganding and independent of where the binding site is occupied. E₀ is set to the negligible value of 10⁻⁵. The total number of free parameters for the global fit of the 16 concatamers is then 15 + 4 + 4 = 23.

To explain the rationale of the analysis, it is useful to skip to our earlier nomenclature. The closed state 0 is then termed C0000, where the four zeros tell that no binding site is occupied. The other 15 closed states in the 4-D hypercube are then given by Cijkl, where i, j, k, and l can be either 0 (empty) or 1 (occupied).

The ligand concentration L has the power of the sum of the indices α = i + j + k + l, corresponding to the number of bound ligands. The factors of disabling by the mutational factor fd_u (u = 1…4) of the four binding sites have the exponents γ1 to γ4, which are also either 0 or 1. An exponent is 1 if the binding site u is mutated and has bound a ligand. Otherwise, the exponent is 0.

This equation was used for both simulating the data points and computation of the CARs by the respective model in the fits.

Since the formula requires a value for E₀ though no detectable P_o was assumed with no ligand bound, we set E₀ generally to the negligible value of 10⁻⁵.

The nlinfit routine of the Matlab software was used to globally fit the CHA models, operating with an iterative least square estimation. The maximum iteration number was set to 100. The routine provides the fit parameters, x_i (E_x, Z_x, fd_x), the covariance matrix, Cov_EZfd, with the elements σ_ij of which the main diagonal gives, with i = j, the values of the variances and standard deviations $σi2$ and σ_i of the fit parameters, respectively. The Z_x allow to compute all respective equilibrium association constants, K_x, by simple ratios (Tables S1, S2, S3, and S4).

These relative standard deviations are useful for comparing the uncertainty of different parameters. Here, two different measures are used to quantify the uncertainties of the parameters in the fit, imprecision and inaccuracy (Benndorf et al., 2022).

The term precision of a parameter x_i specifies its quality, without knowing that the parameter has been determined correctly. Consequently, the term imprecision, ipr_i, denotes a statistical variability of a parameter with respect to its determined mean. Hence, ipr_i is given by σ_i,rel.

Computing 1,000 of the most challenging fits with stochastically varied scaled unitary start vectors lasted about 10 min for the least complex CHA² model and about 300 min for the most complex CHA⁵ model (Fujitsu Lifebook U series; vPRO i7).

Fig. S1 shows the structure of the HA⁴ model in systematic representation and as 4-D hypercube as published earlier. Fig. S2 illustrates by bar graphs how the fit parameters were set for the HA² to HA⁵ model. Fig. S3 shows the inaccuracy and imprecision of the parameters for 1,000 global fits of the CHA² to CHA⁵ model using a start vector with the true values. Fig. S4 provides color-coded correlation matrices for the CHA^x models. Fig. S5 illustrates the effects of mutational factors on the identifiability of E_x for the CHA^x models. Fig. S6 shows the effects of weakening the mutational factors on the identifiability of the Z_x for the CHA^x models. Fig. S7 shows plots of the inaccuracy and imprecision of the parameters for equal SU start vectors with stochastically varied elements. Fig. S8 illustrates successful_pΔMSE fits with different stochastically varied SU start vectors and stochastically varied elements. Fig. S9 provides histograms of successful_pΔMSE fits with different SU vectors and stochastically varied elements. Fig. S10 shows the fit parameters determined with different SU vectors and stochastically varied elements and compares them with the true parameters. Fig. S11 illustrates the inaccuracy and imprecision of the parameters for successful_pΔMSE fits with different SU vectors and stochastically varied elements. Tables S1, S2, S3, and S4 provide the generation of K_x from Z_x in the HA² to HA⁵ model, respectively. In the supplemental text, the algorithms for specifying the parameters for the CHA^x models are described.

The respective values are 3, 7, 15, and 31 (Fig. 3). The n_Zx values also tell directly the number of independent parameters in a fit for ligand binding. The Z_x define all respective K_x (Tables S1, S2, S3, and S4).

A plot of the ratio n_p/n_CARmax as a function of the dimension n_D shows a decrease in the higher dimensions (Fig. 4). Hence, if all possible CARs are included, this suggests that the chances for a reasonable fit are poorest in the 2-D case and largest in the 5-D case, i.e., a fit with seven parameters would be less promising than that with as much as 41 parameters, respectively, although it is a priori not trivial to conclude this because of the different structures of the models. We therefore investigated the outcome of global fits with the CHA² through CHA⁵ model systematically using the maximum numbers of CARs given by Eq. 21.

Comparison of CHA models with different structures and number of ligands requires the use of related theoretical parameters. These theoretical parameters were generated in some relation to the previous parameters for the heterotetrameric CNG channels (Schirmeyer et al., 2021; Benndorf et al., 2022). Their exact derivation is described in the supplemental text at the end of the PDF, and their values are illustrated in Fig. S2.

For E_x, we set the value of the single-liganded channel to 10⁻¹ and that of the fully liganded channel to 10², i.e., we assumed that maximal liganding causes a 1,000-fold increase of E_x with respect to E₁. The intermediate E_x values were set equidistantly on a logarithmic scale between E₁ and the maximum E_x according to Eq. 23.

For fd_x, the situation is grossly analogous. fd₁ was set to 10⁻⁴ and that of the maximally impeded binding site was set to 10⁻², resulting in a maximum difference between the binding sites of 100. Again, we set the intermediate fd_x values equidistantly on a logarithmic scale between fd₁ and the maximum fd_x according to Eq. 24.

For Z_x, the situation is more complex because the number of Z_x, n_Zx, rises exponentially according to Eq. 19 and the options at an equal degree of liganding follow binomial statistics. The generation of the parameters are described by Eqs. 25, 26, 27, 28, 29, 30, 31, 32, 33, and 34.

While the values for the E_x and fd_x differ substantially for different x, quite similar to values found for heteromeric CNG channels (Schirmeyer et al., 2021; Benndorf et al., 2022), we chose a more challenging condition for the Z_x values to differ at an equal degree of liganding by only one order of magnitude. Moreover, they also superimpose significantly for different degrees of liganding.

With these equilibrium constants, we computed families of data points with the CHA², CHA³, CHA⁴, and CHA⁵ model by including the respective maximum numbers of CARs and using Eq. 8 (Fig. 5). The noise was set to 0.25 (see Materials and methods), a value resembling the scatter in reasonable experimental measurements.

According to the set parameters for respective wild type and maximally mutated channels, these “extreme” CARs cover a similar concentration range and the CARs with intermediate degrees of mutations fill the range between the two extreme CARs, thereby developing characteristic shapes.

We first compared the fitted parameters, x_i, with the true parameters, x_i,tr. To increase the robustness of the result, we repeated the global fit 1,000 times, thereby generating the data points for each fit de novo according to Eq. 9. Despite taking the true values as the start vector, a portion of the fits failed, as identified by at least one negative parameter, a physically senseless result. In the following, all remaining non-failing fits with only positive parameters are termed successful_p fits (Fig. 6 A), of which it is a priori not clear if they have run to the correct minimum or to another minimum. The diagram in Fig. 6 A suggests that models with higher numbers of binding sites have more constraints than fits with lower number of binding sites.

From these successful_p fits, the means for each parameter and their standard deviation were computed and plotted together with the respective true parameter in bar graphs (Fig. 7). The mean relative standard deviations σ_i,rel for each parameter of 1,000 fits (ipr_i,1000; Eq. 11) are plotted below each bar graph as a fraction of 1 to illustrate the relative errors. The following conclusions can be drawn: (1) the global fits found with all models for each mean parameter is a minimum near x_i,tr; (2) the errors obtained from these fits were generally fairly small, but those of the last E_x as well as Z₁ and the last Z_x exceeded the others; (3) particularly well determined are the mutational factors of the binding sites, fd_x; and (4) the x_i are apparently best determined with the CHA⁵ model and worst with CHA² model.

To further characterize the global fits, we determined for the successful_p fits the uncertainties of the parameters by two measures, the inaccuracy, the deviation of the fitted parameters from the respective true values, and the imprecision, the deviation of the fitted parameter from the determined mean (see Materials and methods). We first consider the inaccuracy that is only accessible in simulated data as used herein. The inaccuracy for the CHA⁵ model with n_p = 41 parameters is approximately similarly good as that for the CHA⁴ model and CHA³ model with n_p = 23 and n_p = 13 parameters, respectively, and all of these fits seem to be superior compared with the fit with the CHA² model with n_p = 7 parameters (Fig. S3 A). We also compared the imprecision (Fig. S3 B), a measure that would be accessible also for real experimental data. We considered if the somewhat surprising higher goodness of the fits for the channels with the higher parameter numbers can also be obtained from this measure of uncertainty. The diagrams for the four channels confirm the above result, i.e., the channels with more binding sites provide an apparently lower imprecision than the channels with less binding sites.

We next addressed the question of whether the fact that some of the few parameters are apparently less well-determined than the majority of the others can be attributed to special correlations among the parameters. To this end, we computed the matrix of the correlation coefficients from the covariance matrix provided by the fits of all CHA models (see Materials and methods) and plotted them in false colors (Benndorf et al., 2022) to identify conspicuous specific correlations among parameters (Fig. S4). The result was that the parameters with the highest imprecision, the last E_x and Z_x, show a pronounced negative correlation. But there are also a few other relevant negative correlations between the second to last E_x and Z_x as well as some positive correlations for both the CHA⁴ and the CHA⁵ model between high Z_x. Nevertheless, the most relevant result is that non-correlated parameters (colored yellow) dominate for both the CHA⁴ and the CHA⁵ model.

Finally, we computed for the successful_p fits respective overall-measures for all parameters of a fit, of which each is already the mean of n fits as specified above. For the inaccuracy and the imprecision, we calculated the overall means, miac_np,n and mipr_np,n (Eqs. 15 and 12, respectively; Fig. 6, B and C). In the case of the correlation coefficients, we abolished the influence of different signs by taking the absolute values and taking out the ones of the main diagonal, yielding the mean absolute correlation coefficients_n (Fig. 6 D).

These diagrams demonstrate that the chances to determine the parameters in the respective models are best in the models with the higher number of binding sites. To come back to our initial question, the increasing number of constraints by the increased number of CARs in channels with more binding sites overbalances the effect introduced by the increased number of parameters. The most important result, however, is that it seems to be indeed realistic to analyze a complex system with 41 parameters, as given by the CHA⁵ model, if only all CAR combinations with systematically impeded binding sites are available and if selecting the successful_p fits for the analysis.

We first counted the number of successful_p fits (Fig. 8 B). As expected, for all CHA^x models, a weaker mutational effect causes a lower number of successful_p fits. Moreover, at stronger mutational effects (w→10), the CHA⁴ and CHA⁵ models produce better success than the CHA² model. The CHA³ model is intermediate.

When calculating the corresponding means of the imprecision (Eq. 12) over all parameters, it is also lower for the CHA⁴ and CHA⁵ model than the CHA³ model and even much lower than for the CHA² model (Fig. 8 C). Moreover, the results show that with the selected set of parameters, for w a value of 7 or larger is a good idea, resulting for ⁷fd₁ and ⁷fd_max to a 631- and 25-times reduced affinity, respectively.

We next considered how the different degrees of mutational effects affect the E_x and Z_x individually. Regarding E_x, it is evident that there is a pronounced reduction of the imprecision from the CHA² to the CHA⁴ or CHA⁵ model, telling that the lower imprecision (higher precision) of fits with the more complex models, and respectively more parameters, are indeed better determined (Fig. S5). It is also notable that E₁ and E_max are some more uncertain than the respective intermediate E_x. For Z_x, the results are basically similar. Again, there is a pronounced reduction of the imprecision from the CHA² to the CHA⁴ or CHA⁵ model. Also, Z₁ and Z_max are more uncertain than the respective intermediate Z_x (Fig. S6). Together, these results show that the reduced overall imprecision (Fig. 6 C) of fits with the more complex models, and respectively more parameters, is caused by a lower imprecision of both E_x and Z_x.

Our approach so far allowed us to substantiate that the parameter identifiability for the CHA⁵ model was better than that for the CHA² model and at least as good as that for the CHA³ and CHA⁵ model. Besides the fact that this situation does not match the situation of an experimenter, who, of course, does not know any true values, our approach could have been strongly biased by using the true values as a start vector, certainly promoting convergence to the minimum with the true parameters.

To match the requirements of an experimenter, the logic of the approach is now turned around: a set of data points for the respective maximum number of CARs is viewed as given and an approach is developed to analyze the respective data of the CHA² to CHA⁵ model.

To this end, we followed our previous strategy of using SU vectors, varying over a wide range, and vary each of its elements stochastically (Benndorf et al., 2022; Schirmeyer et al., 2022). If sufficiently many successful fits with parameters providing the same minimum can be obtained and the nonsuccessful fits can be reliably eliminated, any bias introduced by a specific start vector can be ruled out.

We varied the SU vector from 10⁻⁶ to 10³, i.e., over nine orders of magnitude, resulting in 19 SU vectors (Fig. 9 A). As a stochastic factor, we used A = 100, i.e., each element of the SU vector was stochastically varied in the range of 10⁻² to 10² of its value. In total, this covers the huge range of 13 orders of magnitude for the stochastically varied SU start vector. The noise factor was again 0.25 and 1,000 fits were performed for n_CAR,max CARs.

We generated one set of stochastically varied data points, as described, and applied the fit routine for the 19 SU vectors à 1,000 fits. We then identified the successful_p fits. Compared with the above approach using the true parameters as start vector, their portion was expectedly very low because the values of the stochastically varied SU start vector must be extremely unfortunate.

Then a second selection algorithm for successful fits out of the successful_p fits was applied.

Using the MSE_x values provided by the fit routine, we identified the minimum MSE_x value, Min(MSE_x) in the 1,000 fits, calculated the difference values ΔMSE_x of each MSE_x to Min(MSE_x) (Eq. 17), and sorted these values starting with the smallest (see Materials and methods). These sorted ΔMSE_x values were plotted for x > 0 on a logarithmic scale (Fig. 9 B). The number of ΔMSE_x values was very low (<<10⁻¹⁰). At a certain x, ΔMSE_x jumped to much higher values >>10⁻¹⁰, indicative of another but much worse fit minimum. These much larger ΔMSE_x values were in no case similar to each other as were the ΔMSE_x values near ΔMSE₁ = 0, suggesting that no other reproducible minimum was present. We then set a threshold at ΔMSE_x = 10⁻¹⁰ to identify all fits in close proximity of ΔMSE₁ = 0. Furthermore, we specified that for a given SU vector at least five fits with ΔMSE_x values below the threshold of 10⁻¹⁰ had to be present among 1,000 fits to accept a true minimum. With these assumptions, we could identify all fits close to the best fit at ΔMSE₁ = 0. These fits are denoted successful_pΔMSE fits in the following. In the plots, the fits with ΔMSE₁ = 0 were omitted. The diagrams show that out of the 1,000 fits, only a small portion but sufficiently many successful_pΔMSE fits were obtained for all four models. The number of successful_pΔMSE fits was either equal or some smaller than that of successful_p fits.

Plotting histograms for the number of successful_pΔMSE fits as a function of the SU vector shows that it was maximal near the SU vector of 10⁻¹ and drops to both larger and smaller SU start vectors to zero (Fig. 10). As a matter of principle, this means that it is possible to identify a range for the SU vectors that yields the majority of successful_pΔMSE fits. Hence, our chosen SU start vectors cover a reasonable range in which a stochastic variation of the elements of the SU start vector by the factor A = 100 identifies a small but reasonable portion of successful_pΔMSE fits. Moreover, as already indicated in Fig. 9 B, the number of successful_pΔMSE fits for all SU vectors decreases for the more complex models.

We next calculated for each model the mean and the standard deviation for all parameters as obtained from the successful_pΔMSE fits with all SU start vectors. These means are again plotted in bar graphs for all parameters together with the true values (Fig. 11). Also the mean relative standard deviations σ_i,rel for each parameter of 1,000 fits (ipr_i,1000; Eq. 11) are plotted below each bar graph as in Fig. 7. All parameters could be determined similarly well as with the analysis above although now there was no bias introduced by any specific start vectors. This result is further confirmed by plotting the inaccuracy and the imprecision (Fig. S7), which both closely match the profiles obtained from the successful_p fits using the true parameters as the start vector.

In summary, our analysis with the stochastically varied SU start vectors used identical data points for the 19,000 fits with all 19 SU vectors, a strategy that could also be run with real experimental data.

As a control, we finally repeated our stochastic approach with different sets of data points for the different SU start vectors, while for the 1,000 fits at each SU start vector, again, identical stochastically varied data points were used. This approach has the advantage that 19 different data point sets were used. If the results are similar to those with identical data points for all SU start vectors, this would confirm that the successful result shown by Figs. 9, 10, and 11; Fig. S7 does not depend on a specific data set. The corresponding Figs. S8, S9, S10, and S11 show a closely similar result to that provided by Figs. 9, 10, 11, and S7, i.e., all parameters are determined, the maximum number of successful_pΔMSE fits was near the SU start vector of 10⁻¹ (Figs. S8 and S9), the parameters match the true parameters (Fig. S10), and the inaccuracy as well as the imprecision (Fig. S11) match those of the corresponding successful_pΔMSE fits shown in Fig. S7 and also those obtained from the successful_p fits shown in Fig. S3. We conclude that our approach is not essentially influenced by a major bias of selecting a specific set of data points.

In this study, we investigated the identifiability of 7–41 free parameters determining the CHA² to CHA⁵ schemes for LGICs containing two to five binding sites, respectively. These parameters define respective larger numbers of equilibrium association constants as specified in Tables S1, S2, S3, and S4. For our comparison of the performance of the CHA² to CHA⁵ models, we included the respective full sets of all possible CARs whose number is given by Eq. 21. The computations showed that the parameters of the most complex CHA⁵ model with 41 parameters are markedly better determined than those of the much simpler CHA² model with only 7 parameters. The CHA³ and CHA⁴ models were in-between. The likely explanation for this apparently counterintuitive result is that the fit constraints rise steeper with the higher complexity than the degrees of freedom introduced by the growing number of parameters.

We like to emphasize again that all computations performed herein were based on the synthetic parameter sets generated by Eq. 23, 24, 25, 27, 29, and 32 in the supplemental text at the end of the PDF. Though the parameter set for LGIC⁴s was related to some extent to the experimental data in heterotetrameric CNG channels (Schirmeyer et al., 2021), it is synthetic, as are the parameter sets for LGIC²s, LGIC⁴s, and LGIC⁵s.

Our strategy is based on three assumptions: (1) a known subunit composition of the channel, (2) a marked reduction of the binding affinity at the binding sites by appropriate mutations, and (3) a global fit of CARs obtained from recording macroscopic currents with coupled models.

To apply our strategy to experimental data of concatameric channels, three assumptions have to be fulfilled: (1) the architecture and function of the concatameric channels must be similar to channels built by non-concatenated subunits; (2) the mutations at the binding sites must markedly reduce but not entirely abolish the binding affinity; and (3) the CARs used for the global fits with coupled models must be of high quality.

Assumption 1 requires that the expected architecture and an approximately unchanged function of the concatenated channels can be demonstrated. Possible side effects of linking the subunits can be easily identified by comparing the CAR of the wt concatamer with the CAR of channels built by individual wt subunits. Regarding the architecture, an incorrect assembly of the mutated subunits in the concatamers has to be excluded which might appear if the subunits of two concatamers build only one functional channel while part of the subunits are “looped out” (Steinbach and Akk, 2011). This phenomenon has been originally reported for voltage-gated K⁺ channels (McCormack et al., 1992) and later also for pentameric ligand-gated nicotinic receptors (Zhou et al., 2003; Groot-Kormelink et al., 2004). A relatively simple functional test would be to take a concatamer containing only wild type subunits and a concatamer containing only mutated subunits, shifting the CAR to much higher concentrations while the slope is unaffected. If coexpressing these concatamers in a 1:1 ratio, channels formed by the two concatamers only would lead to two clearly separated CAR components with unchanged steepness whereas a significant combination of the subunits of the two concatamers within a channel, while looping out others, should lead to a smeared CAR over the whole concentration range. We demonstrated the preservation of the two components for wt and mutated concatamers of homotetrameric CNGA2 channels (Wongsamitkul et al., 2016). Moreover, the consistency of our previous results on concatenated heterotetrameric CNGA4:A2:B1b:A2 channels (Schirmeyer et al., 2021) and homotetrameric HCN2 channels (Sunkara et al., 2018) suggested both a correct architecture and only minor side effects by concatenation. We therefore assume that is generally possible to follow such a strategy. It should also be noted that enforcement of subunits into a desired stoichiometry and position is only one way to define the composition of channels. In case that all subunits differ in their structure and are consistently incorporated in a channel, concatenation is not even required.

Assumption 2 of our strategy includes that the mutations strongly impede the binding sites but preserve functionality, i.e., all mutated binding sites keep their ability to bind a ligand and contribute to activation. We realized this in our simulations by following the previous strategy for tetrameric channels using impeding parameters fd_x for each of the subunits, i.e., the fd_x values were assumed to be unique for each binding site, independent of the activation of the other subunits. Hence, the fd_x factors are helper parameters that are not relevant for interpreting the constants of the wild type channel.

For assumption 3, the global fit with coupled models, we assumed for each binding step individual equilibrium association constants, K_x, but lumped the equilibrium closed–open isomerization constants, E_x, at an equal degree or liganding (Schirmeyer et al., 2021; Benndorf et al., 2022). This significantly reduced the number of parameters. Future simulation studies will show to what extent it is possible to differentiate between the E_x. One possible analysis strategy we already tested is by including data at a second voltage and assuming that each subunit exerts an invariant subunit promotion energy on channel opening (Schirmeyer et al., 2022). Simple additivity of the subunit promotion energies led to individual E_x values for each constellation of liganding. In all fits, we obeyed the principle of microscopic reversibility, which is widely accepted to be physically meaningful though for some ion channels a violation of this principle has been suggested (Schneggenburger and Ascher, 1997; Xu and Meissner, 1998). It should also be noted that we only simulated and analyzed equilibrium constants here which do not tell anything about the dynamic response of the receptors.

We demonstrated that our fit approach can be used for LGICs with two to five binding sites (LGIC² to LGIC⁵) and that an accurate read-out of channel activation is possible.

Among LGIC³s, possible candidate channels for our approach are trimeric purinergic P2X receptors (Kawate et al., 2009; Hattori and Gouaux, 2012) and related ASICs (Jasti et al., 2007) that form between two adjacent subunits three binding sites. For P2X receptors, there are several reports that concatamers are functional (Nicke et al., 2003; Nagaya et al., 2005; Stelmashenko et al., 2012; Gusic et al., 2021) which would allow us to define the position of the mutated binding sites. For ASICs, we are not aware of respective concatamers.

Typical LGIC⁴s with four binding sites are AMPA receptors (Sobolevsky et al., 2009). So far there are no reports about functional concatamers in these receptors, possibly because of the problem that the N-terminus is outside. However, for related NMDA receptors, it has been demonstrated that functional concatamers can be constructed (Schorge and Colquhoun, 2003), which might open the door to applying our strategy. In contrast, tetrameric CNG channels (Li et al., 2017) and hyperpolarization-activated cyclic nucleotide-gated channels (Lee and MacKinnon, 2017) form willingly functional concatamers that can be systematically analyzed (Ulens and Siegelbaum, 2003; Wongsamitkul et al., 2016; Sunkara et al., 2018; Schirmeyer et al., 2021).

Pentameric ionotropic receptors belong mainly to the class of cys-loop receptors. They include the families of nAchRs, GlyRs, and 5HTRs. Many of these receptors are formed by two α-subunits and three other homolog subunits. For example, the classical nAChR of both the electric eel and the electric fish (Changeux et al., 1970; Miledi and Potter, 1971) is composed of two α-subunits, one β-subunit, one γ-subunit, and one δ-subunit (Weill et al., 1974; daCosta and Baenziger, 2013), forming only two bindings sites between the adjacent subunits α and γ as well as α and δ. Therefore, they would be classified herein as LGIC²s. In these channels, the specification of the binding sites need not necessarily be forced by concatenating any subunits because the subunits form on their own channels with defined composition with defined binding sites. Pentameric receptors with five binding sites also exist, e.g., the α7 nAchR (McKay et al., 2007) or the α1 GlyR (Zhu and Gouaux, 2021). In our classification, they would be LGIC⁵s to be described by the CHA⁵ model. It seems to be promising to investigate these or other pentameric receptors with the strategy described herein because the subunits can be linked to form functional concatameric channels (Mazzaferro et al., 2017; Prevost et al., 2021). In summary, there are many relevant candidate channels for our strategy to be applied.

Our approach should also be applicable to other membrane receptors. This includes voltage-gated ion channels (VGICs), for which functional concatamers have been reported e.g., from Shaker subunits (Tytgat and Hess, 1992; Al-Sabi et al., 2010) and Kv2.1 and Kv6.4 subunits (Möller et al., 2020; Pisupati et al., 2020). Here, a promising idea is to affect the voltage-sensor domain of the subunits by appropriate mutations evoking shifts of steady-state activation to more positive voltages. Though measurement of currents in VGICs is easier than that in LGICs, quantitative analysis as described herein would be more complex because each voltage-dependent equilibrium constant requires two parameters, an equilibrium constant at zero voltage and an apparent gating charge. Moreover, our strategy should be easily usable to analyze natural heterotetramers, such as Nav channels (Jiang et al., 2020) and Cav channels (Wu et al., 2015). Furthermore, our analysis could be applied to oligomeric metabotropic receptors as well if a sufficiently valid read-out of receptor activation is available. Advancements in FRET techniques with dimeric glutamate receptors seem to be a promising strategy (Hlavackova et al., 2005; Grushevskyi et al., 2019; Kukaj et al., 2023).

Another intriguing idea is to analyze double mutations at two relevant positions in a channel or receptor, e.g., one mutation at the ligand binding site and a second in the pore gate. If both mutations exert different single effects, a systematic combination of these mutations should allow us to identify the interaction of the subunit with very high precision. However, the number of required combinations might become immediately very high. For example, introducing in a heteroterameric channel a mutation in the binding site and another mutation in the pore gate would result in maximally 16² = 256 combinations. Nevertheless, it is very likely that such a strategy would allow us to specify very subtle interactions of subunits or binding sites in these oligomeric proteins. Notably, this idea is basically an extension of the idea of the double-mutant cycle analysis (Horovitz, 1996).

A further relevant result is that we successfully used a stochastic approach to identify sufficiently many successful_pΔMSE fits out of a vast number of fits by using SU vectors, scanning a range of nine orders of magnitude, and stochastically varying its elements over additional four orders of magnitude. Though the portion of successful_pΔMSE fits was only small compared with the total number of fits, its results regarding accuracy and precision, quantified by inaccuracy and imprecision, respectively, were entirely reliable. Moreover, the scanning of the SU vector over the useful range of nine orders of magnitude allows an investigator to reliably identify the true parameters without prior knowledge of the parameters. For other applications, the range for the SU vectors could be easily extended. Apparently, there is a maximum incidence of the successful_pΔMSE fits in the range of the true parameters, which was in our case ∼10⁻¹ (c.f., Fig. S2). Our stochastic fit approach seems to be useful for any simple or complex fit problem because the investigator does not have to generate specific starting parameters, in particular if there are no reasonable guesses and the number of parameters is large.

We conclude that for LGICs of known composition with two to five binding sites, our strategy to impede, but not abolish, the affinity of defined binding sites by mutations and to globally fit all respective CAR combinations of wild type and impeded binding sites by systems of respective kinetic schemes is a powerful approach to analyze subunit cooperativity of LGICs in great detail. To exclude any bias in the fit arising from specific start vectors, we show that scaled unitary start vectors stochastically varied over an extremely wide range provide reliable fit results. Remarkably, and to some extent counterintuitively, our results show that the identifiability of the parameters in complex models with more binding sites is superior to that in simple models with less binding sites.

The data are available from the corresponding author upon reasonable request.

Joseph A. Mindell served as editor.

This work was supported by the Research Unit 2518 DynIon (project P2 to K. Benndorf) of the Deutsche Forschungsgemeinschaft.

Author contributions: K. Benndorf designed the project, created the simulation and fit software, did the computations, designed the figures, and wrote the manuscript. E. Schulz developed the strategy for using hypercubes in the analysis.