Numerical Solution for the Counterions Distribution in a Hexagonal DNA Lattice within Mean Field Theory Using Finite Element Method

Abstract

DNA has proven to be a promising material in fabrication and construction of complex structures with precise controlled nanoscale features. Most of the systems involving DNA are functioned in an aqueous solution where DNA molecules are strongly negatively charged. Therefore, understanding electrostatics of DNA system is essential for better understanding and design of DNA as a biomaterial and as a biological structure. In this work, the mean-field Poisson-Boltzmann equation for the distribution of mobile ions in a two-dimensional hexagonal lattice of DNA cylinders is solved numerically using finite element method. The weak formulation of the Poisson-Boltzmann equation is derived. The equation is then solved numerically using FreeFem++ scripting language. Our results show that the excess counterions of DNA dominate over the bulk ion concentration for the physiological salt concentration considered, and they condense on the DNA surface leading to very high charge density. The results also demonstrate the strong influence of the entropic confinement of the ions when the distance between neighboring DNA is smaller than 10 nm. This effect cannot be ignored in this case, and should be taken into account in any electrostatic investigation of DNA system.

Top view of a hexagonally DNA lattice (left) and the electrostatic potential iso-surfaces around a DNA molecule approximated as a negatively charged cylinder (right).

1. Introduction

Deoxyribonucleic acid (DNA) has been proved as great potential material in fabrication and construction complex structures with precise control over nano-scale features.^1–3) DNA is an important molecule that carries genetic information encoded for the growth and development of all life forms including viruses.⁴⁾ Due to the unique structural, and physicochemical properties, DNA has many interesting and promising applications in nanomedicine and biomedical devices such as drug delivering, biosensors, synthetic nanopore formation for diagnostics, and gene delivery systems.^1,5,6)

In soft materials, especially those of biomaterials, electrostatic interactions associated with charged objects play an important role.^4,5,7) In aqueous environment, phosphate groups of the DNA backbone lose their protons and DNA becomes negatively charged. Hence, in a water solution DNA is highly charged with the effective density of one elementary (negative) charge every 1.7 Å along its principle axis. Similarly, if we consider DNA as a charged cylinder with a radius of 1 nm, then the surface density is about 1e/nm². Both of these linear or surface charge densities are among the highest values for biological molecules. Therefore, electrostatics dictates many aspects of the structure and functions of DNA’s system, and investigation of electrostatics of DNA containing systems remained one of most frequent topics in recent scientific literature despite a large of bodies of work already published. Understanding of these electrostatics will be the main focus of this work as well. In many biotechnological and biological system, such as DNA packaging in cells, inside viruses, in toroids, DNA molecules arrange in hexagonal bundle^8,9) (see Fig. 1). It would be interesting to understand electrostatics of these hexagonal DNA bundles. This is the main purpose of this work.

Fig. 1

DNA bundle is modeled as a hexagonal lattice with d is the distance between DNA cylinders. Each DNA molecule is modeled as a hard-core cylinder with radius R. The periodic space-filling unit cell (Wigner-Seitz cell) is the hexagonal region. Due to symmetry, only the green domain (1/12 of the cell) is used for numerical finite element method solution. The rest of the cell can be obtained by symmetry. The parallelogram is the basic unit cell of the hexagonal lattice.

Many quantitative investigations of DNA system in particular or charged macroions in general are simulation based, especially when describing the overcharging effect or like charge attraction phenomena. Examples are Monte Carlo simulation of charge systems in planar geometry,^10,11) cylindrical geometry^12–14) and spherical geometry.^15,16) Phenomenological theory based on mean-field theories are also developed.^17–20) Nevertheless, the mean-field Poisson-Boltzmann equation remains the corner stone of many analyzation tools of molecular software packages such as VMD,²¹⁾ MM/PBSA²²⁾ which provides intuitive understanding of electrostatics of these systems at molecular level. Therefore, we will employ this mean-field equation for studying hexagonal DNA lattice.

The Poisson-Boltzmann equation is a very useful analytical approximation with many applications. Because the equation is non-linear partial differential equation, it has closed-form analytical solutions only for a limited number of simple charged boundary conditions.²³⁾ On the other hand, by solving it numerically or within some further approximations or limits, we can obtain the ionic profiles as well as the free energy of complex structures. Therefore, in this work we focus on studying the electrostatics by numerically solving Poisson-Boltzmann equation. The method of choice will be the finite element method which can be advantageous to standard finite difference method in many situations.

The paper is organized as follows. In Section 2, the weak formulation of the Poisson-Boltzmann equation for hexagonal DNA lattice is derived. In Section 3, numerical parameters and system setup are presented. Various results on the concentration of the ions, the ion concentration variation in distance from DNA axis, hexagonal DNA lattice constants are presented and discussed. We conclude in Section 4.

2. Derivation of Poisson-Boltzmann for Hexagonal DNA Lattice for Finite Element Numerical Solution

Although one can find formulation of the mean field Poisson-Boltzmann equation in many standard textbooks for electrolytes, we would like to very briefly derive it here to introduce important notations and to transform it later to a form suitable for weak formulation of the finite element numerical solution.

2.1 Poisson-Boltzmann equation

For deriving Poisson-Boltzmann equation for DNA system in solution, we start from the Poisson equation for the electrostatic potential V(r) which is generated by a charge system:   

\begin{equation} \nabla^{2}V(\boldsymbol{{r}}) = -\frac{\rho_{f}(\boldsymbol{{r}}) + \rho(\boldsymbol{{r}})}{\varepsilon\varepsilon_{0}}. \end{equation} (1)

Here, ρ_f(r) is the fixed charge density (such as the immobile charge of DNA molecules considered in this work) and ρ(r) is the charge density of mobile ions in solution, ε₀ is the electric permeability of free space, and ε is the relative dielectric constant. If we have N ions at the position r_i and charge q_i then:   

\begin{equation*} \rho(\boldsymbol{{r}}) = \sum\nolimits_{i = 1}^{N}q_{i}\delta(\boldsymbol{{r}} - \boldsymbol{{r}}_{i}). \end{equation*}

In Poisson-Boltzmann mean-field theory, instead of explicitly examining the position of each ion over time, we use the average quantities according to Boltzmann statistics where the probability of a state is proportional to $e^{ - \frac{E}{k_{B}T}} = e^{ - \frac{qV}{k_{B}T}}$, where k_B $ \simeq $ 1.38 × 10⁻²³ J/K is the Boltzmann constant and T is the temperature of the system. Specifically, one replaces the charge density with its statistical mean:   

\begin{align} \rho(\boldsymbol{{r}}) &= \rho_{+}(\boldsymbol{{r}}) + \rho_{-}(\boldsymbol{{r}}) \\ &= (\mathit{ez}_{+}n_{0+})e^{-z_{+}eV(\boldsymbol{{r}})/k_{B}T} + (-\mathit{ez}_{-}n_{0-})e^{+z_{-}eV(\boldsymbol{{r}})/k_{B}T}, \end{align} (2)

where e $ \simeq $ 1.6 × 10⁻¹⁹ C is the elementary charge, n₀₊, ez₊ and n₀₋, −ez₋ are bulk concentrations and valence of positive ions (counter-ions of DNA) and of negative ions (co-ions of DNA), respectively. For example, in aqueous solution with NaCl salt, z_± = 1; with MgCl₂ salt, z₊ = 2, z₋ = 1, …. At infinity (in a solution far from a DNA system), the electrostatic potential V(r) and the average charge density ρ(r) must equal zero. This leads to condition   

\begin{equation} z_{+}n_{0+} = z_{-}n_{0-} \end{equation} (3)

Combining Poisson equation, eq. (1) with Boltzmann statistic mean charge density, eq. (3), we obtain the famous mean-field Poisson-Boltzmann equation for DNA system in solution with mobile ions as:   

\begin{align} \nabla^{2}V(r) &= -\frac{1}{\varepsilon\varepsilon_{0}}[\rho_{f}(\boldsymbol{{r}})V(\boldsymbol{{r}}) + \mathit{ez}_{+}n_{0+}e^{-\frac{z_{+}eV(\boldsymbol{{r}})}{k_{B}T}} \\ &\quad - \mathit{ez}_{-}n_{0-}e^{+\frac{z_{-}eV(\boldsymbol{{r}})}{k_{B}T}}]\\ & = -\frac{1}{\varepsilon\varepsilon_{0}}\{\rho_{f}(\boldsymbol{{r}})V(\boldsymbol{{r}}) + \mathit{ez}_{+}n_{0+}[e^{-(\frac{z_{+}}{z_{-}})z_{-}eV(\boldsymbol{{r}})/k_{B}T} \\ &\quad- e^{+z_{-}eV(\boldsymbol{{r}})/k_{B}T}]\} \end{align} (4)

In practice, we consider fixed charges ρ_f as boundary conditions and solve this equation to find the charge density and the electrostatic potential in the solution. Then Poisson-Boltzmann equation becomes:   

\begin{equation} \nabla^{2}V(\boldsymbol{{r}}) = -\frac{\mathit{ez}_{+}n_{0+}}{\varepsilon\varepsilon_{0}}[e^{-(z_{+}/z_{-})z_{-}eV(\boldsymbol{{r}})/k_{B}T} - e^{+z_{-}eV(\boldsymbol{{r}})/k_{B}T}], \end{equation} (5)

provided that the boundary is near a surface with a charge density per σ unit area   

\begin{equation} \hat{n}\nabla V(\boldsymbol{{r}}) = -E = -\frac{\sigma}{\varepsilon\varepsilon_{0}} \end{equation} (6)

where $\hat{n}$ is the normal unit vector perpendicular to the charged surface facing the solution.

Rewrite eq. (5) using dimensionless quantities $\bar{V} = ez_{ - }V/k_{B}T$, we have:   

\begin{align} \nabla^{2}\bar{V}(r) &= -\frac{e^{2}z_{+}z_{-}n_{0+}}{\varepsilon\varepsilon_{0}k_{B}T}[e^{-(z_{+}/z_{-})\bar{V}} - e^{\bar{V}}] \\ &= -\bar{n}_{0+}[e^{-Z\bar{V}} - e^{\bar{V}}], \end{align} (7)

with   

\begin{equation*} Z = \frac{z_{+}}{z_{-}};\bar{n}_{0+} = \frac{e^{2}z_{+}z_{-}n_{0+}}{\varepsilon\varepsilon_{0}k_{B}T} = 4\pi l_{B}Zn_{0+};l_{B} = \frac{e^{2}z_{-}^{2}}{4\pi\varepsilon\varepsilon_{0}k_{B}T} \end{equation*}

where, l_B is called the length Bjerrum length, the length scale at which the electrostatic interaction between two elementary charges is equal to the thermal energy k_BT. In typical aqueous solution at room temperature, l_B $ \simeq $ 0.71 nm.

The boundary condition at fixed charge surfaces in dimensionless quantities   

\begin{equation} \hat{n}\nabla\bar{V} = -\frac{\mathit{ez}_{-}\sigma}{\varepsilon\varepsilon_{0}k_{B}T} = \bar{\sigma} \end{equation} (8)

From now on, we will work exclusively with dimensionless quantities. Without any confusions, one can remove the overline symbol above their notations in the rest of the paper.

2.2 Weak formulation of Poisson-Boltzmann equation for numerical finite element solution

In this work, we use numerical finite element method to solve for the charged and ion density in a hexagonal DNA lattice. The computer simulation process is programmed by using FreeFEM++ software package.²⁴⁾ This is a popular 2D and 3D partial differential equations (PDE) solver that allows us to easily implement our own physics modules using the provided FreeFEM++ language with a large list of finite elements. For this program, one needs the weak formulation of our partial differential equation.

From the point of view of numerical calculation, one can see that Poisson-Boltzmann equation is the result of the minimization of the following free energy functional.   

\begin{align} F(V) &= \int_{\varOmega}dr\left\{\frac{1}{2}|\nabla V|^{2} + \rho_{f}V + \frac{1}{Z}(n_{0+}e^{-ZV} + n_{0-}e^{V})\right\} \\ & = \int_{\varGamma}ds\sigma V_{\varGamma}\\ &\quad + \int_{\varOmega}dr\left\{\frac{1}{2}|\nabla V|^{2} + \frac{1}{Z}(n_{0+}e^{-ZV} + n_{0-}e^{V})\right\} \end{align} (9)

where Ω is the region of space where the ions live, Γ is the surface of DNA and the fixed charge density is ρ_f = σδ(r − r_Γ).

Indeed, a variation, V = V + δV gives   

\begin{align} \delta F & = \int_{\varGamma}\sigma\delta V_{\varGamma} + \int_{\varOmega}\nabla V\nabla \delta V - n_{0+}e^{-ZV}\delta V + \frac{n_{0-}}{Z}e^{V}\delta V\\ & = \int_{\varGamma}\sigma \delta V_{\varGamma} + \int_{\varOmega}\nabla (\delta V\nabla V) \\ &\quad - \delta V\nabla^{2}V - n_{0+}[e^{-ZV} - e^{V}]\delta V\\ & = \int_{\varGamma}(\sigma - \hat{n}\nabla V)\delta V_{\varGamma} + \int_{\varOmega}\nabla(\delta V\nabla V)\\ &\quad - \delta V\nabla^{2}V - n_{0+}[e^{-ZV} - e^{V}]\delta V\\ & = \int_{\varGamma}(\sigma - \hat{n}\nabla V)\delta V_{\varGamma} \\ &\quad + (-\nabla^{2}V - n_{0+}[e^{-ZV} - e^{V}])\delta V \end{align} (10)

The requirement that this variation is zero for all δV immediately leads to the Poisson-Boltzmann equation, eq. (7), with proper boundary condition at the fixed charge surface, eq. (8).

We will use this minimized functional to solve iteratively for V. Starting from a guess V₀, we evolve the solution in time according to   

\begin{align} \frac{\partial V}{\partial t} &= -\frac{\delta F}{\delta V} = (-\sigma + \hat{n}\nabla V)\delta (r - r_{\varGamma}) \\ &\quad + \nabla^{2}V + n_{0+}[e^{-ZV} - e^{V}]. \end{align} (11)

Therefore,   

\begin{align} \frac{\partial F}{\partial t} & = \int_{\varGamma}\sigma V_{\varGamma}' + \int_{\varOmega}\nabla V\nabla V' \\ &\quad + \frac{1}{Z}[n_{0+}e^{-ZV}(-ZV') + n_{0-}e^{V}V']\\ & = \int_{\varGamma}\sigma V_{\varGamma}' + \int_{\varOmega}\nabla[V'\nabla V] - V'\nabla^{2}V \\ &\quad + \frac{1}{Z}[n_{0+}e^{-ZV}(-ZV') + n_{0-}e^{V}V']\\ & = \int_{\varGamma}\sigma V_{\varGamma}' + \int_{\varOmega}\nabla [V'\nabla V] \\ &\quad - V'\left\{\nabla^{2}V - \frac{1}{Z}[n_{0+}e^{-ZV}(-Z) + n_{0,-}e^{V}]\right\}\\ & = -\int_{\varOmega}(V')^{2} \leq 0. \end{align} (12)

Hence the free energy functional F[V] is guaranteed to decrease overtime to its minimum. Note that in the last equality of this equation, the time evolution equation, eq. (11), was used, and that the electrostatic potential and electric fields are assumed to be zero at infinity.

To formulate the weak formulation for numerical calculation at each step, we start from the equation that needs to solved numerically,   

\begin{equation} \frac{\partial V}{\partial t} = \nabla^{2}V + n_{0+}[e^{-ZV} - e^{V}] \equiv \nabla^{2}V + f, \end{equation} (13)

with the boundary condition $\sigma = \hat{n}\nabla V$. Multiply eq. (13) with a test function, ω(r), and integrate over all space we have   

\begin{equation*} \int_{\varOmega}d\varOmega\frac{\partial V}{\partial t}w - \frac{1}{4\pi l_{B}}\int_{\varOmega}d\varOmega w\nabla^{2}\bar{V} = \int_{\varOmega}d\varOmega wf \end{equation*}

  

\begin{align*} &\int_{\varOmega}d\varOmega \frac{\partial V}{\partial t}w - \frac{1}{4\pi l_{B}}\int_{\varOmega}d\varOmega \nabla(w\nabla \bar{V}) \\ &\quad + \frac{1}{4\pi l_{B}} \int_{\varOmega}d\varOmega(\nabla w)(\nabla\bar{V}) = \int_{\varOmega}d\varOmega wf \end{align*}

  

\begin{align*} &\int_{\varOmega}d\varOmega\frac{\partial V}{\partial t}w - \frac{1}{4\pi l_{B}}\int_{\varGamma}w\hat{n}\nabla\bar{V} + \frac{1}{4\pi l_{B}}\int_{\varOmega}dr(\nabla w)(\nabla \bar{V}) \\ &\quad = \int_{\varOmega}d\varOmega wf \end{align*}

which must be satisfied for all ω(r). This is the weak formulation of our equation.

This weak formulation plus the charge neutrality condition, eq. (3), unfortunately, is not enough to obtain unique solution. We put further constrains on the numerical ion density to guarantee the uniqueness of the solution. These are done by making sure at every step that the total number of negative ions is fixed at the value equaling the volume times the expected concentration. The requirement that the whole system is charge neutral consequently leads to the fact that the total number of positive ions in the system must be so that their total charge equals the sum of both DNA charges and negative ions charges. During implementation, these constraints are imposed at the end of each iteration by rescaling the densities. As the solution converges, the rescaling factor approaches the value of 1.

Lastly, for the initial guess for numerical calculation we will let this be Debye-Huckel solution of a charged cylinder at infinite dilution.   

\begin{equation} V_{0}(r) = -\frac{2}{b\kappa_{D}}\frac{K_{0}(\kappa_{D}r)}{K_{1}(\kappa_{D}R)} \end{equation} (14)

with κ_D the inversed Debye-Huckel screening length, and b the Gouy-Chapman length. The former is found by linearizing the Poisson-Boltzmann equation,   

\begin{align} \nabla^{2}V_{0}(r) &= -n_{0}[e^{-ZV_{0}(r)} - e^{V_{0}(r)}] = n_{0}(Z + 1)V_{0}(r) \\ &= -\kappa_{D}^{2}V_{0}(r), \end{align} (15)

so   

\begin{equation} \kappa_{D} = \sqrt{n_{0}(Z + 1)}. \end{equation} (16)

The Gouy-Chapman length is matched by differential with respect to r:   

\begin{equation} \frac{\partial V_{0}}{\partial r} = \frac{2}{b}\frac{K_{1}(\kappa_{D}r)}{K_{1}(\kappa_{D}R)} \end{equation} (17)

which at r = R equals σ. So   

\begin{equation*} b = 2/\sigma. \end{equation*}

It should be noted that this initial guess is for the case of infinite dilution (d → ∞) and small surface charge density (σ → 0). It satisfies the Neumann boundary condition (correct gradient) at the DNA surface, but it does not satisfy other boundary conditions. Nevertheless, with our numerical procedure, all boundary conditions will be satisfied after just one iteration. In fact, any initial guesses for the solution will work. This particular guess is chosen so that we can start from a right-order-of-magnitude initial solution with correct gradient at the DNA surface. This helps with faster convergence of the numerical procedure.

3. Numerical Setup, Results and Discussion

3.1 System parameter setup

For numerical calculation, we assume the DNA molecules in the condensate to arrange in a two-dimensional hexagonal lattice with lattice constant d (see Fig. 1). The DNA axis is parallel to the z-axis. The space-filling periodic Wigner-Seitz cell of the lattice is the hexagonal region in the horizontal (x,y) plane. Due to symmetry, we only need to numerically solve the Poisson-Boltzmann equation for the green region in this Figure, which is 1/12 of the hexagonal cell (see Fig. 1). The individual DNA molecule is modeled as an impenetrable cylinder with radius of R = 1 nm. This corresponds to the radius of real B-DNA molecule (the most common form of DNA molecule). The linear charge density of B-DNA is −2e for every base pairs along the DNA axis, where e = 1.6 × 10⁻¹⁹ C is the elementary charge. Since the nucleotide base pairs distance along the DNA axis is 0.34 nm, the surface charge density of our cylinder is set to σ = −0.936e/nm².

The rest of the Wigner-Seitz cell can be obtained from the green region by mirror or rotational symmetry. The solvent water is treated as a continuous dielectric medium with dielectric constant ε = 78 and the system is assumed to be at room temperature T = 300°K. At the non-DNA boundaries of the green region, the normal derivative of the electrostatic potential is set to zero due to symmetry.

In Fig. 2(a), a typical triangle mesh that we generated for finite element solution of the Poisson-Boltzmann equation is shown. The maximum length of a side of a triangle element is 0.001 nm. The basis element function used is the second order P2 polynomial. The mesh is dynamically adapted to the solution during its converging procedure. This dynamical adaption of the mesh allows one to have finer mesh density at location where the gradient of the solution is highest, leading to more accurate discretization. However, this mesh adaptation is carried out only within the first 100 iterations. After that, the mesh is fixed.

Fig. 2

(a) A typical triangle mesh of the green region in Fig. 1 used for solving the Poisson Boltzmann equation. This mesh is dynamically adapted to the numerical solution during the numerical procedure; (b) The free energy functional F[V], eq. (9), decreases after each iteration as the solution converges after about 1000 steps.

In Fig. 2(b), the typical behavior of the minimizing free energy functional F[V], eq. (9), is plotted as a function of iterations in numerical finite element calculation. The stopping criteria of convergence of free energy functional is chosen to be that the relative change in free energy is less than 10⁻⁷. One can see clearly from this graph that our iterative minimization protocol works well, and the free energy functional indeed decreases overtime to its final value as the eq. (12) dictated. From these graphs, one can see that our solution of the non-linear Poisson Boltzmann equation does converge after about 1000–2000 iterations.

In Fig. 3, a typical numerical solution for the electrostatic potential of the Poisson-Boltzmann for our system is plotted using contour plot. Only one space-filling Wigner-Seitz cell is plotted. As expect, near the surface of the DNA, the isopotential lines are circular. The lowest, most negative potential is at the DNA surface. As one moves away from the DNA, the isopotential lines are more hexagonal to adapt to the Wigner-Seitz cell geometry. Since the normal derivative of the potential at the cell boundary must be zero by symmetry, the isopotential lines cross the cell boundary at perpendicular direction. The highest value of the electrostatic potential is at the corner of the cell, where the distance from the DNA is furthest. This is simply because the DNA molecules are negatively charged.

Fig. 3

A contour plot of the typical electrostatic potential, $\bar{V}(\boldsymbol{{r}})$, obtained from numerical solution to the Poisson-Boltzmann equation, eq. (7), for the hexagonal DNA system. Here is the solution for the case d = 5 nm in 100 mM salt concentration. Only one space-filling Wigner-Seitz cell is shown.

3.2 Ion distribution profile

Let us now investigate the distribution of ions for 1:1 salt at 100 mM, a typical physiological NaCl concentration. The results are plotted in Fig. 4. In these figures, the ion molar concentrations are plotted as a function of the distance from the central DNA axis and along the longest edge of our integration region, from r = 1 nm to $r = \sqrt{3} d/2$.

Fig. 4

(a) Distribution of negative ions molar concentration, ρ₋(r), as function of the distance from DNA surface, for different lattice constant, d, of the hexagonal DNA lattice; (b) Distribution of positive ions molar concentration, ρ₊(r), as function of the distance from DNA surface, for different lattice constant, d, of the hexagonal DNA lattice.

The left figure is for the negative ion (co-ions) concentration, the right figure is for the positive ion (counter-ions) concentration. In all the figures, the concentration profiles are plotted for different lattice constants of the hexagonal DNA lattice.

One immediate conclusion one can draw from Fig. 4(a) is that the negative ion concentration increases as one moves away from the DNA molecules. This is obviously the results of the fact that DNA molecules is negatively charged, which repels the negative mobile ions away. This repulsion has to compete with the entropic confinement of the ions within the DNA lattice. For smaller DNA lattice constant, d, the ions are confined to a smaller volume, therefore the saturation values of the ion concentrations at the mid-point are higher, which is as much as 60% higher at d = 3 nm as compared to the bulk concentration of 100 mM. As d increases to 15 nm, this saturated concentration is only 10% higher than the bulk concentration. Note that the negative ion concentration crosses the bulk value of 100 mM at some distance r^* less than d/2. At this value, the potential is zero (bulk potential). As one moves away from r^* and further into the middle space between DNA, the concentration is higher, meaning the electrostatic potential is “positive”. This should not be confused with the charged of DNA becoming positive (overcharged). The positive potential value here is simply because of the confinement of negative ions which DNA repels and pushes them toward the middle. As d increases further toward large value, the negative ions concentration in the middle space among DNA decreases toward the bulk value of 100 mM. The profile also becomes flatter as one would expect. Qualitatively, one can see that DNA influences the negative ion concentration up to a distance of about 4–5 nm from its surface. One easily identifies this length to be of the order of the well-known Gouy-Chapman length, eq. (18), which for our charged DNA is about 4.19 nm.

Vice versa the density of positive charges, in Fig. 4(b), decreases monotonically with increasing r and approach a saturation value as r approaches the middle region farthest from each DNA, at distance $r = \sqrt{3} d/2$. The concentration of positive ions is always greater than that of negative ions at every value of r (note the scale of the vertical axis in the two figures, Fig. 4(a) and 4(b) are different). Also, the concentration of positive ions is always greater than the bulk value of 100 mM as well. Since both negative and positive ions are of the same magnitude of charge valence, the difference in the concentrations is only due to excess counterions needed to neutralize DNA. Since n₊(r) > n₋(r) at all point r of our integration region, we conclude that the excess counterions of DNA are distributed all over the space. This confirms that the DNA is not overcharged at this concentration at the level of mean-field theory.

The next observation one can see from the positive ion distribution profile is the influence of the entropic confinement of the ions as the inter-DNA distance d decreases. Similar to the negative ion distribution profiles, when d decreases, the concentration saturation value increases as the ions are becoming more confined. Since the bulk concentration of positive ions is 100 mM = 0.1 M, the positive ions in Fig. 4(b) come mostly from the excess DNA counterions. One can see from this figure that for d < 10 nm, the difference in the positive ion distributions at large r is minimal. For both negative and positive ions distribution, entropic confinement is important when DNA molecules are within 10 nm from each other.

To further elucidate the role of entropic confinement, in Fig. 5, the contact ion densities of the negative and positive ions as a function of DNA lattice constant are plotted. Contact ion densities are the densities of the ions right next to the DNA surface. One can see that the concentration of ions, both of negative and positive charges, decreases with increasing lattice constant d and saturated when d is increased above 10 nm. This once again confirms that the entropic confinement leads to the accumulations of both negative and positive ions at the DNA surface. The entropic confinement is very strong for d < 10 nm, where even the negative ions (which is supposed to be repelled from DNA surface) can accumulate at the DNA surface to nearly 30 mM (one third of the bulk value of 100 mM). As d increases, this accumulation decreases toward zero. For positive ions and d > 10 nm, the contact ion density saturates at about 4.8 M. As d decreases, the contact density starts to increase and reaches as high as 6.7 mM at d = 3 nm, a significant 50% increases. This entropic confinement effect is usually neglected in discussion about electrostatics of DNA. For a DNA bundle, our results show that it is a significant effect, and should be taken care of for proper description of DNA electrostatics.

Fig. 5

(a) Contact concentration, ρ₋(R) of negative ions on the surface of DNA, plotted as function of the hexagonal DNA lattice constant, d; (b) Contact concentration, ρ₊(R) of positive ions on the surface of DNA, plotted as function of the hexagonal DNA lattice constant, d.

Going back to the ion distribution profile near DNA in Fig. 4(b) for d ≥ 10 nm, we see that the majority of DNA excess counterions condensed near the DNA surface within a surface distance of about 1 nm. One easily identifies this length with the radius of the DNA cylinder. These agree qualitatively with the Manning condensation picture.²⁵⁾ In Manning theory, no condensation length is given since the line of charge is assumed infinitely thin. For our DNA with finite radius, one can approximate the condensation thickness to be of the same size as the radius of the DNA cylinder. From scaling argument, obviously this is a natural length scale for this condensation region, since below this length, one sits very near the surface of the DNA. Therefore, below this length, the DNA can’t be considered as an infinitely thin line of charges. For our DNA, we also obtained the condensed counterion on the DNA surface to be 4.8 M which is 48 times higher than the bulk concentration of 100 mM. These cannot be obtained from Manning phenomenological theory. Note that these behaviors agree well both qualitatively and quantitatively with results obtained by using Monte-Carlo simulation²⁶⁾ at low concentrations and small counterion valence. However, at high concentration or high counterion valence, the decreasing of negative charge concentration is not observed. This is because of neglect of strong correlation among counterions within mean-field theory.

4. Conclusions

In conclusion, a computational study of the electrostatics of DNA condensation in hexagonal bundle is presented. We derived a weak formulation of the Poisson-Boltzmann equation and solve it numerically using finite element method for ions distribution in the hexagonal DNA lattice. The component concentrations of negative ions and positive ions around DNA are obtained. The dependence of ion concentrations on the inter-DNA distance in the bundle is investigated. The results show strong influence of the entropic confinement of the ions when the distance between neighboring DNA decreases to below 10 nm. This pushes the ion concentrations near the surface of the DNA increase significantly, as much as 50%. This effect should be considered in any proper electrostatic investigation of DNA bundles. For the physiological monovalent salt of 100 mM concentration, our results show that the positive ion population is dominated by the excess counterions of DNA. These results are consistent qualitatively with previous simulation results in Monte-Carlo simulations. Quantitative deviation for high valence counterions needs more investigation.

The Poisson-Boltzmann equation is considered in mean-field theory and it fundamentally neglects fluctuations and correlations. Additionally, the physical parameters of the system such as the finite size of the ions is also neglected. The strength of DNA-DNA short-range attraction can be affected by this parameter.¹⁴⁾ Moreover, in fact the solvent is not a continuous media therefore discrete solvent corrections need to be added.²⁷⁾ These effects go beyond the purpose of this work.

In our computer simulation, we focus on the concentration of ions. It is believed that most important qualitatively behaviors will not change significantly for higher counterion, co-ion valency. However, in the cases of very high concentration of ions in solution, the strong deviations from Poisson-Boltzmann behavior is observed. At high concentrations, electrostatic interactions are strongly screened and the specificity of ions, the structure of the water shell around them (hydration shell) and the ionic finite size and polarizabilities need to take into account. There has been shown that the water shell around ions causes short-range attraction by molecular dynamic (MD) simulations.²⁸⁾ Another interesting effect for high electrolyte concentration is phase transitions and related critical phenomena are reviewed.²⁹⁾

Additionally, there are many other factors that were neglected in our simplified model. One approximation is that in the simulation, the position of the DNA cylinders is straight with infinite bending rigidity. Inside viruses, DNA is bent, and the configuration entropy of the DNA is not necessarily zero, and there is not a perfect hexagonal arrangement of DNA cylinder with fixed inter-DNA distance.³⁰⁾ The relative orientation of individual DNA cylinder is also fixed in our simulation. Allowing for various orientations of DNA charges can lead to a frustration effect similar to spin glass in a two-dimensional Ising model.³¹⁾ The dielectric discontinuity at the DNA surface is also another important factor that was neglected in this model and can have important consequences on the screening of DNA in different electrolytes.³²⁾

The corrections that take into account all these factors are beyond the scope of this paper and needs to be considered in more detail in future work. Nevertheless, we believe that the qualitative contribution of ion concentration discussed in this paper will not be affected.

Acknowledgements

This work is financially supported by the Vietnam National University, Hanoi through grant number QG.16.01.

REFERENCES

• 1)   S.  Surana,  A.R.  Shenoy and  Y.  Krishnan: Nat. Nanotechnol. 10 (2015) 741–747.
• 2)   H.  Chang and  S.-H.  Wu: Mater. Trans. 51 (2010) 2094–2098.
• 3)   H.  Chang and  W.-C.  Tzeng: Mater. Trans. 49 (2008) 1467–1473.
• 4)   C.M.  Knobler and  W.M.  Gelbart: Annu. Rev. Phys. Chem. 60 (2009) 367–383.
• 5)   M.R.  Jones,  N.C.  Seeman and  C.A.  Mirkin: Science 347 (2015) 1260901.
• 6)   K.M.  Abu-Salah,  A.A.  Ansari and  S.A.  Alrokayan: J. Biomed. Biotechnol. 2010 (2010) 715295.
• 7)   B.D.  Mai,  H.T.  Nguyen and  D.-Q.  Hoang: Mater. Trans. 60 (2019) 2499–2505.
• 8)   A.  Leforestier and  F.  Livolant: Proc. Natl. Acad. Sci. U.S.A. 106 (2009) 9157–9162.
• 9)   N.V.  Hud and  I.D.  Vilfan: Annu. Rev. Biophys. Biomol. Struct. 34 (2005) 295–318.
• 10)   G.M.  Torrie and  J.P.  Valleau: J. Chem. Phys. 73 (1980) 5807–5816.
• 11)   A.G.  Moreira and  R.R.  Netz: Eur. Phys. J. E 8 (2002) 33–58.
• 12)   V.  Vlachy and  A.D.J.  Haymet: J. Chem. Phys. 84 (1986) 5874–5880.
• 13)   L.  Guldbrand,  L.G.  Nilsson and  L.  Nordenski: J. Chem. Phys. 85 (1986) 6686–6698.
• 14)   T.T.  Nguyen: J. Chem. Phys. 144 (2016) 065102.
• 15)   L.  Degreve,  M.  Lozada-Cassou,  E.  Sanchez and  E.  Gonzalez-Tovar: J. Chem. Phys. 98 (1993) 8905–8909.
• 16)   L.  Degreve and  M.  Lozada-Cassou: Mol. Phys. 86 (1995) 759–768.
• 17)   M.  Deserno,  C.  Holm and  S.  May: Macromolecules 33 (2000) 199–206.
• 18)   T.T.  Nguyen and  B.I.  Shklovskii: Phys. Rev. Lett. 89 (2002) 018101.
• 19)   A.Yu.  Grosberg,  T.T.  Nguyen and  B.I.  Shklovskii: Rev. Mod. Phys. 74 (2002) 329.
• 20)   R.R.  Netz: Eur. Phys. J. E 5 (2001) 557–574.
• 21)   W.  Humphrey,  A.  Dalke and  K.  Schulten: J. Mol. Graphics 14 (1996) 33–38.
• 22)   S.  Genheden and  U.  Ryde: Expert Opin. Drug Discov. 10 (2015) 449–461.
• 23)  D. Andelman: Handbook of Biological Physics, (Elsevier Science B.V., Amsterdam, 1995) pp. 603–642.
• 24)   F.  Hecht: J. Numer. Math. 20 (2012) 251–266.
• 25)   G.S.  Manning: J. Chem. Phys. 51 (1969) 924.
• 26)   V.D.  Nguyen,  T.T.  Nguyen and  P.  Carloni: J. Biol. Phys. 43 (2017) 185–195.
• 27)   Y.  Burak and  D.  Andelman: J. Chem. Phys. 114 (2001) 3271–3283.
• 28)   E.  Guardia,  R.  Rey and  J.A.  Padro: Chem. Phys. 155 (1991) 187–195.
• 29)   Y.  Levin: Rep. Prog. Phys. 65 (2002) 1577–1632.
• 30)   P.K.  Purohit,  M.M.  Inamdar,  P.D.  Grayson,  T.M.  Squires,  J.  Kondev and  R.  Phillips: Biophys. J. 88 (2005) 851–866.
• 31)   G.M.  Grason: Phys. Rev. Lett. 105 (2010) 045502.
• 32)   A.  Naji,  M.  Ghodrat,  H.  Komaie-Moghaddam and  R.  Podgornik: J. Chem. Phys. 141 (2014) 174704.