Shake it like a polaroid picture: MD in pyrosetta

> This blog post has been unfinished for two years. So I am posting in the hopes it will spur me to finish it.

The score of a pose reflects how good its interactions are in that static arrangement, a static snapshot. However, given some energy several of these interactions may break and a different conformation is seen. The best way to describe what does 1 kcal/mol mean is that it is the typical strength of a hydrogen bond, but this is rather weak... in fact this is also the average collision energy of water molecules at 37°C, because that is the molar Boltzmann constant times temperature (k_BT/N_A). (At that point in the explanation is it paramount to resist the urge to explain that k_BT coincides with the mean of the Boltzmann distribution describing the energy of collisions as per Maxwell–Boltzmann statistics or else you get that glazed look thermodynamics seems to illicit even in folk that aren't hangover students)

Therefore, hydrogen bonds do come apart and together rather frequently and in some cases these dynamic properties result large scale switching. This cannot really be determined from a static score —even the per residue scores aren't an indication of dynamic properties. So how does one do an MD run in Pyrosetta?

There are several great MD tools, such as NAMD or Gromacs, but if one have set up a topology file and dealt with constraints already in Rosetta one really does not want to switch even if it means using a single CPU. Luckily there is a molecular dynamics mover in Rosetta that means one does not have to —if it is a small task.

Solvent

Rosetta uses an implicit solvent model. This means that water is a magical ghostly thing that is less energetically favourable that crystallographic water and that may result in cavities where a real water molecule would be unhappy for being squeezed. As discussed in a previous post, there is a hydration mover and dedicated scorefunctions to allow a hybrid explicit/implict model. However, this is really not enough. However, on a pinch the MD mover does the job nicely.

Cartesian MD

First I ought to say that there seem to be two MD movers. One called CartesianMD and the other MolecularDynamics. The latter does not work for me and is more undocumented. While the former works. This seems to have been created for a homology refinement paper

scorefxn = pyrosetta.create_score_function('ref2015_cart')
md = pyrosetta.rosetta.protocols.md.CartesianMD(pose, scorefxn)

As you can see it is in cartesian geometry and not dihedral geometry (in cartesian the atoms have an absolute x,y,z position, while in dihedral they have a distance and a some angles releative to their neighbours). This means that it is slower than other movers —we won't even mention speed difference against an Nvidia CUDA accellerated multicore NAMD run.

Anyone that has done an MD run knows there are endless settings. CartesianMD mover has several but is lighter.

• md.set_ncyc_premin(0) this is the burn-in, the steps that are discarded at the beginning
• md.set_ncyc_postmin(0) these are the number of steps discarded at the end
• md.set_nstep(500) how many steps to run. If rattle is one, a step is 2 femptoseconds, so 500 steps is a picosecond. Remember that in Python scientific notation always gives floats, so to give an integer one writes int(1e6), alternatively one can use underscores as a thousands separator so 1_000_000

md.set_ncyc_premin(0)
md.set_temperature(310)
md.use_rattle(True)
md.set_store_trj(True)
md.apply(pose)
#md.get_last_move_status()
len(md.dump_poses(pose)), md.kinetic_energy(), md.potential_energy(), pyrosetta.create_score_function('ref2015_cart')(pose), md.n_dof(), md.last_proposal_density_ratio()

ps = md.dump_poses(pose)
pyrosetta.rosetta.core.scoring.CA_rmsd(ps[1], ps[3])
xxx
MD
This is rubbish??
md = pyrosetta.rosetta.protocols.cartesian.MolecularDynamics(pose, pyrosetta.get_fa_scorefxn())
md.doMinimising(pyrosetta.get_fa_scorefxn())
md.doMD(pyrosetta.get_fa_scorefxn(), 100, 200, 200)
print(md)

Testbed

Barnase/barnstar is the quintessential example for protein-protein interactions, myoglobin of a protein with the most annoying ligand possible, and for MD I would argue that the Trp-cage mini protein is a nice and simple testbed as it's only 20 amino acids long, but the NMR gives multiple states:

original = pyrosetta.toolbox.rcsb.pose_from_rcsb('1L2Y') # this is a mess of 720 residues, i.e. the models have been concatenated.
# usable test bed
pose = pyrosetta.rosetta.protocols.grafting.return_region(original, 1,20) 
# ensemble
poses = pyrosetta.rosetta.utility.vector1_core_pose_Pose()
n = [pyrosetta.rosetta.protocols.grafting.return_region(original, 1+i,20+i) for i in range(0, original.total_residue(),20)]
poses.extend(n)