A numerical technique for studying topological effects on the thermal properties of knotted polymer rings
Abstract
The topological effects on the thermal properties of several knot configurations are investigated using Monte Carlo simulations. In order to check if the topology of the knots is preserved during the thermal fluctuations we propose a method that allows very fast calculations and can be easily applied to arbitrarily complex knots. As an application, the specific energy and heat capacity of the trefoil, the figureeight and the knots are calculated at different temperatures and for different lengths. Shortrange repulsive interactions between the monomers are assumed. The knots configurations are generated on a threedimensional cubic lattice and sampled by means of the WangLandau algorithm and of the pivot method. The obtained results show that the topological effects play a key role for shortlength polymers. Three temperature regimes of the growth rate of the internal energy of the system are distinguished.
I Introduction
Topological polymers are a fascinating subject by themselves. Moreover, the interest in this topic is motivated by possible industrial applications, for example because knots in polymer materials are able to influence the viscoelastic properties grosberg (), and due to their relevance in biology and biochemistry, see for instance dna (); katritch96 (); katritch97 (); krasnow (); laurie (); cieplak (); liu (); wasserman (); sumners (). Polymer rings may be found in nature or in artificial materials in the form of knots that can be additionally linked together. In this work we explore the case of the statistical mechanics of a single polymer knot. Various topological configurations are taken into account: The unknot, denoted according to the AlexanderBriggs notation , the trefoil , the figureeight and the knot . In particular, we study the thermal properties of the above knots by computing numerically their internal energy and heat capacities. The most important problem of numerical simulations of polymer knots is to preserve their topology while their configurations are fluctuating randomly. To this purpose, several methods have been developed vologodski (); orlandini (); kurt (); mehran (); pp2 (); metzler (); muthu (). Most of them, are based on the Alexander or Jones polynomials, which are rather powerful topological invariants able to distinguish with a very high degree of accuracy the different topological configurations. The main drawback of these polynomials is that their calculation is very expensive and time consuming from the computational point of view.
Here we propose a new strategy that consists in the following. First, we apply the pivot algorithm to change a given knot configuration. Starting from a seed knot, for instance those given in ref. scharein (), we change it at randomly chosen points by applying the pivot algorithm pivot (). After each pivot move, the difference between the old and new configurations forms a closed loop. Around it, an arbitrary surface is stretched, whose boundary is the closed loop itself. The criterion to reject changes that destroy the topology of the knot is the presence or not of lines of the old knot that cross such surface. If these lines are present, the trial pivot move is rejected. This combination of pivot algorithm and excluded area (PAEA) provides an efficient and fast way to preserve the topology that can be applied to any knot configuration, independently of its complexity. For pivot moves involving a small number of segments the method becomes exact. This technique may be employed in the study of the thermal and mechanical properties of polymer knots. Here we have limited ourselves to the computation of the internal energy and heat capacity in the case in which the monomer interactions are of very short range, but the extension to any other interactions or to the calculations of more complicated physical quantities is possible. The configurations of the unknot, trefoil, figureeight and knot are modeled on a cubic lattice. The sampling of the canonical ensemble is achieved by using the WangLandau algorithm wanglaudau () at different temperatures.
The rest of this paper is organized as follows. In Sec. II, we describe
the sampling and simulation methods of the polymer knots.
Our strategy to distinguish possible violations of the topology
of a given knot after a random change of its trajectory is explained.
The
calculation of the specific
energy and the specific heat
capacity for various polymer topologies and for different polymer lengths is
presented in section III. In section IV we draw our conclusions and
discuss the open problems together with new possible directions of our
research.
Ii Model and Computational method
In the current work, we will study the topological effects on the thermal properties of polymer rings. Four types of polymers will be considered on a 3dimensional cubic lattice, namely the unknot , the trefoil , the figureeight and the knot , see Fig. 1.
A general knot will be denoted with the symbol . To simulate the trivial polymer ring , we generate a loop starting from a random point and performing a selfavoiding walk on the cubic lattice. The walk is stopped when it crosses itself, forming a loop. We generate many loops of this kind and we retain only those with a given length . For nontrivial topologies, we choose a different strategy. We begin from a seed conformation and then we act repeatedly on it with the pivot method pivot () to get random conformations. The elements of the loop that will be changed by the pivot move have a fixed number of segments or . The elements begin at a random point in the knot, let’s say and end up at the point number , or . can range from 0 to ( is the length of the polymer). We check the topology of the configuration after every pivot move, to be sure that the move will not change the topology of the given knot. Usually, this check is achieved by projecting the knot along some direction on a two dimensional plane and the relative Alexander polynomial is computed. This method has been pioneered in refs. vologodski (). Here we choose another strategy, namely on a given knot we perform a pivot move. In that way a new knot is obtained. A schematic picture of the situation can be found in fig. 2. In order to verify that the topology of the knot has not been changed by the pivot move, we check that no lines of cross the area spanned by the small loop including its border. Otherwise, this could imply that a crossing of the lines of has happened. The definition of the area of in a numerical simulation is a difficult problem of computer graphics. However, if the number of segments involved in the pivot moves is small, it is possible to classify all the geometries of and to stretch a suitable surface around . Next, these data, namely all geometries of and the coordinates of the relative surfaces , are supplied to the code, which will then be able to detect the occurrence of possible lines of crossing or the part of its border , see fig. 2 for a schematic representation of what these symbols mean. If these lines are present, the trial pivot move is rejected.
Let us note that the strategy described above could be considered as
an implementation of the topological
constraints based on the Gauss linking number modified to cope with
the present situation, in which the Gauss linking number becomes not
suitable because and share part of their
trajectories. While the Gauss linking number is a weak topological
number,
our method turns out to be
very efficient when
applied to pivot modes with
segment lengths and . In both cases the geometry of the
loop is simple enough, that it is possible to analyze all
situations in which the crossings of lines could violate the topology of and to
provide an exact criterion
for preserving this topology without increasing exceedingly the
simulation time. Only when it may occur
that the trajectory of the old knot crosses the area
an even number of times.
In this case, due to limitations similar to those of the Gauss
linking number,
allowed moves can
be rejected. If this happens very frequently,
the computational time may increase. However,
this is not a serious drawback because
the above procedure of detecting topology changes is
extremely time efficient. The necessary cpu time depends essentially
only on the number of segments , independently on the complexity of
the knot.
In the following, we will restrict ourselves to the cases
and .
In order to calculate the internal energy of the polymer per unit of
length ,
where denotes
the Boltzmann factor with (natural thermodynamic units) and
the heat capacity
(1) 
the density of states of polymers is needed. To that purpose, the WangLandau (WL) algorithm is used for sampling the polymer conformations. In the case of the trefoil, figureeight and knot, before applying the WL algorithm, the system has been equilibrated, since we start from a given seed which is not in a random conformation. The fluctuations of the radius of gyration are employed to judge the degree of equilibration. These fluctuations are computed according to the formula rq ():
(2) 
Here is the sweep time, is the average radius of gyration at time . The equilibrating procedure is stopped after the time dependent fluctuation of the radius of gyration oscillates around a constant. To crosscheck the equilibration of the nontrivial knot configurations, we use also another method, that consists in counting the number of sites in which the direction of the chain is not changed tk (). This number is a random variable, but at equilibrium its variance around the average is small. The results obtained with both techniques, that based on the fluctuation of the gyration radius and that based on the changes of direction of the segments, are in agreement.
Fig. 3, which is related to the latter method, shows that the system equilibrates after a few tens of thousands of pivot moves in the case of the trefoil . Similar considerations are valid for the other knots considered in this work. In total the system is subjected to one million of pivot moves before treating it with the WL algorithm in order to compute the energy and the heat capacity .
Now we briefly introduce the WL algorithm. This method is considered as a selfadjusting procedure for obtaining the density of states. In the present context the states are distinguished by the number of closest contacts between the polymer monomers, where takes positive integer values. Following a previous work pn (), all initial values of the density of states are firstly set to be equal to one, i. e. for all the states. Then the pivot method is used to visit the next trial state . The energy of the trial state will be denoted by , where is the number of the closest contacts of the monomers, while is the interaction energy of each monomer pair. and correspond to repulsive and attractive case respectively. We assume in the following that . The probability of translating from the state to is
(3) 
Once an energy state is visited, the corresponding density of states
of this state will be updated by multiplying a modification
factor wanglaudau (). After the th sweep, the modification
factor
is reduced according to the equation . The above
process is repeated until the modification factor reaches the predefined limit
wanglaudau ().
Of course, at each step the preservation of the knot topology is
checked with the method explained before.
We show as an example the obtained densities of
states for the unknot, the trefoil and the figureeight of length
in fig. 4.
The obtained distribution of density of energy states can now be used to calculate the canonical averages of the observables according to standard relations. As already mentioned, the observables studied in this paper are the internal energy of the canonical system and the heat capacity as a function of the temperature. The expression of the internal energy is given by:
(4) 
where is the number of the closest contacts of polymers and is the density of energy states computed by means of the WL algorithm. The formula of the heat capacity has been already provided in Eq. (1).
The computational method has been checked using the trivial knot as a test. We could repeat the previous results of pnring () in the case of chains with lengths equal to 12, 30 and 50.
Iii results and discussion
In this section we study the thermal properties of the trefoil and figureeight with lengths in the range of segments for the case . The used algorithm is quite fast. Depending on the length, the typical simulation times vary between 15 minutes to 1 hour on a dedicated pc. It can be applied to much more complicated knots. As an example, we have analysed here also the case of the knot.
The results concerning the trefoil and figureeight are presented in fig. 6 and fig. 7 respectively. For comparison, fig. 5 represents the results for with lengths ranging from to . It is found that the growth of the specific energies in fig. 6(a) and fig. 7(a) is characterised by three regions both for the trefoil and the figureeight. At very low temperatures the energy increase is practically zero. The reason is that the temperatures are too low to allow contacts between the monomers. When the energy of the system grows – let’s remember that we are working in energy units, so that we can talk about temperatures and energies interchangeably – the “first energy state” can be reached. This consists in the case in which at least two monomers are in contact. For instance, we expect that when the chain has length , the first energy state occurs in figs. 6 and 7 at some “critical” temperature , while when this temperature should be given by . These values of slightly increase with increasing topological complexity. In general, this proportionality temperature/length is respected in figs. 6 and 7. After the first energy state is reached, the specific energy grows rapidly as a function of the temperature until saturation is reached after a given temperature and the energy increase becomes moderate. The division into three regions can be expected by looking at the Boltzmann factor appearing in the expression of the energy given in Eq. (4):

i) In the range , the factor goes rapidly to , because low temperatures cause few states to be occupied in this situation, so we get .

ii) When , the factor changes sensitively depending on . In this regime more and more states of the given system are getting excited through heat absorption and increases with the increase of the temperature as expected.

iii) Once , increases slower due to the fact that most states have been already occupied in the range . Just few states are left that can be excited in this situation. One could notice that the factor in Eq. (4) increases slowly at high temperatures. When , we have .
It can also be found that in fig. 6(b) and fig. 7(b), the heat capacity has a peak in the range of temperature . The appearance of this peak together with the sigmoidal behavior of the specific energy points out at an ongoing twostate phase transition according to cooper (). The presence of these pseudotransitions related to the lattice geometry have been clarified in ref. wustlandau (); vogel ().
We would also like to show that the heat capacity and energy are affected not only by the temperature, but also by the polymer length and its topological conformation. First of all, the specific energy of the chain decreases with increasing values of . For instance, the quantity for is larger than that of longer chains as it is illustrated in fig. 6 and fig. 7. This is because longer polymers in the presence of shortrange repulsive interactions prefer states which are diluted comparing with the short ones. Probably for this reason, in proteins knots are localized in regions with a very small number of residues dna (). The total internal energy, instead, decreases with decreasing polymer lengths. This is due to the fact that longer polymers can accommodate a larger number of contacts. Let us note that the peak temperatures of the heat capacity are shifting with the length . This is expected since we are computing the specific heat capacity. If we computed the total heat capacity of the system, these peaks would occur at the same temperature.
Concerning the topological effects, we observe that both the energy and
the heat capacity increase with increasing knot complexity. This can
be seen by comparing the plots for various knots of the same length in
figs. 5, 6 and 7.
On the other side, the topological effects are stronger
when polymers are shorter. This is shown in
Table 1 in which the energy differences for
polymers with the same length but different topology are listed.
In the table is the difference between the
unknot and the trefoil defined by
, while
denotes the difference between the trefoil and figureeight. When the
length ranges from to , varies from
to and from to
at . The heat capacity exhibits a similar
behavior. It is thus
pretty clear that
the topological effects become less important when the length
is increasing. This is similar to the so called ’size effect’ in
nanomaterials yani ().
Topological
properties will have a lesser influence in the case of diluted
monomers. We expect that topological effects will disappear if the
polymer is long
enough.
unknot  trefoil  figureeight  

100  0.174  0.247  0.279  0.073  0.032 
200  0.170  0.197  0.210  0.027  0.013 
300  0.167  0.182  0.186  0.015  0.004 
To get a feeling of the size of the involved quantities, we can make the concrete example of the proteins which, as mentioned before, can be found in different topological conformations like the trefoil and the figureeight. These proteins have a length of about 250 residues dna () and their temperature is around (since this is approximately the temperature of the human body). The parameter for a protein is about ljwell (). Based on those data, for one mole, we can roughly estimate the range of corresponding temperatures in our units, which is . This temperature is located in the region .
Before concluding this section, we would like to stress that the used method for preserving the knot topology after a pivot move is very time efficient and can be applied to any knot configuration independently of its complexity. For example, in fig. 8 we provide the data for the specific energy and heat capacity in the case of the knot . As it is possible to check, also for this knot the same conclusions drawn before for are valid. In particular, the topological effects become less important with increasing polymer length and it is confirmed that both the energy and the heat capacity increase with increasing knot complexity.
Iv Conclusions
We have studied the topological effects on the thermal properties of several knot configurations with the help of Monte Carlo simulations based on the WangLandau algorithm and the pivot method. The topology of the knots, that can change after a pivot move, has been preserved by explicitly checking if during the move some parts of the trajectory of the knot have crossed themselves at some points.
The specific energy and the heat capacity of the trefoil, the figureeight and the knot have been calculated for various lengths ranging from to . As a test of our methodology, these quantities have been derived also in the known case of the unknot. In particular, we have been able to reproduce the results of pnring (). Our investigations show that the specific energy and heat capacity for a polymer knot whose monomers are subjected to shortrange repulsive interactions are depending on the length and the topological conformation of the polymer knot. This could be expected because in the presence of repulsive interactions the polymer will prefer “diluted” conformations in which the monomer density is lower. We may thus suggest that the dilution of polymers in a good solvent should be effective in reducing the topological effects. Moreover, for the growth of the specific energy with respect to the temperature it is possible to distinguish three different regimes that have been explained in details in Section II. The zero energy regime appears at extremely low temperatures and it is characterised by the fact that the energy of the thermal fluctuations is below the threshold of the repulsive energy barrier, so that the monomers have not enough energy to get near. After this threshold has been passed, the number of contacts between the monomers starts to increase rapidly until at the end a saturation phenomenon occurs and the energy of the system increases only slowly with increasing temperatures. The plot of the specific heat capacity shows a pseudophase transition which is related to the lattice geometry, as explained in refs. wustlandau () and vogel ().
Our simulations can be generalised under several aspects.
First of all, we limited ourselves to simple short range
interactions, but there are no obstacles to extend our
procedure to describe more realistic polymer systems.
Moreover, we have limited ourselves to pivot moves involving a short
number of segments. No difference between the case and
have been detected, so that all displayed diagrams are related only to the
case . The independence of up to has been confirmed by
preliminary calculations performed using another method
to preserve the topology based on knot invariants.
We hope to be able to report on these new
developments very soon zhaoferrari ().
Acknowledgements.
We are indebted with V. G. Rostiashvili for helpful discussions and suggestions that stimulated the development of the method for distinguishing the topology of the different knot configurations presented in this work. We wish also to thank heartily J. Paturej and T. A. Vilgis for fruitful discussions. The support of the Polish National Center of Science, scientific project No. N N202 326240, is gratefully acknowledged. The simulations reported in this work were performed in part using the HPC cluster HAL9000 of the Computing Centre of the Faculty of Mathematics and Physics at the University of Szczecin.References
 A. Yu. Grosberg Phys.Usp. 40, 12 (1997).
 W. R. Taylor, Nature (London) 406, 916 (2000).
 V. Katritch, J. Bednar, D. Michoud, R. G. Scharein, J. Dubochet, A. Stasiak, Nature 384, 142 (1996).
 V. Katritch, W. K. Olson, P. Pieranski, J. Dubochet and A. Stasiak, Nature 388, 148 (1997).
 M. A. Krasnow, A. Stasiak, S. J. Spengler, F. Dean, T. Koller and N. R. Cozzarelli, Nature 304, 559 (1983).
 B. Laurie, V. Katritch, J. Dubochet and A. Stasiak, Biophys. Jour. 74, 2815 (1998).
 J. I. Sułkowksa, P. Sułkowksa, P. Szymczak and M. Cieplak, Phys. Rev. Lett. 100, 058106 (2008).
 Z. Liu, E. L. Zechiedrich, and H. S. Chan, Biophys. J. 90, 2344 (2006).
 S. A. Wasserman and N. R. Cozzarelli, Science 232, 951 (1986).
 D. W. Sumners, âKnot theory and DNA,â in New Scientific Applications of Geometry and Topology, edited by D. W. Sumners, Proceedings of Symposia in Applied Mathematics, Vol. 45, ÍAmerican Mathematical Society, Providence, RI, 1992, 39.
 A. V. Vologodski, Ì A. V. Lukashin, M. D. FrankKamenetski Ì and V. V. Anshelevich, Zh. Eksp. Teor. Fiz. 66, 2153 (1974); Sov. Phys. JETP 39, 1059 (1975); M. D. FrankKamenetskii, A. V. Lukashin and A. V. Vologodskii, Nature (London) 258, 398 (1975).
 E. Orlandini, S. G. Whittington, Rev. Mod. Phys. 79, 611 (2007); C. Micheletti, D. Marenduzzo, and E. Orlandini, Phys. Reports 504, 1 (2011).
 T. Vettorel, A. Yu. Grosberg and K. Kremer, Phys. Biol. 6, 025013 (2009).
 P. Virnau, Y. Kantor and M. Kardar, J. Am. Chem. Soc. 127 (43), 15102 (2005).
 P. Pierański, S. Przybył and A. Stasiak, EPJ E 6 (2), 123 (2001).
 R. Metzler, A. Hanke, P. G. Dommersnes, Y. Kantor and M. Kardar, Phys. Rev. Lett. 88, 188101 (2002).
 K. Koniaris and M. Muthukumar, Phys. Rev. Lett. 66, 2211 (1991).
 R. Scharein et. al., J. Phys. A: Math. Gen. 43, 475006 (2009).
 N. Madras, A. Orlistsky and L. A. Shepp, Journal of Statistical Physics 58, 159 (1990).
 F. Wang and D. P. Laudau, Phys. Rev. Lett. 86, 2050 (2001).
 C. Vree and S. G. Mayr, New Journal of Physics 12, 023001 (2010).
 T. Kennedy, Jour. of Stat. Phys., 106 (34), 407429 (2002).
 P. N. VorontsovVelyaminov, N. A. Volkov and A. A. Yurchenko, J. Phys. A: Math. Gen. 37, 15731588 (2004).
 N. A. Volkov, A. A. Yurchenko, A. P. Lyubartsev and P. N. VorontsovVelyaminov, Macromol. Theory Simul. 14, 491504 (2005).
 A. Cooper, Thermodynamics of protein folding and stability, in Protein: A Comprehensive Treatise, Vol. 2, A. G. Stamford (Ed.), JAI Press. Inc., 217 (1999).
 T. Wüst and D. P. Landau, Phys. Rev. Lett. 102, 178101 (2009).
 T. Vogel, M. Bachmann and W. Janke, Phys. Rev. E. 76, 061803 (2007).
 Y.N. Zhao, S.X. Qu and Ke Xia, Jour. Appl. Phys. 110, 064312 (2011).
 M. Schlosshauer and D. Baker, Protein Science 13, 1660 (2004).
 Y. Zhao and F. Ferrari, work in progress.