Biomodeling: DNA in a Nanopore
Intro
Nanopore technology is a rapidly emerging field that utilizes nanoscale pores to characterize molecular identities.
For example, DNA/RNA sequencing, peptide sequencing, detection of small molecules, etc.
Materials for nanopore are diverse: solid-state, biological, 2D materials, etc.
References
Kalra, A., Garde, S. & Hummer, G. Osmotic water transport through carbon nanotube membranes. Proc Natl Acad Sci U S A 100, 10175-10180 (2003). https://doi.org:10.1073/pnas.1633354100
Hello, DNA!
Prepare DNA that is infinitely long.
Use our web server: https://yoo.skku.edu/apps/application/dna
If you want to simulate a DNA of finite lengths, you can prepare the moleculetype by using gmx pdb2gmx.
We made the web server because no MD package can generate DNA effectively infinite under PBC.
The web server generates dsDNA of any sequence.
The length of DNA should be close to a multiple of 10.5, which is the periodicity of dsDNA.
10, 11, 21 bp, 32 bp, 42 bp etc.
For example, try "atagcatcatatatataaaag"
AMBER force field has two major models for DNA
BSC0
BSC1
BSC1 is newer than BSC0. But, both are still actively used.
Outputs from web server
amber*.ff: force fields in Gromacs format.
conf.pdb: Structure file of dsDNA.
Confirm the periodicity by drawing periodic images.
topol.top: topology file ready to use.
*.mdp: MDP files for Gromacs.
Moleculetype definitions:
dna.itp: the essential moleculetype definition for you dsDNA.
dna.pbc.itp: contains the connectivity between two ends of strands.
dna.hbonds.itp: contains constraints between complementary pairs. Optional.
dna.enm.itp: contains ENM-like constraints. It can be useful during equilibration. Optional.
Perform MD of dsDNA in solution
The package is ready to run MD.
See dna.pdb file. The CRYST definition will look like this:
CRYST1 60.000 60.000 70.980 90.00 90.00 60.00 P 1 1
This indicates that you have a hexagonal simulation box: numbers indicates x, y, z, alpha, beta, gamma in the crystallographer's convention.
You can adjust the box size in x and y directions. But the z-dimension cannot be changed because it matches the periodicity of dsDNA.
Add water molecules by using gmx solvate.
gmx solvate -cp dna.pdb -o conf.pdb -p topol.top
Then, topol.top file will look like this:
; Include forcefield parameters
#include "amber99bsc0.ff/forcefield.itp"
;#define POSRES
#include "dna.itp"
#include "dna.pbc.itp"
#include "dna.hbonds.itp"
;#include "posre.itp"
;#include "dna.enm.itp"
; Include water topology
#include "amber99bsc0.ff/tip3p.itp"
; Include topology for ions
#include "amber99bsc0.ff/ions.itp"
[ system ]
; Name
loop in water
[ molecules ]
; Compound #mols
DNA_chain_A 1
SOL 7618
Add ions & minimize
gmx grompp
gmx grompp -f mini.mdp -c conf.pdb -maxwarn 40
gmx genion take TPR file as an input.
Add ions using gmx genion
gmx genion -neutral -conc 0.15 -o conf.pdb -p topol.top
-neutral: add ions such that the system becomes neutral.
-conc: add more ions to get 0.15 M.
See topol.top to confirm the changes.
Minimization
gmx grompp -f mini.mdp -c conf.pdb -maxwarn 40
gmx mdrun -c conf.pdb
the minimized structure will be written to conf.pdb
Gromacs mdrun put molecules in a rectangular box by default. If you want to see the hexagonal box, try this:
gmx trjconv -f conf.pdb -o conf.pdb -pbc mol -ur compact
Confirm the minimized system in VMD or pymol.
MDRUN
Copy mini.mdp to grompp.mdp
Modify grompp.mdp like this:
; RUN CONTROL PARAMETERS
integrator = md
; Start time and timestep in ps
tinit = 0
dt = 0.002
nsteps = 500000000
; Temperature coupling
tcoupl = v-rescale
print-nose-hoover-chain-variables = no
; Groups to couple separately
tc_grps = system
; Time constant (ps) and reference temperature (K)
tau_t = 0.5
ref_t = 300
Because we have a dsDNA along the z axis, the box size in z must be decoupled from x and y sizes. Similar to the lipid bilayer systems, this DNA system fits best with semi-isotropic box.
; pressure coupling
Pcoupl = c-rescale
Pcoupltype = semi-isotropic
; Time constant (ps), compressibility (1/bar) and reference P (bar)
tau-p = 1.0
compressibility = 4.5e-5 4.5e-5
ref-p = 1 1
gmx mdrun
gmx grompp -f grompp.mdp -c conf.pdb -maxwarn 40
gmx mdrun
Analysis
Visual inspection
Draw the concentrations of Na and Cl ions as a function of distance from DNA.
This is called the "counterion condensation" around DNA.
For analysis you need to coordinates of center of mass of dsDNA and coordinates of individual ions.
gmx traj is a Gromacs tool that print out text format coordinates from XTC files.
Center DNA in your simulation box.
In the output XTC file from MDRUN, DNA molecules move freely. For the analysis, having dsDNA at the center of hexagonal box is convenient.
gmx trjconv can do it:
gmx trjconv -f traj_comp.xtc -o center.xtc -pbc mol -ur compact -center
Select group for centering: DNA
Select group for output: 0
-f: input XTC
-o: output XTC.
-pbc mol -ur compact: put water and ions in hexagonal box.
-center: chosen molecules will be put at the center of box. By doing this, your trajectory will be DNA-centric.
Confirm center.xtc using VMD or pymol.
To print CM of dsDNA
gmx traj -s topol.tpr -f center.xtc -ox coord_dna.xvg -com
Choose DNA
-ox: Coordinates will be written to coord_dna.xvg file.
coord_dna.xvg file will contain lines like these:
# lines starting with "#" or "@" are comments.
# Each columns are time (ps) x y z (nm)
0 4.49242 2.58475 3.55809
20 4.60075 2.66173 3.73794
40 4.67276 2.70951 3.65495
60 4.6117 2.61924 3.74992
-com: With this option, center of mass coordinate will be written. Otherwise, coordinates of ALL atoms of DNA will be written.
To print coorinates of ions
gmx traj -s topol.tpr -f center.xtc -ox coord_na.xvg
Choose Na
Because we didn't specify -com option, coord_na.xvg will contain coordinates of all Na ions.
0 {x, y, z of ion 1} {x, y, z of ion 2} ...
20 {x, y, z of ion 1} {x, y, z of ion 2} ...
40 {x, y, z of ion 1} {x, y, z of ion 2} ...
60 {x, y, z of ion 1} {x, y, z of ion 2} ...
...
By using coord_dna.xvg and coord_na.xvg, calculate the distance between DNA and each Na ion, using your favorite tools such as numpy, Excel, etc.
Bin the distance to draw a histogram as a function of r.
Average over entire simulation time.
Then, you will see the ion density profile similar to the theory of counterion condensation.
Perform MD of dsDNA in a Graphene Nanopore (Not yet done)
Solvation
Minimization
You can use mini.mdp downloaded from the web server.
gmx grompp -f mini.mdp -r conf.pdb -c conf.pdb -maxwarn 40
gmx mdrun -c conf.pdb -nt 4
-c: output PDB
md.mdp – Pressure coupling
You can use md.mdp downloaded from the web server.
Mostly same options as other simulations.
Pressure coupling should be anisotropic because graphene sheet is anisotropic.
; pressure coupling
Pcoupl = berendsen
Pcoupltype = anisotropic
; Time constant (ps), compressibility (1/bar) and reference P (bar)
tau-p = 1.0
compressibility = 4.5e-5 4.5e-5 4.5e-5
ref-p = 1 1 1
md.mdp – Electric field
We want to apply a constant electric field along the z axis such that ions move through the pore.
The height of the box in the z axis is 7 nm for this example (could be different for you).
I want to apply a potential difference of 1 V.
Then, the electric field will be 1 V / 7 nm = 0.142857 V/nm.
Add the following lines to md.mdp
; Electric fields
; Format for electric-field-x, etc. is: four real variables:
; amplitude (V/nm), frequency omega (1/ps), time for the pulse peak (ps),
; and sigma (ps) width of the pulse. Omega = 0 means static field,
; sigma = 0 means no pulse, leaving the field to be a cosine function.
electric-field-x = 0 0 0 0
electric-field-y = 0 0 0 0
electric-field-z = 0.142857 0 0 0
Perform MD
gmx grompp -c conf.pdb -f md.mdp -r conf.pdb -maxwarn 40
gmx mdrun
Run for about 10 ns.
Load the trajectory to VMD and see ions moving through the pore.
The rate of ion movements is the ionic current.