Kinematic Flexibility Analysis: Hydrogen Bonding Patterns Impart a Spatial Hierarchy of Protein Motion
KKinematic Flexibility Analysis:Hydrogen Bonding Patterns Impart a SpatialHierarchy of Protein Motion
Dominik Budday, ⇤ , † Sigrid Leyendecker, † and Henry van den Bedem ⇤ , ‡ † Chair of Applied Dynamics, University of Erlangen-Nuremberg, Erlangen, Germany ‡ Biosciences Division, SLAC National Accelerator Laboratory, Stanford University,California, Menlo Park, USA
E-mail: [email protected]; [email protected]
Abstract
Elastic network models (ENM) and constraint-based, topological rigidity analysis are two dis-tinct, coarse-grained approaches to study con-formational flexibility of macromolecules. Inthe two decades since their introduction, bothhave contributed significantly to insights intoprotein molecular mechanisms and function.However, despite a shared purpose of theseapproaches, the topological nature of rigidityanalysis, and thereby the absence of motionmodes, has impeded a direct comparison. Here,we present an alternative, kinematic approachto rigidity analysis, which circumvents thesedrawbacks. We introduce a novel protein hy-drogen bond network spectral decomposition,which provides an orthonormal basis for col-lective motions modulated by non-covalent in-teractions, analogous to the eigenspectrum ofnormal modes, and decomposes proteins intorigid clusters identical to those from topologi-cal rigidity. Our kinematic flexibility analysisbridges topological rigidity theory and ENM,and enables a detailed analysis of motion modesobtained from both approaches. Our analysisreveals that collectivity of protein motions, re-ported by the Shannon entropy, is significantlylower for rigidity theory versus normal mode ap-proaches. Strikingly, kinematic flexibility anal-ysis suggests that the hydrogen bonding net-work encodes a protein-fold specific, spatial hi- erarchy of motions, which goes nearly unde-tected in ENM. This hierarchy reveals distinctmotion regimes that rationalize protein sti↵nesschanges observed from experiment and molec-ular dynamics simulations. A formal expres-sion for changes in free energy derived from thespectral decomposition indicates that motionsacross nearly 40% of modes obey enthalpy-entropy compensation. Taken together, ouranalysis suggests that hydrogen bond networkshave evolved to modulate protein structure and dynamics.
Introduction
Coarse-grained modeling techniques can pro-vide significant insight into the dynamic behav-ior and biological function of macromolecules.Elastic Network Models (ENM ) and rigiditytheory based flexibility analysis are two wellknown approaches that have been applied ex-tensively to study molecular motion. Althoughthey are based on physically distinct concepts,both aim to distinguish more flexible regionsfrom compact ones in the molecule, and to pre-dict large-scale, functional motions.ENM (Figure 1, right) approximate the poten-tial energy function V via pairwise harmonic in-teractions V = P ij C ij ( | r ij | | r ij | ) , where C ij denotes the sti↵ness or an exponentially fad-ing weight of the restraint between atoms i, j a r X i v : . [ q - b i o . B M ] F e b ebble DoFbar constraints dihedral/cartesian DoFspring restraintsconstraint graph kinematic tree spring-mass networkTopological rigidity Kinematic flexibility Elastic Network Modelspebble game constraint Jacobian J eigenmodesrigid clusters and floppy modes full explicit motion basisran( J T )null( J ) eigenfrequenciesconstraint perturbationszero-frequency / zero-perturbation Figure 1: Coarse-grained protein modeling via topological rigidity analysis (left) and elastic networkmodels (ENM, right). Topological rigidity analysis decomposes a protein into rigid clusters ofatoms from explicit constraints, quantifying the number of ’floppy modes’ without an explicitmotion basis. By contrast, ENM obtains an explicit motion basis corresponding to eigenmodesof one-dimensional restraints. Our kinematic flexibility approach (center) combines features fromtopological rigidity and ENM, providing motion modes from a spectral decomposition of explicitconstraints corresponding to, for example, hydrogen bonds.located at r i , r j , respectively, with r ij = r j r i and rest length | r ij | . The resulting spring-mass network can be formulated in terms ofdihedral or cartesian degrees of freedom(DoF) and analyzed with classical Hamiltonianmechanics. Diagonalization of the Hessian ma-trix leads to a spectral distribution of the dy-namics in terms of eigenmodes with correspond-ing eigenfrequencies. The simplified energyfunction of ENM places the native structure ata global minimum, circumventing initial mini-mization necessary in traditional normal modeanalysis (NMA). The motions correspond-ing to low-frequency modes are robust, oftenagree with those of more detailed models, and can be functionally relevant. Dihedral-based approaches often correlate better with ex-perimentally observed conformational changesthan cartesian-based methods.
Other vari-ants include Gaussian or anisotropic net-work models.In rigidity theory, macromolecules are modeledas constraint graphs (Figure 1, left), with edges between interacting atoms (vertices) represent-ing covalent and non-covalent bonds. Di↵er-ent types of constraint graphs, such as the bar-joint, the body-bar, or the equivalentbody-bar-hinge graph all share the conceptof assigning a number of pebbles to atomic ver-tices as DoF, and a number of bars to each in-teraction as constraints. Typically, rotatablesingle-covalent bonds and hydrogen bonds re-tain a dihedral degree of freedom, while pep-tide and double covalent bonds are modeledrigid. More constraints, such as additional hy-drogen bonds, increasingly rigidify the graph.The remaining flexibility, i.e., the coordinatedmotion of internal DoF that do not perturbconstraints, also known as floppy modes , can be determined via constraint countingthrough the pebble game algorithm. Con-straint counting on these graphs has been ap-plied to examine conformational flexibility inmacromolecules, probe the e↵ects of lig-and binding, or estimate thermo-stability and entropic measures. Rigidity informa-3ion has also helped to shape successful per-turbation strategies in complex motion plan-ning algorithms.
However, the topologi-cal nature of this approach denies explicit ac-cess to the kinematics of the underlying 3-dimensional molecular structure. Importantly,the exact motion vectors determined by theconstraints (covalent and non-covalent bonds)remain unknown and can only be approxi-mated by randomized perturbations and iter-ative loop-closure algorithms such as ROCK or FRODA. While previous e↵orts tried tocombine pebble game rigidity analysis withENM, these impediments have so far pre-vented a detailed comparison between the twomethods.Here, we overcome these limitations by us-ing a kinematic approach to characterize flex-ibility and rigidity. Our approach explicitlyprovides basis vectors for motions correspond-ing to floppy modes, allowing us to analyzeand compare motions directly from ENM andrigidity analysis. We previously establishedthat topological rigidity and kinematic analy-sis yield an identical decomposition of the pro-tein into rigid clusters for constraint compliantmotions in non-singular configurations. Topo-logical rigidity analysis fails for protein con-formations corresponding to singular kinematicconfigurations, where kinematic analysis givescorrect decompositions. Furthermore, here weexploit that our kinematic analysis can providea basis for motions that progressively perturbconstraints, i.e., motion modes that are inac-cessible to topological rigidity analysis. In theremainder, we distinguish between topologicalrigidity as counting on constraint graphs viathe pebble game, and kinematic flexibility asour new method.Our study makes three main contributions.First, we establish that the hydrogen bond-ing pattern of a protein structure, togetherwith its rotatable bonds, encodes a hierarchyof proteins motions. These orthogonal motion’modes’ are ranked by perturbations of the hy-drogen bond network energy. Second, we showthat this motion spectrum is conserved acrossa large dataset of di↵erent protein structures,but simultaneously reveals fold-specific di↵er- ences. This fold-specific, spatial hierarchy sug-gests that hydrogen bond networks and foldsare designed to modulate protein structure and dynamics. Third, we derive a formal, quali-tative expression for a mode-specific free en-ergy, which displays a near constant regimewhere increasing energetic cost is compensatedby elevated entropy. Thus, our kinematic anal-ysis provides a conceptual model system forenthalpy-entropy-compensation, with a read-ily accessible interpretation of the controver-sial phenomenon. Collectively, our results sig-nify a deep connection between the structure ofthe hydrogen bond network, protein conforma-tional dynamics, and functional motions. Ourkinematic flexibility analysis is implemented inour KGS software, and available from https://github.com/ExcitedStates . Kinematic flexibility analy-sis of proteins
We model macromolecules as a kinematic link-age, with groups of atoms as rigid body ver-tices and rotatable covalent bonds as directededges with a single degree of freedom (Fig-ure 1, center). The linear, branched topol-ogy of monomeric molecules can be representedwith a single kinematic spanning tree, rootedat an arbitrarily selected vertex (pink). Co-valent double bonds, peptide bonds, and di-hedrals amenable to planarity are modeled asrigid. Proline and aromatic rings are also mod-eled rigid, identical to standard topological ap-proaches.
Groups of atoms without inter-nal DoF are joined into rigid bodies. Dihe-dral angles of remaining single-covalent bondsare the DoF of the molecule, connecting twoneighboring vertices. In contrast to topolog-ical rigidity, we distinguish between covalentDoF and non-covalent constraints such as hy-drogen bonds. Without constraints (red, Fig-ure 1 center), our model corresponds to a se-rial, open kinematic chain. Non-covalent inter-actions form closed kinematic cycles that coor-dinate dihedral motion, reducing the number ofindependent DoF. Multi-chain proteins can betreated within the same framework by connect-4ng chains through inter-chain constraints orto a virtual super-root.
We encode hydro-gen bonds as pentavalent holonomic constraintsthat allow only a rotation about the bond axisbut prevent all other relative motion. In otherwords, the distance between the hydrogen andacceptor atom as well as the donor-hydrogen-acceptor and hydrogen-acceptor-base angles areconstrained. This corresponds to a constraintwith five bars in pebble game algorithms, butexplicitly specifies its kinematics. Other non-covalent interactions such as hydrophobic forcesare not considered in this work. For a moleculewith overall n dihedral DoF, we can distinguishbetween free DoF q f T f that do not appear inclosed kinematic cycles, and constrained, cycleDoF q T d , where f + d = n . The m hydrogenbonds define a constraint variety Q on the cycleDoF Q = { q T d | ( q ) = R m } . (1)Di↵erentiating the constraints with respect totime leads to instantaneous (velocity) con-straint equationsd d t = J ˙ q ( = , if ˙ q constraint observing = , if ˙ q constraint violating (2)which characterizes two disjoint, orthogonalsubspaces corresponding to velocities that ob-serve constraints and those that violate con-straints. Here, J R m ⇥ d is the Jacobian ma-trix of the constraints. Constraint observingvelocities span a subspace null ( J ( q )) of dimen-sion d r , with r the rank of J . Note that r p = min( d, m ) due to the rectangularshape of J . The singular value decomposition(SVD, ) JV = U⌃ , (3)where U = [ u , . . . , u m ] R m ⇥ m and V [ v , . . . , v d ] = R d ⇥ d , provides orthonormalbases for the range and nullspace of J . The rect-angular matrix ⌃ = diag( , . . . , p ) R m ⇥ d contains the singular values on the diagonal,where . . . r > r +1 = . . . = p = 0,and corresponding u i and v i are the i th left and right singular vector, respectively. Letran( J T ) = span { v , . . . , v r } =: R , null( J ) = span { v r +1 , . . . , v d } =: N . (4)Then, any ˙ q = [ R , N ] ⇥ ⌫ TR ,⌫ TN ⇤ T can be ex-pressed with ⌫ i the proportion of velocity point-ing along vector v i , respectively.Since u i and v i are orthonormal vectors, i en-codes non-orthogonality between J and V ; thesingular value encodes the norm of constraintperturbation when moving along v i . There-fore, the two subspaces provide physically dis-tinct insights. For motions along the nullspace,the perturbation is zero, encoded by the van-ishing singular value. Thus, null( J ) spans allconstraint-observing motions. As shown pre-viously, these modes are identical to floppymodes from topological rigidity approaches, i.e., they are the only remaining internal DoF ifno constraint perturbation is admitted. Dihe-drals not part of these floppy modes are rigidi-fied and merge adjacent rigid bodies, leading toa rigid cluster decomposition of the molecule,which is a central result from topological rigid-ity analysis. The nullspace is highly sensitive tothe set of non-covalent constraints and carefultuning of energy cuto↵s, e.g., for the inclusionof hydrogen bonds in protein and RNA, is cru-cial to prevent over-rigidification. It is alsonoteworthy that the nullspace dimension scalesnearly linearly with the size of the protein,i.e., the number of floppy modes increases withincreasing dimensions of the Jacobian (Figure2A).The range ran( J T ) provides a spectrum of mo-tions ranked by increasing constraint pertur-bation. Consistent with the notion in penaltymethods, the singular value i v i , i = 1 , . . . , r . We term those motions kinematic flexibility modes . They are not ac-cessible in topological rigidity analysis, sincethey do not respect constraints, but can poten-tially inform on hydrogen bonds that need toadapt their configuration during functional re-arrangements or exchange with solvent. Thus,5
500 1000 1500 2000 cycle DoF d f ( d ) floppy modesfit, r=0.79 max fit, r=0.8constraintsfit, r=0.99 A B l og l og [ kc a l / m o l ] || geometry|||| energy||nullspace limit C Figure 2: The hydrogen bond network spectrum predicts kinematic and energetic hydrogen bondperturbations, shown for an individual example (A, PDB ID 1p5f) with d = 620 and a large,diverse dataset with 183 high-resolution single-chain proteins (B). Modes are plotted in reverseindex order, leading to a monotonous increase in . Geometry and energy changes are computedfrom small perturbations (step size = 1 e
5) along kinematic flexibility modes. The rankingof energy perturbations by kinematic flexibility modes and the near log-linear motion regime areremarkably conserved among all proteins. Motion modes (horizontal axis) in panel B are limitedto R and indices are normalized by r to remove size-related di↵erences in the dataset. (C) Thenumber of hydrogen bond constraints, remaining floppy modes, as well as the magnitude of thelargest singular value vary fairly linearly (correlation coe cient r ) with the number of DoF d inkinematic cycles, which serves as indicator for protein size.the full range of protein conformational dynam-ics encoded by the hydrogen bonding patternis accessible through orthonormal bases span-ning both matrix subspaces. We therefore referto the SVD (3) of the constraint Jacobian ob-tained from the hydrogen bonding pattern asthe hydrogen bond network spectral decomposi-tion . Results
Kinematic flexibility modes hier-archically rank collective hydrogenbond energy perturbation
First, we analyzed how the hydrogen bond net-work spectral decomposition imposes a rangeof collective motions on the protein. To estab-lish a predictive relationship between singularvalues and the cumulative variation in proteinhydrogen bond geometry and energy, we mon-itored geometry and energy changes while tak-ing a small step along the direction of eachindividual singular vector q ,i = v i . (5) The vectors v i are all unit length, and therefore || q ,i || = , i = 1 , . . . , d . Changes in hydro-gen bond geometry follow from || ( q ,i ) || (1),while energies || E HB ( q ,i ) || are evaluated us-ing the Mayo energy potential for each indi-vidual hydrogen bond E HB = D ( ✓ R R ◆ ✓ R R ◆ ) f ( ✓, , ) , (6)with well-depth D , equilibrium distance R ,hydrogen bond donor-acceptor distance R , andangular terms f that depend on the hybridiza-tion state of donor and acceptor (details see SI(11)). We illustrate our findings using the crys-tal structure of 189-residue human DJ-1 protein(PDB ID 1p5f). Figure 2B shows monotoni-cally increasing singular values (blue) which re-versely rank-order the kinematic motion modesalong the horizontal axis. The magnitudes ofthe corresponding cumulative geometric pertur-bations (magenta) and hydrogen bond energies(red) are shown for = 1 e i are re-lated by a scale factor c V = = 1 e
5, i.e., || J q ,i || = c V i (combining (3) and (5)). Fit-ting the singular values curve (blue) to the mag-nitude of geometric perturbations of hydrogenbonds || ( q ,i ) || = c G i (1), we found a scalefactor of c G = 1 e ( . ± . ⇡ , indi-cating that linearization has negligible e↵ectsat small step sizes. For the highly non-linearenergy function ((6) and SI (11)), we found c E = 1 e ( . ± . ⇡ . . Hence, fora su ciently small step size , the hydrogenbond network spectral decomposition predictsthe cumulative hydrogen bond energy pertur-bation along mode i by || E HB ( q ,i ) || = c E i for this example structure. For floppy modes in-side the nullspace, the di↵erence between singu-lar value prediction and energetic cost is larger.This can be explained by sp2 sp2 hybridizedhydrogen bonds that observe small energeticchanges due to a torsional term in the hydrogenbond energy function (SI), which is not presentin the geometric constraint formulation. Thetorsional term becomes more relevant if otherterms such as distance and angles remain small,which is the case inside the nullspace.This hierarchy of protein motions is remark- ably conserved in the protein universe. Wecompiled a diverse benchmark dataset of 183high-resolution, non-redundant crystal struc-tures ranging in size from 30 to 555 amino acidsin length from the PDB. We first observed anexpected, strong linear relationship between thenumber of cycle DoF d and the number of con-straints 5 m in the structures (Figure 2A). Table1 evaluates statistical properties of the dataset.Similar to our single example crystal structureof human DJ-1 protein, singular values spannedmany orders of magnitude. We also observeda large range in the number and distributionof vanishing singular values across the crys-tal structures, a median of 34 nullspace floppymodes, with a maximum of 148 and one (engi-neered) protein with zero floppy modes (PDBID 5eca). While most structures have slightlymore constraints than cycle DoF, i.e., 5 m > d ,it appears that a fairly constant fraction of con-straints is linearly independent, leading to alinear increase in the number of floppy modesover protein size (Figure 2A). This was alsoobserved previously. We therefore analyzedkinematic flexibility modes corresponding tonon-vanishing singular values, i.e., outside ofthe nullspace, for all 183 crystal structures (Fig-ure 2C). For comparison we normalized modes,and we grouped modes into 50 bins per struc-ture. Repeating our analysis above, we ob-tained a surprisingly universal law from fittingthe mean curves( E q HB ) = 1 . . . || E HB ( q ,i ) || . . . = c E (7)with c E ⇡ . independent of mode number,and error of the same order as before. Thissuggests that protein hydrogen bonding pat-terns impart a distribution of orthogonal, coor-dinated motions on the DoF. Kinematic flexi-bility modes identify preferred directions of de-formation for the protein in (hydrogen bond)internal energy landscapes. Note that ( E q HB ) reports on the cumulative energy change; mo-tions along kinematic flexibility modes i do notnecessarily increase all hydrogen bond energies.7able 1: Statistics of the high-resolution dataset with 183 single-chain proteins. m d d r h-bond energy median 163 116 516 34 -422.94 kcal/mol 485.67min 30 23 107 0 -1506.81 kcal/mol 82.94max 555 445 2032 148 -59.88 kcal/mol 1966.29 Instead, for each direction, fluctuations in hy-drogen bond energies can provide compensatorymechanisms, i.e., many could marginally reducein energy to allow a handful to significantly in-crease. Initial constrained minimization of thedataset ensures starting structures near a localminimum of hydrogen bond energy.
Kinematic flexibility modes of ↵ -helices and -sheets To visualize motions across the spectrum andtheir e↵ect on hydrogen bonds, we analyzed an ↵ -helix and a -sheet in detail. The top panelsof Figure 3 show the matrix product JV , accu-mulated over the five constraints per hydrogenbonds, thereby graphically displaying equation(3). The predicted perturbation of individualhydrogen bonds (rows) for motions along eachkinematic flexibility mode (columns) is color-coded, increasing from dark to light. For the ↵ -helix (A), we observed a striking, alternatingpattern of perturbation, suggesting motions aremodulated by concerted hydrogen bond energyfluctuations. We observed a similar, but weakerpattern for the -sheet (B).We then projected the conformational change q,j for degree of freedom j corresponding to se-lected kinematic flexibility modes onto the he-lix and sheet (Figure 3, bottom panels). Theleft-most, purple-colored set of modes in thematrix correspond to local fluctuations (e.g.,A11, B5, indices label column numbers), af-fecting only DoF close to individual hydrogenbonds shown as dots. As the mode number in-creases, DoF are increasingly engaged. For ex-ample, in super-helical twisting (e.g., A38) andstraight bending (e.g., A46) compensating mo-tions of DoF, indicated by alternating blue-redpatterns, mediate moderate constraint pertur-bation. By contrast, high constraint perturba-tions found in compression/tension (A54/A55) are a result of changes in DoF that reinforceeach other. Interestingly, unraveling (e.g., A19)demonstrates significantly less overall perturba-tion than compression/tension (A54/A55), in-dicating an energetically favorable mode for he-lix dissociation. Similarly, the structure of -sheets permits twisting and bending motions(e.g., B11), while shear (e.g., B33) and espe-cially widening/narrowing (e.g., B38) requiressubstantial distortion of hydrogen bonds. Ourfindings agree well with molecular dynamicssimulations on the mechanical response of typ-ical protein building blocks. They revealedhigh strength of -sheet proteins in responseto shear loading, with an elastic modulus ofabout 240pN / ˚A. Also, the simulation resultspredicted an initial linear elastic regime dur-ing tensile loading of ↵ -helical protein domains.Beyond the elastic regime the helix uncoils, re-leasing one turn at a time. It is worth notingthat there are no nullspace floppy modes alongthe backbone in either structure, i.e., topolog-ical rigidity analysis predicts only side-chainflexibility, but full backbone rigidity.Besides the graphical interpretation, these re-sults hint at two important characteristics ofthe motion spectrum. First, the perturbationpatterns confirm distinct motion regimes, sep-arated by perturbation levels and collectivity,i.e., how much of the structure is engaged ineach mode. Second, the hydrogen bonding net-works encode fold-specific motions. We exam-ine these attributes quantitatively in the follow-ing. Spectral distribution of the hydro-gen bond network
To obtain a physical understanding of the dis-tribution of modes across protein structure, weperformed a full spectral decomposition on ourbenchmark dataset. The size distribution of8igure 3: Hierarchical constraint perturbation in an ↵ -helix (A panels) and an anti-parallel -sheet(B panels). The two top panels plot the matrix JV , ranking singular vectors (columns) acrossindividual hydrogen bonds (rows) by cumulative constraint perturbation, increasing from purple toyellow. Bottom panels depict motions along selected kinematic flexibility modes (column indicesgiven). The motion amplitude is exaggerated (step size = 1) for visualization purposes. Colorcoding indicates increasing change in DoF, from blue to red. The original conformation is shownin green.protein structures is reflected in the SVD ofassociated constraint Jacobian matrices J . Itis well-known mathematically that the max-imum singular values monotonically increaseby adding a column (a degree of freedom),while the smallest non-zero singular values de-crease. Physically, a larger protein with moreDoF can be interpreted by a larger lever armthat allows a larger maximum constraint per-turbation, or simultaneously provides more mo-tion capabilities to avoid perturbation. Addinga row (constraint) will increase the minimumsingular value, i.e., there are more constraintsto be perturbed, and the protein is rigidi-fied. Overall, the spectral range becomes size- dependent, as can be seen in the widened dis-tributions near the spectral limits in Figure 2C.Hence, to compare mode densities across di↵er-ently sized proteins, we normalized singular val-ues of each structure by their maximum toobtain [0 ,
1] for all modes across all struc-tures and grouped them by into ten bins perdecimal power.Singular values for the 183 single-chain pro-teins span several orders of magnitude and showwell-conserved, sharp peaks near integer expo-nents (Figure 4A). The two most common kine-matic flexibility modes occur at ⇡ e ⇡ e
3. The peak at lowest values (Figure4A, far left) represents floppy modes, which are9he only modes considered in topological rigid-ity analysis. Their density is least conservedacross the spectrum and shows a linear increaseover protein size (Figure 2A). Thus, while nor-malizing by helps spectral comparison in theregime >
0, the variability in the number offloppy modes where = 0 remains.The spectral distribution relates directly to thesti↵ness of the protein. Formally, this followsfrom defining a cumulative perturbation p c = Z p ( ) d (8)for perturbation-specific probability densities p ( ) following the spectral distribution, as-suming individual modes are enabled at equalprobability. Thus, proteins enriched in high-perturbation modes require overall more en-ergy to access their motion modes, renderingthem sti↵er. Discrete jumps between spec-tral peaks suggest that modes are distributedacross di↵erent sti↵ness regimes. Interestingly,atomic force microscopy (AFM) on single an-tibody proteins also found two distinct elasticregions with a ⇠ Again, theseexperimental findings agree well with our per-turbation analysis.The geometry of hydrogen bonds in our analy-sis critically depends on accurately placed hy-drogen atoms in structure preparation. To ex-clude the possibility of bias at that stage werepeated our analysis using a smaller set of 34structures from the PDB with experimentallydetermined hydrogen atom positions from neu-tron di↵raction experiments (SI). Its spectralsignature, i.e., peak locations and heights, wasvirtually identical, validating the results fromour initial dataset of 183 high-resolution single-chain structures for subsequent analysis. Fur-thermore, spectral analysis of singular valuesfrom a set of random matrices with the samedimensions as the original protein dataset pro-duced a completely di↵erent distribution, whilea set of random matrices with the same spar-sity pattern led to similar distributions. Thisdemonstrates that the hydrogen bonding net-work, together with the kinematic structure of the protein, stores this fold-specific dynamic in-formation.
Motions from the hydrogen bond-ing pattern are spatially dis-tributed
Besides distinct perturbation levels, we ob-served various levels of mode engagement, i.e.,how much of the molecule is involved in a spe-cific motion mode. To measure this collectivity s of mode v i , we compute the exponential of theShannon entropy of its squared components s i = 1 d exp d X j =1 i,j log( i,j ) ! , (9)where i,j = v i,j / P dj =1 v i,j . This normalizationof v i,j is trivial for singular vectors, since theyare unit length by definition, but non-trivial foreigenmodes from ENM or NMA. The secondnormalization by the number of modes d (equalsthe length of v ) reduces size di↵erences acrossprotein structures, which allows s i [0 ,
1] tobe interpreted as the fraction of significant con-tributors in motion mode i .Figure 4B plots collectivity computed overmodes with similar singular values, groupedinto ten bins per decimal power. Nullspacefloppy modes with zero perturbation at thelower end of the spectrum show a relativelylow collectivity compared to medium or high-perturbation modes. From panels A and B to-gether we observe that the most collective mo-tion modes are also most abundant. Collectiv-ity, and thus entropy, follows a near-exponentialregime between 1 e e
3, showing a sim-ilar trend as energetic perturbations in Figure2B. Finally, modes with strongest constraintperturbation are less collective than medium-perturbation modes.
Protein fold classes have a uniqueh-bond network spectral signature
To analyze protein-fold specific di↵erences en-coded by the hydrogen bond pattern, we ex-amined the hydrogen bond spectrum and col-10
12 -10 -8 -6 -4 -2 0, power of 1000.050.10.15 a v e r age p r obab ili t y den s i t y +/ -12 -10 -8 -6 -4 -2 0, power of 1000.050.10.150.2 a v e r age c o ll e c t i v i t y +/ A BC D
Figure 4: Kinematic flexibility analysis with hydrogen bond constraints reveals fold-specific motionsin proteins. (A) Spectrum of normalized singular values for 183 single-chain proteins, with striking,well-conserved peaks. The nullspace dimension (peak at lowest ) varies with protein size. (B)Collectivity of motion from Shannon entropy. Floppy modes with zero perturbation (lowest )show relatively small collectivity compared to modes at higher perturbations. Across a medium-perturbation range, collectivity increases with increasing singular values. (C) Fold-specific, averagespectrum for classes of ↵ -only, -only, ↵ + and ↵/ proteins (stair representation only for bettervisibility). Peak locations are well conserved across folds, while peak heights are shifted to lower from ↵ -only to -only. Floppy modes (lowest ) vary rather with protein size than fold. (D)Average collectivity as in (B) for fold-specific datasets. Again, floppy modes with zero perturbationshow relatively small collectivity; they are indistinguishable across folds.lectivity of kinematic flexibility modes for fourseparate datasets of ↵ -only, -only, ↵ + and ↵/ proteins, ranging from 655 to 1051 struc-tures (SI). Figure 4 shows the mean curves forspectrum (C) and collectivity (D), correspond-ing to the black mean curves of the originalmixed-fold dataset in (A) and (B), respectively.Percentiles follow similar patterns as for theprevious data-set and are omitted for visibil-ity. While the location of spectral peaks isconserved across folds, the peak heights showclear di↵erences. The class of ↵ -only folds hasmore modes at higher singular values, yield-ing them sti↵er to perturbations (8). By con-trast, -only proteins are enriched in smallerspectral modes, rendering them more flexible,while mixed folds ↵ + and ↵/ show interme-diate spectra. The density of zero-perturbationfloppy modes with vanishing singular valuesis highly size-dependent and less informative regarding fold content. Nonetheless, our re-sults predict more floppy modes in ↵ -helicalfolds than -folds, which can most likely be at-tributed to less inter-helical connectivity com-pared to inter- -strand connectivity.Collectivity, computed as before, is generallyhigher for -folds than ↵ -folds, and interme-diate for mixed folds. Interestingly, medium-perturbation modes in -dominated folds aremore collective than the highest-perturbationmodes, indicating localized hinges, e.g. in -turns, that are able to significantly unfold theprotein (cf. Figure 3B40). Note that floppymodes with zero-perturbation, the only acces-sible modes from topological rigidity analysis,show relatively small and very similar collectiv-ity across folds, rendering them more localizedcompared to medium-perturbation modes andless informative regarding fold content.To further evaluate predictive capabilities of11ur method, we analyzed the spectral distribu-tion of a set of four hyperstable, designed pep-tides, each with an NMR bundle of 20 dis-tinct structures (PDB codes 2nd2, 2nd3, 5jhi,5ji4; details in SI). Their spectrum shows aclear shift towards modes with higher relativeperturbation compared to the high-resolutiondataset (Figure 5) and a ⇠ p c (8), identifying their designed constraint pat-tern as a key contributor to increased stability.Collectivity of modes shows trends as the otherdatasets.Figure 5: Spectrum of normalized singular val-ues for four designed hyperstable, constrainedpeptides, each consisting of 20 NMR structures.The spectrum is shifted to higher perturbationmodes relative to the high-resolution dataset. Free energy of modes
The hydrogen bond energy perturbations pre-dicted by the magnitude of singular values canbe formally combined with the conformationalentropy contributions encoded by collectivityinto a dimensionless expression for free-energychanges F related to each mode i F i = i c T s i , (10)with a dimensionless temperature factor c T .For normalized singular values and normalizedcollectivities, the range of F is between -1and +1 when c T = 1. Clearly, (10) is a for-mal abstraction, since only hydrogen bond en-ergy contributes to internal energy. Surpris- ingly, though, our free-energy changes demon-strate how enthalpic and entropic contributionscompensate each other in the overall spectrumof conformational motion. Figure 6 depicts F over normalized motion modes computed fromall 183 structures in the high-resolution dataset.The red curve averages over individual modesgrouped into 100 bins. F roughly levels forabout 40% of modes (Figure 6 inset) with mostfavorable entropy (collectivity) and medium en-thalpic cost. Modes to the right of the intervalhave unfavorably high enthalpic cost; they cor-respond to unfolding (high-perturbation modesin Figure 3). Modes to the left of the inter-val are more localized, with smaller enthalpiccost, but simultaneously less entropic benefit.Nullspace floppy modes are excluded in thegraph; their enthalpic cost is zero, as encodedby vanishing singular values. Thus, associ-ated free energy is identical to their collectivityas depicted in Figure 4B, which is on average F = .
05 and thus turns out less favorablethan the free energy plateau.Figure 6: Dimensionless free energy of modesdemonstrates entropy-enthalpy compensationencoded by the hydrogen bonding pattern. A40% indicates a near constant free-energy levelfor highly collective modes at acceptable en-thalpic cost (inset). Motions to the right corre-spond to unfolding, while motions to the left aremore localized, with smaller entropic benefit.12 omparison with iMOD
Finally, we compared the spectrum and collec-tivity of motions obtained from kinematic flexi-bility analysis with normal modes from iMOD. iMOD is a versatile ENM-based tool to studynormal modes of macromolecules in internalcoordinates, i.e., dihedral angles. AlthoughiMOD’s mass matrix is based on a full atomrepresentation, flexibility is limited to main-chain and angles, while our model maintainsfull side-chain flexibility. We carried out iMODsimulations with the command line settings n x to enable the degree of freedom,and otherwise default parameters. Comparisonof reduced ENM models to full NMA previ-ously revealed that agreement of low-frequencymodes is conserved, but that higher-frequencymodes can di↵er significantly. Figure 7 shows the iMOD spectrum of eigenfre-quencies (A) and collectivity of eigenmodes (B)for our 183 single-chain protein dataset. Simi-larly to our spectrum, eigenfrequency distribu-tions are broadly conserved across structures,confirming results from previous NMA andENM analysis.
Modes with low to mediumeigenfrequency are most abundant. Similar tokinematic flexibility, the most abundant modesare also most collective. Higher-frequency mo-tions with !> cannot be obtainedwith this coarse-grained model and require fullNMA. Zero-perturbation floppy modes in topo-logical rigidity and our kinematic flexibilityanalysis correspond to zero-frequency modes inENM, i.e., modes with vanishing energeticcost. However, zero-frequency modes in ENMare often considered an artifact, as the networkbreaks down into multiple independent ones.These modes are often avoided by increasingthe number of weak interactions, for exampleby increasing cut-o↵ distances. Floppy mo-tions from topological rigidity therefore cannotbe directly accessed with ENM.Zero-perturbation floppy modes from topologi-cal rigidity analysis are less collective, i.e., cor-respond to more local motions compared to low-frequency normal modes or medium kinematicflexibility modes (compare Figure 4B to Figure 7B). This signifies that topological rigidity the-ory based methods tend to overestimate molec-ular rigidity, and study conformational flexibil-ity based on more local motions than ENM orNMA.When we compared fold-specific spectra (Fig-ure 7C) and collectivity (Figure 7D), we ob-served several important di↵erences betweenthe methods. For example, the location andheight of the main spectral peak in iMODis indistinguishable for ↵ -only and -only,and shifted slightly to higher frequencies formixed folds, while kinematic flexibility analy-sis shows stronger di↵erences in peak heightsacross fold types. Interestingly, ↵ -folds showslightly increased density at very low frequen-cies in iMOD, predicting lower sti↵ness com-pared to other folds, in contrast to other estab-lished methods. Our kinematic flexibility anal-ysis, detailed NMA, and analysis of force-displacement curves from MD simulations all predict that ↵ -folds are sti↵er. However,iMOD predicts lower collectivity in ↵ -folds,similar to our approach.A more detailed comparison between the meth-ods, e.g., analyzing direct mode overlap, re-mains di cult, since iMOD is limited to main-chain and dihedral motion, while our kine-matic model maintains full side-chain flexibil-ity, thus demanding lower-dimensional projec-tions. Other ENM or NMA tools are mostlyformulated in cartesian coordinates or are un-available for other users, making a comparisoneven more di cult. Nevertheless, our resultsreveal similarities and di↵erences important forusers of ENM, topological rigidity, and our newkinematic analysis. Discussion
Our new, kinematic approach to rigidity anal-ysis treats hydrogen bonds as a geometric con-straint network. Compared to topological rigid-ity, it extends analyses to constraint perturb-ing motions, providing a full spectral decom-position of motion modes ranked by their cu-mulative hydrogen bond energy perturbations.While analyses of the structural dynamics of13
200 400 600 800 1000[cm -1 ]00.020.040.060.08 a v e r age p r obab ili t y den s i t y +/ -1 ]00.050.10.150.2 a v e r age c o ll e c t i v i t y +/ A BC D
Figure 7: Eigenspectrum and collectivity for normal modes computed with iMOD. (A) Eigenfre-quency spectra of the data-set with 183 single-chain proteins. (B) Collectivity of eigenmodes. (C)Eigenfrequency spectra of four fold-specific data-sets. (D) Fold-specific collectivity of eigenmodes.proteins often implicitly assume that weak hy-drogen bonds disrupt first, our approach doesnot require such assumptions. Instead, the net-work imparts a hierarchy of motions and sti↵-ness regimes on the protein, which modulate theconformational response. The hydrogen bondspectral decomposition is 1) highly conservedin the protein universe, and 2) reveals key fold-specific di↵erences.Kinematic flexibility analysis indicated thatzero-perturbation floppy modes, the motionsobtained from topological rigidity analysis, aremore localized than low-frequency modes fromENM, suggesting an overly rigidified repre-sentation in topological rigidity or alterna-tively, over-connectivity in ENM. Therefore,strict rigidity theory based methods rely onexceedingly local motions to estimate confor-mational flexibility, ligand binding, en-tropy, or thermo-stability compared toENM or NMA. Nonetheless, they show convinc-ing agreement with experimental data. More-over, kinematic floppy modes that observe hy-drogen bond constraints have guided confor-mational transitions and revealed coor-dinated loop motions, often more successfulthan normal mode based methods. Analyses of motions beyond the floppy modes, i.e., in the kinematic flexibility regime, revealedseveral unexpected insights. First, we foundthat hydrogen bond networks determine struc-ture and modulate structural dynamics. Thiscould have important implications for de novo protein design and folding or hydrogels, where recent attention is focused on design-ing hydrogen bonding patterns to create sta-ble interfaces and mediate specificity. Forexample, our procedure revealed a clear shifttoward sti↵er modes for designed, hyperstablepeptides. While hydrogen-bond guided de-signs are often successful structurally, i.e., thecrystal and predicted structure are near identi-cal, it remains a challenge to design dynamic,functional proteins. Our procedure could pre-dict motion modes of hundreds of designed pro-teins and hydrogen bond networks in minutes.Second, the conserved hierarchy and collectiv-ity of protein motions revealed distinct mo-tion regimes, which are often observed in ex-periments. For example, AFM on single anti-body proteins uncovered two elastic motionregimes, separated by a near 4-fold increase insti↵ness, and a regime of plastic deformation.Our spectral decomposition structurally ratio-nalizes the distribution of these motions un-der di↵erent strains. Low strain triggers low-14erturbation modes that are low in collectiv-ity and mostly locally perturb the structure.Higher strain engage highly collective modes to-ward sti↵er motion regimes, explaining elevatedsti↵ness in the second elastic regime measuredwith AFM. Interestingly, our EEC model sug-gests that a fraction of the energetic cost couldbe entropically balanced, rendering the defor-mation elastic and reversible. Any strainbeyond the elastic regime leads to plastic defor-mation, indicated by a sharp increase in hydro-gen bond energy perturbation. Motion modesin this regime likely completely unfold the pro-tein.Third, the fold-specific di↵erences we observein our perturbation analysis suggest distinctfunctional roles of secondary structure, con-firming previous simulations.
We foundmore modes at higher energy perturbationin ↵ -helices than -sheets, yielding helicessti↵er. This is consistent with experimentaldata from low-frequency Raman spectroscopy and eigenfrequencies from detailed NMA. Detailed NMA captures additional di↵erencesbetween folds in high-frequency regimes suchas the amide I, II, or III bands. These higher-frequency di↵erences can also be detected ex-perimentally with infrared spectroscopy andhave been used to determine secondary struc-ture content.
Simplified elastic networkanalysis using iMOD failed to identify this shifttoward sti↵er modes, highlighting a limitationin the simplified force-fields of most ENM. Ouranalysis also predicted that -sheets are highlyresistant to shear motions, which is confirmedby MD simulations. Simulations on the me-chanical response of silk crystalline units undershear loading showed high rupture forces due toe cient force distribution in the -sheet struc-tures. Buehler and Keten further found aninitial linear elastic regime during tensile load-ing of ↵ -helical proteins, after which the helixunravels turn by turn. Our analysis also pre-dicted unraveling as the favored mode of helixdissociation.Our method can provide quick insight intohow rigid clusters, collective motions, and sti↵-nesses shift as constraints are added or re-moved. While we found convincing agreement with experiment and simulation, important lim-itations remain. Our model considers onlyintra-molecular hydrogen bonds, ignoring othernon-covalent interactions. Bond lengths andangles are fixed and only torsional DoF con-tribute to molecular motion. While hydro-gen bonds are important determinants of pro-tein structure and dynamics, hydrophobic ef-fects, electrostatics, solvent, or protein-proteininteractions also modulate structural dynamics.The e↵ects of these interactions could be ex-plored with our method. For example, neutrondi↵raction reveals the position and orientationof hydrogen atoms in waters. Adding water-mediated hydrogen bonds can help understandhow solvation propagates collective motions inprotein cavities or binding pockets.Our new, kinematic flexibility analysis is a ver-satile method, bridging topological rigidity andENM. By way of a novel spectral decomposi-tion of protein hydrogen bonding patterns, itprovides explicit access to collective motionsand free energy of modes, signifying that hy-drogen bonds store intrinsic, fold-specific func-tional motions. These quantitative models andinsights can help improve de novo protein de-sign and folding, help to understand mechano-biology probed by AFM or single-molecule flu-orescence resonance energy transfer (smFRET)at the molecular level, or help interpret exper-imental data such as hydrogen deuterium ex-change.
Acknowledgement
The authors gratefully ac-knowledge financial support to D.B. from theDeutsche Telekom Stiftung. H.v.d.B. is sup-ported by NIH GM123159.
Datasets
The initial high-resolution dataset of 183 pro-teins was obtained through an advanced entity-based PDB search limited to crystallographi-cally collected single-chain protein structuresbetween 30 and 1000 amino acids in length,with a single oligomeric state and model, atmost 50% sequence identity and no free or mod-ified residues ligands. To ensure high validity,we required a resolution of at least 1.5 ˚ A as15ell as R free and R work values below 0.2. Thedataset was then curated with Schrodinger’sMaestro suite filling small gaps, performing aconstrained minimization with an RMSD of 0.3(0.4 in one single case where no minimum wasfound otherwise), and optimizing hydrogen po-sitions for hydrogen bonding. One structure,2 qvk , had to be excluded due to a large miss-ing domain. Water molecules were included inthe Maestro preparation step, but ignored insubsequent KGS analysis.The neutron scattering control dataset was sim-ilarly obtained from the PDB, specifying exper-imental method to neutron di↵raction with ex-perimental data present and limiting the searchto single-chain proteins at 90% sequence iden-tity, leaving 34 structures in the dataset. Struc-tures were also prepared using Maestro, thistime without resampling hydrogens to preserveexperimentally determined positions and orien-tations.The fold-specific datasets were obtained fromthe PDBs SCOP classification, limited tosingle-chain, single oligomeric state proteinswith at most 90% sequence identity, leading to655 ↵ -only, 877 -only, 885 ↵ + , and 1051 ↵/ examples, respectively. Due to the size of thedatasets, we refrained from full Maestro-basedpreparation and added hydrogen atoms withReduce instead.The four constrained peptides were downloadedfrom the PDB via accession codes 2nd2, 2nd3,5jhi, and 5ji4, respectively. Each of their twentyindividual NMR structures was processed likethe fold-specific datasets, by adding hydrogenswith Reduce and successively identifying hydro-gen bond constraints. Hydrogen bond identification andenergy computation
Hydrogen bonds are traditionally identifiedbased on stringent geometric criteria onhydrogen-acceptor distance, acceptor-donordistance, and various angles between them. Wedeveloped our own Python-based implemen-tation using criteria developed in HBPlus, limited to strong hydrogen bonds between anelectro-negative donor and acceptor (oxygen, nitrogen, or sulfur).There are several ways to estimate hydrogenbond energy. Di↵erent approaches includeelectrostatic modeling, quantum-mechanicalmodeling based on electron density, and ex-tensions with the density’s Laplacian. Allthese, however, require detailed experimentaldata and are thus quite complex.Traditionally, the Mayo energy function is usedin topological rigidity analysis (see full text(6)). Its full definition is given by E HB = D ( ✓ R R ◆ ✓ R R ◆ ) f ( ✓, , ) , (11)with R = 2 . A and D = 8kcal / mol are thecorresponding equilibrium distance and well-depth, and R the donor-acceptor distance. An-gular terms depend on hybridization, with thefollowing specifications:sp donor - sp acceptor: f = cos ✓ cos ( . , sp donor - sp acceptor: f = cos ✓ cos , sp donor - sp acceptor: f = cos ✓, sp donor - sp acceptor: f = cos ✓ cos (max [ , ]) , (12)where ✓ describes the donor-hydrogen-acceptorangle, the hydrogen-acceptor-base angle, and the angle between normals of the planes de-fined by donor and acceptor bonds. PDB basedevaluations show that a heuristic-based en-ergy function might produce more realistic en-ergies, especially for hydrogen bonds outsidethe classical geometric configuration favored bythe Mayo potential. Here, however, we areinterested in energetic changes, where we ex-pect only marginal di↵erences between both ap-proaches, such that the commonly used energyfunction in rigidity analysis su ces.For all datasets, hydrogen bonds were identifiedwith an in-house python implementation, usingkinematic criteria as defined in and an energycut-o↵ of / mol using the above ”Mayo”potential. All hydrogen bonds stronger thanthat were included as constraints in our kine-matic model. Other non-covalent interactions16ere not considered. iMOD iMOD is developed and maintained by the Cha-con lab and performs normal mode analysiswith an elastic network model in internal, tor-sional coordinates. For our comparative anal-ysis, we ran iMOD on each structure in ourdatasets, setting option n x to enable the dihedral degree of freedom and chose defaultoptions otherwise. This representation consid-ers all heavy atoms individually, but only main-chain and flexibility in the resulting modes. References (1) Tirion, M. M. Large amplitude elastic mo-tions in proteins from a single-parameter,atomic analysis.
Phys Rev Lett , ,1905.(2) Laman, G. On graphs and rigidity of planeskeletal structures. J Eng Math , ,331–340.(3) Jacobs, D. J.; Thorpe, M. F. Generic rigid-ity percolation: the pebble game. PhysRev Lett , , 4051.(4) Lop´ez-Blanco, J. R.; Garz´on, J. I.;Chac´on, P. iMod: multipurpose normalmode analysis in internal coordinates. Bioinformatics , , 2843–2850.(5) Noguti, T.; G¯o, N. A method of rapid cal-culation of a second derivative matrix ofconformational energy for large molecules. J Phys Soc Jpn , , 3685–3690.(6) Tirion, M. M.; ben Avraham, D.Atomic torsional modal analysis for high-resolution proteins. Physical Review E , , 032712.(7) Durand, P.; Trinquier, G.; Sane-jouand, Y.-H. A new approach fordetermining low-frequency normal modes in macromolecules. Biopolymers , , 759–771.(8) Tama, F.; Gadea, F. X.; Marques, O.;Sanejouand, Y.-H. Building-block ap-proach for determining low-frequency nor-mal modes of macromolecules. Proteins:Struct , Funct , Bioinf , , 1–7.(9) Suhre, K.; Sanejouand, Y.-H. ElNemo: anormal mode web server for protein move-ment analysis and the generation of tem-plates for molecular replacement. NucleicAcids Res , , W610–W614.(10) Levitt, M.; Sander, C.; Stern, P. S. Proteinnormal-mode dynamics: trypsin inhibitor,crambin, ribonuclease and lysozyme. JMol Biol , , 423–447.(11) Na, H.; Song, G. Conventional NMA as abetter standard for evaluating elastic net-work models. Proteins: Structure, Func-tion, and Bioinformatics , , 259–267.(12) Na, H.; Jernigan, R. L.; Song, G. Bridgingbetween nma and elastic network models:preserving all-atom accuracy in coarse-grained models. PLoS computational biol-ogy , , e1004542.(13) Nicolay, S.; Sanejouand, Y.-H. Functionalmodes of proteins are among the most ro-bust. Phys Rev Lett , , 078104.(14) Mendez, R.; Bastolla, U. Torsional net-work model: normal modes in torsion an-gle space better correlate with conforma-tion changes in proteins. Physical reviewletters , , 228103.(15) Bray, J. K.; Weiss, D. R.; Levitt, M. Opti-mized torsion-angle normal modes repro-duce conformational changes more accu-rately than cartesian modes. Biophys J , , 2966–2969.(16) Bahar, I.; Atilgan, A. R.; Erman, B.Direct evaluation of thermal fluctuationsin proteins using a single-parameter har-monic potential. Fold Des , , 173–181.1717) Bahar, I.; Wallqvist, A.; Covell, D.;Jernigan, R. Correlation between native-state hydrogen exchange and cooperativeresidue fluctuations from a simple model. Biochemistry (Mosc ) , , 1067–1075.(18) Atilgan, A.; Durell, S.; Jernigan, R.;Demirel, M.; Keskin, O.; Bahar, I.Anisotropy of fluctuation dynamics of pro-teins with an elastic network model. Bio-phys J , , 505–515.(19) Jacobs, D. J.; Hendrickson, B. An algo-rithm for two-dimensional rigidity perco-lation: the pebble game. J Comput Phys , , 346–365.(20) Thorpe, M.; Lei, M.; Rader, A.; Ja-cobs, D. J.; Kuhn, L. A. Protein flexibil-ity and dynamics using constraint theory. J Mol Graphics Modell , , 60–69.(21) Jacobs, D. J. Generic rigidity in three-dimensional bond-bending networks. JPhys A: Math Gen , , 6653.(22) Whiteley, W. Counting out to the flexibil-ity of molecules. Phys Biol , , S116.(23) Katoh, N.; Tanigawa, S.-i. A proof ofthe molecular conjecture. Proceedings ofthe 25th annual symposium on Computa-tional geometry. 2009; pp 296–305.(24) Fox, N.; Jagodzinski, F.; Li, Y.; Streinu, I.KINARI-Web: A Server for Protein Rigid-ity Analysis. Nucl. Acids Res. , ,W177–W183.(25) Cai, Y.; Thorpe, M. Floppy modes in net-work glasses. Physical Review B , ,10535.(26) Chubynsky, M.; Thorpe, M. Algorithmsfor three-dimensional rigidity analysis anda first-order percolation transition. Physi-cal Review E , , 041135.(27) Jacobs, D. J.; Rader, A. J.; Kuhn, L. A.;Thorpe, M. F. Protein flexibility predic-tions using graph theory. Proteins: Struct, Funct , Bioinf , , 150–165. (28) Pfleger, C.; Radestock, S.; Schmidt, E.;Gohlke, H. Global and local indices forcharacterizing biomolecular flexibility andrigidity. J Comput Chem , , 220–233.(29) Rader, A.; Hespenheide, B. M.;Kuhn, L. A.; Thorpe, M. F. Proteinunfolding: rigidity lost. Proceedings of theNational Academy of Sciences , ,3540–3545.(30) Fulle, S.; Gohlke, H. Analyzing the flex-ibility of RNA structures by constraintcounting. Biophys J , , 4202–4219.(31) Raschka, S.; Bemister-Bu ngton, J.;Kuhn, L. A. Detecting the native ligandorientation by interfacial rigidity: SiteIn-terlock. Proteins: Structure, Function,and Bioinformatics , , 1888–1901.(32) Rader, A. Thermostability in rubredoxinand its relationship to mechanical rigidity. Physical biology , , 016002.(33) Radestock, S.; Gohlke, H. Protein rigid-ity and thermophilic adaptation. Proteins:Structure, Function, and Bioinformatics , , 1089–1108.(34) Rathi, P. C.; Fulton, A.; Jaeger, K.-E.;Gohlke, H. Application of Rigidity The-ory to the Thermostabilization of LipaseA from Bacillus subtilis. PLoS computa-tional biology , , e1004754.(35) Vorov, O. K.; Livesay, D. R.; Jacobs, D. J.Nonadditivity in conformational entropyupon molecular rigidification reveals a uni-versal mechanism a↵ecting folding cooper-ativity. Biophys J , , 1129–38.(36) Gohlke, H.; Ben-Shalom, I. Y.; Kopitz, H.;Pfei↵er-Marek, S.; Baringhaus, K.-H.Rigidity theory-based approximation ofvibrational entropy changes upon bindingto biomolecules. J Chem Theory Comput , , 1495–1502.1837) Ben-Shalom, I. Y.; Pfei↵er-Marek, S.;Baringhaus, K.-H.; Gohlke, H. E -cient Approximation of Ligand Rotationaland Translational Entropy Changes uponBinding for Use in MM-PBSA Calcu-lations. Journal of chemical informationand modeling , , 170–189.(38) Abella, J. R.; Moll, M.; Kavraki, L. E.Maintaining and Enhancing Diversityof Sampled Protein Conformations inRobotics-Inspired Methods. Journal ofComputational Biology , , 3–20.(39) Thomas, S.; Tang, X.; Tapia, L.; Am-ato, N. M. Simulating protein motionswith rigidity analysis. J Comput Biol , , 839–855.(40) Luo, D.; Haspel, N. Multi-resolutionrigidity-based sampling of protein confor-mational paths. Proceedings of the In-ternational Conference on Bioinformatics,Computational Biology and BiomedicalInformatics. 2013; p 786.(41) Lei, M.; Zavodszky, M. I.; Kuhn, L. A.;Thorpe, M. Sampling protein conforma-tions and pathways. Journal of computa-tional chemistry , , 1133–1148.(42) Wells, S.; Menor, S.; Hespenheide, B.;Thorpe, M. Constrained geometric simu-lation of di↵usive motion in proteins. PhysBiol , , S127.(43) Ahmed, A.; Gohlke, H. Multiscale mod-eling of macromolecular conformationalchanges combining concepts from rigid-ity and elastic network theory. Proteins:Struct , Funct , Bioinf , , 1038–1051.(44) Jimenez-Roldan, J. E.; Freedman, R.;Roemer, R. A.; Wells, S. A. Rapid simula-tion of protein motion: merging flexibility,rigidity and normal mode analyses. PhysBiol , , 016008.(45) Budday, D.; Leyendecker, S.; van denBedem, H. Geometric analysis character-izes molecular rigidity in generic and non- generic protein configurations. J MechPhys Solids , , 36–47.(46) Pachov, D. V.; van den Bedem, H.Nullspace sampling with holonomic con-straints reveals molecular mechanisms ofprotein G ↵ s. PLoS Comput Biol , ,e1004361.(47) Pachov, D. V.; Fonseca, R.; Arnol, D.;Bernauer, J.; van den Bedem, H. Coupledmotions in ↵ s conformational en-sembles. J Chem Theory Comput , , 946–956.(48) Golub, G. H.; Van Loan, C. F. Matrixcomputations ; JHU Press, 2012; Vol. 3.(49) Leyendecker, S.; Betsch, P.; Steinmann, P.Energy-conserving integration of con-strained Hamiltonian systems–a compar-ison of approaches.
Computational Me-chanics , , 174–185.(50) Leyendecker, S.; Betsch, P.; Steinmann, P.Objective energy–momentum conservingintegration for the constrained dynamicsof geometrically exact beams. ComputerMethods in Applied Mechanics and Engi-neering , , 2313–2333.(51) Dahiyat, B. I.; Benjamin Gordon, D.;Mayo, S. L. Automated design of the sur-face positions of protein helices. ProteinSci , , 1333–1337.(52) Buehler, M. J.; Keten, S. Elasticity,strength and resilience: A comparativestudy on mechanical signatures of ↵ -helix, -sheet and tropocollagen domains. NanoResearch , , 63–71.(53) Perrino, A. P.; Garcia, R. How soft is asingle protein? The stress–strain curveof antibody pentamers with 5 pN and 50pm resolutions. Nanoscale , , 9151–9158.(54) Br¨uschweiler, R. Collective protein dy-namics and nuclear spin relaxation. TheJournal of chemical physics , ,3396–3403.1955) Bhardwaj, G. et al. Accurate de novo de-sign of hyperstable constrained peptides. Nature , , 329–335.(56) Na, H.; Song, G.; ben Avraham, D. Uni-versality of vibrational spectra of globu-lar proteins. Physical Biology , ,016008.(57) Keten, S.; Buehler, M. J. Geometric con-finement governs the rupture strength ofH-bond assemblies at a critical lengthscale. Nano Letters , , 743–748.(58) Budday, D.; Fonseca, R.; Leyendecker, S.;Bedem, H. v. d. Frustration-guided mo-tion planning reveals conformational tran-sitions in proteins. Proteins: Struct ,Funct , Bioinf , , 1795–1807.(59) H´eliou, A.; Budday, D.; Fonseca, R.;van den Bedem, H. Fast, clash-free RNAconformational morphing using molecularjunctions. Bioinformatics , btx127.(60) Fonseca, R.; Budday, D.; Bedem, H. v. d.Collision-Free Poisson Motion Planning inUltra High-Dimensional Molecular Con-formation Spaces.
J Comput Chem ,JCC25138, in print.(61) Lee, S.; Wang, C.; Liu, H.; Xiong, J.;Jiji, R.; Hong, X.; Yan, X.; Chen, Z.; Ham-mel, M.; Wang, Y.; Dai, S.; Wang, J.;Jiang, C.; Zhang, G. Hydrogen bonds area primary driving force for de novo proteinfolding.
Acta Crystallographica Section D:Structural Biology , , 955–969.(62) Zhang, J.; Wang, N.; Liu, W.; Zhao, X.;Lu, W. Intermolecular hydrogen bondingstrategy to fabricate mechanically stronghydrogels with high elasticity and fatigueresistance. Soft Matter , , 6331–6337.(63) Stranges, P. B.; Kuhlman, B. A compari-son of successful and failed protein inter-face designs highlights the challenges ofdesigning buried hydrogen bonds. ProteinScience , , 74–82. (64) Huang, P.-S.; Boyken, S. E.; Baker, D.The coming of age of de novo protein de-sign. Nature , , 320–327.(65) Boyken, S. E.; Chen, Z.; Groves, B.; Lan-gan, R. A.; Oberdorfer, G.; Ford, A.;Gilmore, J. M.; Xu, C.; DiMaio, F.;Pereira, J. H.; Sankaran, B.; Seelig, G.;Zwart, P. H.; Baker, D. De novo designof protein homo-oligomers with modularhydrogen-bond network–mediated speci-ficity. Science , , 680–687.(66) van den Bedem, H.; Fraser, J. S. Integra-tive, dynamic structural biology at atomicresolution – it’s about time. Nat Methods , , 307–318.(67) Buehler, M. J.; Wong, S. Y. Entropicelasticity controls nanomechanics of singletropocollagen molecules. Biophysical jour-nal , , 37–43.(68) Painter, P.; Mosher, L.; Rhoads, C. Low-frequency modes in the Raman spectraof proteins. Biopolymers , , 1469–1472.(69) Cai, S.; Singh, B. R. A distinct util-ity of the amide III infrared band forsecondary structure estimation of aque-ous protein solutions using partial leastsquares methods. Biochemistry (Mosc ) , , 2541–2549.(70) Yang, H.; Yang, S.; Kong, J.; Dong, A.;Yu, S. Obtaining information about pro-tein secondary structures in aqueous so-lution using Fourier transform IR spec-troscopy. Nature protocols , , 382–396.(71) Xiao, S.; Stacklies, W.; Cetinkaya, M.;Markert, B.; Gr¨ater, F. Mechanical re-sponse of silk crystalline units from force-distribution analysis. Biophysical journal , , 3997–4005.(72) McDonald, I. K.; Thornton, J. M. Satis-fying hydrogen bonding potential in pro-teins. J Mol Biol , , 777–793.2073) Momany, F.; McGuire, R. F.; Burgess, A.;Scheraga, H. A. Energy parameters inpolypeptides. VII. Geometric parameters,partial atomic charges, nonbonded inter-actions, hydrogen bond interactions, andintrinsic torsional potentials for the nat-urally occurring amino acids. The Jour-nal of Physical Chemistry , , 2361–2381.(74) Oliveira, B.; Pereira, F.; de Ara´ujo, R.;Ramos, M. The hydrogen bond strength:New proposals to evaluate the intermolec-ular interaction using DFT calculationsand the AIM theory. Chem Phys Lett , , 181–184.(75) Lifson, S.; Hagler, A.; Dauber, P. Con-sistent force field studies of intermolec-ular forces in hydrogen-bonded crystals.1. Carboxylic acids, amides, and the C:O. cntdot.. cntdot.. cntdot. H-hydrogenbonds. J Am Chem Soc , , 5111–5121.(76) Biegler-k¨onig, F. W.; Bader, R. F.;Tang, T.-H. Calculation of the averageproperties of atoms in molecules. II. JComput Chem , , 317–328.(77) Grabowski, S. J. Ab Initio calculationson conventional and unconventional hy-drogen bonds study of the hydrogen bondstrength. J Phys Chem A , ,10739–10746.(78) Kortemme, T.; Morozov, A. V.; Baker, D.An orientation-dependent hydrogen bond-ing potential improves prediction of speci-ficity and structure for proteins andprotein–protein complexes. J Mol Biol ,326