Brownian Dynamics in IMP

Table of Contents

• Setup
• Representing parts
• Representing interactions
• Representing dynamics
• Running the script
• Visualization
• Next steps

In this tutorial, we show how to setup and run a Brownian Dynamics simulation in IMP.

This tutorial can be followed in several ways:

• Download the files using the "Clone or download" link at the tutorial's GitHub page and use them in conjunction with this text.
• Download the files from GitHub and, using Jupyter Notebook, open the notebook scripts/pbc_simple.ipynb.
• Load the tutorial directly in your browser, courtesy of Google Colaboratory. (This needs no software installed on your machine, but may take a while to load.)

Setup

We need to import IMP modules that we plan to use in the script:

from __future__ import print_function, division
import IMP.atom
import IMP.algebra
import IMP.rmf
import IMP.core
import RMF
import IMP.container
import IMP.display
import GranuleFactory

Next, we can set parameters for the simulation:

RMF_FILENAME = "pbc_simple.rmf"
# I. Parts parameters:
L = 50000 # Length of our bounding box, A
R = 20000 # PBC radius, A
R_NUCLEUS = 10000 # NE radius, A
N_GRANULES = 50 # Number of granules
R_GRANULES = 1500 # Radius of granules, A
N_GRANULE_PATCHES = 6

# II. Interaction parameters:
K_BB = 0.1  # Strength of the harmonic boundary box in kcal/mol/A^2
K_EXCLUDED=0.1 # Strength of lower-harmonic excluded volume score in kcal/mol/A^2

# III. Time parameters:
BD_STEP_SIZE_SEC= 10E-8
SIM_TIME_SEC= 0.050
bd_step_size_fs= BD_STEP_SIZE_SEC * 1E+15
sim_time_ns= SIM_TIME_SEC * 1E+9
RMF_DUMP_INTERVAL_NS= sim_time_ns / 1000.0

Representing parts

To represent parts of the system, we use a number of Python classes provided by IMP:

• IMP.Model A container for all of the system’s parts
• IMP.Particle A particle is the elementary data unit describing a system part in IMP
• IMP.Decorator Dynamic descriptors of the properties of IMP particles
  • IMP.core.XYZR
  • IMP.atom.Hierarchy
  • IMP.atom.Diffusion
• IMP.container Static or dynamic collections of particles, particle pairs, etc.

To see these classes used in a simulation, look at the file pbc_simple.py in the scripts directory. First, we create an IMP Model, and set up the root of a hierarchy of particles:

# Model:
m = IMP.Model()
# Root of parts hierarchy:
p_root= IMP.Particle(m, "root")
h_root = IMP.atom.Hierarchy.setup_particle(p_root)

Next, we can set up a coarse-grained nucleus by calling the create_nucleus function. The nucleus particle is given XYZ coordinates and a radius using the XYZR decorator. We then add it to the existing hierarchy:

def create_nucleus(m, R):
    '''
    Generate a coarse-grained spherical nuclear envelope
    of radius R in model m
    '''
    p = IMP.Particle(m, "nucleus")
    xyzr = IMP.core.XYZR.setup_particle(p)
    xyzr.set_coordinates_are_optimized(True)
    xyzr.set_coordinates([0,0,0])
    xyzr.set_radius(R)
    IMP.atom.Mass.setup_particle(p, 1.0) # fake mass
    IMP.display.Colored.setup_particle(p,
                                       IMP.display.get_display_color(2))
    IMP.atom.Hierarchy.setup_particle(p)
    return p

# Nucleus:
p_nucleus = create_nucleus(m, R_NUCLEUS)
h_nucleus = IMP.atom.Hierarchy(p_nucleus)
h_root.add_child(h_nucleus)

Next, we create a number of insulin granules at random locations in the cytoplasm. This is done using a simple piece of code, that creates similar XYZR particles, in GranuleFactory.py:

# Granules hierarchy root:
p_granules_root= IMP.Particle(m, "Granules")
IMP.atom.Mass.setup_particle(p_granules_root, 1.0) # fake mass
h_granules_root= IMP.atom.Hierarchy.setup_particle(p_granules_root)
h_root.add_child(h_granules_root)

# PBC cytoplasm bounding sphere:
pbc_sphere= IMP.algebra.Sphere3D([0,0,0], R)
# Actual granules:
nucleus_sphere= IMP.core.XYZR(p_nucleus).get_sphere()
gf=GranuleFactory.GranuleFactory(
    model=m, default_R=R_GRANULES,
    cell_sphere=pbc_sphere,
    nucleus_sphere=nucleus_sphere)
for i in range(N_GRANULES):
    granule= gf.create_simple_granule("Granule_{}".format(i))
    h_granule= IMP.atom.Hierarchy(granule)
    h_granules_root.add_child(h_granule)

Representing interactions

Next, we need to tell IMP the nature of the interactions in the system. These are handled by several IMP classes:

• Scores evaluate some property of the system, for example how far a particle is outside of a sphere.
• Restraints apply Scores to a subset of the system.
• Containers describe the subsets that the Restraints use.

# ----- II. System interactions: -----

# Outer bounding box for simulation:
bb = IMP.algebra.BoundingBox3D(IMP.algebra.Vector3D(-L/2, -L/2, -L/2),
                               IMP.algebra.Vector3D(L/2, L/2, L/2))

# Add enclosing spheres for pbc and outer simulation box
bb_harmonic= IMP.core.HarmonicUpperBound(0, K_BB)
pbc_bsss = IMP.core.BoundingSphere3DSingletonScore(bb_harmonic,
                                                   pbc_sphere)
outer_bbss = IMP.core.BoundingBox3DSingletonScore(bb_harmonic,
                                                  bb)
# Restraints - match score with particles:
rs = []
rs.append(IMP.container.SingletonsRestraint(pbc_bsss,
                                            h_granules_root.get_children()))
rs.append(IMP.container.SingletonsRestraint(outer_bbss,
                                            h_granules_root.get_children()))
# Add excluded volume restraints among all (close pairs of) particles:
ev = IMP.core.ExcludedVolumeRestraint(IMP.atom.get_leaves(h_root),
                                      K_EXCLUDED,
                                      10, # slack affects speed only
                                          # (slack of close pairs finder)
                                      "EV")
rs.append(ev)
# Scoring Function from restraints
sf = IMP.core.RestraintsScoringFunction(rs, "SF")

Representing dynamics

Finally we can tell IMP to run Brownian Dynamics. This uses the IMP.atom.BrownianDynamics class. First, we set properties of the simulation, such as the temperature and number of time steps:

# -------- III. System dynamics: --------

bd = IMP.atom.BrownianDynamics(m)
bd.set_log_level(IMP.SILENT)
bd.set_scoring_function(sf)
bd.set_maximum_time_step(bd_step_size_fs) # in femtoseconds
bd.set_temperature(300)

Next, we request output of the trajectory in RMF format (using the add_optimizer_state method):

# -------- Add RMF visualization --------
def convert_time_ns_to_frames(time_ns, step_size_fs):
    '''
    Given time in nanoseconds time_ns and step size in femtosecond
    step_size_fs, return an integer number of frames greater or equal
    to 1, such that time_ns*step_size_fs is as close as possible to
    time_ns.
    '''
    FS_PER_NS= 1E6
    time_fs= time_ns * FS_PER_NS
    n_frames_float= (time_fs+0.0) / step_size_fs
    n_frames= int(round(n_frames_float))
    return max(n_frames, 1)
sim_time_frames= convert_time_ns_to_frames(sim_time_ns, bd_step_size_fs)
rmf_dump_interval_frames= convert_time_ns_to_frames(RMF_DUMP_INTERVAL_NS, bd_step_size_fs)

print("Simulation time {:.1e} ns / {} frames; "
      "RMF dump interval {:.1e} ns / {} frames".format(sim_time_ns,
                                                      sim_time_frames,
                                                       RMF_DUMP_INTERVAL_NS,
                                                      rmf_dump_interval_frames))

rmf = RMF.create_rmf_file(RMF_FILENAME)
rmf.set_description("Brownian dynamics trajectory with {}fs timestep.\n"\
                    .format(bd_step_size_fs))
IMP.rmf.add_hierarchy(rmf, h_root)
IMP.rmf.add_restraints(rmf, rs)
IMP.rmf.add_geometry(rmf, IMP.display.BoundingBoxGeometry(bb))
IMP.rmf.add_geometry(rmf, IMP.display.SphereGeometry(pbc_sphere))
# Pair RMF with model using an OptimizerState ("listener")
sos = IMP.rmf.SaveOptimizerState(m, rmf)
sos.set_log_level(IMP.SILENT)
sos.set_simulator(bd)
sos.set_period(rmf_dump_interval_frames)
bd.add_optimizer_state(sos)
# Dump initial frame to RMF
sos.update_always("initial conformation")

Finally, we run the simulation by calling optimize (this takes a minute or two to run):

# -------- Run simulation ---------
print("Running simulation")
#m.update()
print("Score before: {:f}".format(sf.evaluate(True)))
bd.optimize(sim_time_frames)
print("Run finished successfully")
print("Score after: {:f}".format(sf.evaluate(True)))

Running the script

The script can be run with Python by simply running:

python pbc_simple.py

Visualization

The final RMF trajectory pbc_simple.rmf can be viewed in UCSF Chimera, or in UCSF ChimeraX with the RMF plugin (installed from the Tools menu under 'More Tools').

Next steps

This demonstration uses a very simple representation of the insulin granules. See also pbc_patches.py and pbc_interacting_patches.py in the scripts directory of the GitHub repository for similar scripts that use a more realistic representation.