Central Assumptions of the New Framework

Considering pairs of amino acid positions in a protein family, assume that, for each pair, the two positions either strongly and permanently coevolve with each other or evolve completely independently. Let denote the set of coevolving pairs. Let represent the set of (structural) contact pairs, specifically side chains contacts. Following pioneering studies [13], [14], [16] an intimate relationship has been conjectured between coevolution and side chain contact. The relationship can be stated in terms of the probabilities and that, for some protein family, a random draw from all pairs or from contact pairs, respectively, gives a coevolving pair:(1)This says that the contact pairs tend to be the coevolving pairs. Let be the set of coevolving pairs predicted by some coevolution detector from sequence data . If the detector is useful then conditioning on has similar effect to conditioning on :(2)Supporting the preceding two assertions it has been shown repeatedly [11]–[14], [16], [20], [22], [23], [29]–[32], [35]–[42] that most detectors can predict contact pairs better than random choice, and so(3)Instead of predicting contact pairs to aid de novo prediction of structure, several studies [11], [18]–[25] aimed to detect coevolving pairs given the set of contact pairs assuming that

(4)The new framework was designed towards that aim and takes all above assumptions and findings as a starting point. As Figure 2 shows, depends on a set of parameters , which specifies the identity of the detector (when a single detector is used) or the relative weights of detectors (when multiple detectors are combined). also determines how data are analyzed by a given (set of) detector(s): how classes of pairs are weighted and how the input alignment is filtered (Figure 2A). Therefore, if the protein structure is known, then can be adjusted for optimal prediction of contact pairs. The individual parameters and the optimization problem will be precisely stated later; at this point another possible formulation is given to be consistent with eq. 3:(5)A crucial assumption of this study is that the optimization in eq. 5 improves the detection of coevolving pairs within the set of contact pairs:

(6)Thus the central goal of this work is to find , which uniquely determines (Figure 2B) and ultimately . A key feature of the new framework is that the known structure plays a dual role in the current analysis. First, the structure is required for the optimization of the parameters (Eq. 5, Figure 2B bottom). Second, the structure (or some alternative conformation of that structure) is used to restrict the predicted pairs to the set of contact pairs by taking the intersection (Eq. 6).

Parameters and Procedures of the New Framework

As mentioned above, is a function of the parameter set . Now the question is: exactly what is , and how does it determine together with the data?

In general, a coevolution detector acts as a binary classifier that divides the set of all pairs into and the complementary set of pairs (the “negatives”). Given the input alignment data , the condition for classification of each pair into is that the test statistic of the detector evaluated at exceeds an adjustable threshold :(7)It is practical to constrain the number of predicted pairs at some chosen fraction of all pairs by treating as a monotonically increasing function of . Then, for a given and ,(8)Consequently, controls the true and false positive rate of the detector, which are defined subsequently in eq. 16–17.

The procedure of filtering of an alignment of homologous sequences, in particular phylogenetic type of filtering, aims to remove redundancies that emerge from the statistical non-independence within any collection of homologous sequences. These redundancies pose challenges to all coevolution detectors, especially to those assuming that homologous sequences are statistically independent from each other.

Any type of filter, applied to alignment , permutes sequences in a given order that depends on the filter type . Then the filter removes a certain number of sequences in that order. Therefore, the filtered is determined both by and by the number of sequences that remain in the alignment. It follows that, for a given ,(9)Filtering will be discussed in more detail in Methods: Alignment Filtering.

For all detectors, is known [34], [35], [38] to depend to some degree not only on the coevolution of position and (where ) but also on the overall rate of amino acid substitution at and at . The dependence on substitution rate deteriorates the performance of the detector but can, in theory, be addressed by conditioning on the rates of the pair. Therefore, the new framework incorporates a novel strategy based on the procedure of partitioning into (substitution) rate classes (Figure 2A):(10)The precise definition of will be given later (eq. 20–22), but it may be worth emphasizing at this point that the members of each are position pairs and not single positions. Now a key feature of the new framework is that can be adjusted separately for each and that is defined as the union of the resulting s:

(11) (12)

The vector thus determines every and therefore every . Like its scalar analog , is also a function of , which imposes the constraint(13)(This is the same as the constraint expressed by the second equality in eq. 8, since s are disjoint sets and thus .) The constraint in eq. 13 still allows individual s to vary, which changes the relative size (the weights) of s. In this work the procedure of changing , while requiring eq. 13 to hold, is referred to as class weighting procedure.

Partitioning also allows the filtering of separately for each rate class so that there is a separate parameter for each ,(14)and thus also depends on the vector . Eq. 14 corresponds to the combination of partitioning + class weighting + filtering in case of a general satisfying eq. 8, or to the combination of partitioning + filtering when all s are set to the same value. Note that in this case “combination” refers to procedures and not detectors.

Up to this point a single detector was assumed. Now let be a collection of detectors, and let denote their logical AND combination [43] and the corresponding thresholds (Figure 3A). Then the set of pairs predicted by the combined detector is defined as(15)It is clear that uniquely determines and that, for a given , the constraint allows individual s to vary. For some , the impact of on , relative to that of any other detector , increases with . In other words, the weight of increases in . Therefore, adjusting s relative to each other is referred to as the procedure of detector weighting and is illustrated by Figure 3A.

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Figure 3. Weighted combination of coevolution detector and .

(A) Green and orange dots represent a set of pairs of amino acid positions in a protein family. , where is the set of structural contact pairs (orange) and is the set of structurally distant pairs (green). and are coevolution detectors with statistics and , respectively, which are evaluated separately for each pair. A combined detector uses a pair of thresholds to define the set of predicted pairs (eq. 15). The set of true positives is defined as ; the true positive rate is linearly related to the number of true positives. False positives and the false positive rate are defined analogously but with instead of (eq. 16–17). Even if is fixed, (and thus ) can still vary if and change in the opposite direction. Changing at fixed is called detector weighting. For example, for all 6 thresholds marked by the arrowheads. For the threshold labeled as “equal ” the two detectors are combined in equal weights. “ more ” refers to the weight of relative to . “Only ” means that has zero weight and therefore is the same as using only. “ more ” and “only ” have analogous meanings. Finally, the threshold denoted as characterizes the optimally weighted , which by definition has the highest for each . Black circles in (B) indicate for all 6 thresholds, at , and thus report on the corresponding performance. The optimal clearly outperforms the equally weighted one, which in this case happens to perform precisely as well as “only (their circles overlap). (B-C) Obtaining for all results in receiver operating characteristic curves, which describe the performance of coevolution detectors with respect to theoretical random, and perfect, detectors. Each curve is determined by the parameter set , which includes and therefore the weights on combined detectors. Integrating a curve on yields the area , which is used as a scalar measure of performance (eq. 18, Figure 4, 5). Conditions: ; ; ; protein family  =  ABC-C; optimal phylogenetic filtering.

https://doi.org/10.1371/journal.pone.0036546.g003

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Figure 4. Influence of alignment filtering.

(A) Random filtering and phylogenetic filtering both remove sequences from the unfiltered alignment, which is represented by the large tree a, but result in trees (b and c) that differ in the length of terminal branches (red). Tree b (random filter) is similar to a in containing many extremely short terminal branches that are known to challenge coevolution detectors. In contrast, tree c (phylogenetic filter) lacks short terminal branches. (B) Opposing effects of progressively increasing strength of filtering, which leaves gradually fewer sequences in the alignment. The top graph shows, for the phylogenetic filter, the minimal sequence-sequence distance among all sequence pairs in the filtered alignment. The two lower graphs show performance, measured by , of a coevolution detector for both the phylogenetic and random filter. The first effect, specific to the phylogenetic filter, is a rise of with increasing strength of filtering (decreasing number remaining sequences). This reflects the disappearance of short terminal branches, which in turn improves performance, until a maximum is reached around 250 sequences remaining. The second effect is the deterioration of performance with increasing strength of filtering, since fewer sequences provide less information for the coevolution detector. This effect is clearly seen for the random filter regardless of the number of remaining sequences but it becomes apparent for the phylogenetic filter only with strong filtering. Conditions: detector  =  MIp; protein family  =  ABC-C. Trees were plotted using FigTree v1.3.1 (http://tree.bio.ed.ac.uk/software/figtree/).

https://doi.org/10.1371/journal.pone.0036546.g004

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Figure 5. Optimizing the prediction of coevolving position pairs.

Performance of several coevolution detectors (identified by color keys) characterized by (A) receiver operating characteristic curves and (B) partial area under these curves. Top graph in (B): low specificity (); bottom graph: high specificity (). (above magenta bars) indicates the optimally weighted detector combination CoMapMIp after partitioning, optimal filtering and optimal class weighting (Figure 2). These optimal conditions yield the parameter set (eq. 19), which determines the set of predicted coevolving pairs, presented in Figure 6 and Table 2, 3. These results were obtained from the ABC-C dataset.

https://doi.org/10.1371/journal.pone.0036546.g005

Given a specific detector , if for all other detectors (), then the weight of these detectors vanish. This special case is equivalent to using detector alone and not in combination with other s. Furthermore, in the general case it is straight-forward to combine detector weighting with partitioning + class weighting (Figure 2A). Then each scalar is replaced by a vector so that . This can be further extended with filtering.

In summary, given the parameter , data , a filter type , substitution rate classes and a set of detectors, the collection of parameters uniquely determines the set of predicted pairs in the new framework. Next, it will be discussed how the optimal is actually found, and eq. 5 will be replaced by a closely related formula. This will be followed by detailed information on and .

Optimization Using Structural Information

Let and have the same meaning as before. Let denote the set of contact pairs and the set of pairs for which and are separated by some substantial distance in 3D space, so that and are unlikely to directly interact with each other in any native conformation of the protein. and will be defined in the next subsection; for now assume that these sets are known. The true positive rate (sensitivity) and false positive rate (reverse specificity) are defined, respectively, as(16)(17)As noted after eq. 8, and are functions of , and therefore eq. 16–17, together with eq. 8, shows that makes both and . Likewise, drives both and . In general, for a given detector and . When , the detector is informative with respect to random selection. In contrast, for a theoretical random detector (Figure 3B-C, dashed line).

The receiver operator characteristic curve of a detector is a mapping that associates each with at a fixed (Figure 3B-C). The partial area under the ROC curve is the Riemann-Stieltjes integral of with respect to over the interval :(18)Thus provides a scalar measure of performance at fixed and . The interval restricts below a chosen . Small is desired when high specificity (obtaining low ) is more important than high sensitivity (achieving high ), as in the case of this study. Note that for a random detector.

Let be a relation transforming to such that . In the new framework, the optimal parameter set is defined as(19)replacing the initial formulation of the optimization problem (eq. 5). Thus, for each , a unique is obtained, which is precisely the central goal of this work (eq. 6).

In the present analysis of ABC transporters detectors were employed, and substitution rate classes were used. This gave adjustable parameters under the constraint expressed by eq. 13. In addition to this, filtering at separate for each and provided parameters and so the parameter space had a dimension of . (Note that in Figure 2A the same is used for all .) To reduce , the present work employed a heuristic optimization strategy for eq. 19, whose details are described in Text S1 (see also Figure 3, S1 and S8).

Structural Models and Contact Pairs

The set of contact pairs was defined as those pairs for which the distance separating the C atom of position from that of is less than 8Å in a structure representing the whole protein family. The set of distant pairs was defined by requiring Å. The remaining “intermediate” pairs (Å) were excluded from as in ref. [37] because a large fraction of them may be connected by chains of coevolving contact pairs [40], [42]. Thus was obtained using only and . These sets were derived separately from Sav1866 (PDB: 2HYD) [44] and CFTR (homology model [45]) representing the ABC-B and the ABC-C family, respectively.

includes the collection of optimized thresholds, which determines the set of predicted pairs (eq. 15). Next, a collection of sets of predicted contact pairs was obtained by using , which was derived from a set of structures that correspond to distinct conformations of the same protein. For the ABC-B family, this set contained Pgp in the inward (3G5U [46]) and outward-facing [47] conformation, and for the ABC-C family, CFTR in the inward [48] and outward-facing [45] conformation. Consequently, a small fraction of predicted pairs were contact pairs selectively in some but not other conformations: for these pairs but ().

Amino Acid Substitution Model and Rate Classes

The definition of rate classes requires some discussion on the amino acid substitution model used in this study. The same model also played a role in the estimation of sequence-sequence distances (which were used for alignment filtering, as explained in the next subsection), in the inference of phylogenetic trees and in the evaluation of the coevolution statistic of certain detectors. Sequence-sequence distances and trees were both estimated by maximum likelihood using RAxML v7.0.4 [49].

The substitution of amino acid residues at each position was modeled as a continuous-time Markov process with a distinct transition rate between each pair of amino acids. The transition rates used in this study were those described by the WAG-F- model [50]. In this model, the transition rates are scaled by a specific factor at each position ; the scaling factor is known as the (overall) substitution rate . In other words, the substitution rate is allowed to vary among positions (p.110 of ref. [51]). Note that substitution rate is inversely related to “residue conservation”.

Considering all positions, the collection of rates is a set of independent, identically distributed random variables. The distribution is -type with cumulative density function . Given the number of rate classes of single positions a new random variable, the discretized substitution rate , is defined as

(20)where denotes the floor function. It follows directly from definition eq. 20 that takes values on and has discrete uniform distribution with probability mass function such that ().

This uniform “prior” probability mass function can be updated, for each position , to the “posterior” the maximum likelihood estimate when an alignment and a tree is given. In this study this was done with CoMAP v1.3.0 [19] using the tree inferred from the alignment (which corresponds to an empirical Bayes approach; see p. 114 of ref. [51]). The estimated discretized substitution rate of position is defined as the mode of the posterior distribution :(21)

Given and for each position pair , the class of pairs is defined as(22)where . By the symmetry of the right side of eq. 22, so it can be required that . Then the number of classes of pairs is derived from according to . In this work and so (Figure S2).

The notation can be replaced by using any function that maps each to a unique . The present work uses the simpler notation to refer to a rate class in general (as in eq. 10), and the form to denote a specific class (e.g. ). Similarly, the symbols , and have the same meaning as , (eq. 11–12) and (eq. 14), respectively.

Multiple Sequence Alignments

A set of ABC-B and a set of ABC-C protein sequences were collected from UniProt release 15.8 using HMMER3 [52]. In both the ABC-B and ABC-C family the “full transporter” is composed of two homologous “half transporters”, each of which contains a TMD and an NBD arranged as TMD-NBD (the “-” means that the domains are on the same subunit). But there are important differences between the two families. In in most ABC-B proteins the two halves constitute separate subunits (domain arrangement: TMD1-NBD1 TMD2-NBD2) while in all ABC-C proteins the halves are covalently linked (TMD1-NBD1-TMD2-NBD2). Moreover, in ABC-B proteins the two halves TMD-NBD () are in general identical or very similar to each other but in ABC-C proteins the halves have extremely diverged from each other. For these reasons, the ABC-B sequence set contained half transporters but the ABC-C set contained full transporters.

A separate multiple alignment (Dataset S1 and S2) was made from each set using MAFFT v6.717b [53] from which all gap-containing positions were removed while keeping the remaining positions aligned. The resulting ABC-B alignment contained 1585 sequences, the ABC-C alignment 553 sequences.

Alignment Filtering

For each unfiltered alignment and filter type , a sequence of filtered alignments was generated by removing sequences, where is the number of sequences in . As mentioned above eq. 9, the type specifies the order of removal. The two types used in this work are called phylogenetic filter and random filter (Figure 4). As discussed before, the role of the phylogenetic filter employed in this work is to remove “sequence redundancies” from the alignment. In contrast, the random filter will be used to study how the performance of coevolution detectors depend on the number of aligned sequences.

In case of the random filter, the order of removal is given by a random permutation of sequences. The phylogenetic filter applies a deterministic permutation rule to the alignment before the next sequence is removed and is generated. The rule is to consider the pair-wise evolutionary distance of all sequence pairs , where and . Next, the pair that has the shortest distance is found. Note that this is the most redundant pair according to the distance measure. Next, either or is swapped with producing the new permutation. Removing the first sequence of the new permutation creates and completes the cycle. Thus is decremented by one in each iteration of the cycle.

In terms of a phylogenetic tree, a single cycle is equivalent to finding the pair of tips connected by the shortest distance and stripping away one of these tips (with its terminal branch). As this cycle is repeated, filtering becomes “stronger”, the number of sequences decreases, and the minimal sequence-sequence distance increases in the alignment (Figure 4B top graph).

To save computational time, only a subsequence of alignments were analyzed with coevolution detectors. For , was chosen to be uniformly spaced (within rounding error) between 1 and , whereas was set to corresponding to the unfiltered alignment.

Selected Coevolution Detectors

Three families of coevolution detectors were used in this study: CoMap [19], [38], mutual information (MI) [54] and CAPS [55]. The CoMap family is conceptually related to detectors in ref. [11], [14], [37]. This family contains detectors of the form CoMap--, where is either correlation or compensation; and is either simple, Grantham, polarity, volume or charge [19]. Unlike other s, simple can be combined only with correlation but not with compensation. In this work CoMap-correlation-simple is referred to as CoMap. The mutual information family contains MI [54] and MIp [31]. The CAPS family, closely related to McBASC and other detectors [13], [16], consists of CAPS and CAPS-t, where “t” denotes time correction [55].

The selected detectors strikingly differ in whether, and how, they account for the non-independence of phylogenetically related sequences. CoMap accounts for this non-independence from “first principles”. This detector considers the set of branches of a phylogenetic tree as a sample space on which, for each position , a random variable is defined, whose value is the expected number of substitutions that occurred along a given branch . For each pair the statistic of CoMap is the correlation coefficient between and . In contrast, MIp and CAPS-t uses empirical correction formulas, whereas MI and CAPS assumes statistical independence of sequences.

Another difference among detectors is related to the transition rates of the substitution process, which is intimately related to the physico-chemical similarities between amino acids. CoMap and CAPS allows realistic, heterogeneous rates by utilizing the empirical rate matrix of the WAG-F- model. MI and MIp, however, assume the same rate for all types of transition.

Unfortunately not all detectors could be applied to all alignments. The time complexity of CAPS is , where is the number of sequences in the alignment. This made alignments with intractable for CAPS in the authors’ implementation [55]. Due to a segmentation fault, CoMap v1.3.0 [19] failed to run on alignments with roughly and with many variable positions. For these reasons only MI and MIp were applied to the large () alignments of ABC-B sequences and a few variable positions, whose discretized substitution rate was typically , needed to be removed from the weakly filtered ABC-C alignments (). Consequently the size of certain rate classes, especially that of , was smaller than others.

Extent and Sources of Improvement by Optimization Procedures

Figure 5 summarizes, for the ABC-C data set, contact prediction performance under (magenta, optimal CoMapMIp) or under conditions lacking some or all of the optimization procedures. The receiver operating characteristic curves (Figure 5A) demonstrate that the relative performance under various conditions depends on the false positive rate , or reverse specificity. Consequently, the partial area under these curves reports on the relative performance in a way that depends on the upper limit of integral of with respect to (eq. 18, Figure 5B). For most optimization procedures the relative improvement in performance was greater at high specificity (, bottom bar graph) than at low specificity (, top bar graph). Importantly, is more relevant to the predicted coevolving pairs (next section) because those represent the fraction of all pairs (eq. 8), whose vast majority is not in contact (the structural model contained more distant pairs than contact pairs).

Figure 5 also demonstrates that all optimization procedures contributed to the improved performance under . At , the greatest improvement was effected by the optimally weighted combination of CoMap and MIp, relative to using either of the two detectors alone. For computational efficiency (Text S1) the remaining 9 detectors were omitted from the weighted combination. Discarding these detectors may be justified by the result that they were clearly inferior to CoMap and MIp in performance (Figure 5 and Figure S5 and S6). At low (Figure 5A) and at (Figure 5B) CoMap greatly outperformed even MIp. Despite this, the optimally weighted CoMapMIp performed markedly better than CoMap alone, which demonstrates the utility of weighted combination of detectors.

Figure 3 illustrates the principle of weighted combination of coevolution detector and , and presents performance for different relative weights. The figure takes as an example MIp and CoMap applied to substitution rate class for the ABC-C family and demonstrates that equal weighting is not in general optimal. In this case, the equally weighted failed to induce any improvement in performance (circles in Figure 3B) in comparison with using only. This result highlights the significance of (possibly unequal) detector weighting. As mentioned before, these effect were greater at low (compare Figure 3B to C).

To understand why phylogenetic filtering improved performance (Figure 5), it is useful to recall that this filter type was designed to remove the redundancies induced by closely related sequences, since these redundancies compromise the performance of all coevolution detectors. Figure 4 exemplifies the effects of alignment filtering for MIp; similar results were found for all other detectors (Figure S7 and S8). Comparing tree c to a in Figure 4A shows that strong phylogenetic filtering had a dual effect on the tree representing the alignment: (i) very short terminal branches (which indicate redundancies) disappeared but (ii) relatively few sequences remained in the alignment. The inverse relationship between effect (i) and (ii) was further established by applying the phylogenetic filter at gradually increasing strength (Figure 4B top).

Phylogenetic filtering had a dual effect also on performance (Figure 4B). Weak filtering (when the number remaining sequences was between ca. 300 and 550) improved, whereas strong filtering () deteriorated performance. Both effects were more pronounced at (bottom graph) than at (middle graph).

The dual effect of the phylogenetic filter on both tree and performance suggested that the increase in performance was related to effect (i) on the tree, whereas the decrease in performance to effect (ii). This hypothesis was tested by applying the random filter, which was designed to dissect effect (ii) from (i). In line with this design, strong random filtering did not affect the distribution of the length of terminal branches (tree b, Figure 4A). Performance (dashed lines in Figure 4B), however, deteriorated at increasing rate with respect to the strength of random filtering. This result, in agreement with the above hypothesis, suggests that the rate of performance deterioration by effect (ii) exceeds the rate of performance improvement by effect (i) at strong filtering. Therefore, optimizing phylogenetic filtering (by finding the maximum location ) is equivalent to balancing these two rates (Figure 4B, bottom).

Partitioning position pairs (explained by Figure S2) into 10 substitution rate classes amplified the filtering-induced improvement in performance particularly in the case of CoMap (Figure 5). Consistently, depended on for all detectors, especially for CoMap (see empty circles marking in Figure S8). This dependence is addressed by the combination of filtering and partitioning, which allows the conditioning of on (eq. 14).

Another benefit of partitioning was related to the possibility of weighting classes. Optimal class weighting substantially improved the performance of CoMap, MIp and MI at (Figure 5). The sources of this improvement were clarified by two further results. First, the distribution of the statistic of each detector clearly depended on (Figure S3 and S4). Second, the conditional version of the performance measure was calculated given each (Figure S7, S8 and in particular Figure S9). This uncovered the dependence of performance on substitution rate; the dependence was especially strong for CoMap. In light of these results, the advantage of class weighting is that it removes both types of dependence by conditioning threshold on (eq. 14).

Predicted Coevolving Pairs

When the fraction (eq. 8) of predicted position pairs was set to 0.001, 95 and 344 coevolving pairs were predicted for the ABC-B and ABC-C family, respectively. The roughly 4-fold difference between these numbers was due to neglecting the relatively small asymmetry between the two homologous halves of ABC-B proteins by creating an alignment from half ABC-B transporter sequences (Methods). Thus, for all pairs , both position and was restricted to the same half ABC-B transporter (this restriction was not used for ABC-C transporters, whose halves are greatly asymmetric).

The main focus of this study is not the entire set of predicted pairs but the subset , where is the set of contact pairs observed in a representative structure. For the optimization procedures, was calculated from the outward-facing Pgp and CFTR structures for the ABC-B and ABC-C family, respectively. contained 41 pairs for the ABC-B and 95 pairs for the ABC-C family. For both families the positive predictive value was an order of magnitude higher than the fraction of contact pairs in the set of all pairs. For example, for the ABC-C family whereas . Consequently, the separation between predicted pairs in -helices was distributed in a way that reflected -helical periodicity (Figure S10, Movie S1) [29], [36].

As a corollary of the unequal size of the 10 substitution rate classes together with the weighting of these classes, the size of sets was also non-uniform. Most predicted pairs fell into class (Figure S1), whose definition (eq. 22) asserts either that the discretized substitution rate at position equals 3 and or that and . As expected, relatively variable positions (exhibiting or ) clustered mainly in the 12 transmembrane helices (TM1-TM12), whereas relatively conserved positions ( or ) were typically located in the 4 intracellular loops (ICL1-ICL4) and the two NBDs, particularly at the central dimer interface (Figure 1B). The positions from which predicted pairs were composed tended to cluster also within the TM helices (Figure 1C). The latter finding, however, does not necessarily imply a natural tendency of coevolving pairs to reside in the TM helices. Rather, it can be seen as a consequence of the previous two results that link, via substitution rate, prediction sensitivity to structural localization.

For detailed exploration of the predicted coevolving pairs (Table 1, 2, 3, Dataset S5, S6), the set was considered, where and is the set of contact pairs in the outward and inward-facing conformation, respectively, of Pgp or CFTR. Thus all predicted pairs were included that were in contact in at least one of these two conformations. At the same time, , and(23)were noted, where and is the 3D distance separating pair in the outward and inward-facing conformation, respectively. Therefore, is the change of distance induced by the complete transition from the outward to the inward-facing conformation.

For the pairs of the ABC-B family (Table 1) and for those in the NBDs of the ABC-C family (Table 2 and Figure 6A) the set of interest was further narrowed to(24)where , i.e the set of pairs fulfilling the condition that and are separated by more than 4 positions in the sequence. This constraint removed “obvious” contact pairs, whose distance is constrained by primary rather than secondary to quaternary structure.

For the pairs of the TMDs of ABC-C proteins (Table 3, Figure 6B and Movie S2), a more restrictive condition was used to define the set . This means that the set(25)contains those pairs that were predicted to coevolve, for which was observed to contact in at least one conformation, and for which and localized to distinct TM helices. In this case, the notion of a “TM helix” included the helices of the ICLs since those are contiguous extensions of the sensu stricto TM helices. Figure 1A and 6 show that each of the 4 ICLs contains two helical extensions and a single “coupling helix” [44], and that pairs of ICLs form compact structural units that predominantly interact with a single NBD: with NBD1 (Figure 6A) and with NBD2. These units of 4 parallel helices are hereby termed intracellular bundle 1 and 2 consistently with the interacting NBD.

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Table 1. Coevolving Position Pairs in ABC-B transporters.

https://doi.org/10.1371/journal.pone.0036546.t001

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Table 2. Coevolving Position Pairs in the NBDs of ABC-C transporters.

https://doi.org/10.1371/journal.pone.0036546.t002

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Table 3. Coevolving Position Pairs in the TMDs of ABC-C transporters.

https://doi.org/10.1371/journal.pone.0036546.t003

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Figure 6. Coevolving position pairs in ABC-C proteins.

(A) Labeled residue side chains connected by lines form subset of predicted coevolving position pairs (eq. 24, Table 2) in NBD1, including (E474, R1066) that connects NBD1 to ICL4. Large colored numbers identify helices of the TMDs. Helix H1 of NBD1 is also labeled as in Figure 7. (B) The subset of predicted pairs (eq. 25, Table 3) are indicated in a topological map of the TMD dimer, in which 12 TM helices (large colored numbers), 2 wings and 4 intracellular loops (ICLs) are labeled. The map was obtained by cylindrical projection of the two polypeptide chains of the TMD dimer. Note that TM1-TM3 are shown twice. In both A and B the color of the lines connecting predicted pairs reports on the extent of distance change induced by the modeled outward inward conformational transition (eq. 23). Black: ; purple: ; red: .

https://doi.org/10.1371/journal.pone.0036546.g006

Pairs Involved in Conformational Changes

Comparison of the CFTR structural models in the outward and inward-facing conformation (Movie S2) revealed possible conformational transitions [48], [56]. The most striking change during the inferred outward inward transition was the dissociation of the tight dimer of NBDs, the closure of the outward-facing cleft delineated by the wings (Figure 7A) and the opening of the inward-facing cleft between the intracellular bundles (Figure 7B). While the NBDs and the lower (i.e. proximal to the NBDs) parts of the IC bundles moved as essentially rigid bodies, the upper parts of IC bundles and especially the wings appeared flexible. A prominent component of that flexibility was the translation of some TM helices along their axes relative to other helices.

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Figure 7. Position pairs evolved to regulate conformational transitions.

(A-B) The entire TMD dimer is shown in surface representation, and selected TM helices (identified by colored numbers) are displayed as ribbons. For each NBD, only helix H1 (cf. Figure 6A) is shown, as well as the bound ATP, if present. The outward-facing model conformation is characterized by a cleft between wing 1 and 2 (A, left) while the inward-facing model conformation reveals a cleft between intracellular bundle 1 and 2 (B, right). All labeled position pairs (connected by red lines and represented as spheres) were predicted as coevolving, and represent structural contact in only one of the two conformations, suggesting that these pairs evolved to regulate conformational transitions. Unlabeled pairs (black connecting lines, stick representation) are expected to remain in contact in both principal conformations, implying that they evolved to enhance structural rigidity. (A) (E873, G1003) and (Q179, V260) appear to regulate the opening and closing of the cleft between the wings and that between the IC bundles, respectively, during the conformational change. (B) (C225, P324), (F311, A876) and (L293, I942) might regulate the relative translation of TM4, TM5, TM7 and TM8 along the helical axes.

https://doi.org/10.1371/journal.pone.0036546.g007

These inferred movements during the outward inward transition were quantified by the distance change (eq. 23), whose extent is indicated by the color of the line connecting each pair in Figure 6, 7 and Movie S1, S2, S3. In Table 2, 3, Figure 6, 7 and in the main text below residues and positions are given for human CFTR (UniProt ID: CFTR_HUMAN), whereas homologous positions for 599 other ABC-C proteins can be obtained from Dataset S4. (E873, G1003) and (Q179, V260) stood out among the pairs in (and in fact also in ), for which was relatively large (Å, red lines). The uniqueness of these two pairs was established by the fact that they contributed to the structural contacts between the closed wings and IC bundles, respectively, but were separated by the cleft between the wings/bundles in the opposite conformation (Figure 7A-B, Movie S3).

For the rest of the red pairs in , position resided in the same IC bundle or wing as (Figure 6B). These included (L293, I942) in IC bundle 2, as well as (C225, P324) and (F311, A876) in wing 1. As Figure 7B and Movie 24 illustrate, the separation of these conformation-specific contact pairs was due to the inferred bending and translation of TM helix 5 with respect to TM7 and TM8. TM4 and TM5 was unusual in that they exhibited marked translation relative to each other at their extracellular ends, containing (C225, P324), whereas the same helix pair appeared relatively rigid in ICL2 (see the 4 unlabeled black and purple pairs in Figure 7B). In this regard, ICL4, formed by TM10, TM11 and a coupling helix directly interacting with NBD1, was similar to ICL2 (Figure 6B). Notably, the coupling helix of ICL4 contains R1066, which together with E474 formed the only pair in that links an NBD to any other domain; was relatively small for this pair too.