Structure and folding prediction

1Proteins¶

Let’s consider a 100-residue peptide chain with a unique folded state. If each amino acid had only two available conformations, the number of configurations available to the chain would be $2^{100} \sim 10^{30}$. If the time required to switch between any two configurations was 10^-12 s (a picosecond), and we assume that no configuration is visited twice, it would take approximately 10¹⁸ seconds to explore all the “phase space” (and therefore, on average, to fold): this is close to the age of the universe! In reality, proteins fold on time scales ranging from microseconds to hours. This is the gist of the famous “Levinthal’s paradox”, which implies that the search for the folded structure does not happen randomly, but it is guided by the energy surface of residue-residue interactions.

Proteins can be denatured by changing the solution conditions. Very high or low temperature and pH, or high concentration of a denaturant (like urea or guanine dihydrochloride) can “unfold” a protein, which loses its solid-like, native structure, as well as its ability to perform its biological function. Investigations of the unfolding of small globular proteins showed that, as the denaturing agent (e.g. temperature or denaturant concentration) increases, many of the properties of the protein go through an “S-shaped” change, which is characteristic of cooperative transitions.

Furthermore, calorimetric studies show that the denaturation transition is an “all-or-none” transition, i.e. that the “melting unit” of the transition is the protein as a whole, and not some subparts. This of course applies to single-domain small proteins, or to the single domains of larger proteins. Here “all-or-none” means that the protein exist only in one of two states (native or denaturated), with all the other “intermediate” states being essentially unpopulated at equilibrium. Therefore, this transition is the microscopic equivalent of a first-order phase transition (e.g. boiling or melting) occurring in a macroscopic system. Of course, since proteins are finite systems and therefore very far from the thermodynamic limit, this is not a true phase transition, as the energy jump is continuous, and the transition width is finite.

But how does the denatured state look like? As discovered by using a plethora of experimental methods, which often seemed contradicting each other, the answer depends on the denaturing conditions. Figure 1(a) shows the phase diagram of lysozyme (a globular small protein) for different temperatures and concentrations of a denaturing agent. Increasing the latter leads to a transition to a disordered coil state where essentially no secondary structure is present. By contrast, the effect of temperature is qualitatively different, as the protein partially melts in a state (the molten globule, MG) that retains most of its secondary structure, but it is not solid like and it cannot perform any biological function. A schematic of the difference between the native and molten globule states is shown in Figure 1(b).

The cooperative transition that leads to the molten globule is due to the high packing of the side chains, which is energetically favourable but entropically penalised, since it impedes the small-scale backbone fluctuations that trap many of the internal degrees of freedom of the protein. The liberation of these small-scale fluctations, and the resulting entropy gain, requires a slight degree of swelling, which is why the MG is not much larger than the native state. However, such a swelling is large enough to greatly decrease the van der Waals attractions, which are strongly dependent on the distance, and to let solvent (water) molecules in the core.

In general, not all proteins behave like this, as there exist some (usually small) proteins that unfold directly into a coil, others that form the molten globule under the effect of specific denaturants, and coils under the effect of others, etc. However, the unfolding transition has nearly always an “all-or-none”, cooperative nature.

Generally speaking, the denaturation of a protein is a reversible process, provided that some experimental conditions are met (the protein is not too large, its chemical structure has not been altered by the denaturation process or by post-translational modifications, etc.). The fact that protein folding and unfolding is a reversible process implies that all the information required to build a protein is stored in its sequence, and that its native structure is thermodynamically stable. Therefore, the folding process can be, in principle, understood by using thermodynamics and statistical mechanics. The sequence $\to$ structure hypothesis is sometimes called Anfisen’s dogma, named after the Nobel Prize laureate Christian B. Anfisen, whose studies on refolding of ribonuclease were key to establish the link between “the amino acid sequence and the biologically active conformation” (Anfinsen et al. (1961), Anfinsen (1973)).

Figure 2(a) shows the origin of the “all-or-none” transition observed in protein denaturation in thermodynamic terms. Here I will use the description given by Finkelstein & Ptitsyn (2002) (whence the figure comes), almost word by word. The panel shows the energy $E$ and the entropy $S$ of a protein as a function of its internal density (in some arbitrary units). The energy is at its minimum at close packing, while the entropy S increases with decreasing density. Initially, this increase is slow, since the side chains can only vibrate, but their rotational isomerization is impossible. When the density decreases enough that the latter becomes possible the increase is much steeper. Finally, when free isomerization has been reached the entropy growth becomes slow again.

The nonuniform growth of the entropy is due to the following. The isomerization needs some vacant volume near each side group, and this volume must be at least as large as a $CH_3$ group. And, since the side chain is attached to the rigid backbone, this vacant volume cannot appear near one side chain only but has to appear near many side chains at once. Therefore, there exists a threshold (barrier) value of the globule’s density after which the entropy starts to grow rapidly with increasing volume (decreasing density) of the globule. This density is rather high, about $80\%$ of the density of the native globule, since the volume of the $CH_3$ group amounts to $\approx 20\%$ of the volume of an amino acid residue. Therefore, the dependence of the resulting free energy $F = E - TS$ on the globule’s density ρ has a maximum (the so-called free-energy barrier) separating the native, tightly packed globule (N) from the denatured state (D) of lower density. It is because of the presence of this free energy barrier that protein denaturation occurs as an “all-or-none” transition akin to melting or nucleation, independent of the final state of the denatured protein.

It is useful, when talking about proteins, to also introduce the concept of random heteropolymers, which are polymers with sequences that have not been selected by nature, and therefore have many (degenerate or quasi-degenerate) ground states rather than a single native conformation. Figure 2(b) shows the qualitative difference between a protein and a random heteropolymer. The bottom panels of the figure show the typical energy spectra of the two chains, where each line corresponds to a single conformation. At high energy both spectra are essentially continuous. However, the low-energy, discrete parts look very different: the energy difference between the ground state and the continuous part of the spectrum, where the majority of the “misfolded” states lie, is much larger in a protein than in a random heteropolymer. It is this energy gap that creates the free-energy barrier of Figure 2(a).

We can also define the entropy $S(E)$, which is proportional to the logarithm of the number of structures having energy $E$, and the temperature corresponding to the energy $E$, implicitly defined by the relation $d S(E) / dE = 1 / T(E)$. As sketched in Figure 2(b), the peculiar shape of the protein energy spectrum gives raise to a concave part of $S(E)$, which means that there is a temperature at which the lowest-energy, native conformation coexists with many high-energy (denatured) ones. In this simple picture, this is the melting temperature $T_M$. This temperature should be compared with the so-called “glass temperature” $T_G$, called $T_C$ in the figure, which is the temperature associated to the lowest energy of the continuous band of the spectrum. Above $T_G$ the chain behaves, from the kinetic point of view, like a viscous liquid, while below the glass temperature the behaviour is glassy-like (i.e. non-Arrenhius dependence on temperature, aging, non ergodicity, etc.).

If $T_M < T_G$, the slow dynamics asssociated to the glass state makes it kinetically impossible to reach the ground state, and the chain remains trapped in higher-energy misfolded states. By contrast, if $T_M > T_G$, folding is possible, and results obtained with minimalistic models show that the fraction $T_M / T_G$ is a good proxy for “foldability”: the larger this fraction, the faster the protein folds.

All these concepts are not merely qualitative, but have been grounded in theory by using sophisticated statistical mechanics approaches that have generated the “funnel landscape” folding picture, which is schematically shown in Figure 3(b).

The basic idea is to leverage the concept of frustration, which happens when the interactions between the different parts of a system are such that they cannot be satisfied, from the energetic point of view, all at the same time. If this is the case, then there is no single ground state, but rather many low-energy stable states that, in a many-body system, are uncorrelated from each other and separated by (possibly large) energy barriers. An example is provided in Figure 3(a). The resulting “energy surface” (or landscape) is said to be rugged (or rough), and generates a dynamics where the system sits in a valley for a long time before being able to jump over the barrier and move to a different stable conformation. This is a hallmark of glassiness.

To draw a parallel with proteins, we can consider a random heteropolymer, where the interactions among the different amino acids will be frustrated, blocking the system from finding a single well isolated folded structure of minimum energy. A candidate principle for selecting functional sequences is thus the minimization of this frustration. Of course even in the absence of frustration, there are still energetic barriers on the path towards the native conformation due to local structural rearrangements which still give raise to a rugged landscape, slowing down the kinetics towards the folded state. This scenario has come to be called a folding “funnel”, emphasizing that there is a single dominant valley in the energy landscape, into which all initial configurations of the system will be drawn.

A classic funnel diagram, such as the one shown in Figure 3(b), represents an energy-entropy landscape, with the width being a measure of the entropy, whereas the depth represents both the energy and two correlated order parameters $Q$ and $A$, which are the fractions of native contacts and correct dihedral angles in the protein backbone. Although in reality the landscape is multidimensional, here the projection attempts to retain its main features, such as the barrier heights, which are a measure of the “ruggedness” or “roughness” of the landscape. The typical height of these barriers, together with the energy gap and the number of available conformations, are the three main parameters of the Random Energy Model, which is one of the main theoretical tools in this context.

1.1The Zimm-Bragg model for the helix-to-coil transition¶

As we discussed, polypeptides can, in general, be reversibly denatured (with temperature or pH, for instance). Here we theoretically investigate this phenomenon in the context of secondary structure with the Zimm-Bragg model, which is a toy model for a chain that can exhibit a continuous coil-to-helix transition. Note that it is possible to synthesise polypeptide chains that undergo such a transition, see e.g. Doty & Yang (1956).

Consider a polymer chain composed of $N$ residues. Each residue can be in one of two states, “coil” (c) or “helix” (h), so that a microscopic state of the chain can be expressed as a one-dimensional sequence, like ccchhhhhcchhhchcc. In equilibrium, the probability (or statistical weight) to observe a given chain microstate $c$ is just the Boltzmann factor, $e^{-\beta \Delta F_c}$, where $\Delta F_c$ is the free-energy cost of the microstate, so that the total (configurational) partition function is

Since (free) energies are always defined up to a constant, we define the free-energy of a coil of any length to be zero by definition, so that the statistical weight of a coil residue is always 1. By contrast, forming a sufficiently long helical segment should be advantageous, otherwise no transition would occur. However, we know that in the helices formed by polypeptides, each $j$-th amino acid is hydrogen-bonded to the $(j + k)$-th, where $k = 4$ for a α-helix. As a consequence, since one full helical turn needs to be immobilised before seeing any benefit from the formation of hydrogen bonds, the formation of a helix incurs a large entropic penalty that needs to be paid every time a coil segment turns into a helical one, $f_{\rm init}$. However, once a helix is formed and the initial price has been paid^[1], adding an H contributes a fixed amount $f_h$ to the helix free energy. As a result, the free-energy contribution of a helical stretch of size $n$ is $\Delta F = f_{\rm init} + n f_h$, so that is statistical weight is

where we have defined the cooperativity (or initiation) parameter, $\sigma \equiv \exp(-\beta f_{\rm init})$, and the helix-elongation (or hydrogen-bonding) parameter, $s \equiv \exp(-\beta f_h)$. Note that, since $f_{\rm init}$ is purely entropic, $\beta f_{\rm init} = -S_{\rm init} / k_B$ and therefore σ does not depend on temperature.

Following Finkelstein & Ptitsyn (2016), we define two accessory quantities: $C_i$ is the partition function of a chain of size $i$, where the last monomer is in the coil state, and $H_i$ is the partition function of a chain of size $i$, where the last monomer is in the helix state. With this definition, the partition function of the chain is $Q_i = C_i + H_i$. We can explicitly compute the first terms:

• For $i = 1$, the sequence is either c or h, so that $C_1 = 1$ and $H_1 = \sigma s$.
• For $i = 2$, we have two sequences ending with c, which are hc and cc, and two sequences ending with h, which are hh and ch. Each partition function is given by a sum of the statistical weights of the allowed microstates, so that we have $C_2 = 1 + \sigma s$ and $H_2 = \sigma s^2 + \sigma s$.
• For $i = 3$, there are 8 available sequences, which are built by taking each $i = 2$ sequence and adding either c or h at its end. Each of the resulting sequence will have a statistical weight that is that of the base $i = 2$ sequence, multiplied by 1 if the sequence ends with c, and by $s$ or $\sigma s$ if the sequence ends with h and the preceeding character is h or c, respectively. Summing up the contribution we obtain $C_3 = C_2 + H_2 = 1 + 2 \sigma s + \sigma s^2$ and $H_3 = C_2 \sigma s + H_2 s = \sigma s + \sigma s^2 + \sigma^2 s^2 + \sigma s^3$.

In general the operations carried out to compute the $i = 3$ partial partition functions can be applied to any other value of $i$, meaning that we can write down recursive relationships connecting $H_{i}$ and $C_{i}$ to $H_{i-1}$ and $C_{i-1}$. These relationships can be neatly expressed in matricial form:

It is convenient to define the matrix

so that the the total partition function of a chain of $N$ monomers can be written as^[2]

where the first vector-matrix multiplication on the right-hand side gives $(C_N, H_N)$, which multiplied by the $(1, 1)$ column vector yields $C_N + H_N$. The partition function can be written in terms of the eigenvalues of $\matr M$ if we note that

where $\matr \Lambda = \matr A^{-1} \matr M \matr A = \begin{pmatrix} \lambda_1 & 0 \\ 0 & \lambda_2 \end{pmatrix}$, and $\lambda_1$ and $\lambda_2$ (with $\lambda_1 > \lambda_2$) are the eigenvalues of $\matr M$, which can be found by solving the secular equation, viz.:

whence

Note that the two eigenvalues are linked by the relation $\lambda_1 + \lambda_2 = 1 + s$.

Since they can be of any length, the two eigenvectors $\vec{A}_1$ and $\vec{A}_2$ of a 2x2 matrix each have a single independent element, $a_k$ (where $k = 1, 2$), which is given by

We obtain $a_1 = \lambda_1 - s = 1 - \lambda_2$ and $a_2 = \lambda_2 - s = 1 - \lambda_1$, and therefore $\vec{A}_1 = \begin{pmatrix} 1 - \lambda_2\\1\end{pmatrix}$ and $\vec{A}_2 = \begin{pmatrix} 1 - \lambda_1\\1\end{pmatrix}$. The two eigenvectors are used to construct the diagonalising matrix and its inverse^[3]

We can now write Eq. (7) explicitly in terms of $\lambda_1$ and $\lambda_2$:

From the partition function we can compute any observable we need. The most interesting quantity is the average fraction of residues that are in the helical state h, which can be compared to experiments. First of all, we note from the $C_i$ and $H_i$ we explicitly calculated and from Eq. (5) that the partition function is a polynomial in $s$, and therefore can be written as

where $g_k$ are coefficients that do not depend on $s$, but only on σ and on the multiplicity of each state with $k$ helical residues. Therefore, by definition the probability that the chain has exactly $k$ helical residues is $P_N(k) = \frac{g_k s^k}{Q_N}$, which means that the average number of helical residues is

since $\frac{ds^k}{ds} = k s^{k-1}$ and only $s^k$ depends on $s$. Defining the fraction of helical residues, $\theta_N \equiv \langle n_h \rangle / N$, and exploiting the known property of logarithms we find

Substituting Eq. (13) it is possible to write $\theta_N$ in the following closed, albeit rather cumbersome, form :

For sufficiently long chains, $(\lambda_2 / \lambda_1)^N \ll 1$ and the expression can be greatly simplified:

The model as described and solved takes into account the cooperativity of the transition through the parameter σ. It is instructive to look at the non-cooperative solution, which can be obtained by setting $\sigma = 1$. In this case $\lambda_1 = 1 + s$ and $\lambda_2 = 0$, so that

which is independent of $N$ and describes the simple chemical equilibrium of a $c \rightleftharpoons h$ reaction. Indeed, defining the equilibrium constant $K_{hc} \equiv [h] / [c]$, where $[x]$ is the concentration of $x$, we find that the fraction of residues in the helical state is

which, if compared with Eq. (19), shows that $s = K_{hc}$.

A comparison between the Zimm-Bragg theory and experiments on poly-γ-benzyl-L-glutamate of different lengths are reported in Figure 4. It is clear that there is a strong $N$-dependence that shows the cooperative nature of the transition. Moreover, in fitting the curves it is assumed that σ does not depend on $T$, since it has a purely entropic origin.

1.2The HP model¶

Globular proteins exhibit a unique characteristic compared to other isolated polymer molecules: their native structure is unique, meaning that the path formed by the backbone is essentially fixed, bar some small fluctuations. The question arises: what forces, determined by the amino acid sequence, contribute to this remarkable specificity? Of course, as discussed before, a significant constraint is provided by the fact that native structures are compact. However, even for compact chain molecules, many conformations are possible, with the number of maximally compact conformations increasing exponentially with chain length. Among the various physically accessible compact conformations, there are in principle many different interactions that can impact their thermodynamic stability and thus select the native structure. However, it is now well accepted that the interaction type that most contribute to this conformational reduction is hydrophobicity: given a specific amino acid sequence, the number of compact conformations that it can take is greatly reduced by the constraint that residues with hydrophobic side chains should reside inside the globule, while polar and hydrophilic residues should be on the protein’s surface (see e.g. Finkelstein & Ptitsyn (2002) and Dill (1990)).

A simple model that has been used to demonstrate the importance of hydrophobic interactions in selecting the native structure is the Hydrophobic-Polar (HP) protein lattice model. In this model, proteins are abstracted as sequences composed solely of hydrophobic (H) and polar (P) amino acids. In the HP model, the protein sequence is mapped onto a grid, or lattice, which can be two-dimensional (2D) or three-dimensional (3D). Common lattice types used include square lattices for 2D models and cubic lattices for 3D models, but other choices are possible.

Here I will present the 2D square lattice version of the model^[4]. Each of the $N$ amino acids of the protein is assigned a type: “H” means hydrophobic, and “P” means polar, so that the sequence of the protein is then the ordered list of amino acids (e.g. “HHPHPPPHHP”). The $i$-th amino acid in the sequence is identified by its index $i$ and occupies a point on this grid, and two amino acids cannot occupy the same grid point. Two amino acids $i$ and $j$ for which $|j - i| = 1$ are “connected neighbours” and share a backbone edge representing the peptide bond, which on the square lattice can be either horizontal or vertical. By contrast, two amino acids $i$ and $j$ for which $|i - j| > 1$ that are adjacent in space but not along the sequence are said to be “topological neighbours”. Every pair of HH topological neighbours form a topological contact that contributes an amount of free energy $\epsilon \equiv \beta \Delta F_{HH} < 0$, while the other contacts give $\epsilon_{HP}$ and $\epsilon_{PP}$. In the simplest version of the model $\epsilon_{HP} = \epsilon_{PP} = 0$. In this case, the total (dimensionless) free energy of a conformation is

where $m$ is the number of HH topological contacts. Following Lau & Dill (1989), the lowest-energy (highest-$m$) conformations of a chain with a specific sequence and a length $N$ are called “native”.

The concept of topological contacts is also useful to describe the degree of compactness of a conformation, defined as

where $t$ is the number of topological contacts, irrespective of the type of neighbouring residues, and $t^{\rm max}$ is the largest possible number of topological contacts for a chain of size $N$, which is

where $P_{\rm min}$ is the perimeter of the smallest box that can completely surround the protein. By this definition, (maximally) compact conformations have $\rho = 1.0$.

As an example, Figure 5 shows the single lowest-energy conformation of a given sequence, together with two other compact conformations of higher free energy. It is rather evident that the native conformation has a core that is richer in hydrophobic residues compared to the other two.

In order to be more quantitative, some useful observables can be introduced. In general, the thermodynamic properties of a chain of length $N$ and a given sequence can be computed from the partition function:

where $\Omega_N$ is the number of distinct conformations of an $N$-chain, so that the sum runs over all the conformations. If the length of the chains is not too large^[5], it is possible to numerically perform a complete (exhaustive) enumeration of all possible chain conformations, from which any quantity of interest can be computed as an ensemble average, viz.

The partition function can also be written as a sum over conformations of fixed number of HH topological contacts $m$:

where $m^{\rm max}$ is the maximum number of HH topological contacts and $g(m)$ is the number of distinct chain conformations having $m$. Note that by this definition, the native conformations are those for which $m = m^{\rm max}$. We also define $G(m, u)$ as the number of conformations that have $m$ HH topological contacts and $u = t - m$ non-HH topological contacts, which is connected to $g(m)$ by

Since partition functions are defined up to a multiplicative constant, we multiply $Q$ by $e^{m^{\rm max} \epsilon}$ to obtain^[6]

With this definition it is straightforward to take averages over native conformations only. Indeed, if we take the $\epsilon \to -\infty$ limit, $e^{(m^{\rm max} - m) \epsilon}$ is non-zero only for $m = m^{\rm max}$, and the partition function becomes

We now define some useful average values to characterise the different states. The average compactness over all conformations is

and the average compactness over the native conformations is

where we used the fact that

The average energy is proportional to the average number of HH contacts, which is

Finally, we know already that most of the hydrophobic residues are buried in the core of globular proteins. In the HP model we define the core as the set of residues that are completely surrounded by other residues. Let the number of core residues be $n_i$, and the number of hydrophobic core residues be $n_{hi}$, then the degree of hydrophobicity of a given conformation can be estimated as

and its averages over all and native conformations can be computed as

The ensemble averages of these quantities, as well as of the partition function, which is a measure of the number of available conformations, as a function of $-\epsilon$ for two sequences are shown in Figure 6. Looking at the top panels it is clear that if the HH interaction is strong the preferred conformations of the HPHPPHPPHH sequence have low energy, are few ($Q \to 1$), compact, $\langle \rho \rangle = 1$, and have a purely hydrophobic core, $\langle x \rangle = 1$: this is a “folding” sequence.

By contrast, the bottom panels show a sequence that behaves very differently: the low-energy (native) conformations are not compact and large in number. Looking at the sequence itself, PPPPPHHHHH, it is clear that the hydrophilic “tail” does not make it possible to form a conformation with a hydrophobic core, and the number of HH topological contacts itself cannot be too high. These are all clear signs of a “non-folding” sequence.

Note that the quantities reported in Figure 6 are somewhat correlated: folding sequences tend tend to have few (sometimes one) low-energy compact conformations with hydrophobic cores, whereas non-folding sequences tend not to have any of these attributes.

The HP model (or its extensions) can be used to understand some general or even specific aspects of protein folding because of its simplicity. The 2-letter alphabet and the greatly reduced residue conformations due to the underlying lattice makes it possible to fully explore (at the same time) the conformational space (i.e. the set of all possible internal conformations of a molecule) and the sequence space (i.e. the set of all possible sequences of H and P residues). Both spaces have sizes that grow exponentially with $N$, and therefore the maximum length of chains that can be investigated by exhaustive enumeration is somewhat limited, but there are computational techniques that can be used to circumvent this issue.

1.2.1Folding HP sequences¶

If we now use a stricter definition of a “folding” sequence as a sequence that has a single native conformation, we can draw a comparison between lattice folding sequences and real proteins. First of all, how many folding sequences there are?

Figure 7 shows the distributions of the number of native conformations $g(m^{\rm max})$, i.e. how many sequences have the given degeneracy, for different chain lengths. As $N$ increases, the distribution becomes more and more peaked on 1: far more sequences have a single lowest-energy conformation rather than, say, 10 or 100. These results suggest that hydrophobicity alone is enough to greatly reduce the number of compact configurations into a small set of folding candidate structures (Dill (1990)), and confirms that the hydrophobic force is a major driving force for protein folding (see e.g. Finkelstein & Ptitsyn (2002)).

But to what extent lattice and real proteins are similar? I will present some results reported in Moreno‐Hernández & Levitt (2012). Therein, an exhaustive enumeration of all conformations of chains of lengths up to 25 has been carried out with two HP models:

• In HP100, which is the one described above, only HH contacts contribute to a conformation’s free energy.
• In HP211 topological contacts between H and P and P and P contribute -1 (i.e. $\epsilon_{HP} = \epsilon_{PP} = -1$) to the free energy, while HH contacts contribute $\epsilon = -2$. As a result, compared to the HP100 model, all residues feel an attraction that makes compact conformations more favourable.

Some interesting results pertaining to folding sequences are reported in Figure 8. Figure 8(a) shows, as a proxy for compactness, the normalised radius of gyration $R_g/\min(R_g)$, where $\min(R_g)$ is the smallest radius of gyration among the structures of a given length, for the lattice and real proteins. For real proteins the authors have considered a non-redundant set of 2401 single-domain proteins from the PDB with lengths between 50 and 200 residues long. By comparing the two lattice models, it is evident that the additional attraction between non-hydrophobic residues enhances the compactness in HP211 more, but the variation in compactness, given by the error bars, is very small in both cases, and the compactness itself does not depend on the chain length. The same behaviour is observed in real proteins.

Figure 8(b) shows the percentage of residues that are hydrophobic in the lattice structures and for the same real proteins used for the radius of gyration. Both lattice and real proteins have an average hydrophobicity of $\approx 50\%$, independent of $m$, and the error bars and total range (dashed lines) are comparable in the two cases (although both are a bit larger for the HP model). Folding sequences in the HP model are “protein-like”.

1.2.2Designable HP conformations¶

Now we focus on compact conformations. How many sequences have a given compact conformation as their native (i.e. lowest-energy) state? Does the answer depend on some properties of the selected conformation itself? In the context of the HP model, some of the first answers were provided here with a HP model with $\epsilon = -2.3$, $\epsilon_{HP} = -1$ and $\epsilon_{PP} = 0$, chosen in order to fulfill the following physical constraints, rooted in the analysis of real proteins:

1. Compact conformations have lower energies compared to non-compact ones.
2. Hydrophobic residues should have the tendency to be buried in the core ($\epsilon_{PP} > \epsilon_{HP} > \epsilon$).
3. Different types of monomers should tend to segregate ($2 \epsilon_{HP} > \epsilon_{PP} + \epsilon_{HH}$).

The output of the computation was to enumerate all the native structures of each possible sequence on a 3D 3x3x3 cube and on 2D 4x4, 5x5, 6x5 and 6x6 boxes. The result of this complete enumeration is the list of all possible sequences that can “design” a given structure, or, in other words, that have that structure as their unique native state. The size of this list, $N_s$ is a measure of the designability of a given structure.

Compact structures differ markedly in terms of their designability: there are structures that can be designed by a large number of sequences, and there are “poor” structures that can be designed by only a few or even no sequences. In fact, for $\approx 10\%$ of the conformations there is no sequence that has that structure as its ground state. The majority of the sequences ($\approx 95\%$ out 2²⁷ total sequences in 3D) have degenerate ground states, i.e. more than one compact conformation of lowest energy, which means that in 3D there exist almost 7 million sequences of length 27 that have a unique compact ground state. Moreover, the number of structures with a given $N_s$ value decreases continuously and monotonically as $N_S$ increases. The data for the 3D and largest 2D cases is shown in Figure 9(a), where the long tails of the distributions, with some structures being the ground states of thousands of sequences, highlight the presence of “highly-designable” compact conformations.

Structures with large $N_S$ exhibit specific motifs (i.e. secondary structures) that small $N_S$ compact structures lack. For instance, the 3D compact structures with the 10 largest $N_S$ values contain parallel running lines packed regularly and eight or nine strands (three amino acids in a row), sensibly more than the average compact structure. Figure 9(b) shows the most designable structures in 3D and in 2D, where it is easy to spot the regularity of the “secondary structures” and, for the 3D structure, also the “strands”.

With the simple lattice model is also possible to directly assess the effect of mutations, which in real proteins is particularly important in the context of homologous sequences (sequences related by a common ancestor). Indeed, focussing on highly designable structures and referring to the $N_s$ different sequences that fold into them as “homologous”, it is possible to observe phenomena that are qualitatively similar to those observed in real proteins. For example, sequences that differ by more than half of their residues can design the same structure (see Figure 20, which will be discussed later, for a real-protein example). Looking at the effect of mutations, the rightmost panel of Figure 9(c) shows that some in very designable conformations residues are highly mutable, whereas others are highly conserved, with the conserved sites being those sites with the smallest or largest number of sides exposed to water.

Another important feature of proteins that is also reproduced by the model is the large energy gap between the lowest-energy and first excited conformations, which stabilises the native structures from the thermodynamic point of view. In the lattice model, the average energy gap $\bar \delta_S$ is defined as the minimum energy required to change a native conformation to a different compact structure, averaged over all the $N_S$ sequences that design that native conformation. Figure 10 shows that there is a sudden jump (in 3D) or change of slope (in 2D) at a specific value of $N_S$. Therefore, the conformations that are highly designable also have large energy gaps, making them highly stable from a thermodynamic point of view.

As shown in the rightmost panel of Figure 10, the same behaviour is also observed when a 20-letter (rather than 2-letter) alphabet is used, with the interactions between the residues being modelled using the MJ interaction matrix, which was derived by analysing the contact frequencies of amino acids from protein crystal structures.

Although the ingredients in the simple lattice model are minimal, lacking many of the factors that contribute to the structure of proteins, the results that have been obtained with it contributed to the understanding that there is a link between designability and thermodynamic stability: “From an evolutionary point of view, highly designable structures are more likely to have been chosen through random selection of sequences in the primordial age, and they are stable against mutations” (Li et al. (1996)).

1.3Sequence alignment¶

Simple models are useful to understand the underlying physics of some particular phenomena. However, how can we understand something very specific, like what is the 3D structure of a particular sequence? The simplest way is to look for similarities: if we already have a list of sequence $\to$ structure connections, we can try to look whether the new sequence, for which the 3D structure is unknown, is similar, and to what degree, to another for which the 3D structure is already known. This operation is called “sequence alignment” (SA).

Sequence alignment is a fundamental technique in bioinformatics. The primary goal of sequence alignment is to identify regions of similarity that may indicate functional, structural, or evolutionary relationships between the sequences being compared. In the context of proteins, sequence alignment can be useful for reasons that go beyond the 3D structure prediction:

• It helps predict the function of unknown proteins based on their similarity to known proteins.
• It aids in understanding evolutionary relationships by identifying conserved sequences among different species, allowing the construction of phylogenetic trees.
• It can identify conserved domains that are crucial for the function of proteins.
• It helps in identifying potential drug targets by finding unique sequences in pathogens that differ from the host.
• It assists in annotating genomes by finding homologous sequences, thus providing insights into gene function and regulation.

Turning a biological problem into an algorithm that can be solved on a computer requires making some assumptions in order to obtain a model with relative simplicity and tractability. In practice, for our sequence alignment model we ignore realistic events that occur with low probability (e.g. duplications) and focus on the following three mechanisms through which sequences vary:

1. A (point) mutation, or substitution, occurs when some amino acid in a sequence changes to some other amino acid during the course of evolution.
2. A deletion occurs when an amino acid is deleted from a sequence during the course of evolution.
3. A insertion occurs when an amino acid is added to a sequence during the course of evolution.

There are many possible sequences of events that could change one genome into another, and we wish to establish an optimality criterion that allows us to pick the “best” series of events describing changes between sequences. We choose to invoke Occam’s razor and select a maximum parsimony method as our optimality criterion^[7]: we wish to minimize the number of events used to explain the differences between two sequences. In practice, it is found that point mutations are more likely to occur than insertions and deletions, and certain mutations are more likely than others. We now develop an algorithmic framework where these concepts can be incorporated explicitly by introducing parameters that take into account the “biological cost” of each of these changes.

1.3.1Homology¶

While here my focus is on protein structure, I also want to point out that in bioinformatics one of the key goals of sequence alignment is to identify homologous sequences (e.g., genes) in a genome. Two sequences that are homologous are evolutionarily related, specifically by descent from a common ancestor. The two primary types of homologs are orthologous and paralogous. While they are outside the scope of these notes, I note here that other forms of homology exist (e.g., xenologs, which are genes in different species that have arisen through horizontal gene transfer^[8] rather than through vertical inheritance from a common ancestor, and often perform similar functions in their respective organisms despite their distinct evolutionary origins).

• Orthologs arise from speciation events, leading to two organisms with a copy of the same gene. For example, when a single species A speciates into two species B and C, there are genes in species B and C that descend from a common gene in species A, and these genes in B and C are orthologous (the genes continue to evolve independent of each other, but still perform the same relative function).
• Paralogs arise from duplication events within a species. For example, when a gene duplication occurs in some species A, the species has an original gene B and a gene copy B’, and the genes B and B’ are paralogus.

Generally, orthologous sequences between two species will be more closely related to each other than paralogous sequences. This occurs because orthologous typically (although not always) preserve function over time, whereas paralogous often change over time, for example by specializing a gene’s (sub)function or by evolving a new function. As a result, determining orthologous sequences is generally more important than identifying paralogous sequences when gauging evolutionary relatedness.

1.3.2Definitions¶

Following Kellis (2016), we introduce a simplified version of the alignment problem and make it more complex by adding the required ingredients. First, some definitions:

Sequence

In mathematics (and combinatorics) a sequence is a series of characters taken from an alphabet Σ. For what concerns us, DNA molecules are sequences over the alphabet $\lbrace A, C, G, T\rbrace$, RNA sequences over the alphabet $\lbrace A, C, G, U\rbrace$, and proteins are sequences over the alphabet $\lbrace A, R, N, D, C, Q, E, G, H, I, L, K, M, F, P, S, T, W, Y, V \rbrace$.

Substring

A consecutive part of a sequence. A sequence of length $N$ has $N$ substrings of length 1, $N(N-1)/2$ substrings of length 2, $\ldots$, 2 substrings of length $(N - 1)$. Given the sequence “SUPERCALIFRAGILISTICHESPIRALIDOSO”, “FRAGILI” is a substring, but “SUPERDOSO” is not.

Subsequence

A set of characters of a sequence. The characters do not have to be consecutive, but they should be ordered like in the original sequence. Given the sequence “SUPERCALIFRAGILISTICHESPIRALIDOSO”, “FRAGILI” and “SUPERDOSO” are two subsequences, but “SOSUPER” is not. Formally, given a sequence $S = (s_1, s_2, \ldots, s_N)$, $K = (z_1, z_2, \ldots, s_n)$, with $n \leq N$, is a subsequence of $S$ if there exists a strictly increasing sequence $i_1 < i_2 < \ldots < i_n$ of indices of $S$ such that $S_{i_j} = K_j$ for each $j: 1 \leq j \leq n$.

Sequence alignment

An alignment of two sequences $S$ and $T$, defined over the same alphabet Σ, is a $2 \times L$ table containing either a character from Σ, or the gap symbol “-”. In the table there is no column in which both characters are “-”, and the non-gap characters in the first line give back sequence $S$, while the non-gap characters in the second line give back sequence $T$.

1.3.3Problem formulation¶

We warm up by considering the following simple problem: given two sequences, $S$ and $T$, what is their longest common substring? We take the two following DNA fragments as examples: ACGTCATCA and TAGTGTCA. The problem can be solved by aligning the first characters of the two sequences, and then shifting one, say $T$, by an integer amount $i \in \mathcal{N}$. By computing the number of matching characters for every $i$ we can find the optimal alignment, which in this case is given by $i = -2$:

--ACGTCATCA
--xx||||---
TAGTGTCA---

The first and third rows of the diagram above describe the alignment itself, as introduced earlier, showing that gaps are characters that align to a space, while the central row shows the character matches: here x and | stand for mismatches and matches, respectively. In this example, the longest common substring is GTCA. The algorithm I just described has a complexity of $\mathcal{O}(n^2)$, where $n$ is the length of the shortest sequence.

A more complicated problem is to find the longest common subsequence (LCS) between two sequences. In the language we are introducing, this means allowing internal gaps in the alignment. The LCS between the $S$ and $T$ sequences defined earlier is

-ACGTCATCA
-|-||x-|||
TA-GTG-TCA

In this case the LCS is AGTTCA, which is longer than the longest common substring.

The LCS problem as formulated here is a particular case of the full sequence-alignment problem. In order to generalise the problem, we introduce a cost function that makes it possible to recast it in terms of an optimisation problem. Given two sequences $S$ and $T$, we define $\delta(S, T)$ as the “biological cost” of turning $S$ into $T$. For the case just considered, I implicitly used as a cost function the Levenshtein (or edit) distance, which is defined as the minimum number of single-character edits required to change one string into another. In other words, the problem has been formulated by implicitly considering that all possible operations (mutations, insertions and deletions) are equally likely. In a more generalised formulation, each edit operation is associated to a specific cost (penalty) that should reflect its biological occurrence, as discussed later.

What about gaps? In general, the cost of creating a gap depends on many variables. There are varying degrees of approximations that can be taken in order to simplify the problem. At the zero-th order each gap can be assumed to have a fixed cost, as we implicitly did above. This is called the “linear gap penalty”. An improvement can be done by considering that, biologically, the cost of creating a gap is more expensive than the cost of extending an already created gap. This is taken into account by the “affine gap penalty” method, whereby there is a large initial cost for opening a gap, and then a small incremental cost for each gap extension. There are other (more complicated and more context-specific) methods that can be considered and that we will not analyse further here, such as the general gap penalty, which allows for any cost function, and the frame-aware gap penalty, which is applied to DNA sequences and tailors the cost function to take into account disruptions to the coding frame^[9].

1.3.4The Needleman-Wunsch algorithm¶

Once we allow for gaps, the enumeration of all possible alignments (as done for the longest common substring method) becomes unfeasible. Indeed, the number of non-boring alignments, i.e. alignments where gaps are always paired with characters, in a sequence of size $N$ containing $G$ gaps (where $N \approx G$) can be estimated as

where we used the second-order Stirling’s approximation, $\log(N!) \approx N \log(N) - N + \frac{1}{2}\log (2 \pi N)$. Note that this number grows very fast: for $N = 30$ it is already larger than 10¹⁷. Considering that we also have to compute the score of each alignment, it is clear that the problem is untractacle with a brute-force method. Here is where dynamic programming enters the field.

Suppose we have an optimal alignment for two sequences $S$, of length $N$, and $T$, of length $M$, in which $S_i$ matches $T_j$. The alignment can be conceptually split into three parts:

1. Left subalignment: The alignment of the subsequences $(S_1 , \ldots , S_{i−1})$ and $(T_1 , \ldots , T_{j −1} )$.
2. Middle match/mismatch: The alignment of $S_i$ with $T_j$ (this could be a match or a mismatch).
3. Right subalignment: The alignment of the subsequences $(S_{i+1}, \ldots, S_N)$ and $(T_{j+1} , \ldots, T_M)$.

The overall alignment score is the sum of the scores from these three parts: the score of the left subalignment, the score for aligning $S_i$ with $T_j$, and the score of the right subalignment. Therefore, if the overall alignment is optimal, both the left and right subalignments must themselves be optimal. This follows from a cut-and-paste argument: imagine that one of the two subalignment was not optimal. This would mean there exists another alignment of these subsequences with a higher score. If such a better alignment for the subsequences existed, we could “cut” the suboptimal subalignment from the original alignment and “paste” in this better alignment. This would result in a new alignment for $S$ and $T$ with a higher total score than the supposed “optimal” alignment. This is a contradiction because the original alignment was assumed to be optimal. Since the same argument applies to both subalignments, the overall alignment’s optimality depends on the optimality of both its left and right subalignments. Of course, this is true only if the scores are additive, so that the score of the overall alignment is the sum of the scores of the alignments of the subsequences. The implicit assumption is that the sub-problems of computing the optimal scoring alignments of the subsequences are independent.

Let $F_{i,j}$ be the score of the optimal alignment of $(S_1 , \ldots , S_i)$ and $(T_1 , \ldots , T_j)$. Since $i \in [0, N]$ and $j \in [0, M]$, the matrix $\hat F$ storing the solutions (i.e. optimal scores) of the subproblems has a size of $(M + 1) \times (N + 1)$.

We can compute the optimal solution for a subproblem by making a locally optimal choice based on the results from the smaller subproblems. Thus, we need to establish a recursive function that shows how the solution to a given problem depends on its subproblems and can be used to fill up the matrix $\hat F$. At each iteration we consider the four possibilities (insert, delete, substitute, match), and evaluate each of them based on the results we have computed for smaller subproblems.

We start by considering the linear gap penalty model, and define $d$, with $d < 0$, as the cost of a gap. We are now equipped to set the values of the elements of the matrix. Let’s consider the first row: the value of $F_{0,j}$ is the cost of aligning a sequence of length 0 (taken from $S$) to a sequence of length $j$ (taken from $T$), which can be obtained only by adding $j$ gaps, yielding $F_{0,j} = jd$. Likewise, for the first column we have $F_{i,0} = id$. Then, we traverse the matrix element by element. Let’s consider a generic element $F_{ij}$: this is the cost of aligning the first $i$ characters of $S$ to the first $j$ characters of $T$. There are three ways we can obtain this alignment:

• Align $S_i$ to a gap (i.e. a gap is added to $T$): the total cost is $F_{i-1, j} + d$;
• Align $T_j$ to a gap (i.e. a gap is added to $S$): the total cost is $F_{i, j - 1} + d$;
• $S_i$ and $T_j$ are matched: the total cost is $F_{i - 1, j - 1} + s(S_i, T_j)$, where $s(x, y)$ is the cost of (mis)matching $x$ and $y$.

Leveraging the cut-and-paste argument sketched above, the optimal cost is given by the largest of the three values. Formally,

This is a recursive relation: computing the value of any $F_{i,j}$ requires the knowledge of the values of its left, top, and top-left neighbours. Therefore, the fill-in phase amounts to traversing the table in row or column major order, or even diagonally from the top left cell to the bottom right cell. After traversing the matrix, the optimal score for the alignment is given by the bottom-right element, $F_{MN}$. In order to obtain the actual alignment we have to traceback through the choices made during the fill-in phase. It is helpful to maintain a pointer for each cell while filling up the table that shows which choice was made to get the score for that cell. Then the pointers can be followed backwards to reconstruct the optimal alignment.

The complexity analysis of this algorithm is straightforward. Each update takes $\mathcal{O}(1)$ time, and since there are $MN$ elements in the matrix $\hat F$, the total running time is $\mathcal{O}(MN)$. Similarly, the total storage space is $\mathcal{O}(MN)$. For the more general case where the update rule is more complicated, the running time may be more expensive. For instance, if the update rule requires testing all sizes of gaps (e.g. the cost of a gap is not linear), then the running time would be $\mathcal{O}(MN(M + N)$).

1.3.5Local alignment: the Smith-Waterman algorithm¶

The Needleman-Wunsch algorithm find the best possible alignment across the entire length of two sequences. It tries to align every character from the start to the end of the sequences, which means both sequences are considered in their entirety. This is called “global alignment”, and it is most useful when the sequences being compared are of similar length and are expected to be homologous across their entire length.

Local alignment, on the other hand, focuses on finding the best alignment within a subset of the sequences. It identifies regions of similarity between the two sequences and aligns only those regions, ignoring the parts of the sequences that do not match well. Local alignment is particularly useful when comparing sequences that may only share a segment of similarity, such as when comparing domains within proteins, detecting conserved motifs, or identifying homologous regions in sequences that may not be overall similar^[10], which is why is very useful for the prediction of the 3D structure of proteins (or protein subdomains).

The most used method for local alignment is the Smith-Waterman algorithm, which is a modification of the Needleman-Wunsch algorithm. The key difference between the two lies in how the scoring matrices are constructed and scored. In Needleman-Wunsch, every cell in the matrix is filled to reflect the best global alignment, with the final alignment score found in the bottom-right corner of the matrix. Smith-Waterman, on the other hand, sets any negative scores to zero, which allows the algorithm to “reset” when the alignment quality dips. The highest score in the matrix indicates the end of the best local alignment, which is then traced back to a zero to identify the optimal aligned subsequence. This approach ensures that only the most relevant, highest-scoring local alignments are highlighted. Here is how the Smith-Waterman algorithm looks like in practice:

1. Initialisation: since a local alignment can start anywhere, the first row and column in the matrix are set to zeros, i.e. $F_{0,j} = jd$, $F_{i,0} = id$.
2. Iteration $\forall (i, j)$: this step is modified so that the score is never allowed to become negative but it is reset to zero. This is done by slightly modifying Eq. (37) as follows:
   $$F_{i,j} = \max \begin{cases} 0\\ F_{i - 1, j} + d \\ F_{i, j - 1} + d \\ F_{i - 1, j - 1} + s(S_i, T_j). \end{cases}$$

   (38)
3. Trace-back: starts from the position of the maximal number in the table and proceeds until a zero is encountered.

A visual comparison between the tables generated by global and local alignment algorithms is shown in Figure 12.

1.3.6Affine gap penalty¶

For both the global and local alignment algorithms introduced we have used a linear gap penalty: the cost of adding a gap is constant, regardless of the nature of the aligned character, or of the length of the gap. From the biological point of view, this means that an indel of any length is considered as the result of independent insertions or deletions. However, in reality long indels can form in single evolutionary steps, and in these cases the linear gap model overestimates their cost. To overcome this issue, more complex gap penalty functions have been introduced. As mentioned before, a generic gap penalty function would result in an algorithmic complexity worse than $\mathcal{O}(N^2)$ (where for simplicity I’m considering two sequences of the same length). Let’s see why. Any gap penalty function can be implemented in the Needleman-Wunsch or Smith-Waterman algorithms by changing the recursive rule, which can be done trivially for element $F_{i,j}$ by evaluating terms such as $\max_{0 \leq k \leq i} \lbrace F_{k,j} + d_{i - k} \rbrace$ and $\max_{0 \leq k \leq j} \lbrace F_{i,k} + d_{j - k} \rbrace$, where $d_x$ is the cost of a gap of size $x$. However, this means that the update of every cell of the dynamic programming matrix would take $\mathcal{O}(N)$ instead of $\mathcal{O}(1)$, bringing the algorithmic complexity up to $\mathcal{O}(N^3)$ (for $N = M$ sequences).

However, the computational cost can be mitigated by using particular gap penalty functions. Here I will present the most common variant, which is known as the affine gap penalty. In this model, the cost of a gap of size $k$ is

where $o$ and $e$ are the opening and extension penalties, respectively, and $o > e$. To incorporate affine-gap penalties into the Needleman-Wunsch or Smith-Waterman algorithms in an efficient manner, the primary adjustment involves tracking whether consecutive gaps occur in the alignment. This requires the alignment process to be split into three distinct cases: insertions, deletions, and matches/mismatches. Instead of using a single dynamic programming table as in the linear-gap penalty approach, three separate tables are used: $\hat I$ for insertions, $\hat D$ for deletions, and $\hat F$ for the overall score.

In this setup, each entry in the $\hat I$ table, denoted by $i_{i,j}$, stores the best alignment score when the last column includes an insertion (i.e. a gap in the first sequence), and each entry in the $\hat D$ table, $d_{i,j}$, captures the best score for alignments ending in a deletion (i.e. a gap in the second sequence). As before, $F_{i,j}$ records the overall score. The recursive update rules of the three tables are

Note that with this algorithm each cell update is $\mathcal{O}(1)$, and therefore the overall complexity remains the same ($\mathcal{O}(NM)$ or $\mathcal{O}(N^2)$ for same-length sequences).

1.3.7Substitution matrices¶

How is the substitution matrix $s(x, y)$ determined? One possibility is to leverage what we know about the biological processes that underlie mutations. For instance, in DNA the biological reasoning behind the scoring decision can be linked to the probabilities of bases being transcribed incorrectly during polymerization. We already know that of the four nucleotide bases, A and G are purines, while C and T are pyrimidines. Thus, DNA polymerase is much more likely to confuse two purines or two pyrimidines since they are similar in structure. As a result, a simple improvement over the uniform cost function used above is the following scoring matrix for matches and mismatches:

Here a mismatch between like-nucleotides (e.g. A and G) is less expensive than one between unlike nucleotides (e.g. T and C).

However, we can be more quantitative by taking a probabilistic approach. Here I will describe how the widely-used BLOSUM matrices are built (see also the original paper). We assume that the alignment score reflects the probability that two similar sequences are homologous. Thus, given two sequences $S$ and $T$ of the same length^[11] for which we have an ungapped alingment, we define two hypotheses:

1. The alignment between the two sequences is due to chance and the sequences are, in fact, unrelated.
2. The alignment is due to common ancestry and the sequences are actually related.

Under the rather strict assumption (“biologically dubious, but mathematically convenient”, as aptly put in Eddy (2004)) that pairs of aligned residues are statistically independent of each other, so that the likelihoods associated to the two hypotheses, $P(S, T|R)$ and $P(S, T|U)$, can be expressed as products of probabilities, viz.

where $p_{xy}$ is the likelihood that residues $x$ and $y$ are aligned because correlated, while the product $q_x q_y$ is the likelihood that the two residues are there by chance: their occurrence is unrelated and therefore the likelihood factorises in two terms that account for the average probability of observing those two residues in any protein.

In order to compare two hypotheses, a good score is given by the logarithm of the ratio of their likelihoods (the so-called log-odds score). Therefore, for our case the alignment score is

so that, thanks to the properties of logarithms and probabilities, the overall alignment score is the sum of the log-scores of the individual residue pairs. Considering 20 amino acids, there are 400 such log-scores, which form a 20x20 substitution (score) matrix. Each entry of the matrix takes the form

where λ is an overall scaling factor that is used to obtain values that can be rounded off to integers. Note that the previous definition can be formally substantiated (see e.g. Karlin & Altschul (1990)).

If $p_{xy} > q_x q_y$ (and therefore $s(x, y)$ is positive), the substitution $x \to y$, and therefore their alignment in homologous sequences, occurs more frequently than would be expected by chance, suggesting that this substitution is favored in evolution. These substitutions are called “conservative”. As noted in Eddy (2004), this definition is purely statistical and has nothing directly to do with amino acid structure or biochemistry. Likewise, substitutions with $p_{xy} < q_x q_y$, and therefore negative values of $s(x, y)$, are termed “nonconservative”.

The procedure applied to compute the $p_{xy}$ and $q_x$ values starts with a collection of protein sequences.

1. The sequences are grouped into families based on their evolutionary relatedness. A common source for these families is the BLOCKS database, which contains aligned, ungapped regions of proteins that are highly conserved across members of the family.
2. Within these protein families, highly conserved regions (blocks) are identified. These regions are short segments of amino acid sequences that are believed to be important for the protein’s function and are conserved across different species.
3. Once blocks are identified, the sequences within each block are clustered, choosing a threshold value $V$: two sequences that share more than this percentage of identical amino acids identity are clustered together, and only one representative from each cluster is used to avoid over-representation of very similar sequences.
4. Within each block, the actual counts of amino acid pairs are made.
   • The background frequencies $q_x$ are estimated by computing the overall frequency of each amino acid $x$ in the sequences, not considering any substitutions.
   • To evaluate the $p_{xy}$ terms, for each position in the aligned block, we count how many times one amino acid (say, Alanine) is substituted by another amino acid (say, Glycine) in the aligned positions across the different sequences. For instance, if there is an alignment of 5 sequences, for each single position we make $\binom{5}{2} = 10$ comparisons.

Of course, the final scores depend on $V$: if $V$ is large, protein blocks that are still rather similar will be considered belonging to different clusters, and therefore compared to each other to derive the scores; the resulting matrix will be more sensitive in identifying closely related sequences but less effective for more distantly related sequences (and vice versa for small values of $V$). Common values for $V$ are $45\%$, $62\%$, and $80\%$, which yield the matrices BLOSUM45, BLOSUM62, and BLOSUM80, with BLOSUM62 being the de-facto standard (and default) one.

The BLOSUM62 matrix is shown in Figure 13. First of all, note that substitution matrices are symmetric, because the biological process of amino acid substitution is symmetric: as empirically observed, there is no preferred direction when one amino acid replaces another. Secondly, it is evident that conservative substitutions tend to be those between amino acids that are similar, as made evident by grouping the amino acids according to the classification introduced by Margaret Dayhoff.

Why the values of the matrix diagonal, representing the scores of “substituting” one amino acid with itself, are all different?

1.3.8BLAST¶

The sheer volume of sequence data generated by modern-day high-throughput sequencing technologies presents a significant challenge. Databases now contain millions of nucleotide and protein sequences, each potentially spanning thousands of characters. When comparing a new sequence against these massive databases, traditional pairwise alignment methods like the ones we just discussed become computationally demanding. Indeed, performing a global or even local alignment between a query sequence and every sequence in a large database can require huge computational resources and time, especially as the number of sequences and their lengths continue to grow exponentially.

The need for a more efficient method resulted in the most cited paper of the 1990s, where BLAST (Basic Local Alignment Search Tool), now a critical tool in bioinformatics, was introduced. BLAST operates by finding regions of local similarity between sequences, which is more computationally feasible and faster than aligning entire sequences globally. The algorithm requires a query sequence and a target database. First of all, the query sequence is broken down into smaller fragments (called words or $W$-mers). For each $W$-mer, a list of similar words is generated, and only those with a similarity measure that is higher than a threshold $T$ are retained add added to the final $W$-mer list. The similarity is evaluated by using a substitution matrix (BLOSUM62 is a common choice). Then, the target database is searched for matches to these words, extending the matches in both directions to find the best local alignments. BLAST assigns scores to these alignments based on the degree of similarity via a Smith-Waterman algorithm, with higher scores indicating closer matches.

The BLAST algorithm can be broken down into the following steps^[12]

1. Optional: remove low-complexity region or sequence repeats in the query sequence. Here “Low-complexity region” means a region of a sequence composed of few kinds of elements. These regions might give high scores that confuse the program to find the actual significant sequences in the database, so they should be filtered out. The regions will be marked with an X (protein sequences) or N (nucleic acid sequences) and then be ignored by the BLAST program. To filter out the low-complexity regions, the SEG program is used for protein sequences and the program DUST is used for DNA sequences. On the other hand, the program XNU is used to mask off the tandem repeats in protein sequences.
2. Make a $W$-letter word list of the query sequence. The words of length $W$ ($W$-mers) in the query sequence are listed “sequentially”, until the last letter of the query sequence is included. The method is illustrated in TODO. $W$ is usually 3 and 11 for a protein and a DNA sequence, respectively.
3. List the possible matching words. A substitution matrix (e.g. BLOSUM62) is used to match the words listed in step 2 with all the $20^W$ $W$-mers. For example, the score obtained by comparing PQG with PEG and PQA is respectively 15 and 12 with the BLOSUM62 matrix^[13]. After that, a neighborhood word score threshold $T$ is used to reduce the number of possible matching words. The words whose scores are greater than the threshold T will remain in the possible matching words list, while those with lower scores will be discarded. For example, if $T = 13$ PEG is kept, but PQA is abandoned.
4. Organize the remaining high-scoring words into an efficient search tree. This allows the program to rapidly compare the high-scoring words to the database sequences.
5. Repeat step 3 to 4 for each $W$-mer in the query sequence.
6. Scan the database sequences for exact matches with the remaining high-scoring words. The BLAST program scans the database sequences for the remaining high-scoring words, such as PEG. If an exact match is found, this match is used to seed a possible ungapped alignment between the query and database sequences.
7. Extend the exact matches to high-scoring segment pair (HSP). The original version of BLAST stretches a longer alignment between the query and the database sequence in the left and right directions, from the position where the exact match occurred. The extension does not stop until the accumulated total score of the HSP begins to decrease. To save more time, a newer version of BLAST, called BLAST2 or gapped BLAST, has been developed. BLAST2 adopts a lower neighborhood word score threshold to maintain the same level of sensitivity for detecting sequence similarity. Therefore, the list of possible matching words list in step 3 becomes longer. Next, exact matched regions that are within distance $A$ from each other on the same diagonal are joined as a longer new region. Finally, the new regions are then extended by the same method as in the original version of BLAST, and the scores for each HSP of the extended regions are created by using a substitution matrix as before.
8. List all of the HSPs in the database whose score is high enough to be considered. All the HSPs whose scores are greater than the empirically determined cutoff score $S$ are listed. By examining the distribution of the alignment scores modeled by comparing random sequences, a cutoff score $S$ can be determined such that its value is large enough to guarantee the significance of the remaining HSPs.
9. Evaluate the significance of the HSP score. Karlin & Altschul (1990) showed that the distribution of Smith-Waterman local alignment scores between two random sequences is
   $$p\left( S\ge x \right) =1-\exp \left( -KMN e^{-\lambda x } \right),$$

   (44)
   where $M$ and $N$ are the length of the query and database sequences^[14], and the statistical parameters λ and $K$ depend upon the substitution matrix, gap penalties, and sequence composition (the letter frequencies) and are estimated by fitting the distribution of the ungapped local alignment scores of the query sequence and of a lot of (globally or locally) shuffled versions of a database sequence. Note that the validity of this distribution, known as the Gumbel extreme value distribution (EVD), has not been proven for local alignments containing gaps yet, but there is strong evidence that it works also for those cases. The expect score $E$ of a database match is the number of times that an unrelated database sequence would obtain a score $S$ higher than $x$ by chance. The expectation $E$ obtained in a search for a database of total length $N$
   $$E = K M N e^{-\lambda S}$$

   (45)
   This expectation or expect value $E$ (often called $E$-score, $E$-value or $e$-value) assessing the significance of the HSP score for ungapped local alignment is reported in the BLAST results. The relation above is different if individual HSPs are combined, such as when producing gapped alignments (described below), due to the variation of the statistical parameters.
10. Make two or more HSP regions into a longer alignment. Sometimes, two or more HSP regions in one database sequence can be made into a longer alignment. This provides additional evidence of the relation between the query and database sequence. There are two methods, the Poisson method and the sum-of-scores method, to compare the significance of the newly combined HSP regions. Suppose that there are two combined HSP regions with the pairs of scores $(65, 40)$ and $(52, 45)$, respectively. The Poisson method gives more significance to the set with the maximal lower score $(45 > 40)$. However, the sum-of-scores method prefers the first set, because $65+40=105$ is greater than $52+45 = 97$. The original BLAST uses the Poisson method; BLAST2 uses the sum-of scores method.
11. Show the gapped Smith-Waterman local alignments of the query and each of the matched database sequences. The original BLAST algorithm only generates ungapped alignments including the initially found HSPs individually, even when there is more than one HSP found in one database sequence. By contrast, BLAST2 produces a single alignment with gaps that can include all of the initially found HSP regions. Note that the computation of the score and its corresponding $E$-value involves use of adequate gap penalties.

The $E$-value is the single most important parameter to rate the quality of the alignments reported by BLAST. For instance, if $E = 10$ for a particular alignment with score $S$, it means that there are 10 alignments with score $\geq S$ that can happen by chance between any query sequence and the database used for the search. Therefore, this particular alignment is most likely not very significant. By contrast, values much smaller than 1 (e.g 10^-3 or even smaller) are likely to signal that the sequences are homologous. On most webserver, it is possible to filter out all matches that have an $E$-value larger than some threshold. For instance, on the NCBI webserver, the “expect threshold” defaults to 0.05.

Note that the sequence database used to search for matches is preprocessed first, which increases further the overall computational efficiency. With BLAST, the algorithmic complexity of searching the database for a sequence of length $N$ is only $\mathcal{O}(N)$. An online webserver^[15] to run BLAST can be found here.

1.3.9Multiple sequence alignment¶

Pairwise sequence alignment suffers from some issues that are intrinsic to it, and becomes glaring when trying to align sequences that are distantly related. In particular, all pairwise methods depend in some way or another on a number of parameters (scoring matrix, gap penalties, etc.), and it is hard to tell what is the “best” alignment if the method finds multiple alignments with the same score. Moreover, a pairwise alignment is not necessarily informative about the evolutionary relationship (and therefore about possibly conserved amino acids or whole motifs) of the sequences that are compared.

In order to overcome these issues, a number of multiple sequence alignment (MSA) methods have been developed. MSA extends the concept of pairwise protein alignment to simultaneously align three or more sequences, providing a broader view of evolutionary relationships, structural conservation, and functional regions among a group of proteins. Aligning multiple sequences make it possible to identify conserved amino acids that may be critical for protein function or stability and therefore are biologically significant. This, in turn, can help to predict structural features, functional motifs, and evolutionary patterns across species, providing a fundamental tool to build phylogenetic trees or inform the classification of proteins into families.

In principle, MSA can be carried out with the dynamic programming algorithms we already introduced. For instance, the recursive rule of the Needleman-Wunsch algorithm for globally aligning three sequences $S$, $T$, and $U$, which extends Eq. (37), is

where the dynamic programming “table” is now three-dimensional, and $s = s(a, b, c)$ is the cost function of aligning $a$, $b$, and $c$, which can be either residues or gaps. It is straightforward to see that the dimensionality of the table grows linearly with the number of sequences $M$, and therefore the algorithmic complexity is $\mathcal{O}(N^M)$. The exponential dependence on $M$ makes this approach practically unfeasible when working with real-world examples.

Unfortunately, better-performing exact methods, i.e. methods that provide the global optimal solution by construction, do not exist (yet). As a result, we have to rely on heuristic methods, which only find local minima. One commonly used approach for multiple sequence alignment is called progressive multiple alignment. Here the prerequisite is that we need to know the evolutionary tree, often called the guide tree, that connects the sequences we wish to align, which is usually built by using some low-resolution similarity measure, often based on (global) pairwise alignment. Then, we start by aligning the two most closely related sequences in a pairwise fashion, creating what is known as the seed alignment. Next, we align the third closest sequence to this seed, replacing the previous alignment with the new one. This process continues, sequentially adding and aligning each sequence based on its proximity in the tree, until we reach the final alignment. Note that this is done using a “once a gap, always a gap” rule: gaps in an alignment are not modified during subsequence alignments. These methods have a computational complexity of $\mathcal{O}(M^2)$, which makes it possible to align thousands of sequences.

One of the most commonly used tools to perform MSAs is Clustal Omega, which is a command-line tool, but it is also available as a webserver. Clustal Omega builds the guide tree using an efficient algorithm (adapted from Blackshields et al. (2010)) which has an algorithmic complexity of $\mathcal{O}(M \log M)$, making it possible to generate multiple alignments of hundreds of thousands sequences. An example of MSA is shown in Figure 15.

1.4Threading¶

If pair alignment tools find that a given query sequence is found to share more than 30% of its sequence with another, then it is common to think that a reasonable model for that sequence can be built. By contrast, an alignment yielding a similarity of 20% - 25% could be purely coincidental. In reality, things are more complicated, as it has been shown that proteins with rather high sequence identity could be very different from a structural point of view (Brenner et al. (1998)). This and other results show that a single quantity, such as sequence identity, is not enough to determine the 3D similarity between two proteins, and more numbers (such as the length of the chain or of the well-aligned regions) are required to build reliable models. For instance, it makes sense that for a short 50-residue protein a 40% sequence identity would be required to generate a good match, while 25% may be enough for 250 residues. However, these numbers have a purely statistical value.

A better approach is to go beyond pairwise alignments by using sequence database searching programs such as BLAST which, as we have seen, provide E-values or similar quantities that estimate the reliability of a sequence match by looking at it in the context of the whole library of sequence scores. In addition, more sophisticated BLAST versions (such as PSI-BLAST) make it possible to obtain good matches with less than 20% sequence identity.

Thanks to these tools, and to the ever-growing number of sequences stored in databases, “simple” database searches are often enough to build a model of the query sequence out of a reliable homolog of known structure. However, if no matches are found, or if an independent confirmation is required, we need methods that do not rely on sequence homology. In this case we recast the problem of protein structure prediction as a problem of choice of the 3D structure best fitting the given sequence among many other possible folds. However, what can be the source of “possible” structures?

While an a priori classification is sometimes possible (see, e.g., Murzin & Finkelstein (1988)) or Ptitsyn et al. (1985)), a more practical answer is the Protein Data Bank (PDB) where all the solved and publicly available 3D structures are collected. However, using structures stored in the PDB turns a “prediction” problem into a “recognition” one: a fold cannot be recognized if the PDB does not contain an already solved analog. This limits the power of recognition. Nevertheless, there is an important advantage associated to this procedure: if the protein fold is recognized among the PDB-stored structures, one can hope to recognize also the most interesting feature of the protein, namely its function—by analogy with that of an already studied protein.

Certainly, not all protein folding patterns have been collected in PDB yet; however, it most likely already includes the majority of all the folding patterns existing in nature. This hope, substantiated by Chothia (1992), is based on the fact that the folds found in newly solved protein structures turn out to be similar to already known folds more and more frequently. Extrapolation shows that perhaps about 1500–2000 folding patterns of protein domains exist in genomes, and we currently know more than half of them (including the majority of the most common folds).

To recognize the fold of a chain having no visible homology with already solved proteins, one can use various superpositions of the chain in question onto all examined (taken from an a priori classification or from PDB) 3D folds in search of the lowest-energy chain-with-fold alignment, as sketched in Figure 16. This is called the threading method. When a chain is aligned with the given fold, it is threaded onto the fold’s backbone until its energy (or rather, free energy) is minimized, including both local interactions and interactions between remote chain regions. The threading alignment allows “gaps” in the chain and in the fold’s backbone (the latter are often allowed for irregular backbone regions only). Many different algorithms have been proposed for finding the correct threading of a sequence onto a structure, though many make use of dynamic programming in some form. Indeed, in principle, threading is similar to a homology search; the difference is that only sequences are aligned in a homology search, while threading aligns a “new” sequence with “old” folds.

Being physics-inspired, threading can be a powerful tool for structural biologists and biophysicists. It uses mostly statistics-derived pseudopotentials rather than actual energies, which are based on the contact statistics between amino acids as found in known protein structures^[16].

As with any other structure-prediction model, there are some issues with threading, the main ones being:

1. The conformations of the gapped regions remain unknown, together with all interactions in these regions.
2. Even the conformations of the aligned regions and their interactions are known only approximately, since the alignment does not include side-chain conformations (which may differ in “new” and “old” structures). Estimates show that threading takes into account, at best, half of the interactions operating in the protein chain, while the other half remain unknown to us. Thus, again, the protein structure is to be judged from only a part of the interactions occurring in this structure. Therefore, this can only be a probabilistic judgment.
3. The number of possible alignments is enormous (remember Levinthal’s paradox?), which makes it very hard to sort out all possible threading alignments and single out the best one (or ones). Here there are many methods that can be employed, one of which being the statistical mechanics of one-dimensional systems, and the related dynamic programming techniques.

As an example, Finkelstein & Ptitsyn (2016) shows the structure prediction done for the replication terminating protein (rtp) by threading. Having threaded the rtp sequence onto all PDB-stored folds, Flöckner et al. (1995) used threading with statistics-based pseudopotentials to show that the rtp fold must be similar to that of H5 histone. This a priori recognition, which is shown in Figure 17, turned out to be correct.

However, it also turned out that the alignment provided by threading deviates from the true alignment obtained from superposed 3D structures of rtp and H5 histone. On the one hand, this shows that even a rather inaccurate picture of residue-to-residue contacts can lead to an approximately correct structure prediction. On the other hand, this shows once again that all the mentioned flaws (insufficiently precise interaction potentials, uncertainty in conformations of nonaligned regions, of side chains, etc.) make it possible to single out only a more-or-less narrow set of plausible folds rather than one unique correct fold. Indeed, the set of the “most plausible” folds can be singled out quite reliably, but it still remains unclear which of these is the best. The native structure is a member of the set of the plausible ones, it is more or less close to the most plausible (predicted) fold, but this is all that one can actually say even in the very best case.

Threading methods became a tool for a tentative recognition of protein folds from their sequences. The advantage of these methods is that they formulate a recipe: do this, this and this, and you will obtain a few plausible folds, one of which has a fairly high chance of being correct.

1.5AlphaFold¶

The CASP (Critical Assessment of protein Structure Prediction) contest is a biennial competition established to evaluate the performance of computational methods for protein structure prediction. Since its inception in 1994, CASP has become a benchmark for assessing progress in the field, gathering researchers from across the globe to tackle one of the most complex challenges in computational biology. Participants receive amino acid sequences of proteins for which the 3D structure remains experimentally unknown to them but has been solved by researchers outside the competition. They then attempt to predict these structures using their algorithms, which are later evaluated against the experimentally determined ones. Figure 18 shows the score of the winners of the last 8 contests. As can be seen in the figure, in 2020 there was a breakthrough that made it possible to predict the 3D structure of proteins with a precision close to the experimental one.

The breakthrough has been possible thanks to AlphaFold2 (which I will refer to simply as “AlphaFold” from now on), which is a deep learning model developed by DeepMind, presented in Jumper et al. (2021), designed to predict protein structures with remarkable accuracy. The work behind AlphaFold has been awarded half of the 2024 Nobel Prize in Chemistry. The model’s architecture is based on neural networks, and integrates several advanced techniques from machine learning and structural biology to predict the 3D structure of a protein out of its sequence (i.e., the folding problem). I will briefly describe the internal architecture of Alphafold, and then show how to use it (in a slightly improved version called ColabFold, presented in Mirdita et al. (2022)).

1.5.1Input and preprocessing¶

The overall architecture, as presented in the original paper, is shown in Figure 19. The first part, which handles the user input and is highlighted in red in Figure 19, is a preprocessing pipeline that can be carried out independently of the rest as done, for instance, in ColabFold. First of all, the AlphaFold system uses the input amino acid sequence to query several databases of protein sequences, and constructs a multiple sequence alignment in order to determine the parts of the sequence that are more likely to mutate. The underlying idea is that, if two amino acids (possibly far apart along the sequence) are in close spatial contact, mutations in one of them will be closely followed by mutations of the other, so that the overall 3D structure is preserved. To make an extreme (and somewhat unrealistic) example, suppose we have a protein where an amino acid with negative charge (say, glutamate) is spatially close to an amino acid with positive charge (say, lysine), although they are both far away in the sequence. Most likely, the resulting electrostatic interaction stabilises the structure of the protein. If for evolutionary reasons the first AA mutates into a positively charged amino acid, the second AA will be under evolutionary pressure to mutate into a negatively charged amino acid in order to preserve the electrostatic attraction and therefore the contribution to the stability of the folded protein. The MSA makes it possible to detect this sort of correlations.

AlphaFold also tries to identify proteins that may have a similar structure to the input (“templates”), and constructs an initial representation of the structure called the “pair representation” which is, in essence, a model of which amino acids are likely to be in contact with each other. Finding templates follows a completely different, but closely related principle: proteins mutate and evolve, but their structures tend to remain similar despite the changes. In Figure 20, for example, I display the structure of four different myoglobin proteins, corresponding to different organisms. You can appreciate that they all look pretty much the same, but if you were to look at the sequences, you would find enormous differences. The protein on the bottom right, for example, only has ~25% amino acids in common with the protein on the top left.