Ab initio gene prediction for protein-coding regions.

Motivation: Ab initio gene prediction in nonmodel organisms is a difficult task. While many ab initio methods have been developed, their average accuracy over long segments of a genome, and especially when assessed over a wide range of species, generally yields results with sensitivity and specificity levels in the low 60% range. A common weakness of most methods is the tendency to learn patterns that are species-specific to varying degrees. The need exists for methods to extract genetic features that can distinguish coding and noncoding regions that are not sensitive to specific organism characteristics. Results: A new method based on a neural network (NN) that uses a collection of sensors to create input features is presented. It is shown that accurate predictions are achieved even when trained on organisms that are significantly different phylogenetically than test organisms. A consensus prediction algorithm for a CoDing Sequence (CDS) is subsequently applied to the first nucleotide level of NN predictions that boosts accuracy through a data-driven procedure that optimizes a CDS/non-CDS threshold. An aggregate accuracy benchmark at the nucleotide level shows that this new approach performs better than existing ab initio methods, while requiring significantly less training data. Availability and implementation: https://github.com/BioMolecularPhysicsGroup-UNCC/MachineLearning.

Protein Function Analysis through Machine Learning.

Machine learning (ML) has been an important arsenal in computational biology used to elucidate protein function for decades. With the recent burgeoning of novel ML methods and applications, new ML approaches have been incorporated into many areas of computational biology dealing with protein function. We examine how ML has been integrated into a wide range of computational models to improve prediction accuracy and gain a better understanding of protein function. The applications discussed are protein structure prediction, protein engineering using sequence modifications to achieve stability and druggability characteristics, molecular docking in terms of protein-ligand binding, including allosteric effects, protein-protein interactions and protein-centric drug discovery. To quantify the mechanisms underlying protein function, a holistic approach that takes structure, flexibility, stability, and dynamics into account is required, as these aspects become inseparable through their interdependence. Another key component of protein function is conformational dynamics, which often manifest as protein kinetics. Computational methods that use ML to generate representative conformational ensembles and quantify differences in conformational ensembles important for function are included in this review. Future opportunities are highlighted for each of these topics.

Functional Dynamics of Substrate Recognition in TEM Beta-Lactamase.

The beta-lactamase enzyme provides effective resistance to beta-lactam antibiotics due to substrate recognition controlled by point mutations. Recently, extended-spectrum and inhibitor-resistant mutants have become a global health problem. Here, the functional dynamics that control substrate recognition in TEM beta-lactamase are investigated using all-atom molecular dynamics simulations. Comparisons are made between wild-type TEM-1 and TEM-2 and the extended-spectrum mutants TEM-10 and TEM-52, both in apo form and in complex with four different antibiotics (ampicillin, amoxicillin, cefotaxime and ceftazidime). Dynamic allostery is predicted based on a quasi-harmonic normal mode analysis using a perturbation scan. An allosteric mechanism known to inhibit enzymatic function in TEM beta-lactamase is identified, along with other allosteric binding targets. Mechanisms for substrate recognition are elucidated using multivariate comparative analysis of molecular dynamics trajectories to identify changes in dynamics resulting from point mutations and ligand binding, and the conserved dynamics, which are functionally important, are extracted as well. The results suggest that the H10-H11 loop (residues 214-221) is a secondary anchor for larger extended spectrum ligands, while the H9-H10 loop (residues 194-202) is distal from the active site and stabilizes the protein against structural changes. These secondary non-catalytically-active loops offer attractive targets for novel noncompetitive inhibitors of TEM beta-lactamase.

Victory Tax: A Holistic Income Tax System.

How can an income tax system be designed to exploit human nature and a free market to create a poverty free society, while balancing budgets without disproportional tax burdens? Such a tax system, with universal character, is deduced from the following guiding principles: (1) a single tax rate applies to all income types and levels; (2) the tax rate adjusts to satisfy budget projections; (3) government transfer only supplements the income of households with self-generated income below the poverty line; (4) deductions for basic living expenses, itemized investments and capital losses are allowed; (5) deductions cannot be applied to government transfer. A general framework emerges with three parameters that determine a minimum allowed tax deduction, a maximum allowed itemized deduction, and a maximum deduction defined by income percentage. An income distribution that mimics the United States, and a series of log-normal distributions are considered to quantitatively compare detailed characteristics of this tax system to progressive and flat tax systems. To minimize government dependency while maximizing after-tax income, the effective tax rate (ETR) as a function of income percentile takes the shape of the letter, V, inspiring the name victory tax, where the middle class has the lowest ETR.

JEDi: java essential dynamics inspector - a molecular trajectory analysis toolkit.

BACKGROUND: Principal component analysis (PCA) is commonly applied to the atomic trajectories of biopolymers to extract essential dynamics that describe biologically relevant motions. Although application of PCA is straightforward, specialized software to facilitate workflows and analysis of molecular dynamics simulation data to fully harness the power of PCA is lacking. The Java Essential Dynamics inspector (JEDi) software is a major upgrade from the previous JED software. RESULTS: Employing multi-threading, JEDi features a user-friendly interface to control rapid workflows for interrogating conformational motions of biopolymers at various spatial resolutions and within subregions, including multiple chain proteins. JEDi has options for Cartesian-based coordinates (cPCA) and internal distance pair coordinates (dpPCA) to construct covariance (Q), correlation (R), and partial correlation (P) matrices. Shrinkage and outlier thresholding are implemented for the accurate estimation of covariance. The effect of rare events is quantified using outlier and inlier filters. Applying sparsity thresholds in statistical models identifies latent correlated motions. Within a hierarchical approach, small-scale atomic motion is first calculated with a separate local cPCA calculation per residue to obtain eigenresidues. Then PCA on the eigenresidues yields rapid and accurate description of large-scale motions. Local cPCA on all residue pairs creates a map of all residue-residue dynamical couplings. Additionally, kernel PCA is implemented. JEDi output gives high quality PNG images by default, with options for text files that include aligned coordinates, several metrics that quantify mobility, PCA modes with their eigenvalues, and displacement vector projections onto the top principal modes. JEDi provides PyMol scripts together with PDB files to visualize individual cPCA modes and the essential dynamics occurring within user-selected time scales. Subspace comparisons performed on the most relevant eigenvectors using several statistical metrics quantify similarity/overlap of high dimensional vector spaces. Free energy landscapes are available for both cPCA and dpPCA. CONCLUSION: JEDi is a convenient toolkit that applies best practices in multivariate statistics for comparative studies on the essential dynamics of similar biopolymers. JEDi helps identify functional mechanisms through many integrated tools and visual aids for inspecting and quantifying similarity/differences in mobility and dynamic correlations.

Molecular function recognition by supervised projection pursuit machine learning.

Identifying mechanisms that control molecular function is a significant challenge in pharmaceutical science and molecular engineering. Here, we present a novel projection pursuit recurrent neural network to identify functional mechanisms in the context of iterative supervised machine learning for discovery-based design optimization. Molecular function recognition is achieved by pairing experiments that categorize systems with digital twin molecular dynamics simulations to generate working hypotheses. Feature extraction decomposes emergent properties of a system into a complete set of basis vectors. Feature selection requires signal-to-noise, statistical significance, and clustering quality to concurrently surpass acceptance levels. Formulated as a multivariate description of differences and similarities between systems, the data-driven working hypothesis is refined by analyzing new systems prioritized by a discovery-likelihood. Utility and generality are demonstrated on several benchmarks, including the elucidation of antibiotic resistance in TEM-52 beta-lactamase. The software is freely available, enabling turnkey analysis of massive data streams found in computational biology and material science.

Rigidity and flexibility characteristics of DD[E/D]-transposases Mos1 and Sleeping Beauty.

DD[E/D]-transposases catalyze the multistep reaction of cut-and-paste DNA transposition. Structurally, several DD[E/D]-transposases have been characterized, revealing a multi-domain structure with the catalytic domain possessing the RNase H-like structural motif that brings three catalytic residues (D, D, and E or D) into close proximity for the catalysis. However, the dynamic behavior of DD[E/D]-transposases during transposition remains poorly understood. Here, we analyze the rigidity and flexibility characteristics of two representative DD[E/D]-transposases Mos1 and Sleeping Beauty (SB) using the minimal distance constraint model (mDCM). We find that the catalytic domain of both transposases is globally rigid, with the notable exception of the clamp loop being flexible in the DNA-unbound form. Within this globally rigid structure, the central ß-sheet of the RNase H-like motif is much less rigid in comparison to its surrounding α-helices, forming a cage-like structure. The comparison of the original SB transposase to its hyperactive version SB100X reveals the region where the change in flexibility/rigidity correlates with increased activity. This region is found to be within the RNase H-like structural motif and comprise the loop leading from beta-strand B3 to helix H1, helices H1 and H2, which are located on the same side of the central beta-sheet, and the loop between helix H3 and beta-strand B5. We further identify the RKEN214-217DAVQ mutations of the set of hyperactive mutations within the catalytic domain of SB transposase to be the driving factor that induces change in residue-pair rigidity correlations within SB transposase. Given that a signature RNase H-like structural motif is found in DD[E/D]-transposases and, more broadly, in a large superfamily of polynucleotidyl transferases, our results are relevant to these proteins as well.

Decomposing Dynamical Couplings in Mutated scFv Antibody Fragments into Stabilizing and Destabilizing Effects.

Conformational fluctuations within scFv antibodies are characterized by a novel perturbation-response decomposition of molecular dynamics trajectories. Both perturbation and response profiles are stratified into stabilizing and destabilizing conditions. The linker between the VH and VL domains exhibits the dominant dynamical response by being coupled to nearly the entire protein, responding to both stabilizing and destabilizing perturbations. Perturbations within complementarity-determining regions (CDR) induce rich behavior in dynamic response. Among many effects, stabilizing any CDR loop in the VH domain triggers a destabilizing response in all CDR loops in the VL domain and vice versa. Destabilizing residues within the VL domain are likely to stabilize all CDR loops in the VH domain, and, when these residues are not buried, the CDR loops in the VL domain are also likely to be stabilized. These effects, described by shifts in normal mode characteristics, initiate a propensity for dynamic allostery with possible functional implications in bispecific antibodies.

JED: a Java Essential Dynamics Program for comparative analysis of protein trajectories.

BACKGROUND: Essential Dynamics (ED) is a common application of principal component analysis (PCA) to extract biologically relevant motions from atomic trajectories of proteins. Covariance and correlation based PCA are two common approaches to determine PCA modes (eigenvectors) and their eigenvalues. Protein dynamics can be characterized in terms of Cartesian coordinates or internal distance pairs. In understanding protein dynamics, a comparison of trajectories taken from a set of proteins for similarity assessment provides insight into conserved mechanisms. Comprehensive software is needed to facilitate comparative-analysis with user-friendly features that are rooted in best practices from multivariate statistics. RESULTS: We developed a Java based Essential Dynamics toolkit called JED to compare the ED from multiple protein trajectories. Trajectories from different simulations and different proteins can be pooled for comparative studies. JED implements Cartesian-based coordinates (cPCA) and internal distance pair coordinates (dpPCA) as options to construct covariance (Q) or correlation (R) matrices. Statistical methods are implemented for treating outliers, benchmarking sampling adequacy, characterizing the precision of Q and R, and reporting partial correlations. JED output results as text files that include transformed coordinates for aligned structures, several metrics that quantify protein mobility, PCA modes with their eigenvalues, and displacement vector (DV) projections onto the top principal modes. Pymol scripts together with PDB files allow movies of individual Q- and R-cPCA modes to be visualized, and the essential dynamics occurring within user-selected time scales. Subspaces defined by the top eigenvectors are compared using several statistical metrics to quantify similarity/overlap of high dimensional vector spaces. Free energy landscapes can be generated for both cPCA and dpPCA. CONCLUSIONS: JED offers a convenient toolkit that encourages best practices in applying multivariate statistics methods to perform comparative studies of essential dynamics over multiple proteins. For each protein, Cartesian coordinates or internal distance pairs can be employed over the entire structure or user-selected parts to quantify similarity/differences in mobility and correlations in dynamics to develop insight into protein structure/function relationships.

Statistical Measures to Quantify Similarity between Molecular Dynamics Simulation Trajectories.

Molecular dynamics simulation is commonly employed to explore protein dynamics. Despite the disparate timescales between functional mechanisms and molecular dynamics (MD) trajectories, functional differences are often inferred from differences in conformational ensembles between two proteins in structure-function studies that investigate the effect of mutations. A common measure to quantify differences in dynamics is the root mean square fluctuation (RMSF) about the average position of residues defined by Cα-atoms. Using six MD trajectories describing three native/mutant pairs of beta-lactamase, we make comparisons with additional measures that include Jensen-Shannon, modifications of Kullback-Leibler divergence, and local p-values from 1-sample Kolmogorov-Smirnov tests. These additional measures require knowing a probability density function, which we estimate by using a nonparametric maximum entropy method that quantifies rare events well. The same measures are applied to distance fluctuations between Cα-atom pairs. Results from several implementations for quantitative comparison of a pair of MD trajectories are made based on fluctuations for on-residue and residue-residue local dynamics. We conclude that there is almost always a statistically significant difference between pairs of 100 ns all-atom simulations on moderate-sized proteins as evident from extraordinarily low p-values.

Mutations in Antibody Fragments Modulate Allosteric Response Via Hydrogen-Bond Network Fluctuations.

A mechanical perturbation method that locally restricts conformational entropy along the protein backbone is used to identify putative allosteric sites in a series of antibody fragments. The method is based on a distance constraint model that integrates mechanical and thermodynamic viewpoints of protein structure wherein mechanical clamps that mimic substrate or cosolute binding are introduced. Across a set of six single chain-Fv fragments of the anti-lymphotoxin-ß receptor antibody, statistically significant responses are obtained by averaging over 10 representative structures sampled from a molecular dynamics simulation. As expected, the introduced clamps locally rigidify the protein, but long-ranged increases in both rigidity and flexibility are also frequently observed. Expanding our analysis to every molecular dynamics frame demonstrates that the allosteric responses are modulated by fluctuations within the hydrogen-bond network where the native ensemble is comprised of conformations that both are, and are not, affected by the perturbation in question. Population shifts induced by the mutations alter the allosteric response by adjusting which hydrogen-bond networks are the most probable. These effects are compared using response maps that track changes across each single chain-Fv fragment, thus providing valuable insight into how sensitive allosteric mechanisms are to mutations.

Novel Ricin Subunit Antigens With Enhanced Capacity to Elicit Toxin-Neutralizing Antibody Responses in Mice.

RiVax is a candidate ricin toxin subunit vaccine antigen that has proven to be safe in human phase I clinical trials. In this study, we introduced double and triple cavity-filling point mutations into the RiVax antigen with the expectation that stability-enhancing modifications would have a beneficial effect on overall immunogenicity of the recombinant proteins. We demonstrate that 2 RiVax triple mutant derivatives, RB (V81L/C171L/V204I) and RC (V81I/C171L/V204I), when adsorbed to aluminum salts adjuvant and tested in a mouse prime-boost-boost regimen were 5- to 10-fold more effective than RiVax at eliciting toxin-neutralizing serum IgG antibody titers. Increased toxin neutralizing antibody values and seroconversion rates were evident at different antigen dosages and within 7 days after the first booster. Quantitative stability/flexibility relationships analysis revealed that the RB and RC mutations affect rigidification of regions spanning residues 98-103, which constitutes a known immunodominant neutralizing B-cell epitope. A more detailed understanding of the immunogenic nature of RB and RC may provide insight into the fundamental relationship between local protein stability and antibody reactivity.

Dynamics and thermodynamic properties of CXCL7 chemokine.

Chemokines form a family of signaling proteins mainly responsible for directing the traffic of leukocytes, where their biological activity can be modulated by their oligomerization state. We characterize the dynamics and thermodynamic stability of monomer and homodimer structures of CXCL7, one of the most abundant platelet chemokines, using experimental methods that include circular dichroism (CD) and nuclear magnetic resonance (NMR) spectroscopy, and computational methods that include the anisotropic network model (ANM), molecular dynamics (MD) simulations and the distance constraint model (DCM). A consistent picture emerges for the effects of dimerization and Cys5-Cys31 and Cys7-Cys47 disulfide bonds formation. The presence of disulfide bonds is not critical for maintaining structural stability in the monomer or dimer, but the monomer is destabilized more than the dimer upon removal of disulfide bonds. Disulfide bonds play a key role in shaping the characteristics of native state dynamics. The combined analysis shows that upon dimerization flexibly correlated motions are induced between the 30s and 50s loop within each monomer and across the dimer interface. Interestingly, the greatest gain in flexibility upon dimerization occurs when both disulfide bonds are present, and the homodimer is least stable relative to its two monomers. These results suggest that the highly conserved disulfide bonds in chemokines facilitate a structural mechanism that is tuned to optimally distinguish functional characteristics between monomer and dimer.

Rigidity Emerges during Antibody Evolution in Three Distinct Antibody Systems: Evidence from QSFR Analysis of Fab Fragments.

The effects of somatic mutations that transform polyspecific germline (GL) antibodies to affinity mature (AM) antibodies with monospecificity are compared among three GL-AM Fab pairs. In particular, changes in conformational flexibility are assessed using a Distance Constraint Model (DCM). We have previously established that the DCM can be robustly applied across a series of antibody fragments (VL to Fab), and subsequently, the DCM was combined with molecular dynamics (MD) simulations to similarly characterize five thermostabilizing scFv mutants. The DCM is an ensemble based statistical mechanical approach that accounts for enthalpy/entropy compensation due to network rigidity, which has been quite successful in elucidating conformational flexibility and Quantitative Stability/Flexibility Relationships (QSFR) in proteins. Applied to three disparate antibody systems changes in QSFR quantities indicate that the VH domain is typically rigidified, whereas the VL domain and CDR L2 loop become more flexible during affinity maturation. The increase in CDR H3 loop rigidity is consistent with other studies in the literature. The redistribution of conformational flexibility is largely controlled by nonspecific changes in the H-bond network, although certain Arg to Asp salt bridges create highly localized rigidity increases. Taken together, these results reveal an intricate flexibility/rigidity response that accompanies affinity maturation.

A virtual pebble game to ensemble average graph rigidity.

BACKGROUND: The body-bar Pebble Game (PG) algorithm is commonly used to calculate network rigidity properties in proteins and polymeric materials. To account for fluctuating interactions such as hydrogen bonds, an ensemble of constraint topologies are sampled, and average network properties are obtained by averaging PG characterizations. At a simpler level of sophistication, Maxwell constraint counting (MCC) provides a rigorous lower bound for the number of internal degrees of freedom (DOF) within a body-bar network, and it is commonly employed to test if a molecular structure is globally under-constrained or over-constrained. MCC is a mean field approximation (MFA) that ignores spatial fluctuations of distance constraints by replacing the actual molecular structure by an effective medium that has distance constraints globally distributed with perfect uniform density. RESULTS: The Virtual Pebble Game (VPG) algorithm is a MFA that retains spatial inhomogeneity in the density of constraints on all length scales. Network fluctuations due to distance constraints that may be present or absent based on binary random dynamic variables are suppressed by replacing all possible constraint topology realizations with the probabilities that distance constraints are present. The VPG algorithm is isomorphic to the PG algorithm, where integers for counting "pebbles" placed on vertices or edges in the PG map to real numbers representing the probability to find a pebble. In the VPG, edges are assigned pebble capacities, and pebble movements become a continuous flow of probability within the network. Comparisons between the VPG and average PG results over a test set of proteins and disordered lattices demonstrate the VPG quantitatively estimates the ensemble average PG results well. CONCLUSIONS: The VPG performs about 20% faster than one PG, and it provides a pragmatic alternative to averaging PG rigidity characteristics over an ensemble of constraint topologies. The utility of the VPG falls in between the most accurate but slowest method of ensemble averaging over hundreds to thousands of independent PG runs, and the fastest but least accurate MCC.

Redistribution of flexibility in stabilizing antibody fragment mutants follows Le Châtelier's principle.

Le Châtelier's principle is the cornerstone of our understanding of chemical equilibria. When a system at equilibrium undergoes a change in concentration or thermodynamic state (i.e., temperature, pressure, etc.), La Châtelier's principle states that an equilibrium shift will occur to offset the perturbation and a new equilibrium is established. We demonstrate that the effects of stabilizing mutations on the rigidity ⇔ flexibility equilibrium within the native state ensemble manifest themselves through enthalpy-entropy compensation as the protein structure adjusts to restore the global balance between the two. Specifically, we characterize the effects of mutation to single chain fragments of the anti-lymphotoxin-ß receptor antibody using a computational Distance Constraint Model. Statistically significant changes in the distribution of both rigidity and flexibility within the molecular structure is typically observed, where the local perturbations often lead to distal shifts in flexibility and rigidity profiles. Nevertheless, the net gain or loss in flexibility of individual mutants can be skewed. Despite all mutants being exclusively stabilizing in this dataset, increased flexibility is slightly more common than increased rigidity. Mechanistically the redistribution of flexibility is largely controlled by changes in the H-bond network. For example, a stabilizing mutation can induce an increase in rigidity locally due to the formation of new H-bonds, and simultaneously break H-bonds elsewhere leading to increased flexibility distant from the mutation site via Le Châtelier. Increased flexibility within the VH ß4/ß5 loop is a noteworthy illustration of this long-range effect.

Thermodynamic stability and flexibility characteristics of antibody fragment complexes.

Free energy landscapes, backbone flexibility and residue-residue couplings for being co-rigid or co-flexible are calculated from the minimal Distance Constraint Model (mDCM) on an exploratory dataset consisting of VL, scFv and Fab antibody fragments. Experimental heat capacity curves are reproduced markedly well, and an analysis of quantitative stability/flexibility relationships (QSFR) is applied to a representative VL domain and several complexes in the scFv and Fab forms. Global flexibility in the denatured ensemble typically decreases in the larger complexes due to domain-domain interfaces. Slight decreases in global flexibility also occur in the native state of the larger fragments, but with a concurrent large increase in correlated flexibility. Typically, a VL fragment has more co-rigid residue pairs when isolated compared to the scFv and Fab forms, where correlated flexibility appears upon complex formation. This context dependence on residue- residue couplings in the VL domain across length scales of a complex is consistent with the evolutionary hypothesis of antibody maturation. In comparing two scFv mutants with similar thermodynamic stability, local and long-ranged changes in backbone flexibility are observed. In the case of anti-p24 HIV-1 Fab, a variety of QSFR metrics were found to be atypical, which includes comparatively greater co-flexibility in the VH domain and less co-flexibility in the VL domain. Interestingly, this fragment is the only example of a polyspecific antibody in our dataset. Finally, the mDCM method is extended to cases where thermodynamic data is incomplete, enabling high throughput QSFR studies on large numbers of antibody fragments and their complexes.

Principal component analysis: a method for determining the essential dynamics of proteins.

It has become commonplace to employ principal component analysis to reveal the most important motions in proteins. This method is more commonly known by its acronym, PCA. While most popular molecular dynamics packages inevitably provide PCA tools to analyze protein trajectories, researchers often make inferences of their results without having insight into how to make interpretations, and they are often unaware of limitations and generalizations of such analysis. Here we review best practices for applying standard PCA, describe useful variants, discuss why one may wish to make comparison studies, and describe a set of metrics that make comparisons possible. In practice, one will be forced to make inferences about the essential dynamics of a protein without having the desired amount of samples. Therefore, considerable time is spent on describing how to judge the significance of results, highlighting pitfalls. The topic of PCA is reviewed from the perspective of many practical considerations, and useful recipes are provided.

A case study comparing quantitative stability-flexibility relationships across five metallo-ß-lactamases highlighting differences within NDM-1.

The Distance Constraint Model (DCM) is an ensemble-based biophysical model that integrates thermodynamic and mechanical viewpoints of protein structure. The DCM outputs a large number of structural characterizations that collectively allow for Quantified Stability-Flexibility Relationships (QSFR) to be identified and compared across protein families. Using five metallo-ß-lactamases (MBLs) as a representative set, we demonstrate how QSFR properties are both conserved and varied across protein families. Similar to our characterizations on other protein families, the backbone flexibility of the five MBLs are overall visually conserved, yet there are interesting specific quantitative differences. For example, the plasmid-encoded NDM-1 enzyme, which leads to a fast spreading drug-resistant version of Klebsiella pneumoniae, has several regions of significantly increased rigidity relative to the other four. In addition, the set of intramolecular couplings within NDM-1 are also atypical. While long-range couplings frequently vary significantly across protein families, NDM-1 is distinct because it has limited correlated flexibility, which is isolated within the active site S3/S4 and S11/H6 loops. These loops are flexibly correlated in the other members, suggesting it is important to function, but the others also have significant amounts of correlated flexibility throughout the rest of their structures.

Towards comprehensive analysis of protein family quantitative stability-flexibility relationships using homology models.

The Distance Constraint Model (DCM) is a computational modeling scheme that uniquely integrates thermodynamic and mechanical descriptions of protein structure. As such, quantitative stability-flexibility relationships (QSFR) that describe the interrelationships of thermodynamics and mechanics can be quickly computed. Using comparative QSFR analyses, we have previously investigated these relationships across a small number of protein orthologs, ranging from two to a dozen [1, 2]. However, our ultimate goal is provide a comprehensive analysis of whole protein families, which requires consideration of many more structures. To that end, we have developed homology modeling and assessment protocols so that we can robustly calculate QSFR properties for proteins without experimentally derived structures. The approach, which is presented here, starts from a large ensemble of potential homology models and uses a clustering algorithm to identify the best models, thus paving the way for a comprehensive QSFR analysis across hundreds of proteins in a protein family.