Genet. Mol. Res. 6 (4): 799-820 (2007)
Received August 03, 2007
Accepted September 25, 2007
Published October 05, 2007

ABSTRACT. The construction of a realistic theoretical model of proteins is determinant for improving the computational simulations of their structural and functional aspects. Modeling proteins as a network of non-covalent connections between the atoms of amino acid residues has shown valuable insights into these macromolecules. The energy-related properties of protein structures are known to be very important in molecular dynamics. However, these same properties have been neglected when the protein structures are modeled as networks of atoms and amino acid residues. A new approach for the construction of protein models based on a network of atoms is presented. This method, based on interatomic interaction, takes into account the energy and geometric aspects of the protein structures that were not employed before, such as atomic occlusion inside the protein, the use of solvation, protein modeling and analysis, and the use of energy potentials to estimate the energies of interatomic non-covalent contacts. As a result, we achieved a more realistic network model of proteins. This model has the virtue of being more robust in face of different unknown variables that usually are arbitrarily estimated. We were able to determine the most connected residues of all the proteins studied, so that we are now in a better condition to study their structural role.

Key words: Protein structure, Complex networks, Myoglobins

INTRODUCTION

The principles of protein folding that result in a balance of stability and flexibility, while maintaining its function, are not perfectly understood and have been difficult to exploit for the development of thermo-stabilized proteins. A key mechanism for thermo-stabilization appears to be the optimization of interactions between atoms within a protein. The underlying principles of protein stability and folding have been studied by a variety of analyses of a large number of different protein structures. Some theoretical studies of protein structures and other empirical methods have been used to understand the stability of proteins. Hitherto, different kinds of experiments have been carried out to understand the stability of proteins, and to determine whether any specific residues exist that could be identified as playing a relevant role in the phenomenon (Fersht and Daggett, 2002; Onuchic and Wolynes, 2004).

This study focuses on understanding the principles of protein structure by considering them as complex networks of noncovalent interactions. We adopt a novel approach on how protein structures could be modeled, with the energy of atomic interactions and the spatial occlusions among them playing important roles in determining the characteristics of the network.

In this study, we also showed the need of solvation of the protein before the analysis of the atomic interactions, in order to avoid unrealistic results when no cut-off value is used to limit the distance between the atoms when their interaction is calculated. Concomitantly, the residues exposed to the solvent tend to show an excessive number of links established with other atoms on the surface of the protein. The solvation of the protein mimics more realistically the scope where this kind of interaction happens. Therefore, the number of links among the atoms decreases drastically and the mathematical patterns of the resulting network become more realistic.

COMPLEX NETWORK THEOTETICAL REVIEW

In recent years, some relevant studies have applied the principles of networks for modeling complex systems, such as in epidemiological processes, social processes, metabolic networks, micro-electronic devices, etc. For networks of tens or hundreds of vertices, it is a relatively straightforward matter to draw a picture of the network with actual points and lines and to answer other specific questions about diverse kinds of network structure by examining this picture. There are some (Newman, 2003) advances towards providing answers to questions related to characterizing and modeling the network structure. On the other hand, studies of the structural effects on system behavior are shown to be still incipient.

A central concept when analyzing complex networks is its topology, and two different topological models have been adopted: "small-world" and "power law" models. For instance, the following concept of "small world" topology is adopted in this study: If the number of vertices within a distance r of a typical central vertex grows exponentially with r, and this is true for many networks including the random graph, then the value of the mean geodesic distance between all pairs of vertices l that have a connecting path will increase by log n. The term "small-world effect" has thus taken on a more precise meaning: Networks are said to show the small-world effect if the value of l scales logarithmically or more slowly with network size for a fixed mean degree. Logarithmic scaling can be proved for a variety of network models and has also been observed in various real-world networks. Some networks have mean vertex - vertex distances that increase more slowly than by log n. Bollobás et al. (2001) have shown that networks with power-law degree distributions have values of l that increase no faster than log n/log log n.

A small-world network is one that has a relatively short characteristic path length, L, and a high clustering coefficient, C. Small-world models can be built on lattices of any dimension or topology, but the best-studied case by far is the one-dimensional one. If we take a one-dimensional lattice of L vertices with periodic boundary conditions, i.e., a ring, and join each vertex to its neighbors k or smaller lattice spacing away, we get a system, with Lk edges (Newman, 2003). The small-world model is then created by taking a small fraction of the edges in this graph and "rewiring" them (Newman, 2003). The rewiring procedure involves going through each edge in turn and, with probability p, moving one end of that edge to a new location chosen uniformly at random from the lattice, except that no double edges or self-edges are ever created. Furthermore, these features will be important to characterize the structural aspects of the proteins studied. In the next sections, other important features of complex networks will be explained.

Transitivity and clustering

In the present study, we adopt the concept of clustering C defined as: if a vertex v has kv neighbors, then the maximum number of links between these neighbors is [k_v (k_v - 1)]/2.

The term C_v gives the fraction of these possible links that actually exist, and C is then defined as the average C_v over all vertices v. C is a measure of local clustering, which means that if two vertices X and Y are both connected to a third, Z, then for large C there is a high probability that X and Y are also directly linked to one another.

In terms of network topology, transitivity means the presence of a heightened number of triangles in the network (sets of three vertices each of which is connected to each of the others). It can be quantified by defining a clustering coefficient C (Newman, 2003):

where a "connected triple" means a single vertex with edges running to an unordered pair of others. In simple terms, C is the mean probability that two vertices that are network neighbors of the same other vertex will themselves be neighbors.

The clustering coefficient measures the density relationships in a network. An obvious generalization is to also ask about the density of longer loops, that is, loops of length four and higher. The clustering coefficient may identify densely connected groups of vertices and the "hubs" inside these structures. Therefore, the clustering coefficient may provide evidence of a modular view of the network’s dynamics (Holme et al., 2003; Guimerá and Nunes Amaral, 2005).

Degree distributions

The degree of a vertex in a network is defined as the number of edges incident on (i.e., connected to) that vertex (Newman, 2003). We define p_k to be the fraction of vertices in the network that have degree k. Equivalently, it should be expressed:

Although in theory one has just to construct a histogram of the degrees, in practice one rarely has enough measurements to get good statistics in the tail, and direct histograms are thus usually rather noisy. There are two accepted ways to get around this problem. One is to construct a histogram in which the bin sizes increase exponentially with degree. This method of constructing a histogram is often used when the histogram is to be plotted with a logarithmic degree scale, so that the widths of the bins will appear even. Because the bins get wider as we get out into the tail, the problems with statistics are reduced, although they are still present to some extent as long as p_k falls off faster than k^-1, which it must if the distribution is to be integrable.

An alternate way of presenting degree data is to make a plot of the cumulative distribution function

which is the probability that the degree is greater than or equal to k. Such a plot has the advantage that all the original data are represented. When we make a conventional histogram by binning, any differences between the values of data points that fall in the same bin are lost. The cumulative distribution function does not suffer from this problem. The cumulative distribution also reduces the noise in the tail.

The degree distribution of a network is interesting as a measure to identify the most connected vertices (i.e., the "hubs") of the network which provide the connectiveness among the existing communities inside the network.

Summarizing

All the real network topologies discussed in the literature so far - chains, grids, lattices and fully connected graphs - have varied from a completely regular lattice to a quite random lattice. Newman (2003) exemplifies three different types of lattice: regular, random and scale-free. These simple models allow us to focus on the complexity caused by the nonlinear dynamics of the nodes, without being burdened by any additional complexity in the network structure itself (Newman, 2003). Ignoring in this approach the dynamic aspects, it is possible to turn our attention to describing the structural aspects of the more complex architectures. For instance, these previously explained features of the networks will be helpful in characterizing the structure of the proteins studied.

MATERIAL AND METHODS

Data set

The data set used in this analysis consisted of structures of 12 myoglobins obtained from the protein data bank - PDB (Berman et al., 2000). This non-redundant set of proteins with a resolution better than 1.8 Å was submitted to a simulated annealing procedure at 312 K in order to achieve their physiological arrangement. The sizes of the proteins considered vary from 2300 to 2700 atoms. Traditionally, the three-dimensional folds of proteins have been perceived as a construct based on elements of secondary structure and fold arrangement. An alternative way to conceptualize and model protein structures, used in our study, is to consider the contacts between atoms, that constitute the amino acid residues, as a network of interactions irrespective of secondary structure and fold type.

Hydropathic nature of interactions

Here, the concepts presented by Sobolev et al. (1999), and implemented in STING (Higa et al., 2004), were adopted in order to distinguish the nature of interactions. For a protein, based on PDB file assignment, one must take a pair of atoms and, according to the nature of each atom in that pair and on an analysis of the interatomic distance, a class is assigned to each atom. However, other aspects must be taken into account as described by Sobolev et al. (1999).

Need of solvation

Constituent atoms of the exposed residues at the surface of the proteins tend to show an excessive number of links established with other atoms on the surface of the protein, when the molecule is considered in vacuo and no distance cut-off is used. Aiming at a solution to this problem, a solvation of the protein was done before the analysis of the atomic interactions in order to avoid unrealistic results. The solvation of the protein mimics more realistically the scope where this kind of interaction happens. In order to solvate the molecules analyzed in this study, the package SOLVATE of NAMD (Phillips et al., 2005) was used, using its default parameters.

Problem of occlusions

Physically, it is unrealistic to consider non-covalent interactions between different atoms of the protein without taking into account the occlusion due to the spatially intervening atoms.

In order to address this problem, we developed an algorithm that does not consider unrealistic interactions that could not occur between two atoms having a third (or more) atom(s) occluding one from the other. Calculation of the occlusion between atoms implies the determination of the Euclidian distance and the angle between them.

Let T = { A_i; A_j; A_k } | i ≠ j ≠ k be a set of three atoms as shown in Figure 1.

Let V_ij be the vector linking (A_i;A_j) and let V_ik be the vector linking (A_i;A_k). Let θ be the angle formed by A_j -A_i -A_k. If cos(θ) > 0 and the perpendicular from A_k to V_ij is less than the sum of the van der Waals radii of A_k and A_j, then A_j is taken to be occluded from A_i by A_k.

In this way, it is important to measure how the energy of interaction between A_i and A_j decreases as the atom A_j becomes obfuscated from A_i by A_k. A reasonable approach would be to take the energy of interaction to be proportional to the remaining exposed projected area from A_j "visible" by A_i. This approach is analogous to a "lunar eclipse", where the earth sheds its shadow over the moon. Geometrically, this rationale can be illustrated as shown in Figure 2.,

In order to better explain this concept, let the projections of two circles of radii r_k and r_j, and centered at (0; 0) and (d; 0) that intersect both of them in a region shaped like an asymmetric lens. The equations of the two circles are:

In order to find the area of an asymmetric "lens" in which the circles intersect, we may adopt an expression for the circular segment of radius R’ and triangular height d’ where:

Since there are two "lens" at the sphere intersection, the total area of intersection will be found by performing this calculation twice. For instance, the heights of the two segment triangles are

The total area of the "lens" may then be expressed as

or

The limiting cases of this expression may be checked varying d = 0 until d = r_j + r_k. As can be seen in Figure 3, the exposed area of A_j (A_f ), decreases quite linearly as A_k intersects the space between A_i and A_j. In the scope of this study, this function will be used as an attenuation factor for the energy of interaction, in order to reproduce the interferences caused by the presence of a third atom in an atomic pair interaction.

Construction of the representative graph of protein structure

The RGPSs are constructed from the three-dimensional atomic coordinates of the protein structures obtained from the PDB. Each protein in the data set is represented as a graph consisting of a set of vertices and edges. Each atom in the protein structure is represented as a vertex and these vertices are connected by edges based on the energy of non-covalent interaction between them. The energy of interaction between two atoms is given by:

where E_ij is the total energy associated with the atomic interaction between the atoms A_i and A_j (i ≠ j) such that A_i and A_j are not part of the same amino acid residue. The total energy associated with each atomic interaction has two components: the energy derived from the Coulomb potential (E^C_ij):

in kcal/mol, where: C - proportionality constant; q_i, q_j - charge is in electron-charges; ε_ij - apparent dielectric constant of the medium; r(ij) - distance in Angstroms; and the energy associated with the Lennard-Jones potential (E^LJ_ij ):

where: ε_ij - apparent dielectric constant of the medium between i and j; n - large coefficient (usually 12); m - low coefficient (usually 6); r_eqm(ij) - parameter distance for an interaction between i and j; r(ij) - present distance between i and j; as described in the AMBER98 (Morris et al., 1998) potential.

Initially, let a protein P be formed by a sequence of NR residues {R₁;R₂;R₃; … ;R_NR}. Let a residue R be formed by a sequence of NA atoms {A₁;A₂;A₃; … ;A_NA}. The total number of atoms in a protein N_AP is calculated as:

A protein P may be viewed as formed by a set of N_AP atoms {A₁;A₂;A₃; … ;A_NAP}.

The following procedure is used in order to identify the pairs of atoms (A_i, A_j) that establish an interaction:

where A_f is the exposed area of A_j.

In this procedure no arbitrary cut-off for the distance between the pair (A_i, A_j) is needed, since the most relevant criterion for establishing an interaction is the presence or not of a third atom A_k that could hide, totally or partially, A_j from A_i. Therefore, if A_j is "visible" by A_i there is an interaction and the energy of interaction is calculated.

For practical reasons, we limited our calculations to a sphere of influence with radius of 10 Å in order to avoid an unnecessary amount of computation. On the one hand, the Lennard-Jones potential is usually calculated from this limit and below. On the other hand, we studied other values beyond this limit but all pairwise interactions, where the distance between the atoms was greater than 9 Å, were occluded. In this study, the dielectric constant for calculations focusing on the atoms inside the proteins was fixed at 4.

Analysis of representative graph of protein structures

The networks are analyzed for the distribution of nodes with k links. For each RGPS, the number of nodes N(k) with k edges (links) is evaluated. The value N(k) for all proteins in the data set is taken, and then N(k) versus k is plotted.

In order to calculate some properties of the RGPS, they are represented as an adjacency matrix A^E, where:

• A^E_ij = E_ij,
  if i ≠ j and i and j are not occluded from each other;
• A^E_ij = 0,
  if i ≠ j and i and j are hidden from each other;
• A^E_ij = 0, if i = j

The adjacency matrix A^E is then analyzed using standard graph techniques to identify distinct clusters and the cluster-forming nodes (atoms) in the RGPS. The largest cluster is then identified, and its size is determined for all the RGPSs. The normalized value of the largest cluster size (with respect to the total number of residues in the protein) is plotted for all the proteins in the data set. Specifically, the high contact number atoms (N_AP > 10), and residues, will be referred as "hubs".

Edge distribution profile of atoms and residues of amino acids

For a given atom A_i, the total number of contacts established by it, may be defined as:

For a given residue R_m, the total number of contacts for this residue N_CR(m) is calculated as:

where P is the protein studied, and:

i ≠ j ^ A_i ∈ R_m ^ A_j R_m ^ {A_i;A_j} ∈ P.