Dielectric Spectroscopy of DNA molecules in solution

Abstract

With the abundance of genomic sequence data available nowadays, DNA plays a vital role in any scientific process in the field of biology, from gene therapy to drug discovery. A multitude of DNA binding proteins have been identified and their structures have been solved in complex with the DNA sequence they have preference for to bind. Yet, it remains still very hard to design a protein with engineered DNA binding preferences even when detailed knowledge of DNA–DNA and DNA–ligand structures can nowadays be obtained by high-resolution techniques (DNA structures down to 0.35 $\mathring{A}$ readily available in PDB). A good explanation of this situation, comes from the fact that energetic analysis of protein-DNA complexes, still remains experimentally and theoretically challenging, because it requires quantification of the different contributions to such interactions, in particular the electrostatic energy term.

Electrostatic forces have long been recognized to inherently influence the DNA structure and interactions including DNA bending and folding and DNA–ligand recognition, owing to the high charge density of the DNA molecule backbone, as well as the polar associative interactions between the nucleotide bases. Standard dielectric characterization tools, such as impedance spectroscopy and dielectrophoresis, only yield average values of DNA polarizability in bulk solution that include major secondary structural contributions and DNA-solvent interfacial effects (shielding). Latest experiments, calculate a value of $\epsilon_{DNA} ∼ 8.5$ , which differs substantially both from values measured by the aforementioned techniques as well as from standard theoretical models ( typically assuming DNA to be a low-polarizable medium with $\epsilon_{DNA} ∼ 2–4$ ).

This work, which is laid in the form of a self-imposed Problem Set, aims to understand frequency dependance of the dielectric constant of DNA in solutions, in order to inform modern computational and experimental DNA assays. It starts by diving into the fields of electrostatic polarization and dielectric relaxation with the hope (hold your breath) to formulate a theoretical model for the dielectric spectrum of polar materials, like DNA polymers. It continues with the use of state-of-the-art Molecular Dynamic simulations, in order to test the full-atom available DNA structures against the theoretical model. Lastly, an experimental validation of the theoretical model is proposed and hopefully will be soon realized.

Questions

1. What is the relationship between the complex dielectric constant of a homogeneous and isotropic material and:

   a) conductivity, impedance, permittivity?

   Hint: Model the material as a circuit of an ideal Resistor and Capacitor in parallel.

   b) dipole moment, polarization and susceptibility? Hint: View the material as an ensemble of electrical dipoles and use Maxwell's equations for Electrostatics

2. Derive the frequency dependent dielectric function from dipole moment auto-correlation. Hint: Use the Linear Response Theory and the Fluctuation-Dissipation Theorem to derive the dielectric relaxation constant applying an impulse and measuring the response.

3. Run a Molecular Dynamics simulation to measure the dipole moment auto-correlation function of the Drew-Dickerson dodecamer.

   a) Describe the simulation

   b) what is the contribution of each type of molecules (nucleic acids, ions, water) in the total dipole moment of the system ?

   c) Use the result from problem (3) to plot the spectrum of the complex dielectric constant of the DNA molecule. How does it compare to experimental values in literature ?

4. Run a Molecular Dynamics simulation to measure the dipole moment auto-correlation function of an ensemble of DNA molecules represented as rods with a permanent dipole moment.

Answers

1. a) Conductivity of can be defined as the ability of charges to move freely through a material. If the material has a large number of free electrons, with high mobility, it is called a conductor. Else, i.e. bound electrons or low mobility, it is called a dielectric.

   The equation defining conductance is:

   $$G = \frac{1}{R} = \frac{\sigma A}{L},$$

   where $G$ is conductance, $R$ is resistance, $A$ is the cross-sectional area and $L$ is the length of the conductor. Resistivity is defined as the ability of a material to impede the flow of charge. Consequently, dielectrics are materials with low conductivity and high resistivity.

   Permittivity $\epsilon$ describes how a material behaves in the presence of an electric field.

   $$\epsilon_r(t) = \frac{\epsilon(t)}{\epsilon_0}.$$

   It is directly related to the molecular properties of the material and further on we will associate it with polarizability. For now, let's assume a simple RC model of our material as the figure below.

   

   Conductivity and permittivity can be used to describe a materals ability to store energy and dissipate energy. Assuming a parallel plate capacitor model with a dielectric material inside,

   it's capacitance is given by

   $$C = \frac{\epsilon_0\epsilon_r(t)A}{d}.$$

   When an electric field is applied across a dielectric material placed between the capacitor’s plates, this field causes the dielectric material to change. At the atomic level, an applied electric field can pull the positively charged nucleus of an atom off center from the negatively charged electron cloud. These opposite charges now separated by a short distance establish their own electric field which opposes the applied electric field. This tiny opposing electric field is established by what is known as an induced dipole. The strength of thesedipoles directly affect the storage capacity of the capacitor. The greater the permittivity of the dielectric placed in the capacitor, the greater the storage capacity of the capacitor.

   The permittivity associated with a material is not a static constant, but is dependent upon the rate of change of the direction of the applied electric field. When the applied electric field is abruptly changed, forces acting on the material are shifted which causes the system to adjust to a new equilibrium point. While the electric field changes abruptly, the material's physical properties do not follow instantaneously. This delay explains why the permittivity of a material is dependent upon the frequency of these changes. This frequency dependence can best be described using the RC circuit above and the concept of impedance.

   Impedance is a measure of the voltage to current ratio in an element and oftenexpresses the resistance to flow of current. For a resistor, the impedance would be equivalentto the resistance which is simply defined as the voltage drop divided by the current passing through the element. The impedance of an ideal capacitor is dependent upon the rate atwhich current direction is changed while passing through the capacitor. If current directionnever changes then its impedance is infinite, and if it is changed at a very high rate then theimpedance will be very low. An ideal capacitor's impedance is related to the frequency of the applied electrical field and the capacitance as:

   $$Z_c=\frac{1}{j2\pi fC}.$$

   The impedance of a resistor is given as (inverse of conductance):

   $$Z_R=R=\frac{1}{G}=\frac{L}{\sigma A}$$

   The parallel impedance of the system abovecan be calculated and is given by

   $$Z_{RC}=(j2\pi fC+\frac{\sigma A}{L})^{-1}.$$

   Admittance $Y_{RC}$ is the inverse of impedance and defines the measure of the ease of electrical conduction, or

   $$Y_{RC} = (j2\pi fC+\frac{\sigma A}{L}) = (\frac{j2\pi f\epsilon_0\epsilon_r(t)A}{d} + \frac{\sigma A}{L}) =\frac{j2\pi f\epsilon_0A}{d}(\epsilon_r(t)+\frac{\sigma A}{j2\pi f\epsilon_0})$$

   The permittivity of a material is due to the physical changes induced in the material by the applied electric field. The effects which related to stored energy we will give the symbol $\epsilon'(t)$ and the effects which relate to energy losses we will give the symbol  $\epsilon''(t)$ . The sum of these two components describe the permittivity of a material. The physical properties of the material cannot change instantaneously and because of this they are described as having a magnitude and a phase. Assuming an exponential change we have:

   $$\epsilon_r(t) =\epsilon'(t)-j\epsilon''(t)$$

   where $\epsilon'(t)$ is related to the dielectric dispersion or energy storage component, and $\epsilon''(t)$ is related to the dielectric losses such as heating or friction. Substituting further above we get:

   $$Y_{RC}=\frac{j2\pi f\epsilon_0A}{d}(\epsilon'(t)-j\epsilon''(t)+\frac{\sigma}{j2\pi f\epsilon_0}+\epsilon_\infty).$$

   At low frequencies the conductivity, $\sigma$ , will dominate the admittance of the circuit and similarly the dielectric response. At high frequencies the effects of conductivity will go to zero and the overall permittivity will reach a limit which we will call  $\epsilon_\infty$ .

   b) An electric dipole is formed by fixing two opposite charges a certain distance from each other. This dipole has a moment $\mu$ , which is defined as the amount of charge $q$ of the dipole multiplied by the charge separation vector $\mathbf{d}$ between the distinct charges.

   $$\mathbf{\mu} =q\mathbf{d}.$$

   The SI unit of the dipole moment is $\text{ C}\cdot\text{m}$ . Yet, most often the unit Debye (D) is used, where $1D = 3.336 × 10^{–30} \text{ C}\cdot\text{m}$ .

   By introducing a material to an external electric field, dipoles (created or existing) inside the material, align in such a way as to oppose the applied electric field. This process is called polarization, and happens every time there is a change in the direction of the external electric field.

   When a dipole is polarized, it physically changes the material by shifting different elements (electrons, ions, molecules). These motions can be related to the resistive and capacitive circuit model of part a). The losses incurred by shifting the material, such as heating and vibration, can bemodeled using a resistor. The reoriented dipoles contain a dielectric potential similar to the electrical potential stored in the electric field of a capacitor.

   

   A material's polarization can be broken down into three components (see figure above), all of which are defined by the atomic structure of the molecular ensemble. The smallest component of a molecule's polarization is called electronic polarization, and it is due to slight distortions of the negatively charged electron cloud surrounding the positively charged nucleus. An order of magnitute up lies the ionic polarization, which is due to shifting of charged ions in a material's structure. The third and predominant component of polarization can be found in materials that do not have a symmetrical charge distribution. The asymmetrical charge distribution causes these material's molecules to produce a permanent dipole moment. This permanent dipole moment is aligned in the presence of an external electric field and this alignment produces the greatest effect on the overall polarization of a material.

   The system's (material or molecule) overall polarization is the sum of all the individual polarization vectors and is referred to as a polarization density. The macroscopic polarization $\mathbf{P}_{per}$ (permanent dipole moment) is defined as the overall dipole moment induced by an applied electric field per unit volume. Thus, if $\mathbf{M}=\sum\mu_i$ is the total dipole moment, we have

   $$\mathbf{P}_{per}=\frac{1}{V}\mathbf{M}=\frac{N}{V}<\mu>.$$

   Electronic polarization takes place on a time scale of $10^{–12} s$ because of the low mass of the electron. Atomic polarization takes place also at a slightly longer time scale. These effects can be summarized by an induced polarization $\mathbf{P}_\infty$ . Thus, the complete polarization density function can be written as

   $$\mathbf{P}=\frac{1}{V}\mathbf{M}+\mathbf{P}_\infty=\frac{N}{V}<\mu> + \mathbf{P}_\infty.$$

   Assuming that the dipoles do not interact with each other and that the electric field $\mathbf{E}_{Loc}$ at the location of the dipole is equal to the outer electrical field (no shielding effects), then the mean value of the dipole moment is is given only by the counterbalance of the thermal energy and the interaction energy $W$ of a dipole with the electric field given by $W = –\mu \cdot \mathbf{E}$ .

   $$<\mathbf{\mu}>=\frac{\int{4\pi}{\mu e^{\frac{\mu\cdot\mathbf{E}}{k_BT}}d\Omega}}{\int{4\pi}{ e^{\frac{\mu\cdot\mathbf{E}}{k_BT}}d\Omega}}$$

   For small values of interaction energy of a dipole with the electric field (field strengths $|E|\leq10^6 \text{Vm}^{-1}$ ) compared to thermal energy, the equation reduces to

   $$<\mu>=\frac{\mu^2}{3k_BT}\mathbf{E}$$

   Plugging this to the previous equation yields:

   $$\mathbf{P}=\frac{\mu^2}{3k_BT}\frac{N}{V}\mathbf{E} + \mathbf{P}_\infty,$$

   which gives the polarization density of a material in the presence of a static electrical field (steady state). Deviding the last equation with $\mathbf{E}$ gives:

   $$\epsilon_s = \frac{1}{3\epsilon_0}\frac{\mu^2}{k_BT}\frac{N}{V} + \epsilon_\infty$$

   where $\epsilon_s=\lim_{\omega \to 0}\epsilon'(\omega)$ is the electrostatic dielectric permittivity and $\epsilon_\infty\lim_{\omega \to \infty}\epsilon'(\omega)$ covers all contributions due to electronic and atomic polarization. Yet, this equation starts to fail when applied for polar associating liquids. The reasons are static orientation correlations between molecules, e.g. hydrogen bonding. Thus to find the effective dipole moment one has to take into account the entire molecular assembly. By employing statistical mechanics for an ensemble of such molecules (e..g DNA in solution), the contribution of the dipole momentum to the dielectric function is given by (Frölich & Kremer):

   $$\epsilon_s - \epsilon_\infty = \frac{1}{3\epsilon_0k_BT}\frac{\langle\sum_i\mu_i(0)\sum_j\mu_j(0)\rangle}{V},$$

   where $\langle\sum_i\mu_i(0)\sum_j\mu_j(0)\rangle$ is the static correlation function of the dipole moment fluctuations.

   Dielectric spectroscopy on polar associating liquids measures the effective dipole moment of such assemblies which can be greater or smaller compared to the dipole moment of the single molecule depending on the molecular structure.

   In the case of a time-dependent external electrical field, re-orientation (alignment) of dipoles in a material happens with every change of the field direction. The external field encourages the charged ends of the dipole to the lie parallel with the lines of the applied field, with opposite direction due to attractive forces. Thus, total effective electrical field in the material is less than the applied. This phenomenon is called dielectric displacement, and is given by the time-dependent Maxwell Equations. For small field strengths, the polarization, dielectric displacement and electrical field density functions are related with the following equation:

   $$\mathbf{P}(\mathbf{r},t)=\mathbf{D}(\mathbf{r},t)-\mathbf{E}(\mathbf{r},t)=\epsilon_0(\epsilon_r(t)-1)\mathbf{E}(\mathbf{r},t)=\epsilon_0\chi(t)\mathbf{E}(\mathbf{r},t),$$

   where $\chi_\epsilon(t) = \epsilon_r(t) -1$ is called the system's susceptibility and expresses the ability of a material to be polarized by an external electric field.

   With that equation set, we can now redefine permittivity as the dispersion of an applied electric field due to dipolar elements present in the system. The dielectric dispersion is a frequency dependent response of a medium. This dispersion is causal and can be best described using complex permittivity, which accounts for the effects in both magnitude and phase.

   $$\epsilon(t)^*=\epsilon_0\epsilon_r(t)\cos(2\pi ft+\theta)=\epsilon'(t)+j\epsilon''(t)$$

2. Dielectric Response from Dipole Moment Auto-Correlation

   Dipole moment auto-correlation function takes the dot product between a reference moment $\mathbf{M}(0)$ and a particular point in time $\mathbf{M}(t)$ , i.e.

   $$\mathbf{\Phi}(t) = \langle\mathbf{M}(0)\cdot\mathbf{M}(t)\rangle$$

   As time progresses the dipole moment of a molecule will drift away from its initial position at an exponential rate, with ever decreasing differences over time, till an average separation is achieved. To elaborate further, the dipoles in a solution will try continuously to align with an alternating electric field applied across the solution; the dipoles will rotate back and forth following the changing field. As the frequency of this applied field increases, drag affects inhibit the dipoles ability to orient with the changing field. At a particular frequency the dipoles will no longer be able to align with the field, and this frequency is related to the relaxation time constant $\tau$ .

   This same time constant can be extracted from a system that is not exposed to an external electric field. If a molecule was placed in a system of other molecules and held at room temperature the thermal energy would cause the different molecules to move and collide with each other. If the dipole moment of a polarized molecule was known at a specific point in time, and then evaluated at subsequent time periods, a drift is observed.

   The most simple ansatz to calculate the time dependence of the dielectric behaviour is the assumption that the change of the polarization is proportional to its actual value, or

   $$\frac{d\mathbf{P}(t)}{dt}=-\frac{1}{\tau_D}\mathbf{P}(t)$$

   This equation leads to the aforementioned exponential decay of the correlation function, with a time constant $\tau_D$

   $$\mathbf{\Phi}(t) = \langle\mathbf{M}(0)\cdot\mathbf{M}(t)\rangle \rightarrow A_0e^{-\frac{t}{\tau_D}}$$

   Whether the system is excited by an external electric field or not, at room temperaturea system will be dynamic. At any time $t$ with respect to the reference time, the system will have changed due to thermal disruption and forces felt from other dipoles in the system. With greater separation between the reference time $t_0$ and time $t$ , the directions of the dipoles become less correlated. This loss in correlation can be graphed to show the decaying exponential behavior similar to the decaying exponential seen in the auto-correlation of a dipole.

   Utilizing these correlation functions we can extract the relaxation time constants withthe decaying exponential equation shown in equation The Fourier-Laplace transform isa bilateral Laplace transform where the variable s has been replaced with jω, an it is definedas

   $$\mathcal{L}\{f(t)\}=\int_{-\infty}^\infty{f(j\omega)e^{j\omega t}dj\omega}$$

   Taking the Fourier-Laplace transform of the correlation functions from equations we can extract the susceptibility values of each of the individual components in the system. The susceptibility is extracted using the equation:

   $$\chi(\omega)=\frac{1}{3Vk_BT}\mathcal{L}\{-\dot{\Phi}(t)\}$$

   Subsequently, we can extract the permittivity from the equation in problem 1:

   $$\epsilon(\omega) = \epsilon_0(\chi(\omega)+1)$$

   Experimental measurement of the dielectric constant typicaly limits itself to the overall response of a sample to the applied external electric field. Yet, although the dielectric constant in an inhomogeneous system is a property of the state of the entire system, we have shown above that it is possible heuristically to separate the contributing components. Thereby, we can simultaneously calculate the contributions of a macromolecule and the explicit solvent to the total dielectric constant, performing Molecular Dynamics (MD) simulations and traching the fluctuations of the dipole moments. Yet, an obstacle to such an approach is that the simulation must be on the order of nanoseconds to attain convergence (see refs).

   Physically, this choice amounts to fixing or immobilizing the DNA, as for instance in a chromosome, and considering only the fluctuations. If the DNA is mobile, it, like al charged species, will yield a finite conductance in the long time limit.

   1. There is only one value for the dipole moment at each MD timestep; thus statistics can be accumulated only overtime, not over the number of atoms.
   2. The dipole moment changes with the orientation of the molecules and the rotation of molecules is usually the slowest physical process in conformationally stable molecular systems.

   Nevertheless, computational power nowadays allows for the MD of DNA molecules in the order of microseconds, so this method increasingly favorable.

3. a) Simulation

   Atomistic Molecular Dynamics simulations of the Drew–Dickerson (DD) dodecamer were conducted to compute the spectrum of the complex dielectric constant of DNA based on the equation derived in problem 2. The DD dodecamer (double stranded 5'CGCGAATTCGCG3')  is a prototypic B-DNA molecule, whose sequence-specific structure and dynamics have been investigated by many experimental and computational studies. 

   The simulations, and analysis of them, were performed using the following tools:

   • NAMD 2.11: a freeware molecular dynamics simulation package written using the Charm++ parallel programming model. The latest version supports GPU acceleration, which is great for calculating Ewald corrections of periodic boundary conditions (FFT transforms).
   • VMD 1.9.2: a molecular visualization program for displaying, animating, and analyzing large biomolecular systems using 3-D graphics and built-in scripting (useful plugin: catdcd)
   • Python 2.7.11: well… uhm… nothing to do with the snake

   The procedure followed to perform the simulations is depicted in the following flowchart:

   

   The DD dodecamer atom coordinate file ( $.pdb$ ) was accessed from RCSB with the name 1BNA. It was split into two files, one for each DNA strand using VMD and Tk Console writepdb function.

   A Protein Structure File ( $.psf$ ) was generated with a custom Psf Generation Script ( $.pgn$ ) and the VMD plugin psfgen, combining the coordinates with charges and masses, from a topology file. This file contains any structural and topographical information about nucleic acids and proteins in general. It issupplied by CHARMM and is optimized for nucleic acids. The topology file includes information about atomic masses and the connections in each residue (Adenine, Cytosine, Guanine or Thymine).

   The script to run with psfgen was the following:

    
   package require psfgen # require the psfgen plugin
   topology toppar/top_all36_na.rtf # specify the topology file (here DNA specific CHARMM topology)
   alias residue DG GUA # converting residue names from rcsb to vmd format
   alias residue DA ADE
   alias residue DC CYT
   alias residue DT THY
   segment DNA1 {pdb DD1.pdb} # 5'->3' strand
   segment DNA2 {pdb DD2.pdb} # 3'->5' strand
   coordpdb DD1.pdb DNA1 # load the coordinate files for each strand 
   coordpdb DD2.pdb DNA2
   guesscoord # guess the hydrogen atom coordinates that were omited in the pdb file
   writepsf DD.psf # write the output
   writepdb DD.pdb

   Subsequently, the script was executed through the VMD text interface.

   Since DNA cannot exist in a vacuum (due to the strong nevativecharged backbones) a solvate must be added. Therefore, water is needed in all directions surrounding the DNA. The DNA molecule was simulated in two solvent environments:

   1. a box of TIP3P water molecules with minimum salt ions (NaCl), Periodic Boundary Conditions, and Particle Mesh Ewald correction for long-range electrostatic interactions. (this is the model used in furhter calculations)
   2. a sphere of TIP3P water molecules with minimum salt ions (NaCl) with non-periodic boundary conditions (used as an initial test and benchmark for NAMD/VMD with DNA)

   The solvent environment were created with custom scripts provided in the NAMD tutorial files.

   DNA is strongly negatively charged which can be compensated by adding positive ions. Natrium ions $N a^+$ are chosen as the positive charged ions and are used to neutralize the system. Adding more $N a^+$ ions in combination with $Cl^−$ ions makes it possible to implement salt in the system. The autoionize package was used to add the ions to the system.

   The resulting models are shown below (water molecules are diffused and ions are enlarged for for clarity).

   

   To test the sanity of the above procedure and to also benchmark my laptop and the CPU/GPU node on MD simulations, a full-atom $5ns$ simulation was performed using the water sphere environment. The configuration file is explained below.

    
   #############################################################
   ## JOB DESCRIPTION                                         ##
   #############################################################
   ​
   # Minimization and Equilibration of DD DNA in a Water Sphere
   ​
   #############################################################
   ## ADJUSTABLE PARAMETERS                                   ##
   #############################################################
   ​
   structure          ../common/DD_ws_ion.psf
   coordinates        ../common/DD_ws_ion.pdb
   ​
   set temperature    300
   set outputname     DD_ws_eq
   ​
   firsttimestep      0
   ​
   #############################################################
   ## SIMULATION PARAMETERS                                   ##
   #############################################################
   ​
   # Input
   paraTypeCharmm      on
   parameters          ../toppar/par_all36_prot.prm 
   parameters          ../toppar/par_all36_carb.prm 
   parameters          ../toppar/par_all36_lipid.prm 
   parameters          ../toppar/par_all36_cgenff.prm 
   parameters          ../toppar/par_all36_na.prm
   parameters          ../toppar/toppar_water_ions.str
   parameters          ../toppar/stream/misc/toppar_ions_won.str
   temperature         $temperature
   ​
   # Force-Field Parameters
   exclude             scaled1-4
   1-4scaling          1.0
   cutoff              12.0
   switching           on
   switchdist          10.0
   pairlistdist        14.0
   ​
   # Integrator Parameters
   timestep            2.0  ;# 2fs/step
   rigidBonds          all  ;# needed for 2fs steps
   nonbondedFreq       1
   fullElectFrequency  2  
   stepspercycle       10
   ​
   # Constant Temperature Control
   langevin            on    ;# do langevin dynamics
   langevinDamping     1     ;# damping coefficient (gamma) of 1/ps
   langevinTemp        $temperature
   langevinHydrogen    off    ;# don't couple langevin bath to hydrogens
   ​
   # Output
   outputName          $outputname
   ​
   restartfreq         500     ;# 500steps = every 1ps
   dcdfreq             250
   outputEnergies      100
   outputPressure      100
   ​
   #############################################################
   ## EXTRA PARAMETERS                                        ##
   #############################################################
   ​
   # Spherical boundary conditions
   sphericalBC         on
   sphericalBCcenter   14.798333168029785 21.214048385620117 8.887923240661621
   sphericalBCr1       26.0
   sphericalBCk1       10
   sphericalBCexp1     2
   ​
   #############################################################
   ## EXECUTION SCRIPT                                        ##
   #############################################################
   ​
   # Minimization
   minimize            100
   reinitvels          $temperature
   ​
   run 2500 ;# 5ps

Benchmarking results are summarized in the following table. Obviously, further simulations were run on the CPU/GPU node.

System	Number of CPU cores	GPU Assist	days/ns	Runtime for 5 ns
Macbook Pro 13''	4	No	0.54	12.5 hours (est.)
Intel Phi w/ 2xTeslas K20c	24	Yes	0.076	2 hours (actual)

The following figures show the computation node doing his job.

For the dielectric analysis, the water box model was used due to better solutions for electrostatics using the Ewald corrections. The configuration file is explained below.

 
​
   #############################################################
   ## JOB DESCRIPTION                                         ##
   #############################################################
​
   # Minimization and Equilibration of DD DNA in a Water Box
​
   #############################################################
   ## ADJUSTABLE PARAMETERS                                   ##
   #############################################################
   structure          ../common/DD_wb_ion.psf
   coordinates        ../common/DD_wb_ion.pdb
​
   set temperature    300
   set outputname     DD_wb_eq
​
   firsttimestep      0
   #############################################################
   ## SIMULATION PARAMETERS                                   ##
   #############################################################
​
   # Input
   paraTypeCharmm       on
   parameters           ../toppar/par_all36_prot.prm 
   parameters           ../toppar/par_all36_carb.prm 
   parameters           ../toppar/par_all36_lipid.prm 
   parameters           ../toppar/par_all36_cgenff.prm 
   parameters           ../toppar/par_all36_na.prm
   parameters          ../toppar/toppar_water_ions.str
   parameters          ../toppar/stream/misc/toppar_ions_won.str
   temperature         $temperature
​
   # Force-Field Parameters
   exclude             scaled1-4
   1-4scaling          1.0
   cutoff              12.0
   switching           on
   switchdist          10.0
   pairlistdist        14.0
​
   # Integrator Parameters
   timestep            2.0  ;# 2fs/step
   rigidBonds          all  ;# needed for 2fs steps
   nonbondedFreq       1
   fullElectFrequency  2  
   stepspercycle       10
​
   # Constant Temperature Control
   langevin            on    ;# do langevin dynamics
   langevinDamping     1     ;# damping coefficient (gamma) of 1/ps
   langevinTemp        $temperature
   langevinHydrogen    off    ;# don't couple langevin bath to hydrogens
​
   # Periodic Boundary Conditions
   cellBasisVector1    42.0    0.   0.0
   cellBasisVector2     0.0  44.0   0.0
   cellBasisVector3     0.0    0   47.0
   cellOrigin          15.206551551818848 22.122560501098633 9.116369247436523
​
   wrapAll             on
​
   # PME (for full-system periodic electrostatics)
   PME                 yes
   PMEGridSpacing      1.0
​
   # Constant Pressure Control (variable volume)
   useGroupPressure      yes ;# needed for rigidBonds
   useFlexibleCell       no
   useConstantArea       no
​
   langevinPiston        on
   langevinPistonTarget  1.01325 ;#  in bar -> 1 atm
   langevinPistonPeriod  100.0
   langevinPistonDecay   50.0
   langevinPistonTemp    $temperature
​
   # Output
   outputName          $outputname
​
   restartfreq         500     ;# 500steps = every 1ps
   dcdfreq             250
   xstFreq             250
   outputEnergies      100
   outputPressure      100
   #############################################################
   ## EXECUTION SCRIPT                                        ##
   #############################################################
​
   # Minimization
   minimize            1000
   reinitvels          $temperature
​
   run 2500 ;# 5ps

b) Dipole Moment Analysis

The analysis of the trajectory files generated from the MD was done using VMD and the Dipole Moment Watcher Tool (Extensions->Visualize).

Concerning the calculation of the dielectric constant for the components with net charges (negatively charged backbone), the instantaneous center of mass of the DNA was chosen as the origin to eliminate the conductance component due to the DNA translational motion. This choice is sane, as we are using a model based fluctuation-dissipation to find the dielectric constant, which is based on thermal fluctuations (see end of problem 2).

The number of atoms for each type of component were calculated through the Tk Console of VMD using the following commands:

 
   set dna [atomselect top nucleic]
   set h2o [atomselect top waters]
   set ion [atomselect top ions]
   $dna num # Number of atoms in DNA molecule
   $h2o num # Number of atoms in solvent box
   $ion num # Number of ions

The total volume of the solvation box was:

The volume for each component was estimated as follows. Each atom was assumed to occupy a sphere with an $1 \mathring{A}$ radius. The volume of the ions is a straightforward calculation. For the volume of the water molecules, we assumed that the volume is the sum of the 3 atom volumes. The DNA volume was estimated as the volume of the box minus the other components.

The mean dipole moment was calculated within the following Python code snippet:

 
   import numpy as np
   filename = 'XXX_dipole.dat' # XXX is one of DNA, H2O and ION
   with open(filename, 'rb') as f:
       frames = []
       dm = []
       dm_mag = []
       header = f.readline()
       data = [map(float,filter(None,l.strip('\n').split(' '))[1:]) for l in f.readlines()]
       data = np.asarray(data)
       for i,d in enumerate(data):
           frames.append(i)
           dm.append([d[0],d[1],d[2]])
           dm_mag.append(d[3])
       frames = np.asarray(frames)
       dm = np.asarray(dm)
       dm_mag = np.asarray(dm_mag)
​
   mean_dm = dm_mag.mean()

The results are shown in the table below.

Component	Number of Atoms	Estimated Volume	Mean dipole moment
DNA	$784$	$39.9nm^3$	$23.9184 D$
Water	$4806$	$20nm^3$	$175.7147D$ *
Ions	$24$	$100\mathring{A}^3$	$~0D$
Total	$5614$	$60nm^3$	$1370.3830D$ *

The values denoted by an asterisk in the table should actually be infinite because of the diverging displacements of the ions at progressively increasing times; however, these values are finite in the table because of the finite simulation time of $5.05275 ns$ . That the static dielectric constantis infinite is a quite general conclusion for any system in which the lowest resonance frequency for the charges is zero.

c) Plotting the Dielectric Spectrum

The dielectric spectrum was derived from the dipole moment fluctuations auto-correlation function, implementing the analysis described in problem 2 using the following Python script.

 
import numpy as np
import matplotlib.pyplot as plt
from scipy.optimize import fsolve
​
# Define constants
V_dna = np.float(39.9 * 10e-27) # Volume of DNA (see above for calculation)
k = np.float(1.3806488e-23) # Boltzmann
T = np.float(300) # Temperature of simulation
eps_0 = np.float(8.854187817e-12) # dielectric permittivity of air
adj_db = np.float(3.3356e-30) # Debye to C*m
A = 1.0/(3*V_dna*k*T) # precalculating formula constant
​
# Calculate and plot mean dot product autocorrelation
count = 0
max_lag = 500
data = np.zeros((max_lag,2))
for lag in xrange(max_lag):
    num_points = len(frames) - lag
    sum_dot_prod = 0
    for i in range(num_points):
        sum_dot_prod += np.dot(dm[i],dm[i+lag])
    mean_dot_prod = sum_dot_prod/num_points
    data[count,0] = lag
    data[count,1] = mean_dot_prod
    count += 1
​
freq = np.logspace(6,15,5*max_lag)
​
t = data[:,0]
phi = data[:,1]
​
plt.figure(1)
plt.xlabel('Seconds')
plt.ylabel('Mean-Dot Product of the Auto-Correlation')
plt.title('Auto-Correlation of DNA in a Water Box')
plt.plot(t*0.5e-12, phi)
f = plt.gcf()
f.set_facecolor('white')
​
# Calculate dielectric permittivity spectrum
dir_phi_p = np.diff(phi)/np.diff(t)
dir_phi_n = -1.0*dir_phi_p
​
hum = np.empty((len(t),len(dir_phi_n)), dtype=np.complex128)  #super careful with overflow
chi = np.empty(len(freq), dtype=np.complex128)
​
# Calculating susceptibility
for i, omega in enumerate(freq):
    for time in range(1, len(t)-1):
        hum[time] = A*adj_db*np.exp(1j*omega*t[time]*0.5e-12)*dir_phi_n[time]
    chi[i] = hum.sum()
​
# Calculating permittivity
eps = eps_0*(chi+1)
​
# Plotting
plt.figure(2)
ax1 = plt.subplot(211)
plt.ylabel('Real Permittivity')
plt.title('Dielectric Spectrum of Drew-Dickerson Dodecamer')
plt.semilogx(freq, np.real(eps))
​
ax2 = plt.subplot(212, sharex=ax1)
plt.xlabel('Frequency [Hz]')
plt.ylabel('Imaginary Permittivity')
plt.semilogx(freq, np.imag(eps))
f = plt.gcf()
f.set_facecolor('white')

Due to larger atomic charges (and hence, larger dipole moments), the $\epsilon_{DNA}$ value obtained with CHARMM27 is slightly higher compared with experimental values in literature ( $\epsilon_{DNA} ∼ 8.5$ ), but still in good accordance with the experimental results.

References

1. Cuervo, A. et al. Direct measurement of the dielectric polarization properties of DNA. Proc. Natl. Acad. Sci. U.S.A. 111, E3624–30 (2014).
2. Abeyrathne, C. D., Halgamuge, M. N., Farrell, P. M. & Skafidas, E. An ab-initio computational method to determine dielectric properties of biological materials. Sci Rep 3, 1796 (2013).
3. Fröhlich, H. The extraordinary dielectric properties of biological materials and the action of enzymes. 1–5 (2012).
4. Broadband Dielectric Spectroscopy. (Springer Berlin Heidelberg, 2003). doi:10.1007/978-3-642-56120-
5. M A Young, B Jayaram, A. & Beveridge, D. L. Local Dielectric Environment of B-DNA in Solution:  Results from a 14 ns Molecular Dynamics Trajectory. The Journal of Physical Chemistry B 102, 7666–7669 (American Chemical Society, 1998).
6. Yang, L., Weerasinghe, S., Smith, P. E. & Pettitt, B. M. Dielectric response of triplex DNA in ionic solution from simulations. Biophysical Journal 69, 1519–1527 (1995).
7. Alberto Pérez, F Javier Luque, A.Modesto Orozco. Dynamics of B-DNA on the Microsecond Time Scale. Journal of the American Chemical Society 129, 14739–14745 (American Chemical Society, 2007).
8. Honig B, Nicholls A (1995) Classical electrostatics in biology and chemistry. Science268(5214):1144–1149.
9. Sharp KA, Honig B (1990) Electrostatic interactions in macromolecules: Theory andapplications. Annu Rev Biophys Biophys Chem 19:301–332.
10. Cherstvy AG (2011) Electrostatic interactions in biological DNA-related systems. PhysChem Chem Phys 13(21):9942–9968.
11. Petrov AS, Harvey SC (2008) Packaging double-helical DNA into viral capsids: Struc-tures, forces, and energetics. Biophys J 95(2):497–50
12. Fenley MO, Harris RC, Jayaram B, Boschitsch AH (2010) Revisiting the association of cationic groove-binding drugs to DNA using a Poisson-Boltzmann approach. Biophys J 99(3):879–886
13. Rohs R, et al. (2009) The role of DNA shape in protein-DNA recognition. Nature 461(7268):1248–1253.
14. Privalov PL, Dragan AI, Crane-Robinson C (2011) Interpreting protein/DNA interactions: Distinguishing specific from non-specific and electrostatic from non-electrostatic components. Nucleic Acids Res 39(7):2483–2491.
15. Davis ME, McCammon JA (1990) Electrostatics in biomolecular structure and dynamics. Chem Rev 90(3):509–521.
16. Jayaram B, Sharp KA, Honig B (1989) The electrostatic potential of B-DNA. Bio-polymers 28(5):975–993.
17. Stigter D (1998) An electrostatic model for the dielectric effects, the adsorption ofmultivalent ions and the bending of B-DNA. Biopolymers 46(7):503–516.
18. Lopes PEM, Roux B, Mackerell AD, Jr (2009) Molecular modeling and dynamics studieswith explicit inclusion of electronic polarizability. Theory and applications. Ther Chem Acc 124(1-2):11–28.
19. van der Maarel JRC (1999) Effect of spatial inhomogeneity in dielectric permittivity on DNA double layer formation. Biophys J 76(5):2673–2678.
20. Cherstvy AG (2007) Effect of a low-dielectric interior on DNA electrostatic response to twisting and bending. J Phys Chem B 111(44):12933–12937.
21. Bone S, Lee RS, Hodgson CE (1996) Dielectric studies of intermolecular interactions in native DNA. Biochim Biophys Acta 1306(1):93–97.
22. Tomi c S, et al. (2007) Dielectric relaxation of DNA aqueous solutions. Phys Rev E Stat Nonlin Soft Matter Phys 75(2 Pt 1):021905–021913.
23. Hölzel R (2009) Dielectric and dielectrophoretic properties of DNA. IET Nanobiotechnol 3(2):28–45.
24. Ermolina I, Milner J, Morgan H (2006) Dielectrophoretic investigation of plant virus particles: Cow pea mosaic virus and tobacco mosaic virus. Electrophoresis 27(20):3939–3948.
25. Jin R, Breslauer KJ (1988) Characterization of the minor groove environment ina drug–DNA complex: Bisbenzimide bound to the poly[d(AT)]·poly[d(AT)]duplex. Proc Natl Acad Sci USA 85(23):8939–8942.
26. Barawkar DA, Ganesh KN (1995) Fluorescent d(CGCGAATTCGCG): Characterization of major groove polarity and study of minor groove interactions through a major groove semantophore conjugate. Nucleic Acids Res 23(1):159–164.
27. Pérez A, et al. (2007) Refinement of the AMBER force field for nucleic acids: Im-proving the description of α/γ conformers. Biophys J 92(11):3817–3829.
28. Pérez A, Luque FJ, Orozco M (2007) Dynamics of B-DNA on the microsecond timescale. J Am Chem Soc 129(47):14739–14745.
29. Dans PD, Pérez A, Faustino I, Lavery R, Orozco M (2012) Exploring polymorphisms inB-DNA helical conformations. Nucleic Acids Res 40(21):10668–10678.
30. Foloppe N, MacKerell AD (2000) All-atom empirical force field for nucleic acids:I. Parameter optimization based on small molecule and condensed phase macromo-lecular target data. J Comput Chem 21(2):86–104.
31. MacKerell AD, Banavali NK (2000) All-atom empirical force field for nucleic acids:II. Application to molecular dynamics simulations of DNA and RNA in solution.J Comput Chem 21(2):105–120.

Cool Stuff

Resources

• VMD/NAMD Tutorials
• Role of solvent in DNA simulations using NAMD/VMD
• PSF Generation for DNA molecules
• Introduction to MD simulation of DNA-protein systems
• VMD User's Guide
• Overview of methods to calculate Dielectric Constants
• Scripting in VMD
• NAMD/VMD Case Study: DNA
• VMD Images and Movies Tutorial

Gallery