"In addition to point-like particles, ARBD simulations can contain rigid body objects that have an orientation in addition to position. In ARBD, rigid bodies can be decorated with point particles and grids that move with the rigid body. Besides having a mass and translational damping coefficient, rigid bodies have a moment of inertia and rotational damping coefficient.\n",
"\n",
"Below we illustrate how one can use rigid bodies to simulate a large biological complex, namely the trombone loop of the replication fork. A lot of predetermined files are used in this example. For details, we recommend checking out the source files from the iScience manuscript upon which this is based: https://doi.org/10.1016/j.isci.2022.104264\n",
"\n",
"As a first step, we'll simulate just the DNA in trombone loop to help illustrate what it takes to add the rigid bodies to the system."
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "1c2b4cdb",
"metadata": {},
"outputs": [],
"source": [
"import numpy as np\n",
"\n",
"## Steric potential was obtained from a the cryo-EM density of a viral replication fork\n",
"replisome_grid = 'grids/replisome.dx'\n",
"\n",
"## Positions of restraints for DNA start and end determined from inspection of cryo-EM structure\n",
"restraint1 = np.array((189.5, 199.5, 209.8))\n",
"restraint2 = np.array((153.9, 161.1, 196.6)) "
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "b54d488e",
"metadata": {},
"outputs": [],
"source": [
"## Coordinates here are from a simulation. Originally, a pre-equilibrated \n",
"## SSB-coated DNA strand was bent into a \"U\" shape, and we skip this step\n",
"\n",
"dna_coordinates = np.loadtxt('0-initial/p-coords.dat') # DNA coordinates\n",
"ssb_coords = np.loadtxt('0-initial/ssb.rigid', usecols=list(range(2,2+3+9))) ## SSB coordinates and orientations"
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "6818145a",
"metadata": {},
"outputs": [],
"source": [
"## Create a function to construct the DNA model (we'll repeat this a few times)\n",
"## Optionally simulate the DNA-only model to see how it performs and behaves\n",
"model = make_dna_model()\n",
"model.simulate(output_name='run',\n",
" directory='trombone-no_ssb',\n",
" num_steps=1e5, output_period=1e3,\n",
")"
]
},
{
"cell_type": "markdown",
"id": "7438e38d",
"metadata": {},
"source": [
"You can load this trajectory in VMD to inspect as usual. You should also load `grids/replisome.dx` or `replisome.{psf,pdb}` so you can see the DNA in the context of the replication machinery. Also, please note that the ssDNA model used here represents poly(dT), so it is unable to base pair or form secondary structure, which is the main reason why cells need SSB proteins."
]
},
{
"cell_type": "markdown",
"id": "6d467e3f",
"metadata": {},
"source": [
"## Adding single-stranded DNA binding protein\n",
"\n",
"SSB proteins are required for replicating cells. Their main function is to prevent secondary structure formation in ssDNA exposed during DNA replication, but they also may act as signaling molecules through their disordered C-terminal tails, which have been truncated here. \n",
"\n",
"They prevent secondary structure formation by binding ssDNA with very high affinity (k_off between SSB and dT_70 is so long that it is difficult to measure experimentally; likely days). But this high affinity is apparently at odds with the constantly moving nature of the replication fork (or replisome). The DNA in the lagging strand \"trombone loop\" where most of the DNA-bound SSB is found is in constant flux.\n",
"\n",
"We previously used the IBI method to construct grid-based potentials for our two-site DNA model interacting with SSB (grids/grid-P.dx and grids/grid-B.dx) using 3D target distributions extracted from all-atom simulations of dT_5 fragments binding to SSB (note the arbdmodel.ibi module does not yet provide support for 3D dsitributions; if you need this feature, please reach out to cmaffeo2@illinois.edu).\n",
"\n",
"We also clustered the atoms in the SSBs according to their Lennard-Jones parameters to allow the atomistic interactions between two proteins to be represented by a tractable number of grids containing pre-computed potentials and the sites they act on. In addition, we computed the electrostatic potential around SSB using the Adaptive Poisson Boltzmann Solver. These operations are to be automated in a future release of arbdmodel."
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "04844c8a",
"metadata": {},
"outputs": [],
"source": [
"## Prepare the interactions for the SSBs\n",
"\n",
"kT = 0.578277 # in arbd energy units of kcal/mol\n",
"\n",
"## Each potential is specified with a tuple consisting of:\n",
"## a keyword, a path to a grid file, and an optional scaling factor\n",
"_potentials = [(f'{key}bead',f'grids/grid-{key}.dx', kT) for key in ('P','B')] # IBI potential for P/B beads\n",
"\n",
"## Potentials will act on any point particle with a matching\n",
"## rigid_body_potentials keyword, so let's add those now\n",
"_P.rigid_body_potentials = ['Pbead'] # note that there can be multiple entries\n",
"_B.rigid_body_potentials = ['Bbead']\n",
"\n",
"## Now let's consider SSB-SSB interactions\n",
"_potentials += [(key,f'grids/1eyg.{key}.pot.dx') for key in ('vdw0','vdw1','vdw2')] # approximate SSB-SSB vdW interactions \n",
" moment_of_inertia = [19375322.0, 18923208.0, 13511100.0], # from VMD\n",
" damping_coefficient = [779.075531386, 757.10491996, 735.261090254], # from hyropro\n",
" rotational_damping_coefficient = [3096.68320292, 3027.72621629, 3547.0215724], # from hydroproo\n",
" potential_grids = _potentials,\n",
" charge_grids = _generalized_charges,\n",
" pmf_grids = _pmf_grids,\n",
" )\n",
"\n",
"## Now create SSB rigid body objects\n",
"ssbs = []\n",
"for i,coord in enumerate(ssb_coords):\n",
" rb = RigidBody( ssb_type, name=f'SSB{i}',\n",
" position = coord[:3],\n",
" orientation = coord[3:].reshape((3,3)) )\n",
" ssbs.append(rb)\n",
" \n",
"## Create new DNA model because prior simulate command messes up cache of rigid body objects\n",
"model = make_dna_model()\n",
"\n",
"model.add( Group( 'SSBs', ssbs ) ) # Add the rigid bodies to the model\n",
"\n",
"# model._updateParticleOrder() # We need to call this\n",
"\n",
"model.simulate(output_name='run',\n",
" directory='trombone-ssb',\n",
" num_steps=1e5, output_period=1e3,\n",
" rigid_body_grid_grid_period = 20, # dynamics of rigid bodies is much slower than DNA, so we can recompute grid-grid force and torque less frequently\n",
" gpu=1\n",
")"
]
},
{
"cell_type": "markdown",
"id": "0884cf4b",
"metadata": {},
"source": [
"Now that you've simulated (or are simulating), have a peak at the `trombonde-ssb/output/` directory. Notice that in addition to a .dcd file containing the particle coordinates, there is a (plain text) .rb-traj file that contains the rigid body coordinates and orientations. Unfortunately for us, VMD does not yet enjoy native support for rigid body objects. Our workaround is to load the rigid bodies as single-frame all-atom molecules and take advantage of the `trace` feature of Tcl to update the positions and orientations of those molecules when the dcd frame for the particles changes in VMD. The result is a little slow, but it allows us to visualize and analyze the results of such trajectories using familiar tools.\n",
"\n",
"From a terinal, run `vmd -e load.tcl` or, from the TkConsole of VMD, run `source load.tcl`. Inspect the trajectory, and close VMD when you are done."
]
},
{
"cell_type": "markdown",
"id": "e21e8d86",
"metadata": {},
"source": [
"## Attaching point particles\n",
"\n",
"Besides attaching grids to rigid body objects, we can attach point particles that can be subject to all the usual point-particle potentials. We can use this feature, for example, to confine the center of a rigid body, or to connect multiple rigid bodies together, or to add disordered tails to a structured protein. Note that rigid-body attached particles are not written to the particle dcd and their mass/friction is ignored.\n",
"\n",
"The syntax for attaching a particle is simple: First create a particle and then run `rigid_body_type.attach_particle( particle )`, or you can provide a list of particles with the keyword `attached_particles` when you create the RigidBodyType object. One caveat: the attached particles cannot belong to any group, and you should wait to define any bonded interactions or restraints, since these act on particle index, which will be unique for the (copied) particle later attached to each rigid body. Non-bonded interactions and those that involve grids or rigid body potentials can be set up beforehand. Once you've constructed rigid body objects, each one will have the attribute `rigid_body.attached_particles` which you can use to access the particles for use in bonded interactions.\n",
"\n",
"As an exercise, see if you can apply a constant force to the SSBs using attached particles subject to a grid-based force."
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "9ff1e22f",
"metadata": {},
"outputs": [],
"source": [
"## code goes here :)"
]
}
],
"metadata": {
"kernelspec": {
"display_name": "Python 3",
"language": "python",
"name": "python3"
},
"language_info": {
"codemirror_mode": {
"name": "ipython",
"version": 3
},
"file_extension": ".py",
"mimetype": "text/x-python",
"name": "python",
"nbconvert_exporter": "python",
"pygments_lexer": "ipython3",
"version": "3.8.18"
}
},
"nbformat": 4,
"nbformat_minor": 5
}
%% Cell type:markdown id:6004c298 tags:
# Modeling rigid-body objects
In addition to point-like particles, ARBD simulations can contain rigid body objects that have an orientation in addition to position. In ARBD, rigid bodies can be decorated with point particles and grids that move with the rigid body. Besides having a mass and translational damping coefficient, rigid bodies have a moment of inertia and rotational damping coefficient.
Below we illustrate how one can use rigid bodies to simulate a large biological complex, namely the trombone loop of the replication fork. A lot of predetermined files are used in this example. For details, we recommend checking out the source files from the iScience manuscript upon which this is based: https://doi.org/10.1016/j.isci.2022.104264
As a first step, we'll simulate just the DNA in trombone loop to help illustrate what it takes to add the rigid bodies to the system.
%% Cell type:code id:1c2b4cdb tags:
``` python
importnumpyasnp
## Steric potential was obtained from a the cryo-EM density of a viral replication fork
replisome_grid='grids/replisome.dx'
## Positions of restraints for DNA start and end determined from inspection of cryo-EM structure
restraint1=np.array((189.5,199.5,209.8))
restraint2=np.array((153.9,161.1,196.6))
```
%% Cell type:code id:b54d488e tags:
``` python
## Coordinates here are from a simulation. Originally, a pre-equilibrated
## SSB-coated DNA strand was bent into a "U" shape, and we skip this step
dna_coordinates=np.loadtxt('0-initial/p-coords.dat')# DNA coordinates
ssb_coords=np.loadtxt('0-initial/ssb.rigid',usecols=list(range(2,2+3+9)))## SSB coordinates and orientations
```
%% Cell type:code id:6818145a tags:
``` python
## Create a function to construct the DNA model (we'll repeat this a few times)
fromarbdmodel.polymerimportPolymerSection
fromarbdmodel.ssdna_two_beadimportDnaModel,_P,_B
defmake_dna_model():
## Create spline-based polymer
a=PolymerSection(name='loop',num_monomers=500,# 500 nt fragment
## Our DNA model is conveniently already implemented. Here we add the
## steric replisome potential (based on cryo-EM density) to both
## bead types. Note this model represents poly(dT), and is not ideal
## for study of the replisome
_P.grid=(replisome_grid,500)# scale it up by a big factor
_B.grid=(replisome_grid,500)
model=DnaModel([a],dimensions=(3000,3000,3000),
timestep=20e-6
)
## Restrain DNA ends
end1=model.strands[0].children[0].children[0]
end2=model.strands[0].children[-1].children[0]
end1.add_restraint((10,restraint1))
end2.add_restraint((10,restraint2))
returnmodel
```
%% Cell type:code id:d74b1f0c tags:
``` python
## Optionally simulate the DNA-only model to see how it performs and behaves
model=make_dna_model()
model.simulate(output_name='run',
directory='trombone-no_ssb',
num_steps=1e5,output_period=1e3,
)
```
%% Cell type:markdown id:7438e38d tags:
You can load this trajectory in VMD to inspect as usual. You should also load `grids/replisome.dx` or `replisome.{psf,pdb}` so you can see the DNA in the context of the replication machinery. Also, please note that the ssDNA model used here represents poly(dT), so it is unable to base pair or form secondary structure, which is the main reason why cells need SSB proteins.
%% Cell type:markdown id:6d467e3f tags:
## Adding single-stranded DNA binding protein
SSB proteins are required for replicating cells. Their main function is to prevent secondary structure formation in ssDNA exposed during DNA replication, but they also may act as signaling molecules through their disordered C-terminal tails, which have been truncated here.
They prevent secondary structure formation by binding ssDNA with very high affinity (k_off between SSB and dT_70 is so long that it is difficult to measure experimentally; likely days). But this high affinity is apparently at odds with the constantly moving nature of the replication fork (or replisome). The DNA in the lagging strand "trombone loop" where most of the DNA-bound SSB is found is in constant flux.
We previously used the IBI method to construct grid-based potentials for our two-site DNA model interacting with SSB (grids/grid-P.dx and grids/grid-B.dx) using 3D target distributions extracted from all-atom simulations of dT_5 fragments binding to SSB (note the arbdmodel.ibi module does not yet provide support for 3D dsitributions; if you need this feature, please reach out to cmaffeo2@illinois.edu).
We also clustered the atoms in the SSBs according to their Lennard-Jones parameters to allow the atomistic interactions between two proteins to be represented by a tractable number of grids containing pre-computed potentials and the sites they act on. In addition, we computed the electrostatic potential around SSB using the Adaptive Poisson Boltzmann Solver. These operations are to be automated in a future release of arbdmodel.
%% Cell type:code id:04844c8a tags:
``` python
## Prepare the interactions for the SSBs
kT=0.578277# in arbd energy units of kcal/mol
## Each potential is specified with a tuple consisting of:
## a keyword, a path to a grid file, and an optional scaling factor
_potentials=[(f'{key}bead',f'grids/grid-{key}.dx',kT)forkeyin('P','B')]# IBI potential for P/B beads
## Potentials will act on any point particle with a matching
## rigid_body_potentials keyword, so let's add those now
_P.rigid_body_potentials=['Pbead']# note that there can be multiple entries
## We can make an SSB feel the potential due to another rigid body by
## giving it a matching generalized "charge" grid. This works by
## finding where each charge voxel is located within a corresponding
## potential, obtaining the gradient of the potential, and applying
## a force (and torque) to the rigid body with the density that is
## proportional to the value of the density voxel, along with
## third-law forces and torques to the rigid body with the
## potential. Please note that the "charge" here is not
## necessarily electrostatic and produces no potential of its own.
_generalized_charges=[(key,f'grids/1eyg.{key}.gencharge.dx')forkeyin('vdw0','vdw1','vdw2')]# atoms in each voxel for steric interactions
_generalized_charges+=[('elec',f'grids/1eyg.charge.dx')]# this one is actually meant to represent the electrostatic charge
## Finally, we don't want our SSB to crash into the replication
## proteins, so we add a "PMF" grid that is fixed in space and acts
## on any generalized charge grid that has a corresponding keyword
_pmf_grids=[('vdw0',replisome_grid,10)]# scale it up a bit
```
%% Cell type:code id:0a1b213b tags:
``` python
## Now let's add SSB using predetermined grids
fromarbdmodelimportRigidBodyType,RigidBody,Group
ssb_type=RigidBodyType('SSB',
mass=48856.41410648823,
moment_of_inertia=[19375322.0,18923208.0,13511100.0],# from VMD
damping_coefficient=[779.075531386,757.10491996,735.261090254],# from hyropro
rotational_damping_coefficient=[3096.68320292,3027.72621629,3547.0215724],# from hydroproo
potential_grids=_potentials,
charge_grids=_generalized_charges,
pmf_grids=_pmf_grids,
)
## Now create SSB rigid body objects
ssbs=[]
fori,coordinenumerate(ssb_coords):
rb=RigidBody(ssb_type,name=f'SSB{i}',
position=coord[:3],
orientation=coord[3:].reshape((3,3)))
ssbs.append(rb)
## Create new DNA model because prior simulate command messes up cache of rigid body objects
model=make_dna_model()
model.add(Group('SSBs',ssbs))# Add the rigid bodies to the model
# model._updateParticleOrder() # We need to call this
model.simulate(output_name='run',
directory='trombone-ssb',
num_steps=1e5,output_period=1e3,
rigid_body_grid_grid_period=20,# dynamics of rigid bodies is much slower than DNA, so we can recompute grid-grid force and torque less frequently
gpu=1
)
```
%% Cell type:markdown id:0884cf4b tags:
Now that you've simulated (or are simulating), have a peak at the `trombonde-ssb/output/` directory. Notice that in addition to a .dcd file containing the particle coordinates, there is a (plain text) .rb-traj file that contains the rigid body coordinates and orientations. Unfortunately for us, VMD does not yet enjoy native support for rigid body objects. Our workaround is to load the rigid bodies as single-frame all-atom molecules and take advantage of the `trace` feature of Tcl to update the positions and orientations of those molecules when the dcd frame for the particles changes in VMD. The result is a little slow, but it allows us to visualize and analyze the results of such trajectories using familiar tools.
From a terinal, run `vmd -e load.tcl` or, from the TkConsole of VMD, run `source load.tcl`. Inspect the trajectory, and close VMD when you are done.
%% Cell type:markdown id:e21e8d86 tags:
## Attaching point particles
Besides attaching grids to rigid body objects, we can attach point particles that can be subject to all the usual point-particle potentials. We can use this feature, for example, to confine the center of a rigid body, or to connect multiple rigid bodies together, or to add disordered tails to a structured protein. Note that rigid-body attached particles are not written to the particle dcd and their mass/friction is ignored.
The syntax for attaching a particle is simple: First create a particle and then run `rigid_body_type.attach_particle( particle )`, or you can provide a list of particles with the keyword `attached_particles` when you create the RigidBodyType object. One caveat: the attached particles cannot belong to any group, and you should wait to define any bonded interactions or restraints, since these act on particle index, which will be unique for the (copied) particle later attached to each rigid body. Non-bonded interactions and those that involve grids or rigid body potentials can be set up beforehand. Once you've constructed rigid body objects, each one will have the attribute `rigid_body.attached_particles` which you can use to access the particles for use in bonded interactions.
As an exercise, see if you can apply a constant force to the SSBs using attached particles subject to a grid-based force.
