WASP: a software package for correctly characterizing the topological development of ribbon structures

We introduce the Writhe Application Software Package (WASP) which can be used to characterisze the topology of ribbon structures, the underlying mathematical model of DNA, Biopolymers, superfluid vorticies, elastic ropes and magnetic flux ropes. This characterization is achieved by the general twist–writhe decomposition of both open and closed ribbons, in particular through a quantity termed the polar writhe. We demonstrate how this decomposition is far more natural and straightforward than artificial closure methods commonly utilized in DNA modelling. In particular, we demonstrate how the decomposition of the polar writhe into local and non-local components distinctly characterizes the local helical structure and knotting/linking of the ribbon. This decomposition provides additional information not given by alternative approaches. As example applications, the WASP routines are used to characterise the evolving topology (writhe) of DNA minicircle and open ended plectoneme formation magnetic/optical tweezer simulations, and it is shown that the decomponsition into local and non-local components is particularly important for the detection of plectonemes. Finally it is demonstrated that a number of well known alternative writhe expressions are actually simplifications of the polar writhe measure.

no evidence that the equality in (1) holds if the ribbon is not closed. Therefore, it cannot constrain the interplay between internal twisting (Tw) and global self entanglement of the ribbon's axis (Wr). Many interesting target applications which can be modelled by ribbons are open structures such as DNA molecules subjected to optical tweezer experiments and chromatin fibers.
To overcome this problem, Fuller originally suggested that one could extend the closed ribbon theorem to open ribbons by artificially extending the structure 15 as shown in Fig. 1b. This approach has been popular in the field of DNA modelling 7,8,[17][18][19][20][21] and elastic tube applications [22][23][24] . It is complicated by the fact that the closure will generally contribute to both the writhing and the linking of the composite and hence the extraction of clear geometrical insight from the calculated quantities cannot be consistently achieved 16,25,26 . Alternative approaches have been formulated which approximate the writhing value 7,8,15,17,23,27,28 , however, these will in general not charcterize the writhing 16,25 and hence cannot consistently be used to constrain the topology of the ribbon. Finally, some authors have simply chosen to utilize the definition of the original closed formula without a closure. In this case the Lk measure is no longer a topological invariant and the relationship loses its fundamental topological anchoring (the sum Tw + Wr is no longer fixed).
Previously, Berger and Prior 16 introduced a version of (1) which is applicable to open-ended ribbons. It is based on a definition of Lk (called the net-winding) which is invariant to all deformations of the ribbon that do not permit rotation at its ends. In supercoiling ribbon models it would be equivalent to the natural linking of the DNA plus the applied over-coiling of the structure. The twisting Tw is the same quantity as in the closed case and a definition of the writhing (termed the polar writhe) was based on the difference of these quantities. For closed ribbons this is equivalent to the original Călugăreanu theorem, but, since the net winding is also an appropriate topological constraint for open ribbons, it has a wider range of applications than the original theorem. This framework has been regularly applied in the field of Solar physics since the original paper e.g. [29][30][31][32][33][34][35] and is increasingly being applied in Elastic tube modelling 26,36,37 as well. However, it has not been widely adopted in the biophysical community and artificial closure methods are still readily employed for problems such as characterizing DNA supercoiling.
Here, we seek to introduce polar writhe to the biophysical community and to provide a user-friendly tool for its implementation. The first aim of this article is to introduce the open-ribbon version of Eq. (1) derived in 16 and extended to consider knotting deformations in 26 . Both papers are somewhat technical, necessarily involving mathematical proofs. In this note we illustrate the properties of the polar writhe which underpin this open decomposition though instructive examples of constrained DNA supercoiling simulations. The second aim is to demonstrate how decomposing the polar writhe into local and non-local components, which measure "spring-like extension" deformations and "buckling-type" deformations/supercoiling changes respectively (see Fig. 1c,d), can provide critical additional insight into the curve's geometry through a series of instructive examples. Third, we demonstrate how an extension of (1) as defined in 26 can detect transitions to knotting possible in open ribbon structures even when the ends of the ribbons are constrained from moving. Finally, we introduce the Writhe Analysis Software Package (WASP) which can calculate the components of Eq. (1) in both closed and open ended cases. WASP allows non-expert users to easily apply these calculations to standard curve files and is particularly geared towards usage with molecular dynamics trajectories, supporting common trajectory file types including those generated by popular MD software including AMBER and OxDNA as well as more general MD file types such as PDB files. (c,d) illustrate what is meant by local and non-local writhing. The curve in (c) coils helically at its centre; the local writhe measures this helical coiling along the curve's length. The curve in (d) is knotted/self entangled, that is to say distinct sections of the curve wrap around each other. This is non-local writhing.

Introduction to the polar writhe
In typical applications, Lk is a prescribed quantity (i.e. a fixed applied number of rotations of the structure in plectoneme experiments), and the twist Tw can be calculated as Lk-Wr. Thus, in this section, we focus on the definition of the writhing as the key quantity to be calculated. In 16,26 the definition of writhing was given the name the polar writhe which we label W p to distinguish from the pre-existing closed ribbon definition (the reason for this name will be clarified shortly). In what follows, we introduce the quantity W p through a series of instructive examples.
The polar writhe. The polar writhe 16 is the sum of a local component ( W pl ) and a non-local component The twisted paraboloids shown in Fig. 2a provide an intuitive example of the behaviors tracked by these distinct components. The parabolas are described by the following formula: The parameter θ , the winding angle, determines the number of rotations applied to the parabola (c.f. Fig. 2ai,ii). The height parameter h determines the relative stretching of this coil (c.f. Fig. 2aiii,iv). For this curve www.nature.com/scientificreports/ W pnl = −θ/2π . The factor of 2π indicates that W pnl measures the number of full rotations. The sign is due to the orientation of the curve as we shall explain shortly.
In the top row of Fig. 2a, h is fixed. In (a), θ = π , resulting in a single loop around the middle of the paraboloid. In this case W pnl = −1 . In b there is an additional looping and W pnl = −2 . The local component, which measures the local helical density of the paraboloid, is increased from (a)(i) to (a)(ii) due to additional winding at fixed height. In the bottom row we fix the winding ( W pnl = −4 ) but vary the curve's height h to highlight the effect on W pl . As the height of paraboloid is increased from (a)(iii) to (a)(iv), the overall winding (and hence the non-local writhe) is unchanged. However, as the coil becomes increasingly less tightly wound, the local writhing decreases as its helical density is decreased. Note that in these cases the two components are generally of opposite sign, however, this is not a general property of the pair (W pl , W pnl ) . As shown in 16 , if one chooses h ≈ 0.37 then W p ≈ 0 irrespective of the choice of θ (see Fig. 2aii). We see the added value of the local/non-local decomposition providing information about the increasingly tight helical coiling which is missing from the total sum. Curve splitting. The idea behind the polar writhe calculation is that we split the curve into sections sharing a mutual vertical height. If we have a Cartesian coordinate system (x, y, z) where z is the height, then we split the curve at turning points for which the curve changes direction from upward to downward pointing: ẑ · dx/dz = 0 . Examples of the parabolas (1 turning point, 2 sections) and a locally looped structure (2 turning points, 3 sections) are shown in Fig. 3a,b. In general, we label the subsections i of the curve that span heights z ∈ [z min i , z max i ] as x i . The local writhe quantifies the coiled geometry of each individual subsection and the non-local writhe quantifies the mutual winding of all distinct pairs of sections.
Local polar writhe. If the unit tangent vector of the curve x is T = dx/ds (s being the arclength of the curve), then: The numerator represents the rotation of the unit tangent vector around the vertical direction. It is positive if the curve coils in a right-handed fashion and negative if left-handed. The denominator indicates (i) the rotation is given more weight if the curve is rotating relatively tightly (as with a small h parabola) and (ii) the modulus sign means that the orientation would be the same for a curve whether it points up or down, i.e. it only depends on the local helical chirality of the curve. The fact that this rotation is measured around the ẑ direction, which on a unit sphere of directions is the north pole, and that W pl can be interpreted as an area on the unit sphere bound by the curve T and the pole is the reason for the name the polar writhe (see 16 ). However, this is not a crucial point in what follows and we do not mention it any further.
Non-local polar writhe. Consider the parabola example. The curve is separated into sections x 1 and x 2 (as shown in Fig. 3c). We define an angle � 12 (z) which is the angle made by the vector joining section x 1 to x 2 at a fixed height z as shown in Fig. 3c. W pnl measures the total number of rotations of this angle: www.nature.com/scientificreports/ To ensure that W p forms part of an invariant sum with the twisting, the orientation of the curve sections is accounted for by multiplying by an indicator function σ i with σ i = +1 if the curve section x i is moving upwards and σ i = −1 if it is moving downwards. Additionally, the calculation is counted twice. This is necessary to form part of the invariant sum W p + Tw . Thus for the parabola: For the looped curve shown in Fig. 3b we follow a similar procedure for the pairs (x 1 , The only extra complication is that the curve sections do not share the same mutual z ranges and the winding is only measured for subsections of the curve which share the same mutual z range z ∈ [z min ij , z max ij ] as indicated in Fig. 3d. The non-local writhe for the pair (x 1 , x 3 ) would be: as indicated in Fig. 3e. Note that the product of σ 1 σ 3 here is positive as both sections have the same vertical orientation. The total W pnl of this curve would be: In general, W pnl is just the mutual winding of all subsections of the curve which share a mutual height range. In the general case where the curve x i has n sections (and n − 1 turning points) the non-local writhing is calculated as: Looped local to non-local transition. We consider a second illustrative example: a curve transitioning from a locally helically coiled curve to a looped curve, the first part of plectoneme formation. In this case, the looped curve is the curve shown in Fig. 3b and used as an example above. This transition is described by a set of curves which are equilibria of elastic rod models: Here s ∈ [−5, 5] is the curve's arclength and τ is a parameter which acts to twist the tube as it is decreased, here from 3 to 0.2. As the tube is twisted, its axis transitions from a curve with a helical perversion around its centre, see Fig. 4(i), to a looped curve indicative of the beginning of plectoneme formation as shown in Fig. 4(ii). One can obtain this behaviour by steadily twisting a cable under tension; that they are are equilibria of an elastic tube model is shown in 38 , and thus they represent the deterministic counterpart of a worm-like ribbon model e.g. 7,8 . In Fig. 4 we see its W p values as τ is decreased from 3 to 0.2 in 200 evenly spaced steps. The polar writhe W p is seen to steadily increase as the curve becomes at first increasingly helically coiled at its centre and then increasingly looped. For curves 1-138 all of the contribution to W p is local due to the helical kinking of the curve around its centre ( W pnl = 0). After this, the loop begins to form and the curve develops non-local writhing ( W pnl > 0 ). What is interesting is the dramatic inter-conversion of local and non-local writhing which occurs whilst the sum changes smoothly through this transition, indicating that the heilical coiling becomes distorted to form the loop. Again we see the local/non-local decomposition providing significant extra information which characterizes the geometrical development of the curve.
Non-local writhing detecting plectoneme formation. An example of a highly plectonemed curve is shown Fig. 4. Marked on the curve are 8 clear crossings indicating looped geometry. Since we have seen in the parabola example ( Fig. 2) that loops typically have a W pnl value of ±1 , we should expect the non-local writhing of this curve to have a magnitude of around 8. In fact, for this curve, W p = 8.88 , W pnl = 8.55 , and W pl = 0.33 (all to 3.s.f). Thus, we can infer that for highly plectonemed curves, W p will be largely dominated by non-local writhing and will roughly count the number of loops of the plectoneme.
Finally, for the interested reader, an additional set of example calculations including local/non-local decompositions of example curves can be found in Chapter 4 of 39 and 16 and examples of elastic rod equilibria can be found in in chapter 6 of 26 . Shortly, we will turn our attention to calculating this quantity for the output of numerical supercoiling experiments. www.nature.com/scientificreports/ Self crossing detection. Both the linking Lk (the net winding) and the polar writhe W p change by a a value of ±2 if the ribbon intersects itself (the twist changes continually). This would, for example, be relevant in DNA models for topoisomerase action. In 26,36 , this fact was used to detect whether solution branches for an elastic tube model were physically valid. A series of evolving elastic equilibria were obtained using a continuation method and when the values of W p jumped by ±2 the equilibrium path could be automatically ruled out as not physically accessible (as elastic tubes cannot self intersect). Of course in DNA modelling we typically expect this to be impossible due to localized repulsion of the polynucleotide chains.

Methods
Numerically simulating DNA tweezer experiments. Plectonemic DNA structures were generated via coarse-grained molecular dynamics simulations of magnetic/optical tweezer experiments. 600 base-pair linear DNA helices were constructed with three different net-windings: Lk = 40 (≈ 24 • per turn) (underlinked) , Lk = 60 (≈ 36 • per turn) (torsionally relaxed), and Lk = 80 (≈ 48 • per turn) (over-linked). Their evolution was simulated using the oxDNA2 model 40 . To simulate tweezers experiments, each helix was first aligned such that its axis was along the ẑ direction. Anchor restraints were placed on the bottom five base-pairs of each helix to permanently fix their position and harmonic restraints were applied to the top five base-pairs of each helix to allow the respective base pairs to move only in the ẑ direction and also to prevent the DNA from shedding its winding by simply untwisting. In addition to these restraints, a range of extension forces (0pN, 2pN, 4pN, 6pN, 8pN) were applied to each helix in separate simulations by pulling the top five basepairs in the +ẑ direction. Each helix was simulated for a total production run of 1.2×10 8 timesteps using the oxDNA2 model with a salt concentration of 0.5 M. We only detail three specific cases below for brevity. These cases are chosen as they highlight the critical characteristics of the polar writhe measure as applied to these numerical supercoiling experiments.  27 , were explored as well. Simulations were performed in explicit solvent using AMBER with a modified DNA BSC1 force field to correct for the circularized DNA and Smith and Dang ion modifications for Na and Cl ions [42][43][44] . DNA minicircles were solvated in a TIP3P rectangular water box with a 15.0 Å buffer. Following solvation, the appropriate cosolute (spermine) was placed along with the DNA within the water box and the entire system was neutralized with Na+ ions. Initial minimization was performed on the cosolute and water molecules with a 10 kcal/mol restraint placed on the DNA residues. After the initial minimization, restraints were removed from the DNA and the entire system was allowed to minimize. Following the minimization phase, restraints were reapplied to the DNA residues and the system was heated from 100 to 300 K. Consequent to heating, restraints were progressively removed from the DNA residues over a period of ≈ 4 ns. MD was performed without restraints on the equilibrated system for 100ns for each trial in the NVT ensemble. The number of spermine molecules in each simulation varied between trials to explore different ranges of cosolute concentration.
Obtaining a DNA axis. To calculate W p , it is first necessary to generate an axis for the DNA molecule. In order to properly characterize the geometry of a given DNA structure, it is important to generate an axis curve in a manner such that its total curvature is minimized. Failure to do so may result in the misrepresentation of the DNA geometry due to excess local writhe introduced to the axis curve as a result of the helical nature of the DNA backbone strands. Simplistic axis curve formulations such as the axis curve obtained by taking the midpoint of  www.nature.com/scientificreports/ the C1' atoms located on opposite sides of each base-pair along the helix result in coiled axes that exhibit the same helical periodicity as their encompassing backbone strands. To avoid such misrepresentation, the WASP package uses an implementation of the WrLINE axis curve method developed by Sutthibutpong, Harris and Noy 45 which effectively smooths the contour of the axis curve and eliminates the unwanted "coiling" and helical periodicity exhibited by alternative methods that generate unwanted excess local writhe. Full details regarding the WrLINE formulation can be found in 45 , but key features of the formulation are summarized below maintaining similar notation as in 45 .
Unlike a basic midpoint method in which an axis is obtained by simply calculating the midpoint m i of the pair of C1' atoms ( r C1',A & r C1',B -see Fig. 5) on opposite sides of a corresponding basepair j: for each base-pair along the DNA helix, the midpoint r i of each "dinucleotide" (2 base-pair) step along the DNA helix consisting of the two pairs of C1' atoms on opposing sides of a set of two consecutive base-pairs j & j + 1 on the helix is calculated as: The remainder of the WrLINE method can then be roughly summarized by the following procedure: 1. Given a dinucleotide midpoint r i , calculate the base-pair step twist θ i from the dinucleotide step corresponding to r i 2. Determine the number of base-pairs (2m) required to complete a full helical turn around the base-pair step corresponding to r i (i.e. m base-pairs above r i and m base-pairs below r i ). 3. Calculate a weighting factor w related to the sum of the base-pair twists θ i−1 , θ i+1 , θ i−2 , θ i+2 ... throughout the helical turn around r i . 4. Use r i , m, and w to determine a point on the axis curve h i .
The above procedure is then repeated for each midpoint r i to obtain a full set of axis points h i that represent the DNA helical axis. Note that the summary of the WrLINE method provided above has been included solely for the purpose of familiarizing the reader with the axis curve method utilized in WASP as well as to highlight the intricacies of determining a proper helical axis. The reader is strongly encouraged to refer to 45 for a comprehensive treatment of the WrLINE method.
It is important to note that the WrLINE method is intended to be used primarily with closed DNA structures such as DNA minicircles. As such, the algorithm requires that there exist a full helical turn of DNA around each point r i as specified above. This requirement is problematic in the case of linear DNA helices where a full helical turn of DNA cannot exist around points r i near the ends of the helix. In order to circumvent this issue, two options are available: 1. Treat the DNA structure as though it were closed (joined at the ends) despite the fact that it is, in reality, not. This allows the WrLINE method to be utilized verbatim as outlined in 45 , but will result in a number of undesirable stray axis points generated in arbitrary positions external to the DNA structure. The number of stray axis points will be proportional to the extent of supercoiling of the helix (how many base-pairs are required to make a full helical turn) and can be safely, manually deleted via knowledge of the system/visual inspection. These stray axis points will always constitute the "ends" of the unaltered axis as they are a product of end effects produced by the WrLINE algorithm. This option is the most general, but should be utilized with extreme caution. 2. Start the WrLINE method on the first points r i such that one full helical turn of DNA can be formed around the respective points r i on each end of the DNA helix (satisfying the criteria for the calculating the m parameter 45 ) as determined prior to beginning the axis curve computation. www.nature.com/scientificreports/ At the time of writing, options 1 and 2 are implemented in WASP via the deleteatoms and autodelete arguments respectively. In either case, the total number of resulting axis points will be less than the total number of basepairs when evaluating open DNA structures using the WrLINE method properly.
Calculating W p . The WrLINE algorithm yields a discrete curve x j n j=1 . The algorithm for calculating W p for this curve is as follows: 1. Split the curve into n sections x i by locating its turning points. This is done from the set of tangent vectors T j = x j+1 − x j and use of tricubic interpolation. 2. For each section x i calculate W pl (x i ) using Eq. (4). The WASP package uses a modified version of Simpson's rule 46 . 3. Identify the range of mutual overlap of all pairs of sections (x i , x k ) . The WASP package uses a simple binary search algorithm. 4. Use Eq. (9) to calculate W pnl (a branch cut tracking method is used to track full windings). This is implemented in C++ for the WASP package, but the main routines can be accessed by a Python interface.
In addition, for comparison to existing calculations, we calculate the writhe as defined by: using the algorithm specified in 47 . This is the expression for writhing used in the closed ribbon Călugăreanu theorem (applied here to open curves). As discussed in the introduction, this measure has been used in open ribbon studies to characterise the writhing geometry of the ribbon's axis and we present the results here for comparison. We do not present any results using closures. We highlight that whilst (13) is a valid quantity for open curves, it does not generally form an invariant sum with the twisting Tw and we lose the topological control associated with this invariance. In these plectoneme experiments (for example), this means the sum W p + Tw is exactly equal to the experimentally applied Lk whilst the sum Wr + Tw is not. It was demonstrated in both 16 and 26 that there is always a closure for which the closed ribbon structure will have the exact same writhe value as W p (in fact it is the typical stadium closure). A numerical demonstration of this fact is also detailed in 26 . Further, it was shown that the closure can account for a significant proportion of the calculation, hence obscuring the interpretation. This last fact is an additional reason that the polar writhe local/non-local decomposition is a preferable measure for quantifying the ribbon's evolving geometry in addition to the fact it forms part of an invariant sum for open ribbons whereas Wr does not.

Results
Undertwisted helices ( Lk = 40). Undertwisted helices exhibited relatively modest writhing dynamics throughout their trajectories. In contrast with torsionally relaxed helices, helical strain due to under-twisting results in the formation of small kink structures in the curve as indicated in Fig. 6a. We detail the results of two particularly pertinent cases here: a weak 2pN stretching force and a stronger 8pN force. The 2pN results are shown in Fig. 6b. The values of both writhe measures Wr and W p settle at values of about −1 , a twentieth of the applied (under) rotation. The W p value is consistently lower than the Wr value, although not by a substantial amount. As indicated in the figures, there is no plectoneme formation, however, partially looped sub-sections often form. The local/non-local decomposition of W p represents this mixture of helical (local) distortion and non-local coiling (see Fig. 6c). There are occasional spikes in the non-local writhing (relative to the local value). One such example is shown to arise from a tight loop formation as shown in Fig. 6d.
The 8pN results are shown in Fig. 7a. The additional force further restricts the writhing of the curve and the magnitude of both quantities is typically bound between [−0.5, 0.5] . The value of Wr is typically significantly larger than its W p counterpart. We verified that the ratio |W p − W r |/(|W p | + |W r |) was larger than 0.5 for over a third of the calculations ( > 50% difference). In Fig. 7b we see that for vast the majority of the curve's geometries there is only local writhing, a clear difference from the 2pN case which indicates that the force is acting to restrict loop formation. The non-local writhe values arose when small kinks in the structure developed, although again these are significantly restricted by comparison to the 2pN case. Intuitively speaking, in the context of this example, the local writhe is seen to reflect the vertical spring-like "bobbing" of the DNA in which the helix oscillates between stretching and compressing without any buckling/plectoneme formation occurring. The intermittency of the non-local writhe indicates the formation of buckling points that are not sustained, unlike the plectonemic formations in the following case.
Overtwisted helices ( Lk = 80). Overtwisted helices exhibit dramatic writhing as a result of high torsional strain on the DNA helix. As a result, helices contract over time in the ẑ direction by forming plectonemic structures (Fig. 8a). Time series plots of W p and W r are shown in Fig. 8b for a helix experiencing no extension force. There is a steady increase in writhing up to nearly 10 in both cases as the curve forms plectoneme structures. This is half the applied over-rotation which indicates that there has been a significant conversion of (over) twisting into writhing. We note that the two measures W p and Wr are nearly identical unlike in the under-twisted case. The reason for this is discussed in greater detail in the supplement. It suffices to state here that both assign the same value to the loops in the plectoneme structures which dominate the axis curve's geometry.
The local/non-local decomposition is shown in Fig. 8c. Generally, except in the initial stage of simulation, there is far more non-local than local writhing, again due to the plectoneme formation. The consistent, gradual www.nature.com/scientificreports/ increase of the non-local writhe reflects the formation of plectonemes within the DNA helix. In contrast, the local writhe remains relatively constant, fluctuating only slightly due to the significant stress on the DNA helix causing any spring-like compression of the helix to be quickly converted into plectonemic buckling. In the first curve shown in Fig. 8b, sampled frame 50 has significant helical looping in its upper part. Near the lower end of the curve we see the initial formation of the first plectoneme, this contributes non-local writhing. As indicated in Fig. 8c, there is a reasonable balance of local and non-local writhing reflecting this two part geometry. The looped section has formed a plectoneme by timestep 200 as shown in Fig. 8b. By this stage, the non-local writhing is dominant (Fig. 8c). The final curve shown in Fig. 8b shows two plectonemes, the second developed from one of the small loops in the upper half of the curve at step 200.

DNA minicircles.
To further highlight the utility of W p , we have included sample analyses of two DNA minicircles (Fig. 9). Both trajectories represent 108 base pair Lk = 14 minicircles exposed to a 10 mM spermine cosolute. The minicircle in Fig. 9a has the base pair sequence (GG) 27   www.nature.com/scientificreports/ loops/buckling points within the minicircle structures caused by spermine bridging. Spermine molecules adsorb to the backbone of the DNA minicircles early in the trajectories and gradually bridge together adjacent tracts of the DNA minicircles over time resulting in the formation of loops/buckling points within the minicircle structures as shown in Fig. 9c. The observed bridging dynamics are in accordance with results from previous studies as sequence dependent interactions have been observed between spermine molecules and AT/AA rich groups of linear DNA helices in which spermine molecules were observed to adsorb parallel to the DNA backbone and bridge together adjacent DNA helices 48 . As shown in Fig. 9a, the axis of the GC/GG minicircle forms a figure 8 structure with a smaller loop occasionally forming at a consistent location on the curve. The W p value is generally dominated by the non-local component with variations being linked to the development (and reduction) of the smaller localised loop. In panel (b) it is shown that the AT/AA minicircle quickly forms two clear loops giving a W p value which oscillates around 2. Again there is a small localised loop which variably develops/recedes. At the point t = 49 (marked with a dashed line) we see a relatively large spike in the local writhing corresponding to the loop localising and forming a tight local helical loop.
For both the AT/AA and GC/GG minicircles, the small localised loops that develop/recede within the minicircle structures are caused by spermine molecules that "jump" around the backbone of the DNA in various locations. These spermine molecules intermittently adsorb to the DNA backbone, bridging sections of the minicircle together which ultimately results in the formation of small loops that recede once the spermine molecules detach from the backbone and move to a new location. It is interesting to note that the non-local component of W p serves as a direct indicator of this phenomena, fluctuating directly in response to the aforementioned spermine/ loop formation dynamics.

Writhe star and winding classes
A complication with open ribbons such as extended DNA structures is the possibility of (un)knotting and belttrick type deformations such as those shown in Fig. 10 17 . These are classes of shape changes in which a section of the ribbon structure loops over one of its end points (as in Fig. 10b-d). In such cases, both Lk and W p can change by values of ±2 26 . Thus, when one prescribes a fixed value of Lk, it is only constrained up to an integer, but critically an integer change can only occur when an "over-the-top" deformation occurs. For many applications this will not be possible, for example, DNA structures in restricted environments and magnetic bead experiments where the bead itself is barrier. However, other potential applications of the writhing such as a measure of the development of tertiary structure in protein folding may need to consider such deformations. Prior and Neukirch 26  End angles and pulled-tight topology. We can think of the two ends of an open ended ribbon as being bound between planes. If any section of curve passes through the end of these planes, it makes an angle with the end point as shown in Fig. 11a. This angle forms part of the non-local writhing calculation (Eq. (9)). As the curve loops over the top of the end point, this angle increases (see Fig. 11a-d). This change would affect both W p and Lk. If the curve then comes back down on the other side of the end point then this angle will have gone through a change of 2π and since this is counted twice it leads to a change in ±2 of both W p and L k . Thus, overall the change of ±2 is detected. However, this means that both W p and Lk would be changing continuously (although still such that Lk − W p = Tw ). This change does not happen if the curve is bound between the two end points where Lk is a fixed quantity. For magnetic field applications where the end planes represent the edge of the system, i.e. the sun's surface 49,50 or the edge of the measurement domain in a plasma experiment 51 , evaluating this angle is crucial. However, for plectoneme formation and protein topology type applications it is important that the section of curve leaving the bounding planes is included in the calculation. Prior and Neukirch 26 derived a variant of the Lk = W p + Tw formula with quantities W * p and Lk * which in effect ignore these end angles and just register a jump in both W p and Lk. More precisely, it was shown one could always extend the ribbon with straight sections (with no Tw) such that the ribbon remains bound between two planes. This adds no local writhing or winding and in effect remains a measure of the topology of the original unextended ribbon, except in one crucial aspect: when the curve/ribbon passes over its end point. In this case, the ribbon self intersects and there is a jump in ±2 in both Lk and W p and thus Lk is fixed up to integer multiples of 2. This has the following intuitive interpretation. If we imagine taking ribbon's axis curve end points and pulling them directly apart, then a section of curve passing over the end point marks the point in its deformation where the curve would change its "pulled tight state", i.e. Fig. 10a-f. For those readers familiar with the literature of artificial closures, there would be the same jump when the curve intersects its closure.
In summary, the sum is invariant under all deformations which forbid curves passing over their end points and changes by values of ±2 when these over-the top deformations occur. This integer change can be used to indicate the ribbon's pulledtight topology (the shape obtained by pulling its ends directly apart) has changed.   Fig. 10 are part of a continuous set of deformations. This deformation was discretized into 360 steps. The first 100 steps, Fig. 10a,b, are simple in that the shape of the knotted (trefoil) curve does not change but its size with respect to the straight end extensions is increased; in short, the shape of the curve is effectively unchanged. From curve 101 to 260, the knot is undone as a section of the curve loops over its top end point, Fig. 10b-d. Then, from step 261-360, the curve is pulled straight (Fig. 10d-f). The net effect is the undoing of a knot. The values of W p , W * p and Wr (given by (13)) are shown in Fig. 12. For the first 100 steps, both W p and W * p are identical and unchanged, correctly representing the fact the knot's shape is essentially unchanged and the corresponding net winding would not change (assuming as always that the ends of the ribbon are prevented from untwisting). The Wr calculated via the Gaussian integral (Eq. (13)), however, does change. The reasons for this change are somewhat technical and are discussed in the supplement; it suffices to note here that it is not characterizing the knots unchanging shape and the sum W r + Tw would not be a fixed quantity. Figure 11. Illustrations of the loss in writhing which occurs to a section of the curve's interior looping over one of its end points and the extension used to make this loss a discontinuous jump which tracks the pulled tight topology of the ribbon. (a) An interior section of the curve rises above one of the planes containing the curve's endpoints. The angle represents the difference of two contributions to the W pnl calculation which involves angles made with the curve's end point. In (b,c) the curve section rises higher and the angle increases. In (d) the angle is nearly 2π and the curve section has (just) passed directly over the top of the curve's end. (e) is a ribbon corresponding to one of the curves in Fig. 11; a section of the curve is in the plane above the ribbon's end points. (f) The planar extension to the ribbon is shown. Now the whole ribbon structure is bounded between two planes. The extension is composed of planar curves with no twisting. www.nature.com/scientificreports/ The set of writhe values for configurations 101-260 (the over-the-top unlooping of the knot) highlight the difference between W p and W * p . For the first part of this deformation they are identical until the looped section passes above the end point's plane (as shown in Fig. 11). Then W p begins to measure a changing value due to the end angle measures between this loop and the end point. W * p does not account for this change. Near the end of this un-looping, the curve passes directly over the top of its end point and we see a jump of −2 in W * p where the measure has classed the curve as jumping into a different pulled tight topological class. In the final set of deformations 260-360, the curve is pulled straight; here all three measurements are in rough agreement that there is a steady decrease in writhe. Note, however, that W * p measure registers a relatively smooth change.

Examples.
Overtwisted DNA reconsidered. The discontinuous jump in the W * p measure can lead to jumpy time series if a section of the curve is continually deforming in the region directly above one of the curve's end points. It happens that this is exactly what is occurring in the overtwisted plectoneme formation simulation in the previous section. The comparison of W p and W * p is shown in Fig. 12b. The W * p measure has some significant discontinuities; each one was analyzed and in each case it was found to result from a section of the curve passing over the bottom end of the curve. An example is shown in Fig. 12c,d. A second observation is that W * p is consistently a small amount larger than W p . This is due to a difference between W pnl and W * pnl ( W pl is always the same for both), specifically the ignoring of the end angle depicted in Fig. 11. It is interesting that the molecule appears to deform such that it partly performs the over-the-top loop, a potential means to shed topology, but never completes this deformation, instead forming a second plectoneme loop.
A final mathematical point is that if one only calculated W * p for the given curve set, then, without looking at the individual curves themselves, it is hard to tell if the jump occurs due to self crossing or over the end deformations. The W p calculations, which have no changes on the order of 2 indicate it is not self-crossing but over the top deformation instead. Thus, one can make this distinction solely from the numerical calculations (rather than visulising the curve set). Of course we should have expected this for the DNA model used as repulsive forces prevent self-crossings, but in other scenarios the combination of information provided by both the W p and W * p quantities could be of significant utility.
The previous two examples highlight aspects of the W * p measure which the user should consider when deciding if it is appropriate for their application. The W * p measure was explicitly designed for elastic rod models in which both knotting and belt-trick style deformations clearly occurred: the end state had unambiguously changed. The pulled tight definition is designed to capture this change. The plectoneme simulation cautions that this change may not always be meaningful. In addition, it is also possible that the net change in W * p in some knotting deformation might add up to zero (an equal number of positive and negative jumps) even though the essential knotting entanglement might have changed. The WASP package gives the user the option to calculate both W p and W * p and we believe these examples will provide a guide as to how to interpret these measures.

Conclusions
Writhe is a fundamental measurement for characterizing the topology of ribbon structures and is extensively utilized in multiple fields ranging from biophysics to solar physics to aneurysm detection (to name just a few) [30][31][32]52,53 . However, the most commonly utilized formalism for writhe (as derived by Călugăreanu ) is only applicable to the subset of problems in which the system under study is a closed loop. Although this limitation has been overcome through the development of polar writhe, it is largely underutilized in the biophysical community.
Here, we have sought to rectify that problem by demonstrating the utility of polar writhe in DNA plectoneme and DNA minicircle calculations. Our results show that polar writhe is not only applicable to the analysis of open curves such as those formed by extended DNA structures but that it can be decomposed into local and nonlocal contributions that can provide additional information about DNA topology that cannot be obtained from the Călugăreanu formalism. In particular, we demonstrate how the non-local component of the polar writhe can be used as a potential indicator of plectoneme formation in both linear and circular DNA structures which we believe will be of particular value to the host of studies performed on effects of plectoneme formation in DNA and braided polymers 54, 55 . To aid in the adoption of polar writhe, we have developed a software package, WASP, that can be used to analyze molecular dynamics simulations of DNA. WASP is an open source software tool that is built to analyze trajectories from popular simulation packages. It is our hope that WASP will be utilized to rigorously analyze the growing field of DNA simulations in which changes of writhe are linked to biophysical phenomena (Fig. S1).