Graphical models for protein structure

Let X = {X_b, X_s} be the random variables representing the entire protein structure. X_b can be represented by a set of 3-d coordinates of the backbone atoms, or equivalently, by a sequence of bond lengths and dihedral angles. The probability of a particular conformation x can then be written as:

${\displaystyle p(X=x|\Theta )=p(X_{b}=x_{b})p(X_{s}=x_{s}|X_{b},\Theta ),\,}$

where ${\displaystyle \Theta }$ represents any parameters used to describe this model, including sequence information, temperature etc. Frequently the backbone is assumed to be rigid with a known conformation, and the problem is then transformed to a side-chain placement problem. The structure of the graph is also encoded in ${\displaystyle \Theta }$. This structure shows which two variables are conditionally independent. As an example, side chain angles of two residues far apart can be independent given all other angles in the protein. To extract this structure, researchers use a distance threshold, and only a pair of residues which are within that threshold are considered connected (i.e. have an edge between them).

Given this representation, the probability of a particular side chain conformation x_s given the backbone conformation x_b can be expressed as

${\displaystyle p(X_{s}=x_{s}|X_{b}=x_{b})={\frac {1}{Z}}\prod _{c\in C(G)}\Phi _{c}(x_{s}^{c},x_{b}^{c})}$

where C(G) is the set of all cliques in G, ${\displaystyle \Phi }$ is a potential function defined over the variables, and Z is the partition function.

To completely characterize the MRF, it is necessary to define the potential function ${\displaystyle \Phi }$. To simplify, the cliques of a graph are usually restricted to only the cliques of size 2, which means the potential function is only defined over pairs of variables. In Goblin System, these pairwise functions are defined as

${\displaystyle \Phi (x_{s}^{i_{p}},x_{b}^{j_{q}})=\exp(-E(x_{s}^{i_{p}},x_{b}^{j_{q}})/K_{B}T)}$

where ${\displaystyle E(x_{s}^{i_{p}},x_{b}^{j_{q}})}$ is the energy of interaction between rotamer state p of residue ${\displaystyle X_{i}^{s}}$ and rotamer state q of residue ${\displaystyle X_{j}^{s}}$ and ${\displaystyle k_{B}}$ is the Boltzmann constant.

Using a PDB file, this model can be built over the protein structure. From this model, free energy can be calculated.

It has been shown that the free energy of a system is calculated as

${\displaystyle G=E-TS}$

where E is the enthalpy of the system, T the temperature and S, the entropy. Now if we associate a probability with each state of the system, (p(x) for each conformation value, x), G can be rewritten as

${\displaystyle G=\sum _{x}p(x)E(x)-T\sum _{x}p(x)\ln(p(x))\,}$

Calculating p(x) on discrete graphs is done by the generalized belief propagation algorithm. This algorithm calculates an approximation to the probabilities, and it is not guaranteed to converge to a final value set. However, in practice, it has been shown to converge successfully in many cases.

Gaussian graphical models are multivariate probability distributions encoding a network of dependencies among variables. Let ${\displaystyle \Theta =[\theta _{1},\theta _{2},\dots ,\theta _{n}]}$ be a set of ${\displaystyle n}$ variables, such as ${\displaystyle n}$ dihedral angles, and let ${\displaystyle f(\Theta =D)}$ be the value of the probability density function at a particular value D. A multivariate Gaussian graphical model defines this probability as follows:

${\displaystyle f(\Theta =D)={\frac {1}{Z}}\exp \left\{-{\frac {1}{2}}(D-\mu )^{T}\Sigma ^{-1}(D-\mu )\right\}}$

Where ${\displaystyle Z=(2\pi )^{n/2}|\Sigma |^{1/2}}$ is the closed form for the partition function. The parameters of this distribution are ${\displaystyle \mu }$ and ${\displaystyle \Sigma }$. ${\displaystyle \mu }$ is the vector of mean values of each variable, and ${\displaystyle \Sigma ^{-1}}$, the inverse of the covariance matrix, also known as the precision matrix. Precision matrix contains the pairwise dependencies between the variables. A zero value in ${\displaystyle \Sigma ^{-1}}$ means that conditioned on the values of the other variables, the two corresponding variable are independent of each other.

To learn the graph structure as a multivariate Gaussian graphical model, we can use either L-1 regularization, or neighborhood selection algorithms. These algorithms simultaneously learn a graph structure and the edge strength of the connected nodes. An edge strength corresponds to the potential function defined on the corresponding two-node clique. We use a training set of a number of PDB structures to learn the ${\displaystyle \mu }$ and ${\displaystyle \Sigma ^{-1}}$.

Once the model is learned, we can repeat the same step as in the discrete case, to get the density functions at each node, and use analytical form to calculate the free energy. Here, the partition function already has a closed form, so the inference, at least for the Gaussian graphical models is trivial. If the analytical form of the partition function is not available, particle filtering or expectation propagation can be used to approximate Z, and then perform the inference and calculate free energy.