Subunit promotion energies for channel opening in heterotetrameric olfactory CNG channels.

Cyclic nucleotide-gated (CNG) ion channels of olfactory sensory neurons contain three types of homologue subunits, two CNGA2 subunits, one CNGA4 subunit and one CNGB1b subunit. Each subunit carries an intracellular cyclic nucleotide binding domain (CNBD) whose occupation by up to four cyclic nucleotides evokes channel activation. Thereby, the subunits interact in a cooperative fashion. Here we studied 16 concatamers with systematically disabled, but still functional, binding sites and quantified channel activation by systems of intimately coupled state models transferred to 4D hypercubes, thereby exploiting a weak voltage dependence of the channels. We provide the complete landscape of free energies for the complex activation process of heterotetrameric channels, including 32 binding steps, in both the closed and open channel, as well as 16 closed-open isomerizations. The binding steps are specific for the subunits and show pronounced positive cooperativity for the binding of the second and the third ligand. The energetics of the closed-open isomerizations were disassembled to elementary subunit promotion energies for channel opening, [Formula: see text], adding to the free energy of the closed-open isomerization of the empty channel, E0. The [Formula: see text] values are specific for the four subunits and presumably invariant for the specific patterns of liganding. In conclusion, subunit cooperativity is confined to the CNBD whereas the subunit promotion energies for channel opening are independent.

Introduction

Cyclic nucleotide-gated (CNG) ion channels play essential roles in the signal transduction of olfactory sensory neurons and photoreceptors [1-5]. Natural olfactory cyclic nucleotide-gated (CNG) channels [2,4] are heterotetramers composed of three homologue subunits, 2×CNGA2 (A2), CNGA4 (A4), and CNGB1b (B1b) [6,7]. Each subunit contains in its intracellular C-terminus an own cyclic nucleotide-binding domain (CNBD) [2,4,8]. At heterologous expression, only the A2 subunit can form functional channels on its own [9]. In contrast, the A4 and B1b subunits introduce diverse functional effects to heteromeric channels [10-12], as modulating the sensitivity to cyclic nucleotides [7,13-15], reducing the unitary conductance compared to homomeric A2 channels [7], and mediating the action of Ca-Calmodulin [3, 7, 10–12,16–18]. Despite the fact that the A4 and the B1b subunit cannot form functional channels on their own, they also bind cyclic nucleotides and process this signal to the activation process [8,15]. The underlying complex mechanism is still elusive.

The Monod-Wyman-Changeux (MWC) model [19] has been successfully used to quantify cooperative [20] processes in proteins, thereby adopting the strongly simplifying assumptions of fixed stoichiometric factors and a joint ‚allosteric‘ conformational change of identical subunits. Consequently, the MWC model requires at equilibrium only one constant for ligand association (K), one constant for the allosteric conformational change (E0), and one allosteric factor (f), i.e. the same number of constants as in a simple model containing only one binding step and two allosteric steps (Fig 1A). The MWC model, or derivatives of it, have also been repeatedly used to quantify the activation of various homomeric ion channels [21-24] including cyclic nucleotide-gated CNGA1 channels [25] and CNGA2 channels [26,27] as well as cyclic nucleotide-gated HCN2 channels [28]. To a lesser extent, the activation of heteromeric channels [29,30] has been quantified. For heteromeric channels, however, this concept must be even more vague because the different subunits require individual equilibrium constants for both ligand association (K) and the closed-open isomerization (E), resulting for the case of a heterotetrameric protein in 32 different K and 16 different E values (Fig 1B). To determine from experimental data such an amount of parameters seems to be completely unpromising.

Fig 1

State models for the activation in heteromeric proteins.

(A) Scheme with two closed and two open states (C0, C1, O0, O1) for a theoretical minimum ligand-gated channel. There is one ligand binding step from closed state C0 to closed state C1, specified by the equilibrium association constant K. From each closed state an opening step is possible that is specified by the respective closed-open isomerization constant E and E. (B) Related scheme for a heterotetrameric channel containing 16 closed (C0-C15) and 16 open states (O0-O15), requiring for the closed channel 32 binding steps (K-K) and 16 allosteric opening steps (E-E). (C) Cartoon scheme of the concatamer used in the analysis with the sequence N-A4-A2(1)-B1b-A2(2)-C [31]. The colors of the subunits are used throughout the article. (D) Detailed cartoon scheme for the concatameric heterotrameric allosteric (HA) model used in the analysis, specifying the counting of the 32 K and the 16 E, here summarized at equal degree of liganding. (E) Representation of the scheme in d in a 4D hypercube to ease the analysis.

In an attempt to address this question, we recently constructed a full set of 16 concatameric heteromeric CNG channels in the sequence N-A4-A2-B1b-A2-C with defined wild-type and disabled, but still functional, binding sites (Fig 1C). Concentration-activation relationships (CARs) were subjected to a global fit analysis with 16 intimately coupled state models [31]. To simplify the analyses, all K for the two A2 subunits and all E at equal degree of liganding were assumed to be equal. This allowed us to determine 32 K and 4 E with reasonable precision (S1 and S2 Tables).

However, these simplifying assumptions still bias the results because it is a priori neither clear that the two A2 subunits are functionally equal due to their different neighbors, as recently elegantly shown for related heteromeric CNGA1/B1 channels by structural analyses [32], nor that an equal number of ligands bound to different subunits generate an equal E.

To better understand the molecular machinery underlying the activation gating in olfactory CNG channels, our strategy of concatenation, subunit disabling, and global fitting is highly attractive because this might answer the question why these three types of homologue subunits are so different in their function. In particular, it is not clear which parts of the subunits cause these significant functional differences and which parts do interact to generate the cooperative activation and which parts do not. Furthermore, continuative analyses of our unusually consistent data set of 16 concatamers with systematically disabled subunits provides a chance to fathom in a more general way to what degree of detail such constructs can be exploited to functionally analyze heterooligomeric channels with per se silent subunits.

Results

The HACO model

We aimed to determine for heterotetrameric olfactory channels the full set of 32 equilibrium association constants, K, and 16 equilibrium closed-open isomerization constants, E. To this end we analyzed ion currents of 16 N-A4-A2-B1b-A2-C concatamers (Fig 2A and 2B) [31] differing by the number of disabled, but still functional, binding sites. In contrast to our previous report [31], we assumed here that the E can be disassembled to elementary subunit components. An elementary subunit component is specified by the promotion factor for opening, fp (n = 1…4), shifting the closed-open isomerization of the empty channel, specified by the constant E, more to the open state,

Fig 2

Concentration-activation relationships (CARs) and cooperative affinity increase at each subunit.

Current amplitudes were measured as late currents at the end of the pulses at +100 mV and -100 mV. Stepping the voltage from negative to positive voltages generates an activating component in the time course. (A) Currents of heteromeric channels generated by co-expression of wt subunits with a cRNA ratio A2:A4:B1b of 2:1:1. (B) Currents generated by the concatamer N-A4-A2(1)-B1b-A2(2)-C show a similar activating component of the time course as co-expressed channels. (C) Global fit of 32 CARs with the 29HACO model obtained from 16 concatamers at either +100 mV (filled circles) or -100 mV (open circles). The model consists of 32 intimately coupled HA models as described in the text. For normalization of the data points see Materials and Methods. The data points at +100 mV and -100 mV were obtained from the same patches (n = 5 to 18). The curves at +100 mV and -100 mV of the same concatamer are visually related by colored areas.

For empty subunits, the fp are unity, whereas for occupied subunits the corresponding fp have values >1 and are specific. Hence, the 16 E are generated by 16 combinational products of the 5 parameters E, fp…fp.

For the 16 closed states of a single heterotetrameric allosteric (HA) model (Fig 1D) and using a 4D hypercube (Fig 1E), this results in 15 independent virtual equilibrium association constants, Z (S1 Fig), determining 32 K [31]. The closed states of the remaining 15 models with one through four disabled binding sites were treated accordingly (S2 to S5 Figs). When using for each subunit an own disabling factor, fdn (n = 1…4), the number of free parameters for ligand binding add up to 19. Together, with the five disassembling parameters for all E, the global fit of 16 CARs requires 24 parameters. In the following, this model is termed Heteromeric Allosteric Combinational Opening (xHACO) model in which x specifies the number of free parameters.

We first fitted the 24HACO model to the 16 CARs at +100 mV (S3 and S4 Tables). Compared to the simpler 17HA model [31] (c.f. S1 and S2 Tables), the fit precision deteriorated markedly, as indicated by the relative SD (S6 Fig).

To further increase the fit constraints, we included also experimental data from 16 CARs at -100 mV, obtained from the same patches. All CARs at -100 mV were shifted to higher concentrations with respect to the data at +100 mV (Fig 2C). Similar to homomeric CNG1 and CNGA2 channels [33-35], depolarization-induced activation becomes directly evident by the current time courses when stepping from -100 mV to +100 mV (Fig 2A and 2B). In the HACO model, the voltage-dependence was assigned exclusively to the E and not to the ligand binding at the cytosolic CNBD because it is not positioned in the transmembrane electric field. This results for the 32 CARs in the 29HACO model containing the parameters Z-Z, fd-fd, E, fp-.fp, E, and fp. Signs in suffixes denote the voltages of +100 and -100 mV. The result is that the 29HACO model is suitable to describe the gating in its intricate complexity (Fig 2C and Tables 1 and S5). Moreover, the errors of the parameters are markedly small, even smaller than with the 17HA model (S6 Fig). When fitting the data at -100 mV with the 24HACO model alone (S3 and S4 Tables), the errors were similarly large as at +100 mV only (S6 Fig), indicating that the superior precision of the fit with the 29HACO model arises from the increased constraints when fitting the 32 CARs at both voltages.

Table 1

Fit parameters of the global fit with the 29HACO model at +100 mV AND -100 mV.

The dimensions of Z-Z, Z-Z, Z-Z, and Z are μM-1, μM-2, μM-3, and μM-4, respectively. fd1-fd4, E, fp1+-fp4+, E, and fp are dimensionless.

Specification	No.	Parameter	29HACO model +100 mV AND -100 mV
			value	rel. SD (%)
Virtual association constants	1	Z 1	7.31E-02	6.8
	2	Z 2	3.60E-02	10.1
	3	Z 3	2.66E-02	10.9
	4	Z 4	4.50E-02	11.2
	5	Z 5	1.14E-02	7.0
	6	Z 6	6.84E-03	7.9
	7	Z 7	2.37E-02	7.2
	8	Z 8	7.75E-03	8.1
	9	Z 9	4.35E-03	8.8
	10	Z 10	1.90E-02	7.9
	11	Z 11	5.18E-03	7.6
	12	Z 12	9.98E-03	5.2
	13	Z 13	2.41E-03	12.6
	14	Z 14	7.54E-03	6.7
	15	Z 15	1.80E-02	6.0
Factors for disabling binding	16	fd 1	8.81E-05	7.5
	17	fd 2	8.91E-03	8.3
	18	fd 3	9.82E-05	7.5
	19	fd 4	2.43E-03	9.8
	20	E 0+	2.19E-02	3.7
Factors for promoting opening at +100 mV	21	fp 1+	1.21E+01	3.7
	22	fp 2+	6.32E+00	5.3
	23	fp 3+	9.39E+00	3.6
	24	fp 4+	4.18E+00	4.5
	25	E 0-	5.93E-03	3.9
Factors for promoting opening at -100 mV	26	fp 1-	1.74E+01	3.0
	27	fp 2-	4.57E+00	3.6
	28	fp 3-	7.22E+00	2.8
	29	fp 4-	4.01E+00	3.2
Mean s.e.m.	***	***	***	6.7

Cooperativity in ligand binding

First, the cooperativity for ligand binding at different pre-occupation of the other subunits is considered for the closed channel by both the K values and the derived affinity increase at a subunit, f (Eq 15), relating each K to the respective K of the empty channel (S6A–S6D Table and Fig 3A–3D).

Fig 3

Effects of pre-occupied subunits on ligand binding.

(A-D) Equilibrium association constants for the closed channel, K, and relative factors f, specifying how many times the affinity of a subunit is increased by pre-occupation of the other subunits with respect to the same subunit in the empty channel, for A4, A2(2), B1b, and A2(1). ‘E’ indicates the empty channel. The values correspond to S6 Table.

In the empty channel, K-K differ only moderately by a factor of 2.7 in the sequence K>K>K>K (A4>A2(2)>A2(1)>B1b) and the two A2 subunits are closely similar. Pre-occupation of one to three subunits increases the affinity at the binding sites notably, generating positive cooperativity. The degree of this increase, however, is unique for each binding site, also among the two A2 subunits.

Regarding A4 (Fig 3A), single pre-occupation of either of the three other subunits generates an f of 3.5 to 7.2, among which the two A2 neighbors exert a slightly larger effect than do the opposite B1b. When pre-occupying two subunits, a strong affinity increase on A4 is only caused by two occupied A2 subunits, with f of 31.4. The two A2 subunits are clearly different when occupied in combination with B1b (red arrows). Pre-occupation of A2(2) and B1b together generates a fourfold lower binding affinity than pre-occupation of A2(2) alone, suggesting negative cooperativity. Triple pre-occupation does not further enhance f compared to pre-occupation of both A2 subunits. The situation of the opposite B1b subunit resembles that of A4 (Fig 3C).

Regarding A2(1) (Fig 3B), again the two neighbors A4 and B1b exert a stronger effect than the opposite A2(2) at both single pre-occupation and in combination with a pre-occupation of A2(2). In contrast to A4 and B1b, however, here triple pre-occupation causes the strongest affinity increase on A2(1). Regarding A2(2) (Fig 3D), all effects are qualitatively similar to those of A2(1) but differ quantitatively up to a factor of about 4, indicating that different neighbors of the A2 subunits exert a different function.

Together, this means for the cooperativity by binding (1) a pronounced positive and distinct negative cooperativity among the subunits, (2) stronger effects by pre-occupation of neighbor than opposite subunits, and (3) stronger positive cooperative effects of the pre-occupied A2(1) on A4 and B1b than of the pre-occupied A2(2).

Subunit promotion energies for channel opening

To ease the considerations, all 16 E were translated to Gibbs free energies (S7A and S7B Table) according to

R and T are the molar gas constant and the temperature in K. Thus, the four subunit-specific promotion factors fp are translated to free energies, , which specify the energy contributions of the individual subunits to channel opening. These energy contributions are termed in the following ‘subunit promotion energies for channel opening’. They add in a combinational way to . For the two voltages of +100 mV and -100 mV, this results in two sets of 16 ΔΔG built by and as well as and , respectively. Their absolute values follow the order A4>B1b>A2(1)≅A2(2) at both voltages (S8 Table).

Overall, progressive liganding forces the E increasingly to the respective open state, some stronger at +100 mV than at -100 mV. The additivity of the values is illustrated in Fig 4. At equal degree of liganding the different generate notable differences in ΔΔG.

Fig 4

Components of free energies for 16 closed-open isomerizations.

In the fit with the 29HACO model, subunit-specific promotion factors for opening, fp, were used by which the closed-open isomerization constant of the empty channel, E, is shifted to the open state. The fp values are either unity for the empty subunit or adopt a specific value if the subunit is liganded. Disassembling of the 16 E leads to 16 combinational products of the five parameters at both +100 mV (E, fp…fp) and -100 mV (E, fp…fp). From these values specific subunit promotion energies for channel opening were obtained as Gibbs free energies, and (colored arrows) by using Eq (16) which were subtracted from the and values (magenta bars), respectively, illustrating their additive genesis. Overall, successive liganding forces the open-closed isomerizations increasingly to an open state, some more at +100 mV than at -100 mV. Hence, the different ΔΔG values at equal degree of liganding can be reduced to the energetic contributions of the different subunits.

Energetic landscape of activation

Knowing all constants of the HA model for wild-type channels allowed us to consider the complete energetic landscape of the activation process.

Starting with the closed channel, we translated the f values to Gibbs free energies, ΔΔG, using Eq (15), and computed the free energies of C1-C15 with respect to C0, ΔΔGCx, (Eq 17). At P = 0.1 (0.286 μM cGMP) there is a strong endergonic situation (grey bars in Fig 5 and S9 Table), i.e. binding is unlikely to the subunits of all states, despite the differences among them. At P = 0.5 (1.05 μM cGMP), binding to the subunits is still endergonic, though less than at P = 0.1. The tendency is continued at P = 0.9 (3.85 μM cGMP) where many ΔΔGCx values approximate zero. Only the quadruple liganded closed channel has a slightly exergonic level. Together, the energetic differences among the states at an equal degree of liganding superimpose with the dominating effects of the ligand concentration, shifting the equilibrium to less endergonic conditions.

Fig 5

Energetic landscape of the activation process.

Scheme with Gibbs free energies of the 16 closed (C0-C15) and 16 open (O0-O15) states and the HA model, as shown by the cartoons on the top and bottom, respectively. The bar graphs in-between provide the free energies ΔΔGCx and ΔΔGOx at three ligand concentrations, generating Po of 0.1, 0.5, and 0.9, respectively. The grey bars indicate the closed channel, the orange bars the open channel. In both cases the empty channel (C0, O0) was used as reference. The values and their errors are listed in S9 Table. Overall, an increasing ligand concentration shifts the diverse E values towards the open states while the subunit-specific energetic differences at equal degree of liganding superimpose.

Next consider the free energies of the states O1-O15 in the open channel with respect to O0, ΔΔGOx, which is possible because the respective free energies of the 32 binding constants, , are given by the and the two adjacent (Eqs 18 and 19). Herein, only the results at +100 mV are considered because the effects of voltage are second-tier.

At Po = 0.1 there is still an endergonic situation although less endergonic than for the closed channel (Fig 5 and S9 Table). At Po = 0.5, the performance changes significantly: For single-, double-, and triple-liganded channels, ΔΔGOx values change from slightly endergonic to slightly exergonic, again superimposed by differences among the individual states. ΔΔGO15 for the quadruple-liganded channel is strongly exergonic. At Po = 0.9, the exergonic situation is further enhanced. Only the single-liganded channel generates ΔΔGOx around zero.

Reliability of the parameters

To further consolidate the reliability of the exceptionally complex global fit, we used scaled unitary (SU) start vectors, containing identical elements, and varied them stochastically. If sufficiently many successful fits converging to the same absolute minimum can be identified, any influence arising from specific start vectors could be excluded. We used SU start vectors between 10−6 and 100, varied each of the 29 elements stochastically over four orders of magnitude, and repeated each fit 500 times (see Materials and Methods). In total, we obtained 126 successful fits with the criterion that no parameter is negative. Subsequently, we calculated from the of each fit x the related quantity, C(), using Eq (12). Setting a threshold at C() = 10−7 allowed us to easily separate 114 successful fits of high consistency with C() at the numerical resolution limit and to separate them from 12 fits with much larger minima, which were also inconsistent among each other (Fig 6A and 6B). The parameters of the 114 consistent fits match those determined above (Table 1) perfectly and are highly consistent among the fits, as indicated by standard deviations below 10−3 (S10 Table). Hence, the obtained parameters are considered to be well determined.

Fig 6

Identification of successful fits with SU start vectors.

As start vectors scaled unitary (SU) start vectors, containing identical elements, were used and varied stochastically. The SU start vectors had values between 10−6 and 100. The 29 elements of the SU start vectors were varied with the stochastic factor A = 100R. R is a linearly distributed random number between -1 and +1. With each SU start vector 500 fits were performed, resulting in total in 3,500 fits (see Materials and Methods). Non-successful fits containing negative parameters were immediately eliminated, leaving 126 successful fits for all SU start vectors. (A) Histogram of the 126 successful fits as function of the value of the SU start vector. With SU start vectors of 10−3, most successful and consistent fits were obtained. The two categories were calculated as described in B. (B) Procedure to identify highly consistent fits at the actual minimum of χ2. From of each successful fit x we calculated the related quantity C() by Eq (12) (see Materials and Methods). 114 successful fits with high consistency at the resolution limit could be easily separated from 12 much less consistent outliers by setting a threshold to 10−7.

Discussion

We describe for the activation of heterotetrameric olfactory channels specific subunit promotion energies for channel opening, , for each subunit in the context of the 32 highly cooperative binding events. The combinations of these four specify together with the energy of the empty channel all 16 closed-open isomerizations of the 29HACO model. For differentiating the 16 closed-open isomerizations at each combination of liganding we exploited a weak voltage dependence of the CARs.

Methodologically, our analysis was based on a combination of 16 concatameric channels with the sequence N-A4-A2(1)-B1b-A2(2)-C and systematically disabled, but still functional, binding sites and on a global mathematical fit with 32 intimately coupled state models.

The trustability of our global fit was demonstrated by a stochastic approach. Additional support comes from the correlation matrix of the 29 fit parameters, yielding for the vast majority of elements values near zero (S7 Fig), suggesting that the constraints in our fit are strong. Moreover, although the number of 29 free parameters seems to be enormously high, it is indeed low when comparing it with 124 parameters required when fitting Hill or double Hill functions to the 32 CARs separately (S11 Table).

As in any case of interpretation of data by a model, also our interpretations depend on the specificities of the used 29HACO model. We like to state again that the 29HACO model seems to us highly plausible for a heterotetrameric channel with four binding sites despite its considerable complexity. The key assumptions of this model are:

We assigned the moderate effects of voltage exclusively to the closed-open isomerizations (E) of the channel at the different combinations of liganding. This assumption is highly plausible because the closed-open isomerizations proceed in the channel core subjected to the transmembrane electric field but not in the binding reactions (K) proceeding in the CNBDs located in the cytosol, and thus outside the transmembrane electric field.

For our HA models, underlying the complex 29HACO model, we assumed that the four subunits in the channel exert only a binary closed-open action in a highly cooperative way, leading to a single opening step in our models. This assumption was supported by the observation, that the ligand concentration controls solely the open probability, Po, but not the single-channel conductance (S8–S10 Figs). We therefore did not consider models of the Koshland-Nemethy-Filmer (KNF) type [36,37].

Our approach assumed that of each occupied subunit is independent of the occupancy of the other subunits. Thus, we reduced the complex interaction of the subunits in the closed-open isomerizations into four independent functional modules. This assumption is a priori not without alternative. The success of our fit, however, suggested to us that this assumption is reasonable.

Modular functionality within the protein was also assumed for the disabling factors fd, i.e. they were taken as independent of the occupation of the other subunits [31]. Again, the success of our fit suggested to us that this assumption is reasonable.

Under these model assumptions, the phenomenon of subunit cooperativity would be confined to the level of the CNBDs. This interpretation gets nice support from recent computational studies on the allosteric signaling within the CNBD and the C-linker of structurally related HCN2 channels [38]. The authors applied MD simulations and a rigidity-theory-based approach and identified complex interactions to mediate the cAMP effect. They identified in the CNBD and the C-linker two intersubunit pathways and one intrasubunit pathway. Our analysis leaves open which structures further process the information from the CNBD and the C-linker to the pore but suggests that that this process runs independent in each subunit, i.e. there is no major additional cooperative interaction of the subunits at the level of the channel core.

The overall view is the following: Ligand binding to the different subunits of a heterotetrameric olfactory channel appears with different affinity, independent of the state of the other subunits but highly sensitive to mutagenesis (Fig 7). In contrast, subsequent conformational changes in the tetrameric CNBD modulate the affinity of the CNBDs with high specificity in a cooperative fashion whereas the specific of a subunit, controlling the channel pore, is independent again.

Fig 7

Schematic summarizing the nature of the processes following ligand binding.

Pre-occupied subunits modulate the binding affinity of the other binding sites, presumably at the level of the CNBD, generating pronounced cooperativity. In contrast, both ligand binding itself and the subunit promotion energies for channel opening, , determined herein are independent processes.

Materials and methods

Ethics statement

The procedures had the approval of the authorized animal ethical committee of the Friedrich Schiller University Jena. The methods were carried out in accordance with the approved guidelines.

Molecular biology and heterologous expression

The heterotetrameric concatamers assembling to olfactory CNG channels were obtained by joining the coding sequences of two CNGA2 (accession No. AF126808) subunits, one CNGA4 (accession No. U12623) and one CNGB1b (accession No. AF068572) subunit from the rat. The methods are essentially the same as reported previously [31]. The point mutations R538E (CNGA2), R430E (CNGA4) and R657E (CNGB1b) were introduced via the overlapping PCR technique, yielding all 16 combinations in the A4-A2-B1b-A4 concatamer.

Oocytes were harvested surgically under anesthesia (0.3% 3-aminobenzoic acid ethyl ester) from adult females of Xenopus laevis[35]. Oocytes were digested with collagenase A (3 mg/ml, Roche, Grenzach-Wyhlen, Germany) for 105 min in Ca2+-free Barth´s solution containing (in mM) 82.5 NaCl, 2 KCl, 1 MgCl2, 5 HEPES, pH 7.5. Oocytes of stage IV and V were manually dissected and injected with ~50 ng of cRNA encoding the respective channels. The oocytes were cultured at 18°C for 1 to 3 days in Barth’s solution containing (in mM) 84 NaCl, 1 KCl, 2.4 NaHCO3, 0.82 MgSO4, 0.41 CaCl2, 0.33 Ca(NO3)2, 7.5 TRIS, Cefuroxim, Penicillin/Streptomycin, pH 7.4.

Electrophysiology

For obtaining concentration-activation relationships at equilibrium, ensemble (macroscopic) currents, generated by hundreds to several thousands of channels, were recorded from inside-out patches with the patch-clamp technique [39]. The amplitude of the late current at +100 mV was evaluated. The patch pipettes were prepared from quartz tubing (VITROCOM, New Jersey, USA) on a P-2000 puller (Sutter Instrument, Novato, USA). The outer and inner diameter of the tubing was 1.0 and 0.7 mm. The pipette resistance was 0.5–1.7 MΩ. Both bath and pipette solution contained (in mM): 150 KCl, 1 EGTA, 10 Hepes (pH 7.4 with KOH). Recording was performed with either an Axopatch 200B amplifier (Axon Instruments, Foster City, CA), controlled by the ISO2 hard- and software (MFK, Niedernhausen, Germany), or with an EPC10 amplifier and the Patchmaster software (HEKA GmbH, Lambrecht, Germany). The sampling rate was either 2 or 5 kHz and the internal filter of the amplifier was set to either 1 or 2 kHz. The data at +100 mV were the same as published previously [31].

For single-channel recording, the patch pipettes were also fabricated from quartz tubing. The outer and inner diameter was 1.0 and 0.5 mm, respectively (VITROCOM, New Jersey, USA). The pipette resistance was 5.0–12.0 MΩ. The pipette solution contained (in mM): 150 KCl, 1 EGTA, 5 Hepes (pH 7.4 with KOH). The recording voltage was +100 mV. The data were recorded at +100 mV or 100 mV in the inside-out patch-configuration at different cGMP concentrations. Recordings from wt heteromers were sampled at 20 kHz and filtered to 5 kHz on-line. All other recordings were sampled at 40 kHz and filtered to 10 kHz on-line. For display the data were off-line filtered to 1 kHz by a Gaussian filter. Amplitude histograms were built from 10-second intervals and fitted with the sum of two normalized Gaussian functions from which the open probability, Po, and the amplitude of the unitary current, i, were obtained.

Data analysis

The global fit

In our recent analysis of the Heteromeric Allosteric (HA) model [31], we used two equal A2 subunits and an equal opening constant E at equal degree of liganding. Herein these two limitations were overcome by extending the HA model to the more general Heteromeric Allosteric Combinational Opening (HACO) model.

Similar to the HA model, the 16 closed states Cijkl are arranged as corners of a 4D hypercube (Fig 1E) [31]. Each state is assigned a four-dimensional vector with the binary indices i, j, k, and l that are either 0 or 1.

The 32 edges of the hypercube correspond to the 32 possible equilibrium association constants, K (i = 1,…,32), between neighboring closed states.

For computations we adopted the assumption of microscopic reversibility [40], resulting for a 16-state model in 15 independent K while the others are given by cycles. With respect to state C0000, this results immediately in 15 virtual equilibrium association constants Z (i = 1,…,15) (violet lines in S1 Fig), defining the 32 K by respective ratios as provided by S2 Table. Microscopic reversibility holds automatically.

Now the variable p0000 specifies the equilibrium occupation probability of the empty closed state C0000. Then the occupation probability pijkl of each other closed state Cijkl can be easily determined from the Zijkl according to

The ligand concentration L appears as the power of the sum of the indices, reaching a maximum value of 4. The factors of disabling by mutation fd (u = 1,…, 4) of the four subunits have the exponents a, b, c, d (equal to 0 or 1). An exponent is 1 if the subunit u is mutated and has bound a ligand. Otherwise, the exponent is zero.

In contrast to the HA model in [31], the influence of each subunit on the opening process is treated separately. The closed-open isomerization constant of the empty channel is E. When a ligand binds to subunit u, an additional factor fp (u = 1 … 4) promotes opening. Then the total opening constant Eijkl between the states Cijkl and Oijkl is given by

This results in the probability pijkl of the open state Oijkl to be occupied of

The open probability P (L) of the whole channel is then given by

Because the occupation probabilities of all states add up to one, this allows to eliminate p0000.

For an effective notation, we define the five subtotals S and T (α = 1,…,4)

Then the open probability P (L) of the whole channel is given by

The promotion factors fp of each binding site are model parameters and their fit result delivers an insight into the contribution of each binding site to the opening process.

The squared differences between the calculated and the measured open probabilities P (L) and P (L) are summed up to χ 2.

The summation covers all n = 16 concatamers, at the respective n concentrations. The weighting factors are the reciprocal of the empirical variances of the mean, which have been estimated from measurements in 6 to 18 patches. Minimization of χ2 was performed with the Levenberg-Marquardt algorithm [41], providing the fit parameters and their covariance matrix cov.

Test for fit validity by stochastic variation of the start parameters

To evaluate the validity of our global fit approach with the 29HACO model minimizing χ2, we designed a strategy to avoid any biasing influence of specific start vectors. We used scaled unitary (SU) start vectors, containing identical elements, and varied them stochastically. If sufficiently many successful fits can be obtained and the non-successful fits can be identified, a major source of bias would be removed. We varied SU start vectors between 10−6 and 100 in decade steps and varied each of the 29 parameters stochastically by a linear stochastic factor A = 100R where R is a random number between -1 and +1. This results in a variation of all parameters over a range of 104. For each SU start vector, the fit was repeated 500 times, resulting in total in 3,500 fits. The maximum number of iterations was set to 200.

In a first step, non-successful fits were identified in which at least one parameter ran into a negative value. This criterion was chosen because a negative equilibrium constant is physically nonsense. The number of the remaining 126 successful fits, containing only positive parameters, depended on the SU start vector reaching a maximum at 10−3 (Fig 6A).

In a second step, we considered the distribution of the minima among the successful fits.

To this end, we used the value of each successful fit x and calculated a related quantity of according to

is zero at the minimum and increases at already very subtle changes. The major advantage of using is that it allows to bring very heterogeneous values onto a single logarithmic scale apart from the minimum value = 0 itself. It turned out that of 114 out of the 126 successful fits was either 0 or 2.2073×10−8, the resolution limit of our calculations. If counting all fits with = 0 and = 2.2073×10−8 together ( ≤ 2.2073×10−8), all 126 could be plotted (Fig 6B). The diagram illustrates the enormous consistence of for 114 fits at the resolution limit as well as 12 outliers which were much larger and highly inconsistent. This allowed us to set a threshold at 10−7 and to distinguish the most frequent and consistent fits at the least minimum easily from rare and inconsistent fits at much larger minima (Fig 6B). The consistent fits at the consistent least minimum are considered to provide the result of the global fit. Notably, these minima were determined without any prior knowledge of the parameters. Although these results make the identified parameters very likely to represent the best fit, these results do not finally exclude a better minimum in the 29-dimensional parameter space.

Fits with Hill functions

For comparison, concentration-activation relationships were fitted with IgorPRO 7 (Lake Oswego, USA) by where I denotes the actual current amplitude and Imax the maximum current amplitude at saturating cyclic nucleotide CN specified for each patch. EC50 is the cGMP concentration evoking half maximum current and n the Hill coefficient.

Part of the concentration-activation relationships required the sum of a high (H) and a low affinity (L) component

The notation corresponds to that in Eq (13). A is the fraction of the high affinity component.

Computation of free energies

Free energies specifying ligand binding to a subunit with respect to the binding to this subunit in an otherwise empty channel were determined from the promotion factors f by R and T are the molar gas constant and the temperature in K.

The free energies for the closed open isomerizations E and E, as well as for fp1+-fp4+ and fp, the values at -100 mV and +100 mV, were calculated accordingly by where h is the respective parameter.

For a given ligand concentration L, the free energies of all closed states C1….C15 with respect to C0 were determined by summing up the ΔΔG values, along each chosen pathway, and adding the respective free energies provided by the ligand concentration L, ΔΔG = −RTln (L), according to

The free energies of the 32 binding constants, , for the open channel were calculated from the respective and the two adjacent according to where and are the free energies before and after the ligand binding, respectively. The free energy of the open states was determined in an analogue fashion to that for the closed states

Error propagation

The fit program provides a set of parameters together with their standard deviations and relative errors. For interpretation of the fit results we calculated from these parameters derived quantities Y of interest. The errors of these quantities were determined according to the rules of error propagation as follows:

A quantity Y (column vector Y, …, Y) depending on the parameters p, p, …, p would be given by

To obtain the statistical errors of Y, the covariance matrix of the Y, cov, is calculated from the covariance matrix of the parameters, cov, with the help of the Jacobian matrix J according to using the rules of matrix multiplication.

The elements J of the Jacobian matrix are the partial derivatives of Y with respect to the parameters p.

The covariance matrix cov is quadratic of type (n x n). The Jacobian matrix J is of type (m x n) and the covariance matrix cov is quadratic of type (m x m).

The standard deviations σ of the Y are obtained from the main diagonal elements of the covariance matrix cov according to

In this way we calculated the errors of the association constants K…K for wild-type and mutated subunits as well as the errors of the opening constants E…E for the HACO model.

The error of the change of the Gibbs free energy, ΔΔΔG, results from the error of the factor h (Eq 16) according to following the rules of the differentiation of a logarithmic function. Δh was calculated from the error of the fit parameters using Eq (22) and an appropriate Jacobian matrix (21).

Determination of Po in CARs

The equilibrium open probability, Po, used in the fits was determined from the amplitude of late currents at the end of the pulses to +100 or -100 mV (Figs 2A and 2B). At +100 mV and infinite cGMP, we assumed that all constructs generate a maximum open probability of 0.99, which is based on single-channel recordings of concatenated wild-type channels [31]. Example single-channel recordings of wt concatamers, triple-mutated concatamers, quadruple-mutated concatamers as well as non-concatenated channels are shown in S8 to S10 Figs, respectively. The currents of triple and quadruple mutated concatamers, not reaching saturating currents at +100 mV, were re-scaled by Po determined in corresponding single-channel measurements for +100 mV as described [31]. At -100 mV, Po was calculated from 5 to 15 patches per construct by relating for each concatamer the late current at -100 mV to that at +100 mV, scaled by a correction factor for the different single-channel currents at both voltages and assuming also at -100 mV Po = 0.99 at infinite cGMP. The latter was based on single-channel recordings at -100 mV where at 100 μM cGMP, the non-mutated concatamer 1234 (S8 Fig) generated consistently larger Po than 0.95. The correction factor was determined using the constructs 1234m, 12m34 and 12m34m, which clearly reached saturation. Here, the ration (-current(-100mV)/current(+100mV)) was on average(n = 26) 1.10+/-0.01. This value was validated by 2 independent approaches, delivering comparable values: Un-normalized currents from an individual patch of 4m2m12m were fitted with a Hill function for -100 and +100 mV. The respective maxima and EC50 were free, the Hill parameter was linked between the fits. The ratio of the fitted amplitudes were 1.09+/-0.06. The ratio of single channel amplitudes in a prolonged measurement for 4m2m1m2 was 1.072+/-0.005. In three concatamers (1234, 12m34, 12m34m), Po at -100 mV determined in this way slightly exceeded unity. These values were arbitrarily set to 0.99.

The errors for the parameters and constants are given as relative standard deviation (SD) in percent (rel. SD (%)). For Gibbs free energies the errors are given as SD.

Virtual equilibrium association constants in the 4D hypercube.

(DOCX)

Click here for additional data file.

HA models for four concatamers containing one disabled binding site.

(DOCX)

Click here for additional data file.

HA models for six concatamers containing two disabled binding sites.

(DOCX)

Click here for additional data file.

HA models for four concatamers containing three disabled binding sites.

(DOCX)

Click here for additional data file.

HA model for the concatamer containing four disabled binding sites.

(DOCX)

Click here for additional data file.

Comparison of errors obtained by fits with different models.

(DOCX)

Click here for additional data file.

Color-coded matrix of correlation coefficients for the fit with the 29HACO model.

(DOCX)

Click here for additional data file.

Single-channel activity in wt concatamers N-A4-A2-B1b-A2-C (1234).

(DOCX)

Click here for additional data file.

Single-channel activity in the concatamers 1m2m3m4 and 1m2m3m4m.

(DOCX)

Click here for additional data file.

Single-channel activity in non-concatenated wt channels.

(DOCX)

Click here for additional data file.

Fit parameters of the global fit with the 17HA model.

(DOCX)

Click here for additional data file.

Equilibrium constants derived from the global fit with the 17HA model.

(DOCX)

Click here for additional data file.

Fit parameters of the global fit with the 24HACO model at +100 mV OR -100 mV.

(DOCX)

Click here for additional data file.

Equilibrium constants derived from the global fit with the 24HACO model at +100 mV OR -100 mV.

(DOCX)

Click here for additional data file.

Equilibrium constants derived from the global fit with the 29HACO model at +100 mV AND -100 mV.

(DOCX)

Click here for additional data file.

Influence of pre-occupation on ligand binding for global fit with the 29HACO model at +100 mV AND -100 mV.

(DOCX)

Click here for additional data file.

Closed-open isomerization constants given by the global fit with the 29HACO model at +100 mV AND -100 mV.

(DOCX)

Click here for additional data file.

Subunit promotion energies for channel opening given by the global fit with the 29HACO model at +100 mV AND -100 mV.

(DOCX)

Click here for additional data file.

Free energies of closed and open states obtained by globally fitting the 29HACO model +100 mV AND -100 mV.

(DOCX)

Click here for additional data file.

Parameters determined with stochastically varied SU start vectors.

(DOCX)

Click here for additional data file.

Hill parameters for the 32 CARs at +100 mV and -100 mV.

(DOCX)

Click here for additional data file.

31 May 2022

Dear Dr. Hospital,

Thank you very much for submitting your manuscript "Relating energetic opening momenta to ligand binding in individual subunits of heteromeric olfactory CNG channels" for consideration at PLOS Computational Biology.

As with all papers reviewed by the journal, your manuscript was reviewed by members of the editorial board and by several independent reviewers. In light of the reviews (below this email), we would like to invite the resubmission of a significantly-revised version that takes into account the reviewers' comments.

We cannot make any decision about publication until we have seen the revised manuscript and your response to the reviewers' comments. Your revised manuscript is also likely to be sent to reviewers for further evaluation.

When you are ready to resubmit, please upload the following:

[1] A letter containing a detailed list of your responses to the review comments and a description of the changes you have made in the manuscript. Please note while forming your response, if your article is accepted, you may have the opportunity to make the peer review history publicly available. The record will include editor decision letters (with reviews) and your responses to reviewer comments. If eligible, we will contact you to opt in or out.

[2] Two versions of the revised manuscript: one with either highlights or tracked changes denoting where the text has been changed; the other a clean version (uploaded as the manuscript file).

Important additional instructions are given below your reviewer comments.

Please prepare and submit your revised manuscript within 60 days. If you anticipate any delay, please let us know the expected resubmission date by replying to this email. Please note that revised manuscripts received after the 60-day due date may require evaluation and peer review similar to newly submitted manuscripts.

Thank you again for your submission. We hope that our editorial process has been constructive so far, and we welcome your feedback at any time. Please don't hesitate to contact us if you have any questions or comments.

Sincerely,

Alexander MacKerell

Associate Editor

PLOS Computational Biology

Jason Haugh

Deputy Editor

PLOS Computational Biology

***********************

Reviewer's Responses to Questions

Comments to the Authors:

Please note here if the review is uploaded as an attachment.

Reviewer #1: The MS “Relating energetic opening momenta to ligand binding in individual subunits of heteromeric olfactory CNG channels” by Jana Schirmeyer et al deals with an important topics in the biophysics of olfactory CNG channels. These channels are heteromeric and are formed by two alpha CNGA2 subunits (A2) one CNGA4 (A4) and one CNGB1b. In the present MS, the authors describe the results they have obtained by using the concatemers N-A4-A2-B1b-A2-C concatemers where they have selectively replaced an Arginine in the CN binding site into a Glutamate so that cAMP cannot bind any more.

I like the scope and the experiments here described, but the MS in its present form is rather obscure and even though I am familiar with the topic I was not able to understand it properly. Therefore I suggest to have a mayor revision and here below are my suggestions:

1 - In the introduction it is useful and almost necessary to explain the motivations for building such model and which new insights and/or a better understanding are expected.

2 - there must be a paragraph clearly stating the assumptions on which the model/analysis is based

3 - In the discussion it is necessary to evaluate the conclusions gained by such modelling/computation

The above three suggestions are - in my view - compulsory and now I list some questions for the authors. Indeed they have the data to answer either in a final and definitive way or at least to provide a substantial better understanding for several important issues :

from the data presented in the MS it seems that some channel opening are observed when all or most of the CN binding sites are knocked off by the mutation R/E. Do the authors have some single channel recordings from these concatened and mutated channels? If so I will be delighted to see these recordings.

in long single channel recordings from the CNGA1 monomeric channel it is possible to observe occasionally long closures. Have the authors observed anything similar in their concatenated channels? If so is this observation compatible with the assumptions of their modelling?

In the MS there is a number of words, which are awkward with a not obvious meaning such:

The authors refer to Energetics of subunit opening momenta. The word “momenta” in science is the product of the mass times the velocity and I presume that the authors do not intend this.

The words endergonic and exergonic are often used and their meaning must be defined or – better – more conventional words must be used.

In the legend of Fig.1 it is stated:

“ Scheme with two closed and two open states (C0, C1, O0, O1) for a theoretical minimum ligand-gated receptor containing for the closed channel one binding step (K1) and two allosteric opening steps “

All this is hard to understand and very confusing and leaves the impression of a rather sloppy writing/editing, as if the MS is a preliminary draft and not a finished and well-polished MS

Reviewer #2: Understanding the complexities of heteromeric ion channel gating has been difficult and many electrophysiological studies to measure heteromeric CNG channel currents have been reported. Here, Schirmeyer and authors refine recently generated models to understand the energetic contributions of individual subunits to gating by ligand binding and voltage. This is an interesting question and a unique way to look at allostery and cooperatively in these channels.

The same concatemers from Schirmeyer et al. 2021 and the resulting patch clamp data recorded at +100 mV were used for both Benndorf et al. 2022 JGP and this current manuscript. For the computational model in this manuscript, the authors add patch clamp data for the same constructs recorded at -100 mV. They fit 32 equilibrium association constants and 16 equilibrium isomerization constants in contrast to the 32 equilibrium association constants and 4 equilibrium isomerization constants (Benndorf et al. 2022 JGP). Most of the results expound upon the additional and new fit constants or derivations thereof (free energies). Some of these fit constants are similar while others are more variable between the previous HA (+100 mV) and HACO models (+100 mV and -100 mV).

Major point

Improvement in fitting are less than 5% relative standard deviations (Figure S6) or approximately equivalent to the error of estimated constants (Table 1), which I believe reduces the enthusiasm for the actual enhancement of the HACO model (+100 mV and -100 mV) over the HA model (+100 mV). Could the authors comment on this?

Minor points

Why are equivalent KCl concentrations in the bath and pipette for the patch clamp experiments?

Show single channel activity of natural CNGs at -100 mV to complement Figure S8. It would also be nice to see representative comparisons for mutants between natural and concatemer channels at the different Po's described at both -100 mV and +100 mV.

While these concatemers were used in previous publications to generate models, are the authors concerned with the differences in kinetics between the natural CNG channel and the concatemer and what this could mean for their conclusions (Figure 2A and 2B)?

Can the authors speculate why the constants are so different between +100 mV and -100 mV?

**********

Have the authors made all data and (if applicable) computational code underlying the findings in their manuscript fully available?

The PLOS Data policy requires authors to make all data and code underlying the findings described in their manuscript fully available without restriction, with rare exception (please refer to the Data Availability Statement in the manuscript PDF file). The data and code should be provided as part of the manuscript or its supporting information, or deposited to a public repository. For example, in addition to summary statistics, the data points behind means, medians and variance measures should be available. If there are restrictions on publicly sharing data or code —e.g. participant privacy or use of data from a third party—those must be specified.

Reviewer #1: No: I have not seen neither the code nor the original data, but they could be somewhere and I did not notice

Reviewer #2: Yes

**********

PLOS authors have the option to publish the peer review history of their article (what does this mean?). If published, this will include your full peer review and any attached files.

If you choose “no”, your identity will remain anonymous but your review may still be made public.

Do you want your identity to be public for this peer review? For information about this choice, including consent withdrawal, please see our Privacy Policy.

Reviewer #1: Yes: Vincent Torre

Reviewer #2: No

Figure Files:

While revising your submission, please upload your figure files to the Preflight Analysis and Conversion Engine (PACE) digital diagnostic tool, . PACE helps ensure that figures meet PLOS requirements. To use PACE, you must first register as a user. Then, login and navigate to the UPLOAD tab, where you will find detailed instructions on how to use the tool. If you encounter any issues or have any questions when using PACE, please email us at .

Data Requirements:

Please note that, as a condition of publication, PLOS' data policy requires that you make available all data used to draw the conclusions outlined in your manuscript. Data must be deposited in an appropriate repository, included within the body of the manuscript, or uploaded as supporting information. This includes all numerical values that were used to generate graphs, histograms etc.. For an example in PLOS Biology see here: http://www.plosbiology.org/article/info%3Adoi%2F10.1371%2Fjournal.pbio.1001908#s5.

Reproducibility:

To enhance the reproducibility of your results, we recommend that you deposit your laboratory protocols in protocols.io, where a protocol can be assigned its own identifier (DOI) such that it can be cited independently in the future. Additionally, PLOS ONE offers an option to publish peer-reviewed clinical study protocols. Read more information on sharing protocols at https://plos.org/protocols?utm_medium=editorial-email&utm_source=authorletters&utm_campaign=protocols

28 Jun 2022

Submitted filename: Answers to the Reviewers Rev1_KB3.docx

Click here for additional data file.

11 Jul 2022

Dear Dr. Benndorf,

We are pleased to inform you that your manuscript 'Subunit promotion energies for channel opening in heterotetrameric olfactory CNG channels' has been provisionally accepted for publication in PLOS Computational Biology.

Before your manuscript can be formally accepted you will need to complete some formatting changes, which you will receive in a follow up email. A member of our team will be in touch with a set of requests.

Please note that your manuscript will not be scheduled for publication until you have made the required changes, so a swift response is appreciated.

IMPORTANT: The editorial review process is now complete. PLOS will only permit corrections to spelling, formatting or significant scientific errors from this point onwards. Requests for major changes, or any which affect the scientific understanding of your work, will cause delays to the publication date of your manuscript.

Should you, your institution's press office or the journal office choose to press release your paper, you will automatically be opted out of early publication. We ask that you notify us now if you or your institution is planning to press release the article. All press must be co-ordinated with PLOS.

Thank you again for supporting Open Access publishing; we are looking forward to publishing your work in PLOS Computational Biology.

Best regards,

Alexander MacKerell

Associate Editor

PLOS Computational Biology

Jason Haugh

Deputy Editor

PLOS Computational Biology

***********************************************************

Reviewer's Responses to Questions

Comments to the Authors:

Please note here if the review is uploaded as an attachment.

Reviewer #1: the ms is acceptable in its present version

Reviewer #2: In my opinion, I believe the authors should consider that a 10% standard deviation versus a roughly 12% standard deviation between the different models is probably not that significant. While it is interesting to think about additional parameters and what this can mean in terms of gating, I still do not see a major mathematical improvement from the simpler model with fewer terms. It's perhaps even more puzzling that the improvements in the HACO model (+100 mV and -100 mV) seemingly invoke voltage criteria for channels that are thought to have very weak voltage dependence, as the authors noted. Reviewer comments were addressed adequately.

**********

Have the authors made all data and (if applicable) computational code underlying the findings in their manuscript fully available?

The PLOS Data policy requires authors to make all data and code underlying the findings described in their manuscript fully available without restriction, with rare exception (please refer to the Data Availability Statement in the manuscript PDF file). The data and code should be provided as part of the manuscript or its supporting information, or deposited to a public repository. For example, in addition to summary statistics, the data points behind means, medians and variance measures should be available. If there are restrictions on publicly sharing data or code —e.g. participant privacy or use of data from a third party—those must be specified.

Reviewer #1: Yes

Reviewer #2: Yes

**********

PLOS authors have the option to publish the peer review history of their article (what does this mean?). If published, this will include your full peer review and any attached files.

If you choose “no”, your identity will remain anonymous but your review may still be made public.

Do you want your identity to be public for this peer review? For information about this choice, including consent withdrawal, please see our Privacy Policy.

Reviewer #1: Yes: Vincent Torre

Reviewer #2: No

17 Aug 2022

PCOMPBIOL-D-22-00564R1

Subunit promotion energies for channel opening in heterotetrameric olfactory CNG channels

Dear Dr Benndorf,

I am pleased to inform you that your manuscript has been formally accepted for publication in PLOS Computational Biology. Your manuscript is now with our production department and you will be notified of the publication date in due course.

The corresponding author will soon be receiving a typeset proof for review, to ensure errors have not been introduced during production. Please review the PDF proof of your manuscript carefully, as this is the last chance to correct any errors. Please note that major changes, or those which affect the scientific understanding of the work, will likely cause delays to the publication date of your manuscript.

Soon after your final files are uploaded, unless you have opted out, the early version of your manuscript will be published online. The date of the early version will be your article's publication date. The final article will be published to the same URL, and all versions of the paper will be accessible to readers.

Thank you again for supporting PLOS Computational Biology and open-access publishing. We are looking forward to publishing your work!

With kind regards,

Zsofia Freund

PLOS Computational Biology | Carlyle House, Carlyle Road, Cambridge CB4 3DN | United Kingdom ploscompbiol@plos.org | Phone +44 (0) 1223-442824 | ploscompbiol.org | @PLOSCompBiol