Mechanical Model of Hydrogen Bonds in Protein Molecules

Veiw figure View Figure

Different geometric formation criteria have been developed to predict hydrogen bonds occurrence in the protein molecules ^[19]. In this paper, in order to identify hydrogen bonds we used the suggested geometric criteria by ^[19] shown in Table 1.

Table 1. Geometric criteria [19]

Download as

Veiw figure View Table

2.2. Hydrogen Bond Energy

Hydrogen bond energy depends on the chemistry of the donor and acceptor atoms as well as their orientation. Based on the potential function used in “Dreiding” force field, (^[23]) suggested a function (equation 1) to calculate hydrogen bond energy. In this energy function, constants are: v₀ = 8 kCal/mol^-1and d₀ = 2.8Å.

(1)

Dahiyat energy function (equation 1) considers hydrogen bonds of protein molecules in any solvent. The energy values predicted by this function (due to both positional and angular coordinates) are commonly regarded in by other researchers in reasonable match with experimental data ^[24]. Although, we are using ^[19] geometric criteria and ^[23] energy function for prediction of the stiffness of hydrogen bonds, the proposed method remains valid for any other method of predicting the bond geometry and its energy.

3. Proposed Hydrogen Bond Mechanical Model

Under the application of external forces on the molecule, or internal forces due to shape changes in the functional proteins, hydrogen bonds exhibit flexibility. Hydrogen bonds can tolerate ±0.17 Å change in bond length in average before breaking ^{[25, 26, 27]}. Therefore we modeled hydrogen bonds as flexible mechanical elements (linear springs). Our simulation is valid for small displacements ranging ± 0.12 Å. Considering the energy function of these bonds shown in equation 1 we modeled each individual hydrogen bond with specific geometry as a combination of three linear springs accurately estimating the bond energy change for small changes in the bond geometry. Consequently our model calculates the mechanical properties of the spring such as spring stiffness and spring free length (in three directions) for each individual hydrogen bond. This is the major achievement of this study. Previous mechanical model of hydrogen bonds considered the bond to be flexible solely in bond direction and assumed the same stiffness for all bonds with different geometry ^[28].

3.1. Equivalent Stiffens for Individual Hydrogen Bonds

To calculate stiffness of each bond the energy function of the bond (shown in equation 1) has been used. In this equation the energy is a function of three geometric parameters: donor-acceptor bond length (d), and the two angles (α&β) as shown in Figure 2-A. To simplify the model, angular coordinates (α&β) in energy function are replaced by the positional coordinates l and r, using geometric relations 2 and 3.

(2)

(3)

Figure 2. (A) The Geometry of hydrogen bonds. (B) The equivalent mechanical model of hydrogen bonds

Download as

Veiw figure View Figure

Figure 2 B shows our proposed equivalent mechanical system to model hydrogen bonds. This model consists of three linear springs. The energy of this system should be the same as the energy of the bond itself computed by equation 1. Using Hook’s law the equivalent energy of the spring system is shown in equation 4.

(4)

To fully define the model, we need to calculate the values of the six constants (k_r, r₀, k_l, l₀, k_d, d₀) in the right hand side of the equation 4. In order to obtain these constants, we have set up an optimization problem with the aim of minimizing the fitting error; that is the difference between the bond energy obtained by empirical studies (equation 1) and our proposed model (equation 4).

The optimization problem has the following constraints:

The energy computed at the configuration prior to the small displacements by both methods (equations 1 and 4) should be equal. This condition is implemented in equation 6.

The three spring constants (for the springs in Figure 2) are the slope of graph describing each spring force versus interatomic separation distance ^[29]. These conditions are implemented in equations 7 to 9.

Stiffness and the free lengths of springs should be positive numbers (equations 10 to 15).

The optimization problem (minimizing the difference between the bond energy obtained by empirical studies and our proposed model after small displacement) is therefore summarized as follows:

Minimize M (the objective function)

(5)

Subject to the following constraints:

(6)

(7)

(8)

(9)

(10)

Figure 3. The new energy function calculated verses the LJ function

Download as

Veiw figure View Figure

(11)

(12)

(13)

(14)

(15)

This optimization problem is then solved using Lagrange Multipliers method ^[30]. Figure 3 shows the energy of the hydrogen bond versus the bond length (r) obtained by empirical studies (dashed curve) and our optimized model (solid curve). As shown in Figure 3, the proposed function can properly estimate the empirical energy function for small variations of the bond length.

3.2. Equivalent Stiffness Ellipsoid

Equivalent stiffness ellipsoid Once we have the stiffness of three springs, modeling individual hydrogen bonds, we can calculate bond stiffness at any given direction. To do this, stiffness matrix and consequently stiffness ellipsoid (ellipsoid equation) for each bond is calculated using static equilibrium equations (16, 17 and 18) for the mechanical system shown in Figure 4. In this system link NH is assumed to be fixed and link OC is free to move (from O1C1 to O2C2). The free body diagram of the system is shown in Figure 4-B.

Figure 4. A: Small displacement of atoms result in small change in the length of springs B: Free body diagram of link OC

Download as

Veiw figure View Figure

(16)

(17)

Figure 5. A schematic presentation of stiffness ellipsoid and hydrogen bond

Download as

Veiw figure View Figure

(18)

Rearranging these equations, the stiffness matrix can be determined from equation 19.

(19)

(20)

The maximum and minimum stiffness of the system at point C, and their corresponding directions can be computed from stiffness matrix. The process is as follows. A unit displacement vector at point C is assumed. The corresponding force that is required to cause this displacement is calculated from equation 20. The process of finding the minimum and maximum possible force vectors is then setup as an optimization problem:

It can be shown that these maximum and minimum forces are the eigenvectors of the matrix, with magnitudes of these forces being the square root of the corresponding eigenvalues. The three obtained eigenvalues are the corresponding diameters of the stiffness ellipsoid, and the eigen vectors are ellipsoid’s principal directions.

The stiffness of the bond at any other direction, Ut, can be computed by calculating the intersection of the stiffness ellipsoid and the line passing through the center of the ellipsoid in the Ut direction. Figure 5 shows a schematic presentation of stiffness ellipsoid and its orientation with respect to bond geometry.

4. Results and Model Validations

4.1. For Single Hydrogen Bond

The developed methodology in previous section has been applied on several sample hydrogen bonds. Table 1 provides an overview of the results, detailing the following parameters: (E) the bond energy (kCal/mol), (k_l,k_r&k_d) stiffness (N/m) of each spring, (l₀,r₀&d₀) free length of springs (Å), eigen values of stiffness matrix which are the maximum, mean and minimum stiffness of the bonds (N/m), (θ*) the angle between maximum stiffness direction and hydrogen bond direction (degrees), (k**) the hydrogen bond stiffness in bond (OH) direction (N/m). Table 3 lists stiffness for various bonds obtained from ^[31]. Our reported stiffness for hydrogen bonds in Table 2, ranging from 1.19-5.74 N/m with an average of 3.43N/m and standard deviation of 1.3N/m, are in agreement with the empirical results reported by ^[31] (3-6N/m). It has been shown in Table 2, column 11 which hydrogen bonds’ maximum stiffness direction needs to be rotated 5 to 10 degrees to be aligned with bond direction (OH). On the other hand as shown in columns 8 and 10 of Table 2 the stiffness ellipsoid is so thin because of very small amount of minimum stiffness comparing to the maximum one. These observations lead us to an important conclusion: “hydrogen bonds are much stiffer in the bond direction than in perpendicular direction to the bond”. This achievement combined with the analysis of particular arrangement of hydrogen bonds in α helices (studied in next section) provides mathematical prove for the observation of α helices being stiffer in axial direction than in lateral direction ^[9].

Table 2. Hydrogen bond stiffness for some sample bonds

Download as

Veiw figure View Table

Table 3. Different bonds’ stiffness

Download as

Veiw figure View Table

4.2. Application to Protein Molecules

The secondary structures of protein molecules (α helices and β sheets) are stabilized by hydrogen bonds ^[32]. Based on this fact, one can infer that the stiffness of protein molecules, as well as their strength in different directions, is greatly determined by hydrogen bonds. Therefore to predict the stiffness of protein molecules the individual stiffness of hydrogen bonds has been used.

Figure 6-A illustrates a schematic arrangement of hydrogen bonds in a helices. It is shown that in each turn of a helix there are 3 to 4 hydrogen bonds ^[33]. Each turn has 3.6 residues with 5.4Å pitch along the helical axis ^[32]. The equivalent mechanical model is shown in figure 6-B. In this model, hydrogen bonds in each turn are assumed to be parallel and turns are connected in a serial form. Considering the Hook’s law for springs, the stiffness of the entire a helix can be obtained using equation 21.

Figure 6. A schematic arrangement of hydrogen bonds inα helices

Download as

Veiw figure View Figure

(21)

4.2.1. Calculating stiffness of some sample protein molecules

Now we can use the approach discussed in previous section to calculate the stiffness of several protein molecules. The proteins in this numerical experiment were chosen such that empirical or computed data on their stiffness is readily available in the literature. Secondly, we chose proteins with helical structures because of their specific geometry and arrangement of hydrogen bonds. The helical motif can greatly simplify the evaluation of protein stiffness since hydrogen bonds are placed almost parallel to the helix axes. Table 4 lists three sample proteins with the required specification, their PDB (Protein Data Bank) code, the number of amino acids and the number of detected hydrogen bonds (based on the geometric and energetic criteria reported in ^[21]). The first molecule is a synthetic peptide, cysteine - lysine - cysteine (C3K30C) specifically designed to study hydrogen bonding by ^[6]. Under the experimental conditions ^[6] this synthetic peptide adopts the a helix structure as a result of hydrogen bonding (31 bonds) within the molecule. The peptide has a length of 53Å and a diameter of approximately 15Å. The stiffness of this protein molecule was measured using atomic force microscopy. Each molecule was stretched from the a-helical state into a linear chain. The stiffness was found to vary with molecule displacement. In order to conduct a meaningful comparison between our calculated stiffness of this protein molecule and the empirical measurement reported in ^[6], only small displacements are considered. In particular, we concentrate on deformation of 0.12Å or smaller for the subject protein. With such deformation protein molecule will keep its helical shape and none of the hydrogen bonds breaks. The reported stiffness of this protein molecule varies between 0.3-0.4N/m. Results of our computational analysis applied on this protein molecule are summarized as follows. 31 hydrogen bonds are detected within this protein molecule.

The stiffness of the studied molecule (0.38N/m) is in agreement with the expected range of the stiffness of this molecule ((0.3-0.4N/m)) from empirical data.

The elastic properties of characterized myosin II S2 sub domain using molecular dynamics and normal mode analysis has been studied ^[9]. Tis protein molecule with 1NKN PDB code, is 87 residue long a-helix. We predicted 146 hydrogen bonds within this molecule structure and 0.083N/m stiffness. This results are in good agreement with the stiffness predicted by ^[9] for this protein molecule (0.06-0.08N/m).

Table 4. Sample proteins

Download as

Veiw figure View Table

To study elasticity and strength of secondary structures of protein molecules ^[1] compared the force-extension behavior of α-helices, β-sheets and tropocollagen domains in protein molecules. Their results are provided by using atomistic modeling of nano-mechanical response of the protein molecules at ultra-slow deformation rates. The stiffness predicted for a molecule with helical structure (1GK6) in this study is 0.571N/m. The predicted stiffness for this molecule using our proposed method is 0.52N/m.

Table 5. Stiffness for different protein molecules (N/m)

Download as

Veiw figure View Table