Unveiling Contacts within Macromolecular Assemblies by Solving Minimum Weight Connectivity Inference (MWC) Problems*

Consider a set of oligomers listing the subunits involved in subcomplexes of a macromolecular assembly, obtained e.g. using native mass spectrometry or affinity purification. Given these oligomers, connectivity inference (CI) consists of finding the most plausible contacts between these subunits, and minimum connectivity inference (MCI) is the variant consisting of finding a set of contacts of smallest cardinality. MCI problems avoid speculating on the total number of contacts but yield a subset of all contacts and do not allow exploiting a priori information on the likelihood of individual contacts. In this context, we present two novel algorithms, MILP-W and MILP-WB. The former solves the minimum weight connectivity inference (MWCI), an optimization problem whose criterion mixes the number of contacts and their likelihood. The latter uses the former in a bootstrap fashion to improve the sensitivity and the specificity of solution sets. Experiments on three systems (yeast exosome, yeast proteasome lid, human eIF3), for which reference contacts are known (crystal structure, cryo electron microscopy, cross-linking), show that our algorithms predict contacts with high specificity and sensitivity, yielding a very significant improvement over previous work, typically a twofold increase in sensitivity. The software accompanying this paper is made available and should prove of ubiquitous interest whenever connectivity inference from oligomers is faced.


CONNECTIVITY INFERENCE FROM SETS OF OLIGOMERS
Structural Inference from Oligomers and Contacts-Unraveling the function of macromolecules and macromolecular machines requires atomic level data, both in their static and dynamic dimensions, the latter coding for thermodynamic and kinetic properties (1). However, obtaining even static snapshots of large systems remains a tour de force, so that alternative methods are being developed, based in particular on reconstruction by data integration (RDI) 1 , a strategy aiming at producing models of assemblies using complementary experimental data (2). In its full generality, RDI accommodates both structural and purely combinatorial data (3). The former typically consists of crystallographic (high resolution) structures, electron microscopy maps, and NMR models. The latter comprise information on the composition and copy numbers of subunits, as well as pairwise contacts. Given a large assembly, information on oligomers (i.e. subcomplexes of the assembly) can be obtained by methods such as tandem affinity purification (4) or native mass spectrometry (5,6), and such oligomers can be complemented by information on pairwise protein -protein interactions (7,8). More specifically, oligomers of varying size can be obtained under various experimental conditions. While stringent conditions (e.g. low pH) result in complete dissociation of the assembly so that the individual molecules are identified, less stringent conditions result in the disruption of the assembly into multiple overlapping oligomers. Assembled together, such oligomers can be used to infer contacts within the assembly (5). In the context of RDI, the models stemming from such analysis do not, in general, achieve atomic resolution. They can, however, be used to bridge the gap to atomic level models of subsystems of the assembly under scrutiny (9,10).
Unweighted and Weighted Connectivity Inference: MCI and MWCI-Consider a macromolecular assembly consisting of subunits (typically proteins or nucleic acids). Assume that these subunits are known but that the pairwise contacts between them are unknown or only partially known-in this latter case, the presence or the likelihood of selected contacts is known. Connectivity inference (CI) is the problem concerned with the elucidation of contacts between these subunits, as it ideally aims at producing one contact for each pair of subunits sharing an interface in the assembly (Fig. 1). Note that mathematically, the subunits may be seen as the nodes of a graph whose edges are defined by the contacts. Thus, in the sequel, we use contacts and edges interchangeably.
To address CI, let an oligomer formula be a list of subunits defining a connected component within the assembly. That is, an oligomer formula is the description of the composition of the oligomer, giving the number of copies of each molecule. We define a connectivity inference specification (specification for short) as a list of oligomers. The solution of a CI problem consists of a set of contacts, denoted S in the sequel. This set is called a valid edge set or a solution provided that for each oligomer and also for the whole complex: Restricting the edges from S to the vertices of an oligomer formula yields a connected graph (Fig. 1). In defining valid solutions, two critical questions arise: How many contacts should one seek, and should all edges be treated on an equal footing?
On the Number of Contacts-In the absence of a priori knowledge on the likelihood of individual contacts, a solution is naturally assessed by its number of contacts. Mastering this size is nontrivial since the number of interfaces between subunits in the assembly is unknown. On the one hand, the trivial solution involving all possible edges is uninteresting since it is likely to contain a large number of false positives. On the other hand, one may solve the minimum connectivity inference problem (MCI), namely the variant of CI minimizing the number of contacts used. To do so, observe that minimally connecting an oligomer merely requires a tree, that is a graph whose number of edges is the number of vertices (subunits) minus one. Thus, solving MCI consists of choosing for each oligomer the tree yielding the solution of minimum size. Yet, in doing so, one is likely to generate false negatives. Given these two extremes, one goal of this work is to optimize the number of contacts reported so as to maximize the number of true positives and true negatives-a goal that will be achieved using so-called consensus edges and a bootstrapping strategy.
A Priori Knowledge on Contacts-In a number of cases, the likelihood of a given contact may be known. On the experimental side, various assays have been developed to check whether two proteins interact, including yeast-two-hybrid, mammalian protein-protein interaction trap, luminescencebased mammalian interactome, yellow fluorescent protein complementation assay, co-immuno precipitation, etc (7,8). But information obtained must be used with care for several reasons, notably because expression systems force promiscuity between proteins that may otherwise be located in different cellular compartments and also because affinity purification typically involves concentration beyond physiological levels. On the in-silico side, various interactions attributes can be used, such as gene expressions patterns (proteins with identical patterns are more likely to interact), domain interaction data (a known interaction between two domains hints at an interaction between proteins containing these domains), common neighbors in protein-protein interaction networks, or bibliographical data (number of publications providing evidence for a particular interaction). Here again, these pieces of information have a number of caveats. In particular, structural data from crystallography or mass spectrometry yield a bias toward stable interactions at the detriment of transient ones. For these reasons, strategies computing confidence scores usually resort to machine learning tools trained on the aforementioned data (11) and also (12).
In any case, being able to accommodate such previous knowledge is precisely another goal of the algorithms developed in this work. Given an assembly whose subunits are known but pairwise contacts are not and for which the composition of a number of oligomers in terms of subunits is also known, the problem consists of inferring contacts between subunits. We consider a toy example involving five proteins and three oligomers (three trimers), as seen on panel (A). As additional information, one may enforce and/or forbid contacts, and one may also weight contacts, depending on their likelihood. To connect each oligomer using as few edges as possible, two edges must be chosen, out of three possible (panels (B, C, D)). The minimum connectivity inference consists of finding the overall smallest number of edges such that each oligomer gets connected. Panel (E) shows a solution with four edges (bold edges). Note that these four edges from a subset of all pairwise contacts.
Algorithms-From a computer science perspective, solving CI problems is a hard task (supplemental section 11). Two algorithms targeting such problems have been developed so far. The first one is a two-stage heuristic method (13). First, random graphs meeting the connectivity constraint are generated by incrementally adding random edges. Second, a genetic algorithm is used to reduce the number of edges and also boost their diversity. Once the average size of the graphs stabilizes, the pool of graphs is analyzed to spot highly conserved edges.
The second one is our method solving MCI problems, based on a mixed integer linear program (14). On the one hand, this work delineates the combinatorial hardness of the CI problem and offers two algorithms. Of the particular interest is MILP since it delivers all optimal solutions of a given MCI problem. (As discussed in the supplemental section 11, this algorithm may require exponential time for hard instances, even though this behavior was not observed for the cases processed.) On the other hand, when assessed against contacts seen in crystal structures, the solutions of MILP suffer from two limitations. First, in all solutions, few false negatives are observed at the expenses of selected false positives. On the other hand, since all edges have a unit weight, one cannot favor or penalize some of them.
In this context, this paper makes two improvements. First, we introduce the minimum weight connectivity inference problem (MWCI), which allows computing optimal solutions incorporating a priori knowledge on the likelihood of edges. Second, we present algorithm MILP-W to solve MWCI problems, an algorithm aiming at maximizing the sensitivity and specificity of the set of contacts reported.

MINIMUM WEIGHT CONNECTIVITY INFERENCE: MATHEMATICAL MODEL
Oligomers and Pools of Edges-In solving CI problems, a valid edge set consists of edges such that each of them involves two subunits belonging to at least one oligomer. More precisely, consider an oligomer O i . This oligomer defines a pool of candidate edges equal to all pairs of subunits found in O i . Likewise, the pool of candidate edges Pool E ( ) defined by a set of oligomers is obtained by taking the union of the pools defined by the individual oligomers. Note that one can also consider a restricted set of oligomers involving the oligomers whose size is bounded by an integer s, denoted Յs , the corresponding pool of candidate edges being denoted Pool E ( Յs ). The rationale for using small oligomers is that they favor local contacts. (Note that the extreme case is that of a dimer since the contact seen in a dimer must belong to every solution.) Note also that one can edit a pool of edges to enforce or forbid a given edge in all solutions. For example, if a cryo-electron microscopy map of the assembly is known and two proteins have been located far apart in the map, one can forbid the corresponding contact even though the two proteins appear in a common oligomer.
We now present two ways to solve CI problems. Unweighted Case-In the unweighted case, each edge from the pool is assigned a unit weight so that the weight of a solution is the number of its edges. The corresponding optimization problem is called MCI, and an algorithm solving it, MILP, has been proposed in (14).
Weighted Case-In the weighted case, each candidate edge e from the pool Pool E ( ) is assigned a weight w(e), namely a real number in range [0,1]. This number encodes the likelihood for the edge to be a true contact. Taking G ϭ 1/2 as a baseline (i.e. no a priori on this contact), a value F Ͼ G is meant to favor the inclusion of this edge in solutions, while a value U Ͻ G is meant to penalize this edge.
Unifying the Unweighted and Weighted Cases: MWCI Problems-Depending on how much information is available on candidate contacts, one may wish to stress the number of contacts in a solution, or their total weight. Both options can actually be handled at once by interpolating between the previous two problems. Using a real number ␦ ʦ [0,1], we define a functional mixing the number of edges and their weights, this latter one being favored for values beyond the threshold 1/2. That is, we define the cost of a solution S using two terms respectively corresponding to the number of edges and their weights: Eq. (1) corresponds to the objective of the optimization problem denoted MWCI in the sequel. The following comments are in order: • In using ␣ ϭ 1, which is the strategy used by algorithm MILP (14), the weights play no role, and the interchangeability of edges favors the exploration of a large pool of solutions. • The situation is reversed for small values of ␣. In particular, the conjunction ␣ ϭ 1 and different weights for all edges typically yields a small number of solutions since ties between solutions are broken by the weights. • A null weight does not prevent a given edge to appear in solutions. To forbid an edge, one should edit the pool of candidate edges, as explained above i.e. remove this edge from the pool.
Remark 1-Assume that each edge has a default weight d instead of 1/2. Eq. (1) is a particular case of the following:  (14) and allows enumerating all optimal solutions with respect to the criterion of Eq. (1). The algorithm solves a mixed integer linear program, using constraints imposing the connectivity constraints inherent to all oligomers. Candidate edges are represented by binary variables taking the value 1 when edges belong to a specific solution (14) and 0 otherwise.
More precisely, algorithm MILP-W iteratively generates all optimal solutions and adds at each iteration extra constraints preventing finding the same solution twice. To this end, the method starts with a first resolution of the problem to get an optimal solution, if any. This solution defines a set of edges and the associated value OPT for the criterion of Eq. (1). To check whether another solution matching OPT exists, a new constraint preventing the concomitant selection of all edges from the first solution is added. More formally, the sum of the binary variables associated with the solution just produced is forced to be strictly less than the number of edges in solutions seen so far. The resolution is launched again, and the criterion value is compared with OPT. This process is iterated until the value of the solution exceeds OPT. As noticed earlier, when ␣ ϭ 1, algorithm MILP-W matches algorithm MILP. Therefore, for the sake of clarity, the solution set, consensus solutions, and the associated edge sets are respectively denoted MILP , MILP-W cons , MILP , and MILP-W cons . These notations are summarized in Table I.
To further assess the quality of the solution set (ϭ MILP , MILP-W ), assume that a reference set of contacts E Ref is known. The ideal situation is that where a high-resolution crystal structure is known since then all pairwise contacts can be inferred (15). This reference set together with the pool Pool E ( ) define positive (P), negative (N), and missed contacts (M) (Fig. 2). From these groups, one further classifies the edges of a predicted solution in set into four categories, namely true positive (TP), false positive (FP), true negative (TN), and false negative (FN).
Positives (P) and negatives (N) decompose as P ϭ TP ϩ FN, and NϭTNϩFPss, from which one defines the sensitivity ROC spec and the specificity ROC spec as follows: Note that specificity requires the set N to be nonempty, which may not be the case if Pool E ( )ʚ E Ref .
We also combine the previous values to define the following coverage score, which favors true positives, penalizes false positives and false negatives, and scales the results with respect to the total number of reference contacts (since P might be included into E Ref if the pool size is too small): Note that the maximum value is one and that the coverage score may be negative.

Algorithm MILP-W B -
The focus on consensus edges is quite natural since these may prosaically be seen as the backbone of the connectivity in the assembly. However, alternative edges of significant importance may exist too. To unveil such edges, we preclude one or more consensus edges, so as to trigger a rewiring of the connectivity of solutions, and check which novel consensus edges appear along the way. Implementing this strategy requires two precautions, namely: (i) edges corresponding to dimers must be kept for a solution to be valid and (ii) hindering too many edges may yield a connectivity inference problem without any solution.
More precisely, we start precluding the consensus contacts i.e. the initial consensus contacts MILP-W cons minus the dimers, one at a time from the pool of contacts to be explored. We subsequently report the union of consensus contacts (including the initial consensus contacts which we began with) yielded after all the MILP-W runs. The process can be iterated by precluding two or more contacts at a time. This strategy triggers rewiring of the system to a greater extent, at the risk of inducing more false positives. (See pseudo code in the supplemental Section 11.) When reporting results obtained with this algorithm, the following conventions are used.

MATERIAL: TEST SYSTEMS
We test the performance of the algorithms MILP-W and MILP-W B on the following three systems for which reference contacts for validation are available either coming from crystal structure or from various biophysical experiments such as cryo-EM based reconstruction, cross-linking, and MS/MS dimers. See supplemental Section 8 for the input to the algorithms and the supplemental Section 9 for the reference contacts.
4.1 Yeast Exosome-The exosome involves 10 protein types (Figs. 3 and 4), and 19 oligomers have been reported (13), ranging in size from two to nine (Table II and supplemental Section 8). Note that originally 21 oligomers are reported in (13), including a trivial case of having all 10 protein types and one other oligomer of size 8 that has a duplicate, leaving behind 19 distinct oligomers.
Oligomers up to size five are required to encompass nine out of 10 proteins-the protein Csl4 is present in size nine oligomers only. In terms of contacts, classical interfaces modeling tools (15) applied to the crystal structure yield 26 contacts amid the 10 proteins, and 20 contacts in the assembly depleted of Csl4 (Fig. 4).
The status of Csl4 is interesting since, as discussed in Section 2, local contacts are favored by small oligomers. In the sequel, we therefore consider two settings, namely the full exosome, and the exosome without Csl4. In the former case, all oligomers define a pool MILP-W cons of 45 candidate edges; in the latter, the pool Pool E ͑8͒ contains 36 candidate edges.
4.2 Yeast 19S Proteasome Lid-Proteasomes are protein assemblies involved in the elimination of damaged or misfolded proteins, and the degradation of short-lived regulatory proteins. The most common form of proteasome is the 26S, which involves two filtering caps (the 19S), each cap involving a peripheral lid, composed of nine distinct protein types each with unit stoichiometry (Fig. 9).
Series of overlapping oligomers were formed by mass spectrometry (MS), tandem MS, and cross-linking using BS3. In total, 14 complexes were obtained out of which eight came from MS, MS/MS and six came from cross-linking experiments (16) (Table III and supplemental Section 8). 4.3 Human eIF3-Eukaryotic initiation factors (eIF) are proteins involved in the initiation phase of the eukaryotic translation. They form a complex with the 40S ribosomal subunit, initiating the ribosomal scanning of mRNA. Among them, human eIF3 consists of 13 different protein types each with unit stoichiometry (Fig. 12). The eIF3 complex in this text refers to the human eIF3 unless otherwise stated. A total of 27 complexes were generated from the assembly by manipulating the ionic strength of the solution and using tandem mass spectrometry (17) (Table IV and supplemental Section 8). The subunit eIF3j is labile, and since none of 27 subcomplexes comprise this subunit, we exclude it from the list of protein types, leaving behind 12 protein types.

RESULTS
We first provide results for MILP-W B with ␣ϭ1, namely when all edges have the same weight. In a second step, we illustrate the benefits of using weights.
5.1 Algorithm MILP-W B -As explained in Section 3.3, algorithm MILP-W B works by accumulating consensus contacts from MILP-W and those due to local rewiring as a result of precluding the initial consensus contacts one (or more) at a time. Consequently, our analysis focuses on the sensitivity, specificity and coverage statistics introduced in Section 3.2 in four settings: • (C-1) The statistics for the edge set MILP , which serve as a baseline. The following consistent observations can be made for the three systems.
• In comparing (C-1) against (C-2), the sensitivity decreases since consensus solutions have fewer edges. On the other hand, the specificity increases, indicating a FIG. 4. Yeast exosome with and without Csl4: Contacts between subunits. Each edge corresponds to an interface between two subunits. The two numbers decorating an edge, respectively, refer to the number of atoms involved at that interface, and to the number of patches (connected components) of the interface. Interfaces were computed with the program INTERVOR, which implements the Voronoi model from (15). Note that a given subunit makes from three (e.g. Rrp40) to seven (e.g. Rrp45) interfaces.  Table  8 for the detailed statistics. large number of true negatives or equivalently a small number of false positives-an observation in line with the high scores of edges in consensus solutions.
• In comparing (C-3) against (C-4), the sensitivity increases, while the specificity decreases. The variations observed for sensitivity and specificity actually depend on the number of contacts precluded in the bootstrap procedure. Indeed, precluding more contacts triggers more rewiring, which in turns yields a larger set of true positive edges (increased sensitivity), at the expense of more false positives (decreased specificity). • Finally, for algorithm MILP-W B , one observes that a small number of iterations, typically in the range 1,…,3, is favorable to high coverage. This is because of the aforementioned counterbalance between sensitivity and specificity. In particular, when the bootstrap procedure halts, all contacts from the pool have been used, which entails a large number of true positives-high sensitivity-but also a large number of false positives-low specificity. Thus, the user may choose the risk level (in terms of false positives) he/she is willing to accept, depending on whether the focus is on sensitivity or specificity.
FIG. 6. Yeast exosome: variation of cumulative sensitivity, specificity and coverage score with iteration index. Note that the iteration index also indicates number of contacts forbidden at a time. See supplemental Table 9 for the detailed statistics.  Table 8 for the detailed statistics. Comparison to Previous Work-The statistics just discussed compare favorably to previous work, obtained in particular with the heuristic network inference algorithm (13). We illustrate this fact with the results produced by MILP-W B after one iteration.
On the yeast exosome with Csl4, the sensitivity of MILP-W B is ϳ1.67 (ϭ 0.77/0.46) times that of network algorithm and Cvg. score increases from -0.08 to 0.35 (lines T3 versus T0 in the supplemental Table 8).
The comparison is not possible for eIF3, though, since the previously published contacts were computed manually using experimental information from various other sources (17). See, however, the supplemental Table 13 for the results using our algorithms.
5.2 Algorithm MILP-W-In this section, we illustrate the role of weights stemming from using Eq. (1) with a 1. One naturally expects a benefice in penalizing an edge that is a negative contact and that is thus predicted as a false positive or a true negative. But an improvement of statistics can also happen by merely favoring positive contacts. As an illustration, we consider the yeast exosome without Csl4 (specifications in Table II), using ␣ϭ0.25. We assign a weight of 0.6 to the following three contacts-(Rrp45, Rrp46), (Rrp40, Rrp46)  Table 11 for the detailed statistics.
FIG. 11. Yeast proteasome lid: variation of cumulative sensitivity, specificity, and coverage score with iteration index. Note that the iteration index also indicates number of contacts forbidden at a time. See supplemental Table 12 for the detailed statistics.  Table 13 for the detailed statistics. and (Rrp41, Rrp42), the remaining contacts having the default weight of 0.5.
Upon moving from the instance without weights (i.e. ␣ϭ1) to the instance with weights (i.e. ␣ϭ0.25), we consider the changes in the sensitivity, specificity, and coverage, namely: ⌬ ϭ ͑⌬ROC sens. ,⌬ROC spec. ,⌬Cvg) (Eq. 6) Consider first MILP . One observes four false positives instead of five (supplemental Fig. 15), improving the assessment tuple by (0, 0.06, 0.05). Thus, while none of the contacts has a weight less than 0.5, the relative value of weights has an incidence on the outcome. Consider now MILP cons . The union of consensus contacts has no false positive, and true positives are increased from 9 to 13, improving the assessment tuple by (0.2, 0.06, 0.45) (supplemental Fig. 16). For MILP-WB , the change in assessment tuple for MILP-WB (first iteration) is (0.05, Ϫ0.12, 0). When more contacts are precluded, the trend seen is similar to earlier cases (supplemental Fig. 17).
The reader is referred to the supplemental Section 12 for a thorough assessment obtained upon varying the weights and the value of ␣.

DISCUSSION AND OUTLOOK
By giving access to a list of overlapping oligomers of a given macromolecular assembly, native mass spectrometry offers the possibility to infer pairwise contacts within that assembly, opening research avenues for systems beyond reach for other structural biology techniques. In this context, our work makes three contributions, based on state-of-the art combinatorial optimization techniques.
First, we introduce the minimum weight connectivity inference problem (MWCI), which generalizes the minimum connectivity inference problem, by introducing weights associated with putative contacts. Second, we develop algorithm MILP-W to solve MWCI problems, taking into account a priori biological knowledge on the likelihood of contacts. Third, we FIG. 14. Human eIF3: variation of cumulative sensitivity, specificity, and coverage score with iteration index. Note that the iteration index also indicates number of contacts forbidden at a time. See supplemental Table 14 for the detailed statistics.   also develop algorithm MILP-W B , a bootstrap strategy aiming at enriching the solutions reported by MILP-W. Algorithm MILP-W B accumulates consensus contacts from MILP-W and those arising due to local rewiring as a result of precluding the initial consensus contacts one (or more) at a time. Our algorithms predict contacts with high specificity and sensitivity, yielding a very significant improvement over previous work, typically a twofold increase in sensitivity. Despite the combinatorial complexity of the problems addressed, all runs of algorithm MILP-W terminated within a handful of seconds for all the cases processed in this work. Calculations with algorithm MILP-W B are more demanding, though, since the runtime depends on the combinatorics of the tuples to be precluded. These algorithms raise a number of opportunities and challenges.
In the context of native mass spectrometry, they offer the possibility to test various parameter sets, in particular regarding the number of contacts and their likelihood, and to compare the solutions obtained. More broadly, the ability to take into account confidence levels on putative edges should be key to incorporate scores currently being designed in proteomics in conjunction with various assays.
In terms of challenges, fully harnessing these algorithms raises difficult questions. On the practical side, one current difficulty is the lack of cases to learn from, namely assemblies for which a significant list of oligomers is known, and a highresolution structure has been obtained. Such cases would be of high interest to tune the balance between the aforementioned two criteria (number of contacts and their likelihood). This would also aid in carrying out an in-depth study of incidence of weights on the solutions obtained from MILP-W runs, given true positives and false positives in the pool of contacts. Unfortunately, mass spectrometry studies are typically attempted on assemblies whose high-resolution structure are unknown and are likely to remain so, at least in the near future. On the theoretical side, outstanding questions remain open. The first one deals with the relationship between the set of oligomers processed and the solutions generated. Ideally, one would like to set up a correspondence between equivalence classes of oligomers yielding identical solutions. The ability to do so, coupled to the understanding of which oligomers are most likely generated, would be of invaluable interest. The second one relates to the generalization of our algorithms to accommodate cases where multiple copies of subunits are present. However, the multiple copies complicate matters significantly so that novel insights are called for not only computing solutions, but also representing them in a parsimonious fashion.
In any case, we anticipate that the implementations of our algorithms will prove its interest for the growing community of biologists using native mass spectrometry.
□ S This article contains supplemental material Tables 8 to 14