Anne-Florence Bitbol
email: anne-florence.bitbol@sorbonne-universite.fr

Inferring interaction partners from protein sequences using mutual information

Keywords: Protein-protein interactions, mutual information, inference, sequence data, information theory, paralogs, homologs

Inferring interaction partners from protein sequences using mutual information

Introduction

Most cellular processes are carried out by interacting proteins. Functional protein-protein interactions allow multi-protein complexes to assemble, and ensure proper signal transduction in various pathways. Hence, mapping functional protein-protein interactions is an important fundamental question. High-throughput experiments have recently elucidated a substantial fraction of protein-protein interactions in a few model organisms [START_REF] Rajagopala | The binary protein-protein interaction landscape of Escherichia coli[END_REF], but such experiments remain challenging. An attractive alternative is to exploit the increasingly abundant sequence data in order to identify functional protein-protein interaction partners.

The sequences of interacting proteins are correlated, both because of evolutionary constraints arising from the need to maintain physico-chemical complementarity among amino-acids in contact, and because of shared evolutionary history. The first type of correlations has received substantial interest, both within single proteins and across protein partners, as evolutionary constraints induce correlations between amino acids that are in contact in the folded protein or in the multi-protein complex. Hence, the correlations observed in multiple sequence alignments of homologous proteins contain information about protein structure. Global statistical models allow direct and indirect correlations to be disentangled [START_REF] Lapedes | Correlated mutations in models of protein sequences: phylogenetic and structural effects[END_REF][START_REF] Burger | Accurate prediction of protein-protein interactions from sequence alignments using a Bayesian method[END_REF][START_REF] Weigt | Identification of direct residue contacts in protein-protein interaction by message passing[END_REF]. Such models, built using the maximum entropy principle [START_REF] Jaynes | Information Theory and Statistical Mechanics[END_REF], and assuming pairwise interactions, known in the field of proteins as Direct Coupling Analysis (DCA), have been used with success to determine three-dimensional protein structures from sequences [START_REF] Marks | Protein 3D structure computed from evolutionary sequence variation[END_REF][START_REF] Su Lkowska | Genomics-aided structure prediction[END_REF], to analyze mutational effects [START_REF] Dwyer | Predicting functionally informative mutations in Escherichia coli BamA using evolutionary covariance analysis[END_REF][START_REF] Cheng | Toward rationally redesigning bacterial twocomponent signaling systems using coevolutionary information[END_REF][START_REF] Cheng | Connecting the Sequence-Space of Bacterial Signaling Proteins to Phenotypes Using Coevolutionary Landscapes[END_REF][START_REF] Figliuzzi | Coevolutionary Landscape Inference and the Context-Dependence of Mutations in Beta-Lactamase TEM-1[END_REF] and conformational changes [START_REF] Morcos | Direct-coupling analysis of residue coevolution captures native contacts across many protein families[END_REF][START_REF] Malinverni | Large-Scale Conformational Transitions and Dimerization Are Encoded in the Amino-Acid Sequences of Hsp70 Chaperones[END_REF], to find residue contacts between known interaction partners [START_REF] Weigt | Identification of direct residue contacts in protein-protein interaction by message passing[END_REF][START_REF] Procaccini | Dissecting the specificity of proteinprotein interaction in bacterial two-component signaling: orphans and crosstalks[END_REF][START_REF] Baldassi | Fast and accurate multivariate Gaussian modeling of protein families: predicting residue contacts and protein-interaction partners[END_REF][START_REF] Ovchinnikov | Robust and accurate prediction of residue-residue interactions across protein interfaces using evolutionary information[END_REF][START_REF] Hopf | Sequence co-evolution gives 3D contacts and structures of protein complexes[END_REF][START_REF] Tamir | Integrated strategy reveals the protein interface between cancer targets Bcl-2 and NAF-1[END_REF][START_REF] Dos Santos | Dimeric interactions and complex formation using direct coevolutionary couplings[END_REF][START_REF] Feinauer | Inter-Protein Sequence Co-Evolution Predicts Known Physical Interactions in Bacterial Ribosomes and the Trp Operon[END_REF], and most recently to predict interaction partners from sequence data [START_REF] Bitbol | Inferring interaction partners from protein sequences[END_REF][START_REF] Gueudre | Simultaneous identification of specifically interacting paralogs and interprotein contacts by direct coupling analysis[END_REF]. DCA models lay the emphasis on interactions between residues that are in direct contact in the three-dimensional protein structure. However, correlations in protein sequences also have important collective modes [START_REF] Halabi | Protein sectors: evolutionary units of three-dimensional structure[END_REF][START_REF] Rivoire | Evolution-based functional decomposition of proteins[END_REF], which can arise from functional selection [START_REF] Yan | Architecture and coevolution of allosteric materials[END_REF][START_REF] Wang | Revealing evolutionary constraints on proteins through sequence analysis[END_REF], and additional correlations are due to phylogeny [START_REF] Halabi | Protein sectors: evolutionary units of three-dimensional structure[END_REF][START_REF] Casari | A method to predict functional residues in proteins[END_REF][START_REF] Qin | Power law tails in phylogenetic systems[END_REF]. These contributions are deleterious to the prediction of contacts [START_REF] Qin | Power law tails in phylogenetic systems[END_REF] but not necessarily to the prediction of interacting partners, since a pair of interacting partners may be subject to common functional selection, and may also have a more strongly shared phylogenetic history than non-interacting proteins.

Here, we present an alternative approach to predict interaction partners from sequence data, among the paralogous proteins belonging to two interacting families. In contrast to the previous pairwise maximum entropy-based approaches [START_REF] Bitbol | Inferring interaction partners from protein sequences[END_REF][START_REF] Gueudre | Simultaneous identification of specifically interacting paralogs and interprotein contacts by direct coupling analysis[END_REF], the present method is based on an approximate maximization of mutual information between the sequences from the two protein families. Specifically, we develop a variant of the iterative pairing algorithm (IPA) introduced in [START_REF] Bitbol | Inferring interaction partners from protein sequences[END_REF], where we use mutual information (MI) as a score to maximize, instead of the effective interaction energy from a pairwise maximum entropy (DCA) model. We demonstrate that this mutual information-based algorithm (MI-IPA) performs slightly better than the one (DCA-IPA) introduced by us and colleagues in [START_REF] Bitbol | Inferring interaction partners from protein sequences[END_REF]. Our findings entail that the statistical dependences allowing interaction partner prediction from sequence data are not restricted to the contacting residue pairs revealed by DCA.

Results

We developed an iterative pairing algorithm (MI-IPA) that pairs paralogous proteins from two interacting protein families A and B by approximately maximizing mutual information between the sequences of the two families. Here, we first introduce the information theory-based pairing score we employ, before briefly explaining the steps of the MI-IPA. Next, we present the results we obtained with the MI-IPA. Throughout, we compare the performance of the MI-IPA to that obtained with the DCA-IPA from [START_REF] Bitbol | Inferring interaction partners from protein sequences[END_REF], which infers a pairwise maximum entropy model and approximately maximizes the resulting effective interaction energies. First, we consider the case where the MI-IPA starts with a training set of known protein pairs, obtaining good performance even with few training pairs. Then, we demonstrate that the MI-IPA can make accurate predictions starting without any training set, as would be needed to predict novel protein-protein interactions. Next, we assess to what extent the MI-IPA is successful at maximizing mutual information. We further demonstrate the robustness of the MI-IPA by successfully applying it to several pairs of proteins. Finally, we show how the MI-IPA reveals signatures of protein-protein interactions between two protein families.

A pairing score based on pointwise mutual information (PMI)

Consider an alignment of M concatenated sequences AB of length L, where A is a protein from family A and B is a protein from family B. At each amino-acid site i ∈ {1, .., L}, a given sequence can feature any amino acid (represented by α ∈ {1, .., 20}), or a gap (represented by α = 21). To describe the statistics of this concatenated alignment (CA), we employ the singlesite frequencies of occurrence of each state α at each site i, denoted by f i (α), and the two-site frequencies of occurrence of each ordered pair of states (α, β) at each ordered pair of sites (i, j), denoted by f ij (α, β). These empirical frequencies are obtained by counting the sequences where given residues occur at given sites and dividing by the number M of sequences in the CA. (Note that when computing frequencies from real protein data, it is useful to weight sequences so as to attenuate the impact of biased sampling, and to include pseudocounts, in order to mitigate finite-sample effects, see Methods.) The empirical frequencies constitute estimates of the corresponding probabilities p i (α) and p ij (α, β), and tend toward them in the limit where the number M of sequences in the alignment tends to infinity.

The pointwise mutual information (PMI) of a pair of residues (α, β) at a pair of sites (i, j) is defined as [START_REF] Fano | Transmission of Information: A Statistical Theory of Communications[END_REF]:

PMI ij (α, β) = log p ij (α, β) p i (α)p j (β) ≈ M ≫1 log f ij (α, β) f i (α)f j (β) . (1)
Averaging this quantity over all possible residue pairs yields the mutual information (MI) between sites i and j [START_REF] Cover | Elements of Information Theory[END_REF]:

MI ij = α,β p ij (α, β) PMI ij (α, β) = α,β p ij (α, β) log p ij (α, β) p i (α)p j (β) . ( 2
)
PMI has been used in linguistics to study the co-occurrence of words [START_REF] Church | Word association norms, mutual information, and lexicography[END_REF][START_REF] Role | Handling the impact of low frequency events on co-occurrence based measures of word similarity[END_REF]. Note that in some instances [START_REF] Fano | Transmission of Information: A Statistical Theory of Communications[END_REF][START_REF] Church | Word association norms, mutual information, and lexicography[END_REF] PMI is called MI, and MI is then referred to as the average value of MI. We define a pairing score S AB for each pair AB of proteins as the sum of the PMIs of the inter-protein pairs of sites of this concatenated sequence (i.e. those that involve one site in protein A and one site in protein B):

S AB = L A i=1 L j=L A +1 PMI ij (α i , β j ) , (3)
where we have denoted the concatenated sequence AB by (α 1 , . . . , α

L A , β L A +1 , • • • , β L ), with L A
the length of the A sequence. This score can be computed for a pair AB that is a member of the CA used to estimate the PMI of each residue pair at each site, but also for any other pair AB comprised of the sequences of members of the protein families A and B.

Next, consider a candidate assignment X of M ′ pairs AB, where each protein A is paired with a protein B from the same species, resulting in a CA of M ′ sequences of length L = L A + L B . Again, the pairs in this CA can involve the proteins in the CA of M pairs used to estimate the PMIs, with the same assignment or a different one, or any other pair AB comprised of the sequences of members of the protein families A and B. We define the overall pairing score S X of the assignment X by the average of all pairing scores (see Eq. 3) of the pairs involved:

S X = 1 M ′ AB∈X S AB = 1 M ′ L A i=1 L j=L A +1 AB∈X PMI ij (α i , β j ) = 1 M ′ L A i=1 L j=L A +1 α,β AB∈X α i =α, β j =β PMI ij (α, β) = L A i=1 L j=L A +1 α,β f ′ ij (α, β) PMI ij (α, β) , (4)
where f ′ ij (α, β) denotes the joint empirical frequency of amino acid state α at site i and amino acid state β at site j in X. It corresponds to the number of concatenated sequences in X featuring both amino acid state α at site i and amino acid state β at site j, divided by M ′ .

In the limit of large alignments, the empirical frequencies tend toward probabilities. Besides, if the scored CA is the same as the one used to calculate the PMIs, all frequencies will be the same for both of them (i.e.

f ′ ij (α, β) = f ij (α, β)).
For different CA of proteins from the same families A and B, in the case of an assignment consistent with the CA used to estimate the PMIs, the limiting two-body probabilities will be the same in the two CA. Hence, in these cases, combining Eqs. 2 and 4 yields

S X -----→ M →∞ M ′ →∞ L A i=1 L j=L A +1 MI ij . (5)
For large alignments, maximizing S X thus corresponds to maximizing the sum of the MIs of inter-protein site pairs, which is itself a pairwise approximation of the MI between the sequences of two protein families A and B. A brute-force self-consistent maximization of S X over all possible assignments X of a realistic dataset would result in a combinatorial explosion, since all allowed permutations of pairs would need to be considered, and PMIs would need to be computed for each of them. In practice, since biologically meaningful pairings have to be made within a species, this means that we would need to consider all combinations of all permutations within each species, which already yields prohibitively large numbers of assignments to test. Hence, we propose an algorithm to perform an approximate maximization of S X .

An iterative pairing algorithm (IPA) based on MI

In order to approximately maximize mutual information via the score S X (see Eq. 4), we propose an iterative pairing algorithm (referred to as MI-IPA) inspired by that of [START_REF] Bitbol | Inferring interaction partners from protein sequences[END_REF] (referred to as DCA-IPA), where the effective interaction energy from a global statistical model was approximately maximized.

In each iteration, we first estimate PMIs for all inter-protein residue pairs from a concatenated alignment (CA) of paired sequences. The initial CA, used at the first iteration, is either built from a training set of known correct protein pairs, or made from random pairs, assuming no prior knowledge of interacting pairs. We calculate the pairing scores S AB (see Eq. 3) for every possible protein pair AB within each species, by summing the inter-protein PMIs. Next, within each species, we assign pairs by maximizing the sum of S AB scores in the species (assuming one-to-one specific interactions), thereby maximizing S X (see Eq. 4) over biologically relevant pair assignments, where each protein has a partner within its species. We attribute a confidence score to each predicted pair, by using the difference of scores between the optimal assignment of pairs in the species and the best alternative assignment that does not involve this predicted pair. The CA is then updated by including the highest-scoring protein pairs, and the next iteration can begin. At each iteration, all pairs in the CA are re-selected based on confidence scores (except the initial training pairs, if any), allowing for error correction. More details on each step of the MI-IPA are given in Methods.

The MI-IPA accurately predicts interaction partners from a training set of known partners

As in [START_REF] Bitbol | Inferring interaction partners from protein sequences[END_REF], we use histidine kinases (HKs) and response regulators (RRs) from prokaryotic twocomponent signaling systems as our main benchmark. Two-component systems are important pathways that enable bacteria to sense and respond to environment signals. Typically, a transmembrane HK senses a signal, autophosphorylates, and transfers its phosphate group to its cognate RR, which in turn induces a cellular response [START_REF] Laub | Specificity in two-component signal transduction pathways[END_REF]. Importantly, most cognate HK-RR pairs are encoded in the same operon, so actual interaction partners are known, which enables us to assess performance.

Unless otherwise specified, our results were obtained on a "standard dataset" comprising 5064 HK-RR pairs for which the correct pairings are known from gene adjacency. Each species has on average m p = 11.0 pairs, and at least two pairs (see Methods).

We start by predicting interaction partners starting from a training set of known pairs. As our training set, we pick a random set of N start known HK-RR pairs from the standard dataset. The first iteration of the MI-IPA uses this concatenated alignment (CA) to compute PMIs and score possible pairs. We blind the pairings of the remaining dataset, and use it as a testing set on which we predict pairings. At each subsequent iteration n > 1, the CA used to recompute PMIs contains the initial training pairs plus the (n -1)N increment highest-scoring predicted pairs from the previous iteration (see Methods).

As in the case of the DCA-IPA [START_REF] Bitbol | Inferring interaction partners from protein sequences[END_REF], iterating, and thereby progressively adding high-scoring pairs to the CA, allows us to increase the fraction of pairs that are correctly predicted. This gradual improvement of the TP fraction during the iterations of the MI-IPA is shown in Fig. 1A for different training set sizes N start . The increase of TP fraction is especially spectacular for small training sets. Fig. 1B shows the initial TP fraction, obtained at the first iteration, and the final TP fraction, obtained at the last iteration, versus the size of the training set N start , both for the MI-IPA and for the DCA-IPA [START_REF] Bitbol | Inferring interaction partners from protein sequences[END_REF]. In both cases, comparing the initial and final TP fractions demonstrates the major interest of our iterative approach, through the massive increase in TP fraction, especially for small training sets. Moreover, the final TP fraction depends only weakly on N start : the iterative approach removes the need for large training sets. Both algorithms yield very good performance, and the MI-IPA even outperforms the DCA-IPA in the trickiest case of small training sets. In this limit (N start = 1), the MI-IPA yields 86% true positive (TP) pairs while the DCA-IPA yields 84% TP, while both start from 12% TP at the first iteration. These final TP fractions are strikingly higher than the random expectation of 9%, while such small training sets contain very little information about pairings, as illustrated by the associated low initial TP fraction.

Since the MI-IPA to the DCA-IPA yield similar performance, we asked whether they tend to predict the same correct pairs when starting from the same training set. To assess this, consider a species with m AB pairs, and denote by p (resp. q) the number of pairs correctly assigned by the MI-IPA (resp. by the DCA-IPA) in this species. If the two algorithms made independent predictions, it would correspond to independently and randomly drawing p (resp. q) proteins A among m, to be correctly paired. In this null model, the number of possible MI-IPA assignments that share k correct pairs with the DCA-IPA is q k m-q p-k : k correct pairs are chosen among the q pairs correctly assigned by the DCA-IPA, and the other p -k ones are chosen among the m -q proteins A incorrectly paired by the DCA-IPA. Besides, the total number of possible MI-IPA assignments with p correct pairs is m p . Hence, the probability P (p, k, m, q) that k correct pairs are predicted by both algorithms is given by the hypergeometric distribution:

P (p, k, m, q) = q k m-q p-k m p . ( 6
)
The expectation of k under this distribution is given by In both panels, the standard dataset of HK-RRs is used, and the CA includes N increment = 6 additional pairs at each iteration. All results are averaged over 50 replicates that differ by the random choice of HK-RR pairs in the training set. Dashed lines represent the average TP fraction obtained for random within-species HK-RR pairings.

k = p q m . (7)
Hence, in this fully independent null model, the expectation of the total number of correct pairs predicted by both algorithms can be obtained by summing the expectations in Eq. 7 over all species in the dataset. Moreover, in each species, the observed number k obs of shared correct pairs can be compared to this null model, as well as to the extreme case where all correct pairs that can be shared are shared, yielding min(p, q) shared pairs. Hence, we define the relative excess E of shared predictions by

E = S i=1 k obs,i -k i S i=1 min(p i , q i ) -k i , (8)
where the index i corresponds to a particular species, the sums run over the S species present in the dataset, and the expectations are given by Eq. 7. If E > 0, the two algorithms tend to predict the same pairs more frequently than if their predictions were fully independent. In addition, the maximal value E can take is 1, including the case where predictions from both algorithms are exactly the same. We calculated the average value of E across 50 replicates where both algorithms were started from the same training set, for various N start values, in the same conditions as in Fig. 1. As expected, we found that E increases with N start , as the algorithms share more information to begin with. However, E depends rather weakly on the size N start of the training set, varying smoothly from 53% for N start = 1 to 68% for N start = 2000. This indicates a significant tendency of the two algorithms to make the same correct predictions, even in the case of small training sets.

Great accuracy is maintained in the absence of a training set

In [START_REF] Bitbol | Inferring interaction partners from protein sequences[END_REF], we showed that the DCA-IPA yields very good identification of interacting pairs without any training set, i.e. without any prior knowledge of interacting pairs. Given this previous result, and given the success of the MI-IPA with very small training sets, we ask whether the MI-IPA also makes good predictions in the absence of a training set. To test this, we followed the approach introduced in [START_REF] Bitbol | Inferring interaction partners from protein sequences[END_REF] by randomly pairing each HK with an RR from the same species, and using these 5064 random pairs to train the initial model. At each subsequent iteration n > 1, the CA only contained the (n -1)N increment highest-scoring pairs from the previous iteration (see Methods).

Fig. 2A shows the final TP fraction obtained for different values of N increment , both for the MI-IPA and for the DCA-IPA. In both cases, the iterative method performs best for small increment steps, which highlights again the interest of the iterative approach. Importantly, the MI-IPA performs better than the DCA-IPA for all values of N increment , and requires substantially less small increments than the DCA-IPA to reach the same performance. The low-N increment limit of the final TP fraction is 0.87 for the MI-IPA, versus 0.84 for the DCA-IPA. These values are consistent with those obtained above with a single training pair, N start = 1 (Fig. 1A). We emphasize that the striking TP fraction of 0.87 is attained by the MI-IPA without any prior knowledge of HK-RR interactions. Ref. [START_REF] Bitbol | Inferring interaction partners from protein sequences[END_REF] showed that an important ingredient for the DCA-IPA to bootstrap its way toward high predictive power is that sequence similarity is favored at early iterations, which increases the TP fraction in the CA, because correct pairs have more neighbors in terms of Hamming distance than incorrect pairs. In [START_REF] Bitbol | Inferring interaction partners from protein sequences[END_REF], this was called the Anna Karenina effect, in reference to the first sentence of Tolstoy's novel. The same explanation holds for the success of the MI-IPA starting from no training set. In addition, both with MI and with DCA, when starting from random pairings, signal from actual partners should add constructively, while noise should add incoherently. To confirm that the predictive power arises from correlations and similarities in the data, we ran the MI-IPA on a version of our standard alignment of HK-RR sequences where amino acids at each site (each column) are randomly scrambled, thus removing correlations (shown in green in Fig. 1A). As expected, performance is then the same as for random within-species pairings (dashed line in Fig. 1A).

In Fig. S1A, we investigated the performance of other variants of the IPA in the same conditions as in Fig. 2A. Instead of PMIs, we explored the possibilities of using normalized PMIs

(NPMIs), NPMI ij (α, β) = -PMI ij (α, β)/ log(f ij (α, β)) [32] and covariances, C ij (α, β) = f ij (α, β) -f i (α)f j (β) (recall
that in the mean-field approximation, DCA employs the inverse of this covariance matrix). For small N increment values, we find that NPMIs perform similarly as PMIs, while covariances do significantly worse than both PMIs and DCA scores. At the first iteration, the CA is made of random within-species HK-RR pairs. At each subsequent iteration n > 1, the CA includes the (n -1)N increment top predicted pairs. All results are averaged over 50 replicates employing different initial random pairings. The MI-IPA was also run on an alignment where each column is randomly scrambled. The dashed line represents the average TP fraction obtained for random within-species HK-RR pairings. (B) Final TP fraction obtained by the MI-IPA and the DCA-IPA versus the total number M of HK-RR pairs in the dataset, starting from random pairings. For each M , except that corresponding to the full dataset, datasets are constructed by picking species randomly from the full dataset, and results are averaged over multiple different such alignments (from 50 up to 500 for small M ). For the full dataset (largest M ), averaging is done on 50 different initial random within-species pairings. All results in (B) are obtained in the small-N increment limit.

An important parameter for the performance of the MI-IPA is dataset size. Qualitatively, larger datasets imply more close neighbors, which is favorable to the success of bootstrapping, and they also allow one to estimate MI more accurately, so we expect the MI-IPA to perform best for large datasets. Indeed, Fig. 2B shows that the performance of the MI-IPA increases with dataset size. Moreover, the rise of performance occurs for slightly smaller datasets in the case of the MI-IPA than for the DCA-IPA. With DCA, a sufficiently large dataset is necessary to properly infer the pairwise maximum-entropy based global statistical model at the heart of the method [START_REF] Weigt | Identification of direct residue contacts in protein-protein interaction by message passing[END_REF][START_REF] Marks | Protein 3D structure computed from evolutionary sequence variation[END_REF][START_REF] Morcos | Direct-coupling analysis of residue coevolution captures native contacts across many protein families[END_REF]. While being data-thirsty too, the MI-IPA bypasses this particular need, and thus, it is better suited for partner prediction in smaller datasets. For the complete dataset (23,424 HK-RR pairs, see Methods), both methods reach the same striking final TP fraction of 0.93.

The fact that the MI-IPA often performs better than the DCA-IPA for predicting interacting partners stands in contrast with the fact that DCA substantially outperforms MI for residue contact prediction, and suggests that relevant covariation information is contained in pairs of residues that are not in contact. In order to test this, we ran the MI-IPA starting from random initial pairings on our standard HK-RR dataset, but suppressing all MI contributions from the contacting pairs of residues (defined with a generous 8 Å threshold on the minimum distance between amino acids). For N increment = 6, a striking TP fraction of 0.83 was obtained. This value is only slightly smaller than the 0.87 TP fraction obtained when all residue pairs are included, thus confirming that substantial MI relevant for partnership prediction is present in non-contacting residue pairs.

As above, we investigated the extent to which the MI-IPA and the DCA-IPA tend to make the same correct predictions. Note that this time, the two algorithms do not start from shared information since there is no training set. For the standard dataset of 5064 HK-RR pairs, with N increment = 6, we found that the average excess shared fraction (see Eq. 7) of correct pairs between the two algorithms is E = 53%, i.e. the same as the value obtained in the limit of a very small shared training set, N start = 1 (see above). Recall that this positive value means that the two algorithms tend to make the same correct predictions. Moreover, the bias is 53% of the maximal value it could take. It is interesting to compare this result to the bias toward shared correct predictions across different replicates of the same algorithm. We obtained E = 73% for the DCA-IPA, and E = 94% for the MI-IPA, meaning that predictions from the two different algorithms are less similar than those made by the same algorithm.

Since the two algorithms do make some different correct predictions, we next asked whether this can exploited, with the intuition that pairs that are predicted both by the MI score and by the DCA score will tend to be correct more often than other pairs. First, we considered using both scores within the IPA. We tried a combined IPA which calculates both scores at each iteration and computes separately the two corresponding pair assignments. Our ranking of pairs puts first the pairs contained in both assignments, and these pairs are ordered by decreasing MI-based confidence scores. Next come the other pairs: the assignment predicted using MI is conserved for them, and they are ordered by decreasing MI-based confidence scores. This means that the pairs consistently predicted by both methods are going to enter the CA earlier in the iterative process than those that differ. A very minor improvement was obtained over the MI-IPA using this method (see Fig. S1B). Since DCA typically requires large datasets to be reliable, we also tried combining MI and DCA by using MI at early iterations, when the CA is small, and switching to DCA at later iterations. No improvement over the MI-IPA was obtained using this method (see Fig. S1B).

Next, we tried combining final results from the DCA-IPA and the MI-IPA. In [START_REF] Bitbol | Inferring interaction partners from protein sequences[END_REF], we showed that multiple different random initializations of the DCA-IPA can be exploited to increase the TP fraction. Specifically, when ranking all possible HK-RR pairs by the fraction of replicates of the MI-IPA in which they are predicted ("replication fraction"), in decreasing order, we found a TP fraction of 0.89 among the M best-ranked pairs, a significant improvement over the 0.84 average of TP fractions obtained in individual replicates. (Here, N increment = 6 is used. Recall that M = 5064 is the total umber of pairs in the CA of HK-RRs, see Methods, and thus the number of predicted pairs in each individual replicate of the algorithm.) Applying the same strategy to the MI-IPA also yields a TP fraction of 0.89, higher than that obtained from individual replicates, 0.87. In both cases, 500 replicates of the algorithms, differing only by their initial random within-species pairings, were used to estimate replication fractions. Among the M top pairs thus predicted by the two separate methods, 87% were common, and the TP fraction among those common pairs was 0.97. We also combined the results from the 500 replicates of the DCA-IPA and the 500 replicates of the MI-IPA, and ranked pairs using their overall replication fraction. Then the TP fraction among the M best-ranked pairs is 0.91. Hence, combining final predictions from both methods yields a further improvement of performance.

The MI-IPA reaches near-maximal MI

The MI-IPA approximately maximizes MI between the sequence alignments of two protein families A and B. However, there is not guarantee that it will find the assignment with highest MI. In addition, the score we maximize converges toward the sum of MIs of all inter-protein residue pairs (henceforth called "pairwise MI") only in the limit of large alignments and assignments consistent with the CA used to estimate the PMIs (see Eq. 5). In practice, how well does the MI-IPA approach the goal of maximizing pairwise MI?

In order to answer this question, we now compare the pairwise MI of the pairs assigned by the MI-IPA to that of the actual protein pairs, and to the pairwise MI of the random within-species assignment which is used to initialize the MI-IPA in the absence of a training set. Fig. 3A shows these three quantities as a function of dataset size M , for HK-RR datasets to which the MI-IPA is applied starting from no training set, as in Fig. 2B. We observe a global trend of all computed pairwise MIs to decrease when M is increased. This arises from a well-known finite size effect that occurs when estimating entropies from real datasets [START_REF] Bialek | Biophysics: Searching for Principles[END_REF], and thus affects entropy-derived quantities such as MIs. To illustrate this point, Fig. 3A also shows the pairwise MI of HK-RR datasets where each column of the alignment is randomly scrambled, thus destroying actual correlations while retaining finite-size noise, as well as one-body frequencies. This null model features a similar decreasing trend as the other curves, thus demonstrating that this trend comes from finite-size effects. Note also that the scrambled alignment features pairwise MI values close to those of the initial random assignment, which makes sense because in both cases inter-protein residue pairs are decorrelated by the scrambling. The slightly higher MI of the initial random assignment arises from the fact that random partners are chosen within each species in this case, while complete columns are scrambled in the null model. Apart from the downward trend, a striking observation from Fig. 3A is that the pairwise MIs of the assignments predicted by the MI-IPA are significantly higher than those of the initial random assignments, and close to those of the actual protein pairs. Fig. 3B highlights these points by considering the excess pairwise MI in the actual protein pairs versus the initial and final assignments in the MI-IPA. Interestingly, even for small datasets, where the MI-IPA yields small TP fractions (see Fig. 2B), the pairwise MI of the assignment predicted by the MI-IPA is much closer to that of the actual assignment than to that of the initial random assignment. This suggests that the MI-IPA does a good job at maximizing pairwise MI, even though it does not reach the absolute maximum.

So far, we have used naive estimates of the pairwise MI, employing empirical frequencies instead of probabilities in MI (see Eq. 2), without correcting for the finite-size effect discussed above and visible in Fig. 3A. Various approaches have been proposed in order to correct for the systematic error due to finite-size effects in entropy (and thus MI) estimates. One can use the fact that these finite-size effects have a leading-order contribution of the form MI-MI ∝ K/M , where MI is the actual value of the mutual information and MI is its naive estimate using empirical frequencies, while K is the number of independent values that can be taken by the pair of random variables considered, and M represents dataset size [START_REF] Bialek | Biophysics: Searching for Principles[END_REF]. Estimates MI obtained by subsampling the initial dataset can be linearly (or polynomially, to include higher-order corrections) fitted versus 1/M , yielding the actual MI as the intercept. However, if M is not larger than K, subleading corrections become too important, and this method becomes inaccurate. Here, K can go up to the square of the number of possible states per site, i.e. ∼ 400. Consistently, we found that polynomial fits become poor when considering datasets with less than ∼ 1000 sequences. Fig. S2A presents results obtained using corrections from these fits, for datasets of at least 1000 sequences. A more sophisticated method to reliably estimating the entropies of discrete distributions was introduced in [START_REF] Nemenman | Entropy and inference, revisited[END_REF]. It employs a Bayesian approach and a flat prior on the entropy, and results in a correction of finite-size bias in entropy estimates. This NSB (Nemenman-Shafee-Bialek) correction is successful in undersampled cases [START_REF] Nemenman | Entropy and inference, revisited[END_REF][START_REF] Nemenman | Entropy and information in neural spike trains: progress on the sampling problem[END_REF]. In Fig. S2B, we show results obtained on our datasets using the NSB correction. Employing either of these two corrections of finite-size effects on the MI estimates (see Fig. S2) confirms our previous conclusion. For all dataset sizes, the pairwise MI of the assignment predicted by the MI-IPA is significantly closer to that of the actual assignment than to that of the initial random assignment, thus confirming that the MI-IPA yields results with near-maximal pairwise MI.

Performance of the MI-IPA is robust across various protein pairs

To demonstrate the generality of the MI-IPA, we applied it to several pairs of protein families [START_REF] Ovchinnikov | Robust and accurate prediction of residue-residue interactions across protein interfaces using evolutionary information[END_REF], beyond HK-RRs. First, we considered several protein pairs involved in ABC transporter complexes, which permit the translocation of different substances across cell membranes [START_REF] Rees | ABC transporters: the power to change[END_REF]. We built alignments of homologs of the Escherichia coli interacting protein pairs MALG-MALK, FBPB-FBPC, and GSIC-GSID, all involved in ABC transporter complexes, using the same method as in [START_REF] Bitbol | Inferring interaction partners from protein sequences[END_REF] (see Methods). As in the case of HK-RRs, for each of these pairs of protein families, we worked on subsets of ∼ 5000 protein pairs from species containing at least two pairs. In addition, we considered smaller families of proteins, yielding ∼ 2000 pairs. More specifically, we chose two pairs of proteins involved in enzymatic complexes: PAAH-PAAJ is a pair of proteins involved in the fatty acid β-oxidation multienzyme complex [START_REF] Ishikawa | Structural basis for channelling mechanism of a fatty acid beta-oxidation multienzyme complex[END_REF], and XDHC-XDHA is a pair of proteins involved in the xanthine dehydrogenase complex [START_REF] Dietzel | Mechanism of Substrate and Inhibitor Binding of Rhodobacter capsulatus Xanthine Dehydrogenase[END_REF].

Fig. 4 shows the performance of the MI-IPA and of the DCA-IPA for these protein pairs, starting from no initial training set, for various values of N increment . In all cases, we find that both methods perform very well for small N increment , but that higher performances are reached by the MI-IPA, particularly at larger N increment values. This is consistent with the results obtained with HK-RRs (see Fig. 2A). Note that here, compared to HK-RRs, larger N increment values suffice to obtain good performance. This is due to the fact that the pairing tasks are easier here, since the average number of proteins pairs per species is smaller (see Fig. 4, to be compared to m p = 11.0 for HK-RRs). (D,E) Pairs involved in enzymatic complexes; full datasets of ∼ 2000 pairs. In each case, the mean number m p of pairs per species is indicated, and dashed lines represent the average TP fraction obtained for random within-species pairings. All results are averaged over 50 replicates that differ in their initial random within-species pairings.

These accurate predictions demonstrate the broad applicability of the MI-IPA. Note that HK-RRs interact transiently, while the ABC transporter and enzymatic proteins we considered form permanent complexes, which highlights the generality of the method.

Predicting protein-protein interactions

As with the DCA-IPA [START_REF] Bitbol | Inferring interaction partners from protein sequences[END_REF], an important potential application of the MI-IPA is to predict new protein-protein interactions between two protein families. A signature of interactions based on exploiting multiple different random initializations of the algorithm was introduced in [START_REF] Bitbol | Inferring interaction partners from protein sequences[END_REF], and can be applied to the MI-IPA. Briefly, if two given proteins from the families considered are paired by the algorithm for any initial random within-species pairing, it indicates that these two particular proteins may interact. Hence, if many pairs are in this case, it constitutes a clue that proteins from these two families generally interact. In Fig. 5, the distributions of the fraction of replicates of the MI-IPA in which each particular pair of proteins is predicted ("replication fraction") are shown for two interacting pairs of protein families with similar mean number m p of pairs per species: the subset of HK-RRs homologous to BASS-BASR, and the homologs of the interacting ABC transporter proteins MALG-MALK. A pair with no known interaction, composed of homologs of BASR-MALK, is also considered. For both interacting pairs of protein families (Fig. 5A-B), the distribution of replication fractions is strongly bimodal: it favors values close to 0 (mostly corresponding to wrong pairs), and close to 1 (mostly corresponding to actual cognate pairs). No such strong bimodality is observed for the pair with no known interaction (Fig. 5C), and the observed distribution of replication fractions is closer to the null model constructed by randomly scrambling the amino acids at each site (column) of the alignment, thus removing correlations. Hence, as for the DCA-IPA [START_REF] Bitbol | Inferring interaction partners from protein sequences[END_REF], the distribution of the replication fraction presents a distinctive shape for interacting protein families. Note that the bimodality observed for the interacting pairs is even stronger with the MI-IPA than with the DCA-IPA, consistent with the better predictive power of the MI-IPA. However, the non-interacting pair was even closer to the null model with the DCA-IPA than it is here.

Figure 5: Signature of protein-protein interactions. Distribution of the fraction of MI-IPA replicates where each possible within-species protein pair is predicted as a pair, for three different pairs of protein families. The MI-IPA replicates only differ in their initial random pairings. (A and B) Interacting protein families: BASS-BASR homologs and MALG-MALK homologs. (C) Protein families with no known interaction (BASR-MALK homologs). In all panels, the distribution of replication fraction obtained on actual sequence alignments is shown together with the same distribution obtained by running the MI-IPA on alignments where each column is randomly scrambled (null model). All alignments include ∼ 5000 protein pairs, with average number of pairs per species m p ≈ 5, and each distribution is estimated from 500 MI-IPA replicates, using N increment = 50.

Another possible approach to identifying interacting protein families is to focus on MI scores, with the intuition that MI between interacting protein families should be higher than those between non-interacting families. DCA-based approaches [START_REF] Ovchinnikov | Robust and accurate prediction of residue-residue interactions across protein interfaces using evolutionary information[END_REF][START_REF] Feinauer | Inter-Protein Sequence Co-Evolution Predicts Known Physical Interactions in Bacterial Ribosomes and the Trp Operon[END_REF] to distinguishing interacting protein families from non-interacting ones have employed such strategies, relying on the average value of the DCA interaction scores over the top-ranking inter-protein pairs of residues. If this value is higher than a threshold, the protein families are predicted to interact, and the top-ranking pairs of residues become predicted contacts. In [START_REF] Gueudre | Simultaneous identification of specifically interacting paralogs and interprotein contacts by direct coupling analysis[END_REF], this strategy was combined with a solution of the paralog-pairing problem. In [START_REF] Bitbol | Inferring interaction partners from protein sequences[END_REF], a variant based on detecting outliers in DCA interaction scores was proposed and combined with the DCA-IPA. In Fig. S3, we show the MI of all inter-protein residue pairs, averaged over 500 replicates of the IPA, and then ranked by decreasing MI value, for two interacting pairs of protein families (BASS-BASR and MALG-MALK, see above) and one without known interactions (BASR-MALK). These plots do not feature strong outliers in the case of interacting protein families, contrary to the DCA interaction scores studied in [START_REF] Bitbol | Inferring interaction partners from protein sequences[END_REF]. However, MI values are substantially larger for the interacting pairs of protein families than for the non-interacting one. This suggests that a score based on MI values, e.g. their average across all inter-protein residue pairs, or across the top-ranked ones, could be employed to distinguish interacting from non-interacting pairs of protein families. In order to predict new protein-protein interactions, one would need to set a threshold on this score, for instance by assessing its performance on known interacting and non-interacting protein families, as for the DCA-based approaches [START_REF] Ovchinnikov | Robust and accurate prediction of residue-residue interactions across protein interfaces using evolutionary information[END_REF][START_REF] Feinauer | Inter-Protein Sequence Co-Evolution Predicts Known Physical Interactions in Bacterial Ribosomes and the Trp Operon[END_REF].

The absence of outliers in MI in Fig. S3 for the interacting pairs of protein families likely stems from the fact that interactions between amino acids can cause indirect correlations, which are disentangled by DCA but not by MI. Indeed, the outliers in DCA scores found in [START_REF] Bitbol | Inferring interaction partners from protein sequences[END_REF] mainly correspond to contacting residue pairs. DCA performs better than MI at identifying contacting residue pairs [START_REF] Weigt | Identification of direct residue contacts in protein-protein interaction by message passing[END_REF][START_REF] Marks | Protein 3D structure computed from evolutionary sequence variation[END_REF][START_REF] Morcos | Direct-coupling analysis of residue coevolution captures native contacts across many protein families[END_REF], and the top-ranking MI residue pairs are less likely to be actual contacts than the top-ranking DCA residue pairs (see Fig. S4). Note that the performance of both DCA and MI at uncovering interacting residue pairs can often be improved by the average product correction (APC) [START_REF] Dunn | Mutual information without the influence of phylogeny or entropy dramatically improves residue contact prediction[END_REF][START_REF] Jones | PSICOV: precise structural contact prediction using sparse inverse covariance estimation on large multiple sequence alignments[END_REF][START_REF] Ekeberg | Improved contact prediction in proteins: using pseudolikelihoods to infer Potts models[END_REF]. However, this correction did not substantially improve inter-protein contact prediction in HK-RRs or yield strong outliers in MI, so we did not employ it.

Given the success of the MI-IPA at providing global signatures of interactions between two protein families and at globally pairing interacting paralogs, we can also ask whether MI-based methods can accurately predict partners of isolated proteins. DCA-based methods have been used to predict interacting partners of "orphan" HKs and RRs involved in prokaryotic twocomponent signaling systems [START_REF] Procaccini | Dissecting the specificity of proteinprotein interaction in bacterial two-component signaling: orphans and crosstalks[END_REF][START_REF] Baldassi | Fast and accurate multivariate Gaussian modeling of protein families: predicting residue contacts and protein-interaction partners[END_REF]. These particular proteins are not encoded on the genome close to their cognate partners, contrary to the vast majority of HKs and RRs. Crosstalk between non-cognate HKs and RRs, which exists in a few cases [START_REF] Laub | Specificity in two-component signal transduction pathways[END_REF], has also been studied using DCA [START_REF] Cheng | Connecting the Sequence-Space of Bacterial Signaling Proteins to Phenotypes Using Coevolutionary Landscapes[END_REF][START_REF] Procaccini | Dissecting the specificity of proteinprotein interaction in bacterial two-component signaling: orphans and crosstalks[END_REF]. Given the striking success of the MI-IPA at identifying cognate interaction partners among HK-RRs, we asked whether MI-based scores can predict orphan partnership and crosstalk. Table S1 reports the performance of MI and DCA-based scores at predicting known orphan and crosstalk partners using only sequences. Specifically, PMI scores (see Eq. 1) and DCA-based scores (direct couplings [START_REF] Bitbol | Inferring interaction partners from protein sequences[END_REF]) were calculated from the full alignment of cognate HK-RR pairs (which excludes orphans, see Methods). Next, new alignments including both cognate and orphan proteins, and RRs paired with non-HisKA HKs (see Methods), were considered, focusing on species where orphan and crosstalk partners were reported in the literature, either from in vivo or from in vitro experimental assays. For each HK considered, all RRs in the species were ranked using either MI-based scores (see Eq. 3) or DCA-based scores (effective interaction energies [START_REF] Bitbol | Inferring interaction partners from protein sequences[END_REF]) for partnership with this HK. With both methods, known partners often ranked among the top potential partners, showing that DCA and MI capture well these interactions. Besides, we checked that our DCA-based predictions were similar to previous ones reported in [START_REF] Procaccini | Dissecting the specificity of proteinprotein interaction in bacterial two-component signaling: orphans and crosstalks[END_REF] (see Table S1). Out of the 27 known orphan and/or crosstalk pairs we considered, 7 obtained the same rank among all possible partners with MI as with DCA, 15 a lower rank, and 5 a higher rank (see Table S1). Hence, our MI-based approach performs quite well at predicting orphan and crosstalk partners, but slightly less well than DCA.

Discussion

Here, we introduced a method based on MI to predict interacting partners among the paralogs of two protein families, starting from sequences only. Specifically, our iterative pairing algorithm (MI-IPA) finds an assignment of protein pairs that approximately maximizes the MI between the sequences of the two protein families. We demonstrated that the MI-IPA allows one to predict interacting protein pairs with high accuracy, starting from sequence data only. High performance is obtained even in the absence of an initial training set of known interacting pairs. We also showed how the MI-IPA could be employed to discover new protein-protein interactions from sequence data.

In [START_REF] Bitbol | Inferring interaction partners from protein sequences[END_REF], us and colleagues introduced a similar iterative pairing algorithm (DCA-IPA) that approximately maximized an effective interaction energy between proteins, instead of MI. This effective interaction energy was calculated from a global statistical model, more specifically a pairwise maximum entropy model, approximately inferred from the empirical one and two-body frequencies of the sequence data. Such models, often called Direct Coupling Analysis (DCA), have been successful at predicting amino-acid pairs that are in contact in folded proteins [START_REF] Weigt | Identification of direct residue contacts in protein-protein interaction by message passing[END_REF][START_REF] Marks | Protein 3D structure computed from evolutionary sequence variation[END_REF][START_REF] Morcos | Direct-coupling analysis of residue coevolution captures native contacts across many protein families[END_REF], and have permitted prediction of the three-dimensional structure of proteins from sequence data [START_REF] Weigt | Identification of direct residue contacts in protein-protein interaction by message passing[END_REF][START_REF] Marks | Protein 3D structure computed from evolutionary sequence variation[END_REF][START_REF] Su Lkowska | Genomics-aided structure prediction[END_REF][START_REF] Morcos | Direct-coupling analysis of residue coevolution captures native contacts across many protein families[END_REF][START_REF] Ovchinnikov | Robust and accurate prediction of residue-residue interactions across protein interfaces using evolutionary information[END_REF][START_REF] Ovchinnikov | Large-scale determination of previously unsolved protein structures using evolutionary information[END_REF]. For structure prediction, these global statistical models outperform the use of MI [START_REF] Weigt | Identification of direct residue contacts in protein-protein interaction by message passing[END_REF][START_REF] Marks | Protein 3D structure computed from evolutionary sequence variation[END_REF], and have yielded major progress in the field, despite the promising results obtained by MI-based methods implementing corrections for background MI [START_REF] Dunn | Mutual information without the influence of phylogeny or entropy dramatically improves residue contact prediction[END_REF]. In the specific case of interacting proteins, DCA allows one to simultaneously infer interaction partners and structural contacts between them [START_REF] Gueudre | Simultaneous identification of specifically interacting paralogs and interprotein contacts by direct coupling analysis[END_REF].

Here, we demonstrated that the MI-IPA performs at least as well as the DCA-IPA, and often outperforms it. The MI-IPA is also faster (see Methods) and requires fewer iterations to reach good performance. Our results highlight an interesting difference between the prediction of contacting pairs of residues, where DCA substantially outperforms MI, and the prediction of interacting partners, where MI often outperforms DCA. Crucially, we have achieved an accurate identification of interacting partners without the need to build a global statistical model of the sequence data highlighting effective pairwise interactions between contacting residues, which is at the root of DCA. An important motivation underlying DCA and other related approaches, such as the Bayesian network method of [START_REF] Burger | Accurate prediction of protein-protein interactions from sequence alignments using a Bayesian method[END_REF] and the sparse inverse covariance estimation of [START_REF] Jones | PSICOV: precise structural contact prediction using sparse inverse covariance estimation on large multiple sequence alignments[END_REF], is to disentangle correlations arising caused by direct interactions from indirect correlations due to a chain of couplings [START_REF] Lapedes | Correlated mutations in models of protein sequences: phylogenetic and structural effects[END_REF][START_REF] Burger | Accurate prediction of protein-protein interactions from sequence alignments using a Bayesian method[END_REF][START_REF] Weigt | Identification of direct residue contacts in protein-protein interaction by message passing[END_REF][START_REF] Jones | PSICOV: precise structural contact prediction using sparse inverse covariance estimation on large multiple sequence alignments[END_REF]. While this distinction between direct and indirect correlations is crucial to infer which residues are in contact in the folded protein, its importance is probably reduced in the determination of interacting partners among paralogous proteins. Besides, MI has been extremely successful in determining "specificity residues" crucial in the interactions between HKs and their cognate RRs: these residues are those involved in interprotein residue pairs with highest MIs [START_REF] Skerker | Rewiring the specificity of two-component signal transduction systems[END_REF][START_REF] Podgornaia | Determinants of specificity in two-component signal transduction[END_REF][START_REF] Podgornaia | Protein evolution. Pervasive degeneracy and epistasis in a protein-protein interface[END_REF]. Strikingly, mutating only these few specificity residues has allowed to fully switch the specificity of HK-RR pairs [START_REF] Skerker | Rewiring the specificity of two-component signal transduction systems[END_REF]. In addition, MI has revealed a cluster of co-evolving dynamic residues that are not in direct contact, but that are important to interaction specificity, through their involvement in the conformational arrangement of the active site residues both in the HK and in the RR, as confirmed by mutagenesis and NMR studies [START_REF] Szurmant | Coevolving motions at protein-protein interfaces of two-component signaling systems identified by covariance analysis[END_REF]. Hence, our present results reinforce the conclusion of these earlier studies, showing that MI can accurately reveal the specificity of interacting protein pairs. Finally, the MI between two sequence sites simply measures how much observing a residue at one site tells us about the other. It is entirely agnostic regarding the origin of this statistical dependence, and incorporates contributions from phylogeny [START_REF] Casari | A method to predict functional residues in proteins[END_REF][START_REF] Qin | Power law tails in phylogenetic systems[END_REF] as well as from those arising from functional selection. In addition, the latter contributions include those arising directly or indirectly from structural contacts, but might also involve aspects of protein function other than structural ones, including collective correlations between residues [START_REF] Halabi | Protein sectors: evolutionary units of three-dimensional structure[END_REF][START_REF] Rivoire | Evolution-based functional decomposition of proteins[END_REF][START_REF] Yan | Architecture and coevolution of allosteric materials[END_REF][START_REF] Wang | Revealing evolutionary constraints on proteins through sequence analysis[END_REF][START_REF] Cunningham | Coupling between Protein Stability and Catalytic Activity Determines Pathogenicity of G6PD Variants[END_REF]. This stands in contrast with DCA, which lays the emphasis on interactions between residues that are in direct contact in the three-dimensional protein structure. The fact that the MI-IPA slightly outperforms the DCA-IPA thus constitutes a hint that sources of covariation other than those that maintain structural contacts help to identify interaction partners among paralogs. This might be due to multiple effects. In particular, global functional selection could potentially affect both interacting partners together. For instance, a functionally important mechanical deformation of the complex formed by these partners could be subject to selection, yielding collective correlations that extend in the sequences of both partners [START_REF] Yan | Architecture and coevolution of allosteric materials[END_REF][START_REF] Wang | Revealing evolutionary constraints on proteins through sequence analysis[END_REF]. Besides, interacting partners may also share more phylogenetic history than non-interacting proteins, e.g. if the genes encoding the partners tend to be duplicated and/or horizontally transferred together [START_REF] Alm | The evolution of two-component systems in bacteria reveals different strategies for niche adaptation[END_REF][START_REF] Capra | Evolution of two-component signal transduction systems[END_REF][START_REF] Rowland | Crosstalk and the evolution of specificity in two-component signaling[END_REF]. Then, phylogenetic correlations could aid the prediction of interacting partners among paralogous proteins, despite being deleterious to residue contact identification [START_REF] Qin | Power law tails in phylogenetic systems[END_REF]. Next, it will be interesting to further investigate these various sources of covariation, both functional and phylogenetic. This should be useful for the particular problem of prediction of interacting protein pairs, as well as for the more general understanding of the sequence-function relationship in proteins.

Methods

Here, we first explain the different steps of the MI-IPA and briefly discuss its run time. Next, we describe the datasets used and the way the statistics of these alignments are computed. Matlab implementations of both the MI-IPA and the DCA-IPA on our standard HK-RR dataset are freely available at https://doi.org/10.5281/zenodo.1421781 and https://doi.org/10.5281/zenodo.1421861 respectively. The various sequence datasets discussed here are available as S1 Data at https://doi.org/10.1371/journal.pcbi.1006401.

Iterative pairing algorithm based on mutual information (MI-IPA)

Ref. [START_REF] Bitbol | Inferring interaction partners from protein sequences[END_REF] introduced an iterative pairing algorithm (IPA) to predict interaction partners among paralogs from two protein families. It essentially performs an approximate maximization of the average effective interaction energy between pairs of proteins comprised of one protein of family A and one of family B, and the effective interaction energy is calculated from a pairwise maximum entropy model (see also [START_REF] Gueudre | Simultaneous identification of specifically interacting paralogs and interprotein contacts by direct coupling analysis[END_REF]). Here, we propose the MI-IPA, a variant of the IPA that approximately maximizes MI via S X (see Eq. 4), instead of the effective interaction energy. Importantly, the MI-IPA does not require the construction of a global statistical model of the data, contrary to the algorithms from Refs. [START_REF] Bitbol | Inferring interaction partners from protein sequences[END_REF][START_REF] Gueudre | Simultaneous identification of specifically interacting paralogs and interprotein contacts by direct coupling analysis[END_REF]. Let us now describe each step of an iteration of the MI-IPA, after explaining how the CA is constructed for the very first iteration.

For simplicity, we assume that in each species, there is the same number of proteins of family A and of family B. If this is not the case, an injective matching strategy can be used [START_REF] Gueudre | Simultaneous identification of specifically interacting paralogs and interprotein contacts by direct coupling analysis[END_REF], so in each species, the relevant number of proteins is the minimum of the number of proteins of family A and of proteins of family B.

Initialization of the CA

If starting from a training set of known interaction partners AB, the CA for the first iteration of the IPA is built from the pairs AB in this training set. In subsequent iterations, the training set pairs are always kept in the CA, and additional pairs with the highest confidence scores (see below) are added to the CA.

In the absence of a training set, each protein A of the dataset where we wish to predict pairings is randomly paired with a protein B from its species. All these random pairs are included in the CA for the first iteration of the MI-IPA. Hence, this initial CA contains a mixture of correct and incorrect pairs, with one correct pair per species on average. At the second iteration, the CA is built using only the N increment AB pairs with the highest confidence scores obtained from this first iteration. Now that we have explained the initial construction of the CA, let us describe each step of an iteration of the MI-IPA.

Calculation of pointwise mutual informations (PMI)

Each iteration begins by the calculation of PMI scores from the CA of paired AB sequences. The empirical one-and two-site frequencies, f i (α) and f ij (α, β), of occurrence of amino-acid states α (or β) at each site i (or j) are computed for the CA, using a specific weighting of similar sequences, and a pseudocount correction (see below) [START_REF] Weigt | Identification of direct residue contacts in protein-protein interaction by message passing[END_REF][START_REF] Marks | Protein 3D structure computed from evolutionary sequence variation[END_REF][START_REF] Morcos | Direct-coupling analysis of residue coevolution captures native contacts across many protein families[END_REF][START_REF] Procaccini | Dissecting the specificity of proteinprotein interaction in bacterial two-component signaling: orphans and crosstalks[END_REF]. The PMI of each residue pair (α, β) at each pair of sites (i, j), defined in Eq. 1, is then estimated from these frequencies as

PMI ij (α, β) = log f ij (α, β) f i (α)f j (β) . (9)

Calculation of pairing scores for all possible AB pairs

Having calculated PMIs on the CA, we next turn to the dataset where we wish to predict pairings. The MI-based pairing score S AB of each possible AB pair within each species in the dataset is calculated by summing all the inter-protein PMIs involved, as defined in Eq. 3.

AB pair assignments and ranking by confidence score

We make one-to-one AB pairs within each species in the dataset by maximizing the sum of the scores of all pairs in this species. Considering the matrix of the scores of all possible pairs, this amounts to choosing one element per line and per column, such that the sum of all of them is maximal. This assignment problem is solved exactly and efficiently using the Hungarian algorithm (also known as the Munkres algorithm) [START_REF] Kuhn | The Hungarian Method for the assignment problem[END_REF][START_REF] Munkres | Algorithms for the Assignment and Transportation Problems[END_REF][START_REF] Buehren | Functions for the rectangular assignment problem[END_REF]. Each pair is scored by a confidence score ∆S AB , which is the difference between the sum of the scores of the assigned pairs in the species and that obtained by using the Hungarian algorithm again while this pair is disallowed.

Once pairs are made and confidence scores are calculated in each species, all the assigned AB pairs from all species are ranked in order of decreasing confidence score. Note that in [START_REF] Bitbol | Inferring interaction partners from protein sequences[END_REF], we had used an approximate greedy approach to the assignment problem, where the pair with lowest energy is selected first, and the two proteins involved are removed from further consideration, and the process is repeated until all As and Bs are paired. For the DCA-IPA of [START_REF] Bitbol | Inferring interaction partners from protein sequences[END_REF], the greedy algorithm performed marginally better than the Hungarian algorithm, while in the present MI-IPA, the Hungarian algorithm yields better performance (see Fig. S5). This slight difference may be explained by the fact the PMI scores are most meaningful collectively, as their sum tends to a pairwise approximation of MI (see Eq. 5).

Update of the CA

The ranking of the AB pairs is used to pick those pairs that are included in the CA at the next iteration. Pairs with a high confidence score are more likely to be correct because there was less ambiguity in the assignment. The number of pairs in the CA is increased by N increment at each iteration, and the MI-IPA is run until all the As and Bs in the dataset have been paired and added to the CA. In the last iteration, all pairs assigned at the second to last iteration are included in the CA.

If starting from a training set of AB pairs, the N start training pairs remain in the CA throughout. The As and Bs from all the other pairs in the CA are re-paired and re-scored at each iteration, and only re-enter the CA if their confidence score is sufficiently high. In other words, at the first iteration, the CA only contains the N start training pairs. Then, for any iteration number n > 1, it contains these exact same N start training pairs, plus the (n -1)N increment assigned AB pairs that had the highest confidence scores at iteration number n -1.

In the absence of a training set, all As and Bs in the dataset are paired and scored at each iteration, and all the pairs of the CA are fully re-picked at each iteration based on the confidence score. For any iteration number n > 1, the CA contains the (n -1)N increment assigned AB pairs that had the highest confidence scores at iteration number n -1.

Once the new CA is constructed, the next iteration can start.

Run time

The run time of the MI-IPA strongly depends on the number of iterations made, as well as on dataset size (length of sequences, number of sequences in the dataset). Because the inversion of the correlation matrix is not necessary, the MI-IPA is substantially faster than the DCA-IPA, and this difference is stronger for longer sequences. In practice, for our standard HK-RR dataset, one round of calculation of all PMI scores from sequence data (Eq. 1) takes 7s, while for DCA, the calculation of all pairwise couplings takes 9s. The longest concatenated sequences considered here were XDHA-XDHC (see Fig. 4), with total length L = 911 (to be compared to 176 for HK-RRs): in this case, the calculation of all PMI scores takes 135s, while for DCA, the calculation of pairwise couplings takes 289s. Times are given for a single processor (Intel Core i7) of a standard laptop computer. Regarding the full MI-IPA, run times with all standard datasets used here were under 10 hours.

Comparison to the DCA-IPA Throughout, the DCA-IPA was used as described in [START_REF] Bitbol | Inferring interaction partners from protein sequences[END_REF], with the same parameters (threshold Hamming distance for weighting neighboring sequences θ = 0.3, pseudocount weight Λ = 0.5).

Dataset construction

We use the HK-RR datasets described in [START_REF] Bitbol | Inferring interaction partners from protein sequences[END_REF]. Briefly, the complete dataset was built using the online database P2CS [START_REF] Barakat | P2CS: a twocomponent system resource for prokaryotic signal transduction research[END_REF][START_REF] Ortet | P2CS: updates of the prokaryotic two-component systems database[END_REF], yielding a total of 23,424 HK-RR pairs (known by genome proximity) from 2102 different species. We focused on the protein domains through which HKs and RRs interact, which are the Pfam HisKA domain present in most HKs (64 amino acids) and the Pfam Response reg domain present in all RRs (112 amino acids).

In most of the paper, and as in [START_REF] Bitbol | Inferring interaction partners from protein sequences[END_REF], we focused on a smaller "standard dataset" extracted from this complete dataset, both because protein families that possess as many members as the HKs and RRs are atypical, and in view of computational time constraints. This standard dataset was constructed by picking species randomly, and comprises 5064 pairs from 459 species. In our datasets, we discarded sequences from species that contain only one pair, for which pairing is unambiguous. It allows us to assess the impact of training set size (N start ) without including an implicit training set via these pairs, and it also enables us to address prediction in the absence of any known pairs (no training set).

While we used HK-RRs as the main benchmark for the MI-IPA, we assessed the generality of its performance by applying it to several other pairs of protein families. For these proteins, paired alignments of homologs of known Escherichia coli interacting protein pairs were built using a method adapted from [START_REF] Ovchinnikov | Robust and accurate prediction of residue-residue interactions across protein interfaces using evolutionary information[END_REF], as detailed in [START_REF] Bitbol | Inferring interaction partners from protein sequences[END_REF]. Note that we kept full sequences, without restricting them to Pfam domains, for PAAH-PAAJ and XDHC-XDHA.

Statistics of the concatenated alignment (CA)

Let us consider a CA of paired AB sequences. At each site i ∈ {1, .., L}, where L is the length of an AB sequence, a given concatenated sequence can feature any of 21 amino acid states α.

The raw empirical frequencies f i (α) and f ij (α, β), obtained by counting the sequences where given residues occur at given sites and dividing by the number M of sequences in the CA, are subject to sampling bias, due to phylogeny and to the choice of species that are sequenced [START_REF] Marks | Protein 3D structure computed from evolutionary sequence variation[END_REF][START_REF] Morcos | Direct-coupling analysis of residue coevolution captures native contacts across many protein families[END_REF]. Hence, we use a standard correction that re-weights "neighboring" concatenated sequences with mean Hamming distance per site < θ. In practice θ = 0.15 was found to yield the best results, but the dependence of performance on θ is weak (see Fig. S6A). The weight associated to a given concatenated sequence a is 1/m a , where m a is the number of neighbors of a within the threshold (including the sequence considered) [START_REF] Marks | Protein 3D structure computed from evolutionary sequence variation[END_REF][START_REF] Morcos | Direct-coupling analysis of residue coevolution captures native contacts across many protein families[END_REF][START_REF] Procaccini | Dissecting the specificity of proteinprotein interaction in bacterial two-component signaling: orphans and crosstalks[END_REF]. This allows one to define an effective total number of sequences, M eff = M a=1 1/m a . We also introduce pseudocounts via a parameter Λ [START_REF] Weigt | Identification of direct residue contacts in protein-protein interaction by message passing[END_REF][START_REF] Marks | Protein 3D structure computed from evolutionary sequence variation[END_REF][START_REF] Morcos | Direct-coupling analysis of residue coevolution captures native contacts across many protein families[END_REF][START_REF] Procaccini | Dissecting the specificity of proteinprotein interaction in bacterial two-component signaling: orphans and crosstalks[END_REF] to avoid issues due e.g. to amino-acid pairs that never appear. Indeed, those can yield mathematical difficulties, such as diverging PMI estimates. Note that pseudocounts are widely used in DCA too [START_REF] Morcos | Direct-coupling analysis of residue coevolution captures native contacts across many protein families[END_REF]. The corrected one-body frequencies fi read fi

(α) = Λ q + (1 -Λ)f i (α) , (10)
where q = 21 is the number of possible states per site. Similarly, the corrected two-body frequencies fij read

fij (α, β) = Λ q 2 + (1 -Λ)f ij (α, β) if i = j , (11) fii
(α, β) = Λ q δ αβ + (1 -Λ)f ii (α, β) = fi (α)δ αβ , (12)
where δ αβ = 1 if α = β and 0 otherwise. We found that Λ = 0.15 yields the best performance of the MI-IPA, but the dependence of performance on Λ is weak (see Fig. S6B). Note that this value is lower than the typical value used in DCA (Λ = 0.5) [START_REF] Marks | Protein 3D structure computed from evolutionary sequence variation[END_REF][START_REF] Morcos | Direct-coupling analysis of residue coevolution captures native contacts across many protein families[END_REF]. Note that our demonstration that the average pairing score S X of an assignment tends to the sum of MIs of inter-protein residue pairs (see Eqs. 4 and 5) did not include proximity weightings and pseudocounts. In practice, at each iteration of our algorithm, assignments are made within each species separately, a scale at which it is convenient to just use the sum of pairing scores S AB (see Eq. 3), and find the assignment that maximizes it. This means that the convergence to the sum of MIs is approximate when using proximity weightings and pseudocounts to calculate the frequencies in the CA and estimate PMIs.

Note also that empirical frequencies without weightings or pseudocounts were employed to study how well the MI-IPA maximizes MI (Figs. 3 andS2). The MI of each inter-protein pair of amino acids was evaluated at the final iteration of the MI-IPA, for three different pairs of protein families, including two interacting ones (BASS-BASR homologs and MALG-MALK homologs) and one with no known interaction (BASR-MALK homologs). These MI scores were averaged over 500 MI-IPA replicates with N increment = 50 that differ in their initial random pairings, and then ranked by decreasing value. Datasets used are the same as in Fig. 5: all alignments include ∼ 5000 protein pairs, with average number of pairs per species m p ≈ 5.

Figure S4: Prediction of inter-protein residue contacts. Fraction among the top interprotein residue pairs ranked by MI or DCA scores that are in contact in the experimentallydetermined three-dimensional structure of complexes. Inter-protein residue pairs were scored both by MI (see Eq. 2) and by the Frobenius norm of the coupling strengths from DCA [START_REF] Bitbol | Inferring interaction partners from protein sequences[END_REF], and then ranked using these scores. (A) HK-RR dataset; contacts were determined from the experimental structure with PDB identifier 3DGE [START_REF] Casino | Structural insight into partner specificity and phosphoryl transfer in two-component signal transduction[END_REF]. (B) Dataset of homologs of MALG-MALK; contacts were determined from the experimental structure with PDB identifier 3RLF [START_REF] Oldham | Snapshots of the maltose transporter during ATP hydrolysis[END_REF]. The alignments of ∼ 5000 cognate protein pairs employed in the rest of the paper to test the MI-IPA (see Methods) were used. A generous threshold of 8 Å on the minimum distance between two amino acids was used to define contacts. In both panels, the chance expectation for finding contacts (i.e. the overall fraction of inter-protein residue pairs that are in contact) is shown (dashed lines). Final TP fraction obtained versus N increment for the standard HK-RR dataset, starting from random within-species HK-RR pairs. For both the MI-IPA and the DCA-IPA, two variants of the protein pair assignment strategy are presented. Once pairing scores are computed, pairs are assigned within each species either using the Hungarian algorithm, which maximizes the sum of scores in the species, or using a greedy algorithm, that favors individual pairs with top scores (see Methods). In all the rest of the paper, the Hungarian algorithm is used for the MI-IPA, and the greedy one is used for the DCA-IPA, as in [START_REF] Bitbol | Inferring interaction partners from protein sequences[END_REF]. All results are averaged over 50 replicates employing different initial random pairings, and the dashed line represents the average TP fraction obtained for random within-species HK-RR pairings. MI and DCA scores trained on the full alignment of cognate HK-RR pairs (excluding orphans, see Methods) were used to score all possible HK-RR pairs, including orphans, in species where orphan and crosstalk partners have been reported in the literature. For each HK considered, RRs were ranked using the MI and DCA scores. The table indicates the species, the HK and RR names, the literature references where partnership was reported, the total number of RRs in the species considered, and the rank among all these RRs of the RR considered, both using MI scores and using DCA scores. Finally, we report a comparison of the DCA ranking with the previous one in [START_REF] Procaccini | Dissecting the specificity of proteinprotein interaction in bacterial two-component signaling: orphans and crosstalks[END_REF]: does our ranking give a better (higher) or a worse (lower) rank to the known partner than this previous ranking? Special mentions in the table: (pair): the cognate pair was included, in order to compare it with the crosstalk partner; Worse*: 2 ranks lower (including an inversion between PleD and DivK, which had close scores in [START_REF] Procaccini | Dissecting the specificity of proteinprotein interaction in bacterial two-component signaling: orphans and crosstalks[END_REF]); Worse**: 1 rank lower (including an inversion between PleD and DivK, which had close scores in [START_REF] Procaccini | Dissecting the specificity of proteinprotein interaction in bacterial two-component signaling: orphans and crosstalks[END_REF]).

Species

Figure 1 :

 1 Figure 1: Performance of the MI-IPA for different training set sizes N start . (A) Increase of the TP fraction with the number of iterations during the MI-IPA. (B) Initial and final TP fractions (at the first and last iteration) versus N start for the MI-IPA and for the DCA-IPA of [21].In both panels, the standard dataset of HK-RRs is used, and the CA includes N increment = 6 additional pairs at each iteration. All results are averaged over 50 replicates that differ by the random choice of HK-RR pairs in the training set. Dashed lines represent the average TP fraction obtained for random within-species HK-RR pairings.

Figure 2 :

 2 Figure 2: Performance of the MI-IPA without a training set. (A) Final TP fraction obtained by the MI-IPA and the DCA-IPA versus N increment for the standard HK-RR dataset.At the first iteration, the CA is made of random within-species HK-RR pairs. At each subsequent iteration n > 1, the CA includes the (n -1)N increment top predicted pairs. All results are averaged over 50 replicates employing different initial random pairings. The MI-IPA was also run on an alignment where each column is randomly scrambled. The dashed line represents the average TP fraction obtained for random within-species HK-RR pairings. (B) Final TP fraction obtained by the MI-IPA and the DCA-IPA versus the total number M of HK-RR pairs in the dataset, starting from random pairings. For each M , except that corresponding to the full dataset, datasets are constructed by picking species randomly from the full dataset, and results are averaged over multiple different such alignments (from 50 up to 500 for small M ). For the full dataset (largest M ), averaging is done on 50 different initial random within-species pairings. All results in (B) are obtained in the small-N increment limit.

Figure 3 :

 3 Figure 3: Increase of pairwise MI obtained by the MI-IPA. (A) Estimates of the pairwise MI (sum of MIs of all inter-protein residue pairs) are shown versus the total number M of sequences in the dataset. The different curves correspond to the pairwise MI of the actual set of correctly paired HK-RRs, of the initial random within-species assignment used to initialize the MI-IPA, and of the final assignment predicted by the MI-IPA, as well as of an alignment where each column was scrambled. Both axes have a logarithmic scale. The slope -1, expected from leading-order finite-size effects, is indicated by the dashed line. For each M , HK-RR datasets are constructed by picking species randomly from the full dataset, and results are averaged over 50 different such alignments. (B) Excess pairwise MI of the actual set of correctly paired HK-RRs, compared to the initial random within-species pair assignments and to the final assignment predicted by the MI-IPA. The MI-IPA successfully reduces this excess pairwise MI, thus approaching the pairwise MI of the actual alignment. Same data as in (A).

Figure 4 :

 4 Figure 4: Performance of the MI-IPA for different protein pairs. Final TP fraction obtained without a training set by the MI-IPA and the DCA-IPA versus N increment . (A,B,C) Pairs involved in ABC transporters; datasets of ∼ 5000 pairs extracted from larger full datasets.(D,E) Pairs involved in enzymatic complexes; full datasets of ∼ 2000 pairs. In each case, the mean number m p of pairs per species is indicated, and dashed lines represent the average TP fraction obtained for random within-species pairings. All results are averaged over 50 replicates that differ in their initial random within-species pairings.

Figure S2 :

 S2 Figure S2: Increase of pairwise MI obtained by the MI-IPA: finite-size effect corrections.In both panels, similar results as those in Fig.3Bare reported, but two different finite-size effect corrections to the MI estimates are implemented. Excess pairwise MI of the actual set of correctly paired HK-RRs, compared to the initial random within-species pair assignments and to the final assignment predicted by the MI-IPA, are plotted versus the total number M of sequences in the dataset. (A) Finite-size correction using subsampling of each dataset and fitting versus 1/M (see Results). Third-degree polynomial fitting was employed. (B) NSB correction[START_REF] Nemenman | Entropy and inference, revisited[END_REF][START_REF] Nemenman | Entropy and information in neural spike trains: progress on the sampling problem[END_REF] (see Results). In both panels, naive estimates (see Fig.3B) are also shown for comparison. For each M , HK-RR datasets are constructed by picking species randomly from the full dataset, and results are averaged over 50 different such alignments.

Figure S3 :

 S3 Figure S3: Alternative signature of protein-protein interactions.The MI of each inter-protein pair of amino acids was evaluated at the final iteration of the MI-IPA, for three different pairs of protein families, including two interacting ones (BASS-BASR homologs and MALG-MALK homologs) and one with no known interaction (BASR-MALK homologs). These MI scores were averaged over 500 MI-IPA replicates with N increment = 50 that differ in their initial random pairings, and then ranked by decreasing value. Datasets used are the same as in Fig.5: all alignments include ∼ 5000 protein pairs, with average number of pairs per species m p ≈ 5.

Figure S5 :

 S5 FigureS5: IPA assignment variants. Final TP fraction obtained versus N increment for the standard HK-RR dataset, starting from random within-species HK-RR pairs. For both the MI-IPA and the DCA-IPA, two variants of the protein pair assignment strategy are presented. Once pairing scores are computed, pairs are assigned within each species either using the Hungarian algorithm, which maximizes the sum of scores in the species, or using a greedy algorithm, that favors individual pairs with top scores (see Methods). In all the rest of the paper, the Hungarian algorithm is used for the MI-IPA, and the greedy one is used for the DCA-IPA, as in[START_REF] Bitbol | Inferring interaction partners from protein sequences[END_REF]. All results are averaged over 50 replicates employing different initial random pairings, and the dashed line represents the average TP fraction obtained for random within-species HK-RR pairings.

Figure S6 :

 S6 Figure S6: Impact of the weighting of similar sequences and of the pseudocount.(A) Final TP fraction obtained for the standard HK-RR dataset, starting from random withinspecies HK-RR pairs, with N increment = 6, for various values of the threshold θ of mean Hamming distance per site below which two sequences are considered as neighbors and weighted accordingly (see Methods). Pseudocount Λ = 0.15 was used. (B) Similar graph, but varying the pseudocount Λ (see Methods), at fixed θ = 0.15. In both panels, all results are averaged over 50 replicates employing different initial random pairings.

Table S1 :

 S1 Orphan partnership and crosstalk prediction in HK-RRs by MI and DCA.

		HK usual	HK ordered lo-	RR usual	RR ordered lo-	Ref.	Number	MI	DCA	Comparison
		name	cus name	name	cus name		of RRs	rank	rank	with [14]
	Thermotoga maritima		TM0853		TM0468	[59]	11	2	2
	Synechococcus elongatus	SasA/CikA	Synpcc7942 2114	SrrA	Synpcc7942 2416	[60]	21	2	2
	Synechococcus elongatus	SasA/CikA	Synpcc7942 2114	RpaB	Synpcc7942 1453	[60]	21	6	4
	Synechococcus elongatus	SasA/CikA	Synpcc7942 2114	RpaA	Synpcc7942 0095	[60, 61]	21	7	5
	Synechococcus elongatus		Synpcc7942 0453	Ycf29	Synpcc7942 1860	[60]	21	2	12
	Synechococcus elongatus		Synpcc7942 2242	SrrB	Synpcc7942 0556	[60]	21	8	5
	Synechococcus elongatus	NblS	Synpcc7942 0924	RpaB	Synpcc7942 1453	[60]	21	1	1
	Synechococcus elongatus	NblS	Synpcc7942 0924	SrrA	Synpcc7942 2416	[60]	21	3	4
	Synechocystis	Hik2	slr1147	Rre1	slr1783	[60, 62]	41	2	7
	Myxococcus xanthus	CrdS	MXAN 5184	CrdA	MXAN 5153	[63]	138	1	1
	Myxococcus xanthus	SasS	MXAN 1249	SasR	MXAN 1245	[63]	138	9	3
	Caulobacter crescentus	DivJ	CC 1063	DivK	CC 2463	[64]	46	8	1	Same
	Caulobacter crescentus	DivJ	CC 1063	PleD	CC 2462	[64]	46	13	2	Same
	Caulobacter crescentus	PleC	CC 2482	PleD	CC 2462	[64]	46	2	1	Better
	Caulobacter crescentus	PleC	CC 2482	DivK	CC 2463	[64]	46	3	4	Worse*
	Caulobacter crescentus	CenK	CC 0530	CenR	CC 3743	[14, 64]	46	17	7	Better
	Caulobacter crescentus	DivL	CC 3484	DivK	CC 2463	[14]	46	16	24	Worse
	Caulobacter crescentus		CC 1062	DivK	CC 2463	[14]	46	14	2	Worse**
	Clostridium acetobutylicum		CA C0903	Spo0A	CA C2071	[65]	42	22	12
	Clostridium acetobutylicum		CA C3319	Spo0A	CA C2071	[65]	42	22	12
	Bacillus subtilis	KinA	BSU13990	Spo0F	BSU37130	[14]	35	1	1	Same
	Bacillus subtilis	KinB	BSU31450	Spo0F	BSU37130	[14]	35	6	3	Better
	Bacillus subtilis	KinC	BSU14490	Spo0F	BSU37130	[14]	35	6	1	Same
	Bacillus subtilis	KinD	BSU13660	Spo0F	BSU37130	[14]	35	1	1	Same
	Bacillus subtilis	KinE	BSU13530	Spo0F	BSU37130	[14]	35	2	1	Same
	Streptococcus mutans (pair)	VicK	Smu 1516	VicR	Smu 1517	[66]	15	1	1
	Streptococcus mutans	VicK	Smu 1516	GcrR	Smu 1924	[66]	15	2	2

Acknowledgments

A.-F. B. thanks Ned S. Wingreen for illuminating discussions. A.-F. B. acknowledges the Aspen Center for Physics, which is supported by NSF Grant PHY-1607611, and more particularly the participants of the 2017 Aspen Working Group on "Covariance Analysis of Proteins" for stimulating discussions.

Supplementary figures and tables