This workflow constructs and simulates a variety of protein-containing systems in explicit solvent: apo proteins, protein-ligand complexes, protein-protein complexes, and systems containing small-molecule cofactors. Upon trajectory completion, Rowan automatically provides users with a thorough analysis of the simulated dynamics.
Protein MD simulations are configurable in both time and number of replicas. By default, Rowan runs 10 ns of production MD per replica with 4 replicas, following a brief minimization, heating, and NPT equilibration.
The post-processing features both analysis per-replicas and pooled across replicas. These include, but are not limited to, binder RMSD, receptor isotropic radius of gyration, hydration site analysis, per-frame SASA, etc.
Solvent is removed from downloadable trajectories to keep file sizes manageable by default. However, users who need the full solvated trajectory can choose to retain solvent.
As a representative benchmark, a single-replica run on a dasatinib•ABL kinase system (PDB 2GQG) with 55,517 atoms solvated, parameterized, minimized, equilibrated, and completed a 10-ns production trajectory in 18 minutes total, reaching speeds of >1 µs/day on our standard L40S GPU hardware.
As a practical starting point for structures containing a ligand, it is often useful to retain heterogens, remove existing hydrogens, and then add hydrogens back in Rowan's protein editing GUI. If heterogens are removed, note that preparation tools may remove ions (including structural ones) in the process. Rowan's protein preparation workflow, whose functionality is described here, is very effective for preparation, too.
Protein MD begins from a prepared protein structure. Small molecules, when present, are parameterized with the protein during OpenMM system construction. Individual small molecules can be designated as a part of the binder; any remaining small molecules are retained as part of the receptor structure. SMILES strings must be provided for all small molecules. Protein binders are identified by chain ID.
By default, Rowan uses the AMBER ff14SB forcefield for proteins, TIP3P for water and ions, and OpenFF Sage 2.3.0 for parameterized small molecules. The system is placed in an explicit-solvent box with an 8 Å water buffer by default. The ionic strength of the solution can be adjusted by users; the default is 0 M added salt. Long-range electrostatics are treated using particle-mesh Ewald, with a default nonbonded cutoff of 8 Å.
Bonds to hydrogens are constrained by default. Rowan also applies hydrogen mass repartitioning, with a default hydrogen mass of 3.0 Da, allowing a 4 fs integration timestep under the default settings. Both the timestep and hydrogen mass are configurable. Rowan generally provides warnings in the GUI when users choose settings that could lead to inaccurate results or unstable simulations.
Protein restraints can optionally be applied in binder-containing systems. When a restraint cutoff is supplied, protein Cα atoms within the cutoff distance of any binder atom are held in place by a harmonic force, with spring constant 100 kcal mol−1 Å−2.
Molecular dynamics simulations are performed with OpenMM using a Langevin-middle integrator. The default temperature is 300 K, with a Langevin friction timescale of 1 ps and a 4 fs timestep. All of these parameters can be changed by users.
After system construction, Rowan first performs energy minimization for up to 5,000 iterations. The minimized structure is saved and can be downloaded after simulation.
The minimized system is then gradually heated from 5 K to the target temperature. Heating proceeds in 5 K increments over roughly 50,000 integration steps. After heating, Rowan adds a Monte Carlo barostat and performs NPT equilibration at the requested pressure and temperature. The defaults are 1 atm, 300 K, and 0.5 ns of equilibration.
The equilibrated positions, velocities, and periodic box vectors are saved as a shared initial state for each replica. The requested number of production replicas are then simulated, with a default production length of 10 ns per replica.
Coordinates are saved every 10 ps by default, though the frame saving interval is also configurable. Potential energy is evaluated for each saved production frame.
By default, solvent is removed from the trajectories that are downloadable by the user. Solvent can instead be retained using the advanced options. For binder-containing systems, users can also choose to retain a specified number of water molecules closest to the binder, providing a compact trajectory that preserves the most relevant local solvent.
Before trajectory analysis, trajectories are imaged to account for periodic boundary conditions and aligned to remove global translation and rotation.
For apo systems and small-molecule-binder systems, trajectory alignment uses all protein Cα atoms. In systems with a protein binder, the trajectory is instead aligned using receptor Cα atoms so that motion of the binder relative to the receptor is preserved.
Some reported observables perform their own alignment independently. In particular, protein RMSD uses all protein Cα atoms, including those belonging to a protein binder.
Rowan reports several general metrics for every replica.
Protein RMSD is calculated for all protein Cα atoms relative to the first trajectory frame. This calculation performs its own Kabsch alignment. So, in systems with protein binders, the reported protein RMSD describes the structural evolution of the complete protein assembly rather than inheriting the receptor-only alignment used for other analyses.
Protein RMSF is calculated for every protein Cα atom relative to the mean trajectory structure, providing a residue-level description of conformational flexibility. The mean trajectory structure is obtained by a coordinate-wise average of each atom's position over the entire trajectory.
Isotropic radius of gyration is evaluated for the receptor polymer chains. Binder chains, small molecules, water molecules, and ions are excluded. This produces a measure of receptor compactness that is not distorted by translation of a binder or by stray heteroatoms.
Potential energy is recorded for each saved frame and returned in Hartree.
When the binder consists of exactly one small molecule, Rowan calculates its heavy-atom RMSD after alignment of the surrounding protein. Because the ligand itself is not used for the fit, this metric captures motion of the ligand relative to the receptor.
When the binder consists of exactly one protein binder chain and no small-molecule binder, Rowan instead computes binder RMSD using the protein's backbone N, Cα, C, and O atoms.
A single binder RMSD is intentionally not reported for multi-component binders (e. g. multiple binder chains, multiple small molecules, combined protein-chain and small-molecule binder, etc.) because collapsing the motion of several independent components into one RMSD would be ambiguous.
Solvent-accessible surface area analysis can optionally be enabled at an inputted sampling interval from the saved frames. Rowan reports both total SASA and polar SASA for the receptor. For these calculations, the receptor is defined as the system excluding solvent, buffer ions, binder chains, and binder-side small molecules. Receptor-side cofactors remain part of the receptor selection.
When a binder is present, the same analysis interval also enables MM/GBSA interaction-energy calculations. Rowan constructs receptor and binder MM/GBSA contexts and evaluates single-point MM/GBSA interaction energies on the selected trajectory frames, returning the results in kcal/mol.
Because SASA and MM/GBSA can be substantially more expensive than basic trajectory metrics, these analyses are not evaluated on every frame unless the requested interval corresponds to every saved trajectory frame.
In addition to frame-by-frame trajectory properties, Rowan provides two structures intended to summarize each production trajectory. The average structure is obtained by averaging atomic coordinates over the trajectory after alignment. Because an averaged structure does not generally correspond to a physically sampled conformation, Rowan also identifies a representative medoid frame taken directly from the trajectory. The medoid is the trajectory frame with the smallest total structural distance to the other frames.
Trajectory clustering can optionally be enabled to identify recurrent conformational states. Rowan supports k-means clustering, requiring a number of clusters (default 10 clusters), and greedy clustering, requiring a conformational similarity threshold in Å (default 2 Å). Rowan returns the cluster assignment for each frame together with the corresponding cluster-centroid frame indices.
For systems containing a binder, Rowan performs hydration-site analysis after all replicas have completed. This analysis operates across all replicas and identifies hydration sites associated with the binder. Protein-water, ligand-water, and protein-water-ligand hydrogen bonding is tracked for each site and returned. Unique water counts per hydration site are also reported. Details on how hydration-site analysis is implemented can be found in this blog.