Gaussian network model

Figure 2 shows a schematic view of elastic network studied in GNM. Metal beads represent the nodes in this Gaussian network (residues of a protein) and springs represent the connections between the nodes (covalent and non-covalent interactions between residues). For nodes i and j, equilibrium position vectors, R⁰_i and R⁰_j, equilibrium distance vector, R⁰_ij, instantaneous fluctuation vectors, ΔR_i and ΔR_j, and instantaneous distance vector, R_ij, are shown in Figure 2. Instantaneous position vectors of these nodes are defined by R_i and R_j. The difference between equilibrium position vector and instantaneous position vector of residue i gives the instantaneous fluctuation vector, ΔR_i = R_i - R⁰_i. Hence, the instantaneous fluctuation vector between nodes i and j is expressed as ΔR_ij = ΔR_j - ΔR_i = R_ij - R⁰_ij.

The potential energy of the network in terms of ΔR_i is

${\displaystyle V_{GNM}={\frac {\gamma }{2}}\left[\sum _{i,j}^{N}(\Delta R_{j}-\Delta R_{i})^{2}\right]={\frac {\gamma }{2}}\left[\sum _{i,j}^{N}\Delta R_{i}\Gamma _{ij}\Delta R_{j}\right]}$

where γ is a force constant uniform for all springs and Γ_ij is the ijth element of the Kirchhoff (or connectivity) matrix of inter-residue contacts, Γ, defined by

${\displaystyle \Gamma _{ij}=\left\{{\begin{matrix}-1,&{\mbox{if }}i\neq j&{\mbox{and }}R_{ij}\leq r_{c}\\0,&{\mbox{if }}i\neq j&{\mbox{and }}R_{ij}>r_{c}\\-\sum _{j,j\neq i}^{N}\Gamma _{ij},&{\mbox{if }}i=j\end{matrix}}\right.}$

r_c is a cutoff distance for spatial interactions and taken to be 7 Å for amino acid pairs (represented by their α-carbons).

Expressing the X, Y and Z components of the fluctuation vectors ΔR_i as ΔX^T = [ΔX₁ ΔX₂ ..... ΔX_N], ΔY^T = [ΔY₁ ΔY₂ ..... ΔY_N], and ΔZ^T = [ΔZ₁ ΔZ₂ ..... ΔZ_N], above equation simplifies to

${\displaystyle V_{GNM}={\frac {\gamma }{2}}[\Delta X^{T}\Gamma \Delta X+\Delta Y^{T}\Gamma \Delta Y+\Delta Z^{T}\Gamma \Delta Z]}$

In the GNM, the probability distribution of all fluctuations, P(ΔR) is isotropic

${\displaystyle P(\Delta R)=P(\Delta X,\Delta Y,\Delta Z)=p(\Delta X)p(\Delta Y)p(\Delta Z)}$

and Gaussian

${\displaystyle p(\Delta X)\propto \exp \left\{-{\frac {\gamma }{2k_{B}T}}\Delta X^{T}\Gamma \Delta X\right\}=\exp \left\{-{\frac {1}{2}}\left(\Delta X^{T}\left({\frac {k_{B}T}{\gamma }}\Gamma ^{-1}\right)^{-1}\Delta X\right)\right\}}$

where k_B is the Boltzmann constant and T is the absolute temperature. p(ΔY) and p(ΔZ) are expressed similarly. N-dimensional Gaussian probability density function with random variable vector x, mean vector μ and covariance matrix Σ is

${\displaystyle W(x,\mu ,\Sigma )={\frac {1}{\sqrt {(2\pi )^{N}|\Sigma |}}}\exp \left\{-{\frac {1}{2}}(x-\mu )^{T}\Sigma ^{-1}(x-\mu )\right\}}$

${\displaystyle {\sqrt {(2\pi )^{N}|\Sigma |}}}$ normalizes the distribution and |Σ| is the determinant of the covariance matrix.

Similar to Gaussian distribution, normalized distribution for ΔX^T = [ΔX₁ ΔX₂ ..... ΔX_N] around the equilibrium positions can be expressed as

${\displaystyle p(\Delta X)={\frac {1}{\sqrt {(2\pi )^{N}{\frac {k_{B}T}{\gamma }}|\Gamma ^{-1}|}}}\exp \left\{-{\frac {1}{2}}\left(\Delta X^{T}\left({\frac {k_{B}T}{\gamma }}\Gamma ^{-1}\right)^{-1}\Delta X\right)\right\}}$

The normalization constant, also the partition function Z_X, is given by

${\displaystyle Z_{X}=\int _{0}^{\infty }\exp \left\{-{\frac {1}{2}}\left(\Delta X^{T}\left({\frac {k_{B}T}{\gamma }}\Gamma ^{-1}\right)^{-1}\Delta X\right)\right\}d\Delta X}$

where ${\displaystyle {\frac {k_{B}T}{\gamma }}\Gamma ^{-1}}$ is the covariance matrix in this case. Z_Y and Z_Z are expressed similarly. This formulation requires inversion of the Kirchhoff matrix. In the GNM, the determinant of the Kirchhoff matrix is zero, hence calculation of its inverse requires eigenvalue decomposition. Γ⁻¹ is constructed using the N-1 non-zero eigenvalues and associated eigenvectors. Expressions for p(ΔY) and p(ΔZ) are similar to that of p(ΔX). The probability distribution of all fluctuations in GNM becomes

${\displaystyle P(\Delta R)=p(\Delta X)p(\Delta Y)p(\Delta Z)={\frac {1}{{Z_{X}}{Z_{Y}}{Z_{Z}}}}\exp \left\{-{\frac {3}{2}}\left(\Delta X^{T}\left({\frac {k_{B}T}{\gamma }}\Gamma ^{-1}\right)^{-1}\Delta X\right)\right\}}$

For this mass and spring system, the normalization constant in the preceding expression is the overall GNM partition function, Z_GNM,

${\displaystyle Z_{GNM}={Z_{X}}{Z_{Y}}{Z_{Z}}={(2\pi )^{3N/2}{\Biggl |}{{\frac {k_{B}T}{\gamma }}{\Gamma ^{-1}}}{\Biggr |}}^{3/2}}$

The expectation values of residue fluctuations, <ΔR_i²> (also called mean-square fluctuations, MSFs), and their cross-correlations, <ΔR_i · ΔR_j> can be organized as the diagonal and off-diagonal terms, respectively, of a covariance matrix. Based on statistical mechanics, the covariance matrix for ΔX is given by

${\displaystyle <\Delta X\cdot \Delta X^{T}>=\int \Delta X\cdot \Delta X^{T}p(\Delta X)d\Delta X={\frac {k_{B}T}{\gamma }}\Gamma ^{-1}}$

The last equality is obtained by inserting the above p(ΔX) and taking the (generalized Gaussian) integral. Since,

${\displaystyle <\Delta X\cdot \Delta X^{T}>=<\Delta Y\cdot \Delta Y^{T}>=<\Delta Z\cdot \Delta Z^{T}>={\frac {1}{3}}<\Delta R\cdot \Delta R^{T}>}$

<ΔR_i²> and <ΔR_i · ΔR_j> follows

${\displaystyle <\Delta R_{i}^{2}>={\frac {3k_{B}T}{\gamma }}(\Gamma ^{-1})_{ii}}$

${\displaystyle <\Delta R_{i}\cdot \Delta R_{j}>={\frac {3k_{B}T}{\gamma }}(\Gamma ^{-1})_{ij}}$

The GNM normal modes are found by diagonalization of the Kirchhoff matrix, Γ = UΛU^T. Here, U is a unitary matrix, U^T = U⁻¹, of the eigenvectors u_i of Γ and Λ is the diagonal matrix of eigenvalues λ_i. The frequency and shape of a mode is represented by its eigenvalue and eigenvector, respectively. Since the Kirchhoff matrix is positive semi-definite, the first eigenvalue, λ₁, is zero and the corresponding eigenvector have all its elements equal to 1/√N. This shows that the network model translationally invariant.

Cross-correlations between residue fluctuations can be written as a sum over the N-1 nonzero modes as

${\displaystyle <\Delta R_{i}\cdot \Delta R_{j}>={\frac {3k_{B}T}{\gamma }}[U\Lambda ^{-1}U^{T}]_{ij}={\frac {3k_{B}T}{\gamma }}\sum _{k=1}^{N-1}\lambda _{k}^{-1}[u_{k}u_{k}^{T}]_{ij}}$

It follows that, [ΔR_i · ΔR_j], the contribution of an individual mode is expressed as

${\displaystyle [\Delta R_{i}\cdot \Delta R_{j}]_{k}={\frac {3k_{B}T}{\gamma }}\lambda _{k}^{-1}[u_{k}]_{i}[u_{k}]_{j}}$

where [u_k]_i is the ith element of u_k.

By definition, a diagonal element of the Kirchhoff matrix, Γ_ii, is equal to the degree of a node in GNM that represents the corresponding residue’s coordination number. This number is a measure of the local packing density around a given residue. The influence of local packing density can be assessed by series expansion of Γ⁻¹ matrix. Γ can be written as a sum of two matrices, Γ = D + O, containing diagonal elements and off-diagonal elements of Γ.

Γ⁻¹ = (D + O)⁻¹ = [ D (I + D⁻¹O) ]⁻¹ = (I + D⁻¹O)⁻¹D⁻¹ = (I - D⁻¹O + ...)⁻¹D⁻¹ = D⁻¹ - D⁻¹O D⁻¹ + ...

This expression shows that local packing density makes a significant contribution to expected fluctuations of residues.[5] The terms that follow inverse of the diagonal matrix, are contributions of positional correlations to expected fluctuations.

Equilibrium fluctuations of biological molecules can be experimentally measured. In X-ray crystallography the B-factor (also called Debye-Waller or temperature factor) of each atom is a measure of its mean-square fluctuation near its equilibrium position in the native structure. In NMR experiments, this measure can be obtained by calculating root-mean-square differences between different models. In many applications and publications, including the original articles, it has been shown that expected residue fluctuations obtained by the GNM are in good agreement with the experimentally measured native state fluctuations.[6][7] The relation between B-factors, for example, and expected residue fluctuations obtained from GNM is as follows

${\displaystyle B_{i}={\frac {8\pi ^{2}}{3}}<\Delta R_{i}\cdot \Delta R_{i}>={\frac {8\pi ^{2}k_{B}T}{\gamma }}(\Gamma ^{-1})_{ii}}$

Figure 3 shows an example of GNM calculation for the catalytic domain of the protein Cdc25B, a cell division cycle dual-specificity phosphatase.

Figure 4: Slow modes obtained from GNM calculations are depicted on Cdc2B catalytic domain. A. Plot of the slowest mode. B. Mapping of the amplitude of motion in the slowest mode onto protein structure. The alpha-helix nearby the catalytic site of this domain is the most mobile region of the protein along the slowest mode. Expected values of fluctuations were also highest at this region, as shown in Figure 3. The results are retrieved iGNM server. PDB ID of Cdc25B is 1QB0.

Diagonalization of the Kirchhoff matrix decomposes the conformational motions into a spectrum of collective modes. The expected values of fluctuations and cross-correlations are obtained from linear combinations of fluctuations along these normal modes. The contribution of each mode is scaled with the inverse of that modes frequency. Hence, slow (low frequency) modes contribute most to the expected fluctuations. Along the few slowest modes, motions are shown to be collective and global and potentially relevant to functionality of the biomolecules. Fast (high frequency) modes, on the other hand, describe uncorrelated motions not inducing notable changes in the structure. GNM-based methods do not provide real dynamics but only an approximation based on the combination and interpolation of normal modes.[8] Their applicability strongly depends on how collective the motion is.[8][9]

There are several major areas in which the Gaussian network model and other elastic network models have proved to be useful.[10] These include:

• Spring bead based network model: In spring-bead based network model, the springs and beads are used as components in the crosslinked network. Springs are cross-linked to represent mechanical behavior of the material and bridge molecular dynamics (MD) model and finite element (FE) model (see Figure. 5). The beads represent material mass of cluster bonds. Each spring is used to represent a cluster of polymer chains, instead of part of a single polymer chain. This simplification allows to bridge different models at multiple length scales and improves the simulation efficiency significantly. At each iteration step in the simulation, forces in the springs are applied to the nodes at the center of the beads, and the equilibrated nodal displacements throughout the system are calculated. Different from the traditional FE method for obtaining stress and strain, the spring–bead model provides the displacements of the nodes and forces in the springs. The equivalent strain and strain energy of spring–bead based network model can be defined and calculated using the displacements of nodes and the spring characteristics. Furthermore, the results from the network model can be scaled up to obtain the structural response at the macroscale using FE analysis.[11][12]
• Decomposition of flexible/rigid regions and domains of proteins [13][14][15]
• Characterization of functional motions and functionally important sites/residues of proteins, enzymes and large macromolecular assemblies [16][11][17][18][19][20][21][22][23][24][25][26]
• Refinement and dynamics of low-resolution structural data, e.g. Cryo-electron microscopy[27][28][29][30]
• Molecular replacement for solving X-ray structures, when a conformational change occurred, with respect to a known structure[31]
• Integration with atomistic models and simulations [32][33]
• Investigation of folding/unfolding pathways and kinetics.[34][35]
• Annotation of functional implication in molecular evolution [36][37]

• iGNM: A database of protein functional motions based on GNM http://ignm.ccbb.pitt.edu [39]
• oGNM: Online calculation of structural dynamics using GNM https://web.archive.org/web/20070516042756/http://ignm.ccbb.pitt.edu/GNM_Online_Calculation.htm

• Anisotropic Network Model web server http://www.ccbb.pitt.edu/anm [40]
• elNemo: Web-interface to The Elastic Network Model http://www.sciences.univ-nantes.fr/elnemo/
• AD-ENM: Analysis of Dynamics of an Elastic Network Model http://enm.lobos.nih.gov/
• WEBnm@: Web-server for Normal Mode Analysis of proteins http://apps.cbu.uib.no/webnma/home

• ProDy: An Application Programming Interface (API) in Python, that integrates GNM and ANM analyses and several molecular structure and sequence analyses and visualization tools: http://prody.csb.pitt.edu [41][42]
• HingeProt: An algorithm for protein hinge prediction using elastic network models http://www.prc.boun.edu.tr/appserv/prc/hingeprot/, or http://bioinfo3d.cs.tau.ac.il/HingeProt/hingeprot.html
• DNABindProt: A Server for Determination of Potential DNA Binding Sites of Proteins http://www.prc.boun.edu.tr/appserv/prc/dnabindprot/
• MolMovDB: A database of macromolecular motions: http://www.molmovdb.org/