Irena Roterman1*, Mateusz Slupina1, Leszek Konieczny2
1Department of Bioinformatics and Telemedicine, Jagiellonian University Medical College, ul. Medyczna 7, 30-688 Kraków, Poland
2Chair of Medical Biochemistry, Jagiellonian University Medical College, ul. Kopernika 7, 31-034 Kraków, Poland
*Corresponding authors: Irena Roterman, Department of Bioinformatics and Telemedicine, Jagiellonian University Medical College, ul. Medyczna 7, 30-688 Kraków, Poland.
Received: 06 December 2024; Accepted: 09 December 2024; Published: 06 January 2025
DOI: 10.26502/jbb.2642-91280176
ShareProtein structure prediction remains an unsolved problem despite the significant progress brought by AI. Methods are evolving through reference to the growing resources of protein structure databases. However, the mechanism of the protein folding process remains unrecognised. A model was proposed to explain the incompatibility of the predicted structure with the native one. The reason for this is the influence of the environment, which, in a differentiated manner, directs the folding process to a form appropriate for the environment's characteristics. This is demonstrated using the lactococcin protein as an example. The structure of this protein was proposed for the targeting that the polar water environment brings. The native structure of this protein appears to represent the reconstitution of a distinct external force field with a hydrophobic component absent. This protein is active in the periplasmic space. The effect of this different environment on the folding process and the resulting different final structure was revealed. The results obtained allow a quantitative assessment of the environment affecting the location of energy minima to be introduced into the funnel model. The individual differentiated energy minimum imaged in the funnel model were placed on a scale expressing the specificity of the local environment of the folding protein.
Polar environment; Protein folding; In silico folding; Funnel model
Polar environment articles; Protein folding articles; In silico folding articles; Funnel model articles
An experiment called CASP (Critical Assessment of Structure Prediction [1]), which has been running since 1996 (in two-year cycles), provides an opportunity to assess progress in the field of protein structure prediction based on amino acid sequence. Participants in this competition propose 3D structures - referred to as models - for an amino acid sequence selected from a list of targets (target - a protein with a native structure known to the organisers). To date, this project has not fully answered questions about the mechanism of the protein folding process [2-19]. The ab initio models (models that are not referenced to available structures in the PDB database) currently used are based on the search for an energy minimum using molecular dynamics simulations [20-24]. A method that has recently made significant progress based on artificial intelligence (AI) techniques provides satisfactory results but with certain limitations [25-27]. In fact, numerous models obtained by known computer programs provided in CASP13 with a level of > 90 on the GDT_TS scale (Global Distance Test - Total Score) are accompanied by examples with a low score [28]. Combining the efforts of the leading teams in the WeFold project provided results with a comparable level of correctness for each team separately [29]. The question relates to the reasons for failure when using software with a high degree of correctness in the other models provided. The simplest answer to this question is that if the model used works very well in some cases and fails in others, it means that the proteins differ in the mechanism by which they have been folded. According to the observations obtained with the FOD-M model, the reason for this lies in the varying environmental conditions, which, by actively participating in the folding process, direct its course in a manner specific to that environment. The high degree of structural heterogeneity is clearly a result of the diversity of amino acid sequences. However, a second factor influencing the variation of the folding process is environmental manipulation. In the polar water environment, an adaptation in the form of exposure of hydrophilic residues on the surface with isolation of hydrophobic residues in the centre of the molecule is evident. In a hydrophobic membrane environment, an exposure of hydrophobic residues on the surface with polar residues in the centre is expected (especially in the case of transmembrane ion channels). These two extreme examples of radical change in the specificity of the local environment suggest the presence of a whole spectrum of changes in the level of polarity of the environment with intermediate characteristics.
In the model discussed in this paper, the structure obtained by simulated folding in an environment conducive to the generation of a micelle-like structure (hydrophobic core at the centre of the protein molecule with exposure of polar residues) is analysed against the native structure, which represents a system with a diffuse distribution of hydrophobicity. The two structural forms, while differing, reveal adjustment to different external conditions. The example discussed here allows us to propose a funnel model that takes into account the factor influencing the selection of an appropriate minimum in a multiple minima problem [30]. The main point of the analysis presented here is to introduce the environmental factor (external force field) as a factor that directs the folding process in a manner specific to that environment. Comparison of the structure obtained for polar water environment conditions with the native structure that functions in a periplasmic environment reveals structural differences. The discussed chain folding in a water environment is not able to obtain a native structure, which features a significantly different external environment.
Data
The two structures of the lactococcin 972 protein from Lactococcus lactis subsp. Lactis. One structure is the 2LGN available in PDB ID ID [31]. The other is a model obtained using the UNRES program with an additional module directing the in silico folding process towards the generation of a hydrophobic core and the exposure of polar residues on the surface. In this procedure, optimisation of energy is accompanied by optimisation of the ordering of hydrophobicity towards a micelle-like distribution. A detailed description is provided in [32].
Model applied
A model called fuzzy oil drop in its modified version of FOD-M [33] was used to describe both structures. The model assumes that polypeptide chains made up of bipolar molecules - amino acids - during the folding process in a polar water environment tend towards a micelle-like structure. This means aiming for a concentration of hydrophobic residues in the centre of the molecule with exposure of polar residues on the surface. Such a hydrophobicity distribution can be described by a 3D Gaussian function spanned over the body of the protein:

The size and shape of the protein is expressed by the adjusted values of the σX, σY and σZ parameters, respectively. The theoretical hydrophobicity level determined for a given effective atom calculated using eq. 1 is denoted hereafter as Ti.
The actual distribution may vary, as expressed by the function determining an inter-amino acid interaction. Such interaction depends on the distance between the effective atoms (the averaged position of the atoms comprising a given amino acid) [34]:

Where rij - distance between the positions of the effective atoms and c - cutoff distance adopted from [34] 9 Å. Hr means intrinsic hydrophobicity - any scale may be applied [35]. The observed hydrophobicity level determined for a given effective atom calculated using eq. 2 is denoted hereafter as Oi.
Both distributions after normalisation can be subjected to comparative analysis using Kullback-Leibler divergence entropy [36].

Where the tested distribution P (in the model used here - distribution O), distribution Q - the reference distribution (in the model used here - distribution O).
The value of DKL thus determined cannot be interpreted. Therefore, a second reference distribution R with values Ri=1/N is introduced, where N is the number of amino acids in the polypeptide chain is introduced. Such a distribution is opposite to the T distribution (hydrophobic core), representing a uniform distribution throughout the protein body.
The two values describing the status of a protein can be replaced by a single RD (relative distance) parameter:

A value of RD < 0.5 indicates "proximity" of the O distribution to the T distribution and thus denotes a protein with a hydrophobic nucleus present. RD > 0.5 indicates structuring without the presence of a hydrophobic nucleus.
The aqueous environment is not the only one in which proteins exhibit biological activity. The membrane environment requires reverse distribution in the form of exposure of hydrophobic residues on the surface with a polar centre (mainly for ion channels). In this situation, the function describes the expected distribution of hydrophobicity.

However, it appears that the hydrophobicity distributions identified in membrane proteins show a distribution expressed as the following function:

Where Ti denotes the distribution expressed by the 3D Gaussian function modified by the presence of the function TMAX - Ti to the extent expressed by the parameter K.
Proteins with K=0.0 (distribution compliant with the distribution expressed by the 3D Gaussian function) and K < 0.4 are down-hill, fast-folding and ultra-fast-folding, including type II antifreeze proteins [37]. The vast majority of domains treated as individual structural units show a structure described by low RD and K values [38]. In contrast, membrane proteins show values in the range of 0.9 < K < 3 [39].
A graphical representation of the model in question (FOD-M) is visualised in Figure1.
A- idealised distribution (Gaussian function) - T - blue line, observed distribution O - pink line, distribution R - grey line.
B- the RD value determined for the set of distributions as shown in A is 0.702.
C- the value of the K parameter determined for the minimum value of DKL(O|M) designated for the set of profiles shown in A is 0.5.
D- a compilation of the distributions T (blue), O (pink) and M (cyan) for the designated optimum K.
The value of K is determined as that value of this parameter for which DKL(O|M) takes a minimum (Fig. 1B). This procedure provides the form of the modified T distribution (M distribution) closest to the O distribution.
Protein status is determined by two parameters:
The program used is UNRES, which comes with a micelle-like environment introduction module. The UNRES (UNited RESidues) program is a significant tool in the field of protein structure prediction. A detailed presentation of this program is given in [24, 29, 32].
The analysis presented here aims to propose a model for in silico protein folding that takes into account the diverse environment.
The status of the structures in question: a native protein and a model for an identical amino acid sequence obtained using the in silico method is presented in Tab. 1. The structure of the model was obtained by means of a procedure which, in addition to optimising the internal energy, also includes a procedure which directs the folding towards the generation of a hydrophobic (micelle-like) nucleus [32].
|
LACTOCOCCIN |
RD |
K |
σX |
σY |
σZ |
|
Native form (2LGN) |
0.610 |
0.6 |
9.18 |
6.16 |
5.67 |
|
In silico model |
0.364 |
0.2 |
8.63 |
7.58 |
6.33 |
Table 1: Parameters expressing the status of the discussed structures: RD and K - parameters describing the hydrophobicity distribution in protein, parameters σX, σY and σZ, express the size and shape of molecule (3D Gauss function parameters) .
A compilation of parameter values (Tab. 1) identifies the model as a globular protein with a near spherical shape (3D Gaussian function parameters comparable). The status of the model clearly indicates the presence of a centric hydrophobic nucleus. A low K parameter value indicates shaping directed by the polar environment. The parameter values for the native structure reveal its elongated shape (Sigma parameter values varied). A high RD value indicates a distribution that is rather close to uniform (R distribution). A relatively high K value indicates the proportion of the environment with a significant presence of a hydrophobic agent.
A - set of T, O and M profiles for the value of K=0.6 together with a 3D presentation of the native structure of this protein.
B -T (blue), O (red) and M (green) profiles for the value of K=0.2 together with a 3D presentation of the model structure determined for this protein.
Axis in A and B: vertical - hydrophobicity, horizontal - residues according to sequence.
The comparison of the T, O and M profiles for the native protein suggests significant differences between the T distribution and the O distribution. The optimally fitted M distribution also reveals significant differences to the T distribution.
The high agreement of the T, O and M distributions with a low K value highlights the micelle-like organisation of hydrophobicity in the native form. It is also evident that the globular form of the resulting model is clearly visible, in contrast to the native form showing significant elongation (Sigma parameter values quantify this Table 1.).
The data presented in Tab. 1 allows the results obtained to be presented using the funnel model. The model visualises the multiple minima problem in relation to the energy states. The modification of the funnel model introduced using the FOD-M model is to introduce a quantitative assessment of the factor expressed by the horizontal axis (X-axis). In the visualisations of the model presented so far, the funnel is qualitative in character. In the visualisations of the funnel model to date, the horizontal axis is not described or is limited to qualitative concepts. According to the FOD-M model used, it simply expresses the value of the K parameter, i.e. the variability of the environment preferring different structuring (Figure 2).
In order to complement and expand the model examples, the status of the structures proposed for the T0953s2-D3 target have been summarised.
Using the interpretation based on the FOD-M model, the status of the individual models can be interpreted. A funnel model could also be proposed for this set, taking into account the value of the K parameter on the horizontal axis, which would indicate the contribution of external factors to the process of “folding” the models discussed. The force fields applied to obtain the forms in question can be interpreted as a force field tailored to the conditions expressed by the value of the K parameter.
In the presentation of Fig. S1, revealing the status of the models as expressed by the values of the RD and K parameters, specific groups with common characteristics are distinguished:
On the vertical axis in the funnel model, there is energy value. Here, the value of the RD parameter, which can also be interpreted as a sort of energy state associated with the presence of a hydrophobic core, is given as a factor characterising the status of a given model (model structures indicated by the red oval in Fig. S1). The presence of a hydrophobic core, together with the possible presence of disulphide bonds, is regarded as a factor stabilising the tertiary structure of proteins. Others, lacking or having a deformed hydrophobic core, can be interpreted as less stable. Energy status is not the only criterion for structural correctness. Structural correctness is the assurance of biological activity, which does not necessarily represent an energy minimum. According to the FOD-M model, a stable structure is one with a low internal energy level but adapted to the conditions directing a given folding process.
It can be speculated that the force fields used by the participants providing models for the T0953s2-D3 target in CASP13 are parameterised in a specific, different way favouring a particular type of folding.
Figure 4: Summary of the RD and K parameter values for
A- the models provided in the CASP13 project - T0953s2-D3 target. The colour variation classifies the protein structures described in the text - blue pentagon - the status of target protein
B- the real proteins listed in Tab. 2. The blue squares - membrane proteins.
A very interesting point representing extremely high RD and K values is a structural form that is unfolded to a considerable extent (Fig. 4). Theoretically, a structure of this type would have to be in highly altered (in relation to polar water) external conditions preventing the formation of a packed globular form.
|
PDB ID |
FUNCTION |
RD |
K |
PDB ID |
FUNCTION |
RD |
K |
|
3F6C[40] 3F13[41] 3J39[42] 3S30[43] 3UOE[44] 3T6G[45] 2JG7[46] 4HE8[47] 2ZTG[48] 1W8X[49] 2RIN[50] |
DNA binding hydrolase ribosome prot. hydrolase oxidoreductase signaling oxidoreductase oxidoreductase ligase capsid virus transport |
0.380 0.388 0.419 0.501 0.589 0.569 0.630 0.767 0.606 0.686 0.674 |
0.17 0.18 0.21 0.42 0.50 0.59 0.73 0.82 0.81 0.92 1.02 |
2F1M[51] 2OBD[52] 1LQS[53] 2RHS[54] 3CSL[55] 3SYS[56] 2XA7[57] 3FIP[58] 3DWO[59] 3J0C[60] |
transport transport transport ligase membrane membrane transport transport membrane structural |
0.735 0.763 0.745 0.741 0.761 0.751 0.708 0.751 0.819 0.861 |
1.05 1.21 1.25 1.44 1.56 1.74 1.84 1.85 1.91 2.06 |
Table 2: The set of randomly selected proteins representing the spectrum of K values.
The target T0953s2-D3 (CASP13) discussed in detail in [28] is an example of a set of models representing the full spectrum of environmental conditions (a detailed description of these models is given in Supplementary Materials). A protein target characterised by extremely low RD and K values, representing a micelle-like distribution, turned out to be a target far from the typical parameterisation of the numerical techniques used, which in their majority represent parameterisations corresponding to a status with an RD close to 0.6 [28]. This value is probably due to the fact that the techniques for the force field parameterisation used in the applied models are based on the analogy of structures available in databases including the PDB. According to our analyses (unpublished data), most proteins represent such a level on the K parameter scale in their large numbers.
The funnel model shape proposals available in the literature treat structural variation qualitatively, in contrast to the model proposed here, where the distribution of the energy minimum is attributed to the corresponding value of the parameter characterising the environment for the folding process [61-66].
It is assumed that environmental changes leading to amyloid transformation will enable simulation of this transformation by way of gradually changed environmental specificity by introducing stepwise variation of the K parameter considered in the form of external force field record [67]. This type of analysis has been positively validated using the example of membrane proteins, whose specific form of membrane anchoring has been shown to be varied [39].
The mathematical model, which shall be applied for simulation In Silico of folding process is the multi-object optimisation [68]. The energy, the minimum of which is searched shall take into account the dependence on two functions as shown below:

Where: F(rij) the function expressing the energy status of protein as dependent on two functions: fINT (rij) - internal force field generating the status dependent on the rij distances between interacting atoms according to non-bonding interactions and fEXT(K, rij) - external force field depending also on rij distances and on K value modifying the external conditions (eq. 6) which influences the structuralisation process of protein. The two f functions act contradictory. The minimisation fINT (rij) may cause the increase of fEXT(K, rij) values and vice versa. The method called front Pareto is the one specially oriented on finding the solution represented by the consensus between these two functions [69].
The funnel model expressing the multiple minima problem can be solved by introducing conditions imposed by the external force field. Correct models provided for subsequent targets (CASP project) are obtained for those proteins, for which the force field parameterisation (or other evaluation criteria) is consistent with the environmental conditions. Obtaining the correct structure for this protein via the In Silico route definitely requires the varying environmental influences to be taken into account. It should be noted that the protein in question is active in the periplasmic space. This environment - as this example shows - influences structure formation in a different way to the aquatic environment. This is made apparent by the funnel model form limiting the potential for structural variation to a narrow range of the K parameter reflecting the environmental conditions for the folding protein. The innovation proposed in this analysis is to give a quantitative form to the horizontal axis in the funnel model, which in other forms of this model is defined qualitatively. The conclusion regarding simulation In Silico of protein folding process, it is necessary to treat such procedure as multi-object optimisation. The minimised energy expressing the status of protein shall be treated as the function dependent on two functions: internal force field (non-bonding interactions) and external force field (specific influence of environment expressed by K parameter). The consensus for these two functions shall be found for final structure of folded protein.
Many thanks to Anna Śmietańska and Zdzisław WiŚniowski for technical support. This research was partially supported by the European Union's Horizon 2020 programme under grant Sano No. 857533. and the Sano project carried out within the International Research Agendas programme of the Foundation for Polish Science, co- financed by the European Union under the European Regional Development Fund.
We gratefully acknowledge Poland's high-performance Infrastructure PLGrid ACK Cyfronet AGH for providing computer facilities and support within computational grant “plgdrippy”.
Author contributions
Conceptualisation: IR, LK; Advisory: IR; Calculation : MS
Funding
This research was funded by Jagiellonian University Medical College grant number N41/DBS/001127.
Data availability
All data can be made available via a request addressed to the corresponding author. The program enabling the calculation of RD is accessible on the GitHub platform: https://github.com/KatarzynaStapor/FODmodel and on the platform https://hphob.sano.science
Data presented in Supplementary Materials are attached to this submission - Excel file with data and profiles for proteins described in Tab. 2.
Declarations
Ethics approval and consent to participate
Neither human nor animal data are used in the presented research.
Consent to publication
Not applicable.
Competing interests
The authors declare no conflict of interest.