The electrostatic profile of consecutive Cβ atoms applied to protein structure quality assessment

The structure of a protein provides insight into its physiological interactions with other components of the cellular soup. Methods that predict putative structures from sequences typically yield multiple, closely-ranked possibilities. A critical component in the process is the model quality assessing program (MQAP), which selects the best candidate from this pool of structures. Here, we present a novel MQAP based on the physical properties of sidechain atoms. We propose a method for assessing the quality of protein structures based on the electrostatic potential difference (EPD) of Cβ atoms in consecutive residues. We demonstrate that the EPDs of Cβ atoms on consecutive residues provide unique signatures of the amino acid types. The EPD of Cβ atoms are learnt from a set of 1000 non-homologous protein structures with a resolution cuto of 1.6 Å obtained from the PISCES database. Based on the Boltzmann hypothesis that lower energy conformations are proportionately sampled more, and on Annsen's thermodynamic hypothesis that the native structure of a protein is the minimum free energy state, we hypothesize that the deviation of observed EPD values from the mean values obtained in the learning phase is minimized in the native structure. We achieved an average specificity of 0.91, 0.94 and 0.93 on hg_structal, 4state_reduced and ig_structal decoy sets, respectively, taken from the Decoys `R' Us database. The source code and manual is made available at https://github.com/sanchak/mqap and permanently available on 10.5281/zenodo.7134.

The challenge of deriving the native structure of a protein from its sequence has intrigued researchers for decades 1 . Methods that predict putative structures from sequences are based either on features from databases of known structures (template-based methods) 2-4 or use first principles of atomic interactions (ab initio or de novo methods) 5-7 . Typically, these methods yield multiple, closely-ranked possibilities. Model quality assessment programs (MQAP) that validate accuracy of these predicted structures are used to select the best candidate from the set of predicted structures.
MQAPs can be classified as energy, consensus or knowledge based. Two major sources of errors in energy based methods used for refining or discriminating protein structures are inaccuracies in the force field due to the inherent approximations in equations that model multi-atomic configurations, and inadequate sampling of the conformational space [8][9][10][11][12] . Consensus based methods are based on the principle that structural features that are frequently observed in a population of structures are more likely to be present in the native structure [13][14][15][16] . These clustering methods outperform other MQAP methods 14 and are "very useful for structural meta-predictors 17 ". However, they are prone to be computationally intensive due structure-to-structure comparison of all models 16 , and are of limited use when the number of possible structures is small 18 . Knowledge based methods proceed by deriving an empirical potential (also known as statistical potential) from the frequency of residue contacts in the known structures of native proteins 19,20 . For a system in thermodynamic equilibrium, statistical physics hypothesizes that the accessible states are populated with a frequency which depends on the free energy of the state and is given by the Boltzmann distribution. The Boltzmann hypothesis states that if the database of known native protein structures is assumed to be a statistical system in thermodynamic equilibrium, specific structural features would be populated based on the free energy of the protein conformational state. Applying a converse logic, Sippl reasoned that the frequencies of occurrence of structural features such as interatomic distances in the database of known protein structures could be used to assign a free energy (potential of mean force) for a given protein conformation 21,22 . Furthermore, this statistical potential can be used to discriminate the native structure 23-27 . The proper characterization of the reference state is a critical aspect in applying statistical potentials 23 . In spite of their popularity, the application of such empirical energy functions to predict and assess protein structures are vigorously debated 28,29 . Many MQAP programs perform better when multiple statistical metrics are combined [30][31][32][33] . The paramount importance of obtaining high quality protein structures from sequences using in silico methods can be estimated by the effort invested by researchers every two years 34 to evaluate both structure prediction tools 35 and MQAPs 17, 34,36 .
Here, we propose a novel statistical potential to assess the quality of protein structures based on the electrostatic potential difference (EPD) of Cβ atoms in consecutive residues -EPD profile of sidechain atoms used in assessment of protein structures (ESCAPIST). Previously, we have established that the EPD is conserved in cognate pairs of active site residues in proteins with the same function [37][38][39][40] . The ability of finite difference methods to quickly obtain consistent electrostatic properties from peptide structures provides an invaluable tool for investigating other innate properties of protein structures 41 . We plot the EPD profiles for different atom types (Cα atoms, Cβ atoms and the C-N bond) in consecutive residues from a set of non-homologous protein structures obtained from the PISCES database (http://dunbrack.fccc.edu/PISCES.php) 42 . We proceed to show that the EPD between Cβ atoms in consecutive residues can be used to generate a scoring function that assesses the quality of protein structures. This EPD scoring function is then applied to standard decoy sets from the Decoys 'R' Us database (http://dd.compbio.washington.edu) to establish the validity of our method 43 .

Electrostatic potential difference (EPD) based discrimination
To extract feature values we chose a set of 1000 proteins from the PISCES database with percentage identity cutoff of 20%, resolution cutoff of 1.6 Å and a R-factor cutoff of 0.25 (SI Table 1).
Invariance of the EPD in the C-N peptide bond and between Cα atoms of consecutive residues Adaptive Poisson-Boltzmann Solve (APBS) writes out the electrostatic potential in dimensionless units of kT/e where k is Boltzmann's constant, T is the temperature in K and e is the charge of an electron. The units of EPD are same as that of the electrostatic potential. The EPD of the C-N peptide bond has a Gaussian distribution with mean = 420 EPD units and SD = 55 EPD units ( Figure 1). In the probability distribution for four pairs of amino acids the mean of all pairs of amino acids are the same (Figure 1a). Figure 1b shows the scatter plot for the mean and standard deviation (SD). Thus, the amino acids are indistinguishable using the profile of the EPD of the C-N peptide bond across all protein structures since they have identical mean values and a large variance (SD=~50).
The probability distribution for four pairs of amino acids for the EPD between the Cα atoms of consecutive residues (Figure 2a) have means that are slightly more varied than those for the C-N bond (Figure 1a). In the scatter plot for the mean and SD of all pairs ( Figure 2b) the outliers are pairs that include proline, which have a higher mean, although the magnitude of SD is the same ( Table 1).
Distinctive EPD between Cβ atoms of consecutive residues for certain amino acid pairs In contrast to the results described above, the EPD between the Cβ atoms in consecutive residues in the peptide structure can be used to discriminate different amino acid pairs in the protein structure. The mean EPD of all amino acid pairs are much more varied (Figure 3a). These pairs do not include glycine, which lacks a sidechain. In the scatter plot for the mean and SD, the outliers are pairs that include

Amendments from Version 2
We have modified our discussion to reflect the point we had missed previously regarding the pKa values in accordance with Dr Kammerlin's suggestion.  cysteine (Figure 3b), which have a much higher SD (=~90) as compared to other pairs (SD=~35) ( Table 2), and thus cannot be used for discriminatory purposes.

REVISED
These values are used as a discriminator when choosing the native structure from a set of possible candidates (Table 3). To establish the non-triviality of these values, we also show that the variance of the EPD between these pairs increases with increasing sequence distance. Thus, the EPD between the pairs 'DF' and 'HS' has lesser correlation as the sequence distance between them increases (sample size for each sequence distance is > 30) ( Figure 4). The SD for distance 1 (i.e. consecutive residues) is 29.8 EPD units and 31.8 EPD units for 'DF' and 'HS', respectively -and rises to around 60 EPD units with increasing sequence distance.

Validating using decoy sets
We obtained the score (PDScore) of any given protein structure by comparing the electrostatics of the Cβ atoms based on Table 3. To benchmark model quality assessment programs, we used decoy sets from the Decoys 'R' Us database 43 . We detail our results from some of these decoy sets. Each set has several structures that are supposed to be ranked worse than the native structure. It is seen that the mean is much more varied than the electrostatic potential difference (EPD) for Cα and the C-N peptide bond.  The misfold decoy set has ~20 protein structures, each of which has a correct and an incorrect structure specified (three structures have two incorrect structures: we randomly chose the first) 44 . The PDScore of these proteins were computed (Table 4). Barring three structures (PDBids: 1CBH, 1FDX and 2SSI), the PDScore of the incorrect structure is higher than that of the correct structures.
The hg_structal set has about ~30 proteins. Each protein has 30 structures (including the native structure). Table 5A shows specificity obtained for structures in this decoy set. The average specificity obtained for this decoy set is 0.91 (Table 5A). The decoy set 4state_reduced has ~600 structures for each of the seven proteins. We obtain an average specificity of 0.94 for this decoy set (Table 5B). Similarly, for the ig_structal decoy set we obtain a specificity of 0.93 (Supplementary Table 1).

Discussion
The functional characterization of a protein from its sequence using in silico methods based on the 'sequence to structure to function' paradigm is contingent upon the availability of its 3D-structure. The rapidly developing field of next generation sequencing has exacerbated the bottleneck of obtaining structural data using crystallization techniques 45 . This ever-widening gap has been filled by methods that predict structures from sequences 46 , based either on features from databases of known structures 2-4 or from first principles of atomic interactions 5,6 .
The various sources of error in protein structure prediction have been previously discussed in detail 47 . An incorrect model of a protein structure will result in an inaccurate analysis of its properties 48 . For example, continuum models 49 that compute potential differences and pK a values from charge interactions in proteins 50 are sensitive to the spatial arrangement of the atoms in the structure.
It must be pointed out that other detailed methods using quantum mechanical/molecular mechanical (QM/MM) techniques and doing extensive conformational sampling have been able to determine the side chain pK a values with high accuracy 51 . Accurate structural information is indispensable for in silico methods that extract the electrostatic profile of atoms in the peptide structure 41,52 , and for other methods widely used in pharmacology for drug discovery 53 . Model quality assessment programs (MQAP) that validate the accuracy of predicted structures are thus a critical aspect in the process of modeling a protein structure from its sequence. MQAPs can be classified as energy 8-12 , consensus 13-16 or knowledge based (statistical potential) 21-27 . The state of the art methods for predicting structures 35 and MQAPs 17,34,36 are evaluated by researchers every two years.
Previously, we hypothesized and demonstrated that the electrostatic potential difference (EPD) in cognate pairs in the active site are conserved in proteins with the same functionality 37,40,54 , even when evolution has converged to the same catalytic from completely different sequences 55 . This similarity is observed in structures solved independently over many years and establishes the reliability of the APBS and PDB2PQR implementations 41,56 . We focused on unraveling other electrostatic features that are innate to protein structures.
Here, we first demonstrate that the EPD between the C-N peptide bond and the Cα atoms of consecutive residues are independent of the amino acid type. This is expected, since the distance between these atoms are almost invariant across all structures. The EPD of the C-N bond has a high variance, implying that the backbone accommodates relatively large variations while seeking energetically minimized structures.
The true source of the chemical and structural diversity in protein structures is the side chain atoms. Every amino acid, except glycine, has a Cβ atom that hosts a unique moiety of atoms. Although the reactive groups are different for amino acids, we show that this difference is encapsulated in the backbone Cβ atoms. We first show that different pairs of amino acids have significantly different mean EPD values in side chain Cβ atoms (Figure 3), unlike the EPD of the C-N peptide bond (Figure 1) or the EPD between consecutive Cα atoms ( Figure 2). Further, the variance is much less than in the EPD of the C-N bond. These facts suggested that the EPD between Cβ atoms of consecutive residues in the peptide structure might act as a discriminator. Our hypothesis is based on the insightful Boltzmann law that lower energy conformations are disproportionately sampled, on the thermodynamic hypothesis 57 that the native structure has minimal energy, and the hypothesis that statistical derived features in known protein structures have a Gaussian distribution 21 . We apply our discriminator to standard decoy sets from the Decoys 'R' Us database to establish the validity of the method 43 .
Our work also highlights the unique properties of proline in the protein structure 58 . This is evident from the different magnitude of EPD in consecutive Cα atoms involving proline (Table 1). Another noteworthy aspect is the large variation in EPD in consecutive Cβ atoms involving cysteine (Table 2), demonstrating the unique role cysteine plays in providing flexibility to protein structures, a critical element in the evolution of complex organisms 59 . The discrimination of Cβ atoms also provides a uniform basis for methods that require a single-atom representation of a residue. Such methods depend on a correct parameterization of the reactive atoms 37 , a task further complicated by amino acids such as histidine which has two reactive atoms. For example, the EPD between the negatively charged E and D with respect to the aromatic phenylalanine is -108 and -93 EPD units, in spite of the difference in their reactive atom. Similarly, the EPD between alanine and the three aromatic amino acids (F, W and Y) are -67, -66 and -63 EPD units respectively.
We achieved an average specificity of 0.91, 0.94 and 0.93 on hg_structal, 4state_reduced and ig_structal decoy sets, respectively, taken from the Decoys 'R' Us database. We have previously demonstrated that the fisa decoy set can be discriminated by a distance based discriminator 60 . ESCAPIST does not discriminate the native structure in this decoy set (Table 5C). The physical implication of ESCAPIST results on the fisa decoy set, which has significant RMSD for backbone Cα atoms, needs further elaboration. The input to a finite difference Poisson-Boltzmann (FDPB) analysis is a charge distribution that might be unfeasible due to energy functions other than electrostatics. For example, van der Waals force or the elastic bond length force components might prevent two atoms from being in close proximity. However, if such a physically impossible configuration were presented to a FDPB-based analysis tool, such as APBS 41 , it would still generate an electrostatic potential landscape. Inferences based on this potential landscape would be incorrect due to its physical non-viability. Thus, before invoking the EPD constraints specified here, it is imperative that other spatial constraints that are rarely violated in structures are checked. Possibly for this same reason, MQAPs that combine many features in their scoring functions are superior. Moreover, it should be kept in mind that decoy sets, like most benchmarking sets, are prone to biases 61 and possible errors 31 . In fact, the fisa decoy set has been shown to violate the van der Waals term 61 . To summarize, we present a novel discriminating feature in protein structures based on the electrostatic properties of the side chain atoms. We validated this discrimination in several decoy sets taken from the Decoys 'R' Us database, and achieved high specificities in most decoy sets.

Methods
Our proposed method has two phases. In the learning phase, we process multiple structures to extract the feature values -mean values The EPD between a pair of amino acid is order-independent -for example, the EPD statistics for the pair 'AC' (alanine-cysteine) includes the EPD of both 'AC' and 'CA' (with the sign reversed).
Quality assessment phase Algorithm 2 shows the function AssessEPDQuality() that generates the PDscore for a given protein from the template file generated by the learning phase. The set of proteins Φ AssessmentPhase proteins consists of the native structure P 1 and N-1 decoys structures (Equation 4). Once again, barring x=IgnoreNTerm and y=IgnoreCTerm number of residues from the N and C terminals, the pairwise EPD for consecutive residues are computed. The absolute value of the difference of these values from their corresponding means, if they exist, in the template file is added to generate the absolute score (Equation 5). This score is normalized with the number of residues that have been compared to obtain the final PDscore. In summary, the PDscore encapsulates the average distance of the EPD for the given atom pairs (it may be Cβ, Cα or the C-N bond) of consecutive residues from their mean values. We hypothesize that in the native or a near native structure, the PDscore will be minimized for the EPD of Cβ atoms of consecutive residues, i.e. given a set of proteins Φ proteins consists of the native structure P 1 and N-1 decoys structures, P 1 will have the minimum PDscore (Equation 6). of electrostatic potential difference (EPD) for each amino acid pair. These feature values are applied on query proteins to obtain a score (PDscore) that indicates the deviation of the feature values in the given structure from the 'ideal' values. Thus, a lower PDscore indicates a better candidate.
Learning phase Algorithm 1 shows the procedure LearnFeatures() that extracts the feature values from a set of proteins Φ LearningPhase proteins (Equation 1). We ignore the first x=IgnoreNTerm and last y=IgnoreCTerm pairs of residues in the protein structure to exclude the terminals. For every consecutive pair of residues in the structure, we calculate the EPD (see below for method) between two specified atoms (atomP and atomQ). Both atomP and atomQ are set to Cβ to obtain EPD values for Cβ atoms, while we set atomP to 'C' and atomQ to 'N' in order to obtain the C-N peptide bond EPD values. The mean (Mean learnt value -MLV) and standard deviation (SD) are computed for each amino acid type pair (AAType) in protein (Equation 2), and the mean computed for the globals set of proteins (MLV(AAType x , AAType y )) for each pair of amino acid types (Equation 3). Pairs whose EPD have a SD greater than a threshold value (sdThresh, 50 by default) are ignored. Means for all significant pairs (φ pairMean ) are noted to a file, which is the input to the quality assessment phase.    The quality of protein structure models is assessed by the geometry of adjacent C beta atoms. The approach successfully distinguishes properly folded proteins in most cases. It adds a new way to assess protein models and could be included in the protein model assessment toolbox.
Just a few comments: a) Table 4 lists three of 20 structures where the incorrect one has a lower score. A few comments about structural features of those examples that lead unexpected scores would be useful.
b) It might be preferable to note in the title that the method is applied to computational models of protein structure as a way to distinguish the manuscript from those that deal with quality assessment of experimentally determined structures. c) I did not test the available source code.
I have read this submission. I believe that I have an appropriate level of expertise to confirm that it is of an acceptable scientific standard.
No competing interests were disclosed. We would like to thank you for taking the time and reviewing our paper, and deeply appreciate your positive comments. We do not have the expertise to comment on the structural characteristics that lead to the 3 negative results in Table 4 -these require insights into crystal structures that we do not possess at the present time. Also, we believe that this method could be applied to any structure -computational or experimentally obtained -and thus are leaving the title unchanged.
No competing interests were disclosed. Competing Interests: