Protein location in an md simulation

[Bartłomiej Fliszkiewicz]

Hello
I am pretty new to molecular dynamics and wanted to ask a question about protein location. I do a simulation of a protein in water (a little modification of the tutorial - a merge of TeachOpenCADD and the tutorial) and visualize the results. What struck me is that the protein is at the brink of a cube that is created of water molecules. The visualization image:

why.png (176.33 KiB) Viewed 642 times

20.png (150.62 KiB) Viewed 642 times

The code:

Code: Select all

import openmm
import openmm as mm
from openmm import app
from openmm import unit
import mdtraj as md
import pdbfixer

def prepare_protein(
    pdb_file, ignore_missing_residues=True, ignore_terminal_missing_residues=True, ph=7.0
):
    """
    Use pdbfixer to prepare the protein from a PDB file. Hetero atoms such as ligands are
    removed and non-standard residues replaced. Missing atoms to existing residues are added.
    Missing residues are ignored by default, but can be included.

    Parameters
    ----------
    pdb_file: pathlib.Path or str
        PDB file containing the system to simulate.
    ignore_missing_residues: bool, optional
        If missing residues should be ignored or built.
    ignore_terminal_missing_residues: bool, optional
        If missing residues at the beginning and the end of a chain should be ignored or built.
    ph: float, optional
        pH value used to determine protonation state of residues

    Returns
    -------
    fixer: pdbfixer.pdbfixer.PDBFixer
        Prepared protein system.
    """
    fixer = pdbfixer.PDBFixer(str(pdb_file))
    fixer.removeHeterogens()  # co-crystallized ligands are unknown to PDBFixer
    fixer.findMissingResidues()  # identify missing residues, needed for identification of missing atoms

    # if missing terminal residues shall be ignored, remove them from the dictionary
    if ignore_terminal_missing_residues:
        chains = list(fixer.topology.chains())
        keys = fixer.missingResidues.keys()
        for key in list(keys):
            chain = chains[key[0]]
            if key[1] == 0 or key[1] == len(list(chain.residues())):
                del fixer.missingResidues[key]

    # if all missing residues shall be ignored ignored, clear the dictionary
    if ignore_missing_residues:
        fixer.missingResidues = {}

    fixer.findNonstandardResidues()  # find non-standard residue
    fixer.replaceNonstandardResidues()  # replace non-standard residues with standard one
    fixer.findMissingAtoms()  # find missing heavy atoms
    fixer.addMissingAtoms()  # add missing atoms and residues
    fixer.addMissingHydrogens(ph)  # add missing hydrogens
    return fixer

# protein = prepare_protein("inputs/8aen_Achain.pdb", ignore_missing_residues=False)
protein = prepare_protein("inputs/3poz_cleaned.pdb", ignore_missing_residues=False)

forcefield = app.ForceField('amber14-all.xml', 'amber14/tip3pfb.xml')

modeller = app.Modeller(protein.topology, protein.positions)
modeller.addSolvent(forcefield, padding=2.0 * unit.nanometers)
# modeller.addSolvent(forcefield, padding=1.0 * unit.nanometers, ionicStrength=0.15 * unit.molar)

system = forcefield.createSystem(modeller.topology, nonbondedMethod=app.PME)
integrator = mm.LangevinIntegrator(
    300 * unit.kelvin, 1.0 / unit.picoseconds, 2.0 * unit.femtoseconds
)
simulation = app.Simulation(modeller.topology, system, integrator, openmm.Platform.getPlatformByName('CUDA'))
simulation.context.setPositions(modeller.positions)

with open("protein_water_outputs/topology_with_solvent.pdb", "w") as pdb_file:
    app.PDBFile.writeFile(
        simulation.topology,
        simulation.context.getState(getPositions=True).getPositions(),
        # simulation.context.getState(getPositions=True, enforcePeriodicBox=True).getPositions(),
        file=pdb_file,
        keepIds=True,
    )

print('Minimizing...')
simulation.minimizeEnergy()

with open("protein_water_outputs/topology2.pdb", "w") as pdb_file:
    app.PDBFile.writeFile(
        simulation.topology,
        # simulation.context.getState(getPositions=True, enforcePeriodicBox=True).getPositions(),
        simulation.context.getState(getPositions=True).getPositions(),
        file=pdb_file,
        keepIds=True,
    )

simulation.reporters.append(app.PDBReporter('protein_water_outputs/trajectory2.pdb', 1000))

simulation_time = 20 * unit.picoseconds
timestep = 2.0 * unit.femtoseconds

output_interval = 1000
print('Running simulation...')
simulation.context.setVelocitiesToTemperature(300 * unit.kelvin)
simulation.step(int(simulation_time / timestep))
print('Simulation complete!')

Does it somehow affect the simulation? Intuitively I think that it does not since the water box is periodical, but I wanted to ask.

[Peter Eastman]

It doesn't affect the simulation. See this FAQ entry which discusses periodic boundary conditions.