Running a Basic NERDSS Simulation
After installing NERDSS, let's run a basic simulation to understand its philosophy. We will simulate the reversible reaction A(a) + R(r) <-> A(a!1).R(r!1) in a 3D solution. For a comprehensive guide to using NERDSS, please refer to user guide in the NERDSS repository.
Prepare the .inp file for the simulation
start parameters nItr = 700000 # steps timeStep = 0.1 # us timeWrite = 200 # steps pdbWrite = 10000 # steps checkPoint = 70000 # steps restartWrite = 10000 # steps end parameters start boundaries WaterBox = [939.993,939.993,939.993] # nm end boundaries start molecules A : 1000 R : 1000 end molecules start reactions A(a) + R(r) <-> A(a!1).R(r!1) onRate3Dka = 988.19 # 3D microscopic binding rate, nm^2/us offRatekb = 99.15 # microscopic dissociation rate, s^-1 norm1 = [0,0,1] norm2 = [0,0,1] sigma = 2.0 assocAngles = [nan,nan,nan,nan,nan] bindRadSameCom = 1.1 end reactions
Prepare the .mol files for the molecule types in the simulation
The following A.mol file is needed for the simulation and can be downloaded from here.
name = A # translational diffusion constants D = [20.0,20.0,20.0] # rotational diffusion constants Dr = [0.0,0.0,0.0] # Coordinates COM 0.0000 0.0000 0.0000 a 0.0000 0.0000 0.0000
The following R.mol file is also required for the simulation and can be downloaded from here.
name = R # translational diffusion constants D = [20.0,20.0,20.0] # rotational diffusion constants Dr = [0.0,0.0,0.0] # Coordinates COM 0.0000 0.0000 0.0000 r 0.0000 0.0000 0.0000
Please refer to user guide in the NERDSS repository for the explanation of each parameter.
Use the GUI provided in the ioNERDSS library to generate the .inp and .mol files (optional)
You can generate input files using the GUI provided by the ioNERDSS library. After installing ioNERDSS using pip, start a Python interpreter and run the following command to start the GUI:
import ioNERDSS as ion ion.gui()
Fill out the parameters in the different sections, add each molecule one by one, followed by each reaction. Finally, click on the "Generate" button. Below are some snapshots.
Run the simulation
Make sure your .inp and .mol files are in the same folder. Navigate to that folder.
Run the simulation in the local environment (Generally faster)
If you want to run NERDSS locally, add it to your PATH and start the simulation by running:
./nerdss -f <your-input-filename>.inp
Run the simulation using Docker
If you are using Docker, start the simulation by running:
docker run -e RUN_NERDSS=true -e ANALYZE_OUTPUT=true -p 8888:8888 -v $(pwd):/SIMULATION -it sikaoguo/nerdsstutorial:latest
The simulation will then begin. The standard output is written to output.log. Once it is done, a Jupyter environment with
ioNERDSS installed will be ready for use.
Outputs of the simulation
The file copy_numbers_time.dat stores the time dependence of the copy numbers of all species in the system. Below are the first five lines of this file.
Time (s),A(a),R(r),A(a!1).R(r!1) 0,1000,1000,0 2e-05,990,990,10 4e-05,980,980,20 6e-05,966,966,34
The file histogram_complexes_time.dat contains the time dependence of the complex components. The first ten lines of this file are shown below.
Time (s): 0 1000 A: 1. 1000 R: 1. Time (s): 2e-05 990 A: 1. 10 A: 1. R: 1. 990 R: 1. Time (s): 4e-05 980 A: 1. 20 A: 1. R: 1. 980 R: 1.
The simulation snapshots are stored in the
PDB/ folder and are in the PDB format.
Analyze the output
We will create a plot of monomer A(a) vs Time (s) and compare the simulation result with the ODE theory.
import pandas as pd import matplotlib.pyplot as plt # Read the simulation data from the file sim_df = pd.read_csv('copy_numbers_time.dat', sep=',') # Read the ODE theory data from the file # Using regex to handle one or more spaces as the separator ode_df = pd.read_csv('theoryODE.dat', sep=r'\s+', skipinitialspace=True) # Plotting plt.figure(figsize=(6, 4)) # Adjust the figure size (width, height) in inches plt.plot(sim_df['Time (s)'], sim_df['A(a)'], 'o-', label='Simulation') plt.plot(ode_df['time(s)'], ode_df['A(t)'], 's-', label='ODE Theory') plt.xlabel('Time (s)', fontsize=14) plt.ylabel('A(a)', fontsize=14) plt.title('Monomer A(a) vs Time (s)', fontsize=16) plt.legend(fontsize=12) plt.grid(True, which='both', linestyle='--', linewidth=0.5) plt.tight_layout() # Save the figure in PNG format with a resolution of 300 dpi plt.savefig('comparison_plot.png', dpi=300, bbox_inches='tight') plt.show()
Render movies for the trajectory
You can generate movies using the PDB format with OVITO.