{'fraction_non_alphanumeric': 0.027302466029189735, 'fraction_numerical': 0.011355057876195269, 'mean_word_length': 5.347175084847275, 'pattern_counts': {'":': 1, '<': 2, '<?xml version=': 0, '>': 6, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Rational, Reproducible, and Rigorous Computational Enzymology: The Case of Catechol-O-Methyltransferase

In order to accurately simulate the inner workings of an enzyme active site, it is crucial that not only the reactive species are modeled at the quantum mechanical-level (QM) but also any important surrounding residues, solvent, or coenzymes involved in crafting the microenvironment. As a first step towards removing subjectivity in enzyme model creation, our lab has been developing the Residue Interaction Network ResidUe Selector ( RINRUS ) toolkit to utilize interatomic contact network information for automated residue selection and QM-cluster model generation. Starting from an X-ray crystal structure of catechol-O-methyltransferase (COMT), RINRUS was used to construct several series of QM-cluster models by either incrementally building models residue-by-residue according to metric-based ranking or by adding combinations of residues grouped by common features. The reactant, product, and transition state of the methyl transfer reaction was computed for a total of 508 models, and the resulting free energies of activation and reaction were used to evaluate model convergence. RINRUS -designed models with only 200 – 300 atoms are shown to converge, significantly less than conventional radially expanding models requiring more than 400 atoms. RINRUS will serve as a cornerstone for improved and automated cheminformatics-based enzyme model design. site of catechol-o-methyltransferase, an enzyme which catalyzes the methyl transfer from S adenosyl methionine cofactor to catechol substrates. Our results demonstrate using RINRUS to rationally design small and accurate active site models.


Abstract
In order to accurately simulate the inner workings of an enzyme active site, it is crucial that not only the reactive species are modeled at the quantum mechanical-level (QM) but also any important surrounding residues, solvent, or coenzymes involved in crafting the microenvironment. As a first step towards removing subjectivity in enzyme model creation, our lab has been developing the Residue Interaction Network ResidUe Selector (RINRUS) toolkit to utilize interatomic contact network information for automated residue selection and QM-cluster model generation. Starting from an X-ray crystal structure of catechol-O-methyltransferase (COMT), RINRUS was used to construct several series of QM-cluster models by either incrementally building models residue-by-residue according to metric-based ranking or by adding combinations of residues grouped by common features. The reactant, product, and transition state of the methyl transfer reaction was computed for a total of 508 models, and the resulting free energies of activation and reaction were used to evaluate model convergence.
RINRUS-designed models with only 200 -300 atoms are shown to converge, significantly less than conventional radially expanding models requiring more than 400 atoms. RINRUS will serve as a cornerstone for improved and automated cheminformatics-based enzyme model design.

Significance Statement
The efficiency, accuracy, and replicability of enzyme simulations is often hampered by ad hoc model design. To address this problem, we have developed the Residue Interaction Network ResidUe Selector (RINRUS) toolkit. RINRUS utilizes residue contact networks to automate construction of rational quantum mechanical cluster models. This work examines this problem by computing the reaction kinetics and thermodynamics for 508 models of the active site of catechol-o-methyltransferase, an enzyme which catalyzes the methyl transfer from Sadenosyl methionine cofactor to catechol substrates. Our results demonstrate using RINRUS to rationally design small and accurate active site models.

Introduction
For nearly two centuries, the structure, function, and catalytic power of enzymes have fascinated scientists, with countless studies seeking to understand their underlying mechanisms.
Atomic-scale computer modeling of enzymes is currently a necessary part of the global multibillion-dollar research effort that aids the design of new pharmaceuticals, helps to investigate and engineer protein structure and function, and advances our understanding of the molecular basis of disease (1, 2). The importance of atomic-level simulation of enzymecatalyzed reactions was publicly acknowledged with the 2013 Chemistry Nobel Prize being awarded to Warshel, Levitt, and Karplus, who developed methods to treat the active site of an enzyme with quantum mechanics (QM) and the periphery with classical or "molecular" mechanics (MM) (3).
QM-only (also called QM-cluster), QM/MM, and ONIOM modeling are complementary approaches that have leveraged advancements in quantum mechanical theory and molecular dynamics (MD) to continually increase the ubiquity of computational enzymology (4)(5)(6). As with all forms of modeling, the comparative accuracy of a model to reality is limited by the design of the model and relevant/reliable experimental data. For simulating the active site of enzymes, it is crucial to ensure not only the amino acids directly involved with the reaction are modeled at the QM-level but also any residues, water molecules, ions, and coenzymes sterically and/or electrostatically crafting the active site microenvironment (4,(7)(8)(9). While this is a simple idea in principle, it is far harder in practice to identify rationally which residues must be partitioned into the QM level.
While ad hoc protocols exist for selecting residues for inclusion in QM-level modeling, recommendations are typically ambiguous and generally inefficient (4,7). One of the most common practices is to simply include all residues within a certain radial distance from a point, perhaps the center of mass of substrate(s) or an active-site metal center. While suitable models could be constructed this way, calibration studies have confirmed large spheres (and consequently large models) are needed for convergence of simulated enzyme thermodynamics/kinetics (8,(10)(11)(12)(13)(14)(15)(16)(17)(18). These results are perhaps unsurprising as nature does not enforce any geometric requirement to the design of an enzyme active site. Clearly, inclusion of less important residues unnecessarily increases the computational cost of any model. Attempts to quantify the importance of residues have been performed via a posteriori computations, such as QM/MM thermodynamic cycle perturbations (19,20), linear response functions (21), or Fukui/Charge Shift Analysis (14,22). However, such methods are essentially performing complicated, time-consuming computations on enzyme models in order to decide what model should be used for further complicated, time-consuming computations. Iterating residue selection processes to self-consistency is even more expensive.
Ideally, there would be a computationally inexpensive, a priori approach to enzyme model construction that utilizes structural and chemical data to rationally select residues (or parts of residues) for QM-modeling. As a potential solution for this model creation problem, our lab has been developing the software Residue Interaction Network ResidUe Selector (RINRUS) which computes a contact-based residue interaction network (23,24) and uses the data to identify and rank residues for modeling. Further, RINRUS automatically trims and caps the residues via a rules-based criterion to form appropriate models and generates formatted input files for several popular electronic structure theory packages (see SI for details). Through establishing an automated and rigorous workflow, we envision solutions to several communitywide problems including standardization of enzyme model creation, reducing learning curves for new users, minimizing trial and error using poorly or incorrectly designed models, and improving reproducibility of workflows and published results.
Additionally, enzyme models used to obtain insightful results must be reported in a reproducible manner, the simplest way being the inclusion of Cartesian coordinates for transformative leaps in computational biochemical method development will be severely hindered without ontologies, rigorous calibration efforts, and shared best practices. We argue that RINRUS is the first community-oriented software to facilitate "rational, reproducible, and rigorous" computational enzymology workflows and models.
Ideally, the RINRUS workflow would be capable of identifying a singular or handful of models that best capture the balance between maximizing the number of key residues included to simulate the active site while minimizing the size of the QM-region for computational efficiency.
This leads to questions such as what makes the enzyme model "good"? What easily obtainable metrics might be universal in computational biochemistry for ranking the importance of interatomic/inter-residue interactions? We begin to answer these questions within the context of contact-based residue interaction networks (23,24).
Previous computational studies have shown substantial variation in both ΔG ⧧ and free energies of reaction (ΔGrxn) with respect to QM-cluster or QM/MM model size. Recent results from QM/MM calibration studies suggesting that asymptotic convergence of thermodynamics/kinetics requires radial distance-based QM-regions of 400 -600 atoms (8,18,30). Unfortunately, conventional DFT calculations of 400 -600 atom models are prohibitively expensive for many research groups. The large QM-region size required to study the COMT mechanism also defies conventional wisdom that kinetic/thermodynamic properties should converge quickly as the size of the QM-region grows in a QM/MM partition. Slow convergence behavior of COMT has been attributed to the non-spherical active site, requiring an accurate description of both the Mg 2+ /catechol coordination chemistry and the electrostatic stabilization of the large SAM cofactor (30).
While the paradigm of calibrating expanding QM-regions in a radial distance-based fashion has been established to provide poor convergence for COMT, there is a surprising dearth of exploring alternatives to distance-based active site models in the literature. Using the workflow employed by RINRUS, we present the reaction thermodynamics and free energies of activation for hundreds of QM-cluster models of COMT with the goal of identifying interresidue contact features that lead to accurate and efficient models of COMT. Our cheminformatics perspective will be the first rigorous step towards a translatable and generalized computational enzymology protocol.

Results and Discussion
We began by computing a contact-based residue interaction network ( Figure 1B) for an X-ray crystal structure of human COMT (Protein Data Bank ID 3BWM), which indicated 27 protein residues and 4 crystallographic waters had contact interactions with any fragments central to the catalytic reaction (termed the "seed": SAM, CAT or Mg 2+ ). The residue contacts with the seed were classified into five different types: wide contacts, close contacts, small overlaps, big overlaps, and hydrogen bonding. All QM-cluster models of COMT were constructed using the crystallographic coordinates of these residues and, unless otherwise indicated, trimmed according to the RINRUS workflow (refer to SI). Models were expanded from the seed by one of two general ways: residues were incrementally added based upon a ranking criterion (e.g. distance from the seed, number of contacts with the seed) or groups of residues were added to the seed based upon similar residue features (e.g. type of contacts). Three additional models were also constructed using the QM-region residues from a previous QM/MM study (18) on COMT (see SI). Further details on the model building schemes beyond what will be outlined in the discussion are provided in the SI. In total, the methyl transfer transition state and connecting reactants/products for 508 unique QM-cluster models were computed. 1590 DFT-optimized stationary points were analyzed in this work.

Expansion of QM-cluster models by Ranking of Residues
We will first detail several ways COMT QM-cluster models were incrementally built-up by ranking residues. The first metric is the current paradigm of ranking residues based on their distance to the active site. Though a simple distance metric may seem straightforward, this method can be ambiguous and tricky to replicate without reporting very precise definitions of the radial origin and the thresholds for adding residue fragments or entire residues. Subtle variances in definitions might qualitatively affect which residues or atoms are captured within varying radially expanding models. For this work, 25 models were constructed with RINRUS by incrementally adding residues ranked by the shortest distance from the position of any atom of the seed to the position of any atom of the surrounding residues. The models were expanded until all residues predicted by the contact network were incorporated, encompassing a 3.10 Å expansion from any atom of the seed. In building up to the largest model, two residues (K46 and N92) with no contact interactions with the seed but within the 3.10 Å distance were necessarily included in these distance-based models.
Computed ΔG ⧧ and ΔGrxn are plotted against the distance-based expansion from the seed ( Figure 2A). As the size of the model increases, the predicted ΔG ⧧ converges (ΔG ⧧ within ±2 kcal/mol of the largest distance-based model) with QM-cluster models containing >342 atoms, but the predicted ΔGrxn does not similarly converge even with the largest distance-based models.
Some of the largest distance-based models (containing 444 and 447 atoms) incorrectly predict an endergonic reaction. Qualitatively incorrect thermodynamics corresponds to the addition of the charged residue K46, which as previously noted, is not predicted by RINRUS to have direct contact interactions with the seed. Residues without contact interactions with the seed are not expected to meaningfully contribute to model accuracy or convergence (as seen with the addition of N92 changing ΔG ⧧ and ΔGrxn by <0.2 kcal/mol in the 486-atom model). At best, the addition of non-interacting residues adds unnecessary time to the DFT simulations, but in this case the charged K46 side chain distorts the simulated microenvironment in the absence of counterbalancing residue(s). Growing QM-cluster or QM/MM enzyme models by distance-based rules is well-known to converge slowly, but until now, alternative schemes have not been proposed or at least widely used.
As a step towards rectifying this issue, we next considered the convergence of COMT QM-cluster models constructed by ranking based on the number of contacts each residue has with the seed and incrementally building models from residues with the most contacts to fewest contacts with the seed. ΔG ⧧ and ΔGrxn both converged by the 302-atom model ( Figure 2B). To establish a more appropriate reference value for the accuracy and convergence of RINRUSdesigned models, the average relative free energies of the five largest models was calculated, yielding 11.6 kcal/mol for ΔG ⧧ and -6.9 kcal/mol for the ΔGrxn. While this ranking method fares better at prioritizing residues than distance-based expansion, there are some inherent limitations.
Namely, larger residues with more surface area (e.g. lysine or tryptophan) are more likely to have more contacts with the seed and may bias the ranking compared to smaller residues.
Ranking by number of contacts with the seed also does not weight or quantify electrostatic influences (e.g. charge, hydrogen bonding, and polarity). Nevertheless, even with this nonoptimal metric, constructing models by contact count still yields impressively small, converged models.
Below, we will detail two combinatoric workflows for building models where residues are classified into sets by common contact type. The third method for ranking residues involves ordering residues by the number of times each residue appears in a unique model from the Combinatoric Scheme 2 model sets (see below and SI for details). This ranking inherently favors residues with more than one type of contact interaction. In using this residue ordering, ΔG ⧧ and ΔGrxn are converged when QM-cluster model size is greater than 300 atoms ( Figure 2C), similar to the models designed through ranking residues by total contacts with the seed. The model with the greatest overestimation of ΔG ⧧ and endergonic ΔGrxn (236 atoms) corresponds to the addition of the positively charged residue, K144. The subsequent inclusion of the negatively charged E199 residue places the predicted free energies within qualitative accuracy, reemphasizing the point that particular care in model design must be given towards charged residues and nearby residues that counter their effective charges.

Automation Versus Constructing QM-cluster Models Manually
The RINRUS package is still undergoing rapid development and needs further testing to address broader QM-cluster model design issues such as residue/substrate protonation states, orientation of explicit solvent molecules, and conformational sampling (7,9). While these factors may be manually addressed by the user, doing so places a potential bottleneck in the throughput of QM-cluster model applications. Removing this bottleneck with RINRUS must not come at the cost of producing unreliable models.
In consideration of possible differences between manual and automated model building, models built by ranking residues via their frequency of appearance in Combinatoric Scheme 2 models ( Figure 2C) were reconstructed by-hand by the PI. The models were designed without any guidance beyond the identity of the specific residues in contact with the seed and their ranked order. The results of these "bespoke" models are presented in Figure 2D and are shown to be comparable to the models built by RINRUS in Figure 2C. There is reduced fluctuation in the ΔG ⧧ for the smaller bespoke models versus comparably-sized RINRUS-generated models, likely attributable to manual sampling of conformers, a treatment not done for any of the RINRUSderived models. However, for the models greater than 300 atoms, there is no qualitative difference between the automated and the "by-hand" approach. These results demonstrate how RINRUS, even without carefully attending to residue protonation and conformational sampling, can construct QM-cluster models in a way similar to that by an experienced scientist, but which is founded on a traceable cheminformatic basis and a reproducible, rigorous workflow.

Expansion of QM-cluster Models by Residue Interaction Features
The remaining models were built up from the seed by combining residues with common features, specifically by inter-residue contact type. The contact types contain two pieces of information used in QM-cluster model construction: the section of the residue contacting the seed (classified as either residue main chain, residue side chain, or water) and the contact type (wide contact, close contact, small overlap, big overlap, hydrogen bonding). Models were constructed by taking all combinations of the contact types and, for each combination, building a QM-cluster model using all residues with the specific contact types of that combination. These models represent a combinatoric approach to building-up models by adding groups of residues by common features to the seed (Combinatoric Scheme 1, see SI for details). To further increase the number of models and dataset size, the sets of residues classified by contact types were repartitioned into a second combinatoric approach (Combinatoric Scheme 2, see SI for details), though the generation of these sets is not rigorous or necessarily applicable to other biosystems.
Given the limitations of time and resources, 69 (of 204 possible) models of Combinatoric Scheme 1 and 357 (of 736 possible) models of Combinatoric Scheme 2 have been simulated, representing all unique combination-based models up to at least 320 atoms. As our goal is identifying small, yet accurate, QM-cluster models, the cost of expanding the dataset to include hundreds of additional large models is not expected to lead to substantial improvements in analysis.
In plotting ΔG ⧧ and ΔGrxn of models built through the two combinatoric schemes ( Figure   3A and B), the models are shown to predict a wide range of ΔG ⧧ and ΔGrxn throughout the model sizes. This wide spread originates from differences in model composition rather than models optimizing into unnatural orientations, since the root mean square deviation (RMSD) of unconstrained residue heavy atoms of the geometry optimized reactant state compared to the Xray crystal structure is on average only 0.54 Å for all models ( Figure S4; standard deviation, 0.18 Å). Similar to the models built by ranking residues, there are models with fewer than 300 atoms that yield accurate predictions, affirming that cluster model convergence does not require large > 400 atom models.

Identifying Important Residues
A general clustering of COMT QM-cluster models that predict similar (though not necessarily accurate) free energies is observed in Figure 3 for both combinatoric schemes. This leads to the question of which residues are required to form an accurate model? To more clearly distinguish the clusters of models, the k-means clustering algorithm was used to partition the entire dataset of all unique QM-cluster models into six groups ( Figure 3C) based upon their predicted ΔG ⧧ and ΔGrxn (44). Though an unsupervised method was used to group the models, the identified clusters are reasonable and properly differentiate the models with both converged ΔG ⧧ and ΔGrxn (Cluster 5) from markedly inaccurate models (Clusters 1 and 4) as well as models with converged values for either ΔG ⧧ or ΔGrxn, but not both (Clusters 2, 3, and 6).
The residues that differ among the clusters give insight into which residues have a comparably strong influence on convergence. Tabulating the percent occurrence of each residue within the COMT models of each cluster ( Figure 4, Table S2), eight residues present in >90% of the Cluster 5 models are absent or have a greatly reduced presence in other clusters. For example, in the models of Cluster 1, which systematically overestimate ΔG ⧧ and 65% of which incorrectly predict an endergonic reaction, none contain E199 and only 11% contain M40.
Without these residues, the QM-cluster models are missing the stabilizing hydrogen bonding interactions between E199 and the catechol alongside the weak hydrogen bonding and hydrophobic interactions between M40 and the SAM, resulting in consistently poor accuracy with respect to the converged free energies.
Surprisingly, residues identified as particularly important for convergence are not always localized around the atoms directly involved in the methyl transfer. For instance, E90 (which is present in 99% of the models in Cluster 5 but only in <37% of the models in Clusters 3 and 4) is 10 Å from the catechol, but plays a role in stabilizing and properly orienting the SAM. Other residues apart from the eight illustrated in Figure 4 such as Y68, I91, A118, S119, and H142 are present in >70% of the models in Cluster 5 and appear to play important roles in crafting the active site microenvironment.
With residues crucial for accurate QM-cluster modeling of COMT identified, the next step is to examine contact and classification metrics to see if any were particularly suitable for predicting the relative importance of residues. For the contact classification combinations, there is unfortunately no consistent combination of contact types among the Cluster 5 models for yielding converged models. Using the contact dot frequency ( Figure 2B) as a ranking system proves modestly successful as 9 of the 13 residues present in >80% of the Cluster 5 models have a high frequency of contacts with the seed and would be correctly prioritized. The four residues with low contacts (N41, A67, Y71, A118) are adjacent to high-contact residues and largely have main chain interactions with the seed, explaining the fewer contacts. The general success of using contact dot frequency as a ranking scheme was previously shown in Figure 2B where converged models had 302 atoms as a lower bound. Improvements to this ranking method are warranted (and are under current investigation by our lab), ranging from incorporating additional chemical descriptors to the interatomic contacts (e.g., through Arpeggio), to developing a weighting system to favor certain contact interactions (e.g., hydrogen bonding, polar, aromatic) (45). In the end, RINRUS provides a computationally inexpensive, rational, reproducible means to building enzyme QM-cluster models.

Conclusions
Computational enzymology has made incredible impact in understanding the atomiclevel intricacies of enzyme function. While computational resources and scaling limitations of quantum chemistry are among factors limiting progress in this field, little attention has been given towards how poor or irreproducible model design might be hampering scientific progress.
Many publication-quality enzyme models have been founded on rationale not necessarily suited for modeling non-spherical active sites (e.g. radial distance criterion) or via rationale prone to fallibility (a researcher's chemical intuition). Techniques addressing this problem by identifying important residues a posteriori have been useful but fail to address the true need for a computationally inexpensive a priori method for designing enzyme models.
As a step towards addressing community-wide problems in computational enzymology, we have been developing the RINRUS toolkit to automate the residue selection and construction of QM-cluster models. RINRUS utilizes the cheminformatics of interatomic contact networks as the rationale for identifying active site residues and ranking/classifying them. The catalytic methyl transfer reaction of the human COMT enzyme was simulated with a total of 508 unique models, illustrating how information from RINRUS may be used to build models up from a base structure by either adding residues incrementally via a ranking scheme (e.g., total contacts with the seed) or by adding combinations of groups of residues (e.g., type of contacts). Clusters of models with common predictions of reaction and transition state free energies were compared to identify residues important for accurate simulations of COMT. Ranking residues by the frequency of their contacts with the seed demonstrated particular usefulness, with QM-cluster models with 210 -300 atoms yielding converged thermodynamic and kinetic properties.
Additionally, the methodology employed by RINRUS to identify seed-residue interactions and accordingly trim QM-cluster models favorably comparable to that of the by-hand models created by an experienced computational biochemist.
A major focus of this work has been to quickly converge energetic properties of smaller QM-cluster models to those of a maximally sized QM-cluster model. Further tests of QM-cluster modeling require comparison to well-defined experimentally known quantities. Another subsequent avenue of study will be calibration of Density Functional Theory, one-electron basis set, implicit solvation parameters, empirical dispersion corrections, and other variables of electronic structure theory to truly assess the accuracy of QM-cluster modeling beyond a metric of internal consistency. Recent developments in linear scaling coupled cluster theory suggest ways to incorporate more rigorous "black box" electronic structure theories into the realm of computational enzymology. Additional avenues currently being investigated by our lab include developing improved chemical descriptors and ranking schemes, integrating machine learning into the workflow, and expanding into automated QM/MM modeling construction. A forthcoming publication will describe the RINRUS software package and include thorough tutorials. Public availability and adoption of RINRUS will substantially reducing learning curves for new practitioners of QM-cluster modeling.
Though model design and reproducibility questions have been largely ignored within the greater computational enzymology community, we hope this work will foster self-reflection on the underlying assumptions behind how atomic-level enzyme simulations are derived. The current practices require unnecessarily large models to obtain accurate results, which is limiting progress and is undoubtedly daunting to inexperienced chemists/biochemists interested in contributing to the field. Through the automated workflows provided by RINRUS and its successful results demonstrated in this work, we present the first steps towards discovering and implementing a computationally inexpensive, cheminformatic-based means for constructing reproducible, rational, and rigorous enzyme models. Reproducible workflows in computational enzymology, supported by RINRUS development, will improve openness, data sharing, and facilitate novel cyber-and software infrastructure in biochemistry and biology.

Methods
The initial atomic coordinates for building the models were taken from an X-ray crystal structure of COMT (PDB ID: 3BWM) containing the coenzyme S-adenosyl methionine (SAM) and a 3,5-dinitrocatechol inhibitor coordinated to the active site metal (46). Hydrogens were added to this protein structure using the program Reduce (47), and the two nitro-groups of 3,5dinitrocatechol were replaced with hydrogens to form the catechol (CAT) substrate. The program Probe (48) was used to roll a small spherical probe over the van der Waals surface of this modified structure to identify and classify non-covalent interatomic contact interactions. All QM computations were performed using the Gaussian16 software package (50). The models were geometrically optimized using density functional theory (DFT) with the hybrid B3LYP exchange-correlation functional (51,52). The computations used the 6-31G(d') basis set for N, O, and S (53); the 6-31G basis set for C and H atoms (54); and the LANL2DZ effective core potential and basis set combination for Mg (55). The Grimme D3 (Becke-Johnson) dispersion correction (GD3BJ) was also included (56) along with a conductor-like polarizable continuum model (CPCM) using UAKS sets of atomic radii, a nondefault electronic scaling factor of 1.2, and a dielectric constant of ε = 4 (57,58). Unscaled harmonic vibrational frequency calculations were used to confirm all stationary points as either minima or transition states. Stationary points were found by first pre-optimizing the model to the reactant structure. This preoptimized structure was then used to construct an approximate transition state structure by translating the methyl midway between the sulfur of SAM and the oxygen of CAT and flattening the methyl to a planar configuration. The transition state was optimized, and intrinsic reaction coordinate computations were used to confirm the formal reactant and product minima and calculate reaction free energies. Whether this procedure biases the simulated active site to more strongly stabilize the reactant structure (and whether such a bias would be of any significance) is unknown and an uninvestigated facet of computational enzymology.
The k-means clustering analysis (44) was run through RStudio v.3.6.3 (59). Elbow and gap statistics ( Figure S6) were used to identify a k = 6 for the cluster analysis (60).   Values are proportionally shaded to emphasize differences in residue composition among kclusters. B) Visualization of a hypothetical model containing residues that occur in >85% of Cluster 5 models. Substrates are colored magenta.