Example 2. Simulating the Amyloid-\beta42 monomer.

Example 2. Simulating the Amyloid-\(\beta42\) monomer.#

Set simulation#

This example will demonstrate how to create a two-bead coarse-grained model using the default AICG2+ force field of OpenAICG2. We choose the Amyloid $beta 42$ monomer as the model protein. The Amyloid $beta 42$ monomer is an IDP (intrinsically disordered protein), with most of its structure being unordered. However, when its monomer concentration exceeds the critical concentration, Amyloid $beta 42$ will aggregate and deposit, forming relatively ordered structures. Next, we will use OpenAICG2 to run a Langevin dynamics simulation to calculate its radius of gyration.

# add required packages
import openmm as mm
from openmm import app
from openmm import unit
import mdtraj as md
import json
from openmicron.forcefield.aicgmodel import AICG2Model
from openmicron import utils

Load the parameters required for the simulation, including temperature, friction coefficient, time step, total simulation time, total number of steps, trajectory recording interval, output file name, etc.

# load simulation parameters from json file
simu_params_path = '../input/simulationparams.json'
configure_params = json.load(open(simu_params_path,'r'))
T = configure_params['Temperature']
friction = configure_params['friction']
timestep = configure_params['timestep']
tot_simu_steps = configure_params['total_mdsteps']
report_period = configure_params['report_period']
output_file = configure_params['output_file_name']
platform_type = configure_params['platform_type']
initial_pdb_path = configure_params['initial_pdb']
monomer_psf_path = configure_params['monomer_psf']
box = configure_params['box_vector']
native_info_path = configure_params['native_information']
lambdakh_scale = configure_params['lambdakh_scale']
cutoff_kh = configure_params['cutoff_kh']

Initialize the topology from the PDB and PSF files (if the PDB file includes the CONNECT section, you do not need to load the PSF to complete the topology in the PDB).

# load pdb and psf
pdb = md.load(initial_pdb_path)
psf = md.load_psf(monomer_psf_path)
top = pdb.topology.to_openmm()
top._bonds = []
bonds = psf._bonds
# Refine the bond in topology
redefine_top = utils.RedefineTopology()
redefine_top.redefine_bond(top,bonds)

Use the ParserNinfo in the utils to parse the native information file to obtain the force field parameters and store them in the ParserNinfo class variable. # load parameter in native information file

ParserNinfo=utils.ParserNinfo()
ParserNinfo.get_ninfo(native_info_path)

After initialization, use the AICG2Model class to create an AICG model, then create an OpenMM system based on the model. The specific parameters include the initialized top, use_pbc=True to control whether nonbonded interactions are used, and nonbondedMethod to set the properties of nonbonded interactions, including whether they are truncated and whether periodic boundary conditions are used. After creating the system, we use the append_ff_params function to input the parsed force field parameters from ParserNinfo, and then call add_all_default_ener_function to add all the necessary default energy functions. The parameter oriented_Hbond set to False indicates that orientation-dependent hydrogen bonds are not added.

# create model
model = AICG2Model()
model.create_system(top,use_pbc=True,
                 box_a=box['x'], box_b=box['y'], box_c=box['z'],
                 nonbondedMethod=app.CutoffPeriodic,
                 remove_cmmotion=False)
# input native information
model.append_ff_params(ParserNinfo)
# add all default energy function
model.add_all_default_ener_function(oriented_Hbond=False,cutoffkh=cutoff_kh,
                                 kh_epsilon_scale=lambdakh_scale,temperature=T)

Create a Langevin integrator using OpenMM, and then build the simulation.

# create simulation
integrator = mm.LangevinMiddleIntegrator(T*unit.kelvin,friction/unit.picosecond,timestep*unit.femtosecond)
init_coord = pdb.xyz[0,:,:] * unit.nanometer
model.set_simulation(integrator, platform_name=platform_type,properties={'Precision': 'mixed'},init_coord=init_coord)
model.move_COM_to_box_center(use_pbc=False)
model.simulation.context.setVelocitiesToTemperature(T*unit.kelvin)
model.simulation.minimizeEnergy()
# reporter trajectory and log about simulation
model.add_reporters(tot_simu_steps, report_period,
                 output_traj_name='../output/%s'%(output_file),report_traj_format='dcd'
                 ,report_traj=True,report_state_log=True)

print('Running simulation!!!')
try:
   model.simulation.step(tot_simu_steps)
except:
   print('simulation interruption')
   model.save_system('../output/system_%s.xml'%(output_file))
   model.save_state('../output/state_%s.xml'%(output_file))

Compute the radius of gyration \(R_{g}\)#

The standard polymer physics formula for calculating the radius of gyration is as follows:

\[R_g = \sqrt{\frac{1}{N}<\sum_{i=1}^{N}(r_i - r_{CM})>}\]

In the equation, $N$ refers to the number of particles, :math: r_{i} is the coordinate of particle :math: i, :math: r_{CM} is the center of mass coordinate, and <…> denotes ensemble average.

import numpy as np
import matplotlib.pyplot as plt

_A_to_nm = 0.1
lambdascale = 1.13
traj = md.load('../output/monomer.dcd',top='../output/monomer.pdb')
rg = md.compute_rg(traj,masses=None)
Rg = np.mean(rg)/_A_to_nm
print('radius of gyration:',Rg,'angstrom')
Contents