Physics Biology
ComputerScience,1970-2008-Computer Science-1998.jsonl.gz.gz
Computer Science
3d5a22358a619358f86403b3905432e49562e94f
Mapping of mutation-sensitive sites in protein-like chains

In this work we have studied, with the help of a simple on-lattice model, the distribution pattern of sites sensitive to point mutations ('hot' sites) in protein-like chains. It has been found that this pattern depends on the regularity of the matrix that rules the interaction between different kinds of residues. If the interaction matrix is dominated by the hydrophobic effect (Miyazawa Jernigan like matrix), this distribution is very simple - all the 'hot' sites can be found at the positions with maximum number of closest nearest neighbors (bulk). If random or nonlinear corrections are added to such an interaction matrix the distribution pattern changes. The rising of collective effects allows the 'hot' sites to be found in places with smaller number of nearest neighbors (surface) while the general trend of the 'hot' sites to fall into a bulk part of a conformation still holds.


I. INTRODUCTION
In this paper we study how the choice of a particular Hamiltonian is responsible for the distribution pattern of sites sensitive to point mutations in a Heteropolymeric chain.
As shown in [1] using a very simple model [2], in every optimized [3] sequence there are sites at which point mutations are likely to cause misfolding of the native state (we call them 'hot' sites), while there are other sites at which point mutations have no relevant thermodynamic effect ('cold' sites), and we call intermediate sites 'warm'. As known from both experimental [4] and theoretical studies [1] usually proteins are made of few 'hot' sites while the majority of the other sites are 'cold' . Because of the thermodynamic importance of the 'hot' sites it is of general interest to investigate the principles guiding their positioning along the protein chain.
The model used in this work describes the polymer as a chain of beads in a cubic lattice, interacting through the Hamiltonian where B ij is the contact energy between the two residues situated on the i-th and j-th positions, L is the length of the chain and ∆(r i − r j ) has the value 1 if the i-th and j-th are nearest neighbors and zero if not.
In literature, different choices of the matrix B ij have been used. The results of Ref. [1,2] have been found with a matrix [5] whose elements are distributed according to a Gaussian, with average B = 0 and standard deviation σ = 0.3 (in units of kT = 0.6 kcal/mol). With these matrix elements, for every target structure it is possible to select sequences whose native state is unique, stable and kinetically accessible [6]. It has been found that a reliable condition for a site to be 'cold' or 'hot' is intimately connected to the change in the native state energy caused by a point mutation: the bigger is this difference, the more probable for this site to be a 'hot' site. Hot sites of such sequences are mainly in the bulk sites of the native conformation, but can be on the surface as well, while some bulk sites can be rather insensitive to mutations (Fig.1a).
Repeating the same calculations with a random-generated, Gaussian-distributed set of interaction energies, we have observed similar distribution patterns for 'hot' and 'cold' sites ( Fig.1b).
Another choice of the interaction matrix can be made to take into account explicitly the hydrophobic effect encountered in real proteins. The simplest way is to choose the matrix B ij to be composed of only three different elements, namely B HH , B P P and B HP = B P H , where B HH < B P P < B HP . These elements are responsible for the interaction between hydrophobic (H) and polar (P) residues. Unfortunately it was found [2] that for a given target structure and 'HP' interaction matrix it is difficult to construct optimal sequences for which this target structure would be kinetically accessible. As a rule, the optimization of the sequence puts H-residues in the bulk sites (the sites with the greatest number of nearest-neighbors), and every substitution of an H residue with a P residue causes the chain misfolding (See Figs. 2 and 3 in Ref. [7]).
On the way between the random matrix and the highly regular HP interaction model stands the matrix deduces by Miyazawa and Jernigan in [8]. This matrix contains 210 different elements which can still be grouped into 3 big blocks, according to their hydrophobicity.
In this case it is possible to find sequences for which the native state is both stabile and kinetically accessible. It was also found that, as in the situation with only two kinds of residues, 'hot' sites are invariably the sites with the highest number of contacts (bulk sites, The question again is what are the principles guiding the positioning of 'hot' sites along the protein chain and how a particular form of the interaction matrix can influence the distribution pattern of 'hot' and 'cold' mutation sites. This paper is organized as following.
In the next section we will present a convenient representation of an interaction matrix as a function of 'mixing' parameter β. Variation of the 'mixing' parameter will correspond to a change from highly ordered HP like interaction matrix at β = 0 to highly nonlinear matrix at non-zero values of 'mixing' parameter. A distribution pattern of 'cold' and 'hot' sites will be investigated with the use of these interaction matrices. Next we will address the same question by introducing a parameter of 'randomness' γ which would allow us to investigate the mutation sites distribution pattern with highly ordered HP like interaction matrix at γ = 0 and random Gaussian interaction matrix at large γ. Conclusions will be drawn at the end.

II. LTW PARAMETERIZATION
In the work by Lee, Tang and Wingreen [9] a particular interesting parameterization of semi-experimental MJ matrix [8] was introduced, as a consequence of its regularity. In their work it was shown that elements of the MJ interaction matrix can be very nicely fitted as where q α is ascribed to a monomer of kind α and β is a constant. This parameterization involves 20 parameters, instead of 210 parameters of the MJ matrix, each specifying the strength of a particular residue. Using this parameterization we have found that, for the best fit, all q α are negative and range from −2.3 to 0, while β = −0.42. Furthermore, it has been pointed out that the origin of the additive term is due to the hydrophobic effect, while the second order term is responsible for the segregation of dissimilar residues.
Using this parameterization, we can write the Hamiltonian of the chain in a form [10] where n and q have dimensions equal to the number of monomers in the chain. The ith coordinate of n is the number of nearest neighbors for the i-th monomer and the i-th coordinate of q specifies the strength of the residue in the i-th site. C is the contact matrix for a given conformation. In the above reference it has been shown that the Hamiltonian (3) fits very well the original MJ Hamiltonian and, what is more precious, it is very convenient to handle analytically.
One of the problems that can be solved easily using the above Hamiltonian is, given a target structure (C matrix and n) and the composition of the chain in term of residues, to find the sequence which minimize the energy [10]. Particularly, using the fact that the second order term is usually 2−3 times smaller than the first order term, a first order approximation solution can be found minimizing only the solvent exclusion term. The straightforward way is to choose a sequence so that the vectors q and n are as anti-parallel as possible, keeping the constraint of a fixed number of different kinds of monomers. Knowing that all components of q are negative while all components of n are positive, it is necessary to put the residues with the most negative value of q i in the sites of the target structure with the largest number of nearest-neighbors. The effect is, roughly speaking, to put 'hydrophobic' residues (i.e. low q i ) inside the structure while keeping 'hydrophilic' ones (high q i ) on the surface. The second order term in the Hamiltonian (3) is responsible for a fine tuning of residue distribution, mostly inside the hydrophobic/hydrophilic regions, and causes a further decrease of the sequence energy.

III. MUTATIONS AND HOT SITES
According to [1], we label each mutation with the difference of the native state energy of the wild-type sequence and of the mutated sequence where B 0 ij is the interaction element associated with the wild-type sequence and B α ij with a mutated sequence. It has been shown that the information about the thermodynamic features of the mutated sequence is mostly contained in the value of ∆E loc [6], [11].
The energetic effects of mutations in a given site i are studied introducing the average ∆E loc over the 19 possible mutations in this site, For optimized sequences there are few sites (up to 10%) where ∆E loc is large (∼ ∆, where ∆ is the gap between the native state and the random conformations energy [2]), while for the others ∆E loc is much lower. We call the former sites 'hot' and the latter 'cold'. If a mutation occurs in a 'hot' site, there is a high probability that it fills the gap, eliminating the feature of design and causing misfolding of the chain.
For random interaction matrices whose elements have Gaussian distribution, the behavior of hot sites is the same as indicated in Ref. [1]. Hot sites can be both in the bulk and on the surface of the native configuration (See Fig. 1b) with dominance of 'hot' sites in the bulk.
On the other hand, in a model with only two kinds of residues (H and P), the optimization of the sequence puts H residues in the sites with the highest number of nearest neighbors.
In this case (see Ref. [7]), every site containing a H residue is a hot site.
We are now interested in studying the pattern of hot sites in a model where the interaction matrix contains both features of hydrophobic-hydrophilic separation and randomness or complicated non-linearity.

IV. CONSTRUCTION OF THE MAP OF 'HOT' AND 'COLD' SITES
To investigate the behavior of ∆E loc (i) we use the parameterized Hamiltonian (3). It is straightforward to show that, for a wild-type sequence characterized by q where < q > α = 1/20 20 α=1 q α , is the average of the values q α corresponding to the 20 monomers different from the wild-type. With this assumption we have found < q > α ≈ −1.23.
We shell consider first the case β = 0, which is exactly solvable, and then the consequences of the non-linear term introduced via non-zero 'mixing' parameter β..
In the case β = 0 the Hamiltonian contains only the solvent exclusion term H 0 = i q i n i .
As we discussed earlier, optimization of the sequence, given a structure, consists in putting the most hydrophobic residues into the sites with the greatest number of contacts, keeping with the constraint of fixed number of different residues, so that For short chains and 20 letter code used in the model with β = 0 the optimization procedure described above works reasonably good in constructing optimal sequences while for longer chains it seems that the model is not adequate to ensure a single ground state or a ground state well-separated from the other states. Nevertheless, the case β = 0 is a good starting point to investigate the role of the different terms of Hamiltonian (3).
In this case equation (6) can be written as Defining < q >= L i=1 q i / L i=1 n i and substituting (7) in (8) as an approximation for an optimal sequence, we obtain For typical values of < q >∼ −0.6 and < q > α ∼ −1.26 then The shape of this function is plotted in Fig. 2. It is interesting to compare the value ∆E loc (i) for the optimized sequence with the one for a random sequence. We will identify, in the spirit of Random Energy Model [12], 'hot' sites (as defined in Sect. III) with those sites in which the average impact of mutations is greater than for a random sequence. We define a sequence to be random if there is no correlation between the strength q i at a given site, and its number of nearest neighbors n i , so that where q ranges between −2.3 and 0. The values that ∆E rand loc (i) can assume are comprised between the two straight lines plotted in Fig.2a. While the dependence of ∆E loc (i) for the random sequence is linear, for a selected sequence it is quadratic. In the case of the selected sequence the quadratic behavior of the mutation energy versus the number of closest nearest neighbors induces a sharp distinction between bulk sites (n > ∼ 3) and surface sites. As clear form (Fig. 2a), all bulk sites are 'hot', while surface sites are 'cold'. As consequence, 'hot' sites have a certain degree of symmetry in target structures, i.e. no one of the sites with the same number of nearest neighbors is privileged to the others. The map of the 'hot' and 'cold' sites for a Hamiltonian β = 0 and a 36 monomer target structure is presented on (Fig.   1c).
If β is not zero, but still small in absolute values, it is possible to minimize the energy of the sequence in two steps, first minimizing H β=0 = n q and finding an initial trial optimal sequence q tr , and then re-minimizing the same H β=0 with an effective n ef f = n + β 2 q tr C. As was shown in Ref. [10] this procedure is quite reliable for small β.
Using the same approximation (7) as in the previous section for q tr one finds where N i = j C ij n j is the number of nearest neighbors of the nearest neighbors of site i and can assume values in the range {n i , ..., 4n i + 1}. The expression of ∆E β loc (i), then, takes into account, through the value of N i , sites lying further than the nearest neighbors of i.
As consequence of this, we can observe (Fig. 2b) a broadening of the range of values that ∆E β loc (i) can assume for each n i (we will refer to this broadening as 'energy bands' in our further discussion) can be broken thus allowing some hot sites to be found on surface.

C. Explicit Calculations with MJ parameterized Hamiltonian
To investigate how the 'energy bands' depend upon the strength of the nonlinear contribution in the interaction matrix we have made some explicit calculations, using as target structure the 36mers chain displayed on (Fig. 1). We have first calculated the best values of q α and β to fit MJ matrix according to Eq. 2, finding β = −0.42. Using these values for q α , we varied β in the range (−1.5, 1.5), optimizing each time the sequence with a genetic energy minimization technique. The composition has been kept fixed for all values of β and chosen in such a way as to satisfy q i ∼ n i (condition for optimal composition at β = 0). For the sake of computational convenience, the interaction matrix elements have been rescaled to have zero average and standard deviation equal to 1.0. Then, we have plotted the value of ∆E β loc (i) for each lattice site, as function of β. First, we consider the case β < 0. The raising of different 'energy bands' is shown in Fig.3a. For −0.7 < ∼ β < 0 all the bands lay into four distinct groups. The first two groups, which correspond to the sites with n i = 1 or n i = 2, contribute to cold sites. The other two groups correspond to the sites with n i = 3 or n i = 4 defining the 'warm' and 'hot' sites. This situation is very similar to the case of β = 0, where bands do not overlap and 'hot' sites are exclusively in the bulk. It is not surprising to find the same distribution pattern for the MJ matrix (β = −0.42) (Fig.1d).
and HP-like interaction matrix (β = 0) (Fig.1c). For −2.5 < ∼ β < ∼ −1.0 the nonlinear part of the interaction matrix starts playing its role. As was shown in the previous section the second neighbor shell contribute to the value of ∆E β loc (i). At these values of parameter β we observe that some of the energy bands corresponding to 'warm' (n=3) and 'cold' sites (n=1,2) mix, while the bands corresponding to the 'hot' sites stay well separated from the other bands.
In the case of β > 0 (Fig.3b) the non-linear effect is much more dramatic. For β > ∼ 0.2 all bands start mixing allowing 'cold' sites penetrate the bulk while pushing 'hot' sites on a surface. We can rationalise the different pattern of 'energy bands' at β < 0 and β > 0 by examining the Eq.6 and Hamiltonian 3. Noticing that β 2 qC q is of the order of β 2 q n < q > we can rewrite the Hamiltonian 3 in a form where η contains the collective effects and q n(1 + β 2 < q >) describes a 'renormalized' hydrophobic effect. In the case β < 0, ( β 2 < q >) > 0 and the interaction matrix is largely dominated by the hydrophobic effect, while the collective contribution to the Hamiltonian is not strong enough to substantially mix the 'energy bands'.
For the β > 0 case ( β 2 < q >) < 0 and the renormalized hydrophobic effect becomes comparable or smaller than the collective term in the Hamiltonian. This allows all bands to mix substantially.

D. Explicit Calculations for Random Hamiltonian
Another interesting question which can be addressed is how our conclusions are modified by the addition of a random term in the β = 0 Hamiltonian. To study this problem, we have chosen a Hamiltonian in the form where n q comes from the parameterization (2)  In the inserts of Figs.3a,b and Fig.4 shown the energy of the optimal non-mutated sequence for different values of β and γ, respectively. It is possible to observe a sudden decrease in energy of the optimal sequence as β and γ increases signifying that the addition of non-linear terms into a Hamiltonian allows for a much better energy minimization.

V. CONCLUSIONS
In this work we have considered how the regularity of the interaction matrix influence the distribution pattern of 'hot' sites. It has been found that if the matrix is polarized (dominant hydrophobic effect, β ∼ 0 or γ ∼ 0) this distribution is very simple. All the 'hot' sites can be found at the positions with maximum number of closest nearest neighbors (bulk).
With increasing importance of non-linear terms (β = 0,γ = 0) the distribution pattern changes so that the 'hot' and 'warm' sites can be found in places with smaller number of nearest neighbors (surface) while the general trend of the 'hot' sites to fall into a bulk part of a conformation still holds.
As pointed out above this can be rationalized by noticing that if mixing parameter is different than zero each site starts feeling not only its nearest neighbors but also the more distant sites. This leads to a collective nature of the interactions giving rise to a modified distribution pattern of 'hot' sites. Energy of optimal sequence