Let denote the number of states available to a linear molecule comprised of base pairs experiencing the B-Z transition. In any given state each base pair in the sequence is either a monomer (i.e. in B-form) or part of a dimeric pair with one of its neighbors (i.e. in Z-form). There are two possible arrangements for the first base pair in the sequence. It can be a monomer (i.e. in B-form), in which case there are ways that the rest of the sequence can be arranged. Otherwise, it can be in a dimeric unit with the second base pair (i.e. in Z-form). In this case the disposition of the first two base pairs is determined, so there are ways to arrange the rest of the sequence. Because these two alternatives are mutually exclusive and exhaust the possibilities, it follows thatUsing the fact that and , one can calculate for any length from Eq. (1). This recursion relation together with these initial conditions show that , the st Fibonacci number. An excellent approximation to this number for all is given byThis shows that the number of possible state in the linear B-Z transition grows exponentially with sequence length , but with base equal to the golden ratio , rather than base 2 as holds for transitions in which both states have monomer units.

Next we determine the free energy of each state of this molecule. This is comprised of the energy cost of the transition, plus the energy of accommodating the residual superhelicity. The quadratic free energy associated to residual superhelicity [52] is given byHere [53], [54], is the gas constant, and is the absolute temperature. The energy cost of the transition includes two factors. First, a nucleation energy for each run of Z-DNA is required to form the junctions between the B-form and the Z-form. This has been measured to be approximately 5.0 kcal/mol/junction, so the nucleation energy for each Z-DNA run is kcal/mol [21], [36], [40]. Second, the energy to transform the specified base pairs into the Z-DNA conformation is also needed. The published values of the B-Z transition free energies for all 10 dinucleotide pairs are given in the second and third columns of Table 1 [35]. The bases in each dinucleotide pair must alternate, one in anti and the other in syn conformation. As the values in the table show, the transition energy of a particular dinucleotide depends strongly on whether it is AS (5′ anti 3′ syn) or SA (5′ syn 3′ anti). A Z-Z junction occurs when adjacent dinucleotides have different anti-syn alternations, either (AS)(SA) or (SA)(AS). This violation is energetically costly, as shown in the last column of the Table.

Some of the energy values shown in Table 1 are calculated estimates, while others were measured experimentally [21], [35], [37]–[39]. The energies shown can be evenly divided between the two base pairs involved, whether in a dinucleotide or in a Z-Z junction, to determine the transition free energy associated with each base pair.

It is energetically most favorable for purines and pyrimidines in the Z-form to be in the syn and the anti states, respectively. There is a substantial energy cost for a base pair to have the opposite conformation. Although in principle every dimeric sequence can be driven into Z-form, four (CG, GC, CA=TG, and AC=GT) have substantially lower transition energies than do the others. The equilibrium distribution will be dominated either by the untransformed state or by states in which transition occurs at sequences composed of the energetically most favored dinucleotides.

Since the population of a state at equilibrium decreases exponentially as its free energy increases, the states that are neglected because their energies exceed the threshold are occupied with very low frequencies. If, for example, the threshold is set at kcal/mol, as is used below, the neglected states are occupied no more than times the frequency of the lowest energy state at K. A careful density of states calculation was performed to assess the aggregate influence of the neglected high energy states, and to approximately modify the calculated ensemble averages to account for them [49], [51]. Although this correction could not be performed on the transition probabilities of individual base pairs, this approach showed that all other ensemble averages were accurate to between four and five significant figures. Subsequently, we developed an exact, but very slow algorithm that performs these calculations [50]. By comparing its results with those of SIDD at different energy thresholds it was determined that, when thresholds in the range 10 to 12 kcal/mol were used, the accuracy of all calculated parameters, including the transition probabilities of individual base pairs, exceeds four significant digits [45].

Since there is a substantial nucleation energy kcal/mol associated with opening each run [21], [36], [40], states with numerous Z-runs will be correspondingly less populated. This allows us to impose a cutoff on the number of runs that may occur. In the SIDD analysis of strand separation it was found that three run states were the most that were encountered at any reasonable energy threshold and superhelix density . However, sequences that are energetically most susceptible to the B-Z transition tend to be shorter than the A+T-rich regions that favor denaturation. As shown by the energies in Table 1, there are many very costly types of “imperfections” which may hinder the extension of a Z-DNA run. So in the B-Z transition it can be energetically less expensive to initiate a new run at another favorable site than to extend an existing run into an energetically unfavorable region. Extensive sample calculations performed in tuning the SIBZ algorithm have shown that limiting consideration to states with a maximum of four runs is sufficient for high accuracy at the superhelical densities and sequence lengths of interest in this paper.

One important step in the SIBZ algorithm is the assignment of transition energies to the dinucleotide repeat units in each Z-run. This is complicated because the transition energy associated to each unit can have any of three significantly different values depending on the anti and syn characters of the base pairs in that unit and in its neighbors. Briefly, within each unit one base pair must be in syn and the other in anti. Also there is a significant energy penalty assessed in cases where there are Z-Z junctions with one or both neighbors.

In the SIBZ algorithm we use the following procedure to assign the energetically most favorable anti or syn conformations to all base pairs in a potentially Z-forming run. First, since the unit cell of Z-DNA is a dinucleotide, we allow only an even number of base pairs in any Z-run. We assign all purines in the run to be syn and all pyrimidines to be anti, as the dinucleotides with the lowest transition energies have this character. (See Table 1.) Since in most cases forming a Z-Z junction is more costly than flipping a base pair into its non-favorable conformation, we next change the conformation of a base pair if it has the same anti or syn character as both of its nearest neighbors. This procedure eliminates any repetition of the same conformation (anti or syn) longer than two base pairs. Then we search for quartets of the form AASS or SSAA. In these we flip the central two bases so that an alternating syn-anti character is obtained. Finally, when two bases within a dinucleotide unit have the same conformation, which is not permitted in Z-DNA, we flip the base pair which yields the minimum Z-Z junction energy with its neighbor. There is an ambiguity in this procedure for internal runs of even length at least four, such as ASSSSA, due to the order in which they are flipped. However, the states involved are always high energy because they will have at least one Z-Z junction as well as at least two unfavorable dinucleotides.

To determine the B-Z transition energies of a region, we scan its sequence by dinucleotide units, adding energies of syn-anti or anti-syn pairs according to Table 1. Whenever the alternating anti-syn character is disturbed we add the appropriate Z-Z junction energy. Also, we set the minimum length of a Z-run to be eight base pairs, as shorter regions have not been seen experimentally to form Z-DNA [21]. In this way we assign a B-Z transition energy to each segment in the DNA sequence of even length between 8 and 250 base pairs, which is a reasonable maximum cutoff for a Z-run length.

The algorithmic strategies for finding the lowest energy state, identifying all states that satisfy the threshold condition, and analyzing them to determine equilibrium values of parameters of interest are the same as those developed previously in the SIDD algorithm [45], [49]. We treat circular sequences as described there. A linear sequence is circularized by joining its end to its start with 50 T bases, and the resulting sequence is treated using the same method. In this case we do not report the information for the augmenting segment. These issues and methods have been described previously [45], [49], [51].

To assess its performance characteristics, we used SIBZ to analyze the pBR322 plasmid ( bp) using an energy threshold of kcal/mol. At superhelical density this analysis took 0.12 minutes to run on a MacBook Pro with dual Intel processors. On average there were 23.9 Z-form base pairs and 2.3 runs of transition. A total of 15,829,349 states were found to satisfy the energy cutoff condition . At superhelical density this analysis took 2.25 minutes to run on the same machine. In this case there were averages of 44.2 Z-form base pairs and 3.9 runs of transition, and 1,047,067,293 states satisfied the energy cutoff condition. We find that the execution time is almost constant at superhelix densities , and increases quadratically thereafter. The algorithm scales approximately linearly with sequence length. These performance characteristics suggest that the SIBZ analysis of the complete human genome at would take approximately 12 hours on a 100 CPU cluster of slightly faster (viz. Opteron) processors, if the sequence was partitioned into 5 kb segments that were analyzed individually. A similar analysis at would take approximately ten days.

This competition between different Z-forming regions within a sequence can lead to a rich repertoire of complex, interactive behaviors. We illustrate this with sample calculations on a designed sequence containing two regions susceptible to Z-DNA formation. This sequence consists of 5000 T base pairs, into which we insert two Z-susceptible regions at distant locations. The thymidine background is chosen to insure that only our inserted segments are susceptible to Z-formation. The first insert is a segment while the second is a segment that contains six Z-Z junctions. The segment is less costly to transform, because the Z-Z junctions in the segment are energetically expensive. However because it is shorter, transition at the segment also relieves less superhelicity. We used SIBZ to calculate the probability of transition of each of these regions over a range of negative superhelix densities. Plots of these probabilities as a function of are shown in Fig. 1.

Just beyond the onset of transition where is still small, the superhelical free energy of the untransformed state also is relatively small. Under these circumstance the energy relief afforded by transition is less than it is at more extreme superhelicities. So in this regime the magnitude of the transition energy is the dominant factor in determining which regions transform. This is why the shorter but energetically less costly Insert 1 is the first to transform, as shown in the figure. As increases the superhelical free energy becomes quadratically larger. Now transitions at longer sequences become more desirable because they relieve more superhelical stress energy. Under these circumstances the difference in transition free energy due to the ZZ-junctions becomes less important than the benefit afforded by transforming a substantially longer segment. For this reason a coupled transition-reversion event occurs around , in which transition of Insert 2 is coupled to the reversion of Insert 1 back to B-form. In the range it is energetically too costly for both segments to transform to Z-DNA simultaneously, so such states occur infrequently at equilibrium. Transition of the long Insert 2 has caused substantial relaxation,which decreases the residual superhelicity felt by Insert 1 below the value that would drive it to transform. So at these stress levels the probability of transformation of Insert 1 drops to near zero. As increases beyond the point where Insert 2 has a high probability of being entirely in Z-form, the additional stress accumulates as negative residual superhelicity. When this reaches a sufficient level Insert 1 again transforms to Z-DNA. Beyond both inserts have high probabilities of simultaneously being in Z-form. One sees that there are specific intervals of superhelicity within which 1) neither insert transforms, 2) the first transforms but not the second, 3) the second transforms but not the first, or 4) both inserts transform simultaneously. The transition in this example experiences every logical possibility.

We compared the performance of SIBZ with those of Z-Hunt and Z-Catcher when run on these sequences [35], [47], [48]. The energies from Table 1 were used in all three programs. When Z-Hunt was applied to the pBR322 sequence it found one 15 bp long segment at location 1448 with a Z-score of 2444, and a 17 bp long segment at position 1407 with a Z-score of 1845. As shown in Fig. 2a, SIBZ also finds these two peaks to be the most dominant, with the segment at location 1448 having the higher transition probability. This agrees with the relative Z-score rankings provided by Z-Hunt. However, the relative probabilities of these two regions are not proportional to their Z-score. Moreover, SIBZ documents several other regions that also have significant transition probabilities that Z-Hunt does not identify. One sees that, although the sites found by Z-Hunt agree with the major sites found by SIBZ, the latter provides more information regarding other Z-susceptible regions. The results from SIBZ, because they are expressed as transition probabilities of the fully competitive transition at the assumed superhelix density, are both more precise and more easily interpretable than are Z-scores.

Z-Catcher does not identify any Z-susceptible regions in pBR322 at superhelix density . At superhelix density it finds the two dominant segments at locations 1407 and 1448. It also finds two other Z-segments at positions where no Z-forming potential is seen by the other algorithms. The segment at position 1448 is found to comprise 28 base pairs, which is significantly longer than predicted either by Z-Hunt or by SIBZ. Thus, the predictions of Z-Catcher seem to differ considerably from those of either Z-Hunt or SIBZ.

The behavior of a B-Z transition varies significantly as the negative superhelical density is modified. In general, as is increased, larger numbers of transformed base pairs are required to relieve torsional stresses. If the most energetically susceptible region is sufficiently long, the most favored way to do this would be by extending the transition to encompass increasing amounts of this region. However, the most Z-susceptible regions in natural DNA sequences are usually relatively short, as is seen in these three sequences. In that case what commonly happens is that, as (and hence also the level of imposed stress) increases, successively more energetically costly distinct Z-forming regions are transformed (possibly coupled to reversion of other regions as shown above). In either of these strategies the average number of base pairs adopting the Z-form increases with negative superhelical density. This is demonstrated in Fig. 3, which shows the average number of Z-DNA base pairs as a function of for pBR322 and for X174. One sees that in both sequences there is a threshold for the onset of transition around , and the expected number of transformed base pairs increases approximately linearly thereafter. (We note that the onset of transition in Fig. 1 occurs at a slightly less extreme superhelix density than is seen in Fig. 3. This occurs because according to Table 1 the shorter Z-susceptible insert in the constructed sequence has the lowest possible transition energy. This is not true of the most Z-susceptible site in either X174 or pBR322.)

The Z-forming insert sequences 1 and 2 both have the same length, each containing 16 dinucleotides (i.e. 32 base pairs). This means they have the same ability to relieve superhelical strain. However, Sequence 2 contains a Z-Z junction consisting of GG bases, which breaks the alternating syn-anti pattern and costs an extra 4 kcal/mol to form (see Table 1). So the onset of transition in Sequence 2 is expected to occur at a more extreme superhelix density because it has a higher total B-Z transition energy. Transition at Sequence 3 has two disadvantages relative to the others. It is shorter in length, containing 13 Z-forming dinucleotides instead of 16, and it has two energetically costly “imperfections”. The GA and TC nucleotides in the anti-syn conformation cost 3.4 kcal/mol each, resulting in an additional transition cost of 6.8 kcal/mol relative to a perfect alternating CG sequence. Therefore, the critical superhelix density for complete Z-formation of this region is expected to be substantially higher than for either Sequence 1 or 2.

Fig. 4 shows the probabilities predicted by SIBZ of the inserted sequences transforming to Z-DNA in each of these three pBR322 derivatives, plotted as a function of superhelical density . The curves labeled Sequence 1, Sequence 2 and Sequence 3 refer to the plasmids with those as their inserts.

Experimentally measurements have been made of the critical superhelical density required to flip the entire inserted sequence into Z-DNA [37], [56]. These analyses regarded the transition as two-state, analogous to an “on-off” switch, and determined the critical superhelicity at which these inserted segments were switched on. These critical densities are shown in column 3 of Table 2.

It is difficult to make exact comparisons between our theoretical predictions and these experimental results because they do not involve entirely comparable quantities. The experiments measure a critical for completion of an “on-off”, two-state B-Z transition, while SIBZ calculates the equilibrium probability of the plasmid experiencing any amount of transition. However, the trend observed in the SIBZ results of Fig. 4 closely agrees with experimental data. Transition occurs at the least extreme superhelix density in the plasmid with the Sequence 1 insert, later in the plasmid with Sequence 2 insert, and lastly in the plasmid with the Sequence 3 insert. Moreover, the horizontal displacement between these curves is about 0.0025 between Sequence 1 and Sequence 2, and about 0.008 between Sequence 2 and Sequence 3. These agree closely with the experimentally measured differences.

The first experiment used anti-Z antibodies to isolate regions of Z-DNA formation [29]. They found three restriction fragments from the c-myc gene region that showed antibody reactivity when the gene was being transcribed. These three fragments are all located in the upstream region of the gene in proximity to its three promoters. Fig. 5 shows the SIBZ transition profile of a 5kb region around the c-myc gene that spans these locations. This profile was calculated at superhelical density , a reasonable value for transcriptionally driven superhelicity [6], [7]. The three large, gray horizontal bars labeled Z1, Z2, and Z3 identify the segments that were found in this experiment to contain Z-forming regions [29]. Five of the six sites identified by SIBZ as having the highest Z-forming probabilities occur within these three segments. The peak that is not within any of the three experimentally identified fragments is located around position 2100 in Fig. 5. Perhaps this region was missed because of its proximity to the Z2 fragment.

We are grateful to Sally Madden for her assistance with the computational aspects of this research.