One-Loop Fluctuation Entropy of Charge Inversion in DNA

TL;DR

This paper models charge inversion in DNA using a lattice-gas approach, incorporating fluctuation corrections to mean-field theory, revealing significant effects on entropy especially when dimers are repelled from DNA.

Contribution

It introduces a lattice-gas model with fluctuation corrections to better understand charge inversion and entropy in DNA solutions, advancing beyond mean-field approximations.

Findings

Fluctuation corrections significantly affect entropy calculations.

Charge distributions depend on dimer orientation and concentration.

Model captures correlation-driven charge inversion phenomena.

Abstract

Experiments have revealed correlation-driven behavior of DNA in charged solutions, including charge inversion and condensation. This paper presents calculations of a lattice-gas model of charge inversion for the adsorption of charged dimers on DNA. Each adsorption site is assumed to have either a vacancy or a positively-charged dimer attached with the dimer oriented either parallel or perpendicular to the double helix DNA chain. The entropy and charge distributions of these three species are calculated including the lowest order fluctuation corrections to mean-field theory. We find that the inclusion of the fluctuation terms has a significant effect on the entropy, primarily in the regime where the dimers are repelled from the DNA molecule and compete with the chemical potential in solution.

Equations

$${\mathcal{H}}_{\text{NI}}=\sum _{\ell ,b}\sum _{\gamma }{\epsilon }_{\gamma }{n}_{\ell ,b}^{(\gamma )}.$$

$${\Omega }_G(T,V,\mu )=E-TS-\sum _{\gamma }{\mu }_{\gamma }{N}_{\gamma },$$

$$N_{\gamma }=\sum _{\ell ,b}{n}_{\ell ,b}^{(\gamma )}.$$

$$S=-{\left(\frac{\partial {\Omega }_G}{\partial T}\right)}_{V,\{{\mu }_{\gamma }\}}=k_B{\beta }^2{\left(\frac{\partial {\Omega }_G}{\partial \beta }\right)}_{V,\{{\mu }_{\gamma }\}},$$

$$Z_G=\sum _{\text{configurations}}{e}^{-\beta (E-{\sum }_{\gamma }{\mu }_{\gamma }{N}_{\gamma })}=e^{-\beta {\Omega }_G},$$

$${\Omega }_G=-\frac{\ln Z_G}{\beta },$$

$$S=k_B\ln Z_G-\frac{k_B\beta }{Z_G}{\left(\frac{\partial Z_G}{\partial \beta }\right)}_{V,\{{\mu }_{\gamma }\}}.$$

$$Z_{\text{NI}}=\sum _{\text{configurations}}{e}^{-\beta {\sum }_{\ell ,b}{\sum }_{\gamma }\left({\epsilon }_{\gamma }-{\mu }_{\gamma }\right)n_{\ell ,b}^{(\gamma )}},$$

$$⟨n^{(\gamma )}{⟩}_{\text{NI}}=\frac{1}{{\mathcal{N}}_{\text{sites}}\beta Z_{\text{NI}}}\frac{\partial Z_{\text{NI}}}{\partial {\mu }_{\gamma }}.$$

$$⟨n^{(\gamma )}{⟩}_{\text{NI}}$$

$$Z_{\text{NI}}=\sum _{\text{configurations}}\prod _{\ell ,b}\prod _{\gamma }{\left[e^{-\beta \left({\epsilon }_{\gamma }-{\mu }_{\gamma }\right)}\right]}^{n_{\ell ,b}^{(\gamma )}}.$$

$$a_{\text{NI}}^{(\gamma )}=e^{-\beta \left({\epsilon }_{\gamma }-{\mu }_{\gamma }\right)}.$$

$$Z_{\text{NI}}=\sum _{\text{configurations}}\prod _{\gamma }\prod _{\ell ,b}{\left[a_{\text{NI}}^{(\gamma )}\right]}^{n_{\ell ,b}^{(\gamma )}}.$$

$$\{n_{\ell ,b}^{(\perp )},n_{\ell ,b}^{(\parallel )},n_{\ell ,b}^{(v)}\}=\{1,0,0\},\{0,1,0\},\text{mbox}or\{0,0,1\}.$$

$$\sum _{\begin{array}{c}\text{single-site}\\ \text{configurations}\end{array}}\prod _{\gamma =\perp ,\parallel ,v}{\left[a_{\text{NI}}^{(\gamma )}\right]}^{n_{\ell ,b}^{(\gamma )}}=\sum _{\gamma }{a}_{\text{NI}}^{(\gamma )},$$

$$Z_{\text{NI}}=\prod _{\ell ,b}\left(\sum _{\gamma }{a}_{\text{NI}}^{(\gamma )}\right).$$

$$Z_{\text{NI}}={\left(\sum _{\gamma }{a}_{\text{NI}}^{(\gamma )}\right)}^{{\mathcal{N}}_{\text{sites}}}.$$

$$⟨n^{(\gamma )}{⟩}_{\text{NI}}=\frac{1}{{\mathcal{N}}_{\text{sites}}\beta Z_{\text{NI}}}\frac{\partial }{\partial {\mu }_{\gamma }}{\left(\sum _{{\gamma }^{\prime }}{a}_{\text{NI}}^{({\gamma }^{\prime })}\right)}^{{\mathcal{N}}_{\text{sites}}}.$$

$$⟨n^{(\gamma )}{⟩}_{\text{NI}}$$

$$\frac{\partial }{\partial {\mu }_{\gamma }}{a}_{\text{NI}}^{(\gamma )}=\frac{\partial }{\partial {\mu }_{\gamma }}{e}^{-\beta \left({\epsilon }_{\gamma }-{\mu }_{\gamma }\right)}=\beta a_{\text{NI}}^{(\gamma )}.$$

$$⟨n^{(\gamma )}{⟩}_{\text{NI}}=\frac{a_{\text{NI}}^{(\gamma )}}{{\sum }_{{\gamma }^{\prime }}{a}_{\text{NI}}^{({\gamma }^{\prime })}}.$$

$${\rho }_{\ell b}=\sum _{\gamma }{q}_{\gamma }{n}_{\ell b}^{(\gamma )}.$$

$${\rho }_{\text{NI}}\equiv ⟨{\rho }_{\ell b}{⟩}_{\text{NI}}=\sum _{\gamma }{q}_{\gamma }⟨n^{(\gamma )}{⟩}_{\text{NI}}.$$

$${\sigma }_{\text{NI}}^2=⟨{\rho }^2{⟩}_{\text{NI}}-{\rho }_{\text{NI}}^2,$$

$$⟨{\rho }^2{⟩}_{\text{NI}}=⟨{\rho }_{\ell b}^2{⟩}_{\text{NI}},$$

$${\rho }_{\ell b}^2={\left[\sum _{\gamma }{q}_{\gamma }{n}_{\ell b}^{(\gamma )}\right]}^2=\sum _{\gamma ,{\gamma }^{\prime }}{q}_{\gamma }{q}_{{\gamma }^{\prime }}{n}_{\ell b}^{(\gamma )}{n}_{\ell b}^{({\gamma }^{\prime })}.$$

$$[n_{\ell b}^{(\gamma )}{]}^2=n_{\ell b}^{(\gamma )}.$$

$${\rho }_{\ell b,\text{NI}}^2=\sum _{\gamma }{q}_{\gamma }^2{n}_{\ell b,\text{NI}}^{(\gamma )}.$$

$$⟨{\rho }^2{⟩}_{\text{NI}}\equiv ⟨{\rho }_{\ell b}^2{⟩}_{\text{NI}}=\sum _{\gamma }{q}_{\gamma }^2⟨n^{(\gamma )}{⟩}_{\text{NI}}.$$

$$S_{\text{NI}}$$

Full text

I Introduction

Polyelectrolytes in solution have been shown to produce a wide array of intriguing phenomena in many different systems. One example is DNA, which in the presence of a physiological salt solution (a 0.1 molar solution of NaCl) is usually negatively charged, with a double helix DNA molecule strongly repelling another. However, if DNA is in a dilute salt solution in which a positively charged polyelectrolyte, such as spermine or sperimidine, has been added to the solution, it has been shown to roll up into a tightly-packed torusBloomfield (1991, 1997); Raspaud et al. (1999); Gelbart et al. (2000). In fact, DNA is usually in a very compact state in cells and viruses.

Another example is F-actin, whose fibrils stack in cross-linked rafts when positive alkaline earth ions are in the solutionWong et al. (2003). Also in F-actin, counter-ions form one-dimensional charge density waves that have a periodicity equal to twice the actin monomer spacing, coupling to twist distortions of the oppositely-charged actin filamentsAngelini et al. (2003). These phenomena and others seem to owe their existence to Coulomb interactions between the constituent parts, and the geometry of the underlying structure may also play a vital role.

Determination of the optimal structure of any of these systems requires minimization of the free energy, which involves a competition between the internal energy and entropy. Often in biological systems the changes in entropy are greater than the changes in internal energy. For instance, in the hydrophobic effect, an unfolded protein lowers the entropy by ordering the water molecules, and so the protein prefers to be in the folded stateUdgaonkar (2001); Lodish et al. (2004). In the chloroplast stroma, it has been shown that there is an entropy-driven attraction that determines chloroplast ultrastructure through spontaneous Mg2+-induced stacking of membranesYoshidome (2010). In sickle hemoglobin, the aggregation of monomers into polymers is also entropy driven, with the internal energy and entropy in a delicate balanceCao and Ferrone (1997).

Historically, the most common description of charged solutions has been Poisson-Boltzmann theory, but continued research has indicated that ions in solution can have far more complex and counterintuitive effects than simple charge screening Gelbart et al. (2000). In 1969 Manning proposed Manning (1969a, b) that a portion of the ambient ions condense onto (i.e., attach to) the surface of a charged macro-ion, partially neutralizing the bare charge. This occurs up to the point at which the energy required to condense another ion equals the available thermal energy $k_{B}T$ at temperature $T$, where $k_{B}$ is the Boltzmann constant. The proposed effect, later termed “Manning condensation,” marked a significant departure from linearized Poisson-Boltzmann theory that predicted only exponential screening. In Manning’s treatment of polyelectrolytes, the long chains of charged subunits, the ion condensation was addressed Gelbart et al. (2000) separately from the ions that remain in solution, which were treated using linearized Poisson-Boltzmann theory Manning (1969a, b).

Further refinements in the treatment of these ambient-ion solutions have been motivated in part by surprising effects observed in DNA that cannot be accounted for using Poisson-Boltzmann theory. By increasing the charge on the salt-ions in the DNA solution from $+e$ (monovalent) to $+2e$ (divalent) and higher, it is possible to cause the DNA to undergo a radical structural transition into a variety of new geometries, including rod-like bundles and toroids Bloomfield (1996). The toroidal structures, in particular, have received considerable attention in the biological community because such toroidal packing motifs are employed by spermidine and other molecules to contain their own DNA in small volumes Hud and Vilfan (2005). A variety of techniques, including cryo-transmission electron microscopy Hud and Vilfan (2005) and UV spectroscopy Raspaud et al. (2006), have enabled direct observation of the formation of these toroidal condensates, and similar studies on other biological polyelectrolytes like F-actin Wong et al. (2003) indicate that condensation is a general phenomenon. If the electrical interactions between the polyelectrolyte, such as DNA, and the ambient-ion solution are purely those predicted by the Poisson-Boltzmann description, then two like-charged polyelectrolytes always repel one another, and such condensates will not form. For these condensates to be stable, the net electrical force between like-charged polyelectrolytes must be $\it attractive$, which is incompatible with the Poisson-Boltzmann theory of ions in solution.

Furthermore, under the same conditions of divalent solution-ions, a radical change in the electrical properties of a single DNA molecule is observed. In gel electrophoresis experiments Raspaud et al. (2006) under these conditions, DNA is observed to drift toward the $\it negative$ electrode under the influence of an electric field, indicating that the net charge is $\it positive$. This difference is quantified as a change in sign of the electrophoretic mobility, which occurs for some minimum concentration of solution-ions with valence $\geq 2$. This change in sign cannot be explained by the Poisson-Boltzmann description of charge screening, and, although Manning condensation can reduce the effective charge, it cannot invert it Grosberg et al. (2002). As with polyelectrolyte condensation, charge inversion requires a fuller explanation of the ambient-ion solution than Poisson-Boltzmann theory can provide. Understanding the nature of these phenomena is important not only for the insight into the natural role of DNA in a cell, but also for use in medical techniques, such as gene therapy Luo and Saltzman (2000).

Many of the approaches to understanding the role of polyelectrolytes in biological systems have adopted models of continuum electrostatics, and some of these have focussed on solving the Poisson-Boltzmann equation for a cylindrical model of DNA in the presence of multivalent ionsGelbart et al. (2000). However, it has been recognized that including charge correlations is essential to understanding these systemsShklovskii (1999); Grosberg et al. (2002). Nguyen and Shklovskii Nguyen and Shklovskii (2002a, b) proposed a simple theory of charge inversion that considers the geometrical structure of the polyelectrolyte together with its electrostatic interaction with the substrate. The basic theory could be applied to arbitrary-length polyelectrolytes adsorbed on various surfaces from linear biomolecules to biomembranes. The idea was that there is fractionalization of charge in which a polyelectrolyte can either neutralize the charge or reverse it depending on how it attaches. In that model, a single double-helix strand of DNA was represented by a rigid cylinder with two one-dimensional lattices of negative charges $-e$ in helices around the surface to represent DNA’s double helix of negatively-charged phosphates, as shown in Fig. 1. Since the DNA is modeled to be in a polyelectrolyte solution, each strand is surrounded by positively-charged species that are long and possibly flexible. These adsorbing species have multiple charges that may partially attach to the surface, with excess charge protruding into the solution, as shown in Fig. 2. By this means, since there is excess charge dangling into the solution, the charge on the surface can be not only neutralized but reversed.

The results of these generalized theories indicate that the correlation effects in the solution are also strongly-dependent on the system geometry. Thermodynamically, the distribution of charge in the solution is governed by a competition between energy and entropy. The bound state in which the ions are condensed on the surface of the macroion/polyelectrolyte is energetically-favored, and the continuum state, in which the ions are free to drift in any direction, is entropically-favored. The balance between these two factors in minimizing the free energy has been shown to vary significantly based on the geometry of the macroion Gelbart et al. (2000).

To understand this, consider the electric potential outside spherical, cylindrical, and planar charged surfaces in vacuum. For spheres, the potential decays as the inverse of the distance; for the cylinder, the potential varies as the logarithm of the distance; and for the plane, the potential increases linearly with the distance. Gelbart et al. Gelbart et al. (2000) assert that, in an ambient-ion solution, the entropy of the point ions varies logarithmically with their concentration, and therefore with the distance from the cylinder. Thus they conclude that, in a crude comparison, for spherical geometry, the $r^{-1}$ potential is dominated by the logarithmic entropy; for planar geometry, the entropy is dominated by the linear potential; and for cylindrical geometry, both the energy and the entropy vary logarithmically.Sharp (1995) In this paper, we will incorporate the impact of the ions in solution on the free energy through the chemical potential of these ions, as detailed in Appendix B.

The helical geometry of DNA, then, sits precisely balanced on the fulcrum between energy- and entropy-driven processes under physiological conditions. The geometry of the biomolecule plays an even more crucial role when the highly-charged ions in solution are not point charges but have geometries of their own. Such is the case for DNA immersed in a solution of polyelectrolytes. These polyelectrolytes can be proteins or other fragments of biomolecules which are routinely found in the nucleus Lehninger (1982), so that, again, the central biological processes occur in precisely the most difficult regimes to model.

A special case of the the work of Nguyen and ShklovskiiNguyen and Shklovskii (2002a) studied by Bishop and McMullenBishop and McMullen (2006) considers the case of dimers (polyelectrolytes that are two units in length) on DNA. The choice to model dimers is motivated by the wealth of studies (refs. Fisher (1966) and Adams and Chandrasekharan (2003), among others) in the literature about dimer models in other branches of physics, as well as for geometrical simplicity. Dimers, the shortest polyelectrolytes, have only two possible orientations when adsorbed on DNA, assuming that the length of the dimer is comparable to the helical spacing between charged sites on DNA, as shown in Fig. 3. A dimer lying on the cylindrical surface of the molecule must lie parallel to the helical strands, neutralizing the charge on two adjacent sites. Otherwise, the dimer must adsorb perpendicular to the cylindrical surface, extending one end radially out from the central axis of the cylinder and inverting the charge on a single site. Charge inversion by dimers, then, is quite similar to charge inversion due to two species of point ions: one monovalent species representing parallel-adsorbed dimers and one divalent species representing perpendicular-adsorbed dimers. In this way, Bishop and McMullenBishop and McMullen (2006) modeled the adsorption of dimers in a lattice-gas model as a two-component solution of point ions, allowing the possibility of vacancies, and used field-theoretic methods to describe the thermodynamics of the system. They carried their calculations to a mean-field level of approximation of the inverted charge on DNA, but did not calculate the entropy. Their work confirmed the possibility of charge inversion within this model.

For charge inversion on DNA with polyelectrolytes, “physiological conditions” require incorporating the combined effects of charge correlation, thermodynamic fluctuations, crowding, and geometrical considerations all at once. As we have discussed in this introduction, there has been considerable work in addressing all of these issues. Some approaches treat only the geometry, with no interactionsMaltsev et al. (2006). Others include both geometry and interactions, but use a continuum model that neglects discreteness effectsKunze and Netz (2002). The lattice gas dimer model is unique in its simultaneous consideration of all these effects, and the thermodynamic and geometrical idealizations it does make can be systematically improved.

While the formalism of Bishop and McMullenBishop and McMullen (2006) included both a mean field theory and corrections due to fluctuations or charge correlations, the computed results were obtained only at the mean field level. Those computations gave the charge per site on the DNA helix as a function of the chemical potential, or equivalently the concentration of the polyelectrolyte in solution. While the Nguyen-ShklovskiiNguyen and Shklovskii (2002a) calculation assumed complete filling of the lattice, the Bishop-McMullen approach allowed for vacancies, represented as negatively-charged sites, in addition to the neutral or positively-charged sites arising from dimer adsorption. However, the lattice-gas model, which assumes all sites are equivalent, does not take into account that the parallel dimer occupies two sites. Instead, the occupation of two sites is mimicked by making the binding energy of the parallel dimer twice as large as that of the perpendicular dimer, and we will use the same approach here. In this paper, this simple model is extended to calculate the fluctuation corrections and the entropy of the system in this simple model. The purpose is to determine the importance of the fluctuation terms for inverting the charge and to see whether these terms have a significant impact on the entropy of the system. A preliminary version of this work was in Sievert’s master’s thesisSievert (2007). This paper does not attempt to include all the possible effects of coions, counterions, and the nonzero hard-core radius of the background ions of the solution, as was done in the work by Solis and Olvera de la Cruz Solis and Olvera de la Cruz (2001); Solis (2002). We also do not attempt to study complicated geometric effects of multivalent polyelectrolyte macroions, such as the “ion-bridging” model of Olvera de la Cruz Olvera de la Cruz et al. (1995), which is used as an explanation of the data of Raspaud et al. Raspaud et al. (1998) on spermine-induced aggregation of DNA, nucleosome, and chromatin. However, this approach could certainly be extended to include all of those possibilities, and in this model, some of these effects can considered to be contained in the binding energies of the polyelectrolytes to the surface and the screened Coulomb interaction between charged species.

In Sec. II of this paper, we outline the model of the charged binding sites on DNA and present the computation of the entropy due only to the hard-core repulsion that prevents multiple binding on the same site. In Sec. III we explain the geometry of the double helix and show how it can be represented as a one-dimensional lattice and in Sec. IV derive the form of the potential and determine the orthogonal transformation that diagonalizes it. In Sec. V, we use functional integral techniques to derive the partition function, which uses a Gaussian integral identity to perform the sum over configurations exactly, at the expense of an integral over a new auxiliary field. In Sec VI, we show how the grand canonical potential can be computed order by order, in which the first two terms are the mean field and the one-loop correction to the mean field. In Sec. VII, we examine the saddle point equation that defines the mean field, and this is used to calculate the entropy, the number densities of all species, and the charge density. In Sec. VIII, we find the explicit form for the inverse-Hartree-field-fluctuation propagator. In Sec. IX, we include the one-loop order terms in the calculations to reveal the effects of fluctuations. Finally, in Sec. X, we compare the results of the various approximations. Details concerning the Gaussian integral identity and the chemical potential of dimers in solution are given in Appendices A and B.

II The Noninteracting DNA Lattice Gas

In this section, we will derive the entropy and the average occupations for the three species of vacancies, parallel dimers, and perpendicular dimers for the lattice-gas model of DNA in the absence of electrostatic interactions between the species on different sites of the lattice. Although this could be done using combinatorics, we will use an approach analogous to what be used later when interactions are included, and this will make the procedure for the interacting case clearer. In the lattice-gas model that we consider here, one feature of the adsorption of species must necessarily be neglected. When a dimer is adsorbed parallel to a phosphate chain it occupies two adjacent sites, but the lattice-gas model treats all sites independently and cannot take into account the blocking of an adjacent site by a parallel dimer. This would be an even more complicated issue for longer polymers, and is generally difficult to incorporate. As a first step toward building a comprehensive model of such systems, here we neglect the adjacent blocking effect for parallel adsorption. Nevertheless, we retain the different binding energies $\varepsilon^{(\parallel)}$ and $\varepsilon^{(\perp)}$, since the two-site occupation implies $\varepsilon^{(\parallel)}\approx 2\varepsilon^{(\perp)}$. This allows us to describe adsorption of either species independently for each site. Any configuration of the DNA-dimer complex in this model can then be described by identifying the type of dimer (parallel, perpendicular, or none) adsorbed on each site on the DNA molecule. Such a model resembles the lattice gas model of condensed matter physics Plischke and Bergersen (1989), which treats the ways of distributing particles of different types onto a regular lattice of sites. If the particles on different lattice sites do not interact with one another, then the structure of the lattice does not matter.

It is convenient to consider a vacancy as a third species $\gamma$ of particle, because then we can impose the constraint that each site is singly occupied, either by a parallel dimer $(\gamma=\>\parallel)$, a perpendicular dimer ($\gamma=\>\perp$), or a vacancy ($\gamma=v$). For each our three species of “particles” that can reside on our lattice sites, we define a relative charge $q_{\gamma}$ in units of the magnitude $e$ of the charge of an electron. These relative charges are then $q_{v}=-1$ for vacancies, $q_{\parallel}=0$ for parallel dimers, and $q_{\perp}=+1$ for perpendicular dimers. We will assume that the binding energy $\varepsilon_{\gamma}$ of species $\gamma$ to the lattice depends only on the type of species and not the location. We will be specifying each lattice point by its location $\ell$ along helix $b$, with the pair $(\ell,b)$ specifying that lattice position. Then the quantity $n_{\ell,b}^{\gamma}$ will be the number of particles of species $\gamma$ on the lattice site at $\ell$ on chain $b$. For each chain, $\ell$ extends from $-{\mathcal{N}}$ to ${\mathcal{N}}$, such that the total number of sites is ${\mathcal{N}}_{\text{sites}}=2(2{\mathcal{N}}+1)$, and we employ periodic boundary conditions.

We begin by considering the Hamiltonian ${\mathcal{H}}_{\text{NI}}$ for this noninteracting system (that is, with no interactions between different sites aside from the hard-core repulsion that blocks double-occupancy), which is

[TABLE]

Because we have a system that exchanges particles with its surroundings, specifically the ions in the solution surrounding the DNA, which attach to the surface, we work in the grand canonical ensemble. If the system contains $N_{\gamma}$ particles of type $\gamma$ in equilibrium with its surroundings, with an average internal energy $E$, the grand canonical potential $\Omega_{G}$, which is a function of the temperature $T$, volume $V$, and chemical potentials $\mu_{\gamma}$ for each species $\gamma$, is written asPlischke and Bergersen (1989); Reif (1965),

[TABLE]

where the number of particles of type $\gamma$ is given by

[TABLE]

Then the entropy can be written in terms this potential as

[TABLE]

where $\beta=\frac{1}{k_{B}T}$ and $k_{B}$ is Boltzmann’s constant. In this partial derivative, the volume and all the chemical potentials $\mu_{\gamma}$ of all species $\gamma$ are held constant. It is convenient to define the grand canonical partition function $Z_{G}$ in terms of the grand canonical potential $\Omega_{G}$ as

[TABLE]

where the sum over configurations includes all the accessible states of the system. The grand canonical potential can alternatively be written in terms of the partition function as

[TABLE]

and this allows us to write an expression for the entropy in terms of $Z_{G}$ as:

[TABLE]

For the noninteracting lattice, the average internal energy is represented here by the Hamiltonian (1), $E={\mathcal{H}}_{NI}$, and the noninteracting grand canonical partition function $Z_{\text{NI}}$ becomes

[TABLE]

where we have suppressed the $G$ subscript for simplicity.

Since a vacant site does not really correspond to a particle, we recognize that energy and chemical potential of the vacancy must be related such that $\varepsilon^{(v)}-\mu_{v}=0$. In addition, because the dimers adsorbed on the surface and those in solution are in equilibrium, $\mu_{\parallel}=\mu_{\perp}=\mu_{\text{dimer}}$ is the chemical potential of dimers in solution at the appropriate concentration. We thus need an estimate for $\mu_{\text{dimer}}$. Approximating the dimer as a uniformly-charged cylinder, this value is shown in Appendix B to be $\beta\mu_{\text{dimer}}\simeq 0.79$ at the physiological temperature of $T\simeq 310$ K.

The average occupancy $\langle n^{(\gamma)}\rangle_{\text{NI}}$ for species $\gamma$ per site is found by taking the derivative of this partition function with respect to $\mu_{\gamma},$

[TABLE]

This can be verified by taking the derivative of Eq. (8).

[TABLE]

which is by definition the average number per site of species $\gamma$.

Continuing with the evaluation of the partition function, we note that the exponential of a sum can be written as the product of exponentials, enabling us to rewrite the partition function (8) as

[TABLE]

In order to simplify and to appreciate the physical meaning of this expression, it is useful to define the relative activity Moore (1972) of species $\gamma$ in the noninteracting lattice-gas model, given by

[TABLE]

Note that this is independent of the lattice site $\ell$ or chain $b$. The partition function can then be written in the simple form

[TABLE]

Because we have assumed that there can be only one species per site, parallel dimer, perpendicular dimer, or vacancy, the sum over configurations can now be done over each site separately, where there are three possible configurations

[TABLE]

This is the same as saying that $n_{\ell,b}^{(\gamma)}=1$ for one and only one of $\gamma=\perp,\parallel$, or $v$ and is zero otherwise. This means that

[TABLE]

and so the grand partition function is

[TABLE]

Since every term in the product is the same, the expression in parentheses is simply raised to the power ${\mathcal{N}}_{\text{sites}}$, giving

[TABLE]

At this point, we can easily see that we can obtain the average number per site using our derivative form from Eq. (9) as

[TABLE]

Taking the derivative and substituting the expression for $Z_{\text{NI}}$ in the denominator, we have

[TABLE]

where the derivative of the activity is

[TABLE]

Using this in the expression for the mean site occupancy in Eq. (19) and cancelling ${\mathcal{N}}_{\text{sites}}-1$ factors of the sum over $\gamma$, the mean occupancy for species $\gamma$ is given by

[TABLE]

In Fig. 4, we show the mean occupancies for the three species versus $\varepsilon_{\parallel}$, with $\varepsilon_{\parallel}=2\varepsilon_{\perp}$. Negative $\varepsilon_{\parallel}$ corresponds to an attraction of dimers to the lattice, and we see that parallel adsorption of dimers dominates because of the stronger binding, followed by perpendicular adsorption, with vacancies becoming nonexistent. Positive $\varepsilon_{\parallel}$ corresponds to a repulsion of dimers from the lattice, so that at large $\varepsilon_{\parallel}$, the ordering is reversed, for the same reasons. At small positive $\varepsilon_{\parallel}$, the dimers are repelled from the lattice, but this competes with the chemical potential, which tries to put dimers back onto the lattice. In this low-coverage regime, perpendicular adsorption dominates, while vacancies dominate for large $\beta\varepsilon_{\parallel}$ because the dimers would prefer to stay in solution. For $\varepsilon_{\parallel},\varepsilon_{\perp}\rightarrow 0$, the parallel and perpendicular mean occupancies become the same. Similarly, where $\varepsilon_{\parallel}-\mu_{\text{dimer}}=0$, $\langle n^{(\parallel)}\rangle_{\text{NI}}=\langle n^{(v)}\rangle_{\text{NI}}$ because, as mentioned earlier, $\varepsilon_{v}-\mu_{\text{dimer}}$ is always zero. Also, $\langle n^{(\perp)}\rangle_{\text{NI}}=\langle n^{(v)}\rangle_{\text{NI}}$ where $\varepsilon_{\parallel}=2\mu_{\text{dimer}}$ because that is where $\varepsilon_{\perp}=\mu_{\text{dimer}}$. It is at this point where $\langle n^{(\parallel)}\rangle_{\text{NI}}$ becomes a maximum.

The average charge per site can now be determined by multiplying $q_{\gamma}$, the charge of species $\gamma$, by $n_{\ell,b}^{(\gamma)}$ the occupation of species $\gamma$ on site $\ell$ of chain $b$,

[TABLE]

Since we assume that every site has either a parallel dimer, a perpendicular dimer, or a vacancy with single occupancy, then for a given site and species, $n_{\ell b}^{(\gamma)}$ is either 0 or 1. The average charge density is then the sum the averages of the individual terms,

[TABLE]

This charge density is plotted as the solid blue curve in Fig. 5. Because a parallel dimer has no charge, $q_{\parallel}=0$, a perpendicular dimer has a positive unit charge, $q_{\perp}=+1$, and a vacancy has a negative unit charge, $q_{v}=-1$, the total mean field charge per site $\rho_{c}$ in Fig. 5 is the difference between the curves $\langle n^{(\perp)}\rangle_{\text{NI}}$ and $\langle n^{v}\rangle_{\text{NI}}$ in Fig. 4 and goes to zero where those two curves cross.

A measure of the magnitude of the charge fluctuations is given by the average of the square of the charge on a site minus the square of its average, which is the charge variance $\sigma_{NI}^{2}$, given by

[TABLE]

where

[TABLE]

and

[TABLE]

The product $n_{\ell b}^{(\gamma)}n_{\ell b}^{(\gamma^{\prime})}$ describes a double occupancy of site $(\ell,b)$ by species $\gamma$ and $\gamma^{\prime}$. Since we have required single occupancy, then we must have $\gamma=\gamma^{\prime}$. Also, since $n_{\ell b}^{(\gamma)}$ can only take the values [math] or $1$,

[TABLE]

Then the charge density per site reduces to

[TABLE]

Taking the thermal average of both sides, we have

[TABLE]

In Fig. 5, we have plotted the standard deviation in the charge density, $\sigma_{\text{NI}}$, which is the square root of the charge variance, for the noninteracting model. As we saw in Fig. 4, there are three distinct regions, $\varepsilon_{\parallel}\ll 0$, $\varepsilon_{\parallel}\gg 0$ and $\varepsilon_{\parallel}$ near $\mu_{\text{dimer}}$. These three regions are also reflected in Fig. 5. At large negative $\varepsilon_{\parallel}$, parallel adsorption dominates, which leads to $\rho\approx 0$. At large positive $\varepsilon_{\parallel}$, vacancies dominate, leading to $\rho<0$, and for $\varepsilon_{\parallel}$ near $\mu_{\text{dimer}}$, there is a small window where perpendicular adsorption of dimers dominates, leading to positive values of $\rho$. Correspondingly, the fluctuations, represented by the standard deviation $\sigma_{NI}$, become small when $|\varepsilon_{\parallel}|$ becomes large, and the fluctuations are largest at small $|\varepsilon_{\parallel}|$, when there are comparable numbers of all species. The phenomenon of charge inversion is demonstrated in Fig. 5 because the average charge is positive, indicating that sufficiently many dimers adsorb in a perpendicular configuration to invert the charge on the molecule from negative to positive. The magnitude of charge inversion increases in the weak-binding limit $\varepsilon_{\perp},\varepsilon_{\parallel}\rightarrow 0$.

The noninteracting entropy $S_{\text{NI}}$ can be obtained from the partition function using Eq. (7) as

[TABLE]

Because the derivative of the activity with respect to $\beta$, can be written in terms of its logarithm as

[TABLE]

the entropy can be written in the simple form

[TABLE]

Substituting the mean occupancy for the noninteracting model from Eq. (21) allows us to write the dimensionless entropy per site for the noninteracting model in a simplified form as

[TABLE]

which agrees with the standard result for the entropy of mixing of an ideal solution with species $\gamma=(\parallel,\perp,v)$ Moore (1972).

In Fig. 6, we show the entropy $S_{NI}$ of the noninteracting model as a function of the binding energy $\beta\varepsilon_{\parallel}$, assuming $\varepsilon_{\parallel}=2\varepsilon_{\perp}$. We also show the individual contributions of Eq. (33) to the entropy. The entropy is a maximum when disorder is greatest, and this occurs when the numbers of each of the species are as close to equal as possible, which occurs at $\varepsilon_{\parallel}=\mu_{\text{dimer}}$, the maximum of $\langle n^{(\perp)}\rangle$ in Fig. 4.

III Geometry of the Charged Double Helix of DNA

While stored in the nucleus of a cell, DNA is wrapped compactly both around histone protein complexes and around itself, but on sufficently small scales ($\approx$ 150 base pairs Moukhtar et al. (2007), or 15 turns Lehninger (1982), or $50$ nm Gelbart et al. (2000)), DNA’s dominant geometrical structure is the familiar double-helix structure shown in Fig. 1.

Each strand of the DNA takes the shape of a helix with a characteristic radius $R_{\text{DNA}}=0.946$ nm and pitch angle $\psi=29.6^{\circ}$, as shown in Fig. 7. Here the pitch angle $\psi$ denotes the angle with respect to the $xy$-plane that gives the appropriate altitude per unit circumferential winding; in cylindrical coordinates, $\tan\psi=\Delta z/R_{\text{DNA}}\Delta\phi$, where $\Delta z$ is the distance along the $z$ axis. Each strand of DNA has a “direction”, identified by a particular carbon on the backbone structure, corresponding to the chirality of the helix. In the DNA double-helix, the two strands are antiparallel and therefore have opposite chiralities. As a consequence, the azimuthal angle between the two helices is always a constant, $\Delta\phi=160^{\circ}$. Because this phase shift is not exactly $180^{\circ}$, the chains have unequal separation in the clockwise and counterclockwise azimuthal directions. The larger gap is referred to as the major groove, and the smaller as the minor groove.

When the hydrogen atoms dissociate under physiological conditions, the resulting negative charges ($-e$, where $e$ is the magnitude of the charge of the electron), can be regarded as located at the mean position of the oxygen atoms on the backbone, as shown in Fig. 1. This figure was constructed using geometrical data taken from Bishop and McMullen Bishop and McMullen (2006), which is based on experimental x-ray diffraction data Arnott and Hukins (1969, 1972) as input to the SYBYL molecular modeling programSYBYL (2006). The oxygen atoms occur at regular intervals along each strand, separated by a helix segment of arc length $a^{(1)}=0.684$ nm. It is these sites located at regular intervals along the helix to which dimers will adsorb. These negatively-charged sites do not occur at the same altitudes on both strands, however. Rather, there is a vertical separation $\Delta z=0.023$ nm between corresponding sites on the two strands. With the relative phase of the strands and the vertical separation between sites on those strands taken together, corresponding sites on the two strands may be viewed as connected by a helical segment of arc length $a^{(2)}=3.34$ nm at a pitch angle of $\alpha=0.394^{\circ}$. This geometry is shown in Fig. 7.

When positively-charged dimers approach the DNA molecule, they will be attracted to the negative charges at the sites on the double-helix. We consider dimers with a length comparable to the spacing $a^{(1)}$ between sites on a strand and having positive charges $+e$ at either end. The dimers can then adsorb onto the surface in two possible orientations, either parallel to the strand or perpendicular to the helix axis. (See Fig. 3.) If the dimer adsorbs parallel to the strand, the positive charges from the dimer lie directly over the negative charges on the strand, neutralizing the charge on two adjacent sites. If the dimer adsorbs perpendicular to the surface of the bounding cylinder, one end of the dimer sits atop the site, while the other extends radially outward. This perpendicular adsorption effectively inverts the charge on the site from $-e$ to $+e$. In order to use a lattice-gas model, this geometric constraint is loosened by having the parallel dimer block only a single site. This deficiency can be somewhat compensated by making binding energy of the parallel dimer twice as large, $\varepsilon_{\parallel}\approx 2\varepsilon_{\perp}$.

Note that because the length of the dimer (equal to the same-chain site spacing $a^{(1)}=0.684$ nm) is much smaller than both the cross-chain site spacing $a^{(2)}=3.34$ nm and the vertical separation $\Delta z=3.4$ nm between turns of the helix, other orientations of the dimer are not possible.

The problem thus described is a complex one, but the similarities with the lattice gas models of condensed matter physics provide guidelines for how to proceed. These prescriptions, however, are aimed at the treatment of a periodic crystalline lattice, and, although the DNA sites exhibit helical symmetry, they do not constitute a periodic lattice in the strict sense of the term. However, an appropriate choice of coordinates can take advantage of the helical symmetry, so that, in these new coordinates, the positions of the sites will fall on a regular, one-dimensional lattice.

We will define these coordinates on a cylinder of radius $R_{DNA}$, as shown in Figs. 7 and 8. The first coordinate $x^{(1)}$ traces out a path with pitch angle $\psi$ along a single helical strand, and the other coordinate $x^{(2)}$ traces out a path with pitch angle $\alpha$ that connects corresponding sites on the two strands. Geometrically, a cylinder can be regarded as flat in the sense that it has no curvature. In Fig. 8, we show the way that the cylinder can be cut with scissors and unwrapped so that this lattice can be mapped on a flat surface as shown in Fig. 9. If we define the origin of coordinates $x^{(1)}$ and $x^{(2)}$ to be halfway along the helical path between the partners on the two chains, as shown in Figs. 8 and 9, then the positions of the sites on both strands form a one-dimensional lattice in the coordinates $(x^{(1)},x^{(2)})$.

These coordinates can be written simply in terms of cylindrical coordinates $(\phi,z)$ in matrix form as

[TABLE]

Inverting this gives definitions of the two coordinates as

[TABLE]

and

[TABLE]

With these definitions, the difference in coordinates between adjacent sites on the same strand is $\Delta x^{(1)}=a^{(1)}$, and the difference in coordinates between corresponding sites on the two strands is $\Delta x^{(2)}=a^{(2)}$. That is, $a^{(1)}$ is the distance along the helical path of a single chain from one phosphate ion to the next, and $a^{(2)}$ is the distance along a helical path from a phosphate ion on one chain to its partner phosphate ion on the other chain.

Next we define a lattice index $\ell$, which specifies the cell (altitude on the double-helix) and chain index $b$, which specifies the basis site, where $b=-\frac{1}{2}$ for the “down”($\downarrow$) chain and $b=\frac{1}{2}$ for the “up”($\uparrow$) chain, as shown in Fig. 10. Using these variables, the coordinates can be written as

[TABLE]

where

[TABLE]

Thus, although the sites on the DNA molecule do not constitute a periodic lattice in real space, they do constitute a lattice in an appropriately-defined coordinate space (see Fig. 10). As we will see, however, this choice of coordinates will make the form of the interaction potential more complicated as a result.

The use of $(x^{(1)},x^{(2)})$ instead of the cylindrical coordinates $(\phi,z)$ indicates a more fundamental shift in our description of the DNA double-helix. The two-dimensional surface on which the helices lie is a cylinder of radius $R_{DNA}$, and the helices inherit the cylinder’s geometric properties. The geometry of the cylinder, however, is locally indistinguishable from the geometry of the flat plane. One common consequence of this is that it is possible to smoothly wrap a flat sheet of paper around a cylinder. In contrast, it is not possible to smoothly wrap a sheet of paper around a sphere; this problem is well-known because of the geometrical distortions that occur in flat maps of a spherical Earth. Maps of a cylindrical surface, however, have no such distortions. This geometrical difference is quantified by the Riemannian curvature tensor, which vanishes for both the cylinder and the plane, but not for the sphereMisner et al. (1973). This means that the local geometry of the cylinder behaves in exactly the same way as the local geometry of the plane, so that a helix on a cylinder is geometrically equivalent to a line on a plane.

Our choice of coordinates is simply a map of the cylindrical surface that reduces the double-helix to two parallel lines, as shown in Fig. 10. The result of this map is that we have a linear lattice with two phosphate sites per cell, with the cells labeled from $\ell=-{\mathcal{N}}$ to ${\mathcal{N}}$. All $2{\mathcal{N}}+1$ cells are identical, and the lattice can be imagined to satisfy periodic boundary conditions.

The analogy describing DNA as a ladder wrapped around a cylinder of radius $R_{\text{DNA}}$ then has a true mathematical basis, because the local structure of the double-helix is equivalent to the structure of a “ladder” — a one-dimensional chain with a unit cell containing one site from each strand. If interactions are ignored, then the problem is described simply by this one-dimensional lattice.

IV DNA Interaction Potential

The geometrical structure of the double helix is important in determining the electrostatic energy $U_{\text{int}}$ of the sites (with and without dimers) interacting with one another. The full Hamiltonian must contain these contributions and can be written as

[TABLE]

where the interaction energy $U_{\text{int}}$ is the total interaction energy between all the charges on all the sites. If ${\mathcal{V}}_{(\ell_{1},b_{1}),(\ell_{2},b_{2})}$ is the screened Coulomb energy between a charge $q_{\gamma_{1}}$ at site $\ell_{1}$ on chain $b_{1}$ with another charge $q_{\gamma_{2}}$ at site $\ell_{2}$ on chain $b_{2}$, then $U_{\text{int}}$ can be written as

[TABLE]

where $n_{\ell,b}^{(\gamma)}$ is the number of particles of species $\gamma$ on the lattice site at $(\ell,b)$ and $q_{\gamma}$ is the charge of that species, in units of the magnitude $e$ of the charge of an electron. The factor of $1/2$ is to ensure that we are not double counting when we sum over all lattice sites, and we implicitly exclude same-site interactions $(\ell_{1},b_{1})\neq(\ell_{2},b_{2})$. The total charge on a given site is given by summing all the charges on that site, so that the total charge on the site $(\ell,b)$ is

[TABLE]

The interaction potential can then be written in terms the the total charge on each site as

[TABLE]

where ${\mathcal{V}}$ is the electrostatic interaction energy between the sites. This screened electrostatic energy between charges $q_{\gamma_{1}}$ and $q_{\gamma_{2}}$ at positions $(\ell_{1},b_{1})$ and $(\ell_{2},b_{2})$, (in SI units) is

[TABLE]

where $\epsilon$ is the electric permittivity of the medium between the two charges, and $d_{b_{1}b_{2}}(\ell_{1}-\ell_{2})$ is the straight-line distance between the two charges. Distances along the chain are invariant under translations by a lattice spacing, and so $d_{b_{1}b_{2}}(\ell_{1}-\ell_{2})$ depends only on the difference between the lattice site indices $\ell_{1}-\ell_{2}$. As in the noninteracting lattice-gas model, we have three species of “particles” on the lattice sites, vacancies of charge $q_{v}=-1$, parallel ($\parallel$) dimers of charge $q_{\parallel}=0$, and perpendicular ($\perp$) dimers of charge $q_{\perp}=+1$.

Under physiological conditions, the presence of monovalent salt ions such as Na+ leads to screening of the bare charges. Traditionally, screening is treated by modeling the ions as a continuous density, resulting in the nonlinear Poisson-Boltzmann equation Kittel (1996); Grosberg et al. (2002). The Poisson-Boltzmann equation is not analytically solvable, so it is often further approximated by linearization. The resulting Thomas-Fermi model is analytically solvable and gives the screened Coulomb (or Yukawa) potential between two charges given in Eq. (43) where $\epsilon=78.5\epsilon_{0}$ is the permittivity of water (with $\epsilon_{0}$ the permittivity of free space) and $q_{s}$ is the magnitude of the screening wave vector Ashcroft and Mermin (1976). Various names are ascribed to both the nonlinear and linearized equations, including Debye-Hückel, Thomas-Fermi, and Poisson-Boltzmann, but we will always refer to the Poisson-Boltzmann equation when we mean the nonlinear form and the (linearized) Thomas-Fermi equation when we mean the linearized form.

One must be cautious in using the screened Coulomb potential for screened interactions. The continuous density approximation from which the nonlinear Poisson-Boltzmann equation was derived constitutes a form of mean-field theoryGrosberg et al. (2002), which fails to describe any of the effects due to correlations between the ions such as charge inversion and condensation. Thus we cannot model screening of the DNA molecule by the dimers using the Thomas-Fermi model, or even the Poisson-Boltzmann equation. Instead, we must treat the dimers as individual particles and compute their interactions so that we do not ignore correlation effects.

For the screening due to the monovalent salt, however, the valence involved is small enough that a mean-field treatment may still be accurate Grosberg et al. (2002). Thus, we will treat the interaction energy $U_{\text{int}}$ in Eq. (42) using the screened Coulomb potential for the dimer-dimer interactions, with the screening vector $q_{s}$ given in the Thomas-Fermi model as Ashcroft and Mermin (1976)

[TABLE]

where $\beta=1/{k_{B}}T$ is the inverse temperature, $n$ is the concentration of ions, and $\ell_{B}=\beta e^{2}/(4\pi\epsilon)$ is the Bjerrum length Kunze and Netz (2002). The Bjerrum length is a characteristic length scale equal to the distance at which two proton charges interact with energy $k_{B}T$. It should be noted that these expressions differ slightly in the literature because of various authors’ choices of units for the electrostatics. Here, we use SI units.

By choosing a coordinate system based on the arc lengths around the surface of the cylinder, we have managed to align the sites into a regular lattice, but the potential $\mathcal{V}_{(\ell_{1},b_{1}),(\ell_{2},b_{2})}$ depends on the straight-line distance between the sites rather than the surface arc-length connecting them. We must then express Eq. (43) in terms of the straight-line distance function $d_{b_{1}b_{2}}(\ell_{1}-\ell_{2})$ between the sites at $(\ell_{1}a^{(1)},b_{1}a^{(2)})$ and $(\ell_{2}a^{(1)},b_{2}a^{(2)})$. The straight-line distance in cylindrical coordinates between two points on the surface of a cylinder of radius $R_{\text{DNA}}$ is

[TABLE]

and applying our change of variables (35,36) gives the straight-line distance between sites as

[TABLE]

where $\ell=\ell_{1}-\ell_{2}$ and $b=b_{1}-b_{2}$. Because of the translational invariance of the lattice, the distance in Eq. (LABEL:eq-gendist) depends only on the differences between the lattice and basis indices, and not their values individually. From this relation, we see the symmetry property of the distance,

[TABLE]

Since $b_{1}$ and $b_{2}$ are either $+\frac{1}{2}$ or $-\frac{1}{2}$, $b=b_{1}-b_{2}$ takes the three values $b=0$ for interactions of one site on a single chain with all other sites on the same chain, $b=1$ for interactions between a site on the “up” chain with with all the sites on the “down” chain, and $b=-1$ for interactions of one site on the “down” chain with all the sites on the “up” chain. We will write an up arrow $\uparrow$ will indicate $b_{1,2}=\frac{1}{2}$, and a down arrow $\downarrow$ will indicate $b_{1,2}=-\frac{1}{2}$. Plots of $d_{\uparrow\uparrow}(\ell)/R_{\text{DNA}}$ and $d_{\uparrow\downarrow}(\ell)/R_{\text{DNA}}$ are shown in Fig. 11a. There are wiggles in the curve because, as the chain wraps around the cylinder (see Fig. 3), the distances between lattice sites can become larger or smaller than they would be on a linear chain.

In order to understand the consequences of the screened Coulomb potential (43) depending on the straight-line distances $d$ between sites, we define a new notation for the potential, which is expressed in terms of the difference in lattice sites indices $\ell_{1}-\ell_{2}$,

[TABLE]

This potential decays exponentially as a function of the straight-line distance $d$ (inset of Fig. 11(b)). However, as a function of the lattice index $\ell$ as in Fig. 11(b), there are corresponding “wiggles” in the decay. This occurs because the distance between the sites decreases slightly as the helix completes a turn. These wiggles are present for every turn the helix makes, but the effect on the potential becomes negligibly small after about the first turn. Thus, by using the indices $\ell$ and $b$ to describe the relationships between sites, we have reduced the DNA double-helix structure to a regular one-dimensional lattice with two sites per cell at the cost of an irregular potential due to the geometrical structure of the DNA molecule. The potentials for other combinations of $b_{1}$ and $b_{2}$ can be related to those in Fig. 11(b) through the symmetry relations

[TABLE]

Fourier transforming in the site index $\ell$ block diagonalizes this matrix into $2\times 2$ blocks, corresponding to the chain indices $b_{1}$ and $b_{2}$. Explicitly, the Fourier transform is

[TABLE]

where the elements of ${\mathcal{F}}$ are independent of $b_{1}$ and $b_{2}$ and are given by

[TABLE]

and the dagger denotes the adjoint (complex conjugate transpose). We use periodic boundary conditions, so that the wave vector $k$ appearing here, which is dimensionless, takes the $2{\mathcal{N}}+1$ values

[TABLE]

where

[TABLE]

If ${\mathcal{N}}$ is large, the range of $k$ is essentially continuous from $-\pi$ to $\pi$.

The matrix elements of this Fourier transform are

[TABLE]

Substituting for ${\mathcal{V}}$ from Eq. (48), we have

[TABLE]

Since we regard the chain as having periodic boundary conditions, changing summation indices to $\ell=\ell_{1}-\ell_{2}$ and $\ell_{2}$ gives

[TABLE]

where the sum over $\ell_{2}$ becomes

[TABLE]

We can now write the Fourier transform of $V_{b_{1}b_{2}}(\ell)$ as

[TABLE]

Substituting this back into the matrix, we have

[TABLE]

Thus, the Fourier transform ${\tilde{\mathcal{V}}}$ of ${\mathcal{V}}$ is diagonal in the $k$ indices and contains $2\times 2$ blocks

[TABLE]

The symmetries (49) of the screened Coulomb potential on the DNA lattice are reflected in the Fourier transforms (58) as well:

[TABLE]

Also, because $V_{\uparrow\uparrow}(\ell)$ is an even function of $\ell$, $\tilde{V}_{\uparrow\uparrow}(k)$ is real. The Fourier transforms (58) are plotted in Fig. 12.

As can be seen in Fig. 12, the Fourier transforms $\tilde{V}(k)$ have regions of $k$ that are negative, which is a result of transforming with respect to our helical coordinate system. When the screened Coulomb potential (43) is Fourier transformed with respect to the straight-line distance $d$, the function is positive definite. We have instead Fourier transformed the potential a function of the lattice index $\ell$, resulting in the “wiggles” in Fig. 11(b) caused by the turns of the helix. These deviations from the exponential decay of $V(d)$ result in the deviations here from the positive-definite Fourier transform.

Using the symmetry operations in Eqs. (61a) and (61b) obeyed by the Fourier transforms allows us to simplify the $2\times 2$ blocks $\tilde{V}(k)$ in Eq. (60) to

[TABLE]

This matrix is easily diagonalized, yielding the real eigenvalues

[TABLE]

which are plotted in Fig. 13. The transformation that diagonalizes $\tilde{V}(k)$ is a unitary matrix $\xi(k)$ whose columns are the eigenvectors of $\tilde{V}(k)$. Choosing the $\lambda_{k+}$ eigenvector to be first, this transformation matrix is given by

[TABLE]

This diagonalizes the $2\times 2$ matrix $\tilde{V}(k)$ according to

[TABLE]

We can write the combined process of Fourier transforming ${\mathcal{V}}$ into block-diagonal form and then diagonalizing the $2\times 2$ blocks $\tilde{V}(k)$ as a single unitary transformation ${\mathcal{M}}$ given by

[TABLE]

which completely diagonalizes ${\mathcal{V}}$ such that

[TABLE]

Here $\Lambda$ is a diagonal matrix, with each element along the diagonal given by:

[TABLE]

where $\sigma$ is either $+$ or $-$.

Although we have diagonalized the potential with a unitary transformation, the fact that the transformation matrix is complex means that the charge density in the diagonal basis could also be complex, which is physically undesirable. However, the potential matrix in the position basis is real and symmetric, which means that it should be diagonalizable by a real orthogonal transformation matrix ${\mathcal{W}}$. In fact, the transformation matrix is not unique, since there is a degeneracy in the eigenvalues, since $\lambda_{k,\sigma}=\lambda_{-k,\sigma}$. This means that we can choose eigenvectors that are real by taking two orthogonal linear combinations of the eigenvectors of the unitary transformation $\cal{M}$ for $k$ and $-k$. Those two new eigenvectors, which represent the columns $(k,\sigma)$ for $k>0$ and $k<0$ of the new transformation matrix ${\mathcal{W}}$, are given by

[TABLE]

and

[TABLE]

Since $\xi_{b\sigma}(-k)=\xi_{b\sigma}^{*}(k)$, these can be written in terms of real and imaginary parts as

[TABLE]

and

[TABLE]

When $k=0$, the orthogonal transformation is the same as the unitary transformation, that is

[TABLE]

where

[TABLE]

Here the columns are labeled by the eigenvalue $\sigma=\pm$, and the rows are labeled by the chain $b=-1/2$ for the first row and $b=1/2$ for the second row.

This transformation can be written in a more compact form by introducing the unitary transformation matrix $\mathcal{Z}$, whose elements are defined by

[TABLE]

Then the orthogonal transformation is

[TABLE]

In deriving the thermodynamic quantities for the interacting lattice gas model, it will be convenient to be able to use this orthogonal transformation to diagonalize the potential.

V Partition Function for Interacting DNA Chains

For the interacting lattice gas, the partition function is similar to the noninteracting one in Eqs. (5) and (8). However, here the full Hamiltonian ${\mathcal{H}}$ from Eq. (39) is used, and then the partition function is given by

[TABLE]

The full Hamiltonian, including the interaction term from Eq. (42) written in matrix form, is

[TABLE]

In Sec. IV, we showed that the orthogonal transformation matrix ${\mathcal{W}}$ diagonalizes ${\mathcal{V}}$. In order to streamline the calculation, it will be useful to insert the identity matrix in the form ${\mathcal{W}}{\mathcal{W}}^{T}$, to the left and right of of ${\mathcal{V}}$ in the interaction term of the Hamiltonian, which gives

[TABLE]

where we have grouped factors to show the transformed quantities.

The Hamiltonian, with the interaction term written in the diagonal basis, can then be written as

[TABLE]

where $\tilde{\rho}={\mathcal{W}}^{T}\rho$ is the transformed charge density. We have not transformed the first term in the Hamiltonian to the diagonal basis because it is more convenient later in the calculation to have it in the position basis. Substituting this form into the partition function and writing the resulting expression as two separate exponentials, we have

[TABLE]

We can write the exponential of the sum in the interaction term as a product of exponentials as

[TABLE]

From this expression, we can see that the mean site occupancy for the interacting lattice-gas model can be obtained in the same way as in the noninteracting case (9),

[TABLE]

by differentiating the partition function with respect to the chemical potential.

In the noninteracting case, since all the information about the configuration was contained in $n_{\ell,b}^{(\gamma)}$, we were able to decouple the sum over configurations from the sum over lattice sites because the activity was independent of position, and $n_{\ell,b}^{(\gamma)}$ only appeared in its exponent. This is because only linear factors of $n_{\ell,b}^{(\gamma)}$ were in the exponent in the partition function. While it is still true that $n_{\ell,b}^{(\gamma)}$ contains the configuration dependence, it now appears quadratically in the exponent of the partition function, since $\tilde{\rho}_{k,\sigma}$ contains linear factors of the $n_{\ell,b}^{(\gamma)}$, and the exponent in the interaction term contains $\tilde{\rho}_{k,\sigma}^{2}$. This makes it impossible to decouple the sum over configurations from the sum over lattice sites. However, we can use an integral identity, known as the Hubbard-Stratonivich transformation, to replace the term quadratic in $\tilde{\rho}_{k,\sigma}$ in the exponential with a term linear in $\tilde{\rho}_{k,\sigma}$, meaning that there will only be linear terms in $n_{\ell,b}^{(\gamma)}$. This requires introducing the integral over an auxiliary field $\tilde{\Delta}_{k,\sigma}$. Then the sum over configurations is possible to do exactly, although at the cost of having to integrate over the auxiliary fields.

The integral identity, which is a complicated way of writing unity, as is shown in Appendix A, is given by

[TABLE]

where the integration path is over the real axis from $-\infty$ to $\infty$ when $\lambda<0$ and over the imaginary axis from $-i\infty$ to $i\infty$ when $\lambda>0$. In the previous work by Bishop and McMullenBishop and McMullen (2006), the problem was done in the position basis, and this integral identity, known as the Hubbard-Stratonovich transformation, was written with the matrix version of the potential in the exponent, as given in Negele and OrlandNegele and Orland (1988). This form creates a dilemma when there are negative and possibly zero eigenvalues of the potential, and it is difficult to determine the path of integration, since it changes depending on the sign of the eigenvalues. Zero eigenvalues would make the determinant of the matrix in that formula zero, and one then finds that there is a division by zero in the formula. By transforming to the diagonal basis, all these difficulties are avoided, since there is a separate integral for each eigenvalue, and if a particular eigenvalue is zero, the identity is not used at all.

To use the identity in Eq. (84), we make the identifications that $\lambda=\lambda_{k,\sigma}$, $\tilde{\rho}=\tilde{\rho}_{k,\sigma}$, and $\tilde{\Delta}=\tilde{\Delta}_{k,\sigma}$. Substituting this “1” into the partition function, we have

[TABLE]

Now we see that the two exponentials containing $\tilde{\rho}_{k,\sigma}^{2}$ cancel, as we planned. Of course, we have gained this convenience by introducing the auxiliary field ${\tilde{\Delta}}_{k,\sigma}$, and we will have to do the integral over this variable at a later stage. However, since ${\tilde{\Delta}_{k,\sigma}}$ does not depend on the configuration of the system, we can bring the sum over configurations and the leading exponential factors in the partition function inside the integral, and the expression for the partition function becomes

[TABLE]

It will be useful to define an action ${\mathscr{S}}$ such that the partition function can be written in the form

[TABLE]

where

[TABLE]

and

[TABLE]

Note that the factors of $\beta\lambda_{k,\sigma}$ are included here in the grand differential ${\mathcal{D}}[{\tilde{\Delta}}]$, rather than incorporating them in the action ${\mathscr{S}}[{\tilde{\Delta}}]$, as was done in Bishop and McMullenBishop and McMullen (2006) and in SievertSievert (2007). When these factors appear in the action, they introduce a term of the form $\ln\sqrt{\det(\beta v)}$. This adds a large constant value to the mean field results, which is subtracted out when the fluctuation terms are included. It actually has no real physical meaning and is part of the normalization of the integralAltland and Simons (2010).

We see that by diagonalizing first, our auxiliary fields are in the diagonal basis of the potential. If we transform the terms containing $\varepsilon_{\gamma}$ and $\mu_{\gamma}$ to this basis, they will no longer be diagonal. Therefore, we will leave those terms in the position basis. It will also be convenient to have the term containing $\sum_{k_{,}\sigma}\lambda_{k,\sigma}{\tilde{\rho}}_{k,\sigma}{\tilde{\Delta}}_{k,\sigma}$ in the position basis also. Therefore, the expression for the transformed charge density is

[TABLE]

in terms of the quantities in real space.

Substituting this into the action, we have

[TABLE]

It is convenient to identify the last term in this expression as a self energy of species $\gamma$,

[TABLE]

By transforming $\tilde{\Delta}_{k,\sigma}$ to the position basis as

[TABLE]

this self-energy can also be written in the position basis as

[TABLE]

This will be a useful form to use when we discuss the saddle point, or mean-field, approximation.

Keeping the first term in the position basis and the second term in the diagonal basis, we now write the action in a compact form, combining the self energy term with the first term in the action to obtain

[TABLE]

Analogous to the noninteracting case, we define the activity as

[TABLE]

With this definition, the partition function becomes

[TABLE]

The trace over configurations is now done over each site separately exactly as we did for the noninteracting case, where $n_{(\ell,b)}^{(\gamma)}=1$ for one and only one of $\gamma=\parallel,\perp$, or $v$, and zero otherwise. Thus, performing the trace gives

[TABLE]

and the grand partition function becomes

[TABLE]

where the effective action is

[TABLE]

This is the general result, which is in principle exact. To understand what the Hubbard-Stratonovich transformation has accomplished for us, recall that the interaction energy $U_{\text{int}}$ posed two difficulties: the additional configuration dependence and the coupling between the sites. By introducing auxiliary fields through the Hubbard-Stratonovich transformation (84), we managed to separate these two complications so that one factor $\exp\left[-\frac{\beta}{2}\sum_{k_{,}\sigma}\lambda_{k,\sigma}{\tilde{\Delta}}_{k,\sigma}^{2}\right]$ contains interactions between field fluctuations but no configuration dependence, and another factor $\exp\left[\beta\sum_{k,\sigma}\lambda_{k,\sigma}{\tilde{\rho}}_{k,\sigma}{\tilde{\Delta}}_{k,\sigma}\right]$ that is configuration-dependent (through ${\tilde{\rho}}_{k,\sigma}$) but decoupled. This allowed us to define a modified activity (96) that incorporated all the configuration dependence, making it possible to evaluate the sum over configurations directly, as in the noninteracting case.

What we have done in using the Hubbard-Stratonovich transformation is replace the interaction part of the partition function with its functional integral representation. From the form of Eq. (97), we see that the largest contribution to the partition function comes from the values of ${\tilde{\Delta}}_{k,\sigma}^{2}$ for which the effective action ${\mathscr{S}_{\text{eff}}}$ is stationary, i.e., at a saddle point. Thus a Taylor series expansion of the effective action (100) about the saddle point will yield an order-by-order approximation to the exact partition function (99).

VI Expansion in Powers of the Auxiliary Field

Since it is not possible to evaluate the integral in the partition function of Eq. (99) exactly, we will approximate it by expanding the action in Eq. (100) about a mean field $\Delta_{c}$, and include fluctuations to lowest order about this mean field. This mean field is determined by the saddle point of the integrand in Eq. (99), which is the largest contribution to the integral.

We can write the effective action then as an expansion about this uniform saddle point. As an aid to keeping track of the order of the expansions, we add a multiplicative factor $m$ in front of the effective action. as

[TABLE]

Here the expansion parameter $\delta{\tilde{\Delta}}$ is scaled by $m^{-1/2}$ by writing $\delta{\tilde{\Delta}_{k,\sigma}}\equiv\sqrt{m}[{\tilde{\Delta}}_{k,\sigma}-{\tilde{\Delta}}_{k,\sigma,c}]$. The factors of $m$ cancel in the quadratic term, so that it is of order $m^{0}$.

In evaluating the integral in the partition function, the largest contribution comes from the region near a saddle point, and the terms beyond the saddle point give corrections to the integral. Physically, the saddle point solution corresponds to a mean field approximation, and we begin with that. The saddle point is found by setting the first derivative in the expansion to zero, that is

[TABLE]

The second derivatives evaluated at the saddle point are proportional to the components of a matrix $\tilde{D}^{-1}$, which is the inverse propagator for the fluctuations of the auxiliary Hartree field $\Delta$ that we introduced to decouple the interparticle interactions. Because we want to absorb the normalization factor in the grand differential ${\mathcal{D}}[{\tilde{\Delta}}]$, we have defined each of the elements of $\tilde{D}^{-1}$ by the second derivative of ${\mathscr{S}}_{\text{eff}}[{\tilde{\Delta}}]$ divided by $-\beta\lambda_{k,\sigma}$ as

[TABLE]

We will show in Sec. VIII that, for a spatially uniform saddle point, this matrix is diagonal, that is, ${\mathscr{S}}_{\text{eff}}[{\tilde{\Delta}}]$ is diagonal in the $\delta{\tilde{\Delta}}_{k_{1},\sigma_{1}}$’s, that is,

[TABLE]

Although the Gaussian integral can still be done if $\tilde{D}^{-1}$ is not diagonal, the argument becomes much easier to follow if we assume at this point that $\tilde{D}^{-1}$ is diagonal because then the Gaussian integral can be done straightforwardly.

With this assumption that $\tilde{D}^{-1}$ is diagonal, the expansion of the effective action reduces to

[TABLE]

The grand canonical partition function can then be written as

[TABLE]

This can be written as a product of Gaussian integrals as

[TABLE]

and each Gaussian integral can be evaluated using the same change of variable and path of integration as was used in the Hubbard-Stratonovich transformation in Eq. (84) that originally was used to define the auxiliary fields, and as was discussed in Appendix A. Each Gaussian integral produces a result of the same form, regardless of whether $(\tilde{D}^{-1})_{(k,\sigma),(k,\sigma)}$ is positive or negative. This quantity can never be zero, since we would not have defined an auxiliary field for that case, which is the same as the case in which an eigenvalue of the potential is zero. Therefore, the integral can always be evaluated, and the partition function can be written as

[TABLE]

Note that all the factors of $\beta\lambda_{k,\sigma}$ have cancelled out due to our choices of the form of the effective action ${\mathscr{S}}_{\text{eff}}$ and the inverse Hartree field fluctuation propagator $\tilde{D}^{-1}$. Choosing these quantities carefully allows for the correct normalization of the integralAltland and Simons (2010). The product can now be brought inside the square root, and since $\tilde{D}^{-1}$ is diagonal, this gives the determinant of the matrix inside the square root, and so the partition function assumes the form

[TABLE]

and the partition function can be written as a single exponential,

[TABLE]

The grand canonical potential is then

[TABLE]

If we formally write the expansion in terms of orders in $m$, we have

[TABLE]

The leading term in the expansion is therefore given by the saddle point value

[TABLE]

and what is known as the “one-loop” correction term is

[TABLE]

VII The Saddle-Point or Mean Field Approximation

The first task in this section is to find an equation that determines the value of the auxiliary field at the saddle-point or mean-field value of the the effective action. In order to find this saddle point, we must set the first derivative of the effective action in Eq. (100) with respect to the auxiliary field ${\tilde{\Delta}}_{k,\sigma}$ to zero. That derivative is given by

[TABLE]

The activity $a_{\ell,b}^{(\gamma)}[{\tilde{\Delta}}]$ for species $\gamma$ is given by Eq. (96), and its derivative is

[TABLE]

where from Eq. (92), the derivative of the self-energy is

[TABLE]

Substituting this into the derivative of the action, we have

[TABLE]

Setting this first derivative to zero, we obtain an equation whose solution gives the auxiliary field in the saddle point approximation,

[TABLE]

We see that on the right-hand side of the equation is the transformation to the diagonal basis of the potential, as given in Eq. (93). We could therefore write the auxiliary field in the position basis as

[TABLE]

where $a_{\ell,b}^{(\gamma)}[\Delta]$ is the activity with the self energy in the diagonal basis $\tilde{\Sigma}_{\ell,b}^{(\gamma)}[{\tilde{\Delta}}]$, which by Eq. (94) has been replaced by the self energy in the position basis, $\Sigma_{\ell,b}^{(\gamma)}[\Delta]$. Either of these two equations, Eq. (119) or (120), may be considered the equation giving the saddle-point, or mean-field, value of the auxiliary field, and we will refer to this as the mean field equation. In this expression $a_{\ell,b}^{(\gamma)}[\Delta_{c}]$ is a function of all the $\Delta_{(\ell,b),c}$ elements, and so, for ${\cal{N}}=30$, which we have used in our calculations, there are ${\cal{N}}_{\text{sites}}=2(2{\cal{N}}+1)=122$ coupled nonlinear equations. Therefore, for simplicity, we make the assumption that the mean field is spatially uniform, which means that

[TABLE]

where $\Delta_{c}$ is constant. This is still a nonlinear equation in $\Delta_{c}$, and so this must be solved numerically with an iterative approach.

Before solving the mean field equation, it is useful to have a physical interpretation of the mean auxiliary field $\Delta_{c}$. The first step in this direction is the second task of this section, which is to find the mean site occupancies of species $\gamma$ at the mean field level. This is done as in the noninteracting case, by taking the derivative of the partition function with respect to $\mu_{\gamma}$. At this mean field level, the $\mu_{\gamma}$ dependence is contained in ${\mathscr{S}}_{\text{eff}}[{\tilde{\Delta}}_{c}]$ in Eq. (110), neglecting the term containing the inverse propagator $D^{-1}$. The expression for the mean occupancy can be gotten using Eq. (83) as

[TABLE]

Taking the derivative of the exponential, we can cancel out the partition function in the denominator and we are left with the derivative of the effective action as

[TABLE]

Substituting the effective action from Eq. (100), we have

[TABLE]

The derivative of the activity is

[TABLE]

and so the mean occupancy becomes

[TABLE]

The third task of this section is to relate the mean charge density per site to the mean occupancy. The charge density has the same form as in Eq. (22) for the noninteracting model, so that at the mean field level we obtain

[TABLE]

which, like the mean site occupancy, is independent of the site index because of our choice of spatially uniform mean field. This is equivalent to the equation for the saddle-point value of the auxiliary field because the right-hand side of the equation shows that

[TABLE]

on comparison with Eq. (120)

In the mean field or saddle point equation for the spatially uniform saddle point, the auxiliary field is the same for all lattice sites, and the self energy becomes spatially uniform also and can be written as

[TABLE]

where $\mathbb{S}_{\text{lattice}}$ is a lattice sum that is independent of the choice of lattice site $(\ell,b)$ for a long DNA double helix, and is given by

[TABLE]

For the parameters used in the figures, $\beta\mathbb{S}_{\text{lattice}}\approx 2.71$.

Consequently, at the mean field level, the self-energy simply represents the electrostatic energy of charge $q_{\gamma}$ interacting with the rest of the DNA lattice. This term plays the role of the self-energy $\Sigma$ in many-body quantum mechanics Bishop and McMullen (2006); Negele and Orland (1988); Altland and Simons (2010). Since the self-energy is independent of lattice site, so is the activity, which can be written as

[TABLE]

The simplifications made possible by the assumption of a spatially uniform saddle point allow us to rewrite the mean field equation given by Eq. (120) in an explicit form as

[TABLE]

Again, recall that $\varepsilon_{v}-\mu_{v}=0$, identically, $q_{\parallel}=0$, and $\mu_{\parallel}=\mu_{\perp}=\mu_{\text{dimer}}$. Therefore, when $\rho_{c}=0$, we can see from Eq. (132) that $q_{\perp}e^{-\beta(\varepsilon_{\perp}-\mu_{\text{dimer}})}+q_{v}=0$, and since $q_{\perp}=+1$, and $q_{v}=-1$, it follows that $\varepsilon_{\perp}=\mu_{\text{dimer}}$. Therefore, this is the point at which the charge density $\rho_{c}$ goes to zero. When $\varepsilon_{\parallel}=2\varepsilon_{\perp}$, the charge density goes to zero at $\varepsilon_{\parallel}=2\mu_{\text{dimer}}$, which is the same place it went to zero for the noninteracting model, as shown in Fig. 5.

Returning to the mean site occupancy in Eq. (126), the quantity inside the sum over $\ell$ and $b$ is constant and the sum yields ${\mathcal{N}}_{\text{sites}}$, and so the average occupancy per species becomes

[TABLE]

Using the mean-field value $\rho_{c}=\Delta_{c}$ of the charge density, calculated from Eq. (132), we have calculated the average occupation numbers $\langle n^{(\gamma)}\rangle$ via Eq. (133) using Eq. (131) for the activities. In Fig. 14, we have plotted $\delta\langle n^{(\gamma)}\rangle_{c}=\langle n^{(\gamma)}\rangle_{c}-\langle n^{(\gamma)}\rangle_{\text{NI}}$ , which are the differences in the mean-field occupancies from the occupancies in noninteracting case shown in Fig. 4. The corresponding change in the charge density $\delta\rho_{c}=\rho_{c}-\rho_{\text{NI}}$ is shown as the solid blue curve in Fig. 15.

The shifts in occupancy plotted in Fig. 14 can be primarily understood as a consequence of the electromagnetic “self energy” $\Sigma_{c}^{(\gamma)}=\rho_{c}q_{\gamma}\mathbb{S}_{\text{lattice}}$ of each species $\gamma$ interacting with the total charge $\rho_{c}$ of the lattice. Immediately, this implies that the parallel dimers are hardly modified at all, since they are electrically neutral and not directly affected by the total charge $\rho_{c}$ of the lattice. For the vacancies and perpendicular dimers, which are charged, the mean field corrections plotted in Figs. 14 and 15 act to reduce the overall charge of the lattice. Below $\beta\varepsilon_{\parallel}\approx 2$, the total charge of the lattice is positive, so the mean field corrections penalize the positively-charged perpendicular dimers and enhance the negatively-charged vacancies. Above $\beta\varepsilon_{\parallel}\approx 2$, the total charge of the lattice is negative, and this effect is reversed. Thus, the mean-field corrections tend to reduce (but not eliminate) the charge inversion at weak binding seen in Fig. 4. Note that the biggest effect of this electrostatic self energy is to heavily penalize the “naked lattice” of negatively-charged vacancies which was the favored ground state at large positive $\varepsilon_{\parallel}$ for the non-interacting Hamiltonian.

The fourth task of this section is to determine variance of the charge density from the mean field, similar to Eq. (24) for the noninteracting model, which gives a measure of the importance of spatial fluctuations in the charge density, can be written in the form

[TABLE]

where the average $\langle\rho^{2}\rangle_{c}$ of the square of the charge density, similar to Eq. (29) for the noninteracting model, is

[TABLE]

The variance $\sigma_{c}^{2}$ will be useful in knowing the size of the charge fluctuations in the mean field, and this quantity will appear in the next level of approximation. The difference in the standard deviation $\delta\sigma_{c}=\sigma_{c}-\sigma_{\text{NI}}$ is plotted as the red dashed curve in Fig. 15, where $\sigma_{\text{NI}}$ for the noninteracting model was plotted as the red dashed curve in Fig. 5.

Using a spatially uniform saddle point, which yields a charge density $\rho_{c}$ that is equal at the mean field level to the auxiliary field $\Delta_{c}$ in real space, it is interesting to write down the auxiliary field in the diagonal basis, given by

[TABLE]

From Eqs. (71) and (72), we see that when $k\neq 0$, ${\mathcal{W}}_{(k,\sigma),(\ell,b)}^{T}$ has the $\ell$ dependence of the form $e^{\pm ik\ell}$, which by Eq. (57) sums to zero. The only surviving term is for $k=0$ which when summed over $\ell$ gives $2{\mathcal{N}}+1$, and we have from Eqs. (73) and (136)

[TABLE]

There are two different results corresponding to the two different eigenvalue labels $\sigma=+$ and $\sigma=-$ in Eq. (63). Summing over $b$, these are

[TABLE]

where ${\mathcal{N}}_{\text{sites}}=2(2{\mathcal{N}}+1)$. The interpretation of Eq. (138) can be understood from Eq. (74). The sum over $b$ in Eq. (137) is a sum over the elements rows in Eq. (74), which for $\sigma=+$ sums the mean charge densities on the two chains and for $\sigma=-$ subtracts the charge densities on the two chains.

The fifth task of this section is to calculate the entropy for the spatially uniform saddle point, which is our mean field value of the entropy. This requires the effective action at this saddle point, which becomes

[TABLE]

The grand canonical potential is then

[TABLE]

We can now compute the entropy from the derivative of the grand thermodynamic potential with respect to the inverse temperature $\beta$, as in Sec. II. This derivative becomes much more complicated with the introduction of $U_{\text{int}}$ because the screened Coulomb potential (43) depends on $\beta$ through the screening vector $q_{s}$ (44), where

[TABLE]

This introduces a $\beta$-dependence into both terms of Eq. (139), so that the mean field entropy is given by the more complex expression

[TABLE]

where

[TABLE]

is the mixing entropy of the three species using the mean field occupation numbers and

[TABLE]

arises directly from the self energy $\Sigma\propto\mathbb{S}_{\text{lattice}}$. The occupation numbers $\langle n^{(\gamma)}\rangle_{c}$ used in this expression are those of the mean-field level, given by Eq. (133) and shown in Fig. 14. The mean-field entropy $\bar{S}^{(0)}$, given in Eq. (142), is plotted in Fig. 16, together with $-\langle n^{(\gamma)}\rangle_{c}\ln\langle n^{(\gamma)}\rangle_{c}$ for each species. Also shown are the two constituents $\bar{S}_{c}^{(A)}$ and $\bar{S}_{c}^{(B)}$ of the total entropy. This plot shows that, as with $S_{NI}$ in Fig. 6, the mean-field entropy is greatest for $\varepsilon_{\parallel}$ near $\mu_{\text{dimer}}$ and decreases outside this region. The difference $\delta\bar{S}_{c}=\frac{S_{c}-S_{\text{NI}}}{k_{B}T}$ between this mean-field entropy for the interacting model and the entropy for the noninteracting model is plotted in Fig. 17. Note that this curve has roughly the same shape as that of $\delta\sigma_{c}$ in Fig. 15, with a large increase in disorder at large positive $\varepsilon_{\parallel}$ and small changes elsewhere.

By employing field-theoretic techniques borrowed from quantum mechanics, we have included the interaction energy $U_{\text{int}}$ in the partition function and found a mean-field approximation to the exact partition function. However, these mean-field level calculations assume the same average charge $\rho_{c}$ for all DNA molecules in the same solution. Thus, mean-field calculations would predict a repulsive interaction between two DNA molecules, which would not lead to DNA condensation. Indeed, this repulsive interaction is all that this zeroth-order calculation can predict. It is the fluctuations of the charge from its mean-field value that enable the net attraction observed experimentally; to understand the attractive forces between DNA molecules in solution, it is necessary to extend this treatment by at least one additional order. The next order in the expansion (101) will yield Gaussian-type integrals, which are analytically solvable for the lowest-order fluctuations in the average charge and other parameters. Furthermore, once lowest-order fluctuations have been included, it will be possible to calculate thermodynamic properties of the system from the partition function to meaningful levels, such as the entropy changes due to small fluctuations in the local charge.

VIII Inverse-Hartree-field-fluctuation propagator

In order to calculate the inverse-Hartree-field-fluctuation propagator $\tilde{D}^{-1}$, we need to find the second derivatives of the effective action, evaluated at the saddle point. Since $\tilde{D}^{-1}$ is diagonal, from Eq. (103), we can write

[TABLE]

Since we already calculated the first derivative in Sec. VII, we need one more derivative, now with respect to ${\tilde{\Delta}}_{k_{2},\sigma_{2}}$. Substituting from Eq. (LABEL:eq-FirstDerivSeff)

[TABLE]

Performing the derivative, we have

[TABLE]

Substituting the derivative of the activity from Eqs. (116) and (117), we have

[TABLE]

The quantity in square brackets is just the constant $\sigma_{c}^{2}$, the variance of the charge density, given in Eq. (134). Since ${\mathcal{W}}$ from Eqs. (71), (72), and (73), is an orthogonal matrix,

[TABLE]

and this means that the matrix is diagonal. With these substitutions, the inverse-Hartree-field-fluctuation propagator becomes

[TABLE]

This is the inverse propagator for the fluctuations of the auxiliary Hartree field that we introduced to decouple the interparticle interactions, expressed in the diagonal basis.

The inverse propagator $(\tilde{D}^{-1})_{(k_{1},\sigma_{1}),(k_{2},\sigma_{2})}$ in the diagonal basis is plotted in Fig. 18 for the two eigenvalues $\lambda_{k\pm}$, shown in Fig. 13. Notice that the shapes are those of the eigenvalues, but shifted upward. They both have minima at $k=\pi$, with the one corresponding to $\lambda_{k+}$ having its maximum at $k=0$ and the one corresponding to $\lambda_{k-}$ having maxima near $k=0$. The corresponding plots of the propagator $\tilde{D}_{(k_{1},\sigma_{1}),(k_{2},\sigma_{2})}$ are shown in Fig. 19, also in the diagonal basis. The propagator increases with $\beta\varepsilon_{\parallel}$ near $k=\pm\pi$, and this can be attributed, at least in part, to the increase in $\sigma_{c}$, shown in Fig. 15, because of the dependence shown in Eq. (148). Although the auxiliary field propagator normally describes the behavior of fluctuations of a physical field like the charge, its definition given by Eq. (LABEL:eq-partitonfctproduct) shows that the field here have been rescaled by the square root of the eigenvalue. This was discussed earlier in the context of entropy in the paragraph after Eq. (89), and so this form is convenient for the entropy calculation, our focus here, but its interpretation is not quite so direct.

IX Inclusion of Lowest-Order Fluctuations

Going beyond mean field to one-loop order requires use of the full partition function given by Eqs. (108) and (110). This has consequences for the site occupancies, the average charge per site, the charge variance, and the entropy. The expression for site occupancy has the same form in terms of the partition function as for the noninteracting and mean field cases in Eqs. (9) and (83), and can be written explicitly for the one-loop order of approximation as

[TABLE]

which is obtained from Eq. (113) with $m$ set equal to one. In that expression, $m$ was the parameter that was introduced in Sec. VI as an aid to identifying the various orders in the expansion. This can be written as the sum of two terms,

[TABLE]

where $\langle n^{(\gamma)}\rangle_{c}$ is the occupancy of species $\gamma$ at the mean field level, given by Eq. (133) ,

[TABLE]

and $\delta\langle n^{(\gamma)}\rangle$ is the correction to that mean site occupancy due to fluctuations, which is given by

[TABLE]

with the inverse propagator $\tilde{D}_{(k_{1},\sigma_{1}),(k_{2},\sigma_{2})}^{-1}$ from Eq. (150) written in matrix form as

[TABLE]

where $I$ is the identity matrix and $\Lambda$ is the (diagonal) matrix of eigenvalues of the potential. Since $\tilde{D}^{-1}$ is diagonal, its determinant is given by the product of its diagonal elements as

[TABLE]

The logarithm of this determinant then gives the sum over logarithms of diagonal elements as

[TABLE]

The only $\mu_{\gamma}$ dependence is contained in $\sigma_{c}^{2}$, and its derivative is given by

[TABLE]

so that the derivative of $\ln(\det D^{-1})$ is given by the coefficient

[TABLE]

Then the mean fluctuation in occupancy $\delta\langle n^{(\gamma)}\rangle$ becomes

[TABLE]

These changes in the mean field occupancies at this one loop order are plotted in Fig. 20. These changes are not large, and the fluctuation corrections amount to no more than about 10% of the total. The most dramatic effects are seen in the region of large positive $\varepsilon_{\parallel}$, with the corrections to the vacancies and perpendicular dimers going in the opposite direction to the corrections due to the mean field in Fig. 14. This can be interpreted as a relaxation of the rigid penalties imposed by the electrostatic self energies at the mean field level. An enhancement of parallel dimers is also seen across the whole range of $\varepsilon_{\parallel}$ that peaks around $\varepsilon_{\parallel}=\mu_{\text{dimer}}$.

The mean charge density at this one-loop-correction level of approximation is given by the same form as Eqs. (22) and (127) for the noninteracting and mean field approximations and is given by

[TABLE]

where now $\langle n^{(\gamma)}\rangle_{\text{tot}}$ is taken from Eqs. (152) and (160). and the charge density becomes

[TABLE]

where

[TABLE]

The change in the charge density $\delta\rho_{\text{tot}}$ is shown as the solid blue curve in Fig. 21. The most notable features are a negative contribution to the total charge above $\beta\varepsilon_{\parallel}\approx 2$, reflecting the addition of more negatively charged vacancies due to the fluctuations, and a positive contribution to the total charge around $\varepsilon_{\parallel}\approx\mu_{\text{dimer}}$. Interestingly, this positive contribution is not associated with the increase of parallel dimer occupancy, since parallel dimers are electrically neutral. Rather, it is due to the fact that, although both the negatively charged vacancies and positively charged perpendicular dimers are both reduced around $\varepsilon_{\parallel}\approx\mu_{\text{dimer}}$, the vacancies are decreased more.

The variance in the charge is given by

[TABLE]

The average of the square of the charge density is given by

[TABLE]

and is written explicitly as

[TABLE]

where $\langle\rho^{2}\rangle_{c}$ is given by Eq. (135), and $\delta\langle\rho^{2}\rangle_{\text{tot}}$ works out to be

[TABLE]

The charge variance can then be written as

[TABLE]

where the fluctuation correction is

[TABLE]

The standard deviation can then be written as

[TABLE]

The change $\delta\sigma_{\text{tot}}$ in the standard deviation of the charge densities due to fluctuations are given in Fig. 21, where the fluctuations account for changes typically of the order of 10% of the mean field value.

Having obtained results for the effects of fluctuations on the mean occupancies and charge densities, we now consider the corresponding corrections to the entropy. For that, we need the grand thermodynamic potential, which was expanded in Sec. VI, as Eqs. (112)-(114), with the mean field and lowest order fluctuation terms given by

[TABLE]

The corresponding expression for the interacting entropy is

[TABLE]

where the mean field entropy $S_{c}$ was derived and the results presented in Sec. VII in Eqs. (142)-(144). The fluctuation contribution to the entropy $\delta S_{\text{tot}}$ is given by

[TABLE]

which requires performing another derivative.

To evaluate this contribution to the entropy $\delta S_{\text{tot}}$ in Eq. (173), we return to Eq. (157) and substitute into the fluctuation contribution to the entropy from Eq. (173), which after performing derivatives gives

[TABLE]

where

[TABLE]

and

[TABLE]

The derivative of the eigenvalue $\lambda_{k,\sigma}$ of the potential is given by

[TABLE]

where $\tilde{V}_{\uparrow\uparrow}(k)$ is the Fourier transform of the interaction potential of a chain with itself and $\tilde{V}_{\uparrow\downarrow}(k)$ is the Fourier transform of the interaction potential between the two chains. Both of these potentials depend on the screening vector $q_{s}$, and the only $\beta$ dependence of the potential is contained in $q_{s}$, which is proportional to $\sqrt{\beta}$. The derivative we need is then given by Eq. (141). As a consequence, the derivative of the Fourier transform of the potential is

[TABLE]

where the derivative of the potential in the position basis is given by

[TABLE]

Substituting and using the symmetry relations in Eqs. (47) and (49), the derivative of the Fourier transform of the potential is

[TABLE]

Using these relations, the derivative of the eigenvalues becomes

[TABLE]

where ${\text{Re}}[\tilde{V}_{\uparrow\downarrow}(k)]$ and ${\text{Im}}[\tilde{V}_{\uparrow\downarrow}(k)]$ are the real and imaginary parts of $\tilde{V}_{\uparrow\downarrow}(k)]$.

Eq. (174) also requires the derivative of the variance $\sigma_{c}^{2}$, where $\sigma_{c}^{2}$ is given by Eq. (134). This contains the average charge density from Eq. (127) and the average of the square of the charge density from Eq. (135). This derivative is

[TABLE]

where

[TABLE]

[TABLE]

[TABLE]

and

[TABLE]

Here $[\cdots]_{\text{av}}$ is the weighted average of the specified quantity over the distribution given by the mean-field occupancies $\langle n^{(\gamma)}\rangle_{c}$ of the three species.

It is interesting to note that the combinations $[\rho^{2}{\cal{E}}]_{\text{av}}-\langle\rho^{2}\rangle_{c}{\cal{E}}_{c}$ and $[\rho{\cal{E}}]_{\text{av}}-\rho_{c}{\cal{E}}_{c}$ can be written in form similar to that of the zeroth-order entropy. These forms are

[TABLE]

and

[TABLE]

The fluctuation entropy contribution $S_{f}^{(A)}$ is easy to calculate, given the eigenvalues $\lambda_{k,\sigma}$. Eqs. (177)-(LABEL:eq-rhoepsavg) simply collect the results that allow us to calculate $S_{f}^{(B)}$. The results for the total fluctuation entropy $\delta S_{\text{tot}}$ and the two contributions to it are shown in Fig. 22 for $\varepsilon_{\perp}=\frac{1}{2}\varepsilon_{\parallel}$. Like the mean field entropy, the fluctuation entropy becomes larger in magnitude as the binding becomes weak. However, the fluctuation contribution is actually negative in this region, suggesting that the fluctuations are leading to increasing order. These corrections to the entropy are typically an order of magnitude larger than the difference between the mean field and noninteracting values of the entropy.

X Discussion

Entropy forms an important contribution to the free energies of many biological systems. Here we summarize the results for the entropies and particle densities for a lattice gas model representing the adsorption of molecules on double-stranded DNA. The entropies and particle densities are presented at three successive levels of approximation, first, for the noninteracting system then at the mean-field level, and finally with the inclusion of fluctuations to one-loop order. Those fluctuation corrections result from correlations induced by the electrostatic forces among the dimer molecules themselves and with the electrically charged DNA substrate.

Perhaps the most transparent descriptors of the DNA system’s behavior are the average occupation numbers $\langle n^{(\parallel)}\rangle$, $\langle n^{(\perp)}\rangle$, and $\langle n^{(v)}\rangle$ representing the thermal average of the numbers of parallel-adsorbed dimers, perpendicular-adsorbed dimers, and vacancies, respectively, on each site. These particle densities are proportional to the probability that a given site would be found to be occupied by a particular species. These site-occupation numbers are shown in Fig. 23 for the three levels of approximation. The short-dashed curves are for the noninteracting level of approximation $\langle n^{(\gamma)}\rangle_{\text{NI}}$, the long-dashed curves for the mean field level $\langle n^{(\gamma)}\rangle_{c}$, and the solid curves for the inclusion of fluctuations $\langle n^{(\gamma)}\rangle_{\text{tot}}$. The corresponding charge densities and standard deviations of the charge density are shown in Fig. 24.

The site occupancies exhibit similar qualitative behavior in all three levels of approximation. Due to the energetics of binding $\varepsilon_{\parallel}=2\varepsilon_{\perp}$, parallel dimers occupy the greatest fraction of sites for all negative binding energies $(\varepsilon_{\parallel},\varepsilon_{\perp}<0)$, followed the perpendicular dimers and then the vacancies, at least for the parameters used in these calculations. At weak binding ($\varepsilon_{\parallel},\varepsilon_{\perp}\rightarrow 0$) in the noninteracting system, the difference between these energies becomes negligible, and both $n^{(\parallel)}$ and $n^{(\perp)}$ approach the same value of $40\%$. When the binding energy $\varepsilon_{\parallel}$ or $\varepsilon_{\perp}$ increases and reaches the chemical potential $\mu_{\text{dimer}}$ of the dimers in solution, they are repelled from the surface of the DNA and go into solution, leaving a “naked” lattice of negatively-charged vacancies, with the parallel dimers being more repelled than perpendicular ones.

In contrast to the case of noninteracting sites, the inclusion of interactions between sites gives a pronounced gap of about $25\%$ by which parallel adsorption exceeds perpendicular adsorption in the weak binding limit. This enhancement of parallel adsorption in the mean-field approximation results from the inclusion of the charge on the perpendicular dimers in an average or mean field sense, which penalizes the addition of positively-charged perpendicular dimers when the total charge on the DNA lattice is positive. Fluctuations reduce this effect slightly. Recall that the perpendicular dimers are positively charged because one end is not attached to a negative binding site, and this means that it is energetically unfavorable to have a second one nearby. Parallel binding, on the other hand, neutralizes the negative charge on the sites, so that a DNA double strand completely covered in parallel-adsorbed dimers would be electrically neutral. The decrease in perpendicular-dimer occupancy caused by this electrostatic repulsion of like charges leads to an increase in sites occupied by vacancies. Vacancies are negatively charged and are energetically favored when the perpendicular dimers provide a positively-charged region. There may, in fact, be correlations resulting from the lower electrostatic energy of a configuration in which sites alternate between perpendicular dimers and vacancies, i.e., between positively and negatively charged sites, in a Wigner-lattice-like state similar to those described in the literature Nguyen and Shklovskii (2002a); Grosberg et al. (2002); Gelbart et al. (2000).

The average charge $\rho_{c}$ per site depends only on the perpendicular dimers and vacancies since the parallel dimers neutralize sites, so that the charge density is written simply as

[TABLE]

and this connection can be seen in the plots of charge density $\rho$ in Fig. 24, when compared with Fig. 23. For all approximations, the average charge is positive at negative and slightly positive binding energies, indicating charge inversion. This is a result of the occupation numbers shown in Fig. 23, because parallel binding, which neutralizes the charge, increasingly dominates in the strong binding limit ($\varepsilon_{\parallel}=2\varepsilon_{\perp}\ll 0$). Similarly, the decrease in the magnitude of the charge inversion from the noninteracting to the mean-field approximations can be explained by the enhancement of parallel adsorption and reduction of perpendicular adsorption due to like-charge repulsion. The fluctuation corrections to the site occupancies are not large when compared with the mean-field values. They only seem significant at positive binding energies where the site occupancies are shifted in the direction of the noninteracting system. Of course, the site occupancies are single-particle properties and do not measure spatial correlations.

The behavior of the entropy is shown in Fig. 25, in which the entropy per site is plotted versus $\beta\varepsilon_{\parallel}$ for all three levels of approximation. The three curves are practically indistinguishable until the binding energy becomes quite weak, indicating little effect of correlations induced by electrostatic repulsion on the entropy in that region. More pronounced differences are seen in the region of $\varepsilon_{\parallel}\approx\mu_{\text{dimer}}$ and more positive binding energies ($\varepsilon_{\parallel}=2\varepsilon_{\perp}\gg 0$), where the system becomes more disordered, and where the dimers tend to leave the surface and enter the solution. The fluctuation corrections to the entropy are pronounced in that region, and tend to reduce the entropy, implying that the electrostatic correlations tend to increase the order in the lattice gas model. Surprisingly, this decrease in entropy at large positive $\varepsilon_{\parallel}$ due to fluctuations is so significant that it overcompensates for the entropy generated at the mean field level, leading to a total entropy $S_{\text{tot}}$ that is even lower than in the noninteracting case $S_{\text{NI}}$.

The form of the noninteracting entropy comes from the lattice gas model, in which parallel adsorption is treated as “local,” only occupying one site. When the geometrical blocking effect is included, in which parallel-adsorbed dimers occupy two sites, we would expect the entropy to differ significantly from the $n\ln n$ result even at the noninteracting level. It is unclear how large a contribution to the full entropy the blocking effect provides, but it is certainly lower for dimers than for longer polymers, for which the nonlocality is greater. Because of this, it will be important in future work to develop a way of incorporating the blocking effect into these types of field-theoretic models.

Dimers are certainly the simplest example of polymers, although biological polymers are typically considerably longer. There is extensive work done on the dimer model, beginning with the work of Fisher Fisher (1961); Fisher and Stephenson (1963). Some of this work has been motivated by the analogy between dimer models and quantum spin systems Fisher (1966). In this context, Fisher has derived expressions for the contributions of hard-core crowding effects to the partition function, including the geometrical blocking effect. Fisher’s results include an expression of the noninteracting partition function for dimers on a 1D linear lattice. In a future publication, we plan to adapt Fisher’s approach to study the statistical consequences these blocking effectsBaker III (2017).

In addition, it has been suggested that methods can be adapted from analogous problems in particle physicsAdams and Chandrasekharan (2003) and condensed matter physicsSyljuasen and Sandvik (2002). This involves using numerical codes like those of Adams and ChandrasekharanAdams and Chandrasekharan (2003) that are used for simulations in lattice quantum chromodynamics.

Modeling the behavior of DNA and other biologically active polymeric molecules under physiological conditions is an area in which quantitative calculation is particularly difficult. This is because it is necessary to accommodate simultaneously the influences of geometry, electrostatic forces, and charge correlations. All of these play roles of varying significance in determining the properties of these molecules. Despite these challenges, considerable progress has been made in extending our physical understanding of these systems, and some surprising new physics, such as charge inversion and condensation, have been discovered in the process.

The work presented here focuses on electrostatic effects and especially on the role they play in the entropy, with the overall goal of determining the entropy as a contribution to the free energy and as a driving force in biochemical reactions. For the DNA-based model that we consider, we find that, particularly in the low-coverage regime, a uniform mean field theory gives changes $\sim 20\%$ from the usual noninteracting entropy term. We also find, though, that the fluctuation contributions are both considerably larger and in the opposite direction. This establishes that the fluctuations are significant, and the opposite sign of their contribution suggests that they lead to additional order in the system.

Appendix A A Gaussian Integral Identity

In this appendix, we prove the integral identity that we use in the derivation of the partition function $Z_{G}$ for the interacting model. The identity is

[TABLE]

The integration path is over the real axis from $-\infty$ to $\infty$ when $\lambda<0$ and over the imaginary axis from $-i\infty$ to $i\infty$ when $\lambda>0$. This integral can be evaluated by noticing that the exponent is a factor times a perfect square,

[TABLE]

Substituting, the integral becomes

[TABLE]

Now it is clear why the path of integration must be chosen differently for positive versus negative $\lambda$, since it must be chosen such that the integral converges. If $\lambda=0$, the identity fails, but in that case it is clearly unnecessary. If $\lambda<0$, we can write $\lambda=-|\lambda|$, and the integral converges for $\tilde{\Delta}$ ranging over the real axis from $-\infty$ to $\infty$. Substituting $y=\tilde{\Delta}-\tilde{\rho}$, with $y$ ranging from $-\infty$ to $\infty$ the integral can be evaluated as

[TABLE]

For $\lambda>0$, the situation is slightly more complicated. In that case, in order to have a decaying Gaussian in the integrand, we must integrate $\tilde{\Delta}$ over the imaginary axis from $-i\infty$ to $i\infty$. This is accomplished by substituting $y=-i\left(\tilde{\Delta}-\tilde{\rho}\right)$, for $y$ ranging from $-\infty$ to $\infty$. In this case, we have

[TABLE]

where the $i$ the numerator cancels the $\sqrt{-1}$ in the denominator. This proves the identity for positive $\lambda$. Therefore, we have seen that this identity works for positive and negative values of $\lambda$, and, although it does not work $\lambda=0$, it is not needed in that case.

Appendix B Chemical Potential of Free Dimers

The chemical potential of a solution of free dimers is the derivative of the Helmholtz free energy of the dimers $F_{\text{dimer}}$ with respect to the number of dimers $N_{\text{dimer}}$, with temperature and volume held constant,

[TABLE]

Assuming that $F_{\text{dimer}}$ is due only to electrostatic interactions, in a dilute solution the free energy is just the number of dimers times the electrostatic self energy $E_{\text{dimer}}$ of a single dimer. Then the chemical potential becomes

[TABLE]

and so we need the electrostatic energy of a single dimer. Suppose we model a dimer as a cylinder of cross-section radius $R$ and length $L$ of total charge $Q$, as shown in Fig. 26. An exact solution to this problem has been provided by CiftaCiftja (2012), and the derivation proceeds as follows. The volume charge density of the cylinder is given by

[TABLE]

and the Coulomb self-energy is given by

[TABLE]

where $\epsilon$ is the dielectric permeability of the medium surrounding the cylinder. In evaluating the integrals, the length to radius ratio is defined as

[TABLE]

and the exact result, as calculated by CiftaCiftja (2012), is given by

[TABLE]

where the function ${{}_{p}F_{q}}(a_{1},a_{2},\ldots,a_{p};b_{1},b_{2},\ldots,b_{p};z)$ is the generalized hypergeometric function whose series expansion is

[TABLE]

where the variable $z$ can be either real or complex and $(c)_{n}$ is the Pochhammer symbol defined as

[TABLE]

Here $\Gamma(z)$ is the gamma function and $(c)_{0}=1$. For a long, thin rod, the logarithm term would dominate, but a dimer approximated by a cylinder, pictured in Fig. 26, does not fall into this regime.

As shown in the picture, for the dimer, the radius of the dimer is $R=2a$, where $a$ is the lattice spacing along a DNA helical chain. The length of the dimer is $L=2a$, and so the length ratio is $\xi=\frac{L}{R}=4$. The charge is $Q=2e$, where $e$ is the magnitude of the charge of an electron, and the chemical potential becomes

[TABLE]

For the parameters given in the text, the chemical potential for the dimer, in units of $\beta=\frac{1}{k_{B}T}$ at $T=300K$ is given by $\beta\mu_{\text{dimer}}=0.79$. It is interesting to note that approximately the same value was obtained by Nguyen and ShklovskiiNguyen and Shklovskii (2002b) and Bishop and McMullenBishop and McMullen (2006) assuming that the dimer is an infinitely thin rod of length $L$ but that there were upper and lower screening cutoffs. This exact solution avoids the use of those cutoffs.

Acknowledgement

We are grateful to Dr. Kevin R. Ward, Director of Research, Department of Emergency Medicine, University of Michigan, for introducing us to the role of charge in biological systems and to the practical consequences that could be achieved by controlling the charge distribution on biological molecules.