Contributions of topological polar-polar contacts to achieve better folding stability of 2D/3D HP lattice proteins: An in silico approach

1.   Introduction

Nature teaches that there is a clever selection of copolymer sequences of biological molecular monomers aimed to overpass physicochemical challenges. These are well-known as foldamers in the biological sciences and could be generated from aminoacids, nucleic acids or sugars [1]. Such physicochemical challenges are those needed for achieving key functions as: structural shaping & templates, chaperones/transporters, molecular crowders & recognition motifs, phase transfer agents or surfactants, antioxidants, ionic/molecular reservoirs, detoxifying agents, catalysts or inhibitors, among other biophysicochemical functionalities. These tasks are molecularly fulfilled because folded structures were kinetic/thermodynamic balanced where particular structures corresponds to specific functions [2]. This is an intense research field in the proteins domain, being all-atom approaches the most complex case, and perhaps the hydrophobic-polar (HP) Dill's model represents the simplest case, that is why it's entitled “coarse grained model”. All these deal with a classic biological problem called protein folding paradigm stating that: the characteristic or native three-dimensional structure of a given proteic sequence is completely determined (encoded) by the aminoacid sequence itself [3]–[5], very related to the “protein folding backbone approach” proposed by Rose et al [6]. Nevertheless, both the all atom and the coarse grained models [7], have such high complexity that represent NP-complete systems from the mathematical viewpoint, thus they are clustered into “computationally intractable” set of problems [8].

So no matter of which level of “coding” is the protein folding problem treated its complexity gave insights to many approaches to treat it, in parts, or as a whole, experimental of theoretically. In this line, there is the HP Dill's model, being a simple two letter coding, in which H is hydrophobic and P is polar, both aminoacid moieties, embedded in a more complex proteic sequence oversimplifying the natural 20-aminoacid coding [4],[5],[9],[10]. More sophisticated models are three letter code (3LC), or four letter code (4LC), e.g., the 4LC named HPNX-model is constructed from the split of the polar charged monomers as positive (P), negative (N), and neutral (X) [11],[12]. These diversity of models, ranged from simple lattice models [13]–[17], followed by a myriad of intermediate models, to finally reach the all atom classical molecular dynamics simulations, or even its quantum counterpart [18]–[22]. All these computational approaches have become cornerstone trials intended to understand such complexity, where authors have stated: computing has been able to provide ab initio abstract/conceptual hints or motifs leading to other derived theoretical or experimental developments. This latter aids even in cases where none evidence, nor experimental nor theoretical, is available [20]. But returning to basics the HP-model, no matter of its simplicity, cleverly states that hydrophobic interactions among aminoacids [23] represent one of the principal driving forces that yield native states in proteins. All the latter occurring in such a fast folding as required by natural or anthropogenic molecular processes. Hence, a lot of information and research is yet to be surveyed and discovered, even in the simplistic HP-folding approach, as will be seen forthcoming.

It is well known that the non-covalent interactions, both polar and nonpolar, play an important role during protein folding, mainly those formed by hydrophobic contacts, but as well as hydrogen bonds, aromatic interactions, and salt bridges, for example [24]. Therefore, only considering a hydrophobic core, there is not considered enough structural information nor the total interactions of the real system. It has been reported in the literature that hydrophobic contacts contribute 60 ± 4% and hydrogen bonds 40 ± 4% to protein stability [25],[26], hence both interactions are important factors that contribute in stabilizing the native folded state of proteins. Therefore, the polar counterpart of the protein system scored by polar-polar (P···P) interactions, should be taking into account.

Because of the nature of the side chains of aminoacids, physicochemical differences among them are present, developing important changes in polarity, size and conformation. These properties are responsible to tune the packing of aminoacids comprised in proteins. Nevertheless, itself hydrophobicity is not enough to achieve the accurate folding. Therefore, its polar counterpart accounts for the rest of molecular interactions between aminoacids, also molding the outer protein structure. Where, polar residues could strongly interact between them by P···P contacts developing mainly electrostatic interactions, e.g. of short range as salt bridges, zwitterionic contacts, etc. [27]. Indeed, salt bridges and zwitterionic contacts have resulted very important to enhance protein folding stability and they are generally formed in their outer/polar surface [28],[29]. Moreover, polar aminoacids tend to accommodate towards this polar media attracting water molecules and forming hydration cores surrounding polar protein surface [28],[29]. Some of these latter are the main contributions to justify taking into account the P···P interactions trying to enhance the simple HP model. Some authors argue that the inclusion of P···P interactions also favor structure compactness, as happens in real proteic systems [30]. Particularly, Kumar & Nussinov clearly state “While the hydrophobic effect is the major driving force in protein folding, electrostatic interactions are important in protein folding, stability, flexibility, and function” [29]. One major enhancement is that these augmented HP models are known to develop protein-like features undoubtedly indicating that lattice models include in these simple metaheuristics the fundamental physicochemical principles of proteins [30].

With all the latter in mind, in this current research the main goal is to test another very simple metaheuristic to include P···P interactions into the HP lattice model. In our work group we have explored the study and design of simple HP protein structures considering: 1) the hydrophobic-hydrophobic (H···H) interactions and 2) the restrictions imposed by the intracellular space as an osmolyte effect or molecular crowding [31]–[33]. In this work we are now proposing to survey the effect of also including P···P interactions to develop better folding strategies. In this new proposal of augmented HP model, the formation of the hydrophobic core is permanently prioritized, and later on the formation of the polar substructures, in contrast with other augmented HP models where this important requisite is not always maintained [30]. This by considering more contributions and physicochemical properties present in these simplistic HP protein systems. For this purpose, were tested 11 HP sequences, some that were difficult to fold in previous works and some other newly designed, in two types of lattices, 2D-square and 3D-cubic. These sequences were selected/designed taking into account including H···H as well as P···P contacts tracking the effect of these two contributions in the sequence folding procedure.

2.   Materials and methods

2.1.   A lattice model including H···H and P···P interactions

As can be seen below in expressions (1) to (2) describing Model 1 (Dill's model) [34], the HP-model is abstracting the protein interactions by labeling the aminoacids as H or P. Where backbone structures will be energetically chosen by optimizing the amount of -H···H- contacts due to a hydrophobic effect. Nevertheless, substructures based on -P···P- interactions, such as hydrogen bonds and salt bridges should also be considered for structure prediction. Hence, extending this traditional HP model to a variant where P···P interactions are also considered in the optimization process conducted to the design of a new approach named Model 2, were we use a variant of the convex function, as shown in the expressions (3) to (6) below. Here, this variant of convex function has been employed to tune the weight between H···H and P···P interactions. That is, each H···H interaction is assigned the value α–1 while each P···P interaction is assigned the value -α, with 0 ≤ α ≤ 1.

Model 1: Only the H···H interactions are optimized (Dill's model):

where e_ij is defined as:

Model 2: Inspired by the convex function, where both H···H and P···P interactions are optimized:

where:

i. h_ij counts the hydrophobic interactions defined as:

ii. p_ij counts the polar interactions defined as:

iii. α, 0 ≤α ≤1, determines the weight given to each type of interaction. Alpha values very close to 0 favor the hydrophobic collapse observed in Dill's model (expression (1)).

iv. Expression (3) is subject to the following restriction: the best hydrophobic core obtained (optimum number of H···H interactions) must be conserved.

v. The term score is defined here as a double conditional in function F (expression (3)) that i) like Dill's model, generates the best score in the HH core, as the maximum amount of H···H interactions and ii) maximized the outer P···P contacts, without breaking/disassembling the HH core.

Note that in expression (3), when α = 0, Model 2 is simplified to Model 1. And the score term will be used to rank folding results.

One of the principal aims in this contribution is testing α values in Model 2 that conduct to the formation of structures with better local minima, preserving hydrophobic core and now folding also PP substructures.

2.2.   The assembly and optimization of 2D-HP structures & methodology

HP-model states that any aminoacid sequence (S_i) could be transcribed and defined as expression (7):

where “n” is the length of the chain.

As shown in Figure S1 (in Supplementary), the HP transcribed sequence could then be folded in a simplistic 2D or 3D-lattice, following discrete movements between the neighboring cells that conform the lattice. Particular details of this are given in Supplementary Material.

From the transcribed HP-sequence, a random population of HP 2D/3D-structures is generated in a 2D/3D-lattice, following the algorithm explained and illustrated in Figure S2 (in Supplementary). Subsequently, on this random population of structures, an optimization process begins executed by an evolutionary algorithm [31],[32] whose fitness function is given precisely by Model 2. The main characteristics of the evolutionary algorithm used are listed in Table S1 (in Supplementary). Moreover, the folding methodology was carried out for the in silico experiments on our Evolution bioinformatics platform, see further details in Supplementary.

3.   Results

3.1.   Structural design: from HP-sequences to HP simplistic structures

In a previous work [33], we designed a set of simple HP-sequences and the corresponding optimal 2D-HP folding for each of these. The latter allowed us to know in advance the optimal score expressed only by the number of H···H interactions, that characterized each of the expected structures. On that research, the optimization of these 2D-HP structures was carried out employing an evolutionary algorithm whose fitness function was given by Model 1, only developing HH structures, as described above. In the present work, we return to some of these HP-sequences designed to fold in 2D-square lattice model, and plus we incorporated new HP-sequences as well as the corresponding expected 2D or 3D-structures. And moreover, we use Model 2, now including also P···P interactions to pursue a better folding score, but also to verify the contributions of H···H contacts in its own HH substructure, but now also exhibiting P···P contacts in order to form the missing PP substructure. Table 1 lists the characteristics of the HP-sequences in which we will test the contribution of P···P contacts in the formation of the corresponding expected structure. The 2D-sequences are S₁-S₈, meanwhile the 3D-structures are S₉-S₁₁.

3.2.   Role of the computational tool and the methodological approach in obtaining the results

The folding of the sequences proposed in Table 1 and, therefore, the generation of the expected 2D/3D-structures, was carried out as described in Supplementary (employing parameters stated in Table S2) using the Evolution computational tool (http://bioinformatics.cua.uam.mx/site/) [35],[36] and guided by the methodological approach described in Supplementary (see Figures S2 & S3). Evolution is a bioinformatics platform developed by our work group [31]–[33], where its functionality has recently been enhanced from Model 1, to now implement optimization based on Model 2.

It is necessary to mention that the values of α (in Model 2) are the result of a preprocessing phase carried out by executing many batches of experiments a priori, finding that α values greater than 0.3 led rupture of hydrophobic core in the majority of structures, thus banning the hydrophobic collapse. When α is higher it moves away from the Dill's model, hence taking a bet for this concept of hydrophobic collapse as main driving force, and that it should be maintained but not as a unique contribution. The Model 2 optimization normalizes between 0 and 1 the importance of ponderation or weighting given to both H···H and P···P interactions in the folding process. So that the sum of these weights always equals to 1. Moreover, α values equal or lower than 0.3 means that the formation of the hydrophobic core is only receiving 70% or less of the overall weight, in the expected 2D/3D-structures. And this value is just enough weighting to ensure that hydrophobic core is preserved, reason why this weighting scheme still provides the required hydrophobic collapse in structures.

3.3.   Analysis of the in silico experiments

Also, in Table 1 are gathered the results of S₁-S₁₁ folding, applying both Model 1, where α = 0, and Model 2, where α ≠ 0. As an example, in Figure 1A is observed that topological H···H contacts are 9 for S₁, therefore applying Model 1, the optimal fitness (OF) is −9 and the value of best fitness found (BFF) in experiments resulted also −9, but all P···P contacts are not optimized (Figure 1B). When considering Model 2, we tested different α values and the best value of α (α_B), in this case 0.3, is the one that produced the expected optimal structure also folding PP substructure (Figure 1C). Therefore, applying Model 2, the HH fitness (F_HH) was obtained by multiplying the total expected hydrophobic contacts (f (h)) by (α_B-1), that is, f (h) (α_B-1) = 9 (−0.7) = −6.3. And the PP fitness (F_PP) was obtained by multiplying the total expected polar contacts (f (p)) by (−α_B), this is f (p) (−α_B) = 21(−0.3) = −6.3. Finally, the BFF is equal to the sum of F_HH and F_PP, that is (−6.3) + (−6.3) = −12.6, precisely matching this value with that of OF. And in this same way all the other experiments were carried out and BFF & OF values were obtained.

4.   Discussion

4.1.   Structural analysis of folded S₁-S₁₁

Figure 1 shows the behavior of sequence S₁, taking account Model 1 (Figure 1B) and its respectively modification using Model 2 (Figure 1C). S₁ is a H₁₆ hydrophobic core, nevertheless, the terminal P₁₁ & P₂₀ branches were stipulated to be one larger and one shorter, in order to survey their correct fold and its own PP substructure. This should cause a fine tune in the proposed α value to preserve the hydrophobic core but also to fold a new polar core in a neighboring fashion. Another possible effect with this is how the H···H core degeneracy behaves with this new weighting scheme. It was easy to observe that Model 1 results only developed well folded hydrophobic core, obtaining the expected nine H···H contacts, but few P···P contacts were produced. Model 2 usage with α value of 0.3, the hydrophobic core is maintained and a compact PP substructure with 21 topological contacts was developed, hence this value of α = 0.3 was considered to be optimum (Figure 1C).

The 2D-HP structure expected as a result of the optimal folding of the S₆ sequence is illustrated in Figure 2A. There, the expected structure corresponds to a 5 × 5 hydrophobic core with short P₄ branches, positioned between each two rows of HH core. This 2D-HP structure was proposed by us [33], where its folding was studied using both Model 1 and a simulated intracellular medium.

The resulting folding of S₆ sequence in 2D lattice is illustrated in Figure 2B-D. As can be seen in Figure 2B, when sequence folding is guided by Model 1, optimal conformation of hydrophobic core is reached (OF_H···H = −20), however, only 2 of the 4 expected polar substructures achieve proper conformation. Note that the suitable conformation of a polar substructure is characterized by a single polar contact, which corresponds to a P···P score equal to −1. Adequate conformation observed of these two polar substructures is due to restrictions imposed by 5 × 5 hydrophobic core formation, there is no other way to achieve this. Figure 2C shows optimal conformation of hydrophobic core with OF_H···H = −20, folded using α = 0.2, also achieving conformation of the 4P substructures, although one of them does not show expected orientation. One of the goals of this experiment was to achieve OF_P···P = −8 (see Figure 2C and Table 1), but only OF_P···P = −6 was obtained at least in the batches of experiments carried out. But another goal was also to achieve the conformation of the 4P substructures and reach local OF_P···P = −4, i.e., -1 for each one of the 4P₄ substructures, which was well accomplished. It should be pointed that for α = 0.05 and α = 0.1 it was also possible to reach the best possible conformation shown in Figure 2C. However, for α > 0.2, the expected 2D-HP optimal structure was not achieved (Figure 2D).

The S₈ sequence should produce the complex 2D-HP structure illustrated in Figure 3A as a result of its folding. There, the expected secondary structure corresponds to a 5 × 5 hydrophobic core included in a 7 × 7 polar-nonpolar structure. Figure 3B illustrates S₈ folding guided by Model 1, showing an imperfect hydrophobic core, since OF_H···H = −24 and BFF_H···H = −22 (Table 1). This result shows both complexity of expected 2D-HP structure and the need of restriction/boundary conditions provided by considering P···P contacts contribution. Figure 3C shows optimal conformation of HH square 5×5 with α = 0.05, achieving OF_H···H = −24, at the same time that surrounding polar perimeter reached expected OF_P···P = −8. However, acquired optimal HP-structure was not exactly the expected in the surrounding polar substructure (see Figure 3A), being a degenerate optimal. For α > 0.05 not optimal nor degenerated-optimal structures were found, see one example in Figure 3D.

For the 2D folded sequences (S₁-S₈), S₈ sequence was the one that produced the most complex structure. Hence, when Model 2 was used, folding of S₈ sequence was precisely the one that required a lower α value to yield optimized structure. Seeming that at higher complexity of the hydrophobic core of the optimal structure, the required α value should be lower. Also, with α values higher than optimum, expected HH core is not achieved but misfolded structures occurred.

The folding of the 2D HP S₂, S₃, S₄, S₅, and S₇ sequences, using Model 1 and Model 2, is provided in the Supplementary, Figures S4-S8, respectively.

The S₁₁ sequence (see Table 1) is a very interesting case to fold in 3D-space. Looking at its primary structure (P₄H₁₆P₄H₁₆P₄) we could think that optimal should be a 3D-structure made up of two 4 × 4 hydrophobic layers, with 16H residues each; and 3PP substructures, each with P₄ residues, two located at the end of each 4 × 4 hydrophobic layer and one as a connecting unit between two 4 × 4 hydrophobic layers. However, finding the optimal folding does not only mean reaching the maximum number of H···H contacts, with Model 1, or H···H/P···P contacts, with Model 2. This process also involves maximizing the compactness of resulting 3D-structures, as we will see below.

The folding of S₁₁ is shown Figure 4A using Model 1, and in Figure 4B-D using Model 2. In Figure 4A, optimal folding of hydrophobic core is achieved, with OF_H···H = −34 (see Table 1). Note that P···P contacts were not obtained. Also optimal folding of hydrophobic core was obtained by compacting structure into 3HH and not 2HH layers, resulting in the most compact structure found by optimization algorithm. As shown in Figure 4B, with α = 0.1, the optimal conformation of the hydrophobic core is preserved while the polar substructures also reach their best folding, achieving expected 8 P···P contacts (Table 1). Note that each volume shares two rows between them as a boundary condition. With values higher than optimal α = 0.1, conformation of hydrophobic core was not achieved nor PP substructure, e.g. Figure 4C-D show suboptimal structures with α = 0.2 & 0.3, respectively.

The folding of the 3D HP S₉ and S₁₀ sequences, using Model 1 and Model 2, is provided in the Supplementary, Figures S9 and S10, respectively. On the other hand, degenerated structures are diminished using Model 2, due to conformational restraints imposed by PP substructure formation. Tested structures in 2D & 3D correctly folded with α = 0.05–0.3 and α = 0.1, respectively. Values of α beyond these yielded misfolded structures.

4.2.   Summary discussion

The summarized results of the folding of S₁ to S₁₁ sequences, using Model 1 and Model 2, are shown in Table 2 and Figure 5. Here, it is necessary to mention that for each sequence, 10 batches of experiments were run per approach, each batch consisting of 10 trials. Table 2 provides for each sequence the following features: 1) best size of hydrophobic core, 2) best number of polar contacts, and 3) the best total number of contacts), using both Model 1 and Model 2. The folding tendency for features 1), 2) and 3) is shown in Figure 5 captions (A), (B), and (C), respectively.

As can be seen in Figure 5A, the expected hydrophobic core was reached for all eleven sequences when Model 2 was used. Only one miss was found for Model 1, failing to reach the expected hydrophobic core of the S₈ sequence. This fact could be justified as a consequence of the role played by P···P contributions correcting in such a way the folding of the hydrophobic core. For instance, when analyzing some of the expected HP 2D structures (see Figures 1A, 2A & 3A), it could be seen that the expected 2D P substructures impose restrictions on the degrees of freedom of the hydrophobic elements, and as a result contributing to the formation of the expected hydrophobic core if the sequence allows the latter. Figure 5B shows the P···P folding tendency, therein it could be noted that when Model 2 was used the expected polar substructure was achieved for ten of the eleven sequences, where the most contrasting examples for this polar effect are for S₁, S₉ & S₁₁. Only in the case of sequence S₆, Model 2 failed to achieve the polar optimum substructure (Figure 2C). Lastly, in Figure 5C can be observed the global folding tendency, that is the joint contribution of H···H and P···P interactions, for sequences S₁ to S₁₁ using Model 1 and Model 2 clearly showing that Model 2 is in all cases better than Model 1, again where the most contrasting examples are for S₁, S₉ & S₁₁.

5.   Conclusions

Here were tested 11 HP-sequences, S₁-S₈ in 2D-square lattice and S₉-S₁₁ in 3D-cubic lattice using two folding approaches, a) Dill's model, named Model 1, and b) a model inspired by convex function, named Model 2. Last model is heuristically aimed to weight H···H (Dill's model) and also P···P contacts, to gather more structural information in order to reach better folding solutions in any given HP-sequence. In Model 2, H···H interactions were tuned as α-1 and P···P interactions as -α, and HP folding in all cases was more successful than Model 1. When α values above 0.3 are employed to fold (high P···P weighting) this started to ban H···H contacts, and misfolding occurred. There were needed α values very close to 0.3 to optimize longer or shorter sequences, with low difficulty in their folded structures. In comparison, more complex 2D-sequences, required much lower α values of 0.05-0.1 to achieve the accurate folding, being the length of the sequence not as important as its structural complexity. The 3D-sequences, of higher dimensionality and complexity, also required α = 0.1 to improve folding results. Moreover, α parameter might be included into a multi-objective optimization scheme, in such a way that the same algorithm should be able to find both the best α value and the optimal folding at once.