CoEVFold suite: user friendly pipelines to visually represent protein coevolution

7 8 Multiple sequence alignment (MSA) data underlies current principles in protein folding and 9 protein-protein interaction prediction, from which large language models (LLMs) in tandem 10 with protein datasets, can predict protein structure. However, what is missing are user-friendly 11 tools that enable researchers to predict and demonstrate coevolution - the principal input 12 which these MSAs infer. Here we present tools to identify and visualize coevolution, through 13 a pipeline (CoEVFold) that uses basic direct coupling algorithms derived from GREMLIN and 14 alignment of sequences from MMSEQs2. The pipeline generates a visual representation of 15 coevolution for a single protein but can also represent coevolution of homomeric or 16 heteromeric protein complexes, as well as coevolution within protein networks. The input for 17 this pipeline can be an amino acid sequence, or user input protein structures from Alphafold 18 their own files or the PDB database. In validation of CoEVFold’s capabilities, and utilising 19 proteins from known prokaryotic and eukaryotic model systems (Bacillus subtilis, Escherichia 20 coli and Saccharomyces cerevisiae) , as well as phage proteins, CoEVFold predicts 21 coevolution between proteins known to interact, proteins known to oligomerise, and 22 coevolution in proteins known to be part of a protein complex. Collectively, these suite of tools, 23 named ‘CoEVFold suite’, have broad applicability and provide a useful toolkit to those 24 interested in dissecting protein-protein interactions and networks. 25 26 Availability: The code is available online at 27 https://colab.research.google.com/drive/1MSSvNTq7KZ4Lr0XTz89vUuK-J3xOTzwS?usp=sharing and 28 Github. https://github.com/MishterBluesky/CoEVFold 29 30 Supplementary information: Supplementary data is available via Figshar e and supplementary 31 materials. 32 33

34 35 Coevolution is the evolution of one entity relative to another. It ensures continued function of complex 36 systems, allows for environmental adaptation, and establishes new functions. With respect to amino 37 acids in proteins, coevolution often occurs between two amino acids to enable folding and 38 stabilisation of the protein. It is generally theorised, and has been shown by the success of learning 39 .CC-BY-ND 4.0 International licenseperpetuity. It is made available under a preprint (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in The copyright holder for thisthis version posted January 26, 2026. ; https://doi.org/10.64898/2026.01.26.701017doi: bioRxiv preprint .CC-BY-ND 4.0 International licenseperpetuity. It is made available under a preprint (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in The copyright holder for thisthis version posted January 26, 2026. ; https://doi.org/10.64898/2026.01.26.701017doi: bioRxiv preprint .CC-BY-ND 4.0 International licenseperpetuity. It is made available under a preprint (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in The copyright holder for thisthis version posted January 26, 2026. ; https://doi.org/10.64898/2026.01.26.701017doi: bioRxiv preprint .CC-BY-ND 4.0 International licenseperpetuity. It is made available under a preprint (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in The copyright holder for thisthis version posted January 26, 2026. ; https://doi.org/10.64898/2026.01.26.701017doi: bioRxiv preprint .CC-BY-ND 4.0 International licenseperpetuity. It is made available under a preprint (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in The copyright holder for thisthis version posted January 26, 2026. ; https://doi.org/10.64898/2026.01.26.701017doi: bioRxiv preprint .CC-BY-ND 4.0 International licenseperpetuity. It is made available under a preprint (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in The copyright holder for thisthis version posted January 26, 2026. ; https://doi.org/10.64898/2026.01.26.701017doi: bioRxiv preprint .CC-BY-ND 4.0 International licenseperpetuity. It is made available under a preprint (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in The copyright holder for thisthis version posted January 26, 2026. ; https://doi.org/10.64898/2026.01.26.701017doi: bioRxiv preprint .CC-BY-ND 4.0 International licenseperpetuity. It is made available under a preprint (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in The copyright holder for thisthis version posted January 26, 2026. ; https://doi.org/10.64898/2026.01.26.701017doi: bioRxiv preprint .CC-BY-ND 4.0 International licenseperpetuity. It is made available under a preprint (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in The copyright holder for thisthis version posted January 26, 2026. ; https://doi.org/10.64898/2026.01.26.701017doi: bioRxiv preprint .CC-BY-ND 4.0 International licenseperpetuity. It is made available under a preprint (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in The copyright holder for thisthis version posted January 26, 2026. ; https://doi.org/10.64898/2026.01.26.701017doi: bioRxiv preprint .CC-BY-ND 4.0 International licenseperpetuity. It is made available under a preprint (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in The copyright holder for thisthis version posted January 26, 2026. ; https://doi.org/10.64898/2026.01.26.701017doi: bioRxiv preprint .CC-BY-ND 4.0 International licenseperpetuity. It is made available under a preprint (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in The copyright holder for thisthis version posted January 26, 2026. ; https://doi.org/10.64898/2026.01.26.701017doi: bioRxiv preprint .CC-BY-ND 4.0 International licenseperpetuity. It is made available under a preprint (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in The copyright holder for thisthis version posted January 26, 2026. ; https://doi.org/10.64898/2026.01.26.701017doi: bioRxiv preprint .CC-BY-ND 4.0 International licenseperpetuity. It is made available under a preprint (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in The copyright holder for thisthis version posted January 26, 2026. ; https://doi.org/10.64898/2026.01.26.701017doi: bioRxiv preprint models, that a highly co-evolved amino acid pair within a protein are commonly in physical proximity 40 (1). Thus, coevolution between amino acids within a protein can effectively reveal its 3D landscape, or 41 rank-based coevolution amino acid contact map (1). In addition, coevolution can be used to discover 42 protein self-interactions (homomeric) to assign oligomeric state, or to discover protein-protein 43 interactions between two or more different proteins (heteromeric). The predictive power of 44 coevolution, something which relies on the genetic history of the organism therefore, is great. In the 45 absence of experimental data, it alongside structural tools could be used to solve or work towards 46 solving problems which as of yet are impossible or arduous. 47 48 The GREMLIN direct coupling analysis algorithm, established by the Baker lab(1,2) is a powerful 49 coevolution predictor. It was used as the one of the bases towards powerful protein structure 50 prediction tools such as PyRosetta(3); However many coevolution tools have lost support in favour of 51 AI structural prediction tools that give visual representation of the protein tertiary structure. One of the 52 major drawbacks of AI-based tools is the missing information on how proteins folds (secondary and 53 tertiary structure) or protein-protein-protein interactions are derived and modelled. This is important 54 information to understand validity and relatedness of protein structures. This is where coevolution can 55 become an important feature for rationalising modelled protein structures, identifying core structural 56 features and identifying evolutionary relationships. 57 58 Coevolution describes the degree to which two molecular elements whether amino acid residues, 59 domains, or entire proteins change in tandem in a predictable way (4,5). Coevolution occurs at the level 60 of species and prophage, between two co-dependent proteins, or two co-dependent amino acids that 61 reinforce a structure. The simplest coevolution to quantify is on the amino acid level; An illustrative 62 example of intramolecular amino acid contact prediction using coevolution is provided in Figure 1A. A 63 sequence alignment of a hypothetical protein from seven species indicates an Arginine 3 to Aspartate 64 10 pairing (Figure 1) such that when Arginine evolves to Glutamate (E), a reciprocal amino acid change 65 takes place to allow for stabilisation and folding (Figure 1B). If this occurs in the same protein across 66 multiple species, coevolution is likely occurring. The ability to predict coevolution between amino acids, 67 (i.e. likelihood of amino acid Y changing, dependent on amino acid X) can be determined as a 68 coevolution ‘score’ that is calculated by a direct coupling algorithm, which discerns the likelihood of 69 reciprocal change : some amino acids pairs may not always change across different species and 70 therefore harbour lower coevolution, whereas others will change , thus producing a highly coevolving 71 pair. This score can be represented in a 2D grid as a contact map (Figure 1C), and these pairings 72 provide information which build a 3D structure of the protein (Figure 1D). Indeed, until recently some 73 proteins structures were originally solved using coevolution , such as RodA , involved in cell wall 74 remodelling in Escherichia coli (4). Coevolution predictions can also be applied to protein -protein 75 interactions, to provide a interaction network. These are powerful insights for coupling protein sequence 76 to biology, as shown for interactions of RodA with cognate partner protein PBP2(6). 77 78 Coevolutionary analysis informs on both intramolecular interactions required for folding monomers (Fig 79 1B), on protomer structure as well as interactions gene networks , and in some cases intermolecular 80 .CC-BY-ND 4.0 International licenseperpetuity. It is made available under a preprint (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in The copyright holder for thisthis version posted January 26, 2026. ; https://doi.org/10.64898/2026.01.26.701017doi: bioRxiv preprint interactions between protomers of a homomeric complex (7). However the differentiation of 81 intermolecular interactions can be difficult to isolate due to protomers often residing in a shared 82 biochemical environment which require parallel changes similar to coevolution coupling (5). Thus, while 83 coevolution can be used as an additional layer of evidence for protein-protein interactions, it is not 84 always representative of a direct intermolecular contact. 85 86 Figure 1. Basic concepts in coevolution 87 (A) Multiple sequence alignment of a protein in seven species, illustrating reciprocal changes in Arginine 88 3 (R- green) with Lysine 10 (K-green) to ‘stabilise’ a protein (B) Contact interface as observed in 2D, 89 if coevolution is strong, to infer a direct physical connection, (C) contact map of coevolutionary scores 90 (D) Hypothetical overall 3D structure built up of many physical connections informed by coevolution. 91 92 Based on the above principles and with this caveat in mind, we have designed a series of python Jupyter 93 notebooks, called ‘ CoEVFold suite’ that form an easy-to-use pipeline to study coevolution and use 94 structure to further inform on this coevolution. In these notebooks, we have allowed users to input their 95 own protein structures from protein folding models, or experimentally determined structures, which can 96 be used as a basis for the plotting coevolution between proteins in homomeric or heteromeric states 97 and establishing degrees of coevolution between protein complexes. Finally, we utilise coevolution to 98 determine the likelihood of protein-protein interactions within proteins acting within the same biological 99 pathway. 100 101 102 103 104 105 A B CD A T R L T L T A K K D G P C C A T R L T L T A K K D G P C C A T R L T L S A K K D G P C C A T E L T L T A K Y D G P C C A T E P T L I A K Y D G P C C A T E L T L T A K Y D G P C C A T S V T L T A K T D G P C C coevolution 3D coevolution annotation C D .CC-BY-ND 4.0 International licenseperpetuity. It is made available under a preprint (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in The copyright holder for thisthis version posted January 26, 2026. ; https://doi.org/10.64898/2026.01.26.701017doi: bioRxiv preprint

106 107 We have produced three tools that exploit coevolution for the prediction of protein-protein interactions: 108 1) CoEVFold - an easy to use coevolution tool; This tool allows users to identify coevolution between 109 two or more distinct proteins of interes t; 2) CoEVFold:4D, a tool based on CoEVFold to identify 110 coevolution between protomers of homooligomers and 3) CoEVMapper, which predicts protein-protein 111 interactions beyond the protein to itself and a single target creating a network, formed from a target set 112 of proteins within a protein-network. 113 114 Identification of coevolution within two or more distinct proteins in a hetero-oligomer: CoEVFold 115 116 Protein complexes underpin many of the catalytic or regulatory functions of cells from across all walks 117 of life. Understanding assembly and interactions as well as the ev olution of such complexes is an 118 essential part in discriminating mechanistic function . To identify coevolution between two or more 119 distinct proteins, CoEVfold first relies on MMSEQ2 (8) to automatically an alignment of amino acids 120 based on similarity to a user protein of interest. It then uses direct coupling analysis learning algorithm 121 (GREMLIN)(1), to identify scores for each residue-residue contact and plots this in a grid similar to the 122 illustrations of Figure 1A,(5,7). The protein X to protein Y intermolecular coevolution is then plotted upon 123 a user input predicted multiprotein complex. 124 (9) 125 Looking at a specific example, of coevolution of bacterial cell envelope protein FtsW with cell division 126 partner penicillin binding protein B (P bpB) ( Figure 2A-C) (6,9,10). PbpB homologues have been 127 experimentally shown to interact with FtsW , through their transmembrane heli ces, a nd AlphaFold 128 predictions are consistent with this experimental data (6,9,10). However, if a user adds in other 129 transmembrane-containing proteins , these will often interact with FtsW at the same site , due to 130 hydrophobic compatibility, thus making pTM scores alone unreliable. Through coevolution mapping with 131 CoEVfold we can generate coevolution scores and map them to the transmembrane region Figure 2C. 132 This provides discrimination between true binding partners and complexes assembled primarily out of 133 charge compatibility. In another example, we plot coevolution of the multiprotein yeast VO type ATPase 134 complex (11). This multiprotein unit must assemble in a very precise manner to enable its dynamic and 135 essential function. The complex shows inter-protein coevolution between chains in its heteromeric state, 136 in particular sharing coevolution with its alpha helical chain encoding genes (Figure 2D-E) such as 137 chain E that forms the central barrel with subunits C, D and to L. This shows distinctive inter-protein 138 coevolution, visible on the coevolution plot (Figure 2D) as well as on the 3D representation (Figure 2E) 139 with interactions indicated by interchain connectors. Each subset of interactions are shown in separate 140 colours. 141 .CC-BY-ND 4.0 International licenseperpetuity. It is made available under a preprint (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in The copyright holder for thisthis version posted January 26, 2026. ; https://doi.org/10.64898/2026.01.26.701017doi: bioRxiv preprint 142 Figure 2. Schematic representation of the CoEVFold pipeline. 143 (A) CoEVFold window, MMSEQ2 alignment and a3m production based on query sequence – B subtilis 144 FtsW-PbpB >50% sequence coverage. Red – Low coverage/conservation, Blue/Green – High 145 coverage/conservation (B) Coevolution of FtsW-PbpB performed using pipeline on transformed 146 alignment. Indicating regions of intermolecular and intramolecular coevolution Dark Blue – coevolution 147 pbpBftsWA B C Coevolution Threshold Filter <12A contact MMSEQ2 alignment D E A BCDE-LMNO high, White – no observed coevolution Z score min=0, max=3 (C) user input structure FtsW-PbpB 148 predicted by Alphafold, inter-protein coevolution has been indicated using red dashed lines = 149 Coevolution at user thresholds (12A, 50% coverage, Chain indicated above performed using 151 pipeline on transformed alignment indicating regions of intermolecular and intramolecular coevolution 152 Dark Blue – coevolution high, White – no observed coevolution Z score min=0, max=3 (E) Query pdb 153 input (PDB ID 7TA0) Cryo-EM structure of Bafilomycin A1 bound to yeast VO V-ATPase Dashed lines 154 represent coevolution links at user thresholds (30A, Z score > 2.5). Each colour represents a unique 155 inter-chain quaternary link subset. 156 157 Identification of coevolution between homo-oligomers: CoEVFold:4D 158 159 Protein oligomers are also important in building our understanding of the functionality of cells. In our 160 second tool: CoEVFold:4D, we look for homo-oligomers. Other groups have also used coevolution to 161 achieve this, however, have not provided tools for other groups to do this (3,4), in our tool, we hope to 162 identify coevolution for single proteins in their quaternary, ‘homomeric’ state. CoEVFold:4D works on 163 the idea of uncovering coevolution that isn’t explained by a monomeric structure. By comparing the 164 ‘coevolutionary 2D contact map’ with physical contact maps from the amino acid backbone of respective 165 user inputs from high confidence monomer structures and multimers. These are coevolutionary contacts 166 which are predicted in the 2D map (Figure 3B) but aren’t present in the monomeric state (Figure 3C) 167 are evaluated for attribution to new inter-protein links between homomeric forms of the protein (Figure 168 3E). Dense packing of otherwise unexplained coevolved residues neatly points to a quaternary 169 interaction in that space, thus a tool able to point to these regions in 3D space would allow resolution 170 of homo-oligomer coevolution (Figure 3G-H). 171 172 This tool is inherently bias and reliant on structural predictions of the monomer. However, it can either 173 be used to confirm an existing predicted multimer or be used to identify regions of contact. We will 174 explore both a guided ‘confirmatory’ quaternary evolution, and an unguided search for unexplained 175 coevolution. There is a caveat here in both modes, coevolution in addition to new interactions could be 176 due to conformational or phylogenetic coevolution noise(11,12). Not all coevolution can be explained 177 by intra-protein and inter-protein interactions (1,5,7). Regions of unexplained coevolution could 178 inadvertently cause users to believe a quaternary interaction is occurring when none is present. These 179 interactions may occur through shared signalling pathways, shared interaction proteins, pH changes, 180 and a conserved local environment, as well as input of aligned sequences. This can be minimised using 181 a combination of coevolution and predicted structures and assemblies. However, these observations 182 are inherently more informed than structural predictions alone. 183 184 Figure 3. Homomeric inter-protein interaction methodology CoEVFold 4D 185 (A) MMSEQ2 alignment and a3m production based on query sequence – B subtilis motor domain 186 SpoIIIE 50%> coverage Red – Low coverage/conservation, Blue/Green – High coverage/conservation 187 (B) Coevolution performed using pipeline on transformed alignment. Dark Blue – coevolution high, 188 White – no observed coevolution Z score min=0, max=3 (C) User input monomer/lower order multimer 189 and higher order multimer, creating contact map (including intermolecular links). (D) Comparison of 190 monomer contact map with coevolution, showing monomer/low order unexplained coevolution (E). 191 (F/G) Refinement of search area to coevolution only present in input higher order multimer, red- 192 multimer explained evolution, purple – monomer explained coevolution, (H) Optional thresholds of 193 multimer only explained contacts are displayed in a prepared structural .pse file.on input multimer. 194 195 MMSEQ2 alignmentCoevolution (GREMLIN)Input monomer Input multimer Monomer Contact map Multimer Contact map Refine coevolution comparison Isolate Multimer explained CoevolutionPlot on input multimer structure A B D Isolate unexplained coevolution Coevolution Threshold Filter <12A contact C E F G H We applied our tools on a set of multimeric proteins, to confirm that coevolution pairs not explained by 196 monomeric structures can be caused by higher order multimerization of the protein. Using methods as 197 previously described in B. subtilis sporulation genes (5). We took an existing known tetramer of E. coli, 198 MutS (13), and compared the monomer coevolution with known crystal structure of the dimer. 199 Unexplained coevolution of MutS beyond the dimer was identified in overlaid contact maps (Figure 4 200 bottom centre). A major cluster of unexplained coevolution was located at the protein C-terminus. In 201 vitro mutagenesis studies have indeed confirmed the C-terminus to be involved in oligomerisation 202 (Figure 4) (14). Residues Aspartate 835 and Arginine 840, located at the C-terminus of MutS reduced 203 higher order multimerization (14), correspondingly these residues displayed a high unexplained 204 coevolution score. 205 206 207 208 Figure 4. Unexplained Coevolution searching yields similar residues to in vitro searches 209 (top left) MutS MMSEQ2 clustered sequences above 50% coverage Red – Low coverage/conservation, 210 Blue/Green – High coverage/conservation, (top centre) coevolution predicted by tool, Dark Blue – 211 coevolution high, White – no observed coevolution Z score min=0, max=3 (top right) user input dimer 212 12Å amino acid contact map black (bottom left) Overlay onto coevolution (bottom centre) Unexplained 213 coevolutionary contacts -Red (bottom right) MutS dimer predicted by AlphaFold. Black box, mutated 214 residues shown by in vitro experiments to be responsible for tetramer formation. 215 216 We then used this technique to re-identify higher order structures of proteins, in a bias fashion, the 217 second mode of CoEVFold:4D. We tested the known multimers; SpoIIIE’s motor domain (B. subtilis 218 inferred by FtsK homologue (14), MlaD (Escherichia coli)(15) and Tequatrovirus T4 phage tail protein 219 4UXG(16). Unexplained coevolution contacts unique to contact maps of the multimers were identified 220 (Figure 5). Coevolution interaction clusters that are only present in the multimer, again agreed with 221 known interfaces in structures. with a coevolution cut-off of 2.5 or more deviations above median 222 coevolution (unadjusted Z score) and monomer contact map cut-off threshold of 12Å between backbone 223 residues, thus allowing removal of contacts explained by monomers. 224 225 226 227 Figure 5. Multimer coevolution and ipTM comparisons, SpoIIIE, MlaD, Tail Fibre 228 (A) SpoIIIE (FtsK homologue) Bacillus subtilis motor domain, (left) MMSEQ2 clustered sequences with 229 more than 50% coverage Red – Low coverage/conservation, Blue/Green – High 230 coverage/conservation, (centre left) Coevolution, Dark Blue – coevolution high, White – no observed 231 coevolution Z score min=0, max=3 (centre right) Multimer only coupling contacts (red) , monomer (blue) 232 (right) Contacts in 3D to self between chains above Z score 1 (B) MlaD Escherichia coli monomer 233 contrast to crystal PDB ID, (C) Enteroccus Phage tail PDB ID 4UXG (D) AlphaFold3 ipTM score of 234 each protein oligomer at user input oligomeric states. 235 236 When comparing contact maps generated from monomeric protein models, small differences in 237 predicted structures at the tertiary level such as backbone bond angles and lengths can create false 238 0 0.2 0.4 0.6 0.8 0 5 100 0.2 0.4 0.6 0.8 1 0 500.10.20.30.40.50.60.7 0 5 10 MlaD Tail FibreSpoIIIE(motor)D SpoIIIE(motor) MlaD Tail Fibre ipTMscore Oligomeric stateOligomeric stateOligomeric state A B C contacts. Subtly different monomers or an array of monomers in alternative conformations can produce 239 a larger interaction surface for self-contact, thus resulting in unexplained contacts before true 240 ‘quaternary contacts’ can be discerned. To avoid this significant larger clusters indicative of a new 241 protein interface, should be considered first. Therefore, as a final step of inter-molecular interaction 242 selection for biased quaternary contact prediction, coevolution links are plotted in 3D using a generous 243 12 angstrom or higher distance threshold to eliminate conformer errors. The distance threshold defines 244 the minimum and maximum distance between the backbone of chains in the hypothetical multimer 245 meaning that only interactions between subunits within the threshold are annotated, and not potential 246 differences in tertiary interactions, where interactions may be extremely far apart i.e. above 30 247 Angstroms Å. 248 249 To ensure coevolution analysis is identifying only real interfaces within oligomers ‘Template modelling’ 250 (TM) scores of the inter-protein interactions and in intra-protein interactions; ipTM or pTM respectively 251 should be used as an additional guide, where available. We compared ipTM scores of three proteins at 252 different oligomeric states MlaD, SpoIIE and Phage tail protein of a virus (17), and found that the highest 253 ipTM was observed at the correct multimerization state (Figure 5D), corresponding with coevolution 254 data. MlaD had the highest ipTM in both hexamer and heptameric forms, the tail phage protein as a 255 trimer, and SpoIIIE motor as a heptamer, with coevolution couplings at this interface across model’s 256 indicative of multimerization (highlighted in Figure 5 zoom boxes). The combined use of coevolution 257 with the structural predictions allowed for accurate re-production of existing crystal structures or cryo-258 EM structures, which cannot be gleaned from a single method alone. This suggests that in the absence 259 of any experimental data both tools need to be combined to accurately predict protein structure and 260 oligomeric states. 261 262 In order to investigate oligomer prediction of larger complexes, the combined prediction approach was 263 also applied to Bacillus subtilis protein SpoIIIAG (Figure S1) (18). This assembly was not predictable in 264 the state-of-the-art AlphaFold3 using public accounts due to character limits. However, coevolution 265 software is not limited in the same manner. The AlphaFold ipTM suggested the asymmetric hexameric 266 intermediate (Figure S1A) was the best fit oligomeric state for the SpoIIIAG. However, when the 267 coevolution contacts and symmetry restraints are considered the 30-mer becomes the favoured 268 assembly, consistent with experimental data (Figure S1B vs S1D). Thus, highlighting that using 269 coevolution to analyse the interaction surface of these larger assemblies combined with predicted 270 structures of smaller ‘symmetrical units’ provides the most biologically relevant outcome for protein 271 structure prediction. 272 273 Mapping networks of protein interactions using coevolution: CoEVMapper 274 275 In addition to modelling protein structures and predicting oligomeric states, coevolution provides a 276 powerful tool for resolving protein interaction networks. Experimentally, it is incredibly challenging to 277 identify a complete network of protein interaction partners within a realistic time scale. Often there are 278 many candidates of potential interactors, and limited genetic phenotypic or biochemical data to validate 279 the interaction, other than co-occurrence of genes. STRING (20,21) provides an excellent resource in 280 gene co-occurrence and neighbourhood which can be applied to discern the likelihood that two genes 281 interact (19,20). However, it does not take into account the product proteins and the role that slight 282 mutations within genes can make in perturbing protein interactions. Further analysis with coevolution 283 data can resolve these features and identify new pathways that are undetected using gene information 284 alone. 285 286 A key issue with coevolution is that it is difficult to compare at the scale of genomes due to the 287 computational power required to look at permutations. We wanted to leverage coevolution, similar to 288 how STRING can leverage co-conservation and gene neighbourhoods, to create gene cluster 289 mappings. The idea is simple, we would use the framework of CoEVFold to create an alignment using 290 MMseqs2 (Figure 6A), and a direct coupling contact map (Figure 6B), then rank whole protein inter-291 protein interactions by the summed coevolution scores. Due to the box sizes and RAM limitations of 292 large contact maps (4000 x 4000 amino acids) which crash typical servers available to users on 12GB 293 RAM python notebooks. We reduced the contact map size and the map was altered to be a 294 representation of a bootstrapped, randomly sampled coevolution between a set of genes (set to 800 by 295 800 grids). Each square within the grid was representative of a region in the protein. Finally, the 296 coevolution of each protein relative to sampled proteins was compared to respective intra-protein 297 coevolution, to identify potential interaction partners. A ‘high coevolution region’ correction was applied 298 using average product correction (apc), meaning proteins which are intrinsically well conserved or 299 overrepresented in sequences will show to what degree their coevolution is intra and inter protein 300 explained (Figure 6C). 301 302 Our pipeline instigated in CoEVMapper works well for multiprotein systems between 600 and 8000 303 amino acids in size. Tested examples included the Division and Cell Wall cluster of Bacillus and its 304 closest clades. Currently the system is limited by MMseqs2 inputs as the alignment provider, although 305 this too could be bootstrapped for extremely large queries in the future, or alternatively other filtering 306 methods such as FilterDCA(21). CoEVMapper resolved interactions with peptidoglycan modification 307 genes of Bacillus subtilis and its relatives, a system well studied experimentally. We find that DivIB and 308 SpoVE-PbpB interact, FtsA and FtsZ preferably interact, and the MreB Mur-ligase cytoskeletal proteins 309 interact with MreC among all three groups, data congruent with the literature(22). We also analysed a 310 eukaryotic gene cluster (Supplementary Figure 2) using CoEVMapper to compare the coevolution of 311 the RNA polymerase from Homo sapiens. Several known core interactions were successfully replicated, 312 particularly those of the well-studied set of proteins (PCF11, CLP1, SSU72, SYMPK, CSTF1, CPSF2, 313 SCAF8 and POL2RA). This highlights that the CoEVMapper method can be applied to study interaction 314 networks across a range of species. Calculating a network however, promiscuity correction could be 315 preferred, as the network interactions can vary greatly dependent on the conservation and therefore 316 level of data in ‘self-coevolution’ a gene has eg: MurE and PbpB are more highly conserved than other 317 proteins. However, this in itself provides information on how some proteins evolve (Figure 6C). This is 318 applied to create a map of hierarchy among coevolution interactions (Figure 6D). 319 320 321 322 Figure 6. Mapping gene networks using coevolution, CoEVMapper 323 (A) MMSEQ2 alignment of input DCW cluster fasta, Red – Low coverage/conservation, Blue/Green – 324 High coverage/conservation, (B) – compressed coevolution (50 bootstraps) Dark Blue – coevolution 325 high, White – no observed coevolution Z score min=0, max=3 (C) (left) Normalised interaction score to 326 self (centre) abundance of coevolution above threshold Z score compared to row (right) Weighted 327 A B C D interaction score (D) – Network of coevolution interactions, line thickness indicating coevolution signal. 328 PbpB coevolution with Mur ligases and SpoVE, FtsA linking with FtsZ. 329 330 331 Coevolution mapping of gene networks needs to be done with caution. The process of coevolution 332 refinement of raw data, means that the signal is ‘normalised’ and therefore higher signals are quenched 333 and lower signals in local regions are amplified. This bias can be seen with PbpBs relative propensity 334 to show coevolution interactions with other proteins compared to itself as quite highly scoring, whereas 335 MreB the inverse is true. Bias can be minimized by limiting Z scores to above a user threshold which 336 we have set in the default parameters to above a Z score of 1 or standard deviations above the normal 337 of one and counting the number of interactions and average. Future work will target optimizing these 338 parameters within the program to streamline the network mapping and reducing the risk of user induced 339 bias. Nonetheless coevolution maps from CoEVMapper provide a useful tool for quickly identifying 340 protein interaction networks for hypothesis generation and experimental rationalization. It also provides 341 a platform that streamlines data generation for introduction into other prediction tools such as Alphafold. 342 In our software we provide both weighted and unweighted versions, dependent on the average Z score 343 above the GREMLIN threshold, across lanes. The best approach is not yet clear, and a focus of future 344 work. For example, clade specific use of coevolution and coverage as well as identify cut-offs may be 345 able to get around the biasing of extremely well conserved proteins. 346 347 CONCLUSION 348 349 We applied established coevolution prediction methods to generate 2D contact maps of genes of 350 interest. Using a combination of MMseqs2 and the GREMLIN direct coupling analysis algorithm 351 coevolution profiles can be readily generated, providing a powerful tool for inferring biology when 352 combined with protein structure modelling tools or experimental data. There are a variety of applications 353 for coevolution, three of which we have incorporated into pipelines; CoEVFold, CoEVFold:4D and 354 CoEVmapper. These tools effectively map and validate protein interaction interfaces both within 355 homomeric and heteromeric complexes and finally provide a hypothetical gene network visualiser. Our 356 structural interface pipelines in ‘CoEVFold’ predict in vivo and in vitro results similar to previous tools, 357 both for homomeric multimerization and heteromeric multimerization. Coevolution is well established in 358 studying protein complexes; however, our goal was to develop a streamlined pipeline that, together with 359 the examples we provide, encourages researchers to explore both heteromeric and homomeric 360 interactions in the context of coevolution and is usable for the unspecialised investigation. By doing so, 361 we aim to reduce the “black box” nature of structural interaction predictions, giving users tools to verify 362 or provide supporting evidence for potential multimeric links within the context of coevolution. Finally, 363 we highlight that coevolution could have unprecedented potential for mapping larger scale gene-cluster 364 wide interactions. Ultimately, we hope that our software grounded in established algorithms will help 365 shed light on how complex structures predicted by new AI tools are established and serve as another 366 means of validating their accuracy, whilst also facilitating new hypothesis generation using genetics as 367 a baseline. 368 369 METHODS 370 371 Prediction of protein structures 372 373 All protein structures and multimers listed, unless otherwise stated as an existing pdb structure were 374 created using submissions to the AlphaFold3 server with no modifications. B subtilis SpoIIIAG was 375 predicted using residues 40-189. The B subtilis SpoIIIE motor domain was predicted using residues 376 353-789. 377 378 CoEVFold and CoEVFold 4D 379 380 CoEVFold works in the following manner. First the user opens a Jupyter notebook explorer in Google 381 Colab, and navigates to CoEVFold, similar to the interface of the widely used ColabFold(23). Then the 382 user uploads a protein structure, or a protein model with two or more proteins interacting in a 383 hypothetical orientation. The query sequences are then derived from the input pdb file and aligned by 384 MMSEQs2 using user defined custom parameters. We recommend a >50% coverage threshold and 385 95% similarity, as these ensure that the majority of the protein sequences used are similar to the protein 386 of interest and that coevolution is not dominated by a specific clade. After alignment, MMSEQs2 387 provides an .a3m alignment file of the input amino acid sequence, with the sequence of all amino acids 388 (independent of protein number) conjoined for each species. The conjoined multiple sequence 389 alignment is then subject to direct coupling analysis by a Markov Random Field as originally described 390 by GREMLIN (1). In other words, this coevolution pair algorithm is applied to determine amino acids 391 that co-evolve, and a coevolution score between 0 (no contact) and 1 is determined between the 392 residues in the conjoined sequence. The regions of coevolution can be mapped onto the 3D structure 393 of the protein, with the threshold Z score controlling coevolution displayed. This Z score is based on 394 average intramolecular interactions, as well as intermolecular interactions in the case of heteromers. 395 We recommend based on existing data for inter-protein interactions between heteromers, to use a Z 396 score three deviations above the mean score to reduce false hits. A higher number of sequences input 397 originally will correspond with a higher confidence in the Z scores estimated. For example, consider a 398 protein complex, composed of 4 molecules of protein A and 2 molecules of protein B, obtained by 399 AlphaFold3 or RosettaFold. If coevolution is present at the subunit interfaces, CoEVfold will primarily 400 display coevolution interactions between the chains of each protein in the output .pse file (Figure 2). A 401 higher concentration of high scoring coevolution contacts increases the likelihood of the interaction 402 having occurred over a longer period of evolutionary history. In addition, CoEVFold provides a contact 403 map and a list of the most coevolved residues, which can aid in experimental design such as site-404 directed mutations to probe oligomeric state. 405 406 In CoEVFold, any interprotein interactions which have coevolution and are in contact are included as a 407 connection dependent on the user input coevolution and back-backbone angstrom thresholds. Default 408 parameters; GREMLIN score above 3 and backbone-backbone contact threshold of 15Å. The top-409 ranking coevolution connections only visible in the multimer are noted and listed as an output, along 410 with their coevolution score, and distance dashes measurements are drawn in PyMOL, to create a pse 411 output of the overall structure explained by coevolution. This shows which amino acids have a high 412 coevolution score and are in contact dependent on user threshold. CoEVFold also notes the reliability. 413 Eg if the number of sequences provided is less than the total amino acid length, the score is unreliable. 414 Eg: if 100 sequences are used to find the coevolution of a 300 a.a total length protein complex. 415 CoEVFold is limited by the contact map matrix size, therefore any interactions which involve complexes 416 in excess of 1000 amino acids may cause the colab to crash, this is a limitation for coevolution discovery. 417 A user could use higher GPU limits possible on A100 GPUs, to assess large complexes such as the V-418 O type Atpase visualised in Figure 2. 419 420 In CoEVFold 4D, coevolution is scored for each amino acid pair (users are only able to give an input 421 sequence) and users may provide a pdb of the monomer or other lower order multimer (for example 422 from AlphaFold database or a known monomer dimer) and optionally that of the hypothetical higher 423 order multimer, please note: These files must be trimmed to have the same amino acid length as the 424 input sequence. The 12Å carbon backbone to backbone map is derived from the monomeric and 425 multimeric maps. The backbone map is then used as a mask for identification of monomeric and 426 multimeric coevolution interactions. Any residue-residue contacts which have coevolution and are in 427 contact in each structure are included as a connection dependent on the user input coevolution and 428 back-backbone angstrom thresholds. Default parameters GREMLIN score above 3 and backbone 429 backbone contact 15Å. The top-ranking coevolution connections only visible in the multimer are noted 430 and listed as an output, along with their coevolution score, and distance dashes measurements are 431 drawn in PyMOL, to create a .pse output of the overall structure explained by coevolution. This shows 432 which amino acids have a high coevolution score and are in contact dependent on user threshold. 433 434 There are two sub-modes within this homomer search mode on CoEVFold: 4D. 435 436 a. Confirmatory quaternary evolution – Here the user provides the input sequence, 437 which is aligned using MMseqs2 (Figure 3A) and coevolution is calculated as 438 previously described (Figure 3B). Then the user provides both a monomeric or lower 439 order .pdb file of the protein, as well as a hypothetical higher order multimeric .pdb file 440 predicted by a structural prediction tool such as AlphaFold or TrRosetta (Figure 3C/F). 441 The program finds user thresholded coevolution only explained by the contact maps of 442 the multimer, and plots these on to the structure, saving the result as a PyMOL readable 443 pse file. (Figure 3G). In multimer verification and analysis mode, the contact map of the 444 multimer is also compared (Figure 3H). The coevolution score is plotted on the same 445 contact map for both inputs for plotting graphically, red indicating multimer interactions 446 and blue monomer interactions. 447 b. Unguided quaternary evolution – Here the user provides only a low order or 448 monomeric .pdb file of the protein; They can search for unexplained coevolution which 449 might be explained by higher order structures, visualised on a 2D contact map (Figure 450 3E). The 12Å contact map of the monomeric protein is then used to filter out only 451 coevolution not explained in this 3D map, this gives a set of residues and their scores, 452 which are yet to be explained, which if there is a multimer, especially in areas of high 453 clustering, may be where a multimer forms. 454 455 CoEVMapper 456 457 CoEVMapper uses a user input sequence separated by ‘:’ and gene names separated by ‘-‘ to establish 458 an input of genes and their respective sequences. It then simplifies these to below a 600.a.a total 459 sequence using random sampling and applies this over a user set number of bootstraps to identify 460 coevolving regions CoEVMapper relies on paired-unpaired coevolution. Finally, the coevolution is 461 compared, using user settings. The code is available on Colab 462 https://colab.research.google.com/drive/1MSSvNTq7KZ4Lr0XTz89vUuK-J3xOTzwS?usp=sharing., 463 and Github; https://github.com/MishterBluesky/CoEVFold/tree/main 464 465 466 ACKNOWLEDGEMENTS 467 468 We would like to thank Sergei Ovchinnikov for his helpful insights into his Github page, and those 469 mentioned in the ‘Beerware license’ of the GREMLIN algorithm. ‘Sergey Ovchinnikov and Peter Koo, 470 as well as Hetu Kamisetty’ who wrote the first GREMLIN algorithm for coevolution searching. We would 471 also like to thank Phillip Stansfeld for his helpful insights and guidance and Melissa Webby for her 472 internal review of the work. This work has been supported by the BBSRC grant BB/X008533/1 awarded 473 to C.D.A.R. We would like to thank the Rodrigues lab for their help in editing the manuscript and testing 474 of the tools, in particular Daniel Dunbar. 475 476 REFERENCES 477 478 479 1. Balakrishnan S, Kamise2y H, Carbonell JG, Lee S, Langmead CJ. Learning genera=ve 480 models for protein fold families. Proteins: Structure, Func=on, and Bioinforma=cs. 481 2011 Apr 25;79(4):1061–78. 482 2. Ovchinnikov S. Protein structure determina=on using evolu=onary informa=on. 483 Thesis. 2017; 484 3. Baek M, DiMaio F, Anishchenko I, Dauparas J, Ovchinnikov S, Lee GR, et al. Accurate 485 predic=on of protein structures and interac=ons using a three-track neural network. 486 Science (1979). 2021 Aug 20;373(6557):871–6. 487 4. Sjodt M, Brock K, Dobihal G, Rohs PDA, Green AG, Hopf TA, et al. Structure of the 488 pep=doglycan polymerase RodA resolved by evolu=onary coupling analysis. Nature. 489 2018 Apr 28;556(7699):118–21. 490 5. Kilian M, Bischofs IB. Co-evolu=on at protein–protein interfaces guides inference of 491 stoichiometry of oligomeric protein complexes by de novo structure predic=on. Mol 492 Microbiol. 2023 Nov 30;120(5):763–82. 493 6. Nygaard R, Graham CLB, Belcher Dufrisne M, Colburn JD, Pepe J, Hydorn MA, et al. 494 Structural basis of pep=doglycan synthesis by E. coli RodA-PBP2 complex. Nat 495 Commun. 2023 Aug 24;14(1):5151. 496 7. Quadir F, Roy RS, Halfmann R, Cheng J. DNCON2_Inter: predic=ng interchain contacts 497 for homodimeric and homomul=meric protein complexes using mul=ple sequence 498 alignments of monomers and deep learning. Sci Rep. 2021 Jun 10;11(1):12295. 499 8. Steinegger M, Söding J. MMseqs2 enables sensi=ve protein sequence searching for 500 the analysis of massive data sets. Nat Biotechnol. 2017 Nov 16;35(11):1026–8. 501 9. Mohammadi T, Van Dam V, Sijbrandi R, Vernet T, Zapun A, Bouhss A, et al. 502 Iden=fica=on of FtsW as a transporter of lipid-linked cell wall precursors across the 503 membrane. EMBO Journal. 2011; 504 10. Cho H, Wivagg CN, Kapoor M, Barry Z, Rohs PDA, Suh H, et al. Bacterial cell wall 505 biogenesis is mediated by SEDS and PBP polymerase families func=oning semi-506 Autonomously. Nat Microbiol. 2016; 507 11. Gandarilla-Pérez CA, Pinilla S, Bitbol AF, Weigt M. Combining phylogeny and 508 coevolu=on improves the inference of interac=on partners among paralogous 509 proteins. PLoS Comput Biol. 2023 Mar 30;19(3):e1011010. 510 12. Rodriguez Horta E, Weigt M. On the effect of phylogene=c correla=ons in coevolu=on-511 based contact predic=on in proteins. PLoS Comput Biol. 2021 May 24;17(5):e1008957. 512 13. Jiang Y , Marszalek PE. Atomic force microscopy captures MutS tetramers ini=a=ng 513 DNA mismatch repair. EMBO J. 2011 Jul 20;30(14):2881–93. 514 14. Jean NL, Rutherford TJ, Löwe J. FtsK in mo=on reveals its mechanism for double-515 stranded DNA transloca=on. Proceedings of the Na=onal Academy of Sciences. 2020 516 Jun 23;117(25):14202–8. 517 15. Wotherspoon P , Johnston H, Hardy DJ, Holyfield R, Bui S, Ratkevičiūtė G, et al. 518 Structure of the MlaC-MlaD complex reveals molecular basis of periplasmic 519 phospholipid transport. Nat Commun. 2024 Jul 30;15(1):6394. 520 16. Granell M, Namura M, Alvira S, Kanamaru S, Van Raaij M. Crystal Structure of the 521 Carboxy-Terminal Region of the Bacteriophage T4 Proximal Long Tail Fiber Protein 522 Gp34. Viruses. 2017 Jun 30;9(7):168. 523 17. Abramson J, Adler J, Dunger J, Evans R, Green T, Pritzel A, et al. Accurate structure 524 predic=on of biomolecular interac=ons with AlphaFold 3. Nature. 2024 Jun 525 13;630(8016):493–500. 526 18. Zeytuni N, Hong C, Flanagan KA, Worrall LJ, Theiltges KA, Vuckovic M, et al. Near-527 atomic resolu=on cryoelectron microscopy structure of the 30-fold homooligomeric 528 SpoIIIAG channel essen=al to spore forma=on in Bacillus sub*lis. Proceedings of the 529 Na=onal Academy of Sciences. 2017 Aug 22;114(34). 530 19. Szklarczyk D, Morris JH, Cook H, Kuhn M, Wyder S, Simonovic M, et al. The STRING 531 database in 2017: Quality-controlled protein-protein associa=on networks, made 532 broadly accessible. Nucleic Acids Res. 2017; 533 20. Szklarczyk D, Gable AL, Lyon D, Junge A, Wyder S, Huerta-Cepas J, et al. STRING v11: 534 Protein-protein associa=on networks with increased coverage, suppor=ng func=onal 535 discovery in genome-wide experimental datasets. Nucleic Acids Res. 2019; 536 21. Muscat M, Croce G, Sar= E, Weigt M. FilterDCA: Interpretable supervised contact 537 predic=on using inter-domain coevolu=on. PLoS Comput Biol. 2020 Oct 538 9;16(10):e1007621. 539 22. Graham CLB, Newman H, Gille2 FN, Smart K, Briggs N, Banzhaf M, et al. A Dynamic 540 Network of Proteins Facilitate Cell Envelope Biogenesis in Gram-Nega=ve Bacteria. Int 541 J Mol Sci. 2021 Nov 27;22(23):12831. 542 23. Mirdita M, Schütze K, Moriwaki Y , Heo L, Ovchinnikov S, Steinegger M. ColabFold: 543 making protein folding accessible to all. Nat Methods. 2022 Jun 30;19(6):679–82. 544 545 546 547 548 549