BACKGROUND
plasmid dna  is a family of dna molecules widely used in life sciences, more specifically in gene therapy research.
these molecules are produced inside host cells in a supercoiled conformation , which is the desired conformation for therapy purposes.
however, such molecules can lose their original conformation in the production and purification processes, assuming more relaxed or even linear conformations, owing to thermodynamic changes .
one of the main challenges for researchers is to find optimal thermodynamic conditions for plasmid dna therapeutic application without losing its supercoiled conformation or, at least, minimizing the occurrence of relaxed or open dna molecules.
for many years, computational methods based on laboratory experimental data have been proposed to model and simulate the dynamic behavior and conformational changes in pdna molecules under certain conditions.
the monte carlo  method has generally been accepted as a reliable tool for simulation purposes, and is seen as the standard.
this iterative method tries to minimize the elastic energy of the molecule in each iteration step of the simulation process, testing the probability of acceptance of each new trial.
the goal is to make the molecule converge to an equilibrium state after performing as few iterations as possible, i.e., maximizing the acceptance ratio of the trials without compromising the effectiveness and reliability of the simulation.
to simplify the simulation process, each plasmid dna molecule is reduced to a linear skeleton  with equal sized segments that represents the topological conformation of the molecule.
random deformations are then applied to this skeleton, generating new trial conformations, which are either accepted or rejected.
interestingly, the essence of the method used to randomly generate each new trial, referred to as the standard crankshaft move, has remained the same for many years, with its origins dating back to the early 1960s , more specifically in the context of lattice polymer chains.
this move was later adapted for simulation of flexible molecules like dna using mc methods.
however, the standard crankshaft move has a very low acceptance ratio of trials, i.e., many trials are rejected.
moreover, it can present very unnatural behavior as it features very sudden motions along large portions of the molecule.
to enhance the efficiency of mc moves, biasing was found to be a solution  <cit> .
however, as earl and deem noted  <cit> , biasing a deformation move implies that the probability of moving from one state to another is no longer symmetric; consequently, the acceptance rule used must be altered to maintain the detailed balance.
in this paper, we present a new unbiased move for plasmid dna, whose skeleton is a closed polyline.
this move not only preserves the size of each segment and its connectivity, but is also very effective in maximizing the acceptance ratio of the trials and stabilizing the molecule, thereby allowing steady, gradual temperature changes during the simulation.
our method also generates natural and realistic animations that can be used in real-time simulation and visualization.
related work
in this section, we briefly review the mc methods in computational biology and chemistry, as well as the generative methods for dna conformations that form the core of these mc methods.
monte carlo simulations
the mc simulation method is one of the most important methods used in dna simulations.
this method, which was originally presented by metropolis et al.
<cit> , generates dna conformations combining energy calculations, random conformational changes, and statistics.
frank-kamenetskii et al.
<cit> , vologodskii et al.
<cit>  and lebret  <cit>  were the first to use an mc method to present numerical results of the probability of the occurrence of knots on pdna.
frank-kamenetskii and vologodskii also presented valuable information on dna torsional rigidity  <cit> .
a few years later, vologodskii et al.
used mc simulations to study the conformational and thermodynamic properties of dna molecules with physiological levels of supercoiling  <cit> .
vologodskii also included a chapter on “monte carlo simulation of dna topological properties” in the book “topology in molecular biology”  <cit> , and with rybenkov, they reviewed how conformational properties of dna catenanes can be studied using mc simulations  <cit> .
gebe et al.
<cit>  presented an mc algorithm to simulate supercoiling free energies in unknotted and trefoil knotted inextensible circular chains with finite twisting and bending rigidity, while marko et al.
<cit>  made use of mc simulations to study the relationship between the amount of twisting in dna molecules and its supercoiling.
kundu et al.
used an mc algorithm to explain denaturation characteristics in a supercoiled plasmid and calculate the probability of denaturation for each base pair at different supercoiling degrees  <cit> .
in their work on the relationship between knots and supercoiling, cozzareli et al.
used an mc simulation procedure to generate an equilibrium set of conformations  <cit> .
mc simulations have also contributed to the understanding of the interplay between base-pair stacking interaction and permanent hydrogen-bond constraints in supercoiled dna elasticity  <cit> .
based on the fact that atomic force microscopy has generated images of supercoiled dna confined to a surface, which affects conformational properties such as twist and writhe, fujimoto and schurr modified an existing program, developed to perform mc simulations of supercoiled dna in solution, flattening the dna to simulate the effect of deposition on a surface  <cit> .
fujimoto and schurr also presented a method to estimate torsional rigidities of weakly strained dna  <cit> .
burnier et al.
used mc calculations to identify a mechanism by which topoisomerases can keep the knotting level low  <cit> .
more recently, olson et al., in their paper “how stiff is dna?”, used mc simulations to understand the behavior of a long, double-helical polymer in the tight confines of a cell and in the design of novel nanomaterials and molecular devices  <cit> .
generative methods of dna conformations
it has generally been accepted that supercoiled, i.e., the self twisting of the double stranded molecule over itself, is the desired conformation for pdna molecules  <cit> .
thus, it is necessary to measure the supercoiling of a given molecule.
one of the most important quantitative measures of closed circular dna supercoiling is the linking number  of the two dna strands, which is an integer corresponding to the number of double-helical turns of the molecule.
there are several methods for calculating lk, with one of the most widely used involving the computation of two very important geometrical properties of closed circular dna molecules: twist  and writhe .
tw features the coiling of the two dna strands around the axis of the helix, while wr is a measure of the coiling of the helix axis in space  <cit> .
thus, the main result is:
in our implementation, we used klenin and langowski’s computation method  to calculate wr
<cit> .
owing to the nature of three-dimensional closed polylines, knots can occur in some pdna conformations.
this is not a desirable feature, i.e., each closed circular dna molecule must remain unknotted during the simulation process, keeping its original topology even if supercoiling occurs.
knot detection methods must be used during simulation to reject possible knotted conformations.
we adopted harris-harvey’s knot detection algorithm, which uses the alexander polynomial to detect the existence of knots  <cit> .
this algorithm is based on the predicate that if two knots have different alexander polynomials, then the knots are topologically distinct.
thus, because the alexander polynomial of an unknotted closed circular dna molecule is equal to one, all conformations for which this polynomial does not equal one must be rejected during the simulation.
each trial conformation must be generated in such a way that the size and connectivity of each segment of the dna chain do not change.
a major deformation method used to displace vertices of the dna chain was introduced by klenin et al.
<cit> .
this method, which is just a biased crankshaft move, starts by randomly choosing two vertices vm and vn.
then, all the vertices  are rotated a randomly selected angle θ around the axis defined by the line connecting vm and vn, as shown in figure  <dig>  furthermore, following klenin et al.
<cit> , the value of θ is uniformly distributed over a certain interval, and must be continuously adjusted during the simulation to guarantee that about half the steps are accepted.figure  <dig>
crankshaft motion.
to increase the acceptance ratio of the simulation trials another type of motion has been proposed in the literature.
this improvement performs a sub-chain translation, which is usually referred to as reptation motion
<cit> , and is illustrated in figure  <dig>  first, two vertices vi and vj are randomly chosen.
then, the sub-chain between vi and vj is translated by one segment length along the chain contour.
the segment that was immediately after vj is also translated to fill the gap between vi and vi+ <dig>  this motion suggests movement analogous to a snake slithering and, hence, the name reptation motion.
other types of motion can be adopted if the metropolis microscopic reversibility requirement is satisfied, i.e., if the probability of each trial conformation is the same as that of the reverse movement  <cit> .figure  <dig>
reptation motion.
visualization of dna conformational changes over time is also important as part of the entire simulation of dna behavior.
this is usually performed only when the entire simulation procedure ends and is typically done by assembling the dna atoms along the dna axis.
interestingly, a more efficient dna assembly algorithm was presented to allow the visualization of all the steps of the simulation procedure in real-time  <cit> .
in this method, each dna nucleotide is represented by a three-dimensional building block, allowing the assembly of the entire molecule faster, but in a realistic way.
in geometric terms, each of the four building blocks featuring dna nucleotides is a gaussian isosurface, which was previously generated by an algorithm that triangulates molecular surfaces  <cit> .
methods
the deformation algorithm presented in this paper uses a linear skeleton  with equal sized segments, henceforth called the dna skeleton.
before introducing the core of the method itself, we explain how the dna skeleton can be created for use by the deformation algorithm.
initial conformation of the dna skeleton
the dna skeleton can assume any closed unknotted conformation.
the simplest of these conformations is a completely relaxed circular conformation.
besides, the length of each segment of the dna skeleton corresponds approximately to  <dig> base pairs of the double helix  <cit> .
that said, the first step of the algorithm is to determine the number of segments of the dna skeleton ensuring around  <dig> base pairs per segment.
assuming that the dna has a sequence of n base pairs, we want to find an integer s denoting the number of segments of the dna skeleton.
we define two integer parameters min and max, respectively, as the minimum and maximum numbers of base pairs that are admissible per segment, such that min < 30 <max.
then, for each integer value i,min ≤ i ≤ max, we calculate the corresponding si=round.
finally, we adopt s=si as the number of segments of the dna skeleton that minimizes |n-|.
once we have the number of segments s, we just have to build a regular polygon with s sides inscribed in a circle.
from the number of base pairs, we can infer the approximate perimeter of the circle, as well as the corresponding radius r, from which we obtain the first vertex of the skeleton at p0=.
then, we apply s successive rotations to p <dig> about the origin to obtain all the vertices of the dna skeleton of the initial relaxed conformation; the rotation angle is given by α=2π/s.
note that, although the initial conformation is circular, the methods for dna assembly and deformation apply to any initial conformation.
skeleton deformation algorithm
assuming we have a three-dimensional closed polyline pk representing the dna skeleton, we need to deform this polyline to obtain a new polyline pk+ <dig> with the same number of equal sized segments as pk, but without loss of its connectivity.
let s be the number of segments of pk and {vi}, i= <dig> …,s- <dig>  its vertices.
we choose a random vertex vm, 0≤m≤s- <dig> as our mobile vertex, i.e., the vertex with the most motion freedom in the current trial conformation.
any movement of the mobile vertexvm implies movement of its closest neighbors vm- <dig> and vm+ <dig>  called semi-mobile vertices .
the remaining neighbors vm- <dig> and vm+ <dig> are fixed vertices because they do not move in the deformation.
thus, in each deformation step, only three vertices will be displaced: vm, vm- <dig> and vm+ <dig> figure  <dig>
mobile vertex
v
m
can be displaced randomly in the intersection of three spheres,
n
m
,
s
m- <dig>  and s
m+2
.
however, vm cannot be freely displaced .
in the first instance, vm moves within the sphere nm centered at vm with radius r= <dig> Δ, where Δ= <dig> Å is the distance between two consecutive base pairs.
more specifically, the new position of vm is found randomly within the region resulting from the intersection of the three spheres, nm, sm- <dig> and sm+ <dig>  the latter two spheres with radius 2l are centered on the fixed vertices vm- <dig> and vm+ <dig>  respectively, where l is the length of each segment of the dna skeleton.
note that the optimal value  <dig> Δ for r was found experimentally and based on the rate of successful trial conformations in the first attempt, about  <dig> percent, though this rate remained high for a value of r up to  <dig> Δ. the small radius r of sphere nm ensures that the transition from pk to pk+ <dig> occurs without noticeable jumps.
as noted above, the new position of vm was found randomly within the region nm∩sm-2∩sm+ <dig> , but no explanation of this random procedure was given.
in fact, to calculate the new position of vm, we first convert its cartesian coordinates  to spherical coordinates  relative to vm- <dig>  where d is the distance between vm- <dig> and vm.
next, we randomly generate a new position for vm as , where Δd∈ and Δθ,Δϕ∈.
it is clear that the displacement of the flanking vertices vm- <dig> and vm+ <dig> depends on the previous movement of vm.
here we focus on the computation of the new position of vm- <dig> since the new position of vm+ <dig> can be calculated similarly.
for this purpose, we also convert the cartesian coordinates of vm- <dig> to spherical coordinates  relative to vm- <dig>  where l is the radius of the three spheres sm, sm- <dig>  and sm- <dig> centered on vm, vm- <dig>  and vm- <dig>  respectively.
moving vm- <dig> to a new position must be done in such a way that its distance l to vm- <dig> and vm remains unchanged.
in other words, the new vm- <dig> must lie on the circumference resulting from the intersection of the two surfaces bounding sm and sm- <dig> .
if Δd= <dig>  the new position of vm- <dig> relative to vm- <dig> is given by ; otherwise, the new location of vm- <dig> is , where Δψ is the angle of the angular motion of vm- <dig> on sm- <dig> resulting from the translational displacement Δd of vm along the line defined by vm and vm- <dig> .
we compute Δψ by rearranging the equation that describes the reciprocal motion of the piston with respect to the crank angle as follows :
<dig>  figure  <dig>
displacement of vertices v
m
and
v
m-1
.
translational piston move of vertex v
m
translates into a rotational move of v
m-1
.
note that applying the translational displacement Δd to vm before the rotational motions Δθ and Δϕ means that Δθ= <dig> in eq.
; otherwise, Δθ≠ <dig>  in summary, moving vm implies a translational and two rotational motions relative to vm- <dig> expressed in spherical coordinates.
this causes vm- <dig> to rotate accordingly on the sphere centered at vm- <dig>  with part of this rotational motion determined by the translational displacement Δd of vm.
these types of moves  satisfy the principle of microscopic reversibility  <cit> , although this is not critical for our purposes because the simulation procedure is only used to locate energy minima.
as noted by mauri in  <cit> , for a conservative n-body system, as in the case of a dna molecule, microreversibility stems from the invariance of the equations of motion with respect to time reversal, i.e., every microscopic motion reversing all particle velocities also results in a solution.
this leads to the so-called principle of detailed balance  <cit> , which states that under stationary conditions  each possible transition from one conformation to another balances itself with the reversed transition in time.
in other words, the probability of obtaining trial conformation pk+ <dig> if the current conformation is pk must be equal to the probability of obtaining trial conformation pk if the current conformation is pk+1
<cit> .
with this in mind, and having calculated the constrained position of vm- <dig> as a consequence of the move of vm, we need to determine its new position after rotating it randomly about the axis defined by vm- <dig> and vm.
it is clear that the old and new locations of vm- <dig> lie on the circumference resulting from the intersection of spheres sm- <dig> and sm.
likewise, we find the new position of vm+ <dig> after rotating it randomly about the axis defined by vm+ <dig> and vm.
interestingly, these two rotations can be seen as two particular crankshaft rotations.
finally, it is worth noting that the deformation algorithm described above can also be used in other biochemical systems such as internal coordinate models of cyclic peptides, as well as in some mechanical problems related to articulated arms and chain moves.
in fact, this algorithm can be used to randomly deform any regular polygon  in two or three dimensions with guaranteed preservation of connectivity and the length of segments.
dna assembly algorithm
for a realistic visualization of closed circular dna simulations in real-time, we combine the new deformation algorithm described above with the dna assembly algorithm introduced by raposo and gomes  <cit> .
this dna assembly algorithm uses four three-dimensional building blocks representing dna nucleotides , namely, adenine , cytosine , thymine , and guanine .
each building block is a pre-triangulated isosurface generated by a triangulation algorithm for molecules  <cit> .
the assembly procedure for nucleotides can be thought of as the operation of wrapping helicoidal dna backbones around cylinders along the dna skeleton.
the algorithm iterates over n base pairs of the dna, assembling a single base pair bibi per iteration.
each iteration comprises four distinct stages:
generation of geometric instances for nucleotidesbiandbi.
considering that there are only four possible base pairs, given a nucleobase ni of a dna strand, two geometric instances of nucleotides must be generated, the first for the building block bi and the second for the corresponding building block bi.
positioning of the base pairbibion planez= <dig>  base pair bibi is positioned as if it was the first base-pair of the dna molecule, that is, it is placed perpendicular to the segment that starts at the origin in such a way that its center coincides with the origin.
alignment of base pairbibiwith the plane perpendicular to segment i. note that this alignment involves a rotation about the origin of the coordinate system.
translation of base pairbibito the plane perpendicular to segment i. finally, because all geometric transformations  are performed around the origin, building blocks bi and bi must be displaced to their correct positions relative to the midpoint of the corresponding segment i of the dna axis.
dna building blocks:  adenine,  cytosine,  guanine, and  thymine.
it is important to note that this dna assembly algorithm does not take into account the sharp kinks that may occur at the junctions of the conformation segments, as shown in figure  <dig>  nevertheless, a possible solution to this problem is the smoothing procedure proposed by kummerle and pomplun  <cit> .
for detailed information about the dna assembly algorithm, the reader is referred to  <cit> .figure  <dig>
detail of dna segments assembled around a random piece of the skeleton.
monte carlo simulation
mc simulations are iterative methods based on the concept of elastic energy of closed circular dna and on stochastic parameters and calculations aimed at converging to the energetic and thermodynamic equilibrium of the molecule.
the main principle is to perform random dna deformations and check whether the resulting new conformations should be accepted or rejected according to energy changes and acceptance probability.
more specifically, a random deformation of the dna is accepted if it reduces the elastic energy of the molecule or has some probability of occurring.
in the experiments and results presented in this paper, we used the same mc simulation method and parameters as those used in  <cit> , where elastic energye is calculated as
<dig>
here eb is the bending energy given by
<dig>
where kb is the boltzmann constant, t is the temperature, α= <dig>  is the bending constant, and Θi is the angular displacement between the directions of segments i and i+ <dig>  torsional energyet is given by
<dig>
where c denotes a constant parameter known as the torsional rigidity, l is the total length of the chain, and wr is the writhe of the skeleton.
the linking number difference Δlk in  is the difference between the linking number lk of the dna molecule and that, lk <dig>  of its relaxed dna conformation
<dig>
where - <dig> ≤σ≤- <dig>  is the superhelix density parameter  <cit> .
for calculating writhe wr, we used the method proposed in  <cit> , more specifically, method 2b  <cit> .
this method is based on the principle that writhe can be calculated as the difference between linking number lk and twist tw:
<dig>
this method for computing wr uses an auxiliary chain close enough to the dna skeleton, and with as many segments si as the dna skeleton.
then, considering that ri is the initial point of segment si, we can find the directional writhe as follows:
<dig>
where
<dig>
for the computation of , we must check whether segments si and sj cross, i.e., whether their projections π and, π onto plane z= <dig> intersect.
we used lamothe’s algorithm to check whether the projections of these segments intersect  <cit> .
in turn, and following klenin and langowski  <cit> , twist is given as
<dig>
where pi denotes the vector, normal to both si- <dig> and si, z denotes the z-th component of vector pi vector, and
<dig>
where u is the unit vector in the z axis direction.
then, using the results of  and , we get the final writhe number:
<dig>
once we know how to perform the necessary energy calculations, we can apply the mc method to each iteration of the simulation to obtain a new dna conformation from a random deformation of the dna skeleton.
then, we calculate the energy ei+ <dig> of the new candidate conformation, and compare it with the energy, ei, of the previous conformation.
the new conformation i+ <dig> is accepted if ei+1<ei or
<dig>
where tm is the temperature of the experiment and ρ is a random value between  <dig> and  <dig>  <cit> .
knots detection
it is important to note that knots can occur when random deformations are applied to dna conformations.
because this is not desirable, i.e., dna supercoiling must occur without generating knots, we must check for the existence of knots and reject the deformation if we find one or more knots.
to optimize the performance of the method, this checking procedure is done before the mc acceptance test, avoiding unnecessary energy calculations.
for knot detection we used the method of harris and harvey  <cit> .
in this method, based on the principle that two knots are topologically distinct if they have distinct alexander polynomials, the dna skeleton is converted to a knot data structure, and its alexander polynomial is computed and compared with the alexander polynomial of the circle, which is a trivial knot.
if these two polynomials are different, the dna skeleton contains at least one non-trivial knot.
RESULTS
to evaluate the effectiveness and performance of our deformation method when applied in mc simulations, we performed a set of experiments comparing our method with two types of dna chain moves, namely, the standard crankshaft move and the biased crankshaft move.
the standard crankshaft move is a randomly chosen move.
in fact, the ends, vm and vn, of each sub-chain are randomly chosen, as is the case with the rotating angle θ of the sub-chain around the line that passes through its ends .
that is, the standard crankshaft move does not adjust the size of the sub-chain nor the rotation angle in any way.
on the other hand, as in the deformation method introduced by klenin et al.
<cit> , the biased crankshaft move used here adjusts only the rotation angle.
recall that this type of biased moves is a way of enhancing the efficiency of mc moves, because it allows us to choose moves with a higher acceptance ratio  <cit> .
through these experiments we aim to demonstrate that our method generates a smoother and more controlled deformation, which leads to more consistent and even faster convergence to molecular energy equilibrium.
experimental setup
three experiments  were performed to compare the proposed method with two classic methods, namely, the standard crankshaft move and biased crankshaft move.
we used a setup based on kummerle and pomplun’s work  <cit>  for the puc <dig> plasmid dna molecule.
all the three experiments were performed using the same closed circular dna sequence with  <dig> base pairs  and exactly the same conditions and mc simulation parameters for the three monte carlo moves under comparison, that is: kb= <dig> -23; α= <dig> ; c=3×10-28; and σ=- <dig> .
however, we performed experiment a at a constant temperature of  <dig> k, while experiments b and c were performed progressively reducing the temperature from  <dig> k to  <dig> k.
finally, it should be noted that the mc simulations were performed on an 64-bit windows  <dig> laptop computer with an intel i <dig>  <dig>  ghz cpu,  <dig> gb ram and an nvidia geforce gt  <dig> mx  <dig> gb graphics card.
experiment a: puc <dig> with constant temperature
in experiment a, we performed an mc simulation with  <dig>  steps using the puc <dig> closed circular dna sequence at a constant temperature of  <dig> k as in  <cit> .
this experiment was replicated using:  the standard crankshaft move,  the biased crankshaft move as described in  <cit> , and  the proposed method.
in the particular case of the biased crankshaft move, after an exhaustive optimization procedure with more than  <dig>  steps, we came to the conclusion that the crankshaft rotation angle should initially be in the range  , as shown in figure  <dig>  furthermore, as expected, this angle decreases over time figure  <dig> figure  <dig>
optimization of initial angle range for  <dig>  steps.
experiment a: crankshaft rotation angle for  <dig>  steps at a constant temperature of  <dig> k.
we used two measures to compare the efficiency of the three methods:  elastic energy equilibrium and  acceptance ratio of trials.
the graphs of elastic energy for the three methods are shown in figure  <dig>  where we can see that the average elastic energy for each method remains approximately the same over time.
nevertheless, on average, the elastic energy of the proposed method is slightly higher than that of the standard crankshaft method, which in turn is higher than the energy associated with the biased crankshaft method.figure  <dig>
experiment a: elastic energy for  <dig>  steps at a constant temperature of  <dig> k.
on the other hand, the acceptance ratio of trials was evaluated for each slice of  <dig>  steps from a total of  <dig>  steps , according to the acceptance condition .
the acceptance ratio was steadily higher for the proposed method, always remaining above  <dig>   accepted steps for each slice of  <dig>  steps, i.e., an average acceptance ratio around 45%, and achieving higher ratios around 50% in the second half of the experiment.
on the contrary, the acceptance ratio for the standard crankshaft move was always under 30%, and even lower in the first  <dig>  steps of the experiment.
with respect to the biased crankshaft move, the average acceptance ratio was slightly above 30%, but far below the results obtained using the proposed method.
this indicates that our new method generates trials with much higher probabilities of being accepted under the mc simulation conditions at a constant temperature, i.e., it minimizes the number of trial rejections, and avoids useless computations.figure  <dig>
experiment a: acceptance trials for slices of  <dig>  steps from a total of  <dig>  steps at a constant temperature of  <dig> k.
during experiment a, we also noted that the crankshaft moves generated a few hundred conformations that were rejected owing to the existence of knots.
on the contrary, our new method did not produce any knots at any time during the  <dig>  simulation steps because the deformation is done smoothly and without conformational jumps.
this concurs with the fact that, despite dna molecules in living cells being long and compactly coiled, they rarely get knotted  <cit> , which suggests that supercoiling inhibits dna knotting.
experiment b: puc <dig> with variable temperature
experiment b also involved  <dig>  mc steps.
this experiment was also replicated for each of three methods analyzed in this paper.
in the particular case of the biased crankshaft move, once again, after an exhaustive optimization procedure with over  <dig>  steps, we concluded that the initial crankshaft rotation angle should be in the range   because, in experiment a, the temperature is not constant.
adjustments to this rotation angle range during the experiment are shown in figure  <dig> more specifically, the temperature decreases with energy, i.e., if the average elastic energy of a  <dig> -step interval is higher than that of the previous slice of  <dig>  steps multiplied by a  <dig>  factor, the temperature is also multiplied by a  <dig>  factor.
in fact, the temperature decreased progressively from  <dig> k to  <dig> k. as expected, the closer the method converges to the energy equilibrium, the greater is the decrease in temperature.
as in experiment a, the acceptance ratio of the proposed deformation method was always higher than those of the classic deformation methods .figure  <dig>
experiment b: crankshaft rotation angle for  <dig>  steps with temperature varying between  <dig> k and  <dig> k.
experiment b : acceptance trials for slices of  <dig>  steps from a total of  <dig>  steps with temperature varying between  <dig> k and  <dig> k.
it was a somewhat surprising to observe how the energy decreased during the simulation.
as shown in figure  <dig> , when using the proposed method, the elastic energy of the molecule converged sooner and more consistently to equilibrium.
as shown, we achieved energy equilibrium after approximately  <dig>  mc steps, while the crankshaft moves only stabilized after  <dig>  steps.
besides, the standard crankshaft move generated a number of very slight energy jumps, i.e., the energy did not decay as consistently as in the proposed method.
however, the energy level at equilibrium was the same for all three methods, approximately  <dig> ×10- <dig> no less meaningful was the temperature decay during this experiment.
as presented in figure  <dig> , when using the proposed method, the mc temperature decreased more rapidly and in a more consistent way, i.e., the graph for the new method is much smoother with the advantage of reaching the equilibrium temperature sooner.
figure  <dig> shows the final result of experiment b.figure  <dig>
experiment b:  elastic energy during a  <dig> -step experiment with temperature varying between  <dig> k and  <dig> k;  temperature decaying during a  <dig> -step experiment.
experiment b: puc <dig> after  <dig>  mc steps with temperature varying between  <dig> k and  <dig> k.
in summary, we can say that the proposed deformation method requires fewer simulation steps to achieve energy equilibrium, largely owing to its high acceptance ratio.
experiment c: average displacement
in experiment c, we set out to measure the amount of deformation of plasmid dna caused by each type of mc move.
this was accomplished by computing the average displacement of the dna skeleton vertices for each accepted trial.
in this experiment, we only considered the standard crankshaft move and the move proposed in this paper.
taking into account vertices vi,vi+ <dig> …,vi+n that are displaced during a simulation trial, we determined the distances di,…,di+n between the new positions and the previous positions of these vertices, and straightforwardly computed the average displacement given by /.
finally, we considered the accumulated displacement for a slice of steps as the sum of the average displacements of the accepted trials of that slice.more specifically, we performed a  <dig> -step simulation for each of the methods, namely, the standard crankshaft move and the proposed move.
as shown in figure  <dig> , the new move generates much smaller average displacements than the standard crankshaft move.
besides, from figure  <dig>  we can see that the new move generates displacements right from the start of the simulation, whereas the standard crankshaft move starts to produce displacements later.
this can be explained by the high acceptance ratio of the new method, as well as its more steady deformations.moreover, from figure  <dig>  we conclude that smaller displacements in each trial do not mean there will be smaller accumulated displacements.
the accumulated displacements of the new method form a logarithmic curve, while the curve of the standard crankshaft move is clearly exponential ).
this means that in the new method, the closer we get to the point of energy equilibrium, the shorter is the displacement toward a stable conformation, i.e., the average displacement for each accepted trial converges to zero as the conformation converges to the equilibrium.
this is not the case in the standard crankshaft move, as illustrated by the accentuated variations in average displacement of the trials in figure  <dig> .figure  <dig>
experiment c:  average displacement per  <dig> steps for puc <dig>  <dig>  steps;  average accumulated displacement for puc <dig>  <dig>  steps.
discussion
as mentioned above, we used the same molecule , the same conditions/parameters, and the same mc simulation method in experiments a, b, and c. for comparison purposes, each experiment was performed using three different deformation methods:  the standard crankshaft move,  the biased crankshaft move , and  the proposed method.
as expected, the acceptance ratio of trials for the proposed method is higher than that for either of the crankshaft moves.
the acceptance ratio of the new method is almost always greater than 40%, and even reaches more than 60% at certain times.
moreover, unlike the crankshaft moves, the acceptance ratio for the proposed method is very high from the very first steps of the simulation.
more specifically, in a scenario with decreasing temperature, we obtained an acceptance ratio of more than 60% for the proposed method compared with 5% for the standard crankshaft move and 20% for the biased crankshaft move in the first  <dig>  mc steps.
with respect to elastic energy, the experiments also show that as the temperature decreases the new move achieves better performance than either of the crankshaft moves.
in fact, we noted that elastic energy tends to its equilibrium point not only in a smoother and more natural way, but also more quickly with fewer mc steps.on the other hand, with regard to the average displacement of vertices in each trial, which provides the deformation measure of each tentative conformation, we noted that, as expected, the proposed move produces smaller deformations than either of the crankshaft moves.
however, the accumulated displacement of the proposed move is actually greater than those of both crankshaft moves in the first  <dig>  or so simulation steps ).
this high acceptance ratio in the initial simulation steps means that the proposed move generates a much more consistent deformation, the behavior of which obeys a logarithmic curve instead of the exponential curve that describes the accumulated deformation of each of the two crankshaft moves considered in this paper.
CONCLUSIONS
the crankshaft rotation method is the most common move found in the plasmid dna simulation methods for generating new dna conformations.
recall that this classic method first selects two random vertices of the dna skeleton, and then all the segments between these two vertices are rotated around the axis defined by them.
this move is not very effective because many trials are rejected by the mc method.
in addition to its low acceptance ratio, this method can generate unnatural movements with large portions of the dna molecule displaced at once, unless the relevant parameters are appropriately adjusted.
in this paper, we introduced a new move for plasmid dna through mc simulations.
in each iteration of the simulation, only one vertex and its two closest flanking vertices are subjected to the motion procedure.
thus, for each new trial, a single vertex is randomly chosen and then randomly displaced to a point within a small neighborhood.
to maintain connectivity of the dna chain, as well as the size of its segments, the two flanking vertices are also displaced but in a less free way.
thus, only three vertices are displaced in each new trial.
interestingly, considering that our algorithm generates small deformations in the transition from one dna conformation to another, we can conclude that it can be applied not only in the simulation of dna coiling, but also in real-time visualization.
we have already done this by employing the dna assembly algorithm that uses gaussian surfaces as geometric representations of nucleotides, as mentioned at the end of section ‘related work’.
in the future, we intend to incorporate a smoothing mechanism into our dna algorithm like, for example, that presented in  <cit> .
this will enable our algorithm to produce even more realistic simulations, eliminating the occurrence of slightly sharp corners like those shown in figure  <dig>  we also intend generating deformations that depend on the dna’s stiffness, which varies according to the sequence of nucleotides.
this will mean greater deformations on more flexible segments and smaller deformations on less flexible segments of dna.
finally, our ultimate goal is to be able to replicate in silico, and visualize what happens to plasmid dna during the production and purification processes in laboratory experiments.
competing interests
the authors declare that they have no competing interests.
authors’ contributions
the authors were jointly responsible for developing the algorithm and writing the manuscript.
both authors read and approved the final manuscript.

 plasmid dna is a family of dna molecules widely used in life sciences, more specifically in gene therapy research. these molecules are produced inside host cells in a supercoiled conformation, which is the desired conformation for therapy purposes. however, such molecules can lose their original conformation in the production and purification processes, assuming more relaxed or even linear conformations, owing to thermodynamic changes. one of the main challenges for researchers is to find optimal thermodynamic conditions for plasmid dna therapeutic application without losing its supercoiled conformation or, at least, minimizing the occurrence of relaxed or open dna molecules. for many years, computational methods based on laboratory experimental data have been proposed to model and simulate the dynamic behavior and conformational changes in pdna molecules under certain conditions. the monte carlo method has generally been accepted as a reliable tool for simulation purposes, and is seen as the standard. this iterative method tries to minimize the elastic energy of the molecule in each iteration step of the simulation process, testing the probability of acceptance of each new trial. the goal is to make the molecule converge to an equilibrium state after performing as few iterations as possible, i.e., maximizing the acceptance ratio of the trials without compromising the effectiveness and reliability of the simulation. to simplify the simulation process, each plasmid dna molecule is reduced to a linear skeleton with equal sized segments that represents the topological conformation of the molecule. random deformations are then applied to this skeleton, generating new trial conformations, which are either accepted or rejected. interestingly, the essence of the method used to randomly generate each new trial, referred to as the standard crankshaft move, has remained the same for many years, with its origins dating back to the early 1960s, more specifically in the context of lattice polymer chains. this move was later adapted for simulation of flexible molecules like dna using mc methods. however, the standard crankshaft move has a very low acceptance ratio of trials, i.e., many trials are rejected. moreover, it can present very unnatural behavior as it features very sudden motions along large portions of the molecule. to enhance the efficiency of mc moves, biasing was found to be a solution. however, as earl and deem noted, biasing a deformation move implies that the probability of moving from one state to another is no longer symmetric; consequently, the acceptance rule used must be altered to maintain the detailed balance. in this paper, we present a new unbiased move for plasmid dna, whose skeleton is a closed polyline. this move not only preserves the size of each segment and its connectivity, but is also very effective in maximizing the acceptance ratio of the trials and stabilizing the molecule, thereby allowing steady, gradual temperature changes during the simulation. our method also generates natural and realistic animations that can be used in real-time simulation and visualization. related work in this section, we briefly review the mc methods in computational biology and chemistry, as well as the generative methods for dna conformations that form the core of these mc methods. monte carlo simulations the mc simulation method is one of the most important methods used in dna simulations. this method, which was originally presented by metropolis et al., generates dna conformations combining energy calculations, random conformational changes, and statistics. frank-kamenetskii et al., vologodskii et al. and lebret were the first to use an mc method to present numerical results of the probability of the occurrence of knots on pdna. frank-kamenetskii and vologodskii also presented valuable information on dna torsional rigidity. a few years later, vologodskii et al. used mc simulations to study the conformational and thermodynamic properties of dna molecules with physiological levels of supercoiling. vologodskii also included a chapter on “monte carlo simulation of dna topological properties” in the book “topology in molecular biology”, and with rybenkov, they reviewed how conformational properties of dna catenanes can be studied using mc simulations. gebe et al. presented an mc algorithm to simulate supercoiling free energies in unknotted and trefoil knotted inextensible circular chains with finite twisting and bending rigidity, while marko et al. made use of mc simulations to study the relationship between the amount of twisting in dna molecules and its supercoiling. kundu et al. used an mc algorithm to explain denaturation characteristics in a supercoiled plasmid and calculate the probability of denaturation for each base pair at different supercoiling degrees. in their work on the relationship between knots and supercoiling, cozzareli et al. used an mc simulation procedure to generate an equilibrium set of conformations. mc simulations have also contributed to the understanding of the interplay between base-pair stacking interaction and permanent hydrogen-bond constraints in supercoiled dna elasticity. based on the fact that atomic force microscopy has generated images of supercoiled dna confined to a surface, which affects conformational properties such as twist and writhe, fujimoto and schurr modified an existing program, developed to perform mc simulations of supercoiled dna in solution, flattening the dna to simulate the effect of deposition on a surface. fujimoto and schurr also presented a method to estimate torsional rigidities of weakly strained dna. burnier et al. used mc calculations to identify a mechanism by which topoisomerases can keep the knotting level low. more recently, olson et al., in their paper “how stiff is dna?”, used mc simulations to understand the behavior of a long, double-helical polymer in the tight confines of a cell and in the design of novel nanomaterials and molecular devices. generative methods of dna conformations it has generally been accepted that supercoiled, i.e., the self twisting of the double stranded molecule over itself, is the desired conformation for pdna molecules. thus, it is necessary to measure the supercoiling of a given molecule. one of the most important quantitative measures of closed circular dna supercoiling is the linking number of the two dna strands, which is an integer corresponding to the number of double-helical turns of the molecule. there are several methods for calculating lk, with one of the most widely used involving the computation of two very important geometrical properties of closed circular dna molecules: twist and writhe. tw features the coiling of the two dna strands around the axis of the helix, while wr is a measure of the coiling of the helix axis in space. thus, the main result is: in our implementation, we used klenin and langowski’s computation method to calculate wr. owing to the nature of three-dimensional closed polylines, knots can occur in some pdna conformations. this is not a desirable feature, i.e., each closed circular dna molecule must remain unknotted during the simulation process, keeping its original topology even if supercoiling occurs. knot detection methods must be used during simulation to reject possible knotted conformations. we adopted harris-harvey’s knot detection algorithm, which uses the alexander polynomial to detect the existence of knots. this algorithm is based on the predicate that if two knots have different alexander polynomials, then the knots are topologically distinct. thus, because the alexander polynomial of an unknotted closed circular dna molecule is equal to one, all conformations for which this polynomial does not equal one must be rejected during the simulation. each trial conformation must be generated in such a way that the size and connectivity of each segment of the dna chain do not change. a major deformation method used to displace vertices of the dna chain was introduced by klenin et al.. this method, which is just a biased crankshaft move, starts by randomly choosing two vertices vm and vn. then, all the vertices are rotated a randomly selected angle θ around the axis defined by the line connecting vm and vn, as shown in figure furthermore, following klenin et al., the value of θ is uniformly distributed over a certain interval, and must be continuously adjusted during the simulation to guarantee that about half the steps are accepted.figure crankshaft motion. to increase the acceptance ratio of the simulation trials another type of motion has been proposed in the literature. this improvement performs a sub-chain translation, which is usually referred to as reptation motion, and is illustrated in figure first, two vertices vi and vj are randomly chosen. then, the sub-chain between vi and vj is translated by one segment length along the chain contour. the segment that was immediately after vj is also translated to fill the gap between vi and vi+ this motion suggests movement analogous to a snake slithering and, hence, the name reptation motion. other types of motion can be adopted if the metropolis microscopic reversibility requirement is satisfied, i.e., if the probability of each trial conformation is the same as that of the reverse movement.figure reptation motion. visualization of dna conformational changes over time is also important as part of the entire simulation of dna behavior. this is usually performed only when the entire simulation procedure ends and is typically done by assembling the dna atoms along the dna axis. interestingly, a more efficient dna assembly algorithm was presented to allow the visualization of all the steps of the simulation procedure in real-time. in this method, each dna nucleotide is represented by a three-dimensional building block, allowing the assembly of the entire molecule faster, but in a realistic way. in geometric terms, each of the four building blocks featuring dna nucleotides is a gaussian isosurface, which was previously generated by an algorithm that triangulates molecular surfaces. methods the deformation algorithm presented in this paper uses a linear skeleton with equal sized segments, henceforth called the dna skeleton. before introducing the core of the method itself, we explain how the dna skeleton can be created for use by the deformation algorithm. initial conformation of the dna skeleton the dna skeleton can assume any closed unknotted conformation. the simplest of these conformations is a completely relaxed circular conformation. besides, the length of each segment of the dna skeleton corresponds approximately to base pairs of the double helix. that said, the first step of the algorithm is to determine the number of segments of the dna skeleton ensuring around base pairs per segment. assuming that the dna has a sequence of n base pairs, we want to find an integer s denoting the number of segments of the dna skeleton. we define two integer parameters min and max, respectively, as the minimum and maximum numbers of base pairs that are admissible per segment, such that min < 30 <max. then, for each integer value i,min ≤ i ≤ max, we calculate the corresponding si=round. finally, we adopt s=si as the number of segments of the dna skeleton that minimizes |n-|. once we have the number of segments s, we just have to build a regular polygon with s sides inscribed in a circle. from the number of base pairs, we can infer the approximate perimeter of the circle, as well as the corresponding radius r, from which we obtain the first vertex of the skeleton at p0=. then, we apply s successive rotations to p about the origin to obtain all the vertices of the dna skeleton of the initial relaxed conformation; the rotation angle is given by α=2π/s. note that, although the initial conformation is circular, the methods for dna assembly and deformation apply to any initial conformation. skeleton deformation algorithm assuming we have a three-dimensional closed polyline pk representing the dna skeleton, we need to deform this polyline to obtain a new polyline pk+ with the same number of equal sized segments as pk, but without loss of its connectivity. let s be the number of segments of pk and {vi}, i= …,s- its vertices. we choose a random vertex vm, 0≤m≤s- as our mobile vertex, i.e., the vertex with the most motion freedom in the current trial conformation. any movement of the mobile vertexvm implies movement of its closest neighbors vm- and vm+ called semi-mobile vertices. the remaining neighbors vm- and vm+ are fixed vertices because they do not move in the deformation. thus, in each deformation step, only three vertices will be displaced: vm, vm- and vm+ figure mobile vertex v m can be displaced randomly in the intersection of three spheres, n m, s m- and s m+2. however, vm cannot be freely displaced. in the first instance, vm moves within the sphere nm centered at vm with radius r= Δ, where Δ= Å is the distance between two consecutive base pairs. more specifically, the new position of vm is found randomly within the region resulting from the intersection of the three spheres, nm, sm- and sm+ the latter two spheres with radius 2l are centered on the fixed vertices vm- and vm+ respectively, where l is the length of each segment of the dna skeleton. note that the optimal value Δ for r was found experimentally and based on the rate of successful trial conformations in the first attempt, about percent, though this rate remained high for a value of r up to Δ. the small radius r of sphere nm ensures that the transition from pk to pk+ occurs without noticeable jumps. as noted above, the new position of vm was found randomly within the region nm∩sm-2∩sm+, but no explanation of this random procedure was given. in fact, to calculate the new position of vm, we first convert its cartesian coordinates to spherical coordinates relative to vm- where d is the distance between vm- and vm. next, we randomly generate a new position for vm as, where Δd∈ and Δθ,Δϕ∈. it is clear that the displacement of the flanking vertices vm- and vm+ depends on the previous movement of vm. here we focus on the computation of the new position of vm- since the new position of vm+ can be calculated similarly. for this purpose, we also convert the cartesian coordinates of vm- to spherical coordinates relative to vm- where l is the radius of the three spheres sm, sm- and sm- centered on vm, vm- and vm- respectively. moving vm- to a new position must be done in such a way that its distance l to vm- and vm remains unchanged. in other words, the new vm- must lie on the circumference resulting from the intersection of the two surfaces bounding sm and sm-. if Δd= the new position of vm- relative to vm- is given by; otherwise, the new location of vm- is, where Δψ is the angle of the angular motion of vm- on sm- resulting from the translational displacement Δd of vm along the line defined by vm and vm-. we compute Δψ by rearranging the equation that describes the reciprocal motion of the piston with respect to the crank angle as follows: figure displacement of vertices v m and v m-1. translational piston move of vertex v m translates into a rotational move of v m-1. note that applying the translational displacement Δd to vm before the rotational motions Δθ and Δϕ means that Δθ= in eq.; otherwise, Δθ≠ in summary, moving vm implies a translational and two rotational motions relative to vm- expressed in spherical coordinates. this causes vm- to rotate accordingly on the sphere centered at vm- with part of this rotational motion determined by the translational displacement Δd of vm. these types of moves satisfy the principle of microscopic reversibility, although this is not critical for our purposes because the simulation procedure is only used to locate energy minima. as noted by mauri in, for a conservative n-body system, as in the case of a dna molecule, microreversibility stems from the invariance of the equations of motion with respect to time reversal, i.e., every microscopic motion reversing all particle velocities also results in a solution. this leads to the so-called principle of detailed balance, which states that under stationary conditions each possible transition from one conformation to another balances itself with the reversed transition in time. in other words, the probability of obtaining trial conformation pk+ if the current conformation is pk must be equal to the probability of obtaining trial conformation pk if the current conformation is pk+1. with this in mind, and having calculated the constrained position of vm- as a consequence of the move of vm, we need to determine its new position after rotating it randomly about the axis defined by vm- and vm. it is clear that the old and new locations of vm- lie on the circumference resulting from the intersection of spheres sm- and sm. likewise, we find the new position of vm+ after rotating it randomly about the axis defined by vm+ and vm. interestingly, these two rotations can be seen as two particular crankshaft rotations. finally, it is worth noting that the deformation algorithm described above can also be used in other biochemical systems such as internal coordinate models of cyclic peptides, as well as in some mechanical problems related to articulated arms and chain moves. in fact, this algorithm can be used to randomly deform any regular polygon in two or three dimensions with guaranteed preservation of connectivity and the length of segments. dna assembly algorithm for a realistic visualization of closed circular dna simulations in real-time, we combine the new deformation algorithm described above with the dna assembly algorithm introduced by raposo and gomes. this dna assembly algorithm uses four three-dimensional building blocks representing dna nucleotides, namely, adenine, cytosine, thymine, and guanine. each building block is a pre-triangulated isosurface generated by a triangulation algorithm for molecules. the assembly procedure for nucleotides can be thought of as the operation of wrapping helicoidal dna backbones around cylinders along the dna skeleton. the algorithm iterates over n base pairs of the dna, assembling a single base pair bibi per iteration. each iteration comprises four distinct stages: generation of geometric instances for nucleotidesbiandbi. considering that there are only four possible base pairs, given a nucleobase ni of a dna strand, two geometric instances of nucleotides must be generated, the first for the building block bi and the second for the corresponding building block bi. positioning of the base pairbibion planez= base pair bibi is positioned as if it was the first base-pair of the dna molecule, that is, it is placed perpendicular to the segment that starts at the origin in such a way that its center coincides with the origin. alignment of base pairbibiwith the plane perpendicular to segment i. note that this alignment involves a rotation about the origin of the coordinate system. translation of base pairbibito the plane perpendicular to segment i. finally, because all geometric transformations are performed around the origin, building blocks bi and bi must be displaced to their correct positions relative to the midpoint of the corresponding segment i of the dna axis. dna building blocks: adenine, cytosine, guanine, and thymine. it is important to note that this dna assembly algorithm does not take into account the sharp kinks that may occur at the junctions of the conformation segments, as shown in figure nevertheless, a possible solution to this problem is the smoothing procedure proposed by kummerle and pomplun. for detailed information about the dna assembly algorithm, the reader is referred to.figure detail of dna segments assembled around a random piece of the skeleton. monte carlo simulation mc simulations are iterative methods based on the concept of elastic energy of closed circular dna and on stochastic parameters and calculations aimed at converging to the energetic and thermodynamic equilibrium of the molecule. the main principle is to perform random dna deformations and check whether the resulting new conformations should be accepted or rejected according to energy changes and acceptance probability. more specifically, a random deformation of the dna is accepted if it reduces the elastic energy of the molecule or has some probability of occurring. in the experiments and results presented in this paper, we used the same mc simulation method and parameters as those used in, where elastic energye is calculated as here eb is the bending energy given by where kb is the boltzmann constant, t is the temperature, α= is the bending constant, and Θi is the angular displacement between the directions of segments i and i+ torsional energyet is given by where c denotes a constant parameter known as the torsional rigidity, l is the total length of the chain, and wr is the writhe of the skeleton. the linking number difference Δlk in is the difference between the linking number lk of the dna molecule and that, lk of its relaxed dna conformation where - ≤σ≤- is the superhelix density parameter. for calculating writhe wr, we used the method proposed in, more specifically, method 2b. this method is based on the principle that writhe can be calculated as the difference between linking number lk and twist tw: this method for computing wr uses an auxiliary chain close enough to the dna skeleton, and with as many segments si as the dna skeleton. then, considering that ri is the initial point of segment si, we can find the directional writhe as follows: where for the computation of, we must check whether segments si and sj cross, i.e., whether their projections π and, π onto plane z= intersect. we used lamothe’s algorithm to check whether the projections of these segments intersect. in turn, and following klenin and langowski, twist is given as where pi denotes the vector, normal to both si- and si, z denotes the z-th component of vector pi vector, and where u is the unit vector in the z axis direction. then, using the results of and, we get the final writhe number: once we know how to perform the necessary energy calculations, we can apply the mc method to each iteration of the simulation to obtain a new dna conformation from a random deformation of the dna skeleton. then, we calculate the energy ei+ of the new candidate conformation, and compare it with the energy, ei, of the previous conformation. the new conformation i+ is accepted if ei+1<ei or where tm is the temperature of the experiment and ρ is a random value between and. knots detection it is important to note that knots can occur when random deformations are applied to dna conformations. because this is not desirable, i.e., dna supercoiling must occur without generating knots, we must check for the existence of knots and reject the deformation if we find one or more knots. to optimize the performance of the method, this checking procedure is done before the mc acceptance test, avoiding unnecessary energy calculations. for knot detection we used the method of harris and harvey. in this method, based on the principle that two knots are topologically distinct if they have distinct alexander polynomials, the dna skeleton is converted to a knot data structure, and its alexander polynomial is computed and compared with the alexander polynomial of the circle, which is a trivial knot. if these two polynomials are different, the dna skeleton contains at least one non-trivial knot.  to evaluate the effectiveness and performance of our deformation method when applied in mc simulations, we performed a set of experiments comparing our method with two types of dna chain moves, namely, the standard crankshaft move and the biased crankshaft move. the standard crankshaft move is a randomly chosen move. in fact, the ends, vm and vn, of each sub-chain are randomly chosen, as is the case with the rotating angle θ of the sub-chain around the line that passes through its ends. that is, the standard crankshaft move does not adjust the size of the sub-chain nor the rotation angle in any way. on the other hand, as in the deformation method introduced by klenin et al., the biased crankshaft move used here adjusts only the rotation angle. recall that this type of biased moves is a way of enhancing the efficiency of mc moves, because it allows us to choose moves with a higher acceptance ratio. through these experiments we aim to demonstrate that our method generates a smoother and more controlled deformation, which leads to more consistent and even faster convergence to molecular energy equilibrium. experimental setup three experiments were performed to compare the proposed method with two classic methods, namely, the standard crankshaft move and biased crankshaft move. we used a setup based on kummerle and pomplun’s work for the puc plasmid dna molecule. all the three experiments were performed using the same closed circular dna sequence with base pairs and exactly the same conditions and mc simulation parameters for the three monte carlo moves under comparison, that is: kb= -23; α=; c=3×10-28; and σ=-. however, we performed experiment a at a constant temperature of k, while experiments b and c were performed progressively reducing the temperature from k to k. finally, it should be noted that the mc simulations were performed on an 64-bit windows laptop computer with an intel i ghz cpu, gb ram and an nvidia geforce gt mx gb graphics card. experiment a: puc with constant temperature in experiment a, we performed an mc simulation with steps using the puc closed circular dna sequence at a constant temperature of k as in. this experiment was replicated using: the standard crankshaft move, the biased crankshaft move as described in, and the proposed method. in the particular case of the biased crankshaft move, after an exhaustive optimization procedure with more than steps, we came to the conclusion that the crankshaft rotation angle should initially be in the range, as shown in figure furthermore, as expected, this angle decreases over time figure figure optimization of initial angle range for steps. experiment a: crankshaft rotation angle for steps at a constant temperature of k. we used two measures to compare the efficiency of the three methods: elastic energy equilibrium and acceptance ratio of trials. the graphs of elastic energy for the three methods are shown in figure where we can see that the average elastic energy for each method remains approximately the same over time. nevertheless, on average, the elastic energy of the proposed method is slightly higher than that of the standard crankshaft method, which in turn is higher than the energy associated with the biased crankshaft method.figure experiment a: elastic energy for steps at a constant temperature of k. on the other hand, the acceptance ratio of trials was evaluated for each slice of steps from a total of steps, according to the acceptance condition. the acceptance ratio was steadily higher for the proposed method, always remaining above accepted steps for each slice of steps, i.e., an average acceptance ratio around 45%, and achieving higher ratios around 50% in the second half of the experiment. on the contrary, the acceptance ratio for the standard crankshaft move was always under 30%, and even lower in the first steps of the experiment. with respect to the biased crankshaft move, the average acceptance ratio was slightly above 30%, but far below the results obtained using the proposed method. this indicates that our new method generates trials with much higher probabilities of being accepted under the mc simulation conditions at a constant temperature, i.e., it minimizes the number of trial rejections, and avoids useless computations.figure experiment a: acceptance trials for slices of steps from a total of steps at a constant temperature of k. during experiment a, we also noted that the crankshaft moves generated a few hundred conformations that were rejected owing to the existence of knots. on the contrary, our new method did not produce any knots at any time during the simulation steps because the deformation is done smoothly and without conformational jumps. this concurs with the fact that, despite dna molecules in living cells being long and compactly coiled, they rarely get knotted, which suggests that supercoiling inhibits dna knotting. experiment b: puc with variable temperature experiment b also involved mc steps. this experiment was also replicated for each of three methods analyzed in this paper. in the particular case of the biased crankshaft move, once again, after an exhaustive optimization procedure with over steps, we concluded that the initial crankshaft rotation angle should be in the range because, in experiment a, the temperature is not constant. adjustments to this rotation angle range during the experiment are shown in figure more specifically, the temperature decreases with energy, i.e., if the average elastic energy of a -step interval is higher than that of the previous slice of steps multiplied by a factor, the temperature is also multiplied by a factor. in fact, the temperature decreased progressively from k to k. as expected, the closer the method converges to the energy equilibrium, the greater is the decrease in temperature. as in experiment a, the acceptance ratio of the proposed deformation method was always higher than those of the classic deformation methods.figure experiment b: crankshaft rotation angle for steps with temperature varying between k and k. experiment b: acceptance trials for slices of steps from a total of steps with temperature varying between k and k. it was a somewhat surprising to observe how the energy decreased during the simulation. as shown in figure, when using the proposed method, the elastic energy of the molecule converged sooner and more consistently to equilibrium. as shown, we achieved energy equilibrium after approximately mc steps, while the crankshaft moves only stabilized after steps. besides, the standard crankshaft move generated a number of very slight energy jumps, i.e., the energy did not decay as consistently as in the proposed method. however, the energy level at equilibrium was the same for all three methods, approximately ×10- no less meaningful was the temperature decay during this experiment. as presented in figure, when using the proposed method, the mc temperature decreased more rapidly and in a more consistent way, i.e., the graph for the new method is much smoother with the advantage of reaching the equilibrium temperature sooner. figure shows the final result of experiment b.figure experiment b: elastic energy during a -step experiment with temperature varying between k and k; temperature decaying during a -step experiment. experiment b: puc after mc steps with temperature varying between k and k. in summary, we can say that the proposed deformation method requires fewer simulation steps to achieve energy equilibrium, largely owing to its high acceptance ratio. experiment c: average displacement in experiment c, we set out to measure the amount of deformation of plasmid dna caused by each type of mc move. this was accomplished by computing the average displacement of the dna skeleton vertices for each accepted trial. in this experiment, we only considered the standard crankshaft move and the move proposed in this paper. taking into account vertices vi,vi+ …,vi+n that are displaced during a simulation trial, we determined the distances di,…,di+n between the new positions and the previous positions of these vertices, and straightforwardly computed the average displacement given by /. finally, we considered the accumulated displacement for a slice of steps as the sum of the average displacements of the accepted trials of that slice.more specifically, we performed a -step simulation for each of the methods, namely, the standard crankshaft move and the proposed move. as shown in figure, the new move generates much smaller average displacements than the standard crankshaft move. besides, from figure we can see that the new move generates displacements right from the start of the simulation, whereas the standard crankshaft move starts to produce displacements later. this can be explained by the high acceptance ratio of the new method, as well as its more steady deformations.moreover, from figure we conclude that smaller displacements in each trial do not mean there will be smaller accumulated displacements. the accumulated displacements of the new method form a logarithmic curve, while the curve of the standard crankshaft move is clearly exponential ). this means that in the new method, the closer we get to the point of energy equilibrium, the shorter is the displacement toward a stable conformation, i.e., the average displacement for each accepted trial converges to zero as the conformation converges to the equilibrium. this is not the case in the standard crankshaft move, as illustrated by the accentuated variations in average displacement of the trials in figure.figure experiment c: average displacement per steps for puc steps; average accumulated displacement for puc steps. discussion as mentioned above, we used the same molecule, the same conditions/parameters, and the same mc simulation method in experiments a, b, and c. for comparison purposes, each experiment was performed using three different deformation methods: the standard crankshaft move, the biased crankshaft move, and the proposed method. as expected, the acceptance ratio of trials for the proposed method is higher than that for either of the crankshaft moves. the acceptance ratio of the new method is almost always greater than 40%, and even reaches more than 60% at certain times. moreover, unlike the crankshaft moves, the acceptance ratio for the proposed method is very high from the very first steps of the simulation. more specifically, in a scenario with decreasing temperature, we obtained an acceptance ratio of more than 60% for the proposed method compared with 5% for the standard crankshaft move and 20% for the biased crankshaft move in the first mc steps. with respect to elastic energy, the experiments also show that as the temperature decreases the new move achieves better performance than either of the crankshaft moves. in fact, we noted that elastic energy tends to its equilibrium point not only in a smoother and more natural way, but also more quickly with fewer mc steps.on the other hand, with regard to the average displacement of vertices in each trial, which provides the deformation measure of each tentative conformation, we noted that, as expected, the proposed move produces smaller deformations than either of the crankshaft moves. however, the accumulated displacement of the proposed move is actually greater than those of both crankshaft moves in the first or so simulation steps ). this high acceptance ratio in the initial simulation steps means that the proposed move generates a much more consistent deformation, the behavior of which obeys a logarithmic curve instead of the exponential curve that describes the accumulated deformation of each of the two crankshaft moves considered in this paper.  the crankshaft rotation method is the most common move found in the plasmid dna simulation methods for generating new dna conformations. recall that this classic method first selects two random vertices of the dna skeleton, and then all the segments between these two vertices are rotated around the axis defined by them. this move is not very effective because many trials are rejected by the mc method. in addition to its low acceptance ratio, this method can generate unnatural movements with large portions of the dna molecule displaced at once, unless the relevant parameters are appropriately adjusted. in this paper, we introduced a new move for plasmid dna through mc simulations. in each iteration of the simulation, only one vertex and its two closest flanking vertices are subjected to the motion procedure. thus, for each new trial, a single vertex is randomly chosen and then randomly displaced to a point within a small neighborhood. to maintain connectivity of the dna chain, as well as the size of its segments, the two flanking vertices are also displaced but in a less free way. thus, only three vertices are displaced in each new trial. interestingly, considering that our algorithm generates small deformations in the transition from one dna conformation to another, we can conclude that it can be applied not only in the simulation of dna coiling, but also in real-time visualization. we have already done this by employing the dna assembly algorithm that uses gaussian surfaces as geometric representations of nucleotides, as mentioned at the end of section ‘related work’. in the future, we intend to incorporate a smoothing mechanism into our dna algorithm like, for example, that presented in. this will enable our algorithm to produce even more realistic simulations, eliminating the occurrence of slightly sharp corners like those shown in figure we also intend generating deformations that depend on the dna’s stiffness, which varies according to the sequence of nucleotides. this will mean greater deformations on more flexible segments and smaller deformations on less flexible segments of dna. finally, our ultimate goal is to be able to replicate in silico, and visualize what happens to plasmid dna during the production and purification processes in laboratory experiments. competing interests the authors declare that they have no competing interests. authors’ contributions the authors were jointly responsible for developing the algorithm and writing the manuscript. both authors read and approved the final manuscript. 
plasmid dna molecules are closed circular molecules that are widely used in life sciences, particularly in gene therapy research. for that purpose, we have adopted a dna assembly algorithm that uses nucleotides as building blocks. monte carlo methods have been used for several years to simulate the conformational behavior of dna molecules. in this paper, we present a new algorithm for the deformation of plasmid dna molecules for monte carlo simulations. we also present the results of three experiments comparing our deformation move with the standard and biased crankshaft moves in terms of acceptance ratio of the trials, energy and temperature evolution, and average displacement of the molecule. in each iteration these simulation methods randomly generate a new trial conformation, which is either accepted or rejected according to a criterion based on energy calculations and stochastic rules. compared with both crankshaft moves, our move generates simulation trials with higher acceptance ratios and smoother deformations, making it suitable for real-time visualization of plasmid dna coiling. these simulation trials are generated using a method based on crankshaft motion that, apart from some slight improvements, has remained the same for many years. the move underlying our algorithm preserves the size and connectivity of straight-line segments of the plasmid dna skeleton. plasmid dnasimulationmonte carloissue-copyright-statement© the author 2014
 plasmid dna molecules are closed circular molecules that are widely used in life sciences, particularly in gene therapy research. monte carlo methods have been used for several years to simulate the conformational behavior of dna molecules. in each iteration these simulation methods randomly generate a new trial conformation, which is either accepted or rejected according to a criterion based on energy calculations and stochastic rules. these simulation trials are generated using a method based on crankshaft motion that, apart from some slight improvements, has remained the same for many years.  in this paper, we present a new algorithm for the deformation of plasmid dna molecules for monte carlo simulations. the move underlying our algorithm preserves the size and connectivity of straight-line segments of the plasmid dna skeleton. we also present the results of three experiments comparing our deformation move with the standard and biased crankshaft moves in terms of acceptance ratio of the trials, energy and temperature evolution, and average displacement of the molecule. our algorithm can also be used as a generic geometric algorithm for the deformation of regular polygons or polylines that preserves the connections and lengths of their segments.  compared with both crankshaft moves, our move generates simulation trials with higher acceptance ratios and smoother deformations, making it suitable for real-time visualization of plasmid dna coiling. for that purpose, we have adopted a dna assembly algorithm that uses nucleotides as building blocks. keywords plasmid dnasimulationmonte carloissue-copyright-statement© the author 2014 