Cooperative "folding transition" in the sequence space facilitates function-driven evolution of protein families
Akira R. Kinjo
































akinjo@protein.osaka-u.ac.jp

Institute for Protein Research, Osaka University, 3-2 Yamadaoka, Suita, Osaka, 565-0871, Japan



In the protein sequence space, natural proteins form clusters of families which are characterized by their unique native folds whereas the great majority of random polypeptides are neither clustered nor foldable to unique structures.
Since a given polypeptide can be either foldable or unfoldable, a kind of “folding transition” is expected at the boundary of a protein family in the sequence space. By Monte Carlo simulations of a statistical mechanical model of protein sequence alignment that coherently incorporates both short-range and long-range interactions as well as variable-length insertions to reproduce the statistics of the multiple sequence alignment of a given protein family, we demonstrate the existence of such transition between natural-like sequences and random sequences in the sequence subspaces for 15 domain families of various folds.
The transition was found to be highly cooperative and two-state-like.
Furthermore, enforcing or suppressing consensus residues on a few of the well-conserved sites enhanced or diminished, respectively, the natural-like pattern formation over the entire sequence. In most families, the key sites included ligand binding sites. These results suggest some selective pressure on the key residues, such as ligand binding activity, may cooperatively facilitate the emergence of a protein family during evolution. From a more practical aspect, the present results highlight an essential role of long-range effects in precisely defining protein families, which are absent in conventional sequence models.



  
  protein folding molecular evolution protein design sequence analysis Monte Carlo simulation
  














§ INTRODUCTION

Natural proteins can be classified into families based on their sequence similarity<cit.>. This is considered to be primarily a consequence of molecular evolution: proteins evolved from a common ancestral protein share similar sequences. However, evolution alone does not account for the existence of relatively well-defined (domain) families that are distributed rather discretely than continuously in the sequence space <cit.>. A key to understanding the family distribution is protein folding. As observed in protein structure classification databases <cit.>, each protein family corresponds to a unique three-dimensional fold, suggesting the existence of physical constraints imposed on protein sequences during the evolutionary process to maintain the fold <cit.>. While protein structures can tolerate great many mutations to the extent that proteins with little sequence similarity can share the same fold, residue conservation patterns reflect the structural context of protein sequences. This fact has long been exploited in protein structure prediction in the form of position-specific scoring matrices <cit.> and, more recently, direct-coupling analysis and related methods <cit.>.

Under a given physiological condition, a polypeptide is either able or unable to fold into some unique structure. This suggests the existence of a “folding transition” at the border between an “island” of a protein family and the “sea” of random polypeptide, that is analogous to the folding transition of a protein molecule in the conformational space <cit.>. It should be noted, however, that there are many families in the sequence space so that a sequence moving in the sequence space may fall into any one of these families. This is in contrast to protein folding in the conformational space where there is usually only one unique native structure for a given protein sequence. Furthermore, a (structural) domain, rather than a whole protein sequence, should be considered as a unit of folding as a particular domain may be found in different proteins in combination with other, different, domains.
Therefore, the system in which the analogy of protein folding holds should be limited to the vicinity of each protein domain family rather than the entire sequence space. In the following, we focus on the folding transition in a sequence subspace around a given protein domain family. Although biologically important, intrinsically disordered proteins <cit.> are excluded from the present study for the following two reasons. First, the analogy of the folding transition may not apply to those proteins. Second, it is difficult to obtain reliable and comprehensive multiple sequence alignments for this class of proteins <cit.>, which are required for parameter estimation of the statistical model employed in the present study. 

There have been a number of theoretical and computational studies on subjects related to the sequence space such as foldability and design <cit.>, molecular evolution within an island <cit.> and between islands <cit.>, or the size and/or distribution of islands in the sequence space <cit.>.
On the contrary, relatively little attention has been paid to the transition between an island (a set of sequences belonging to the same family) and the sea (the set of sequences that do not belong to the family) apart from a few exceptions. In the context of protein design, <cit.> theoretically predicted the existence of a “folding transition”.  In a study of hierarchical evolution of protein fold families based on a simple model of the evolutionary selection by native stability using a generic contact potential <cit.>, <cit.> observed a sharp transition at a certain design temperature. In neither of these studies, however, the nature of the transition was investigated further.
Characterizing the folding transition in the sequence space may help understand essential features that constitute a protein family and possible evolutionary trajectories that may have led to the emergence of a protein family. It also has practical importance in identifying new family members and designing new proteins. 

In the following, we investigate the folding transition in the sequence subspaces for 15 protein domain families including all-α, all-β, α/β and other folds by performing extensive Monte Carlo (MC) simulations of the modified lattice gas model (LGM) of protein sequence alignment <cit.> that coherently integrates long-range interactions and variable-length insertions. Using the LGM, the existence of a sharp two-state transition between natural-like sequences and random sequences in the sequence subspace is demonstrated. Furthermore, the nature of the transition is examined in detail by analyzing residue distribution of each site along the transition as well as by performing virtual “mutation” experiments.



§ THEORY

We briefly summarize the theory to the extent that is necessary for understanding the present study. For more details about the formulation of the LGM as well as the algorithms for MC simulations and parameter optimization, refer to the previous papers <cit.>.
The LGM for a given Pfam <cit.> family consists of N “core” sites and N-1 “insert” sites which respectively correspond to the “match” states and “insert” states of the Pfam profile hidden Markov model (HMM) of length N,
excluding the N- and C-terminal insert states. Exactly one of the 21 residue types (including the “delete” symbol) can exist on each core site whereas arbitrarily many (including zero) residues out of the standard 20 residue types can reside at each insert site.
The array of core and insert sites are connected via “bonds” (solid arrows in Fig. <ref>) that reflect the linear polypeptide structure. A pair of model sites connected via a bond are called a “bonded pair” in the following. Between core sites more than 2 residues apart along the sequence, there may be interactions (dashed lines in Fig. <ref>) based on a representative native structure of the family. Two sites are defined to be interacting if the residues aligned to those sites are in contact in the corresponding representative native structure. Two residues are defined to be in contact if any non-hydrogen atoms in those residues are within 5Å. Interactions are defined only between core sites for simplicity. Interacting core sites are referred to as “non-bonded” pairs in the following.


Let 𝐚 = a_1⋯ a_L be an amino acid sequence of L residues and an LGM ℳ of length N consist of core sites O_1,⋯, O_N and insert sites I_1 ,⋯, I_N-1. An alignment between the sequence 𝐚 and the model ℳ is represented as a sequence of pairs of a model site (core or insert) and a residue: 𝐗 = X_1⋯ X_L_𝐗 where
L_𝐗 is the length of the alignment and each X_i is a pair such as (S,a) with S∈{O_i}_i=1,⋯,N∪{I_i}_i=1,⋯,N-1 and a ∈{a_1,⋯, a_L}.
For example, given an amino acid sequence, say , and a model of length N=6 (Fig. <ref>), one of many possible alignments is represented as
𝐗 = X_1⋯ X_9 =
(O_1,𝙺)
(O_2,𝙲)
(O_3,𝙵)
(I_3,𝙿)
(I_3,𝙳)
(I_3,𝙶)
(O_4,-)
(O_5,𝚅)
(O_6,𝚆). Note there are multiple occurrence of the insert site I_3 whereas other insert sites are completely absent in this particular alignment.
Since there may be any number of residues at each insert site, the alignment length is variable.

Based on this representation of sequence alignment, the energy function of alignment 𝐗 is defined as

    E(𝐗)    =   
      -∑_k=1^L_𝐗-1J(X_k,X_k+1)
      -∑_(k,l) ∈𝒯K(X_k,X_l) 
            -∑_k=1^L_𝐗μ(X_k)

where J and K are short-range and long-range interaction parameters, respectively, μ's are chemical potentials, and 𝒯 indicates the set of all the interacting non-bonded pairs. The short-range interactions act only between bonded pairs of residues that are consecutive in the alignment (i.e., between X_k and X_k+1 in Eq. <ref>). The long-range interactions act between residues that are
aligned to interacting non-bonded pairs. The chemical potentials are so called because they are used to control the residue densities of each site. Only J and K parameters constitute intrinsic energy, and they are to be determined from a given (observed) multiple sequence alignment (MSA) of the family sequences (see below).

We assume that the probability P(𝐗) of obtaining an alignment 𝐗 is given by the Boltzmann distribution:

    P(𝐗) = exp[-E(𝐗)/T]/Ξ[T]

where T is the “design (selection) temperature”<cit.> in energy unit and Ξ[T] is the partition function 

    Ξ[T] = ∑_𝐗exp[-E(𝐗)/T].

Here, the summation is over all possible alignments (𝐗) of all possible amino acid sequences with the model. Since the alignment length can vary, this ensemble is considered to be a grand canonical ensemble.

Given the parameters J, K and μ, we can sample (alignments of) sequences according to the probability distribution Eq. (<ref>) by running grand canonical Monte Carlo (MC) simulations <cit.>.

We define the standard condition as the system with μ(X_k) = 0 for all X_k (i.e., no perturbation), and the natural condition as the standard condition with T = 1. 
The parameters J and K are determined iteratively so that the pair number densities of residues for bonded and non-bonded pairs over the samples produced by MC simulations under the natural condition match those observed in the given MSA <cit.>. This optimization corresponds to minimizing the following free energy function under the natural condition (T=1 in particular):

    F[T] = -max_J,K(⟨E|_⟩obs - Ω[T])

where ⟨E|_⟩obs is the average energy of the observed MSA and Ω[T] = -TlnΞ[T] is the grand potential. It can be easily shown that this is equivalent to the principle of maximum entropy <cit.>.

The design temperature in the LGM controls the mutation rate uniformly over the entire sites. The higher the temperature, the higher the mutation rate. The temperature of the natural condition T=1 defines the unit of energy and temperature. In other words, the temperature for each domain family is in a relative unit defined by the natural condition (unlike the absolute scale as in <cit.>). The mutation rate here means the number of accepted mutations per attempted mutation, and mutations are attempted at a constant time interval. The relationship between the present “energy” and physical folding free energy has been recently established by <cit.>.



§ RESULTS




 §.§ Sequence pattern formation is accompanied by a cooperative two-state transition

We have determined the energy parameters for each of the 15 families of relatively small domains listed in Table <ref> by using MC simulations following the procedure described previously <cit.>. 

Using the optimized parameters, we run a multicanonical MC simulation <cit.> for each family from which the specific heat,

    C(T) = ⟨E^2|_⟩T - ⟨E|^⟩2_T/T^2,

was computed as a function of temperature (Fig. <ref>A). For all the families but HATPase_c, a single peak in the specific heat was observed, indicating the existence of a transition. An increase in the specific heat of HATPase_c at higher temperatures was due to excessively long insertions. Depending on the family, the peak temperature (T_m; see Table <ref>) ranged between T=1.009 (for DnaB) and T=1.033 (for zf-C2H2).  The energy distribution obtained from grand canonical MC simulations at T_m (Fig. <ref>B) showed a clear two-state transition between low energy and high energy states.


In order to confirm if the low-energy sequences are indeed similar to the natural sequences, we compared the profile HMM of each family against the set of sequences generated by the MC simulations. Out of 10,000 sequences generated at each temperature ranging from T=0.95 to T=1.05, a profile HMM search was performed (using the hmmsearch program <cit.>) to identify the sequences that significantly matched the Pfam profile HMM with full sequence E-value < 0.01 (Fig. <ref>C). The hit ratio (the fraction of sequences significantly similar to the family's profile HMM) shows the clear two-state transition between natural-like sequences and non-natural-like (≈ random) sequences around the transition temperature T_m for each family, and the sequences generated at lower temperatures were indeed more similar to natural sequences. Since these natural-like sequences are expected to fold into the native fold of the domain family but non-natural-like sequences are not, this transition can be regarded as a “folding transition” in the sequence subspace.
The hit ratios at T=1 (the natural condition) ranged from 73.3% (for zf-C2H2) and 92.5% (for HATPase_c), indicating there were a non-negligible fraction of non-natural-like sequences at the natural condition. In particular, the C2H2-type zinc finger domain (zf-C2H2) consisting of only 23 core sites exhibited a large fraction of non-natural-like sequences at lower temperature.
For all the families examined, the hit ratio was a decreasing function of temperature. For later reference, we defined a “disordering temperature”, T_D, of each family as the lowest temperature at which the hit ratio was less than 1% (column “T_D” in Table <ref>).



 §.§ Comparison with a model without long-range interactions: An example


The long-range interations incorporated in the LGM are necessary (though not sufficient) for the model to exhibit the cooperative transition.
Conventional sequence models, however, do not include such long-range interactions. In other words, they are essentially one-dimensional models. It is a well-known fact that one-dimensional systems do not exhibit cooperative (phase) transitions at finite temperatures <cit.>. Although theoretically trivial, it may be instructive to demonstrate the absence of the transition in such a system. Within the framework of the LGM, we can simply discard the long-range interations and train a model without the K parameters in the energy function (Eq. <ref>). In this “J-only” model, the J parameters for short-range interactions can be obtained exactly <cit.>. We trained a J-only model for the SH3 domain family (Table <ref>) and compared its specific heat, hit ratio and energy distributions with those of the original “J+K” model which includes long-range interactions (Fig. <ref>).

The specific heat of the J-only model monotonically increased with temperature (Fig. <ref>A) which is in sharp contrast with the J+K model. The upper limit of the temperature range for the J-only model was set to T=1.17 because it was only at this temperature the hit ratio became less than 1% (Fig. <ref>B), which is significantly greater than the corresponding disordering temperature for the J+K model (T_D = 1.05; Table <ref>). Accordingly, the hit ratio decreases rather slowly as the temperature increases. The energy distributions of the J-only model at varying temperatures are unimodal and the mode energy value continuously increases with increasing temperature (Fig. <ref>C). This indicates that there is no clear boundary between natural-like sequences and random sequences. 
These results for the J-only model clearly demonstrate the absence of cooperative transition in the model without long-range interactions. It it expected that other conventional sequence models also behave similarly.

In the following, we exclusively use the J+K models.


 §.§ Less conserved sites order at higher temperatures, better conserved sites order at lower temperatures

Next, the transition was characterized in terms of residue distribution at each site over a temperature range.
The “naturality” of each core site of simulated sequences at a given temperature was evaluated by the Kullback-Leibler divergence <cit.> D_O_i(T):

    D_O_i(T) = ∑_a=1^21⟨n_(O_i,a)|_⟩Tlog⟨n_(O_i,a)|_⟩T/⟨n_(O_i,a)|_⟩obs

which measures the difference between the residue distribution ⟨n_(O_i,a)|_⟩T of the site O_i at temperature T from the observed residue distribution ⟨n_(O_i,a)|_⟩obs. An example for the SH3 domain is provided in Fig. <ref>C. Lower divergence values imply more natural-likeliness.
By construction, the divergence is nearly 0 for all the sites at T = 1, but the values varied largely among different sites for higher temperatures. In general, better conserved sites (c.f., Fig. <ref>A) have larger divergence at high temperatures. In the following, we use the D_O_i(T_D), i.e., divergence at the disordering temperature T=T_D, as a measure of conservation. As expected, sites with a larger contact number <cit.> (i.e., the number of non-bonded contacts a site makes with other sites; see Fig. <ref>B for example) tend to be more conserved as can be seen in their correlations (Fig. <ref>A).




Ordering of each site can be quantified via the derivative of the divergence with respect to temperature (c.f., Fig. <ref>D):

    C_O_i(T)    =   dD_O_i(T)/dT
       =   1/T^2∑_a=1^21[⟨n_(O_i,a)E|_⟩T - ⟨n_(O_i,a)|_⟩T⟨E|_⟩T]
          ×log[⟨n_(O_i,a)|_⟩T/⟨n_(O_i,a)|_⟩obs]

which we call “differential divergence” in the following.
Note that this is a weighted sum of covariance between the number density and total energy where the unnormalized and signed weights are log[⟨n_(O_i,a)|_⟩T/⟨n_(O_i,a)|_⟩obs]. Thus, the peaks of this quantity provide the site-specific transition point at which a large change in the residue distribution is accompanied with a large change in total energy.
Precisely, the site-wise transition temperature T_m,O_i is defined by

    T_m,O_i = argmax_T < T_D C_O_i(T)

(see Fig. <ref>E for an example).
For some sites, this temperature cannot be defined because of the absence of a unique peak for T < T_D (due to lack of non-bonded interactions; this was observed for several sites of the fn3 family).
The range of the site-wise transition temperatures are plotted along with T_m and T_D in Fig. <ref>B. The site-wise temperatures all appeared slightly above the global transition temperature T_m, indicating a global (higher order) ordering process between the minimum T_m,O_i and T_m.
The site-wise transition temperature was found to be negatively correlated with site conservation D_O_i(T_D) (Fig. <ref>C) although the correlation was not always significant. Thus, as temperature lowers, local ordering starts from less conserved sites, but these sites have little global effects (smaller values of C_O_i(T_m,O_i); see Fig. <ref>E, green line), then it proceeds to increasingly better conserved sites with increasingly greater effects (larger change in total energy). Then, higher-order ordering proceeds down to the global transition temperature, T_m. 
Interestingly, a nearly perfect correlation (R ≈ 1) was found between {C_O_i(T_m)} and {D_O_i(T_D)}. Thus, the transition of each site appears almost completely determined by the conservation patterns.



 §.§ Suppressing key residues inhibits sequence pattern formation under the natural condition: “Alanine” scanning


Although the above analysis of divergence maps hinted some differential roles of individual sites
in determining natural-like sequences, the site-wise transition temperature T_m,O_i did not vary
much across the sites. 
In order to explore the role of each site more directly, we performed virtual site-directed mutation experiments by manipulating the chemical potential μ.
First, we performed “alanine” scanning experiments. That is, for each core site O_i, the chemical potential for alanine (“𝙰”) residue (or glycine “𝙶” if the consensus residue was alanine)  was set to a large positive value (namely, μ(O_i,𝙰) = +20) while other μ's were set to 0. This effectively fixes alanine residue at the site O_i while allowing other sites to freely adopt residues that are consistent with the perturbation. By performing MC simulations at T=1 and applying the hmmsearch program to the sampled sequences, the hit ratio h_O_i^ALA (the fraction of natural-like sequences) of the “alanine mutant” was determined (c.f., Fig. <ref>C, magenta line). This hit ratio was compared to the hit ratio h_0(T=1) of the natural condition (c.f., Fig. <ref>C). Namely, we evaluated the “mutants” in terms of the following enrichment value:

    En_O_i^ALA = h_O_i^ALA / h_0(T=1).


While most of the sites were not largely affected by alanine mutations (the enrichment ranging between 0.5 and 1.5), a few sites exhibited enrichment of less than 0.5 (Fig. <ref>A), that is, the fraction of natural-like sequences for these few mutants was less than a half of that under the natural condition.
The hit ratio of mutants {h_O_i^ALA}_i=1,⋯,N and a measure of conservation of each site {D_O_i(T_D)}_i=1,⋯,N showed negative
correlations (Fig. <ref>C; R < -0.50 for most families), indicating better conserved sites are affected more severely by the alanine mutation. This result suggests that fixation of some key residues is essential for the global pattern formation of natural-like sequences.




 §.§ Enforcing key residues promotes sequence pattern formation under disordering condition: “Enforcedly conserved residue” scanning


Conversely to alanine scanning, we can enforce the consensus (the most dominant) residue a_O_i^* of each core site O_i by setting μ(O_i,a_O_i^*) = +20 (and all other μ's were set to 0) at the disordering temperature T = T_D. With this “enforcedly conserved residue (ECR)” mutation at each site, the hit ratio h_O_i^ECR was calculated (c.f., Fig. <ref>C, green line) and the enrichment with respect to the standard hit ratio at T_D

    En_O_i^ECR = h_O_i^ECR / h_0(T=T_D).

was compared (Fig. <ref>B).
While majority of sites exhibited enrichment of less than 10 fold, a few sites showed more than 40-fold enrichment.
The hit ratio for ECR mutants are even more correlated with site conservation (Fig. <ref>C; R>0.87 for all the families).
Thus, enforcing a few well-conserved residues can enhance the sequence pattern formation over the entire sequence under otherwise disordering conditions.





§ DISCUSSION

The LGM for studying the “folding transition” in the sequence subspace around a given protein family is analogous to Gō-like models for studying the folding transition in the conformational space for a given native structure <cit.>. One major conceptual difference is that while there is a unique native structure as the global energy minimum for a Gō-like model, there is an ensemble of sequences as the global free energy minimum for an LGM (Eq. <ref>). This is because in the latter the parameters are optimized using a set of natural sequences at a finite temperature (T = 1). This may be considered as a consequence of the asymmetric one-to-many relationship between a protein fold and family sequences.
The design of the LGM assumes a funnel-like free energy landscape of the sequence subspace analogous to that in the conformational space <cit.>. Based on simple exact models, it has been suggested that the free energy landscape of the sequence subspace within a family of sequences sharing the same ground-state conformation is indeed funnel-like (“superfunnel”) <cit.>. The present LGM extends the superfunnel concept to include sequences slightly beyond the family boundary.

The parameters of the LGM are determined so that statistics of simulated sequences match that of an observed MSA. This may imply that the sequences sampled by the LGM are biased towards natural sequences. Nevertheless, based on computational protein design experiments, it has been suggested that the volume of sequence space optimal for a protein structure is restricted to a region around the natural sequence<cit.>. Therefore, the ensemble of sequences sampled by the LGM under the natural condition is expected to be a good approximation for the ensemble of all the sequences, both natural and artificial, compatible with the protein family.

The existence of a rather sharp two-state transition implies the existence of a clear boundary between family members and non-family members. This transition is made possible by the inclusion of the long-range (non-bonded) interactions in the LGM. Conventional sequence models such as profile HMMs do not (or cannot) incorporate such long-range effects because they are essentially one-dimensional (1D) systems, and therefore, they cannot exhibit such transition <cit.>. Thus, the discrimination between family and non-family members based on the conventional sequence models is necessarily fuzzy (Fig. <ref>C), inevitably leading to false positives and false negatives. In fact, this situation may have already biased our understanding of protein families. For example, the C-terminal 3_10 helix and β strand (β5 in <cit.>) of the SH3 domain are completely missing in the corresponding Pfam profile HMM (c.f., grey regions in Fig. <ref>), but these regions not only comprise an integral part of the SH3 fold, but also contain a functionally important tyrosine residue (Fig. <ref>) involved in peptide binding <cit.>. More generally, there are many partial domains found in the MSAs of Pfam and other databases that are shorter than 50% of the domain model length, and most of them are considered to be alignment and annotation artifacts <cit.>. The present results suggest that these alignment artifacts are an inherent property of conventional sequence models.
A notable exception to the conventional sequence models is MRFalign <cit.> based on essentially the same principle as the LGM <cit.>. Based on the present discussion, the success of MRFalign in remote homology detection is explained by its ability to set a clear boundary for family members owing to the inclusion of long-range correlations.


In this study, we have used one of the conventional methods, namely the profile HMM, for quantifying natural-likeliness in terms of the “hit ratio”. However, due to the limitation of profile HMMs just pointed out, such a measure may be inaccurate. In designing artificial WW domain proteins, Socolich et al. <cit.> showed that sequences designed without pairwise correlations were unable to fold whereas all the designed foldable proteins were designed with pairwise correlations (by statistical coupling analysis <cit.>). It has been shown that other statistical models incorporating “direct” long-range correlations can discriminate foldable sequences from unfoldable ones to a very good accuracy <cit.> for the artificially designed WW domains <cit.>. The LGM also incorporates the direct long-range correlations in the form of interacting non-bonded pairs so its energy value is expected to be a good discriminator of foldable proteins.
It is of interest and in principle possible to experimentally measure the foldability of the sequences generated by the LGM.

In the context of molecular evolution, the present results suggest that if some selective pressure is imposed on one or a few key sites, the entire sequence may evolve rather quickly to form a pattern compatible with the family characterized by a unique fold. 
How can such selective pressure act on residues that are to be conserved?
According to recent studies, it appears that protein families have been naturally (and positively) selected so as to stabilize their native folds <cit.>. In fact, most of the key residues listed in Table <ref> comprise the structural core of native folds, suggesting the structural importance of these residues. However, a specific native structure per se is unlikely to be subject to natural selection unless it accompanies some advantageous phenotype <cit.> through molecular functions. Visual inspection of the native structures in Table <ref> and their homologs revealed that the key residues often included those involved in intermolecular interactions such as binding of small molecules, polypeptides/proteins or polynucleotides, or phosphorylation (Table <ref>, residues marked in boldface). In the case of the SH3 domain (Fig. <ref>), for example, the key residues Y6, W32 and P47 (consensus residues according to Fig. <ref>A) comprise the “xP pockets” essential for the proline-rich peptide binding activity <cit.>. Thus, although it has been suggested that functionally important residues are not necessarily optimal for structural stability <cit.>, some of them play an essential role in cooperatively shaping the overall pattern of family sequences and hence the native fold <cit.>. This may also help explain the observation that most of the structural motifs for ligand binding are confined within single protein families <cit.>. 



§ MATERIALS AND METHODS



 §.§ Data preparation

All the 15 domain families studied here were selected from the “top 20” list of the largest Pfam families (except for the SH3 domain) that had average sequence length of less than 120 residues. Only those of the “domain” type were used (i.e., “family” and “repeat” were discarded). These domain families consist of various folds.
Multiple sequence alignments (MSA) and profile hidden Markov models (HMMs) were downloaded from the Pfam database (version 30.0). The MSAs were based on the “uniprot” set when available or on the largest available representative set, otherwise (Table <ref>). Aligned sequences that contained more than or equal to 10% deleted residues on core sites were discarded. A representative structure of each family was assigned based on resolution and alignment quality using the PDBj Mine2 relational database <cit.> integrated with the SIFTS <cit.> resource (Table <ref>).



 §.§ Monte Carlo simulations

In each MC simulation of an LGM model of length N, one sweep consists of N and N-1 attempted mutations for randomly selected core and insert sites, respectively. For the algorithmic details of the MC simulations of the LGM, refer to <cit.>.
For the results given in Figs. <ref>B, C, <ref> and <ref>, 20 grand canonical MC simulations (under the conditions described in the text) were performed for 10^6 sweeps after 5× 10^4 sweeps of equilibration, and 10,000 sequences (alignments) were saved every 100 steps for each of the 20 trajectories. For each set of 10,000 sequences resulted from each trajectory, a sequence database was created by removing the “delete” symbols from the sequences, against which the hmmsearch <cit.> program was applied with the Pfam profile HMM with the E-value cutoff of 0.01 in order to identify significantly similar sequences. From the 20 runs, the average hit ratio and the standard deviation were obtained. 

For the results given in Figs. <ref>A, <ref>, and <ref>, first, the density of states was determined by the Wang-Landau method <cit.> by using the energy bins with bin width of 0.2 energy units. The lower bound of the energy range was fixed at a sufficiently low value, and the upper bound was determined by tentative grand canonical MC simulations at a high temperature T = 1.2. The Wang-Landau iteration was terminated when the density of states was determined to the precision better than 10^-7.
Next, using thus determined density of states, a multicanonical MC simulation <cit.> was performed until all the energy bins are filled with at least 2,000 counts (sequences were saved at every sweep). From this trajectory, the specific heat and residue distributions at various temperatures were obtained by the reweighting method <cit.>.



§ ACKNOWLEDGEMENTS

The author thanks Ken Nishikawa, Motonori Ota and Gautam Basu for valuable comments, Eugene I. Shakhnovich for providing some references.














§ REFERENCES

elsarticle-harv 




