NPanalysis - An analysis tool for two-component nanoparticles

NP(nanoparticle)analysis is a tool to analyze nanoparticle properties in all-atom or coarse-grained molecular dynamics systems. The scripts are written for DNA-PEI nanoparticles, but can be used for other two-component nanoparticles. It can currently only read Gromacs .gro and .xvg files generated by Gromacs to calculate different properties:

1. Connection matrix of molecules
2. Number of molecules performing different roles
3. Cluster analysis to calculate different molecules in a nanoparticle
4. Calculating the size of the nanoparticle in terms of number of molecules of one type (DNA or PEI)
5. Calculating size distribution of nanoparticles
6. Charge of nanoparticles as a function of size
7. Make nanoparticles whole across periodic boundaries
8. Radius of gyration, and hydrodynamic radius of nanoparticle as a function of time
9. Radius of gyration, and hydrodynamic radius of nanoparticle as a function of size
10. Calculate transition diagram and principal transition diagram of dynamic processes and plot them.

Published work

This toolbox is based on published work:

[1] Mahajan S., and Tang T., Polyethylenimine–DNA Ratio Strongly Affects Their Nanoparticle Formation: A Large-Scale Coarse-Grained Molecular Dynamics Study, J. Phys. Chem. B, 2019, 123, 45, 9629-9640, DOI: 10.1021/acs.jpcb.9b07031

[2] Mahajan S., and Tang T., Polyethylenimine–DNA Nanoparticles under Endosomal Acidification and Implication to Gene Delivery, Langmuir 2022, 38(27), 8382-8397, DOI: 10.1021/acs.langmuir.2c00952

Installation

NPanalysis can be installed as a module using

python setup.py install --user

Quickstart

For details check Tutorials.

Prepare Gromacs files

First make molecules whole across periodic boundaries and output .gro files:

gmx trjconv -f md_1.xtc -s md_1.tpr -pbc whole -o Whole/DP.gro -sep

Output the mass data:

gmx dump -s md_1.tpr > tpr.dump

Prepare and read constants

Generate a constants file 'constants.dat':

ndna = 27 # number of DNA 
npei = 270 # number of PEI
adna = 154 # number of DNA atoms
apei = 13 # number of PEI atoms
Qdna = -22 # Charge of DNA
Qpei = 3 # Charge of PEI
sdna = 1 # first atom of DNA
spei = 4159 # first atom of PEI
contact_dist = 0.5 # contact distance threshold
pbc = xyz

Read constants from the file and pickle it:

import NPanalysis as NPa
NPa.pickle_constants('constants.dat')

Generate molecule ID index or read it:

NPa.gmx.gen_index_mol('molndx.pickle')
NPa.gmx.read_index('index.ndx')

Read and pickle mass information:

NPa.gmx.pickle_mass('tpr.dump')

Calculate connection matrix

This can be calculated from any gromacs strucrture files, molecules does not need to be whole.

NPa.connMat.gro2connected(inGRO='Whole/DP')

Calculate PEI roles:

Average over 50 timesteps and output to 'PEI_roles.dat' and 'PEI_roles2.dat'. get_roles() calculates the number of PEIs in free, peripheral, bridging PEIs. get_roles2() calculates the average number of bridging PEIs between a pair of bridged DNAs and total number of bridged DNAs.

NPa.connMat.get_roles(50,'PEI_roles.dat')
NPa.connMat.get_roles2(50,'PEI_roles2.dat')

Determining clusters (nanoparticles)

Calculates clusters in 'cluster.pickle': connection matrix has to determined before this calculation.

NPa.cluster.gen_cluster()

Calculate cluster properties: Number average size, weight average size, number of nanoparticles as a function of number average size, charge of nanoparticles as a function of number average size.

NPa.cluster.gen_avgsize(50,'avgsize.dat')
NPa.cluster.gen_ncNP_s(3,'num_charge_NPs_size.dat')

Make nanoparticles whole

Input gromacs files must have whole molecules, connection matrix and clusters has to be calculated.

NPa.gmx.make_NPwhole(inGRO='Whole/DP',outGRO='NPwhole/DP')

Radius of gyration and hydrodynamic radius

Calculate radii for each cluster at different timestep: radius of gyration is stored as square. The input gromacs must have whole nanoparticles. Clusters has to be known.

NPa.radius.calc_Rh_Rg(Rh_pickle='Rh.pickle',Rg2_pickle='Rg2.pickle',inGRO='NPwhole/DP')

Calculate average radii: the average is performed over every 50 timesteps. Square root of average and standard error is taken for radius of gyrtion.

NPa.radius.gen_rad_avg('Rh.pickle',50,'avg_Rh.dat')
NPa.radius.gen_rad_avg('Rg2.pickle',50,'avg_Rg.dat',sqrt=True)

Calculate average radii per number average size: the average is calculated over timestems 200 throug 251 (not included). Square root of average and standard error is taken for radius of gyration.

NPa.radius.gen_rad_per_size('Rh.pickle','Rh_size.dat',200,251)
NPa.radius.gen_rad_per_size('Rg2.pickle','Rg_size.dat',200,251,sqrt=True)

Calculate principal transition diagram

Calculate principal transition diagram from connection matrix. In the code, "GPT" stands for "graph of principal transitions", which is identical to principal transition diagram.

NPa.cluster.gen_GPT(main_mol=0, GPT_out_pickle='GPT.pickle', UniqNSGs_out_pickle = 'UniqNSGs.pickle',
    connected_pickle='connected.pickle', time_pickle = 'time.pickle')

Plot the principal transition diagram using the unique nanoparticle IDs.

NPa.cluster.plot_GPT(GPT_pickle='GPT.pickle', show=True, edge_label=True, time_pickle='time.pickle')

Calculated nanoparticle structure graph from each unique ID, where number of bridges between each pair of main molecule (for example DNA) is evaluated. Below, NSG stands for "nanoparticle structure graph" and "weighed" or "w" indicates that number of bridges between main molecule is evaluated.

NPa.cluster.get_weighted_NSGs(main_mol=0, wNSG_pickle='wNSG.pickle', GPT_pickle='GPT.pickle', 
    UniqNSGs_pickle='UniqNSGs.pickle', time_pickle='time.pickle', connected_pickle='connected.pickle')

To generate publication quality images for NSG, you can use an interactive plot (See Tutorial 3). UID refers to the unique ID. Use the script below to determine the representation for all starting UIDs (or ending UIDs in case of dissociation).

NPa.cluster.plot_GPT_NSG(UID, GPT_pickle='GPT.pickle', wNSG_pickle='wNSG.pickle', 
    time_pickle='time.pickle', NSG_pos_pickle='NSG_pos.pickle', labels=True, ref_UIDs=None)

After the representation of all reference UIDs have been determined (using the command above), PNG of all NSG can be determined using the following command,

NPa.cluster.plot_GPT_NSG_all(GPT_pickle='GPT.pickle', wNSG_pickle='wNSG.pickle', 
    time_pickle='time.pickle', NSG_pos_pickle='NSG_pos.pickle', labels=False, arrow=True, rnd_off=4, 
    figname='GPT-NSG', dpi=200, ref_UIDs= UID_list)

Calculate transition diagram

The script is written such that one has to calculate the principal transition diagram before calculating the transition diagram. The advantage is purely asthetic. We can plot the transition diagram and highlight the principal transitions in the plot. See Tutorial 4 and ref [2] for details.

To calculate the transition diagram referred to as GT (graph of transitions), use the following python command,

NPa.cluster.gen_GT(main_mol=0, GT_out_pickle='GT.pickle', UniqNSGs_out_pickle = 'UniqNSGs.pickle', 
    connected_pickle='connected.pickle', time_pickle = 'time.pickle')

To plot the GT use the following python command,

NPa.cluster.plot_GT(GT_pickle='GT.pickle', GPT_pickle='GPT.pickle', time1=0, time2=0.15, 
    node_pos_pickle='GT_nodepos.pickle', iters=2000, show=True,  show_gt_edge=True)

iters = 0 plots the GT using neato networkx node positions. For a value greater than zero, updates the position of nodes in a molecular dynamics fashion.

Contributors

1. Subhamoy Mahajan, PhD student, Mechanical Engineering, University of Alberta, Edmonton, Canada (will graduate by Sep, 2021)
2. Tian Tang, Professor, Mechanical Engineering, University of Alberta, Edmonton, Canada

Funding Sources

v1.1 and v1.2: This research was funded by Natural Science and Engineering Research Council of Canada through Tian Tang. My research was also funded by several scholarships, Mitacs Globalink Graduate Fellowship (2016-2019), RR Gilpin Memorial Scholarship (2019), Alberta Graduate Excellence Scholarship (2020), Sadler Graduate Scholarship (2020).