Recent advances in user-friendly computational tools to engineer protein function

Abstract Progress in technology and algorithms throughout the past decade has transformed the field of protein design and engineering. Computational approaches have become well-engrained in the processes of tailoring proteins for various biotechnological applications. Many tools and methods are developed and upgraded each year to satisfy the increasing demands and challenges of protein engineering. To help protein engineers and bioinformaticians navigate this emerging wave of dedicated software, we have critically evaluated recent additions to the toolbox regarding their application for semi-rational and rational protein engineering. These newly developed tools identify and prioritize hotspots and analyze the effects of mutations for a variety of properties, comprising ligand binding, protein–protein and protein–nucleic acid interactions, and electrostatic potential. We also discuss notable progress to target elusive protein dynamics and associated properties like ligand-transport processes and allosteric communication. Finally, we discuss several challenges these tools face and provide our perspectives on the further development of readily applicable methods to guide protein engineering efforts.


Introduction
Proteins have been studied intensively already for several decades to reap immense benefits through their applications in green industry, biomedicine, sustainable agriculture and other areas [1][2][3][4][5][6][7]. Prominent development of methods from genetic engineering and molecular biology has laid the foundation mimics natural evolution by introducing a large number of mutations to be screened or selected [14,15]. Both engineering approaches profited extensively from progress in technology and scientific knowledge, enabling their systematic application, which led to numerous designs of proteins with improved function, solubility and stability. Their inherent requirements limit both methods. Directed evolution relies on a highthroughput system capable of evaluating large libraries of generated protein variants. In contrast, rational design requires a profound knowledge of the investigated protein and/or intensive computer simulations enabling the precise design of mutants.
The recent trend is to use both these approaches in unison, which is termed semi-rational engineering or focused directed evolution, to overcome the primary restrictions of the two approaches. In this strategy, rational components and computer predictions are used to prioritize the most promising protein sites for mutagenesis and frequently also implement restrictions on the diversity of introduced mutations to the most viable ones [16][17][18][19][20][21], resulting in smart-and-small mutant libraries with a large fraction of functional variants. With the tremendous advances in computer technology, availability of protein structure models and mathematical methods, computational tools have become indispensable for the semi-rational engineering process [22][23][24][25][26]. Such approaches continuously support the successful delivery of proteins adopted for use in various biotechnologies [27][28][29][30][31][32]. Following continuous efforts to introduce increasingly sophisticated computational approaches into semi-rational engineering, a plethora of new tools with distinct purposes and uses are being developed and released each year. This unceasing addition of available tools has positive and negative implications. On one hand, there is a tool available for almost any particular task; on the other hand, the resulting diversity of offered tools can be overwhelming for new users and active practitioners alike.
To limit the never-ending literature search and guide researchers toward appropriate tools, we have critically surveyed structure-based computational tools dedicated to protein engineering that have emerged between 2016 and 2019 (Table 1). Tools published before this period have already been thoroughly reviewed [33][34][35][36]. We have focused on recent additions to the software toolkit of user-friendly and readily applicable approaches for altering protein function that can be employed by a broad spectrum of researchers, whereas more advanced tools and methods for protein engineering and design relying on the utilization of intensive computation and expertise have been reviewed elsewhere [37][38][39][40][41][42][43]. Also, we have not covered tools for the evaluation of protein stability or solubility, as those have been reviewed too [44,45]. As bioengineering research can comprise various strategies leading to the selection of the most promising candidates with improved function of interest, we consider tools for the following sequential stages of the engineering process: (i) hotspot identification for sitesaturation mutagenesis ( Figure 1A), (ii) in silico mutagenesis to evaluate the effects of mutations and prioritize promising variants ( Figure 1B) and (iii) analysis of results to guide further engineering efforts ( Figure 1C). Additionally, we discuss an integrative computational workflow that aims at providing complete computational support for protein engineers.

Tools for hotspot identification
The initial step in a bioengineering project is to identify promising and relevant positions or regions of the protein to mutate. These sites, hotspots, are often located in or near structurally or functionally important regions of the protein, increasing the likelihood that mutations at these positions impact the protein's properties. Various strategies utilizing protein structure for hotspot identification can be applied depending on the respective target property, e.g. focusing on the proximity of binding sites [67], ligand-transport pathways [68], flexible regions [69], allosteric networks [70] or conversely structural voids [71]. Nonetheless, the majority of such tools rely on an analysis of a multiple sequence alignment (MSA) to obtain insights into the importance of the individual positions in homologous proteins. Recently developed tools such as visualCMAT [46], PDB2Graph [47], STRESS [48] and AlloSigMA [49] aim to tackle the challenging prediction of hotspots involved in allosteric communication and others, including PPI3D [50] and DisruPPI [51], focus on hotspots which govern protein-protein interactions.

Allosteric hotspots
Analysis of correlated or co-evolving residues has been the tradition in protein structure modeling to predict direct or indirect coupling between pairs or groups of residues [72][73][74]. This method has been employed to improve the quality of predicted protein structures and protein-protein complexes [75][76][77][78]. Identified co-evolving residues have also become a frequent target of engineering aiming at stability [79] and allosteric regulation [80] ( Figure 2A). The visualCMAT web-server is a recent tool focused on the analysis and identification of correlated or co-evolved hotspots [46]. This tool uses an MSA and 3D structure to assess the correlation between residues based on mutual information. Since prediction quality is dependent on the quality of the MSA, the server recommends an integrated Mustguseal web-server [81] for the preparation of the MSA, which reduces the minimal input to a PDB file. A particular limitation is that the server can evaluate only one chain at a time. As a result, correlated or co-evolved residues are identified and mapped onto the structure. These residues are then differentiated into two categories based on a predefined distance cutoff: physically interacting residues and long-range interactions. These long-range interactions may be indicative of allosteric pathways that would be difficult to identify without evolutionary information. Additionally, the server can perform a binding site prediction by identifying pockets with Fpocket [82] and ranking them by the number of present correlated residues. Outputs comprise the list of correlated pairs of positions with the corresponding statistics, the cumulative statistics enumerating involvement of a given position in all pairs found, and their structural visualization. As a test case, the authors predicted correlated residues for the FecA protein from the Porins superfamily [46]. They showed that mutations in some of the identified residues notably altered the transport function of the protein. Owing to a large number of proposed hotspots, the experimental data for validation were available only for a small fraction of the hotspots, excluding some of the highest-ranking ones [46]. Therefore, the lack of data leaves the full extent of visualCMAT applicability to be established yet.
An alternative approach to identify allosteric interactions without the need for a high-quality MSA is to consider a protein structure as a graph of residue-residue interactions ( Figure 2B). In a nutshell, protein residues are converted to nodes, often positioned on their Cα atoms, and the edges among these nodes are drawn based on distance cutoffs representing various interresidue interactions [83]. Finally, the network topology can be analyzed to reveal structurally and functionally relevant connections among residues. Many tools like NetworkAnalyzer [84] or RINerator [85] have been developed for graph-based analyses using different and often incompatible file-formats to store their results, restricting further analyses of generated graphs to specific software, e.g. RINalyzer [85]. To provide a bridge between the most common formats, a user-friendly PDB2Graph toolbox has been developed [47], which, however, depends on proprietary Matlab software. The tool produces an undirected, coarse-grained, distance-based graph that can be exported to different graph formats including Cytoscape [86], Pajek [87] and UCINET [88]. Furthermore, various centrality indices identifying residues critical for protein structure and function can be calculated. Calmodulin, phage T4 lysozyme, Barnase and Ribonuclease HI were used as test cases to assess the applicability of the PDB2Graph tool [47]. Many experimental mutations coincided with the residues identified based on centrality indices, suggesting this method can indeed be of service to protein engineers. Unfortunately, the sparsity of the experimental data did not allow systematic validation of the tool, similarly to visualCMAT.
Taking into consideration a protein not only as a network of connected nodes but also considering protein dynamics ( Figure 2C), the structurally identified essential residues (STRESS) tool aims at disclosing allosteric hotspots on the protein surface as well as in its interior [48]. For the surface hotspots, the STRESS tool employs a modified binding leverage approach [89], which was previously implemented in the SPACER web-server [90]. This method combines Monte Carlo (MC)-based ligand docking with normal mode analysis (NMA), which is a computational approach that approximates the local dynamics of a system by a harmonic motion. Outlining the principle of the STRESS method, a simplified representation of a ligand consisting of four-beads is used as a probe to identify putative binding sites on the protein surface. Then, the deformability of these sites is predicted from the 10 lowest frequency normal modes generated by a coarse-grained representation (Cα atoms) NMA provided by the Molecular Modeling Toolkit (MMTK) [91]. Finally, the putative sites are scored and ranked according to their deformability to estimate the degree that the bound ligand would interfere with predicted conformational change. In contrast to SPACER, the STRESS tool markedly reduces a large number of identified putative binding sites by considering all heavy atoms of a protein during docking and applying automatic thresholding. Additionally, STRESS also combines the above described NMA approach with a network analysis in which a protein is represented as a network of interacting residues to expand the scope of analysis to residues critical for communication along a given allosteric pathway (buried allosteric hotspots). Within this network, each edge is weighted depending on the correlation between the movements exhibited by the corresponding interacting residues during NMA. Such a weighted network is subdivided into communities using the Girvan-Newman algorithm [92], and residues critical for interconnection between these communities are detected according to their highest betweenness. For surface hotspot predictions, a list of ranked putative sites with their scores and constituent residues is produced, whereas for interior hotspots, the identity of the critical residues is reported. The applicability of the method was evaluated by its authors on 12 well-studied proteins, in which an average of 55% of known binding sites were identified correctly [48]. Further, the relevance of the identified hotspots was supported by their significantly higher evolutionary conservation in comparison to the non-critical residues calculated on a large dataset of more than 1000 proteins.
Likewise, the AlloSigMA web-server uses NMA aiming not only at the identification of hot-spots but also enabling evaluation of the effects of ligand binding or a mutation on allostery [49] ( Figure 2C). This method recognizes four structurally relevant states: (i) unbound/wild type (WT), (ii) bound to a ligand, (iii) mutated and (iv) bound and mutated. In each of these states, allosteric free energy is calculated for each residue using the lowfrequency normal modes derived from the Cα-representation of protein provided by the MMTK package and employing a previously developed structure-based statistical mechanical model of allostery [93]. Two scanning approaches can be used in order to detect allosteric hotspots: mutation-based scanning of selected regions or a whole protein, or binding-based scanning with a small probe to triplets of residues. Alternatively, allosteric effects originating from ligand binding, mutation or both combined can be evaluated. It is important to pinpoint that since the method is based on a coarse-grained representation, only two types of mutations can be considered. 'UP' mutations, by which the method emulates a mutation to a bulky residue, resulting in stabilizing effect on the local contact network, and 'DOWN' mutations, which models alanine/glycine substitutions, resulting in destabilization of the contact network. Aside from interactive visualization, the server provides PDB files with allosteric energies in the B-factor columns. The underlying method was applied to design allosteric mutations aiming at improving the activity of insulin-degrading enzyme toward amyloid β peptide [94]. Out of five constructed single-point mutants, three mutants showed up to 50% increased overall efficiency, while the other two mutants exhibited decreased efficiency [94]. Clearly, more extensive benchmarking is still required to obtain a robust estimate of the method predictive performance.
Whereas all four tools were successfully applied for the identification of hotspots, their quantitative comparison is prevented by the lack of suitable systematic benchmark datasets. Hence, the differences among tools can only be appreciated from a user perspective, considering their speed, ease-of-use and the nature of delivered results. PDB2Graph is the fastest tool, closely followed by visualCMAT, while STRESS is notably slower (Table 1 and see Supplementary Table 1 available online at https: //academic.oup.com/bib). AlloSigMA is by far the most time-demanding tool evaluated in our review, and its allosteric scanning is restricted to a maximum of 2000 residues per analysis. On the other hand, AlloSigMA is the only tool capable of delivering quantitative evaluation of the effect of the mutation on allostery and is much appreciated by the community (Table 1 and see Supplementary Table 2 available online at https: //academic.oup.com/bib). Regarding userfriendliness, only the STRESS tool does not provide interactive analyses of results, and its installation might turn out to be a bit complicated due to its dependence on somewhat older Python 2 modules. Similarly, readers interested in using PDB2Graph will have to procure proprietary Matlab software or consider using older tools like NetworkAnalyzer or RINerator that rely on the open-source Cytospace platform [86].

Protein-protein interaction hotspots
Interactions of proteins with other proteins are fundamental characteristics of most biological processes such as substrate recognition, metabolism, signaling, pathogenic recognition, protein activation and inactivation. The involved residues have received much attention describing their potential for disrupting or enhancing activity, gaining knowledge about interactions or guiding protein structure prediction [95,96]. Interface residues were initially identified based on the distance between the interacting partners. Recently, approaches based on the Voronoi diagram tessellation have been used instead. Advantages of the Voronoi diagram include robust descriptions of the curvature and connectivity of these interfaces, allowing an unambiguous definition of the interface boundaries and providing a direct way to calculate the contact area [97,98]. Although several tools applying Voronoi tessellation to interface analysis have been published in recent years [97][98][99][100][101], many of them are obsolete or require a substantial degree of programming expertise.
As a user-friendly alternative focusing on protein-protein interactions, the PPI3D web-server was released [50]. This web-server uses a curated local database to retrieve and perform analyses. PPI3D assesses and differentiates three types of interactions: protein-protein, protein-peptide and domaindomain interactions using the weighted Voronoi tessellation implemented in Voronota [102], the correctness of which have been thoroughly tested on more than 90 000 structures from protein databank and almost 30 000 predicted protein structures of various qualities [102]. Initially, PPI3D uses Voronota to identify inter-atom contacts, which are later grouped into interresidue contacts. To execute the analysis, the web-server accepts three types of submissions: (i) a single sequence to predict protein and peptide-binding sites, (ii) two sequences to find all possible interactions and (iii) a PDB-ID code to identify all interactions in the entry. Importantly, the search for interfaces is not performed on the query protein alone but employs information on homologs too. The results are summarized in the form of a table containing links to a detailed description of interactions, such as interface residues, their contact area and type of interaction.
When searching for interface hotspots, the objective is typically to improve the binding affinity or avoid mutating those residues rather than destabilizing the protein complex. However, in some cases, disruption is desirable, which is the purpose of the DisruPPI software [51]. To obtain hotspot residues and disruptive mutations, DisruPPI assesses the stability of each monomer and the interactions between them thus enabling disruption of the binding but maintaining or improving the stability of the monomers ( Figure 3A). It is important to note that the interface region for analysis must be specified by the user, for which tools like PPI3D can be utilized. Whenever a query sequence is submitted, the software performs a search for homologs in a non-redundant database, the homologs are then aligned, and the interface residues are assessed based on conservation statistics. If the number of homolog sequences is not sufficient, structural modeling is employed to generate likely variants for assessment. In this case, the modeling starts by mutating each residue on the interface to the remaining amino acids, not present in an MSA. Mutations that destabilize the protein according to FoldX [103] or Rosetta [104] are discarded from the candidate pool. With the MSA of homologs, the software selects the most promising candidate sites and takes into consideration their conservation and hydrophobicity. Once candidates are selected, each is tested to obtain disruption and stability scores using a modified INT5 score as the metric [105]. Finally, DisruPPI identifies the lowest energy variant with the highest disruption score. The authors successfully employed the method in three experimental cases: the Hen Egg Lysozyme (HEL) with two anti-HEL antibodies, the HIV-1 glycoprotein gp120 with the cellular CD4 receptor and a red fluorescent protein (RFP) from Discosoma sp. These studies resulted in the following: (i) mutations disrupting anti-HEL antibody binding at levels of 27% and 59% of the WT, (ii) the most disruptive mutation of HIV-1-CD4 complex and (iii) five mutations disrupting RFP oligomerization while preserving stable monomers [51].

Interaction hotspots in ligand-transport tunnels
Residues governing the efficient exchange of cognate ligands between binding sites buried in protein cores and the bulk solvent have recently been recognized as potent hotspots for the engineering of a wide range of properties like activity, selectivity or stability [106]. However, understanding of ligand access and egress pathways is mostly reserved for methods based on molecular dynamics (MD) simulations [107]. (green and cyan cartoon) are analyzed to identify interaction hotspots (gold and pink sticks in the zoomed-in region). (B) Transport tunnels (blue, red and green objects) leading from the buried cavity of a target protein to the bulk solvent are delineated and explored by molecular docking of ligand to find protein residues that significantly contribute to the energy barrier of the ligand transport through the tunnel.
CaverDock software was developed as an alternative to expensive analyses of ligand transport using MD simulations [58,59]. It requires pre-computed tunnels for a given macromolecule by software such as CAVER 3.02 [108]. These tunnels are then systematically explored by docking the ligand with a modified AutoDock Vina tool [109]. To automate these calculations, CaverDock has been integrated into a Caver Web 1.0 server [60], which provides automatic calculation of tunnels for a protein of interest that is coupled with follow-up analyses of ligand transport through the user-selected tunnel ( Figure 3B). To start a Caver Web 1.0 calculation, the user selects a starting point for tunnel calculation based on detected pockets, the position of ligands, manual selection of residues or Cartesian coordinates. At this stage, detailed information on computed tunnels including their visualization can be accessed. Further, the user can specify tunnel(s) for the examination of transport of ligands, for which the molecule of interest must be uploaded in any format supported by the Open Babel tool [110], provided as an accession code to ZINC15 database [110], or drawn in an interactive window. As a result, the energy profile for each tunnel-compound pair is estimated. The presented functionality of the Caver Web 1.0 platform relating to ligand transport examination is a unique capability compared to alternative tools for tunnel detection [60]. Importantly, the applicability of CaverDock for protein engineering was verified by a detailed computational study of the transport of toxic pollutant 1,2,3trichloropropane in two variants of haloalkane dehalogenase featuring several advanced MD simulation methods [111]. For this model system, CaverDock was able to pinpoint similar hotspots as the sophisticated MD simulations, confirming its potential for the engineering of ligand transport [111].

Tools for predicting the effects of mutation
Once hotspots have been selected, the next step is to assess the effects of particular mutations at the site on the target properties of the protein. With this in mind, we have grouped tools based on the type of evaluated feature, starting with protein interactions, followed by flexibility and electrostatics. In this section, we present 10 tools that can be useful for prioritizing particular mutations for inclusion into a smart library based on their predicted effects. MutaBind [52], iSEE [53] and mCSM-PPI2 [54] can be used to estimate the effects of mutation on protein-protein interactions. mCSM-NA [55], PremPDI [56] and mCSM-lig [57] predict mutational effects on interactions of the protein with nucleic acids and ligands. The DynaMut server aims to rapidly evaluate changes in protein dynamics and stability after mutation [61]. The study of modifications in electrostatic potential upon mutation is the domain of the Mutantelec tool [62] and the analysis of electrostatic structures of proteins (AESOP) library [63]. Finally, the BioStructMap integrates data from various sources to help harness the knowledge for the next round of engineering [66]. The majority of the tools discussed in this section rely on machine-learning to derive predictive models from training datasets, while their quantitative performance is evaluated on the testing datasets. To compare performance of tools quantitatively, we summarized Pearson correlation coefficients (PCC) and root-mean-square errors (RMSE) achieved on these datasets as well as the main dataset parameters ( Table 2).

Effects of mutations on protein interactions
Proteins interact with each other, nucleic acids or small chemical compounds. As such, they are involved in all essential processes in any living cell. Recognizing the importance of protein-protein interactions, numerous tools have been developed to predict the impact of mutations in residues forming these interactions [124][125][126].
To expand a set of publicly available methods capable of quantitatively predicting the effects of single-point mutations on binding energy, the MutaBind web-server was developed [52]. This tool uses a consensus of multiple linear regression and a model trained with random forest (RF), both based on 1925 mutations of 80 protein-protein complexes from the SKEMPI database [112] to calculate the changes in binding free energy upon a single-point mutation ( G). The web-server relies on six physicochemical descriptors: Van der Waals interactions, polar solvation energies, unfolding free energies, solvent accessible surface area, conservation score and the ability of proline to introduce constraints on the protein backbone. The prediction performance of the tool was evaluated in a leave-one-complexout cross-validation, attaining a notable PCC of 0.68, and RMSE of 1.41 kcal.mol −1 [52]. Additionally, a separate model was provided for analyses of protease-inhibitor complexes, for which the authors observed improved performance (0.76 PCC and 1.48 kcal.mol −1 RMSE) on a subset of their training dataset (862 mutations on 16 complexes) [52]. Two de novo designed influenza inhibitors complexed with hemagglutinin used as targets (T55 and T56) in the 26th round of the CAPRI prediction experiment [120] with about 1000 mutations each were employed for independent evaluation. With PCC of 0.56 and 0.37 and RMSE of 2.58 and 4.27 kcal.mol −1 for the T55 and T56 targets, respectively [52], MutaBind scored better than three already well-established methods (BeAtMuSiC [126], FoldX [103] and molecular mechanics Poisson-Boltzmann surface area (MM/PBSA) [127]), showing appreciable improvement in PCC by 0.16 over the second bestperforming method. The relatively low predictive power might be caused by the use of de novo proteins forming an interface that is not commonly present in the training dataset of native proteins, the lack of respective crystal structures and the use of enrichment value as a measurement of binding affinity changes. Also, we would like to highlight that due to the limited number of highly stabilizing mutations present in the SKEMPI and CAPRI datasets, the prediction accuracy for such mutations cannot be adequately estimated and is presumably less reliable. MutaBind web-server can process up to 16 single-point mutations per run, delivering the binding affinity change, confidence level and structure model of each mutant. A similar approach has been adopted by the iSEE method [53], utilizing an RF approach to evaluate G at proteinprotein interfaces based on changes in structure-derived solvent accessible surface area, van der Waals, Coulomb and solvation energies and evolutionary information. A subset of the DACUM database [113], consisting of 1102 mutations from 57 protein dimers, was used for method training. In the authors hands, iSEE exhibited 0.80 PCC and RMSE of 1.41 kcal.mol −1 during 10-fold cross-validation [53]. Unfortunately, when using a blind dataset of 487 mutations from 56 protein complexes derived from the new SKEMPI 2.0 database [114], iSEE exhibited a PCC of only 0.25 and RMSE of 1.32 kcal.mol −1 , which is comparable to the performance of the other three state-of-the-art tools evaluated by the authors [53], i.e. FoldX, mCSM [128] and BindProfX [129]. The drop in the ranking performance observed for all four tools is in line with relatively small G values present in the dataset, making their evaluation rather sensitive to errors of experimental measurement [53]. Nonetheless, iSEE achieved good correspondence with experimental data on the classification of mutations in the MDM2-p53 complex. It is important to note that iSSE is provided as an R-model only and requires non-trivial data as its input, i.e. the 3D structure of both the WT and mutant complexes, eight energy terms, and evolutionary information in the form of position specific scoring matrix, which represent a marked barrier to widespread application of this tool.
The most recent addition to the toolbox for evaluation of the effects of mutations on protein-protein interaction is the mCSM-PPI2 web-server [54]. A characteristic feature of the mCSM-PPI2 is the use of a graph-based structural signature to represent the environment of the WT residue [128], where the residue environment is described as a graph, with the atoms as nodes and interactions among them as edges. In this scheme, atoms are assigned physicochemical-based pharmacophore types such as hydrophobic, positive, negative, acceptor, donor, aromatic, sulfur and neutral. Aside from the graph-based signatures, seven other types of features are employed in the model: pharmacophore changes due to mutation, structural and sequential residue environment in the WT protein, the nature of the mutant and WT residues, evolutionary information, non-covalent interaction network metrics, energetic terms and atomic fluctuations. Using these features, the mCSM-PPI2 server was trained using the ExtraTrees method on a derivative of the SKEMPI 2.0 database which consists of 4169 experimental variants and their hypothetical reverse mutations from 319 different complexes, giving a total of 8338 single-point mutants. During a leave-onecomplex-out cross-validation, the server achieved a PCC of 0.75 and RMSE of 1.30 kcal.mol −1 and showed an RMSE of 2.55 and 4.06 kcal.mol −1 for mutations of T55 and T56 targets from the 26th round of CAPRI competition [54]. Based on these results and Kendall scores, mCSM-PPI2 is on par with MutaBind and significantly better than the battery of 25 methods including FoldX, BeAtMuSiC, mCSM and MM/PBSA [54]. The mCSM-PPI2 server provides two modes of operation: (i) evaluating the effects of specific mutations defined by the user or (ii) assessing the mutation effects on the interface region by alanine scanning or saturation mutagenesis. For single-point mutations, the server presents the predicted change in binding affinity together with the visualization of the mutation in an interactive NGL viewer [130]. In the interface evaluation mode, an overview of all identified interfaces is provided, with each interface linked to a results page containing a summary with individual mutants listed and available for further exploration, including hotspot identification.
While there is no single study benchmarking all three reviewed tools for predicting the effect of mutations on protein-protein interactions, the results on the independent datasets from the CAPRI competition (Table 2) indicate that both MutaBind and mCSM-PPI2 are comparably accurate and among the best quantitative predictors available. Regarding iSEE performance, this tool reached similar performance to MutaBind in the recent study [131] performed on the testing dataset of iSEE ( Table 2). From the user point of view, mCSM-PPI2 offers a convenient option of automated mutation scanning of whole interfaces, and its results are presented in a far more interactive manner. On the contrary, the iSEE method requires the user to precompute all non-trivial input data manually. As the oldest among reviewed tools for proteinprotein interaction analysis, Mutabind is the most established in the community (Table 1 and see Supplementary Table 2 available online at https: //academic.oup.com/bib). Finally, our testing indicates that mCMS-PPI2 can perform the complete computation in the shortest time, followed by Mutabind (Table 1 and see Supplementary Table 1 available online at https: //academic.oup.com/bib).
Similarly to protein-protein interactions, interactions of proteins with nucleic acids are the basis of many crucial cellular processes such as replication, repair, recombination, transcription, translation and gene expression regulation. However, the development of rapid predictive methods has had much less success given a notably different physicochemical nature of these interactions, i.e. the prevalence of polar interactions that are often much harder to model precisely [132], and the limited availability of experimental structures as well as affinity data.
Benefiting from the recent release of high-quality data for protein-nucleic acid interactions in the second version of ProNIT database, the mCSM-NA web-server was developed. The tool uses the previously described graph-based approach, adopting additional descriptors for the atoms of nucleic acids, which are divided into three categories: a phosphate group, sugar and nitrogenous base [55]. From this representation, interactions between the protein and nucleic acid atoms were encoded as the graph-based structural signatures and served as an input for Gaussian process regression to train the predictive model. As a training dataset, the ProNIT database composed of 222 mutations from 28 complexes of protein-dsDNA, 42 mutations on six protein-ssDNA complexes and 67 mutations from five protein-RNA complexes was used [115]. On the whole training dataset, a PCC of 0.70 was achieved during 10-fold cross-validation, constituting a small improvement over its generalist predecessor, the mCSM tool, having PCC of 0.67 [128]. Interestingly, PCCs of 0.54, 0.85 and 0.75 were attained for dsDNA, ssDNA and RNA complexes when considered separately, respectively [55]. The authors also performed a blind test on 79 mutations from 14 protein-RNA complexes from the study by Barik and co-workers [121], in which mCSM-NA achieved a PCC of 0.56, raising to 0.68 when considering only mutations in the proximity of the RNA [55]. The user can submit up to 20 single-point mutations per run to mCSM-NA web-server to obtain the predicted change in binding affinity and the effect on protein stability. Importantly, predictions concerning protein-dsDNA interactions should be approached more cautiously, given their notably lower PCC on the training dataset.
An alternative approach to assess the effect of mutations on protein-DNA interactions was adopted by the PremPDI webserver [56], relying primarily on the interaction energy terms computed with the MM/PBSA method. The web-server uses the FoldX tool to model a structure of mutation, which is then, together with the WT structure, shortly energy minimized with NAMD software [133] using CHARMM36 force field [134]. Minimized complexes are then analyzed with the CHARMM package [135] to calculate differences in polar solvation, Van der Waals and electrostatic interaction energies, the number of hydrogen bonds, solvent accessible surface areas, which are further supplemented by information if the mutation occurs on the protein-DNA interface, the length of the protein chain and the pairwise statistical potential for protein folding obtained from the AAindex database [136]. Subsequently, the predictive model is trained by multiple linear regression on a PremPDI dataset, which was compiled from ProNIT and dbAMEPNI [116] databases and training dataset of SAMPDI tool [117], comprising 219 mutations from 49 protein-DNA complexes. A leave-one-complex-out cross-validation yielded PCC of 0.63 and RMSE of 0.95 kcal.mol −1 . The authors evaluated the performance of PremPDI, SAMPDI and mCSM-NA tools on subsets of training datasets overlapping among the pairs of tools, showing that on these datasets, Prem-PDI performs similarly to mCSM-NA and markedly outperforms SAMPDI. The PremPDI server provides an interactive calculation setup, enabling up to 16 single-point mutations to be specified, for which the prediction of the effect of the mutation and mutant 3D structure is generated.
Again, the two reviewed tools exhibited similar predictive performance for the effect of mutations on interactions of proteins with DNA. The main difference comes from the user experience and server applicability. Firstly, mCSM-NA can evaluate complexes containing RNA unlike PremPDI. The machinelearning predictions by mCSM-NA are several-fold faster than the demanding physics-based predictions implemented in PremPDI, especially for larger biomolecular systems (Table 1  and see Supplementary Table 1 available online at https: //academic.oup.com/bib). Additionally, the mCSM-NA web-server offers information on the effect of a mutation on protein stability and interactive online visualization. The overall advantages of the mCSM-NA tool are well reflected in its preferential use by the researchers (Table 1 and see Supplementary Table 2 available online at https: //academic.oup.com/bib).
Due to the importance of protein-ligand binding affinities to drug design and enzymology, their predictions have been the focus of numerous tools, including those employing deeplearning algorithms like DeepDTA [137] or K DEEP [138]. However, these methods have not been developed for predicting the impact of single-point mutations on binding affinity.
The recently assembled large-scale database of the effects of mutations on the affinity of protein-ligand complexes, Platinum [118], provided a foundation to develop the mCSM-lig web-server [57] capable of such prediction. This web-server is the last of the series of three tools based on the same graphbased approach that, aside from graph-based descriptors of interactions, employs physicochemical descriptors of the ligand, which are complemented by several other features, including a predicted change in stability, depth of evaluated residues, and experimental affinity of WT protein to the ligand. Two separate models for regression and classification of mutations were derived by using a Gaussian process for predictive regression and RF for classification tasks. For model training, a subset of the Platinum database was used, composed of 763 mutations from more than 200 protein-ligand complexes. During the leave-one-complex-out cross-validation on this dataset, the tool reached a PCC of 0.63 and RMSE of 2.06 kcal.mol −1 [57]. Moreover, the authors tested the capability of mCSM-lig to discriminate resistance-mutation profiles of different drugs binding to the same site. Evaluation of the binding of chemotherapeutics Imatinib, Nilotinib and Dasatinib to human ABL-kinase with this tool identified over 75% of resistance-causing mutations correctly [57]. In the same way, mutations affecting the binding of Efacirenz and Rilpivirine to HIV-1 reverse transcriptase were identified with more than 80% sensitivity [57]. Application of mCSM-lig to identify mutations improving the binding affinity of fluorescein was tested on the FluA protein. The tool achieved appreciable PCC of 0.69 with experimental data, indicating an aptness for protein engineering tasks [57]. The mCSM-lig webserver requires a structure of a protein-ligand complex, and the ligand affinity to the WT protein to deliver a prediction of the effect of mutation and the mutant 3D structure.

Effects of mutations on protein dynamics
Mutations can also affect the dynamics of proteins, altering their flexibility or rigidity [42,139,140]. The effect of mutations on protein dynamics can be studied with MD simulations, a powerful method that is computationally demanding and requires considerable expertise to perform. Therefore, quick and ready-to-use methods such as NMA have been developed as an alternative.
To facilitate the use of NMA for evaluating effects of mutation on protein dynamics and stability, DynaMut web-server [61] employs two NMA tools Bio3D [141] and ENCoM [142]. Protein dynamics is combined with the graph-based signature described for mCSM-PPI2, mCSM-NA and mCSM-lig tools to represent the WT structure, stability evaluation of individual conformations via DUET [143], and structural descriptors of the environment of the mutated residue like solvent accessible surface area, residue depth and secondary structure. Using the RF algorithm, the consensus predictor was trained on a dataset derived from the ProTherm database [119], containing 4594 experimental and hypothetical reverse mutations of 131 proteins. During 10-fold cross-validation on the training dataset, DynaMut achieved a PCC of 0.67 and RMSE of 1.31 kcal.mol −1 , performance was further confirmed with a blind dataset containing 702 forward and reverse mutations (PCC of 0.70 and RMSE of 1.45 kcal.mol −1 ) [61]. With two types of analysis, the web-server can be used to study the dynamic nature of a protein and the effect that point mutations have on its flexibility. DynaMut provides a detailed comparison of the interactions, flexibility and deformability of residues in the WT and mutant protein. 3D structure of the protein is supplemented by predicted energy change in DynaMut, and presented along with three other methods: mCSM, DUET [143] and SDM [144]. Here we would like to add that NMA is known to have, in some instances, a limited sensitivity to a mutation, unless it is responsible for a substantial conformational change [145].

Effects of mutations on protein electrostatic potential
The electrostatic potential of proteins can affect protein adsorption, ligand binding or thermodynamic stability of a protein complex [146][147][148][149][150]. Different software packages that calculate the electrostatic potential of proteins are available [151][152][153]. Apart from tools targeting protein stability such as the pStab web-server [154], there are not many tools that assess the effects of mutations on this vital feature.
To overcome this deficiency, the Mutantelec web-server was developed to enable evaluation of the effects that mutations have on electrostatic potential [62]. Additionally, the effects of phosphorylation of serine, threonine and tyrosine can be assessed, broadening the type of analysis that can be performed. The web-server uses Modeller [155] to optimize mutant structures, PDB2PQR [156] to assign atom parameters and APBS [151] to calculate electrostatic potential. As a result, the difference in electrostatics for each residue is represented as a histogram. Additionally, the web-server allows for download of 3D structures and electrostatic maps to visualize in PyMOL (Schrödinger, LLC., USA, http: //www.pymol.org/) or VMD [157]. The p53 protein was analyzed with Mutantelec to explain the effects of Arg249 mutations shown to inactivate the protein function [62].
Similar intentions have prompted a reimplementation of the AESOP framework [158] into a more accessible Python 2 library and web-server [63]. The AESOP library was shown to be applicable to study of the electrostatic similarity of protein families [159,160], perform an alanine scanning of ionizable amino acids to identify possible electrostatic hotspots [161,162], and assess the effects of single-point mutations on the free energy of association [163,164]. Analogously to the Mutantelect workflow, AESOP makes use of Modeller [165] to generate and optimize mutant structures and PDB2PQR and APBS software packages to calculate electrostatic potentials. The capabilities of this webserver are currently limited to alanine scanning. As a result, the AESOP library and web-server generate PDB files of mutants, information on the energy change caused by the mutations, and several predefined structural visualizations.

Data integration
Engineering efforts do not inevitably end with characterizing properties of the altered protein since information about the mutants can be stored and further reused. This way, an additional cycle of design can be performed iteratively, combining experimental and computational stages. For more convenient analyses, data coming from different sources can be mapped onto the 3D structure of proteins, and there are web-services capable of mapping conservation [166], coevolution [167], biochemical and biomedical annotations [168]. All those services are, however, limited to displaying only one property at a time.
Hence, the BioStructMap tool was developed [66], allowing users to map any sequence-associated function that returns a numeric value onto a 3D protein structure. The BioStructMap package includes pre-defined functions to analyze data, such as mapping polymorphic hotspots, amino acid propensity scales, Tajima's D index, nucleoside diversity and customized data aggregation. As an input, the package requires a sequence alignment, a PDB file, and a reference sequence matching the sequence of the PDB. The output includes residue-values, as a Python dictionary, which can be mapped to the PDB file on the B-factor column or returned as a text file. Although simple, the most attractive feature of the BioStructMap package is the capability of customization, the only requirement is that the function should return a numerical value. For instance, the user could use numerical results coming from other tools like MutaBind, AESOP or alanine scanning and map those values onto the 3D structure. The resulting PDB file can be analyzed with PyMOL to detect hotspots for the mapped property. BioStructMap is also offered as a web-server that supports analyses with the pre-defined functions.

Integrative platforms for protein engineering
The typical engineering workflow often covers a range of structural and evolutionary analyses, predictions of the effects of mutation on protein function and stability, and their visualization. Execution of such analysis requires bioengineers to search for suitable computational tools and transfer data from one tool to another, which is frequently complicated by incompatibilities among these tools. To relieve bioengineers of such hurdles, integrative services gather different methodologies for hotspot identification and mutation analyses and streamline the flow of data into automated workflows.
The HotSpot Wizard web-server is a platform developed with the primary purpose of finding protein engineering hotspots by combining structural and sequential analyses [169]. The second version of the service notably expanded the original focus on functional hotspots by providing access to the other three design strategies: stability by structural flexibility, stability by sequence consensus or correlated hotspots [64]. Functional hotspots are identified based on the protein active site and transport tunnels. At the same time, highly conserved residues that are often indispensable for catalysis are avoided in order not to compromise enzyme function. This method of hotspot selection is recommended when substrate specificity is the target property. Regarding stabilizing hotspots, highly mobile residues are identified by their thermal B-factors. Alternatively, an MSA of homologous proteins can be analyzed to find residues differing from the prevalent variant at any given position in homologous proteins. Both strategies were shown to be useful when a protein with improved thermal stability is desired [170,171]. Lastly, the correlated hotspot approach is based on the MSA that is scrutinized by seven methods for coordinated changes in the sequence and combined to a Z-score based consensus. Considering the vast number of identified correlated hotspots, further analyses using alternative approaches are advisable to obtain more focused predictions, unless the modulation of allostery is the engineering target. The main innovation of the current third version of the platform was the integration of an automated protein structure modeling workflow to overcome the applicability limit of previous versions [65]. Additionally, the Rosetta and FoldX tools can now be utilized for the evaluation of the effect of particular mutations at hotspots of interest. As the input, a protein sequence, PDB file or PDB-ID code is required. Additionally, some of the most influential parameters controlling calculations of integrated tools can be modified in advanced settings. Once the calculations are finished, results for each strategy can be accessed from the navigation panel. Further, the effects of selected mutations can be evaluated, used to develop an optimized smart library and generate nucleotide sequences employing the codon usage of a target organism. Overall, the HotSpot Wizard is a user-friendly and well-established service for hotspot identification and library design that has been available for over 10 years and is kept up to date, with the latest version published in 2018. It has been successfully verified in seven engineering studies reported so far (see Supplementary Table 3 available online at https: //academic.oup.com/bib).

Conclusions, perspectives and challenges
The reviewed computational tools can provide valuable guidance for protein engineering toward tailoring a wide range of target properties. Notably, the applicability of the protein engineering toolbox has been expanded from targeting traditionally engineered regions of proteins, such as binding sites and transport-pathways, to addressing even more elusive targets including protein dynamics and allostery. Here, most tools resorted to approximating the dynamics by rapid methods based on NMA to minimize overall execution time and required computing resources. However, we expect that other approaches including the geometric-constraints-based method, tCONCOORD [172], or the perturbation-based methods, L-RIP or RIPlig [173], will be integrated as an alternative or used alongside NMA to overcome the intrinsic limitations of each method and to provide complementary insights into the structure-dynamics-function relationships of engineered proteins.
There is a clear trend of protein engineering software becoming more user-friendly, commonly featuring an attractive webserver interface. However, we have noted only a handful of integrative efforts to provide meta-servers where the user can employ different approaches to protein design. Even the current one-stop-shop platforms like HotSpot Wizard do not yet offer comprehensive approaches. In particular, they do not allow analysis of possible interactions of the target protein with other molecules, leaving such burdensome tasks to their users. Therefore, an essential step toward complete workflows is to integrate tools for the prediction of macromolecular binding sites and evaluation of actual protein interactions with cognate binding partners, ultimately creating complete pipelines which will be even more attractive for protein engineers.
Since a significant fraction of the reviewed tools employs various methods of supervised machine learning, the quality of the available datasets for development and testing is of utmost importance. We have seen an example of the role that a more extensive dataset can play for the accuracy of the predictive methods for the effect of a mutation on proteinprotein interactions. In that case, the release of the SKEMPI 2.0 database, which more than doubled the number of annotated mutations and almost tripled the diversity in protein-protein complexes, enabled the derivation of more accurate predictive models adopted by the mCSA-PPI2 tool. However, even this database is notably biased by the strategies and research objectives traditionally applied in the protein engineering community [114], for example, there is a prevalence of singlepoint mutations (most frequently to alanines), the majority of mutations being located at the interface (mostly at its core), and the distribution of observed effects of mutations being shifted toward destabilization. This lack of proportionality hampers the development of more accurate and general predictive methods. The impact of such biases on the development of machine learning models for protein engineering and some necessary steps toward resolving these issues have been discussed elsewhere [174].
Another aspect related to the experimental datasets is their sparsity. For most proteins, information on the effects of mutations is available for a few sites only, and even for those that were mutated, often not all variants were tested. Such situations markedly complicate the performance evaluation of tools for hotspot identification and mutation prioritization. The data sparsity also makes any rigorous comparison among different methods almost impossible since experimental data rarely cover the top hits proposed by individual tools. As a remedy, we suggest turning to datasets from systematic mutagenesis projects like the dataset of 13 massively mutated proteins compiled from the literature and patents [175] and, more recently, datasets originating from deep mutation scanning techniques that provide almost complete mapping of the mutation landscape [176,177]. Also, we would like to highlight another initiative, ProtaBank repository [178], which introduces a standardized format for storing and reporting data from protein engineering studies to facilitate accurate comparisons among them. Its utilization has the potential to overcome some of the challenges discussed here.

Key Points
• Computationally supported rational and semirational protein engineering constitute widely accepted and efficient approaches. These approaches become more accessible to the broad scientific community due to the development of user-friendly tools, 18 of which are reviewed here.
• These novel tools enable identification of engineering hotspots promising for the modification of protein function based on allosteric communication, interactions with other proteins or ligands.
• Mutations can be evaluated for their impact on protein interactions with other macromolecules, protein dynamics or an electrostatic potential that can be estimated by user-friendly software covered in this review.
• Identification of hotspot residues and prediction of the effect of their mutation can be easily performed via integrative platforms like HotSpot Wizard. This platform was successfully applied in multiple design experiments leading to the identification of enhanced variants showcased here.

Supplementary data
Supplementary data are available at Briefings in Bioinformatics online.