"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",
"The positions of attached particles are defined relative to the rigid body’s own frame, not the lab frame. To attach a particle, first create the particle and then use `rigid_body_type.attach_particle(particle)`. Alternatively, you can provide a list of particles using the `attached_particles` keyword when creating the `RigidBodyType` object.\n",
"\n",
"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."
