Page 187

identity of the base pair. In practice, a >> b_GC> _bAT.Because the initiation free energy a is large, low-energy states tend to have only a small number of runs of strand separation. Because bGC is larger than b_AT, these runs tend to be in A+T-rich regions. The free energy of interstrand twisting within separated regions is quadratic in the local helicity of the deformation, with coefficient denoted by C. The free energy of residual superhelicity has been measured experimentally to be quadratic in that deformation, with coefficient K. Combining these contributions (and allowing the interstrand twisting to equilibrate with the residual superhelicity), the free energy G of a state is found to depend on three parameters: the number n of separated base pairs, the number n_ATof these that are A·Ts, and the number _{r of runs of separation:}

                  (7.6)

The energy parameters in this expression, a, b_AT, b_GC, C, and K, all depend on environmental conditions such as salt concentration and temperature. The values of the b's are known experimentally under a wide variety of conditions (Marmur and Doty, 1962; Schildkraut and Lifson, 1968). However, values for the other parameters are not so well understood. These parameters must be evaluated before the methods can yield quantitatively accurate results. We will do this by fitting these parameters to actual experimental data.

Analysis of Superhelical Equilibria

To calculate the equilibrium strand separation behavior of superhelical DNA molecules, we proceed as follows. First, the DNA sequence is analyzed and key information needed for later stages is stored. This step need be done only once per sequence. Next, the linking difference a and environmental conditions are specified, which sets the energy parameters and determines the free energy associated with each state. The state having minimum free energy under the given conditions

Page 188

is found from the free energy expression and the sequence data. Then an energy threshold q is specified, and all states i are found that have free energy exceeding the minimum G_min by no more than this threshold amount. Three inequalities occur, one each for n, n_AT, and r. Together the satisfaction of all three inequalities provides necessary and sufficient conditions that a state satisfy the energy threshold condition. For every set of values n, n_AT, and r satisfying these inequalities, all states with these values are found from the sequence information. This is a very complex computational task. The number of states involved grows approximately exponentially with the threshold q. In cases where r > 1, care must be taken to verify that a collection of r runs having the requisite total length and A+T-richness neither overlap nor abut, but rather are distinct. An approximate partition function Z_cal is computed from this collection of low-energy states to be

                                        (7.7)

By focusing only on the low-energy states, approximate ensemble average (that is, equilibrium) values are computed for all parameters of interest. These may include the expected torsional deformation of the strand-separated regions, expected numbers of separated base pairs, of separated A·T pairs, and of runs of separation, the ensemble average free energy , and the residual superhelicity.

The most informative quantities regarding the behavior of the molecule are its destabilization and transition profiles. The transition profile displays the probability of separation of each base pair in the molecule. The separation probability p(x) of the base pair at position x is calculated from equation (7.5) using parameter z_x, where z_x = 1 in states where base pair x is separated and z_x = 0 in all other states. This calculation is performed for every base pair in the sequence. The transition profile displays p(x) as a function of x.

The destabilization profile is the incremental free energy needed to induce separation at each base pair. To calculate this quantity, let i(x) index the states in which the base pair at position x is separated. Then the average free energy of all such states is

Page 189

                               (7.8)

To determine the destabilization free energy G(x), we normalize by subtracting the calculated equilibrium free energy:

                                                    (7.9)

Base pairs that require incremental free energy to separate at equilibrium have G(x) > 0, while base pairs that are energetically favored to separate at equilibrium have G(x) £ 0. This calculation is performed for each base pair in the molecule, and the destabilization profile plots G(x) versus x. Examples of these profiles are given in Figure 7.1.

Although individual high-energy states are exponentially less populated than low-energy states at equilibrium, they are so numerous that their cumulative contribution to the equilibrium still may be significant. The next step in this calculation requires estimating the aggregate influence of the states that were excluded from the above analysis because their free energies exceeded the threshold. This involves estimating the contribution Z(n,n_AT,r) to the partition function from all states whose values n,n_AT, and r do not satisfy the threshold condition. Here

                (7.10)

where M(n,r) is the number of states with n separated base pairs in r runs, which for a circular domain is

                                         (7.11)

The only part of Z(n,n_AT,r) not amenable to exact determination is P_n,r(n_AT), the fraction of (n,r)-states that have exactly n_ATseparated A·T base pairs.

Page 190

Figure 7.1
The transition (top) and the helix destabilization (bottom) profiles are shown for the
circular pBR322 DNA molecule at  a = -26 turns. The promoter (P) and terminator
(T) (the sites that control the start and end, respectively, of gene expression) of the
beta-lactamase gene are indicated. These results were calculated using the
energetics found by the method described in the text. Reprinted (bottom
panel), by permission, from Benham (1993). Copyright © 1993 by the
National Academy of Sciences.

The estimation of the influence of the high-energy states is done in two steps. First, growth conditions are found for Z(n,r), the contribution to the partition function from all states having n total separated base pairs in r runs. Assuming that the distribution of A+T-richness among (n,r)-states is approximately the same as that for (n,r + 1)-states, then the ratio

                             (7.12)

Page 191

is monotonically decreasing with r. One can find an such that T=p<1. Then for all r ³ we have Tr £r, so that

                          (7.13)

By a similar line of reasoning one also can find a value , above which the aggregate contribution to the partition function again is bounded above by a convergent geometric series:

                             (7.14)

In practice, low values ( » 0.1) for the series ratios r and s occur at reasonably small cutoffs ( » 8, » 150 for a molecule of N = 5,000 base pairs under reasonable environmental conditions).

The contribution of the intermediate states having n £ and r £ but not satisfying the threshold requires estimating p_n,rn_AT, the fraction of n,r-states that have exactly n_ATseparated A·T base pairs. Although in principle one can compute this quantity from the base sequence, for molecules of kilobase lengths it is feasible to compute only the exact distribution of A+T-richness in _{r=I run states having n £ nmax} » 200. Experience has shown that high accuracy is obtained by calculating p_n,2n_AT exactly, and using p_n,2n_AT as an estimate of p_n,rn_AT for r> 2. Once the sequence information needed in this step has been found (a calculation that need be performed only once per molecule), the performance of the rest of this refinement is computationally very fast.

These results are used to estimate the contribution _neg to the partition function from the neglected, high-energy states:

                                      (7.15)

(Here the carat marks denote approximate values.) Any parameter z that depends only on n, n_AT, or r also can have its previously calculated