Enzyme Similarity Networks

There is a crescent use of enzymes in multiple industries and sciences, ranging from materials and fuel synthesis to pharmaceutical and food production. Their applicability in this variety of fields depends not only on their biochemical function but also on their physicochemical properties. In the present work, we describe how the coincidence methodology can be employed to construct similarity networks of seventy well-studied enzymes of the Glycoside Hydrolase Family 13 and to identify communities of physicochemically related enzymes. More specifically, each of the selected enzymes is mapped into a network node, while the links between pairs of enzymes are determined by the coincidence similarity between selected physicochemical features of interest. The obtained networks have modularity and number of isolated nodes optimized respectively to two parameters involved in the coincidence methodology, resulting in highly modular networks. In order to investigate the effect of the considered physicochemical features on the enzymes relationships, the coincidence-based method also is applied to create a meta-network, in which the enzymes similarity networks obtained by the combination of every possible feature becomes nodes of a feature combination network, and the coincidence similarity between those networks defines the respective links. The obtained feature combination network systematically and comprehensively indicates the impact of the selected physicochemical features on enzyme similarity. Several interesting results are reported and discussed, including the identification of subgroups of enzymes with similar physicochemical features within catalytical classes, providing important information for the selection and design of enzymes for targeted biotechnological applications.


Introduction
Several of the most complex real-world systems can be found in biology, extending from micro to macroscopic time and spatial scales. Examples at micro and mesoscopic level include the immune system, the basal metabolism, and the nervous system. At more macroscopic levels, we have food webs and ecology.
A great deal of the biological complexity stems not only from the vast number of elements that are typically interacting, but also from the diversity of types of these entities. Of exceptional scientific relevance are the systems that involves proteins, which appear in biological systems in an impressive diversity of sizes and types, with multiple functionalities and physicochemical peculiarities.
Proteins are polymers of amino acids, each chosen among 22 possible types. Each amino acid has specific electrical, chemical, and other physical properties that dictate how they behave and how the protein folds in 3D space. The protein and all its resulting features on the 3D space -including but not limited to its geometry, topology, volume, stability, and solubilitywill depend on its amino acid composition and will enable its biochemical function [1]. Thus, the vast number of possibilities to assemble amino acid chains is directly connected to the numerous possible activities that proteins can perform [1].
An activity of particular importance for proteins is their capacity to act as catalysts. This activity is crucial for the existence of life, as almost all biochemical reactions in living systems depend on proteins to occur at the rates needed to maintain biological systems [1]. This role is so important that proteins that act as catalysts are specially called enzymes.
There are four perspectives from which the catalytic activity of enzymes can be better appreciated. The first is that enzymes increase the rates of very specific chemical reactions by selectively interacting with substrates, stabilizing transition states, and releasing chemical products [1]. The second is that the same enzyme can have more than one catalytic activity. For instance, the neopullulanase enzyme of Geobacillus stearothermophilus can participate in both the converting of pullulan in panose and the converting of cyclomaltodextrin in linear maltodextrin [2]. The third perspective is that the same enzyme type can exist in several distinct organisms, extending along with both related and unrelated phylogenetic levels. For example, pullulanase enzymes can be found in Bacillus subtilisand and Thermococcus kodakarensis, two bacteria species, as well as in Oryza sativa and Spinacia oleracea, two plant species [3]. The fourth perspective is that enzymes can be applied outside of their source organism to catalyze reactions of biotechnological interest [4].
The catalytic activity of an enzyme stems from two principal aspects: (i) its 3D geometric shape; and (ii) the physical and chemical properties along its chain of amino acids [1].The 3D geometric shape is very important, as the interactions between the enzyme and its substrates occur at specific active sites in the enzyme structure [1]. For instance, the active site of an enzyme can be a concave and positively charged cavity that can interact with a convex and negatively charged substrate, making it stable and prone to a chemical reaction such as hydrolysis. However, although there are multiple computational and experimental techniques to study enzyme structure, this still is an arduous task.
For this reason, the study of the activity of an enzyme and its possible applications in biotechnological systems usually begins with the characterization of its physical and chemical properties [4]. Subsequent analysis involves organizing a set of enzymes of particular interest according to their mutual physicochemical similarity. This is of particular interest because enzymes with the same catalytic activity and with similar physicochemical properties can be synthesized, purified, and applied in similar conditions, thus being interchangeable in their biotechnological applications [4].
Over the last decade, multiple efforts have been carried out to classify enzymes based on their functional characteristics. The main one is the creation of the Enzyme Commission (EC) number system by the International Union of Biochemistry and Molecular Biology (IUBMB) [5]. This commission labels enzymes with a four-part EC number based on their catalytic activities. For example, the neopullulanase enzyme of Geobacillus stearothermophilus mentioned above receives the EC numbers 3.2.1.135 for its neopullulanase activity and 3.2.1.54 for its cyclomaltodextrinase [3]. These four sequential numbers have different meanings: the first indicates the main catalytic class of the catalysis (e.g., 3 is for hydrolases); the second is for the subclass based on the target chemical structure (e.g., 3.2 is for glycosylases -targets glycosyl groups); the third is for further information on the chemical reaction (e.g., 3.2.1 is for glycosidases -targets O-and S-glycosyl groups); and the fourth is a serial number for the specific reaction [4,6]. This classification method is very convenient for scientifical and biotechnological purposes, as enzymes are promptly grouped with other proteins with the same biochemical applications, irrespective of their original species. The EC numbers system can be easily accessed on the ExplorEnz database, available at https://www.enzyme-database.org/ [6].
Other important efforts on curating and classifying proteins are: the BRENDA database (https://www.brenda-enzymes.org), with nomenclature, functional, kinetic, organism-related, structural and experimental information on enzymes [7]; the BioCatNet database system (https://biocatnet.de/), focused on catalytic and amino acid sequence information of enzymes for engineering purposes [8]; the carbohydrate-active enzymes database (CAZy) (http://www.cazy.org/), specialized on the organization of genomic, structural and biochemical features of carbohydrate-active enzymes (CAZymes), as well as their classification into protein families [3]. These resources can assist biotechnology researchers to both prospect and engineer enzymes with better features for their intended application [9].
The knowledge available on those databases could be better explored by identifying communities of enzymes with similar properties within the groups assembled through phylogenetic, structural and catalytic analysis. For this purpose, in the present work, we apply the coincidence methodology [10] to obtain enzyme similarity networks based on physicochemical data. The coincidence similarity was developed as an extension of the Jaccard similarity index capable of being applied to real values [11,12] and able to incorporate information about the relative interiority (or overlap, e.g..g. [13]) between two datasets. When compared to more traditional similarity indices as the cosine similarity and correlation, the coincidence index has been found to allow more strict quantification of the pairwise similarity between two datasets [12,14].
The coincidence methodology has been proposed [10] as a means to translate datasets of entities characterized by respective features into networks presenting enhanced connectivity details and modularity. It has been applied with encouraging success to several types of data [15,16] and interesting problems including city characterization [17] and motifs identification [18].
In the present work, we apply the coincidence methodology to translate seventy well-studied enzymes of biotechnological interest into similarity networks. Each node of the similarity network corresponds to an enzyme, and each connection between pairs of nodes are specified by the pairwise coincidence similarity between their respective physicochemical featuresnamely number of amino acids, molecular weight, theoretical isoelectric point (pI), instability index, aliphatic index, and grand average of hydropathicity (GRAVY) [19].
The translation procedure involves two parameters, namely α and T , respectively controlling the relative contribution of negative and positive feature values and the threshold for removing weaker connections [10]. Different properties of the similarity networks can be optimized while systematically varying α and T [10,18]. In order to find subgroups of connected entities using similarity networks, one important property that can be optimized is the modularity. However, often the created network with the parameter configuration leading to maximum modularity is also characterized by a large number of isolated nodes. In order to identify the parameter configurations leading to high (but not necessarily optimal) modularity while simultaneously minimizing the number of isolated nodes, we adopted a supervised combined optimization index reflecting both these sought properties.
The enzymes selected for this work are CAZymes of the glycoside hydrolase family 13 (GH13) [3]. They cover four types of catalytic activities [3], each with particular biotechnological interests: hexosyltransferases (EC 2.4.1.-), applied in the production of sugars with novel properties, both in material and food sciences [20,21,22]; phosphate α-maltosyltransferases (EC 2.4.99.16), enzymes of potential pharmaceutical interest [23,24,25]; glycosidases, enzymes that hydrolyse O-and S-glycosyl compounds (EC 3.2.1.-), applied in starch and baking industries [26]; and intramolecular transferases of maltodextrin, sucrose or maltose (EC 5.4.99.-), used for the production of valuable sugars, such as trehalose [27]. In the translation of these enzymes into similarity networks, we calculated the combined optimization index based on the EC number of each enzyme, which was tuned by controlling the α and T parameters. The proposed method led to enzyme similarity networks with both a small number of isolated nodes and high levels of catalytic similarity within the obtained clusters.
The results typically obtained in supervised classification are known to depend critically on the set of features adopted to characterize the data elements of interest (i.e. individual enzymes). In the present work, we employ the feature analysis methodology described in [16] as the means for studying, in a systematic and complete manner, the contribution of each of the considered physicochemical features on the obtained enzyme similarity networks. More specifically, each of the possible combinations of physicochemical features is considered for obtaining optimal enzyme similarity networks, and the similarity between each of these networks is then gauged by using the coincidence index, resulting in a new network that expresses in a comprehensive and systematic manner the effects of considering different sets of features on the construction of enzyme similarity networks. Therefore, one of the contributions of the approach reported in this work is the inspection of how the adopted physicochemical features, calculated based on the properties of each amino acid that forms the enzymes, can for groups of enzymes within and among the four original types of catalytic activities. In this sense, a good discrimination between enzymes of the four original types obtained by using the coincidence approach would indicate, among other implications, that these measurements can be used as subsidies for further characterizing the physicochemical space of each enzyme classes.
Several relevant contributions and results are reported and discussed in this work. Combined with the potential of the coincidence methodology for providing more strict quantification of the similarity between two mathematical structures, the use of the composed optimization index favoring simultaneous maximization of the modularity and minimization of the number of isolated nodes yielded markedly uniform groups of connected enzymes (respectively to their EC number) while minimizing the number of isolated enzymes. In this work, this is achieved by controlling the contribution of the features with the same or opposite signs on the resulting similarity as provided by the α parameter. Remarkably, the optimal configuration could not have been achieved using the standard coincidence formula without resorting to the α parameter. These results provide an important subsidy for the selection and design of enzymes for biotechnological applications based on both their catalytic activity and physicochemical profile.
The current work begins by presenting the adopted enzyme data and basic concepts and methods to be applied, including the coincidence similarity index, network modularity and its optimization, as well as the methodology for systematic analysis of the effect of the adopted features on the obtained protein similarity networks.

Materials and Methods
In this section are presented the model enzymes under study, their respective physicochemical features, the coincidence similarity and its application to translated datasets into networks, the concept of network modularity and its optimization, as well as the description of the systematic methodology for studying the effect of the adopted features on the obtained enzyme similarity networks.

The Model Enzymes
The seventy enzymes considered for this work were selected from the CAZy database, which has been used as a tool for the study of CAZymes since 1998 [3]. This database categorizes CAZymes into five classes, each class with several families of structurally-related proteins [3]. One family of great biotechnological interest and with extensive experimental data available is the GH13, of the glycoside hydrolases class, composed of an evolutionary diversified set of enzymes that hydrolyze glycosidic bonds [28].
The CAZy database has 129333 assigned enzymes under the GH13 family, from which only 147 are both characterized and have their structure available [3]. In this work, we will refer to these 147 enzymes as "well-studied".
Furthermore, the GH13 family has 30 different assigned EC numbers [3]. To simplify this collection of catalytic activities, we shorten the EC number to its third digit (e.g., the EC number 3.  Table 1.
Among the 147 well-studied enzymes, we selected 70 model enzymes based on their evolutionary diversity, reported cases of biotechnological interest, and quality of the data available on CAZy [3] and UniProt [29] databases. Additional information on the adopted model enzymes is in Supplementary Material 1. The distribution of the model enzymes into the four enzyme types, the description of their catalytic activity, and their potential biotechnological uses are presented in Table 1. Table 1: Distribution of the 70 model enzymes into the four types of enzymes under study. The model enzymes were selected from a pool of 147 well-studied CAZymes from the GH13 family. Each enzyme type corresponds to a simplified EC number relative to a set of similar catalytic activities. The asterisk indicates that all the well-studied enzymes available of that specific type have been used in this work.  [27] We collected the amino acid sequences of all chosen model enzymes in Uniprot [29] by using the accession numbers provided by CAZy [3]. The amino acid sequences were used to calculate six physicochemical features important for the synthesis, purification and application of enzymes, namely: 1. Number of amino acids -It is one of the most fundamental features of a protein. Since it is directly proportional to the length of its respective gene, this feature plays a role in assembling the right expression system, which can have restrictions based on gene length [30]. Furthermore, the number of amino acids is important for size exclusion chromatography, a technique that separates proteins based on their capability to infiltrate and be retained in a porous gel medium, which depends on protein size and weight [31,4]. The target proteins can be purified by selecting a gel with compatible pore size, thus the purification of enzymes with similar sizes via molecular exclusion chromatography can be performed using the same chromatography equipment.
2. Molecular weight -Just like the previous, this is an important feature of size-exclusion chromatography [31,4]. Furthermore, the molecular weight can be exploited in ultrafiltration, a size-exclusion technique that consists in passing the solution with the target protein through a semipermeable membrane that retains molecules with a minimal molecular weight [4]. Additionally, enzymes can be retained in semipermeable membrane to be used in membrane reactors [4]. Therefore, enzymes with similar molecular weights could be purified or used on membrane reactors using the same semipermeable membranes.
3. Theoretical pI -The pI of an enzyme is the pH when its net charge is 0. This property is important for setting the experimental conditions of ion-exchange chromatography. This type of chromatography is used to separate and purify a target enzyme. It consists in passing the solution that contains the desired enzyme in a matrix with selected ion-exchange groups [31,4]. If the target enzyme is in a solution with a pH above its pI, it will have a negative surface charge and will be trapped in a matrix with a positively charged ion-exchange group (anion-exchange chromatography), while if the target enzyme is in a solution with a pH below its pI, it will have a positive surface charge and will be trapped in a matrix with negatively charged ion-exchange group (cation-exchange chromatography) [31,4]. Therefore, proteins with similar pI could be purified from a solution using a similar chromatography apparatus. Such principle is also used in bifunctional or mixed-mode chromatography [4] 4. Instability index -Is a fundamental property used to predict the stability of the enzyme in vitro [32,19]. Enzymes with an instability index lower than 40 are considered stable [32], and lower values of this index are expected to be related to more stable proteins. Enzymes with similar instability index should demand similar experimental and industrial treatments during their purification and application -e.g., suitable temperature to perform purification processes, possibility to perform denaturation fractionation, use of protease inhibitors in solution to avoid the destruction of the desired enzyme [33].
5. Aliphatic index -This feature corresponds to the relative volume occupied by aliphatic side chains, and it indirectly indicates the thermostability of the enzyme [19,34]. Similar to the previous index, enzymes with the same aliphatic index should have similar responses to temperature, thus being an important characteristic to evaluate the suitability of an enzyme for industrial applications.
6. Grand average of hydropathicity (GRAVY) -It indicates the average hydropathicity, i.e. the hydrophilic or hydrophobic behavior, of the amino acids that composes the enzyme [19,35]. A positive GRAVY value indicates that the protein is hydrophobic and a negative value that it is hydrophilic [19,35]. Enzymes with similar GRAVY values will have a similar behavior toward the water, thus enabling the use of the same technical apparatus for purification through hydrophobic/hydrophilic interaction chromatography [4,31,36]. Additionally, enzymes with similar GRAVY will have similar solubility to different solvents, therefore increasing the possibility to apply either of them in the same solution conditions.
The calculations were performed with the ProtParam tool, available on the ExPASy server (https://web.expasy.org/protparam/) [19]. Although several experimental methods are available to estimate the physicochemical properties of enzymes, this procedure was used to standardize the data acquisition for this large dataset.

The Coincidence Similarity
Quantifying the pairwise similarity between two mathematical structures represents one of the most essential operations in science and technology. While the Jaccard [37] and other related indices (e.g. [13,38,39]) have been often adopted for quantifying the similarity between two sets (categorical or binary data), concepts such as the Pearson correlation and cosine similarity have been widely employed for measuring the similarity between two realvalued feature vectors (e.g. [40,41,42,43]).
The coincidence similarity index has been proposed as an extension of the Jaccard index in the sense of being applicable to real-valued data and being able to reflect the relative interiority between the datasets [11]. The extension of the Jaccard index to real-valued data has been accomplished through multiset theory, suitably adapted to take into account negative values [11,44,45]. The Jaccard index has been shown [11] to be unable to take into account the relative interiority, or overlap (e.g. [13]) between two sets. This can be addressed by multiplying the Jaccard index with the interiority index, therefore yielding the coincidence similarity index.
The potential of the coincidence index for more strict quantification of the similarity between two generic mathematical structures has been verified respectively to several comparative approaches (e.g. [14,45]) and applications (e.g. [18,17]).
In addition, the real-valued coincidence similarity index can be adapted, with the introduction of a parameter α, so as to allow the contributions of pairs of features with the same or opposite signs to be taken with different weights respectively to the resulting coincidence values [12,11,45].
A method for translating datasets, with data elements characterized by respective features, into respective detailed networks has been reported [10] which is based on assigning a node to each of the data elements, while the pairwise links between nodes correspond to the coincidence similarity between the respective standardized (or not) features. This method has been verified to allow networks characterized by a marked level of interconnection detail as well as by enhanced modularity [10]. As discussed in [45], the coincidence index naturally performs the comparison relatively to the magnitudes of the vectors, being little susceptible to perturbations or noise on single features.
Only two parameters are involved in the coincidence methodology, namely the above mentioned α as well as an overall threshold T which is eventually adopted in case binary connections are required. By varying α and T , the coincidence approach to translating dataset into networks can be optimized respectively to sought properties such as network modularity, number of components, among many alternative possibilities.
In its parameterless version, the coincidence similarity index between two vectors X and Y can be defined as corresponding to the product between the respective Jaccard and interiority indices, namely: The real-valued coincidence can be immediately obtained as: The interiority (or overlap, e.g. [13]) index, in its real-valued version [12], can be expressed as: where S is the support of the vectors, and The real-valued Jaccard index can be defined [12,11] as: with 0 ≤ J R (X, Y ) ≤ 1. Now, the parametrized real-valued coincidence index can be written as: where: The parameter α, with 0 ≤ α ≤ 1 provides an interesting manner to control the relative contributions of the pairwise features combinations with the same or opposite signs onto the resulting similarity values. In particular, for α > 0.5, the pairwise combinations of features with the same sign will have larger relative contribution to the resulting index than those with opposite signs. By varying α, it becomes possible to optimize the obtained coincidence similarity values with respect to several indices, such as the modularity of the resulting networks [10].
It can be verified that when α = 0.5, the parameterized real valued coincidence becomes identical to its parameterless version.

Translating Datasets into Networks by Using the Coincidence Methodology
As proposed in [10], the coincidence index can be effectively employed to translate datasets into complex networks. Each of the N data elements x i is characterized by M respective features j, 1 ≤ j ≤ M , organized into respective feature vectors, which are then standardized respectively to each of the features so as to yield null means and unit variance. The standardization is performed respectively to each feature j of each data element i as: where µ i and σ i are the means and standard deviation along the values of the feature j.
The coincidence similarity is then calculated respectively to each pair of data elements, yielding a respective adjacency matrix which can be thresholded by T (the coincidence values of the entries larger or equal to T are kept). The parameter α can be set so as to optimize some property of interest, such as the network modularity or the number of isolated nodes.
Provided the parameters α and T are set in a suitable manner, this simple methodology has been found [10,14,15,17,18] to allow networks to be obtained that are characterized by enhanced level of details and modularity. These enhanced features are a direct consequence not only of the more strict ability of coincidence to quantify similarity [12,45], but also of the possibility to control, by setting the parameter α, the relative contribution of features with the same or opposite sign on the resulting similarity value.

Network Optimization
Conceptually, a community or module of a network is a set of nodes that are more intensely interconnected with one another than with the remainder of the network [46,47]. Networks that have communities, which are said to be modular, appear frequently in real-world systems, including social networks, scientific collaboration, as well as several biological networks, including protein interaction and, very probably, also enzyme similarity [48,49].
In order to quantify how modular a network is given a proposed division criterion, a possibility is to calculate the network modularity [47,46,50]. The modularity quantifies the difference between the number of edges that connect nodes within and among each candidate module. There are several ways in which the modularity of a network can be defined, but, in this work, we adopt the equation [50]: Where −1/2 ≤ Q(C) ≤ 1, C is a proposed community structure, c is a module, m is the number of nodes in the network, E(c) is the number of edges within the c module, and deg(v) is the total degree of the vertex v [50]. This simple equation can be used to optimize the modularity of a network: the first term should be maximized with modules that have many edges among their nodes, while the second term should be minimized with the module having nodes with a low degree. An equivalent and computationally efficient approach to calculate this equation [47] was used in this work.
Although there are several manners of defining communities on a network based on its properties [46,51], the community structure for biological networks usually has a biological meaning transcending its topology. In the particular cases of enzymes, this has to do with their catalytic activity. In our case, the enzyme similarity networks will be divided into four modules respective to the four enzyme types under study. Therefore, this proposed community structure is expected to be preserved meaning that, in this work, we aim to determine suitable interconnections between the network nodes which lead to a significant agreement with the original four modules. This will be approached by using only the coincidence method while varying the α and T parameters. In our analysis, the modularity Q(α, T ) can be seen as a function of these two parameters.
An unwanted consequence of the optimization of the equation 10 is that networks with a low quantity of edges and with multiple isolated nodes may have high values of modularity. Although this is not necessarily a problem, the proposed coincidence method to create enzyme similarity networks are intended to be used to create modular networks -thus with biologically relevant information -with a reduced number of isolated nodes -thus integrating most of the enzymes under analysis. Therefore, we propose to combine the increase of the modularity Q(α, T ) with the minimization of the number of isolated nodes I into a single combined optimization index F [52]: The index F (α, T ) makes possible to both minimize the relative number of isolated nodes I(α, T )/m and maximize the modularity Q(α, T ) of the network (considering the biologically relevant community structure C), solely by varying (α, T ).
In brief, in the present work, given any set of the 6 physicochemical features under study, we apply the coincidence methodology to create 121 dif-ferent networks between the 70 adopted model enzymes, each network corresponding to a pair of (α, T ), with α and T assuming values in 0.0, 0.1, 0.2, ..., 0.9, 1.0. Then, we estimate the optimization index F for each network considering the division of the nodes into four enzyme types (modules), and we define the optimal network as the one with (α, T ) that yielded the minimal calculated value of F.

Feature Analysis
Before real-world or abstract patterns can be recognized or classified, they need to have their most relevant properties translated into respective mathematical measurements, or features. While there is no definite method for selecting which are the best features relatively to each pattern recognition problem, it becomes of particular relevance to devise means for studying the effect of specific combinations of features into the obtained result(e.g. [53]). Of particular relevance is the potential impact of different feature choices on the performance of the implemented approaches respectively to each particular dataset.
In the present work, each of the 70 adopted model enzymes is characterized by 6 physicochemical features. In order to investigate how possible combinations of these features impact the obtained enzyme similarity network, we adopt the coincidence-based approach reported in [10]. This approach involves obtaining optimal similarity networks for each of the possible feature combinations and then mapping them into nodes of a feature combination network, in which the links between any two nodes are determined by the coincidence between the adjacency matrix of the similarity networks they represent.
The feature combination network expresses in a quantitative and effective manner the effect of each feature in the obtained similarity networks. This type of network has been found to be strongly modular respective to specific sets of features [10,17,18,16], so that each of the identified communities of features can be taken as a set of related models for describing the dataset. A community that is strongly connected within itself indicates particularly homogeneous models, characterized by intense synergy between the features involved.
In the present work, we consider the optimal enzyme similarity network for each of the possible combinations of the six physicochemical features under study. In order to provide a direct indication of the effect of the considered feature combinations on the obtained enzyme similarity networks, the respective optimization index is also shown in the resulting feature combination network.

Results and Discussion
In this section, we derive the optimal enzyme similarity network from the adopted six physicochemical features respectively to the adopted optimization index (which needs to be minimized) and investigate the impact of the combinations between the six features on the resulting networks. We start by presenting, in Figure 1, the two-axes Linear Discriminant Analysis (LDA, e.g. [53]) projection of the 70 model proteins considering the adopted six physicochemical features. The LDA consists of a supervised approach in which the original data elements are projected into a lowerdimensional space to maximize their separation based on the respective intra and intercluster distance matrices (e.g. [53]).
It can be readily observed from this result that the LDA projection is unable to effectively cluster the four types of enzymes, being mostly devoid of clustering tendency. This result is partially due to LDA being a lowerdimensional projection of six-dimensional data, data that represents physicochemical and non-catalytic characteristics of the enzymes. Therefore, the already expected weak interrelationship between different EC groups in a six-dimensional space may have been weakened by a substantial loss of information. In addition, LDA is based on correlation-related similarity, which has been found to be less strict than the coincidence approach (e.g. [14,12,45]).
Next, we construct 63 optimal enzyme similarity networks respectively to the adopted combined optimization index F(α, T ), each one of them respective to different possible combinations of the six physicochemical features under study. This optimization took place over the two-dimensional space defined by the parameters α and T . The 63 optimal networks were gathered in a feature combination network in Figure 2. Figure 2(a) depicts, in respective heatmap colors, the modularity values obtained for each of the 63 combinations of the six physicochemical features. Figure 2(b) shows the same feature combination network, but with heatmap colors indicating the resulting number of isolated nodes for each feature combination, and Figure 2(c) presents the respective combined optimization index values in heatmap colors. Figure 2: Feature combination network, the network of coincidence similarities between the 63 enzyme networks respective to each possible combination of the physicochemical features. Each vertex corresponds to an optimal coincidence graph of enzymes created based on the labeled physicochemical features, while the edges represent the coincidence index between pairs of graphs. The groups of connected vertexes are labeled as Feature Groups (FG). Each vertex is colored according to the (a) modularity, (b) the number of isolated nodes, and (c) the combined index of their respective enzyme similarity graph. The label numbers represent the following physicochemical features: 1 -number of amino acids, 2 -molecular weight, 3 -theoretical pI, 4 -instability index, 5 -aliphatic index, and 6 -GRAVY. This network was obtained for α = T = 0.5.
A large connected component can be readily identified in Figure 2, corresponding to enzyme similarity networks that are similar to one another, as well as some smaller groups and isolated nodes. Each of these groups can be understood as corresponding to a possible model of the enzyme similarities, in the sense that the nodes within each of these groups will lead to enzyme similarity networks that are similar [16]. Distinct connected components, therefore, provide alternative models. Each of these models can be summarized in terms of the respective hub [17].
The feature groups FG9 are worthy of special attention. This group corresponds to the largest obtained connected group, composed of combinations encompassing all the six adopted physicochemical features, which also happens to include one of the largest modularity values. The minimum observed combined optimization index F is to be found within this connected group -respectively to the features 1, 3, 4, 5, and 6. Figure 3 presents the number of occurrences of each of the six physicochemical features within each of the three largest feature groups, namely FG22, FG6, and FG9. These histograms indicate that all the six features are almost equally frequent in FG9. A similar situation, except for a substantially less frequent presence of the first feature, can be observed for FG22, while FG6 contains mostly the three first features as well as one occurrence of the fifth feature. This shows that these six physicochemical characteristics are mostly interchangeable for creating enzyme networks with similar connections, at least when a minimum of three of them are considered in the calculations. Figure 3: Histograms indicating the frequency of occurrence of each of the adopted six physicochemical features within the main models FG22, FG6, and FG9. The optimal enzyme similarity network was identified within FG9. Table 2 presents the five feature combinations (nodes in the network of Figure 2) leading to optimal enzyme similarity networks with the lowest observed values of the combined optimization index. It is presented the α and T values considered to achieve each optimal enzyme similarity network, as well as its respective modularity, number of isolated nodes, number of components, combined index, and feature group (FG). Except for the fifth entry in this table (features 1, 4, and 6), all other configurations involve at least four out of the six adopted physicochemical features.
According to the combined optimization index method, the best enzymatic similarity network achieved was constructed considering the features 1, 3, 4, 5, and 6, with α = 0.9 and T = 1.0. Next, we proceeded to analyze this combination of features now with multiple values of α and T , in order to study how the resulting enzymatic network behaves with the alteration of these two parameters. Figure 4(a) depicts the modularity values, in heatmap colors, of the enzyme networks obtained for several values of α and T considering the features 1, 3, 4, 5, and 6. The respective number of isolated nodes are shown in Figure 4(b), and the combined optimization index in (c).  Interestingly, the region of largest modularity corresponds to a diagonal ridge extending roughly from the bottom-left to the upper-right portions of the graph in Figure 4(a). This indicates that larger modularity values tend to occur for mutually similar, though not identical, values of α and T . It is worth noticing that this modularity ridge coincides with an abrupt transition (a 'cliff') of the number of isolated nodes taking place along the same region in Figure 4(b). This fact tends to imply the highest modularity values to involve an intermediate number of components (orange points in (b)). The combined optimization index, shown in Figure 4(c), corresponds to a smooth ridge.
Of particular relevance is the observation that the optimal configuration was obtained for α = 0.9, being therefore otherwise impossible to be identified by using the parameterless coincidence index equation (Eq. 5), which would keep α fixed ad 0.5. This result corroborates the importance of resourcing to the α parameter while translating datasets into respective coincidence networks. Figure 5 shows the scatter plot of the number of isolated nodes in terms of modularity obtained for the enzyme networks of multiple α and T configurations using the features 1, 3, 4, 5, and 6. The optimal combined index configuration is shown in purple, corresponding to α = 0.9 and T = 1.0. Other combinations of special interest are also highlighted in the figure. It is interesting to observe a relatively well-defined non-linear (V-like) relationship between the two optimization parameters.
The enzyme similarity network yielding the optimal (minimum) combined index is presented in Figure 6.
Several relatively large modules can be identified in Figure 6, most of which are highly uniform in the sense of involving mainly enzymes that catalyze the same type of reaction. Several of the obtained groups involve strongly interconnected enzymes, indicating respectively high uniformity of physicochemical properties. Interestingly, some of the obtained communities, such as PG5 and PG14, consist of a string of connected nodes, which therefore indicates transitive similarity between the involved enzymes, with the enzymes corresponding to the two extremities being the least similar among the other enzymes in that group. Also of particular interest is the fact that all the four enzyme types yielded more than one related connected component, which indicates the presence of subgroups of physicochemical characteristics Figure 6: The enzyme similarity network with the optimal combined index, obtained with the features 1, 3, 4, 5, 6. The nodes are labeled with the respective Uniprot accession numbers of the enzymes. The connected components are labeled as Protein Groups (PG). Most of these components are strongly interconnected, except for PG5 and PG14, which correspond to chains implied by transitivity of the coincidence pairwise interactions. Of particular interest is the good homogeneity of enzyme types within each identified connected component, as well as the fact that those four types have yielded separated clusters, suggesting the presence of subcategories within types of enzymes.
within the original enzyme types.
From the biological perspective, the relatively good agreement between the obtained groups and the four original enzyme types with the formation of pure type PGs -such as PG15, PG11, PG12, and PG7 -substantiate the potential of the physicochemical features 1, 3, 4, 5, and 6 in reflecting, as well as providing subsidies for predicting and enhancing, the catalytic activities adopted for the definition of the four original enzyme types.
Two biochemical explanations can possibly justify the observed relationship between physicochemical properties and enzyme types. First, these five physicochemical features are directly calculated from the amino acids that compose the protein sequences [19], which will also dictate the folding of the enzyme in its final three-dimensional structure and, consequently, its catalytic activity [1]. Second, the physicochemical characteristics of an enzyme are the result of its evolution, which is strongly dependent on its biological function, i.e., its catalytic activities [54]. Although enzymes with similar catalytic activities should not necessarily have similar physicochemical characteristics [54], the existence of this type of relationship is not unexpected. As shown in [55], physicochemical properties can be used, with the help of machine learning algorithms, to classify proteins into their correct families. Figure 7 illustrates the enzyme similarity networks obtained for several equally spaced values of the parameter α, more specifically α = 0.7, 0.8, 0.9, 1, with T fixed at 1.0. The protein similarity network obtained for α = 0.7 is composed almost exclusively of individual, disconnected nodes, except for a few groups involving from two to four proteins. As α is progressively increased, the nodes become more and more interconnected, while the combined index reaches its maximum at α = 0.9, then decreasing for α = 1.0. At the same time, the number of connected components decreases monotonically with the increase of α. Even though the combined index values obtained for α = 0.9 are smaller than the peak value, the respectively obtained networks can still be observed to contain mostly uniform groups of nodes, that are progressively interconnected. It is interesting to keep in mind that, as α increases, with new links being incorporated into the network, the previous connections are maintained. Figure 8 depicts the range of values for each of the five considered physicochemical features characterizing the PGs with at least four enzymes: PG5, PG11, PG12, PG14, and PG15.
According to Figure 8, the enzymes of each group have quite a similar number of amino acids, thus having similar protein and gene sizes. Therefore, the enzymes of a given group could be purified with the same size exclusion chromatography equipment, while different groups would require different chromatography gels. Because the enzymes of the groups PG5 and PG11 have pI values varying from neutral to slightly acidic conditions, these proteins could be purified through both cation-and anion-exchange chromatography, depending on altering the pH of the solutions to acid or basic, respectively. However, the enzymes of groups PG14, PG12, and, particularly, PG11, because they have lower values of pI, should be more easily purified with anion-exchange chromatography.
Furthermore, by choosing enzymes of the group PG14 or PG15, the  biotechnologist will have enzymes with low instability index and high aliphatic index, thus having enzymes with the highest predicted stability for uses in industrial processes. However, those enzymes might have a slightly hydrophilic behavior, thus they could present problems with solubility. Alternatively, enzymes of the PG12 group might have more stability problems, but a higher hydrophilic behavior. Overall, the above analysis illustrates how the proposed methodology could help assemble groups of technically related enzymes for biotechnological uses while taking into account specific enzyme properties and requirements.

Concluding Remarks
Proteins, which include enzymes, have been largely acknowledged to constitute the building blocks of life. Their essential role lies in the fact that they implement different types and levels of functionality across species and biological systems, including many functionalities useful for biotechnological applications. As a consequence, the study of not only the catalytic activities of proteins but their physicochemical properties and applicability outside its original species, as well as how the proteins can be grouped into classes are ubiquitous in biological, biophysical, and biochemical sciences.
The present work has addressed the classification of enzymes based on their physicochemical properties relying on the following three key aspects: (i) the coincidence method is employed to map from several combinations of enzyme features to respective similarity networks; (ii) this mapping is optimized with respect to an optimality index related to the modularity and number of isolated nodes of the resulting networks; and (iii) a coincidencebased approach is applied in order to characterize in a systematic and comprehensive manner the impact of the enzyme features on the resulting optimal similarity networks.
Several interesting results are reported and discussed, including the iden-tification of enzyme similarity networks characterized by relatively high modularity levels and small number of isolated nodes obtained from six physicochemical properties of the considered enzymes. The several combinations of the six features under study were systematically evaluated, and the optimal enzyme similarity network was obtained while considering features 1, 3, 4, 5 and 6, leaving out only the molecular weight. This result suggests that, among the adopted physicochemical features, the five selected measurements have more pronounced potential for characterizing subgroups of physicochemically related enzymes within the four types of enzymes.
Another particularly interesting result concerns the importance of using the parameters α and T while optimizing the properties of the networks obtained by the coincidence methodology. More specifically, it has been shown that the optimal enzyme network could not have been otherwise obtained with the parameterless (i.e., with α = 0.5) version of the coincidence similarity. At the same time, the mapping of the modularity and number of isolated nodes in terms of α and T provided interesting insights suggesting that the latter property defines a cliff that nearly coincides with the main diagonal characterized by α = T , while the maximum modularity configurations tending to occur along the respectively defined ridge. Therefore, the maximum modularity tends to take place at positions in the domain (α, T ) where the number of isolated nodes undergoes an abrupt variation, which reminds a phase transition.
The described methodology and results pave the way for several interesting prospects for further related studies. For instance, the reported approach can be directly applied to other important enzyme and protein families, as well as to other types of biological molecules and interactions. It would also be interesting to consider additional physicochemical and catalytic properties, as well as measurements reflecting the geometrical structure of the considered molecules. Furthermore, this methodology could be incorporated as a complementary tool in prospecting molecules of biotechnological and pharmaceutical interest or in the rational design of proteins and ligands to improve their physicochemical properties.