Structural Rigidity and Protein Thermostability in Variants of Lipase A from Bacillus subtilis

Understanding the origin of thermostability is of fundamental importance in protein biochemistry. Opposing views on increased or decreased structural rigidity of the folded state have been put forward in this context. They have been related to differences in the temporal resolution of experiments and computations that probe atomic mobility. Here, we find a significant (p = 0.004) and fair (R 2 = 0.46) correlation between the structural rigidity of a well-characterized set of 16 mutants of lipase A from Bacillus subtilis (BsLipA) and their thermodynamic thermostability. We apply the rigidity theory-based Constraint Network Analysis (CNA) approach, analyzing directly and in a time-independent manner the statics of the BsLipA mutants. We carefully validate the CNA results on macroscopic and microscopic experimental observables and probe for their sensitivity with respect to input structures. Furthermore, we introduce a robust, local stability measure for predicting thermodynamic thermostability. Our results complement work that showed for pairs of homologous proteins that raising the structural stability is the most common way to obtain a higher thermostability. Furthermore, they demonstrate that related series of mutants with only a small number of mutations can be successfully analyzed by CNA, which suggests that CNA can be applied prospectively in rational protein design aimed at higher thermodynamic thermostability.


Body-and-bar networks
In body-and-bar networks, atoms are considered bodies with six degrees of freedom, and each bar between two bodies removes one degree of freedom. Depending on the strength of an interaction between two atoms, a constraint can be modeled as any number of bars between one and six with six bars completely freezing the motion between two atoms [1,2].
As done previously [3][4][5][6], covalent bonds were modeled with five (single bonds) and six (peptide and double bonds) bars, whereas hydrophobic tethers and hydrogen bonds (including salt bridges; together referred to as hydrogen bonds here) were modeled with two and five bars, respectively. A modified version of the potential by Mayo and coworkers [7] as described in ref. [8] was used to calculate hydrogen bond energies EHB; hydrogen bonds with energies lower than a certain cutoff Ecut were included in the network (see "Thermal unfolding simulation" section in the main text for details). Hydrophobic constraints were considered between pairs of carbon and/or sulfur atoms according to a Gaussian probability function depending on the distances between the atoms (dij), the sum of their van der Waals (dvdw) radii (C: 1.7 Å; S: 1.8 Å), and the full width at half maximum Dcut (eq. S1; see ref. [3] for details). (S1)

Local and global rigidity indices
From a thermal unfolding trajectory, one can calculate both residue-level (local) and overall (global) rigidity characteristics of a protein [9]. Local indices can be used to investigate specific questions regarding the stability and activity of a protein. In the present study, the stability map rcij introduced by us in ref. [5] was used to characterize the local rigidity of BsLipA and to understand the influence of mutations. A stability map is derived by identifying "rigid contacts" between two residues i and j that are represented by their Cα atoms. A rigid contact exists if the two residues belong to the same rigid cluster. During a thermal unfolding simulation, stability maps are then constructed in that, for each residue pair, Ecut (or, equivalently, a temperature derived from the relationship T = f(Ecut) described in refs. [4,5]) is identified at which a rigid contact between these residues is lost. In that respect, the stability map is a two-dimensional itemization of the local rigidity index as detailed in ref. [9]. When filtered such that only rigid contacts between residues that are 5 Å apart from each other (measured as the distance between the closest atom pair from the two residues) are considered, a neighbor stability map results. This map helps focusing on short-range residue contacts that can be directly modulated by mutagenesis with the aim to stabilize them for improving the overall stability of a protein.
In addition, the (local) percolation index pi introduced by us in ref. [9] was used to characterize thermal unfolding pathways of the BsLipA structures. The percolation index monitors the percolation behavior (i.e., the loss of rigidity when diluting the constraint network) of a biomolecule on a microscopic level and so allows identifying the hierarchical break-down of the giant percolating cluster during a thermal unfolding simulation. The giant percolating cluster is the largest rigid cluster present at the highest Ecut value (i.e., at the lowest temperature at the beginning of a thermal unfolding simulation) with all constraints in place. More technically, pi monitors the Ecut at which a bond segregates from the giant percolating cluster during a thermal unfolding simulation. For a Cα atom-based representation, the lower of the pi values of the two backbone bonds is considered.
Global indices help identifying phase transition temperatures Tp at which a network switches from being largely rigid to largely flexible. Previously, we showed that Tp identified by a modified cluster configuration entropy Htype2 [4,9] can be used for predicting the thermodynamic thermostability of and identifying structural weak spots in a protein [4][5][6].
The cluster configuration entropy has originally been introduced by Andraud et al. [10] as a morphological descriptor for heterogeneous materials and is adapted from Shannon's information theory. Htype2 monitors the degree of disorder in the realization of a given network state: As long as a network is dominated by a very large rigid cluster, Htype2 tends to be low because there are only a few configurations of a system with a large rigid cluster possible; Htype2 increases when larger rigid clusters break down in smaller clusters. The Htype2 versus T curve obtained from a thermal unfolding simulation was fitted with a double sigmoid [11] as done previously [6], and the temperature Tp was identified as the inflection point of the sigmoid with the larger difference in the asymptote values. This way, in most cases, a late transition involving the final decay of the giant percolating cluster is identified as Tp [5].        (Table ) I) of all BsLipA variants but the two outliers, wild type and 6B, (blue) and II) of only the two outliers (red). A normal kernel function with an optimal smoothing parameter [12] at each data point was used for calculating the PDFs.    similarly, for hydrophobic tethers, only hydrophobic atoms (C, S) whose van der Waals spheres are within 0.35 Å were considered.