To store the contact types and their conditions, we used a dictionary containing all heavy atoms from the 20 standard amino acids, as defined in (Sobolev et al., 1999; Silva et al., 2019; Fassio et al., 2020; Barroso et al., 2020; Pimentel et al., 2021; Dos Santos et al., 2022). All 20 standard amino acids had their heavy atoms classified by the following characteristics, in binary form (Table 1): tendency to contribute to hydrophobic effects, belonging to aromatic groups, having positive charge, having negative charge, capability of donating or accepting electrons. The full atom classification table is available in the Supplementary Table S1.

The possible contact types are: hydrogen and disulfide bonds; hydrophobic effects; attractive, repulsive, and salt bridge interactions; and aromatic stackings. This dictionary also contains the conditions needed for the contact (e.g., to form an attractive interaction, the atoms must be differently charged), and the range of Euclidean distances, in angstroms, for the contact to occur (Table 2).

The full PDB protein archive, in ‘.cif’ format, was obtained using in-house scripts to query and download entries directly from the PDB API. First, a script was used to query the API for entries containing “Protein” as an exact match from the parameter “entity_poly.rcsb_entity_polymer_type”.

To avoid rate limits and overwhelming the server, queries had a 1 s delay from one another, and only 25,000 IDs were obtained at a time. Then, a second script was used, together with the Biopython Bio. PDB module (Cock et al., 2009), to download all files that matched the IDs gathered in the first step. All files were downloaded between July 4th and 10 July 2024.

To serve as a comparison to our method, the Biopython package (Cock et al., 2009), largely used in bioinformatics, was used. The Bio. PDB module contains tools to parse a. pdb or. cif file, as well as the NeighborSearch (NS) class, which is useful in interatomic contact determination.

We used an in-house Python script to perform an all-atom neighbor search of 6Å radius, the maximum distance for contacts defined in our dictionary. Then, the neighbors were filtered based on their distance and physicochemical properties relative to the parent atom. Redundant comparisons were excluded; for example, if atom “ a ” was identified as a neighbor of atom “ b ,” the comparison was performed only once, preventing the redundant evaluation of “ b ” as a neighbor of “ a .” The code for the NS implementation is available at the Supplementary GitHub repository.

The contacts obtained contained the following information: chain, residue number, and parent atom name; chain, residue number, and neighbor atom name (i.e., the atomic pair making the contact); type of contact; and distance between the two atoms.

To analyze the PDB protein archive and obtain the maximum distances matrix used in the rest of this work, we first devised a Static Cutoff (SC) implementation, where the C α cutoff distance was fixed. Akin to Biopython, proteins are treated as Python objects, containing chains, residues, and atoms. The package includes a customized. pdb/.cif parser, optimized to rapidly extract only the information relevant for contact determination. This makes it more efficient and lightweight than general-purpose parsers, which are typically designed to support a broader range of structural analysis tasks. By default, the parser considers the following criteria: a) Only the first model of each protein is considered (in the case of proteins experimentally resolved by NMR); b) Only atoms with occupancy ≥ 0.50 are considered; c) Water molecules, hydrogen atoms, non-standard residues, nucleic acids (DNA and RNA), and metallic coordination are not considered.

After parsing, the protein object is passed to a contact calculation script, where the C α distances for each pair of residues are obtained, and filtered based on the fixed cutoff. Centroids of aromatic rings were calculated using all atoms belonging to the ring, and the calculation of normal vectors and angles was performed using the Python NumPy library.

The atoms from the residues that are in range to interact are then compared to the dictionary previously described, based on their distance to each other, and their physicochemical properties. Finally, the contacts are returned in a custom object containing all their information, similar to the NS method.

Throughout the processing of the complete PDB archive using the SC method, the maximum distances (across all proteins in the PDB) between the C α atoms of each amino acid pair were stored in a distance matrix. Upon completion of the processing, this distance matrix was then used to update the static cutoff point employed in the SC method, generating specific values for each amino acid pair.

The distance matrix D = [ d i j ] n × n is a square matrix of size n × n , where n represents the number of standard amino acids. Each entry d i j corresponds to the maximum distance between the C α atoms of the amino acids at positions i and j (e.g., d 11 represents an Alanine pair, and d n n represents a Valine pair): D = d 11 d 12 ⋯ d 1 n d 21 d 22 ⋯ d 2 n ⋮ ⋮ ⋱ ⋮ d n 1 d n 2 ⋯ d n n , (1) where each d i j represents the maximum Euclidean distance between the C α atoms of the amino acids at positions i and j .

A total of 210 distance values were obtained, representing each possible residue pair and excluding redundancies (e.g., Ala-Val is the same as Val-Ala) (Equation 2). P = n n − 1 2 + n , (2) where P is the number of non-redundant distance pairs, and n is the number of standard amino acids. In this case, as n = 20, then P = 210.

The COC α DA tool is implemented similarly to the previous implementations (NS and SC). A schematic representation of the implementation is presented in Figure 2.

First, the script “main.py” is executed via the command line, along with the required parameters, which are processed by the script “argparser.py” (1). The parameters include the mandatory file paths (wildcards are accepted), an optional binary flag for generating an output file in. csv format (default = no), an optional parameter to parallelize file processing in batches across any combination of available CPU cores, and an optional parameter to use custom contact distances defined by the user instead of those defined in Table 2 (using the “contact_distances.json” configuration file). All parameters are fully explained using the ‘-h’ or ‘–help’ flags.

When users specify custom contact cutoffs greater than the default (6Å), the static C α –C α distance matrix is extended accordingly. This is done by computing an epsilon ( ϵ ) value, which is the difference between the user-defined cutoff and the default maximum cutoff. If ϵ > 0 , it is added to all distance values in the original matrix. This adjustment ensures that the pruning based on C α distances remains valid even under relaxed contact definitions, preserving consistency with the originally computed matrix while accommodating user-defined thresholds.

The files are then processed by the script “parser.py” (2), which utilizes the previously defined classes to create objects representing proteins, their chains, residues, and atoms. Finally, the script “contacts.py” receives the objects generated by the parser (3). The scripts “conditions.py”, which stores the conditions required for contact detection, and “distances.py”, which stores the values obtained from the distance matrix, are used to compute the contacts.

If the user does not specify the. csv output file parameter, only a summary is displayed in the terminal, containing the protein name, residue count, number of contacts, and processing time in seconds. If the parameter is used, in addition to the summary, a. csv file is generated with detailed information about each detected protein contact. The columns in the output file are organized as follows: Chain 1, Residue Number 1, Residue Name 1, Atom Name 1, Chain 2, Residue Number 2, Residue Name 2, Atom Name 2, Distance, Contact Type. An example output file for PDB ID 101M, as well as a PyMOL (Schrödinger, LLC) visualization script to help users quickly explore the results, are available at the Supplementary GitHub repository.

Two datasets were selected to benchmark our results and compare them to other competitors ( n dataset1 = 896 and n dataset2 = 215,716). The first (D1) is a modified gold-standard set of enzyme superfamilies (Brown et al., 2006), with 365 unique entries ranging from 194 to 6,208 residues. For a more balanced comparison, we split all chains in different files, and treated them separately. The new modified dataset contains 896 entries, ranging from one to 994 residues (available on the Supplementary GitHub repository). The second dataset (D2) includes all PDB proteins with less than 10,000 residues, covering approximately 99.2% of all protein entries. This dataset contains 215,716 unique entries, ranging from three to 10,000 residues.

To ensure fairness and eliminate biases, all benchmarks were conducted simultaneously on a server with the following specifications: NVIDIA A100 GPU, 768 GB RAM, and a 128-thread AMD Ryzen Threadripper 5995WX processor. To prevent memory overload and parallelization issues, each process was executed on an individual core. Due to its size, dataset D2 was divided into nine batches of approximately 25,000 files each, with each batch being processed independently on separate cores.

Although multithreading and batch processing is available for all implementations (NS, SC, and COC α DA) using the Python module “concurrent.futures”, each core handled only a distinct batch to maintain consistency in the results. The total processing time for each entry was defined as the sum of the file reading, parsing, contact detection, and output generation times.

In total, 217,454 PDB entries were downloaded in. cif format, totaling approximately 450 GB (The full ID list is available on the Supplementary GitHub repository). Proteins ranged from three (PDB IDs: 1Q7O, 8DDG, 8DDH) to 503,221 (PDB ID: 8GLV) modeled amino acid residues. To obtain the values for the distance matrix, we processed all the downloaded files using a fixed C α distance cutoff of 21Å for all pairs of residues (SC). This value is comfortably above the maximum distance between the C α of a pair of arginines, the biggest residues by length, that are able to have contacts between their side-chain atoms (considering a maximum contact distance of 6Å, as per Table 2). To confirm this, we compared an all-atom approach (i.e., comparing every atom of the protein against each other, without cutoffs) to the SC approach using D1, and no contacts were missed (Supplementary Table S2).

Using the SC implementation and the 217,454 entries downloaded from the PDB, over 211 million amino acid residues and 819 million contacts were processed and identified. Along this process, we stored the maximum C α distances for every pair of the 20 standard amino acids, and after merging redundancies, we obtained 210 values in a symmetric distance matrix (Figure 3; Equation 1). The full distance table is available in the Supplementary Table S3.

As the distance matrix is color-coded based on the value of the maximum C α distance, we can quickly spot the minimum and maximum values obtained. For the lowest value encountered, we found a pair consisting of an alanine and a glycine residue, both present in the HD chain of PDB ID 6QCM, with a distance of 7.65Å between their C α ’s (Figure 4A). This is expected, as alanines and glycines are two of the smallest amino acid residues, differing only by a single CH 3 group in the side chain of the alanine, while glycine has a hydrogen atom in its side chain. However, even with this difference, the presence of the CH 3 group on the side chain of the alanine does not impact the distance between their C α atoms, only contributing to the chirality of the alanine residue. For these two residues, only hydrogen bonds are possible, as the main chain atoms are only capable of donating (main chain nitrogen) or accepting (main chain oxygen) hydrogen atoms.

On the other hand, a pair of two arginine residues, from chains G and H of PDB ID 3X0Y, represents the highest C α distance encountered, of 20.46Å (Figure 4B). This result is also expected and consistent with the fixed distance of 21Å used in the SC approach, once again demonstrating that the fixed cutoff was appropriate to yield no missed contacts. The contact itself is a repulsive interaction between two NH 2 atoms from the residues side-chains, at a distance of 5.9Å.

If the user wishes to apply custom contact distances instead of those defined in Table 2, the optional ‘-d’ flag can be used, specifying the desired values in the “contact_distances.json” configuration file. This flag extends the C α distance values in the distance matrix to incorporate the user-defined distances, ensuring that no contacts are omitted.

Use-case scenarios for modifying the default cutoff distances include adopting alternative minimum or maximum values reported in the literature, such as 6Å for aromatic stacking (Fassio et al., 2022), or 5Å for hydrophobic effects (Bickerton et al., 2011), instead of the default 5Å and 4.5Å, respectively. Another common scenario involves exploratory analyses using step-wise cutoff variations (da Silveira et al., 2009; Vangone and Bonvin, 2015).

With the maximum possible C α distances properly established for all amino acid residue pairs, we updated with the new values the “distances” dictionary from the SC implementation (Section 2.4), which before was fixed at 21Å for all residue pairs. To the joint implementation of the SC method with the distances updated from the distance matrix, we gave the name COC α DA.

By using tightly-defined, pair-specific cutoff distances, we can further enhance the search space pruning compared to using a single fixed distance threshold. Analysis of the distribution of maximum C α distances reveals that most amino acid pairs exhibit distances well below 21Å, suggesting that introducing pairwise-specific cutoffs is justified and efficient (Supplementary Figures S1 and S2). Moreover, given that dictionary lookups in Python operate with linear average-case complexity (O(1)), storing and querying 210 unique cutoff values introduces negligible computational overhead relative to using a single fixed value.

To improve efficiency, COC α DA employs a stepwise filtering strategy that first evaluates residue-level proximity based on C α –C α distances. An initial coarse filter excludes residue pairs exceeding the global maximum cutoff (20.46Å, for Arg–Arg pairs), allowing early removal of clearly non-interacting pairs. Remaining pairs are then subjected to a more stringent, pair-specific cutoff comparison using the distance matrix. Only those pairs that satisfy both criteria proceed to atomic-level evaluation to determine whether a contact is present. This tiered approach reduces the number of atomic comparisons required, helping to balance computational cost with contact detection accuracy. A schematic illustration of this process is provided in the Supplementary Figure S3.

The residue-pair-specific distance matrix used in COC α DA captures the full range of C α –C α distances observed across known protein structures that are compatible with side-chain atomic contacts. Because these thresholds are grounded in the physical constraints that govern residue-residue contacts, the method is broadly applicable regardless of the overall fold, resolution, or origin of the structural model. Any contact that can realistically occur must still conform to these spatial constraints, ensuring generalization across diverse structural contexts.

While the method is robust to structural variability at the backbone level, the accuracy of contact detection may be influenced by uncertainty in side-chain atom positions, particularly in lower-resolution or flexible regions. In such cases, performance near cutoff boundaries may be affected, which is an inherent limitation of any deterministic cutoff-dependent method.

To benchmark COC α DA against other approaches used in the literature, we selected the following: all atoms against all atoms (AllAtoms, used in (Pimentel et al., 2021)), Arpeggio (Jubb et al., 2017), Arpeggio CLI ³ (Jubb et al., 2017), Biopython Neighbor Search (NS, used in (Bickerton et al., 2011)), and Static Cutoff (SC). Other methods, like nAPOLI (Fassio et al., 2020), STING Contacts (Mancini et al., 2004), and PICCOLO (Bickerton et al., 2011), were not available at the time of search and were not updated recently, so they were not considered.

Both Arpeggio versions were too slow to process even small proteins, as our tests showed processing times of approximately 5 and 23 min for a single 1,000 residue protein (PDB ID 6RTH) for Arpeggio CLI and Arpeggio Web, respectively. For comparison, the same protein was processed in 0.62s using COC α DA. This can be due to several factors, but we believe that the explanation lies mainly in server load (for Arpeggio Web), and the several external libraries and computing time that are needed to run the more complex analysis (for both versions). In this way, since a large-scale analysis would not be feasible due to the large processing times, both versions (Web and CLI) were disregarded in the subsequent analyses.

For the AllAtoms approach, a new implementation was developed in Python, in order to incorporate all the previously defined definitions and constraints, and for the NS method, a custom implementation was created using Biopython (Cock et al., 2009), since the PICCOLO tool is currently unavailable. Thus, for D1, the following methods were compared: AllAtoms, NS, SC and COC α DA.

The first dataset contains 896 entries, ranging from one to 994 residues. The initial choice of a smaller dataset was made to include the AllAtoms approach, which is significantly slower than the other three (but still considerably faster than both versions of Arpeggio), which would make a large-scale analysis unfeasible.

In Figure 5, it is possible to see that the AllAtoms approach (orange) rapidly explodes in a quadratic curve compared to the three others, which maintain rather linear calculation times up to 1,000 residues. Once again, no contacts of any type were missed in any of the approaches (Supplementary Table S2), but the AllAtoms approach was removed from further analysis because of its performance. Comparing the faster approaches, SC (yellow) obtained calculation times 1.5x faster on average than NS (magenta), while COC α DA (cyan) showed calculation times 3.8x faster on average, obtaining the same contacts.

As the results from the first, small dataset showed a significant difference in processing times between the 3 fastest approaches, we then moved to D2, which contains 215,716 unique entries, ranging from three to 10,000 residues, making approximately 99.2% of the PDB protein archive. The choice to remove entries above 10,000 residues was made due to the nature of those entries, which are mostly protein complexes, containing several copies of each unique chain. This makes them not suitable for contact analysis directly, requiring some kind of pre-processing, like splitting only the unique chains or working with each individual protein present in the complex separately. This can also be true for entries below 10,000 residues, but we believe that this slice correctly represents the diversity of experimentally resolved protein structures.

Figure 6 shows the results for D2 when comparing the COC α DA (cyan), SC (yellow), and NS (magenta) approaches. It is possible to see that COC α DA performs better for all proteins, averaging approximately 6x faster times than NS and 2.5x faster times than SC. The SC approach performs better than NS in proteins below 5,000 residues, equal between 5,000 and 7,000 residues, and worse above 7,000 residues. However, since the vast majority of PDB entries fall within the smaller size range (approximately 97.2% of unique entries have 6,000 or fewer residues), where the performance gain of COC α DA and SC over NS is most pronounced, the high density of small proteins in the dataset significantly skews the overall average, leading to the reported 6x and 2.5x improvement.

Outliers were considered as entries that had a processing time ± 5 times the Standard Deviation for each approach, with less than 1% of entries removed. After outlier removal, it is possible to see that COC α DA has a consistent time vs. size distribution, while the other two approaches have more variation. This can be due to the tight and precise definition of cutoff distances for COC α DA, which speeds up a lot of the computation, while also limiting structural variations, like between globular and fibrillar proteins.

The protein size range of 0–1,000 residues is noteworthy, as shown in detail in the upper left corner of Figure 6. At this range, we can see the distribution pattern observed in D1, while also identifying several divergent entries in the NS approach. The divergent spike is composed exclusively of Nuclear Magnetic Resonance (NMR) resolved entries, which are usually deposited as several individual models of the same protein. Due to the nature of Biopython native parsing, all the models need to be parsed even if only the first one is of interest, unless the user creates specific functions for this purpose, thus deviating from the original implementation of the library. As these NMR entries are small, the parsing time of several NMR models outpaces the contact calculation time of the first one, leading to a spike in processing time. This does not occur in the COC α DA and SC approaches, as the customized parser handles only the selected model in the file (the default value is always the first model).

Computing interatomic contacts is inherently a quadratic problem (O ( n 2 ) ) because it requires calculating the distance between every pair of atoms in a protein, such as in the AllAtoms approach. However, sophisticated data structures, such as k -d trees, can be employed to avoid calculating distances between atoms/residues that are too far apart. This approach theoretically reduces the computational space by pruning irrelevant comparisons, leading to practical reductions in computation time and typically logarithmic ( O ( l o g n ) , in the average of cases) or linear (O ( n ) , in the worst of cases) complexity (Bentley, 1975; Berg et al., 2008). However, these complexity values refer only to the operations of search, insertion, and deletion of elements in the trees. For the construction of a new tree—for example, for a new protein—the processing is substantially greater, corresponding to a log-linear complexity ( O ( n l o g n ) ) in the best-case scenario, for balanced trees (Wald and Havran, 2006; Brown, 2015).

In the case of small entries, like most of the protein structures, the memory overhead associated with the allocation and creation of the tree usually does not outweigh the computational gains in search, insertion, and deletion operations. Thus, in addition to the high implementation complexity of k -d trees, the results may turn out to be worse than those obtained by ‘brute-force’ algorithms, provided that the latter are properly implemented. These gains are primarily concentrated in small proteins, which significantly contribute to the average observed speedup of 6x when comparing COC α DA to NS, even though the performance benefit becomes less pronounced on bigger proteins.

Our approach in COC α DA fits as an improvement of classical ‘brute-force’ algorithms (such as AllAtoms) for the calculation of interatomic contacts in proteins. This is because, theoretically, all atoms of the protein are checked at least once, but the precise definitions of cutoff distances for each residue pair significantly reduce processing time, without the additional cost of using complex data structures, as is the case with k -d trees.

In this study, we chose to evaluate the complexity of various algorithms empirically, by comparing standard methods commonly used in the structural bioinformatics community with the COC α DA method. These different methods were tested with inputs of increasing sizes (where n represents the number of residues, which have on average eight atoms each), and we analyzed the resulting fitted curves with real datasets.

The curve fittings of the three approaches against the second dataset demonstrate that both COC α DA ( R 2 quadratic = 0.99, Equation 3) and SC ( R 2 quadratic = 0.97, Equation 4) exhibit quadratic growth trends, while NS shows a linear growth trend ( R 2 linear = 0.97, Equation 5). This results demonstrate, experimentally, the nature of the contact identification functions, which are the most time-consuming operations. COC α DA and SC act basically as heavy improvements of ‘brute force’ algorithms, having a time complexity of O ( n 2 ) , while the NS contact identification function operates with a time complexity of O ( n ) , leading to its linear growth pattern. f n = 1.35 × 1 0 − 7 n 2 + 5.04 × 1 0 − 4 n − 6.36 × 1 0 − 3 ; (3) g n = 1.20 × 1 0 − 7 n 2 + 1.60 × 1 0 − 3 n − 9.18 × 1 0 − 2 ; (4) h n = 2.37 × 1 0 − 3 n + 7.94 × 1 0 − 2 , (5) where f ( n ) , g ( n ) , and h ( n ) are the best-fitted functions for COC α DA, SC, and NS, respectively, and n is the number of residues.

However, in the case of the NS approach, the memory overhead in tree construction is evident, especially for small inputs. It is possible to observe, in Figure 5, that for proteins of up to 200 residues even the most computationally expensive approach (AllAtoms) achieves lower processing times than NS. This can also be seen in the range that includes small entries in D2, although with less detail due to the massive number of points (Figure 6, detail).

The analysis of the coefficients from the obtained equations is another way to explain the poorer performance of NS, even though it grows linearly compared to the quadratic growth of the other approaches. Starting with the quadratic coefficients ( a ) , it can be noted that both are practically insignificant ( 1.35 × 1 0 − 7 for COC α DA and 1.20 × 1 0 − 7 for SC). This means that for small n values, as in D2, the quadratic growth is extremely slow, as can be seen from the slight curve in the data.

As for the linear coefficients ( b ) , there is a significant difference between COC α DA ( 5.04 × 1 0 − 4 ) and the other two approaches ( 1.60 × 1 0 − 3 for SC and 2.37 × 1 0 − 3 for NS), with NS showing the highest value. This difference, along with the low quadratic coefficients, explains how NS’s linear growth can be less efficient than the quadratic growth of the other approaches. Furthermore, as the SC and NS terms are close, between 6,000 and 7,000 residues the curves invert, representing the point where the quadratic growth of SC becomes more influential.

Finally, regarding the constant coefficients ( c ) , again there is a marked difference between the three approaches, with SC having the smallest value ( − 9.18 × 1 0 − 3 ) , followed by COC α DA ( − 6.36 × 1 0 − 3 ) , and NS ( 7.94 × 1 0 − 2 ) , the only approach with a positive c value. This result highlights the memory overhead associated with k-d tree construction, which causes even small inputs to have a baseline processing time higher than the actual contact calculation time.

As previously shown, the NS approach (using k-d trees) has a linear time complexity in the worst case. Meanwhile, COC α DA exhibits quadratic growth, yet shows lower processing times for all analyzed entries. By equating the two equations (Equation 3 for COC α DA; Equation 5 for NS), we find that COC α DA would only yield worse results for proteins with approximately 14,000 residues (Equation 6). f n = h n , when  n ≈ 14.000 , (6) where f ( n ) e h ( n ) are the best-fitted functions to represent COC α DA and NS data, respectively, and n is the number of residues.

However, only slightly more than 2,000 entries in the PDB have a size greater than 14,000 residues, representing less than 1% of all protein structures. Furthermore, all of these entries correspond to protein complexes or repetitions of the same protein, as previously discussed. Therefore, for all practical purposes of contact detection in proteins, the COC α DA approach demonstrates the best temporal performance compared to other methodologies in the literature.