Protein intrachain contact prediction with most interacting residues (MIR)

Lattice geometry

We model a protein as a chain of evenly spaced C_α atoms placed on a lattice [14, 45]. We define a lattice unit (lu) to be 1.7 Å. Hence, C_α atoms are connected by vectors of the form (2, 1, 0). These vectors are 5 lu in length, which corresponds to 3.8 Å, the mean distance between adjacent C_α atoms in proteins. This results in 24 immediate neighbor positions for each point in the lattice. This represents the intersection of a 4×4-segmented cube with a sphere of radius 3.8 Å (5 lu), as shown in Figure 10.

Figure 10

Vectors resulting from intersection of lattice with sphere at origin 0.0.

Our model does not take into account the presence of side chains; therefore, the required separation is modeled with a 3.8 Å minimum distance requirement. On the basis of chain geometry, we limit the angle between some C_αs at positions i, i+1, and i+2 in a sequence by requiring the distance between them to be from 4.1 to 7.2 Å (or from 6 to 18lu). This corresponds to angles from 66° to 143° [14, 43], which are closer to the real angles in α-helices and β-strand conformations than previous cubic lattice methods [5]. This is demonstrated in Figure 11. Here, some residue i is fixed at (0, 0, 0); we then show all 24 possible positions for residue i+1 (black vectors). For each i+2 residue, after excluding the occupied origin, there is a choice of 23 possible vectors: for clarity, only one (0, 1, 2) (the green and red vectors) is shown. Red vectors are those that violate the distance (angle) restriction and will not be permitted by the method.

Figure 11

Subset of i to i+2 residue positions.

Angle restriction: vectors parallel or producing sharp corners, thus violating the angle constraint, are shown in red. The sixth invalid position is not visible because it overlaps with the origin.

To initiate the simulations, 100 different starting conformations (or models) within the lattice are used. Figure 12 displays a sample of these models as a comprehensive plot. These starting conformations were taken from offline computations for chains of 1100 residues. The only requirement is that these randomly computed seed conformations have some level of non-compactness [14, 45]. Starting from the first backbone residue located at position (0, 0, 0), the first n positions in the seed model will be used for an input model with n residues in its sequence. When placing an input protein into the lattice for the first time, each residue in the protein is positioned based on the respective residue in the seed model, i.e., the first residue of the given protein is assigned the position of the first residue in the seed model, the second residue of the given protein is assigned the position of second residue in the seed model, and so on. If the input sequence is shorter than the seed model, then the positions of the final residues in the seed will not be used.

Figure 12

Left: first five initial models. Right: all 100 initial models.

As will be discussed later, a number of possible perturbations may be applied to each of the residues in the model. The model is stored as a set of relative vectors between C_αs, representing their distance relative to the previous C_α. This is useful as we may, before knowing that a move is valid, change the position of a single residue without having to update every other residue that follows it in the chain. For the purposes of energy and neighbor calculations, a displacement vector must be calculated for the absolute position of each residue. We create new conformations by working along the relative positions of residues and accumulating their positions in space, in a neighbor-to-neighbor fashion. Hence, every residue following the one being changed will be translated by the distance between the initial and final positions of the perturbed residue. The resulting models are verified by checking that no non-adjacent residues have come into close proximity (5 lu) with other residues (which should be disallowed by our move set) and that no invalid bond angles have formed. If these checks fail, the perturbation that proposed the new conformation will be abandoned.

At the end of each simulation, the algorithm internally produces a model that has collected all of the perturbations that were accepted, based on an energy criterion described later. It then uses the topology of the model to determine the NCN count at each residue. Two of these resulting models are shown in Figure 13 for the PDB code 1asu. Both of the MIR models follow the rough pattern of the starting model; however, one can observe the same globular regions starting to form at the same places (e.g., ends). In previous works, these were called proto fragments [5].

Figure 13

Models resulting from the first two simulations for the protein with PDB code 1asu.

The blue models are the results, and the green model is seed model used.

Energy model

Although the protein chain is only modeled as a sequence of C_αs, the effect of side chains is included in energy terms associated with each pair of residue interactions. We assume that inter-residue energies are significant when the distance between them is between 3.8 and 5.88 Å (5 to 12 lu). The lattice model requires a minimum distance of 3.8 Å for side-chain separation. For this work, we use the distance-independent statistical pair potential by Miyazawa and Jernigan [4, 46]. This is a 20×20 symmetric matrix where solvent effects are implicit.

We take ER(i) to be the energy at the i^th residue in a sequence and calculate it as

ER(i)=∑j>i+1j<i−1 orEI(i, j),

where the energy interaction EI is calculated according to the distance between residues (i, j) and the energy matrix PE corresponding to the residue-residue energy interaction [4, 14] that is a function of the type of residue (one of the 20 amino acids). We use type(i) to denote a function from residue index to residue type. Let dist(i,j) compute the Euclidean distance between two points on the lattice.

EI(i, j)={PE (type(i), type(j)), if 5≤dist(i, j)≤12 0, otherwise

Monte Carlo algorithm

The core of the MIR algorithm is a Monte Carlo simulation where possible perturbations (moves) are repeatedly applied to an existing conformation and accepted on the basis of standard Metropolis criterion. In the following sections, we detail the concrete MIR implementation. The algorithm is implemented in Fortran 90 and is used in a Linux environment. All random numbers are computed with the “Keep It Simple, Stupid” method [47]. Each time the algorithm is run, the same initial seed value is used; together with the precomputed random initial confirmation, this enables reproducible results. The Monte Carlo simulation is used to generate a total of 100 models.

The limiting number of Monte Carlo steps for each simulation is given by:

MClimit (L)={106,106(L50)2, if L<50    otherwise

where L is the sequence length. During a simulation, we record the state (snapshot) of residue interactions every MChop=MClimit104 steps for a minimum of every 10² states. The same process is followed for each simulation: first, we calculate the MClimit value. We then calculate how often an interaction snapshot should be taken. The overall initial energy of the model is then computed. From this point, we run the main simulation function to create new conformations until the model has performed MClimit number of steps. Let E be the sum of all residue energies for an entire model. For each change to the existing model, we compute ∆E=E_new–E_existing.

A new conformation is accepted if ΔE<0; otherwise, it may be accepted with probability b=e−ΔERT where R is the gas constant and T is the temperate such that RT=1.5, corresponding to the optimized test previously performed [5]. According to [6] acceptance happens with probability: b1+b. Regardless of its acceptance, each attempted conformation represents an MC step. For each position, the number of NCNs is periodically recorded (every MChop steps). After the simulations have completed, we calculate the NCN count for each residue as an average based on the number of snapshots (104). Let N(i, j) be the total number of times i and j are non-covalent neighbors over all snapshots. Thus:

NCN(i)=1 104∑j<i−1 or j>i+1N(i, j)

Residues with NCN of at least 6 are marked as MIRs; any residues with NCN of no more than 2 are marked as LIRs. To reduce statistical fluctuations that can produce successive positions attributed as an MIR, but without physical meaning, a smoothing procedure is implemented on the web server of SPROUTS. On the basis of a Pascal algorithm, it produces a smoothed distribution of NCNs, and the maxima are then considered as SMIRs.

Simulation step

A model is evolving while we have not reached MClimit and can still select perturbations, which change energy between some residues. When we begin processing the model, we randomly select an unseen residue in the sequence to perturb. Call it residue i. We calculate the angle between residues i–1, i, and i+1. We limit the angle between these residues by restricting the distance between them from 4.1 to 7.2 Å (or from 6 to 18). On the basis of the angle between the residues adjacent to i, a number of different perturbations (or moves) may be possible. Below, we will consider one of the six groups of angle perturbations that are possible. In the illustration of Figure 14, we simplify the possible placements by assuming that residue i-1 is located at the origin of the coordinate system and that residue i is located at (0, 1, 2). The following perturbations were originally presented by Skolnick and coworkers [6–8]. We use the term perturbation vectors to refer to the use of the difference vectors defined by the 24 possible positions between covalently bonded residues in the lattice.

Figure 14

Possible perturbations for residue i at position (0, 1, 2).

Residue i-1 is always located at (0, 0, 0), which is the bottom intersection of the bright colored lines. Residues i+1 are located at the ends of the top black vectors. Green are 10, brown is 14, and orange are 18.

For a simple example of a perturbation, consider those perturbations of length 10, 14, 18. These three distances each define a family of vectors (along some rotational symmetry) that can each serve as perturbation of the model. For these so-called corner moves, the perturbation takes the form of an exchange of the relative vectors between residues i-1 to i with i to i+1 (see Figure 1; the continuous and dashed lines are swapped). Figure 14 shows the results of these perturbations to an arbitrary position for residue i (midpoint of bold lines). The nine positions possible, after removing unnatural bond angles, for residue i+1 are indicated in black. For each possibility from these nine initial positions, the dashed green lines indicate the associated result of the perturbation. Notice that for each case, residues i–1 and i+1 do not change position.

For distances of 0, 6, 8, 12, 16, the case is more complex: the algorithm locates multiple position perturbations. For 0, i.e., the N-terminal residue, the algorithm selects any of the 23 permutation vectors that do not overlap with the model (the so-called end move). For other distances, the algorithm also randomly selects one of the generated perturbation vectors that do not overlap with the model. However, it will have fewer valid choices than for the residues at the limits. The process of randomly selecting a perturbation continues until one that generates a valid model is found. When a position is problematic (i.e., it generates a model with residue overlap or unnatural angles), it is marked as seen so that it is not retried during the particular simulation step.

If the new model is valid, the change is applied and the energy value is calculated for the residue that was moved. This local energy result is passed to the energy acceptance function to probabilistically determine if this new model should be accepted. If the model is invalid or is not accepted, then we reset to the previous model. In this case, one begins the process of seeking a residue to perturb again.