Computational methods in the study of self-entangled proteins: a critical appraisal

As mentioned before, knot theory defines the topological state of closed curves. By common experience we know that, while a knot on an open string can be undone by proper manipulation, this is impossible if the ends of the string are attached to form a loop (see figures 1(a) and (b)). Any possible deformation in space of this closed loop preserves its topological state. This suggests that the knotted state of a closed curve can be operationally defined by means of spatial deformations.

Zoom In Zoom Out Reset image size

Let us mathematically represent the knotted closed string of figure 1(b) as a closed curve embedded in the three-dimensional euclidean space $X\subset \mathbb{R}^{3}$, as shown in figure 1(c). In mathematics, all the possible continuous deformations of X in space are called homotopies (see e.g. [2] for further details). To define the topological state of X we need to restrict the class of transformations to those homotopies that prevent the curve from passing through itself, named isotopies. This is still not sufficient, since X has no thickness, any entanglement hosted by the curve can be continuously reduced to a single point, transforming the curve into an un-knotted loop, also called trivial knot. This leads to the definition of an ambient isotopy (AI), which deforms X through the continuous transformation of its embedding space, in this case $\mathbb{R}^{3}$. The action of an AI on X does not change its topological state, thus the AIs are the mathematical analogous of the string manipulation mentioned before.

The topological state of a curve is defined through AIs as an equivalence class, named knot type. Two curves that can be transformed into each other by AIs belong to the same knot type, namely they are topologically equivalent. This definition is related to the concept of knot complement, namely $K=S-N(X)$, where S is a compact region embedding the curve and $N(X)$ is a tubular region that indicates the neighboring space of X. Indeed, K defines the knot type, as equivalent knots have homeomorphic complements [73]. For example, all the curves that can be transformed into a circle belong to the trivial knot type, and are said to be unknotted. There exist infinitely many possible knot types, the known ones being classified in catalogs [74, 75].

From this general definition of knots, we shall restrict to those knot types that can describe a physical objects, which are named tame knots. Tame knots, as demonstrated in [2], can always be represented by closed, non-intersecting and finite polygonal curves, called polygonal knots. An example of a knot type excluded from this definition is the so-called wild knot, shown in figure 2, which features an infinitely recursive character and cannot represent a physical knot.

Zoom In Zoom Out Reset image size

Moreover, since most of the techniques reviewed in the following have been conceived for the topological analysis of polymers, it will be sometimes useful to directly refer to a bead-chain structure rather than to an abstract polygonal curve.

In the process of topological classification, a schematic visualization of knots, capable to underline differences between knot types, provides a significant help. A natural, compact way to visualize a curve embedded in $\mathbb{R}^{3}$ is to project it on a plane. The mapping of a three dimensional curve on a plane naturally generates singular points, in which more than one point of the curve is mapped. The relevant projection for knot identification is the so-called regular projection, in which the singularities are ordinary double points, also called crossing points (CP). In CPs two strands of the curve overlap with different tangents, as shown in figure 3. It can be demonstrated that physical knots can always be represented through a regular projection [2]. To preserve the topological information in the two-dimensional projection it is also necessary to assign the information of relative depth of the crossing branches in each double point. This is normally done as shown in figure 3, by interrupting the branch that underpasses the CP. The resulting representation is named knot diagram.

Zoom In Zoom Out Reset image size

Knot diagrams can be extremely complex, in particular if they represent a collapsed globular structure, such as that of proteins. However, a complex knot diagram can be transformed to reach a simpler representation, without changing its topological state. Intuitively, the transformations that allow one to reduce the complexity of a knot diagram are planar projections of AIs. These consist in all planar deformations of the diagram that do not affect the CPs, plus the three Reidemeister moves (RM) [77]. The latter ones act on the number of CPs of a diagram without affecting its topology (a pictorial representation of RMs is reported in figure 4). These transformations generate the class $\mathcal{D}(X)$ of topologically equivalent diagrams of X. One can thus use Reidemeister moves to reduce the number n of double points in a diagram down to the minimum number $\tilde{n}=\min\nolimits_{\mathcal{D}(X)} n$, named crossing number, obtaining the so-called minimal diagram. We underline that the minimal diagram is generally the simplest representation of a knot, but sometimes it hides specific knot properties (for further insights on this we refer to specific knot theory books, such as [2, 65]).

Zoom In Zoom Out Reset image size

The crossing number $\tilde{n}$ of a knot diagram is a topological invariant, namely a property that depends only on the topological state of the curve and not on a specific three-dimensional realization or planar projection. More precisely, $\tilde{n}$ is a weak invariant, as different knots can have the same crossing number. The known knots are classified using the notation $\tilde{n}_{i}$, in which i indexes different knots with the same crossing number. The trivial knot is classified as 0₁, while the simplest known knot, namely the trefoil knot, is referred to as 3₁. According to these conventions known knots are classified in tables such as the one reported in figure 5(a), including all knot types found in protein structures. A single curve can also form composite knots, that contain more than one knot connected as in figure 5(b). The components of a composite knot are called factor knots, while prime knots cannot be decomposed in two or more non-trivial factors (the trivial knot is always a factor of another knot). Knot tables, as that reported in figure 5(a) usually include only prime knots.

Zoom In Zoom Out Reset image size

A further topological specification is the chirality, or handedness, of a knot. A knot is chiral if there is no AI that can map it to its mirror image. The two mirror images of a chiral knot are named enantiomers. The simplest example of chiral knot is the 3₁, while the 4₁ is an example of a-chiral knot.

For decades mathematicians have searched for topological invariants which can distinguish between different knots, and at the same time be efficiently computed also for very complex embeddings of curves. As mentioned before, the crossing number, on which the knot classification is built, is the simplest topological invariant. However its calculation requires by definition the determination of the minimal knot diagram, and it is therefore a complex operation. Nowadays many invariants have been discovered, ranging from Gauss winding number [78], to the knot group [79, 80].

The most widely used invariants for identifying knots in polymers are the polynomial invariants, that are constructed defining combinatorial rules that apply to any knot diagram, not requiring the construction of the minimal diagram. These invariants include the Alexander polynomial [78], Jones polynomial [81], and HOMFLY polynomial [82] (from the names of its co-discoverers Hoste, Ocneanu, Millett, Freyd, Lickorish and Yetter). A detailed review of these invariants is beyond the purposes of present work, the interested reader can refer to [2, 80]. As an example on how these polynomials are constructed, we report the procedure for the calculation of the Alexander polynomial $\Delta(t)$ [66, 83]:

• (i)  
  First, an arbitrary orientation is assigned to the diagram. This orientation defines the sign of the crossings according to the convention represented in figure 6(b).
• (ii)  
  Following the orientation, a progressive number $x=1\ldots n$ is assigned to the crossing points of the diagram, starting from an arbitrary point.
• (iii)  
  Also the n arcs, the diagram branches separated by underpasses, are numbered progressively following the orientation.
• (iv)  
  An $n\times n$ matrix M is defined. Each row is associated to a crossing point, while each column is associated to an arc. The nonzero elements of the xth row correspond to the three arcs that meet in the xth crossing point. Let us define that in x the ith arc overpasses the two consecutive j  and k arcs. The xth row elements are then defined according to the following rules:

  • if the crossing is positive: $M(x,i)=1-t$, $M(x,j)=-1$ and $M(x,k)=t$,
  • if the crossing is negative: $M(x,i)=1-t$, $M(x,j)=t$ and $M(x,k)=-1$,
  • if i  =  k, or i  =  j , $M(x,j)=1$ and $M(x,k)=-1$.
• (v)  
  An arbitrary minor of M is selected to obtain a $(n-1)\times(n-1)$ matrix, the so-called Alexander matrix.
• (vi)  
  The determinant of the Alexander matrix is computed to obtain the Alexander polynomial $\tilde{\Delta}(t)$.

Zoom In Zoom Out Reset image size

The $\tilde{\Delta}(t)$ obtained in this way is not a topological invariant, but all the $\tilde{\Delta}(t)$'s corresponding to the same knot (computed from different diagrams) differ by a factor $\pm t^{m}$, with $m\in\mathbb{Z}$. Therefore one can obtain the so-called irreducible Alexander polynomial $\Delta(t)=\pm t^{m}\tilde{\Delta}(t)$, where the sign and $m\in\mathbb{Z}$ are chosen so that the lowest order term of $\Delta(t)$ is a positive constant. The irreducible Alexander polynomial is an actual topological invariant. In the following, as it is customary in the literature, we will drop the 'irreducible' and simply indicate it as Alexander polynomial. In figure 6 an example of the algorithm application is displayed.

It is clear from this example that the calculation of polynomial invariants depends on the number n of crossing points associated to the particular projection considered, thus making it desirable to simplify the knot diagram before the calculation, in order to operate on a smaller—if not the smallest possible—number of crossings. It can be shown that Jones polynomial calculation time scales as 2ⁿ, while Alexander polynomial calculation time scales as (n  −  1)³ [84]. For its scaling properties, and for its relatively simple implementation, the Alexander polynomial is largely used in the literature of protein topology (see e.g. [24, 85–87]). However, since Alexander polynomial cannot distinguish between two enantiomers, a more general invariant, such as the HOMFLY polynomial, has to be computed when the chirality information is required [87]. We also stress that polynomial invariants cannot always distinguish among different knots, for example the Kinoshita–Terasaka knot, with crossing number $\tilde{n}=11$, has $\Delta(t)=1$, the same as the trivial knot [2] (while HOMFLY polynomial can differentiate among the two). However, of all knots with $\tilde{n}<11$ only six cannot be differentiated by the Alexander polynomial. Thus $\Delta(t)$ represents a valuable tool to discriminate between relatively simple knots, which are more likely to occur in nature.

In polymer and bio-polymer physics, knot theory is usually applied to analyze globular configurations, as for example in the study of proteins. In such cases the polymer is represented by a dense polygonal curve, whose projections may have a large number of CPs n, even if the topology of the curve is trivial. This makes the calculation of the Alexander polynomial, or of other polynomial invariants, extremely demanding. To overcome this limitation the analyzed curve should be modified to reduce n, without changing its topology. In other words one has to implement an algorithm that mimics AIs.

An effective algorithm that performs this curve smoothing, or rectification, was proposed Koniaris and Muthukumar in [84]. In this method the polygonal curve is modified by progressively analyzing the triplets of consecutive vertexes. If no curve segment crosses the triangle formed by the considered triplet, then its central vertex is removed, otherwise the algorithm moves to the next triplet. The computation time of this procedure scales as $M\log(M)$, where M is the initial number of vertexes of the polygonal. This calculation allows to minimize the crossings n in any possible projection of the curve, significantly accelerating the computation of polynomial invariants. Algorithms based on the same principle of curve smoothing have lately been proposed in [27, 88]. In many cases curve rectification allows also the visual detection of knots, so it can be considered itself as a knot finding technique [27]. Currently, curve smoothing algorithms are included in all the common procedures for finding knots in proteins [85, 87, 89, 90].

A crucial property of physical knots is their size, namely the length of the shortest curve portion hosting the entanglement. In knotted polymers the size of the knots can deeply modify the equilibrium and dynamics properties of the polymeric chain (see e.g. [91, 92]). The qualification of knot size and localization on the backbone (or depth) requires the topological analysis of all possible portions of the considered chain, or arcs, that is all the combinations of consecutive vertices of the polygonal curve [27, 30, 56, 93]. Once an arc is selected, a closure operation is performed, and then its topology is assessed. The results of such an operation can be collected by means of knot matrices, or fingerprints [56, 94, 95], which allow to effectively visualize the size and location of knots [96]. As shown in figure 7, each entry of the knot matrix indicates (via color or shading) the topology of an arc, whose end vertices are indicated by the row and column indexes. The smallest portion of the chain that embeds a specific topology is the so-called knot core, while the two chain segments excluded by the knot core are called knot tails. Recently, an alternative approach was proposed to display the topology of a knotted closed chain and all its subchains [97], named disk matrix. Based on a longitude-latitude map, the disk matrices are meant to reproduce the periodicity of circularized chain, but can provide also useful insights on the topology of open chains such as proteins.

Zoom In Zoom Out Reset image size

This first model was introduced in the framework of small globular protein folding by Clementi et al [172], in Onuchic group. The polypeptide is represented by a chain of N spherical beads centered at the $C_{\alpha}$ atom positions $\mathbf{R}=\mathbf{r}_{1},\ldots,\mathbf{r}_{n}$. The potential energy of the model is given by equation (6), in which the different potential contributions are defined as follows. Let $\mathbf{d}_{i}=|\mathbf{r}_{i+1}-\mathbf{r}_{i}|$ be the Euclidean distance vector between two consecutive beads, and d_i its magnitude. The bond energy between consecutive beads is:

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \label{Ubond} U_{\mathrm{bond}}=\sum_{i=1}^{N-1}k_{\mathrm{b}}(d_{i}-d_{0i})^{2}, \nonumber \end{align} \tag{ 7 }$$

where d_0i is the native distance between the ith and the i  +  1th residues. The backbone stiffness is represented by restraining the bending and dihedral angles of the protein chain to their native values, namely:

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \label{UbbA} U_{\mathrm{bb}}=U_{\mathrm{bending}}+U_{\mathrm{dihedral}}, \nonumber \end{align} \tag{ 8 }$$

where $U_{\mathrm{bending}}$ is a harmonic potential involving triplets of consecutive beads:

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \label{Ubend} U_{\mathrm{bending}}=\sum_{i=1}^{N-2}k_{\theta}(\theta_{i}-\theta_{0i})^{2}, \nonumber \end{align} \tag{ 9 }$$

in which $\theta_{i}$ is the ith angle, formed by beads i, i  +  1 and i  +  2, and $\theta_{0i}$ is its respective native value. The dihedral term is given by:

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \label{Udih} U_{\mathrm{dihedral}}=\sum_{i=1}^{N-2}\left[k_{1\phi}(1+\cos(\phi_{i}-\phi_{0i}))+k_{3\phi}(1+\cos(3\phi_{i}-3\phi_{0i}))\right], \nonumber \end{align} \tag{ 10 }$$

where $\phi_{i}$ is the ith dihedral angle of the chain, i.e. the angle between two intersecting planes, one containing i, i  +  1 and i  +  2 beads and the other containing i  +  1, i  +  2 and i  +  3 beads. $\phi_{0i}$ is the respective native value. Equation (8) has been used also in other protein models, employing standard, non-specific equilibrium angles (see e.g. [160]).

The $U_{\mathrm{nat}}$ term involves 12–10 Lennard-Jones pair potentials:

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \label{Unat1210} U_{\mathrm{nat}}=\sum_{i<j}^{N-1}C_{ij}\epsilon\left[5\left(\frac{\sigma_{ij}}{r_{ij}}\right)^{12}-6\left(\frac{\sigma_{ij}}{r_{ij}}\right)^{10}\right], \nonumber \end{align} \tag{ 11 }$$

where the contact map C_ij is equal to 1 for residues in contact in the native state, and 0 otherwise. Since $U_{\mathrm{bb}}$ regulates the interaction among consecutive quadruplets of beads, pairs of sequential distance $|i-j|\leqslant3$ are excluded from the native contact map (i.e. C_ij  =  0). The length scale $\sigma_{ij}$ is equal to the distance r_0ij between the ith and j th residues in the native state, so that the minimum of the potential is at $r_{ij}=r_{0ij}$. The non-native term $U_{\mathrm{non-nat}}$ contains only a repulsive contribution:

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \label{Unon} U_{\mathrm{non-nat}}=\sum_{i<j}^{N-1}\tilde{C}_{ij}\epsilon\left(\frac{\sigma}{r_{ij}}\right)^{12}, \nonumber \end{align} \tag{ 12 }$$

where $\tilde{C}_{ij}=1-C_{ij}$ is the complement of the contact map, and includes all possible residue pairs except the native contacts. In equation (12) the length scale is constant, typically $\sigma=4$ $\mathring{\rm A}$. A common choice for the parameters in Clementi's model is: $\newcommand{\e}{{\rm e}} k_{\mathrm{b}}=100\epsilon$, $\newcommand{\e}{{\rm e}} k_{\theta}=40\epsilon$, $\newcommand{\e}{{\rm e}} k_{1\phi}=1.0\epsilon$, $\newcommand{\e}{{\rm e}} k_{3\phi}=0.5\epsilon$, where $\newcommand{\e}{{\rm e}} \epsilon$ is the energy unit corresponding to the depth of the native contact well in equation (11). Since its publication, this model was further developed, for example by implementing the shadow contact map [169], or by introducing a Gaussian-shaped well potential for the native contacts [173]. All these methodologies are available to the public by means of the SMOG [174] and SMOG2 [175] interfaces.

Clementi's description has been widely used in protein folding, and it is widespread also in self-entangled protein studies. To provide few examples, it has been employed to investigate the folding of MJ0366 [105, 171], the smallest knotted protein, or even more complex systems, such as the laboratory engineered 2ouf-knot [176] and the 5₂-knotted Ubiquitin C-terminal Hydrolases [177], and other kind of backbone entanglements, including complex lassos [111, 113, 114] and links [116, 178].

A widely used Gō-like description in the framework of self-entangled protein folding has been initially proposed by Cieplak and Hoang to study universality of protein folding times [179]. The main feature of this model is that the $U_{\mathrm{bb}}$ potential is implemented as a four body term that favors the native chirality of the backbone:

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \label{Uchi} U_{\mathrm{bb}}=\sum_{i=2}^{N-2}\frac{1}{2}\kappa \epsilon \left(\chi_{i}-\chi_{0i}\right)^{2}. \nonumber \end{align} \tag{ 13 }$$

$\chi_{i}$ indicates the chirality associated to the residue i, that is defined, in terms of the distance vectors $\mathbf{d}_{i}$:

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \label{chirality} \chi_{i}=\frac{(\mathbf{d}_{i-1} \times \mathbf{d}_{i})\cdot\mathbf{d}_{i+1}}{d^{3}_{0i}}. \nonumber \end{align} \tag{ 14 }$$

$\chi_{0i}$ in equation (13) indicates the value of the native conformation. This form of $U_{\mathrm{bb}}$ was first proposed in [180], and it is a simpler and more numerically efficient version of the original potential proposed in [179]. The coupling constant $\kappa$ is usually chosen equal to 1. The chirality potential of equation (14) plays an analogous role as the angular terms in Clementi's model, but it has in general a weaker action. It can be demonstrated [168] that equation (14) is approximately equivalent to the dihedral term in equation (8). The bending stiffness is here unrestrained, so that the chain is less stiff than in the previous case. To compensate for this, native contacts between i and i  +  3 residues are included in the contact map. A further difference with previous description consists in the use of a 12–6 Lennard-Jones potential for the native contacts:

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \label{Unat126} U_{\mathrm{nat}}=\sum_{i<j}^{N-1}C_{ij}4\epsilon\left[\left(\frac{\sigma_{ij}}{r_{ij}}\right)^{12}-\left(\frac{\sigma_{ij}}{r_{ij}}\right)^{6}\right], \nonumber \end{align} \tag{ 15 }$$

where the $\sigma_{ij}$ are chosen so that the distance of the energy minimum ($2^{1/6}\sigma_{ij}$) corresponds to the native distance r_0ij.

In the framework of self-entangled proteins, Cieplak model has been used to investigate the folding of the smallest knotted protein [49], of deeply knotted proteins YibK and YbeA [48, 57], and of protein dimers [122, 177], and to explore thermal unfolding [181, 182]. Moreover, as discussed in the following (section 4.4), this description has been employed in combination with methods to simulate protein dynamics under specific conditions, such as the presence of interface environment [182], the translocation through a proteasome channel [183], or co-translational folding [48].

Cieplak's model has been also extensively employed to simulate protein stretching experiments, in which single molecules are manipulated by means of atomic force microscopy [185] or optical tweezers [186], and stretched pulling two selected residues apart (see figure 11(A). The resulting force-extension (F  −  d) diagrams show complex peak patterns (an example is shown in figure 11(B), that can reveal useful information about the structure and stability of the analyzed proteins [187, 188]. MD simulations represent a useful tool for the interpretation of these experimental results [189], generating theoretical F  −  d diagrams that can be reconducted to the experimental ones to clarify the nature of the different force peaks. The stretching force is generated by constraining the position of a selected residue, and by imposing a moving harmonic potential to a second residue. Both constant velocity and constant force pulling procedures have been implemented. An interesting aspect of this simulation protocol is that the simulated velocity of stretching can be set orders of magnitude larger than in experiments, still being able to extrapolate to the experimental timescales. As a result, the computational time required for stretching simulations is relatively shorter and allows also the use of atomistic resolution MD (see section 4.3). Nonetheless, CG models such as Cieplak's have shown to be useful for large surveys of stretching simulations, that could help generalize the mechanical properties of proteins [190, 191]. In the field of self-entangled protein studies mechanical stretching represents a crucial tool of analysis, being capable of revealing possible intermediate conformations, or providing information on the stability of backbone entanglements [63]. Cieplak's model has been widely employed to simulate the stretching of self-entangled proteins [58, 59, 181, 192–194] and protein complexes [122, 177].

Zoom In Zoom Out Reset image size

Comparisons between the two presented variants of CG Gō-models can be found in the literature (see e.g. [168, 195] and, in the realm of entangled folding, [49]), showing e.g. the lower stability of chirality potential with respect to the angular potential of equation (8), or indicating the different ranges of folding temperatures for the two models.

An example of off-lattice Gō-model that has been used in combination with MC sampling is that proposed by Prieto et al in [196, 197]. In this description, the protein is represented as a chain of hard spheres centered on the $C_{\alpha}$ atoms. The new conformations are generated by MC moves, which are constrained by the imposition of excluded volume of the spheres, and by fixed bond distances d  =  3.8 $\mathring{\rm A}$. The energy of a configuration is then determined by native and backbone interactions, which are both described by a harmonic well potential, namely:

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \label{Uwell} U_{\mathrm{bb}+\mathrm{nat}}=\sum_{i<j}u_{w}(r_{ij}), \nonumber \end{align} \tag{ 16 }$$

where the pair interaction is given by:

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \label{Uwell2} u_{w}(r_{ij}) \left\{\begin{array}{@{}ll@{}} w_{ij}[(r_{ij}-r_{0ij})^{2}-a^{2}], & {{\rm if~}} r_{0ij}-a<r_{ij}<r_{0ij}+a \nonumber \\ 0, & {{\rm otherwise}} \end{array}\right., \nonumber \end{align} \tag{ 17 }$$

in which a  =  0.6 $\mathring{\rm A}$, r_0ij is the native distance of the residues i and j  and w_ij defines the contacts. w_ij corresponds to the native contact map but, since there are no angular interactions, also contacts separated by two or three bonds are here considered. In order to preserve the native chirality, when contacts between i and i  +  3, a sign is assigned to r_ij. This description has been employed to assess the effects of local flexibility and steric confinement on the knotted folding of MJ0366, in comparison with the results of a lattice model [198].

The models presented up to here are characterized by sequentially homogeneous interaction, meaning that parameters such as the LJ energy scale $\newcommand{\e}{{\rm e}} \epsilon$, or the angular stiffnesses $k_{\theta}$, $k_{1\phi}$ and $k_{3\phi}$, are independent of the residue index. In these descriptions the sequence information is all represented by the native structure information. This can be an useful simplification to reproduce protein folding only on an approximate, qualitative level. However, in some cases, the lack of interaction specificity can prevent the folding, in particular when the native conformation exhibits a complex topology [48]. In order to improve the effectiveness of protein modelling, many choices of heterogeneous native-contact potentials have been explored, some examples can be found in the systematic study of [193]. One of the simplest choices is that of rescaling the native contacts corresponding to cysteine-cysteine bridges, increasing the depth of the potential well so that the rupture of the bond becomes unlikely [193]. This method has been used in the study of self-entangled proteins for which cysteine bridges are topologically relevant, such as complex lasso structures[111, 113–115] and protein dimers [177]. In some other cases cysteine bridges are treated equally to peptidic bonds (e.g. equation (7)) [59, 111, 113, 114].

An example of how interaction heterogeneity can be crucial for knotted folding can be found in [199], where the atomic interaction-based coarse grained (AICG) model, initially proposed to study allosteric proteins [200], is applied to self-entangled proteins. The AICG approach is built on the potential energy of Clementi's model, in which the coupling constants are made residue dependent. This means that AICG has heterogeneous bond energy stiffnesses ($k_{\mathrm{b}}\to k_{\mathrm{b}i}$), angular stiffnesses ($k_{\theta}\to k_{\theta i}$, $k_{1\phi}\to k_{1\phi i}$ and $k_{3\phi}\to k_{3\phi i}$), and native contacts ($\newcommand{\e}{{\rm e}} \epsilon\to\epsilon_{ij}$). The first step to construct the AICG model is the definition of the relative strength of native contacts. The coupling constant is decomposed as $\newcommand{\e}{{\rm e}} \epsilon_{ij}=\epsilon_{\mathrm{av}}w_{ij}$, where $\newcommand{\e}{{\rm e}} \epsilon_{\mathrm{av}}$ describes the average strength of residue-residue contacts, and w_ij is the relative weight. The w_ij's are computed by means of all-atom (AA) implicit solvent MD in the native state, defining an energy decomposition protocol to retrieve CG contact energies from atomistic interactions [200, 201]. One can either perform an all-atom native state simulation for each protein to be studied, or use a linear regression model proposed in [200], trained over an ensemble of different proteins. After the w_ij's are set, the remaining free parameters of the potential are tuned. The multiplicity of bending and angular parameters is first reduced, assigning each residue to a specific secondary structure such as $\beta$-strands or $\alpha$-helixes, to which a unique set $\{k_{\mathrm{b}},k_{\theta},k_{1\phi},k_{3\phi}\}$ is associated. The reduced set of parameters is then optimized iteratively by matching the average equilibrium fluctuations of CG MD with those of AA MD. Again, one can perform this tuning for any specific system, or use a set of trained parameters obtained from an ensemble of optimizations.

In [199] the AICG approach has been applied by using a backbone potential differing from equation (8), named flexible local potential, constructed from a library of angular distributions via Boltzmann inversion [202], with the aim of reproducing more realistic backbone flexibility. This variant of the model, named AICG2 has been employed to study the folding of the engineered knotted protein 2ouf-knot, of the smallest knotted protein MJ0366, and of the deeply knotted YibK. In the first two cases the AICG2 model obtained a much larger propensity of folding with respect to the homogeneous Gō-model, demonstrating that interaction specificity can play a crucial role in entangled folding. In the case of YibK the folding was instead inaccessible, suggesting that crucial interactions in this system are still neglected at the level of this description.

The models presented until now are purely native-centric, meaning that non-native interactions (the $U_{\mathrm{non-nat}}$ term of equation (6)) are described only by excluded-volume contributions, implicitly assuming that these interactions play a negligible role in determining the folding dynamics. This represents a radical approximation, which is nonetheless supported by accurate MD simulation studies [203]. On the other hand, multiple computational studies have shown that non-native interactions might play a relevant role [204], especially when entangled proteins are considered. As seen earlier, the effects of structural mutations on the lattice model of [154] demonstrated that moving away from the purely native-centric picture can improve the folding propensity of knotted lattice proteins.

Non-native interactions were proven to be crucial in entangled folding by means of a CG off-lattice model in the pioneering work of Wallin et al [155]. In this work the folding of the deeply knotted protein YibK was investigated by means of the previously discussed Clementi Gō-model, equipped with an extra non-native attractive potential term. Let us decompose the non-native contact matrix (see equation (12)) as $\tilde{C}_{ij}=R_{ij}+A_{ij}$, where R_ij and A_ij represent repulsive and attractive contacts, respectively. The non-native potential is then given by:

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \label{Uwal} U_{\mathrm{non-nat}}=U_{\mathrm{steric}}+0.8\epsilon\sum_{i<j}^{N-1}A_{ij}\mathrm{e}^{(r_{ij}-\sigma_{\mathrm{nn}})^{2}/2}, \nonumber \end{align} \tag{ 18 }$$

where $U_{\mathrm{steric}}$ is given by equation (12) and $\sigma_{\mathrm{nn}}=4.0$ $\mathring{\rm A}$ is the attractive contact length. In [155] these contacts are chosen ad-hoc within residues that are involved in the entanglement formation. The inclusion of this attractive term has a striking effect on the folding propensity of the model: while the bare native-centric model was unable to reach the correct topology in all the presented simulations, with non-native interaction the model could always reach the native conformation. These results supported the idea that non-native contacts can regulate the kinetic accessibility of the entangled state, and increase the folding probability of the protein.

Further studies on the impact of non-native interactions in the folding of entangled proteins were obtained by means of a more physically detailed model in [53, 156, 157]. This model, presented in detail in [156] is based on a Gō-like CG description, that can be described by equation (6). Apart from the common harmonic bond potential (equation (7)), this description features heterogeneous angular and native contact interactions [205, 206]. The angular potentials do not include native angles as in equations (9) and (10), but depend on the secondary structure associated to the respective residues. The bending term is given by a double-well potential:

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \label{Ubend2} U_{\mathrm{bending}}=\sum_{i=1}^{N-2}-\frac{1}{\gamma}\ln\left[\mathrm{e}^{-\gamma k_{\alpha}(\theta_{i}-\theta_{\alpha})^{2}-\gamma\epsilon_{\alpha}}+\mathrm{e}^{-\gamma k_{\beta}(\theta_{i}-\theta_{\beta})^{2}}\right], \nonumber \end{align} \tag{ 19 }$$

where $\theta_{\alpha}=92^{\circ}$ and $\theta_{\beta}=130^{\circ}$ are equilibrium values for helical and extended chain respectively, and $k_{\alpha}$,$\newcommand{\e}{{\rm e}} \epsilon_{\alpha}$ and $k_{\beta}$ are specific parameters, typical values are reported in [207]. The dihedral term is a generalization of equation (10):

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \label{Udih2} U_{\mathrm{dihedral}}=\sum_{i=1}^{N-2}\sum_{n=1}^{4}k_{n\phi}\left[1+\cos(n\phi_{i}-\delta_{n})\right], \nonumber \end{align} \tag{ 20 }$$

where $k_{n\phi}$ and $\delta_{n}$ depend on the secondary structure associated to the two central residues of the ith dihedral. The employed values are indicated in [205]. The native contacts are described by the following heterogeneous term:

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \label{Unat12106} U_{\mathrm{nat}}=\sum_{i<j}^{N-1}C_{ij}\epsilon_{ij}\left[13\left(\frac{\sigma_{ij}}{r_{ij}}\right)^{12}-18\left(\frac{\sigma_{ij}}{r_{ij}}\right)^{10}+4\left(\frac{\sigma_{ij}}{r_{ij}}\right)^{6}\right], \nonumber \end{align} \tag{ 21 }$$

where $\newcommand{\e}{{\rm e}} \epsilon_{ij}$ is set as the specific hydrogen bond energy, if a native hydrogen bond is present. Otherwise, $\newcommand{\e}{{\rm e}} \epsilon_{ij}$ is chosen proportional to the quasi-chemical contact potentials of Miyazawa and Jernigan [208], properly renormalized to the hydrogen bond energy scale. As shown in figure 12, the pair term in equation (21) qualitatively differs from the other LJ-like potentials in that a small energy barrier must be overcome to establish the contact, mimicking desolvation energy cost.

Zoom In Zoom Out Reset image size

On top of this native-centric potential, a non-native term and a long-range electrostatic term are introduced. As in the previous case, the non-native contact matrix is divided in attractive and repulsive contributions, defining the following non-native potential term:

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \label{UnnSk} U_{\mathrm{non-nat}}=\sum_{i<j}^{N-1}[A_{ij}U_{A}(r_{ij})+R_{ij}U_{R}(r_{ij})+U_{\mathrm{el}}(r_{ij})], \nonumber \end{align} \tag{ 22 }$$

where U_A and U_R are the attractive and repulsive contact potentials, and $U_{\mathrm{el}}$ describes the long-range electrostatics. Attractive contacts interact via a heterogeneous 12-6 LJ potential:

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \label{UnnA} U_{A}(r_{ij})=4|\epsilon_{ij}|\left[\left(\frac{\sigma_{ij}}{r_{ij}}\right)^{12}-\left(\frac{\sigma_{ij}}{r_{ij}}\right)^{6}\right], \nonumber \end{align} \tag{ 23 }$$

and repulsive contacts interact with the following functional form (displayed in figure 12):

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \label{UnnR} U_{R}(r_{ij})= \left\{\begin{array}{@{}ll@{}} 4\epsilon_{ij}\left[\left(\frac{\sigma_{ij}}{r_{ij}}\right)^{12}-\left(\frac{\sigma_{ij}}{r_{ij}}\right)^{6}\right]+2\epsilon_{ij}, & {{\rm if~}} r_{ij}<2^{1/6}\sigma_{ij} \nonumber \\ -4\epsilon_{ij}\left[\left(\frac{\sigma_{ij}}{r_{ij}}\right)^{12}-\left(\frac{\sigma_{ij}}{r_{ij}}\right)^{6}\right], & {{\rm if~}} r_{ij}\geqslant 2^{1/6}\sigma_{ij} \end{array}\right.. \nonumber \end{align} \tag{ 24 }$$

The length-scale of the contacts is $\sigma_{ij}=(\sigma_{i}+\sigma_{j})/2$ where $\sigma_{i}$ and $\sigma_{j}$ are the van der Waals radii of the interacting residues. The interaction strength is set as $\newcommand{\e}{{\rm e}} \epsilon_{ij}=\lambda(e_{ij}-e_{0})$, where e_ij are negative values taken from the Miyazawa and Jernigan residue interaction matrix [208] and $\lambda$ and e₀ are free parameters, obtained in [207] by fitting calculated binding affinities with experimental values. The average non-native interaction strength resulting from equations (23) and (24) is approximately 1/10 of native one. Finally, the electro-static term is given by Debye–Huckel potential:

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \label{Uel} U_{\mathrm{el}}(r_{ij})=\frac{q_{i}q_{j}}{4\pi\epsilon_{0}D r_{ij}}\exp\left(-\frac{r_{ij}}{\xi}\right), \nonumber \end{align} \tag{ 25 }$$

in which q_i and q_j are the residue charges, $\newcommand{\e}{{\rm e}} \epsilon_{0}$ is the dielectric constant in vacuum, D is the relative dielectric constant of water and $\xi\sim10$ $\mathring{\rm A}$ is the Debye screening length. While the resulting model is definitely richer of physical ingredients than the usual Gō-models, non-native interactions introduce frustration, slowing down the diffusion of the protein along the folding funnel. For this reason, in [53, 156] the sampling is performed via MC, applying a set of local moves that allow to retrieve dynamic properties with a significantly lower cost than with MD [129, 209].

By comparing the results obtained by this model with and without equation (22), it had been possible to enlighten the key role of non-native interaction in favoring early knot formation in carbamoyltransferases [156], and in determining the knotted folding of MJ0366 [53].

As mentioned before, in the study of entangled protein folding one of the issues of Gō-like potentials is the untimely formation of native contacts, that can entrap the chain in a misfolded configuration. These bonds need then to be ruptured in order to reach the native topology. This backtracking mechanism amplifies the simulation time required to reach the folded state, resulting in very low folding probabilities. Nonetheless, it is believed that topologically complex proteins have evolved to fold in a reproducible and efficient way, so their sequence should encode sufficient information to avoid kinetic traps and misfolded configurations. Building on this idea, Najafi and Potestio have proposed a CG off-lattice description, dubbed elastic folder model (EFM), which is constructed to fold along optimal pathways, that is in the most reproducible and rapid way.

To achieve this the EFM employs a typical CG $C_{\alpha}$ description of the protein, in which native contact potentials are absent and the folding is governed by the back-bone angular interactions alone. The energy can be written as:

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \label{Uefm} U_{\mathrm{tot}}=U_{\mathrm{bond}}+U_{\mathrm{steric}}+U_{\mathrm{bb}}, \nonumber \end{align} \tag{ 26 }$$

where the Kremer–Grest model [210] is used to represent peptidic bonds and steric interactions. The former are described with the finitely extensible non-linear elastic potential:

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \label{Ufene} U_{\mathrm{bond}}=-\sum_{i=1}^{N-1}\frac{k_{\mathrm{FENE}}}{2}\left(\frac{R_{0}}{\sigma}\right)^{2}\ln\left[1-\left(\frac{r_{i,i+1}}{R_{0}}\right)^{2}\right], \nonumber \end{align} \tag{ 27 }$$

in which $\newcommand{\e}{{\rm e}} k_{\mathrm{FENE}}=30.0\epsilon$ is the interaction strength parameter, in units of an energy constant $\newcommand{\e}{{\rm e}} \epsilon$, $R_{0}=1.5\sigma$ and the length-scale $\sigma=3.8$ $\mathring{\rm A}$ is the typical extension of a peptidic bond. The Weeks–Chandler–Anderson interaction [211] is used for the steric part:

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \label{Usteric} U_{\mathrm{steric}}=\sum_{i<j}^{N}U_{\mathrm{WCA}}(r_{ij}), \nonumber \end{align} \tag{ 28 }$$

with:

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \label{Uwca} U_{\mathrm{WCA}}= \left\{\begin{array}{@{}ll@{}} 4\epsilon\left[\left(\frac{\sigma}{r_{ij}}\right)^{12}-\left(\frac{\sigma}{r_{ij}}\right)^{6}\right]+\epsilon\qquad &{{\rm if}}\ r<2^{1/6}\sigma \nonumber \\ 0\qquad &{{\rm otherwise}}. \end{array}\right. \nonumber \end{align} \tag{ 29 }$$

The structural information of the protein is included in the back-bone angular potential, which has the same functional form used by Clementi model (equation (8)). The driving force of the folding is thus provided by the bending and torsion terms. The idea of folding pathway optimality is implemented by using heterogeneous angular stiffnesses $k_{\theta i},k_{1\phi i},k_{3\phi i}$, which are tuned through a stochastic optimization approach, aimed at maximizing the folding propensity of the model. This optimization procedure is based on iteratively tuning the parameters $k_{\theta i},k_{1\phi i},k_{3\phi i}$ building on the results of several folding simulations. In [55] this is done by a MC sampling of the parameter space, in which a mutation of the stiffness is accepted or rejected according to its effect on the folding probability. This technique has been used to investigate the folding routes of two small knotted proteins [55], observing a qualitative agreement with the results of more detailed models. The optimization procedure has been recently improved via an evolutionary, parallelized strategy that allows a more efficient exploration of the parameter space [115].

During decades of computational protein studies distinct CG descriptions of the polypeptide chain have been proposed, naturally increasing the detail starting from the simple $C_{\alpha}$ representation. However, mainly for computational requirements, these models have not been extensively employed in the field of entangled proteins.

An example of a more detailed description of the protein chain is found in [212], where the thermal unfolding of YibK is simulated with different techniques. Among these, wo CG models were employed, both increasing the detail of $C_{\alpha}$ representation by describing further elements of the amino-acid chain. The first one, described in [213], is a MC approach which includes side-chain and $C_{\beta}$ interaction sites. The Hamiltonian contains heterogeneous backbone and native contact interactions, that combine general potential forms for secondary structures, native structure information, and statistical potentials built over a database globular proteins [214]. The second model as well includes a side-chain representation, and employs interaction potentials parametrized over known protein structures and ab initio calculations [215–217]. This second description is employed with MD sampling.

A second example of detailed CG in knotted proteins can be found in [218], where the folding of tRNA methyltransferase, presenting a deep native trefoil knot, is addressed. Actually only the entangled fragment of the protein chain is considered in the presented calculations. The used representation is inspired from the associative memory Hamiltonian used in structure-prediction studies [219, 220], that describes the structure retaining the position of $C_{\alpha}$, $C_{\beta}$ and O atoms. The energy function is the sum of a back-bone term, including self-avoiding and stiffness interactions, and a heterogeneous gaussian Gō-like term for $C_{\alpha}-C_{\alpha}$, $C_{\alpha}-C_{\beta}$ and $C_{\beta}-C_{\beta}$ pairs. The detailed form of the Hamiltonian can be found in [218, 220]. Another peculiar aspect of this work is the use of structure prediction methods to investigate the kinetics of the model. In these techniques the sampling of the folding space is obtained by generating a complex network of the possible local energy minima, using a combination of MC moves, minimum energy path, and energy optimization algorithms [221–223]. The transition rates across this network of conformations are then computed using discrete path sampling method [224], and the global kinetics of the model is inferred.

Many of the studied proteins constitute functional homo-dimers. In such cases it can be relevant to simulate the cooperative folding of two protein specimens, and observe their interactions and dimerisation process [173]. As previously mentioned, in [105] the folding of MJ0366 was simulated via an AA Gō-model. In this paper, the authors also investigated the possible dimerisation along or after the folding process, modelling the native contacts between two MJ0366 monomers via Gō-potential, with the same interaction strength and functional form (Gaussian) employed for the intramonomer native contacts. The monomers were held together by an harmonic potential acting on their centers of mass. Varying the curvature of this potential, the effect of different crowding on folding and dimerisation was probed. The interaction between protein dimers was simulated also in [122, 178, 246], where the study of both folding and stretching was addressed, this time using the CG Gō-model of Cieplak. Once again, the contacts between monomers were treated like the intra-chain native contacts. Using the same approach, in [122], the folding of known homo-dimeric knotted proteins, such as MJ0366 and YibK, was simulated, assessing the impact of the dimerisation potential on the process.

Cieplak Gō-model has been employed to simulate also other dimeric structures, forming protein links. In [122], the authors propose a four-terminal pulling simulation protocol, that detects the intertlinking among homo and hetero-dimer structures. This technique has been applied to more than 10⁴ n-meric structures, detecting about 9% of entangled complexes. Four of these linked dimers have been studied in [178], assessing their folding and thermodynamic properties. In general the full dimerisation is found to occur as a late step, subsequent to the proper folding of the monomers.

Another aspect that plays a crucial role in protein folding in vivo, is the translation [52]. Indeed, the folding dynamics can take place already while the polypeptide chain is emerging from the ribosome, and this possibility can completely alter the picture of spontaneous folding adopted in standard simulations. This 'co-translational' folding can be particularly relevant for self-entangled proteins, as the sequential character of topology formation can be regulated by the time-scale of translation. This possibility is addressed in [48], where a simple model of on-ribosome protein folding is proposed. In this model the protein is represented by means of the Cieplak CG description, introduced before, while the ribosome is represented as a plane that generates a uniform potential:

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \label{ribo} U_{\mathrm{r}}=\frac{3\sqrt{3}}{2}\epsilon\left(\frac{\sigma_{0}}{z}\right)^{9}, \nonumber \end{align} \tag{ 30 }$$

where z is the distance from the plane and $\sigma_{0}=4\times(2){}^{-1/6}$. The protein chain emerges from a fixed coordinate on the plane, giving birth to a new $C_{\alpha}$ residue every t_w, starting from the N-terminal. This simple representation of the process aims at modelling the sequential generation of the polypeptide, and the excluded volume determined by the ribosome. This technique was used to simulate the folding of knotted proteins [48, 49], observing an improvement in the folding propensity. A simpler, but relatable study was performed in [158], where one terminal of the lattice protein chain was constrained to a chemically inert plane (determining only excluded volume). This method can provide a rough approximation of the entropic limitations present in nascent folding, but also in single molecule experiments. An alternative, more detailed, approach to simulate co-translational folding was recently proposed in [50], where the folding of the deeply knotted protein Tp0624 is simulated by modeling the ribosome channel with a cylindrical 10 $\mathring{\rm A}$ tunnel with a funneled exit located on a planar wall. Initially, the stretched protein is contained in the tunnel, and a constant force is applied to its residues in order to push them through the channel exit. The protein description adopted is the Clementi model, with Gaussian native contacts. The interaction of the wall with the residues can be repulsive (WCA potential, equation (29)) or attractive (LJ 12-6, equation (23)). By assigning attractive interactions to a pre-selected set of residues this scheme could promote the formation of a loop on the ribosome wall at the channel exit. Starting from this configuration the authors demonstrated that the tying of the deep knot could be strongly favored by the ribosome action.

An external factor that can play a key role in protein dynamics is the interaction with an air-water interface [247]. The vicinity of such interfaces had been also found to influence the topology of proteins [182], e.g. favoring non-native entanglements, or untying shallow knots. The method used here is again based on Cieplak's CG Gō model, where the air-water interface is effectively described by a field of forces coupled to the hydropathy of amino acids. The force acting on the ith residue is defined as follows:

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \label{awf} F^{\mathrm{wa}}_{i}=q_{i}A\frac{\exp\left(-z^{2}_{i}/2W^{2}\right)}{\sqrt{2\pi}W}, \nonumber \end{align} \tag{ 31 }$$

where q_i is the hydropathy index associated to the amino-acid [248], z_i its coordinate on the axis perpendicular to the interface, and $\newcommand{\e}{{\rm e}} A=10\epsilon$ and W  =  5 $\mathring{\rm A}$ are energy and length-scale of the interaction.

It is well known that the folding of proteins can be assisted, and accelerated, by the bacterial chaperonin GroEL-GroES, a complex of two large heptameric units that form a compartment capable of accomodating a folding polypeptide [249]. The encapsulated protein exhibits a faster folding rate, which can be the result of an active action of the chaperonin, via e.g. destabilisation of misfolded configurations, and a passive action, via the steric confinement determined by the compartment. The GroEL-GroES complex is of interest also in the field of self-entangled folding, as it has been shown to assist the folding of knotted proteins such as YibK and YbeA [8, 43]. To better understand the action of the chaperonin, computational techiques that simulate folding under spatial confinement have been proposed [250]. Of interest for topologically complex folding is [198], where the confinement effect was studied by means of both lattice (see section 4.1) and off-lattice CG Gō models (see section 4.2.3). The methodology simply consists in limiting the space accessible to the protein chain to a cubic or spherical cavity, the size of which is chosen accordingly to the typical chaperonin size (few nm of radius). Moreover, in [177, 251], the chaperonin cage was modeled by means of a cylindric repulsive cavity, following [250]. The spatially confined folding of different entangled proteins was here studied using Clementi's CG Gō model. Overall, the steric effect of the cage is found to enhance the capability of folding of the protein models, reducing the folding time and promoting different routes with respect to bulk simulations.

Crucial biological processes, such as mithocondrial import, or degradation, require proteins to be translocated through nanopores of 12–14 $\mathring{\rm A}$ of minimum diameter. This size is too narrow to accomodate a native folded structure, and the protein needs to be unfolded by the action of unfoldases, that employ energy from ATP hydrolysis to perform this operation. It has been proposed that topologically complex proteins could prevent such a process, as it was shown that a mechanically tightened protein knot has a size comparable or larger than the pore diameter [28, 252]. It is therefore of interest to understand how the cellular machineries that operate protein translocation can cope with the presence of knots, untying the chain before it jams the pore channel.

In simulations this is attempted by means of an external repulsive potential that models the pore channel, and by applying a pulling force on one extremity of the protein model, to drive its translocation through the pore potential. Once again, because of computational time requirements, the preferential protein descriptions for these simulations, are simple $C_{\alpha}$ representations of self-entangled proteins. As a first example, in [253], specifically designed knotted polypeptides were represented by the model of [160, 166], and pulled by a constant force directed along the axis of the pore. The latter was represented by a cylindrical repulsive potential of length $L\approx50$ $\mathring{\rm A}$ and radius $\rho\approx6.5$ $\mathring{\rm A}$. Different topologies were tested, showing that the knotted polypeptide can still translocate through the pore, but the rate of the process is significantly reduced.

A similar strategy was employed in [254], where models of real knotted proteins such as YibK, YbeA and MJ0366 were tested. Cieplak Gō-model was used for the representation of the protein, and a cylindrical repulsive potential modeled the pore, while a constant pulling force was used to drive the translocation dynamics. To mimic the transport into mitochondria, which is mediated by a short polypeptide sequence attached to the protein end, a 10-bead unstructured chain was added to the end of the protein to be pulled through the pore potential. The effective radius of the modeled pore was similar to the previous case, $\rho\approx7$ $\mathring{\rm A}$, but the results demonstrated how deep protein knots can get stuck at the pore entrance. In [255] the same techniques and test systems were employed, but a periodic pulling force was applied to the protein terminals. This cyclic force, which is more realistic in reproducing an ATP-hydrolisis activated behaviour, could resolve the translocation jam by letting the knot slide off the chain.

A different pore-model has been introduced in [256], to represent the translocation into the proteasome. In this model, the entrance ring of the proteasome is described as a repulsive torus with major and minor radius of 13 and 6 $\mathring{\rm A}$ respectively, placed at the end of a cylindric channel of radius $\rho\approx7.5$ $\mathring{\rm A}$, which mimics the proteasome chamber. These two elements share the longitudinal axis, determining a cylindrical channel with a smooth, funneled entrance of 7 $\mathring{\rm A}$ radius (the torus hole). Also here, proteins were described by means of Cieplak structure-based model, while the pulling of the protein inside the channel was performed either with a constant force or a constant velocity protocol. In [183] this methodology was used to simulate translocation of knotted globular proteins such as YibK and MJ0366 and transiently knotted polyglutamine tracts, assessing the dependence of the translocation process on protein topology, pulling force and pore model. Also a more rugged pore entrance was tested, by substituting the torus with 12 overlapping spheres of 6 $\mathring{\rm A}$ of radius, arranged in circle. This model was also dynamic, implementing the temporary shrinking of the entrance spheres, three at a time. This behaviour is meant to reproduce the allosteric transformations of the proteins that compose the proteasome entrance ring, providing a more realistic description of the biological process.