The Effect of a Knot on the Thermal Stability of Protein MJ0366 : Insights from Molecular Dynamics and Monte Carlo Simulations

We use the approach of topological solitons developed by Antti J. Niemi Niemi (2014). It is based on representing the protein as a partially linear chain consisting of $C_{\alpha}$ atoms and links that connect them (backbone). The positions of all $C_{\alpha}$ atoms in the protein determine the positions of all other atoms due to rigidity of the peptide group.

This curve is considered in terms of discrete Frenet coordinates. In each residue, the Frenet frame (Fig. 2) consisting of three vectors $(\vec{t},\vec{b},\vec{n})$ is built:

{dgroup}

$$\vec{t}_{i}=\frac{\vec{r}_{i+1}-\vec{r}_{i}}{|\vec{r}_{i+1}-\vec{r}_{i}|},$$

(1)

$$\vec{b}_{i}=\frac{[\vec{t}_{i-1},\vec{t}_{i}]}{|[\vec{t}_{i-1},\vec{t}_{i}]|},$$

(2)

$$\vec{n}_{i}=[\vec{b}_{i},\vec{t}_{i}],$$

(3)

where $\vec{r}_{i}$ are coordinate vectors of $C_{\alpha}$ atoms.

The angles of curvature $\kappa_{i}$ and torsion $\tau_{i}$ (Fig. 2) act as new coordinates. They can be calculated as follows: {dgroup}[noalign]

$$\kappa_{i}=\arccos(\vec{t}_{i+1},\vec{t}_{i}),$$

(4)

$$\tau_{i}=\mathrm{sign}([\vec{b}_{i-1},\vec{b}_{i}],\vec{t}_{i})\arccos(\vec{b}_{i+1},\vec{b}_{i}).$$

(5)

The free energy function is based on the Abelian Higgs model.

$$F=\sum_{i=1}^{N-1}\left(\kappa_{i+1}-\kappa_{i}\right)^{2}+\sum_{i=1}^{N}\left(\lambda(\kappa_{i}^{2}-m^{2})^{2}+d\kappa_{i}^{2}\tau_{i}^{2}\\ -b\kappa_{i}^{2}\tau_{i}-a\tau_{i}+c\tau_{i}^{2}\right)+\sum_{i>j}^{N}V(|\vec{r}_{i}-\vec{r}_{j}|),$$

(6)

where the terms with coefficients $a$ and $b$ stand for conservation of charges in discrete non-linear Schrodinger equation hierarchy called helicity and momentum, this terms provide parity breaking and make protein right-handed chiral. Terms with coeficients $c$ and $d$ together constitutes the Kirchhoff energy of an elastic rod. The last term is hard-core Pauli repulsion to avoid self-crossings.

This energy function possesses solutions in the form of topological solitons, and the protein model is built as a composition of solitons (multisoliton). Each soliton describes a region of the protein consisting of parts of regular structures ($\alpha$-helices and $\beta$-strands) and an intermediate structure between them (loop).

The sets of coefficients $(\lambda,m,a,b,c,d)_{k}$ are individual and uniform for each soliton. Based on these coefficients, the native state can be calculated as a configuration satisfying to the minimum of (6):

{dgroup}

$$\frac{\partial F}{\partial\kappa_{i}}=2(2\kappa_{i}-\kappa_{i-1}-\kappa_{i+1})+4\kappa_{i}(\kappa_{i}^{2}-m^{2})+2d\kappa_{i}\tau_{i}^{2}{-2b\kappa_{i}\tau_{i}=0,}$$

(7)

$$\frac{\partial F}{\partial\tau}={2d\kappa_{i}^{2}\tau_{i}-b\kappa_{i}^{2}-a+2c\tau_{i}=0.}$$

(8)

Taking into account that solitons usually aggregate several amino acid residues, the number of parameters needed to describe the protein structure is often less than the amount of amino acids that make up the protein.

As a starting point for modeling the protein structure, we use a structure taken from RCSB Berman (2000). Using a simulated annealing algorithm Pincus (1970), the soliton parameters are chosen such that the configuration corresponding to the minimum free energy corresponds to the initial structure taken from the RCSB data base with a minimum RMSD value.

The obtained structure can be used for simulation of the unfolding-folding process using Monte Carlo methods. To do this we use Glauber algorithm Glauber (1963); Bortz et al. (1975). In comparison to Metropolis-Hastings algorithm it allows us to choose perturbing site randomly regardless if it was chosen before.

At each Monte Carlo step $n$ we independently change either the torsion angle or the curvature angle at a randomly chosen site $i$:

{dgroup}

$$\tau_{i}\rightarrow\tau_{i}+r\sigma_{\tau},$$

(9)

$$\kappa_{i}\rightarrow\kappa_{i}+r\sigma_{\kappa},$$

(10)

where $r$ is a normally distributed around $r=0$ random number with dispersion $\Delta r=1$ and $\sigma_{\tau}=\sigma_{\kappa}=0.01$ is proportionality factor. And then compute the corresponding change $\Delta F(\kappa_{i},\tau_{i})$ in the free energy. The obtained configuration is being accepted with probability

$$P=\frac{1}{1+e^{\Delta F/T_{G}}},$$

(11)

where $T_{G}$ is the Glauber temperature factor that allows to simulate temperature changing. It is important to note that it is Monte Carlo temperature and it can not be concided with the physical temperature measured in Kelvin directly. The relation can be found by the methods of renormalization Krokhotin et al. (2013); Peng et al. (2016).

It has been demonstrated Martinelli and Olivieri (1994a, b) that Monte Carlo simulation using equation (11) converges exponentially to the Gibbsian distribution. Therefore, this approach should offer a statistically significant representation of protein folding dynamics near equilibrium during adiabatic temperature changing.

The simulation process consists of three stages: heating (unfolding), thermalization, cooling (folding). During heating and cooling temperature is increasing (decreasing) by exponential law with growth factor defined by required amount of steps and temperature range:

$$T_{G}^{\prime}=T_{G}b,$$

(12)

$$b=ln\frac{T_{max}}{T_{min}}N^{-1},$$

(13)

where $N$ is number of steps required for each stage. This choice allows to guarantee adiabatic nature of the process.

Knots as topological objects are defined only on closed loops while proteins mostly are not closed. To avoid this different techniques can be applied.

To detect whether protein is knotted or not we used Knot_pull—python package Jarmolinska et al. (2019) by Aleksandra Jarmolińska. Knot_pull refines and streamlines the given structure with topological awareness, ensuring that the inherent chain does not intersect with itself (Fig. 3). Smoothing continues until no further modifications are possible (or the current arrangement has been previously examined), any unnecessary beads are eliminated, retaining only the essential topological points. Based on this structure, the corresponding Alexander-Briggs Alexander and Briggs (1926) notation is calculated.

This approach is more suitable than others based on the calculation of polynomial invariants for the determination of knot type Lua (2012); Tubiana et al. (2018); Dabrowski-Tumanski et al. (2018b, 2020). To detect knot type these approaches close the loop randomly on infinite radius sphere (Fig. 4), such closure repeats many times and gives only the probability of observing given knot type.

The protein MJ0366 was selected for investigation due to its status as the smallest known knotted protein. Its relatively small size and knotted topology make it a particularly interesting subject for studying the structural and functional implications of protein knots, providing a simpler model to understand the biological significance and stability of these complex protein structures.

The MJ0366 protein is a hypothetical protein from the archaeon Methanocaldococcus jannaschii Thiruselvam et al. (2017), which is a hyperthermophilic, methane-producing microorganism found in extreme environments such as hydrothermal vents.

M. jannaschii is a member of the Archaea domain, a group of prokaryotes that often thrive in extreme environments. It grows optimally at temperatures above $85\celsius$ and can reduce carbon dioxide to methane using hydrogen, a process known as methanogenesis. This species is used as a model organism to study extremophilic life and early evolutionary processes due to its ancient evolutionary lineage and unique biochemical pathways.

The KnotProt 2.0 package Dabrowski-Tumanski et al. (2018b) detects the knot of type $3_{1}$ whose core consists of 71 residues, from 6 to 77 amino acids (Figure 6). The knot loop in MJ0366 consists mostly of helices, that probably makes it more tough and thus more stable.

Although the exact biological function of MJ0366 remains hypothetical, knotted proteins, in general, are known for their mechanical and thermal stability. Some knotted proteins have been implicated in various biological processes, including those related to diseases like Parkinson’s disease.

The study of MJ0366 sheds light on the structural complexity of knotted proteins and their potential stability in extreme conditions.

For calculations model 2EFV from PDB database was used.

Knotted protein 2EFV consists of 82 amino acids which form 2 $\beta$-strands and 5 $\alpha$-helices. It means that there are 8 loops and the structure must be describe by at least 8 individual solitons. But $\beta$-strands because of their complicated geometry require additional solitons in bending positions. Finally, we construct the native structure of the protein using 12 solitons. It allows us to reach resolution of 0.52 Å RMSD that is less than reported resolution 1.90 Å of experimental data.

The results of constructing are summarized in Figure 5. In panel a) overlay of 3d models experimental (green) and soliton (purple) structures is represented. Panel b) shows how the angles of curvature and torsion behave in experimental data and in the soliton model.

We performed simulations of the unfolding-folding process and inspected 3d models of the protein.

We found that 2EFV does not fold back to its native state once knot was untied after slow enough heating to the temperature of self-avoiding random walk (i.e. denaturation). It becomes clear from the graphs of radius of gyration $R_{G}$ (Figure 7). In the panel c) enlarged area around $T_{G}=10^{-12}$ preciding to unfolding process is shown. The radius of gyration undergoes a shrinkage that is connected to unknotting of the molecule Lundgren and Niemi (2009). Meanwhile in panel d) which corresponds to the same area, but during the cooling process, shrinkage is totally absent.

We studied the temperature range in which protein can be folded back to the native state. To do that, we heated the protein to the required temperature ($N_{steps}$), did a short thermalization ($N_{steps}/10$) and then cooled it back to $T_{G}=10^{-17}$. We present probability distributions of RMSD after heating and cooling to a certain temperature, the results are represented in Figure 8.

We checked all obtained configurations by Knot_pull and found that all and only configurations which belong to the first peak must be considered as knotted ones. Other peaks are products of misfolding to local energy minimums.

The series of simulation was provided. We studied the dependences of the percantage of knotted configurations after heating and cooling from the maximum temperature in the heating-cooling process and the amount of Monte Carlo steps in the process. As follows from trajectories analysis the unknotted configuration can not be tied back, and after the cooling if the configuration is still knotted, it means that it was not untied. So, we can interpret this graph as percent of untied configurations.

The results obtained were fitted by the sigmoid function

$$f(T_{G})=\frac{a}{(1+e^{-b(\log T_{G}-\log T_{\mathrm{crit}})})}.$$

(14)

The parameter $T_{\mathrm{crit}}$ in (14) can be interpreted as as a critical temperature for particular heating velocity at which protein denatures irreversibly. We plot extracted from fit values of critical temperature and fit it with linear regression. The critical temperature at infinitely slow heating turns out to be

$$\log_{10}T_{\mathrm{crit}}=-11.18\pm 0.05.$$

Thus, we see that 2EFV is relatively thermally stable, especially at high heating velocities. We propose that the high thermal stability of 2EFV is connected to the knot in its structure. To prove that we start simulations not from native state but from unknotted deeply cooled configuration. The unknotted configurations we chose from a huge variety of configurations cooled after unknotting. We made several calculations starting from different configurations closest to the native state having RMSD in a range from 3 to 5 Å  and radius of gyration from 12.20 to 12.30 Å.

It is clear from the graphs that the unknotted configurations starts unfolding earlier for two orders of magnitude of Glauber temperature than the native knotted configuration. Also, it does not shrink, and we can conclude that shrinking before unfolding is connected with the unknotting of the structure.