Summary
In this experiment, the model is a very important part of our experiment. Due to the consideration of biological safety, we can't use Ureb of Helicobacter pylori in practical experiments, so we tested the binding effect of nano-antibody and Ureb of Helicobacter pylori by molecular docking and other means. Through computer simulation, we deeply analyzed the potential inhibitory effect of our nano-antibody on Helicobacter pylori and established a good model foundation for subsequent experiments.
Molecular docking
Goal:
The binding effect of
antigen and nano-antibody was
simulated by
molecular docking method, and then the nano-antibody with the
strongest binding effect on
Helicobacter
pylori was screened out.
Experimental method:
The nucleic acid sequences
of three nano-antibodies (Nb-6,
Nb-human,
Nb-ScFv) were obtained by literature search, and then the
nucleic acid sequences were translated
into
protein sequences by Snapgene. Because there is no suitable
analytical structure for molecular
docking at
present, ProteinX is used to simulate its structure, and Rosetta
Relax is used to optimize its
energy, and
the pdb format file is output for subsequent
calculation.
Firstly, HDOCK was used for
rigid docking, and the initial
conformation
of the two proteins was adjusted. Then, based on the rigid
docking results, RosettaDock was used
for
flexible docking, and the selection of the results depended on
the score module built in
Rosetta. The
bonding between the two proteins was analyzed by Ligplot
software, and the interaction binding
energy of the
two proteins was analyzed by the Interface_analyzer module of
FoldX. The conformational display
of protein
after docking was carried out by Pymol software.
Result:
The 3D position of Ureb
protein after interaction with three
antibodies
is shown in the figure. Among them, the green protein is UreB
protein, the blue protein is Nb6
protein
(Figure1A), the pink protein is Nb-Human protein (Figure1B) and
the purple protein is Nb-ScFv
protein
(Figure1C), and the bonding situation among the residues has
been marked in the figure, in which
the red
dotted line represents hydrogen bond interaction, the yellow
dotted line represents salt bridge
interaction,
and the number near the dotted line of the interaction bond is
the bond length of the
interaction bond in
Amy. Through the comprehensive analysis of these interactions,
we can screen out the
nano-antibody protein
with the strongest binding ability.
Figure 1
Interaction between Ureb protein
and
nano-antibody protein. A. Schematic diagram of the
interaction between UreB protein and Nb-6
protein, in
which UreB protein is shown in green and NB6 protein is
shown in blue. B. Schematic diagram
of interaction
between UreB protein and Nb-Human protein, in which UreB
protein is shown in green and
NB-human protein is
shown in pink. C. Schematic diagram of interaction
between UreB protein and Nb-ScFv protein,
in which UreB
protein is shown in green and Nb-ScFv protein is shown
in purple.
Using the Interface analyzer
module of FoldX, the binding
energies of
UreB protein with three proteins are shown in the following
table(Table 1), among which UreB
protein has a
strong binding tendency with Nb-SvFc protein and Nb-Human
protein (protein-protein interaction
binding
energy is less than -25 kcal/mol, which means that there is a
strong binding tendency). The
interaction of
three pairs of proteins can be further observed in combination
with subsequent molecular
dynamics
simulation.
Table 1 The
binding energies of UreB protein
with three
proteins
Molecule dynamics
Goal:
The binding ability of
antigen and nano-antibody was further
explored by
molecular dynamics simulation.
Experimental method:
Based on the above docking,
the protein complexes obtained
as the
initial structure were simulated by all-atom molecular dynamics,
and the simulation was carried
out by
Gromacs 2023.3 software. The protein was described by AMBER14SB
protein force field[1,2]. The
pdb2gmx
subroutine is used to add hydrogen atoms to the system, the
truncated cubic TIP3P solvent box is
added at
the distance of the system 10[3], and Na+/Cl- is added to the
system to balance the system
charges. Finally,
the topology and parameter files for simulation are
output.
Molecular dynamics
simulation is carried out by Gromacs
2023.3
software[4], and the simulation time is 100ns. Before the
simulation, the energy of the system
was
optimized, and the energy was minimized (canonical ensemble) by
mdrun command and steepest
descent method.
The initial step size was set to 0.01 nm, and the maximum
allowable force was 1000 kJ /mol•nm.
After the
energy optimization of the system is completed, the NVT
(Isothermal Isomer) family simulation of
the system
is carried out at a constant volume and heating rate of 100 ps,
so that the temperature of the
system slowly
rises from 0 K to 310.15 K, and the solvent molecules are
further evenly distributed in the
solvent box.
Then, 100 ps NPT ensemble simulation was carried out by
Berendsen constant pressure machine, and
the
pressure of solvent and complex system was balanced, so that the
system pressure rose to 1 bar.
In MD
simulation, all hydrogen bonds involved are constrained by LINCS
algorithm, and the integration
step is 2
fs. The electrostatic interaction is calculated by
(particle-mesh Ewald) PME method[5], and the
cutoff value
is set to 1.2 nm. The non-key interaction truncation value is
set to 10, which is updated every
10 steps.
The simulated trajectory is periodically eliminated, and the
subsequent analysis of RMSD, RMSF,
Rg and the
number of hydrogen bonds is carried out.
Calculation of binding free energy of
MMGBSA
The binding free energy
between proteins was calculated by
MM/GBSA
method using gmx_MMPBSA software[6-8]. The specific formula is
as follows:
In the above formula，stands
for internal energy、stands for
van der Waals
action andstands for electrostatic interaction. The internal
energy includes bond
energy（Ebond）、angular
energy（Eangle）and torsion energy（Etorsion）；and are collectively
referred to as solvation free
energy.Among
them, GGBis polar solvation free energy、GSAGSA is nonpolar
solvation free energy. For，this paper
uses the GB
model developed by Nguyen et al [9 ]to calculate（igb = 2）.The
nonpolar solvation free energy（GSA
is
calculated based on the product of surface tension (γ) and
solvent accessible surface area
(SA)，GSA= 0.0072
× SASA[10]。Entropy change is neglected in this study because of
its high computational resources
and low
accuracy[6].
Result:
RMSD (Root Mean Square Deviation)
Figure 2
RMSD Curve of Ureb and Three
Proteins in
Simulated Docking Process.
RMSD curve shows the changes
of root mean square
displacement of three
proteins relative to UreB protein in the simulation process and
is used to evaluate the
stability of binding
between proteins in this study. From the graph analysis, it can
be seen that the RMSD curve of
Nb6 protein
relative to UreB protein shows a certain fluctuation trend in
the first 80ns of the simulation.
After 80ns,
the RMSD curve gradually converges and vibrates around 0.8nm
until the end of the simulation.
The RMSD curve
of Nb-ScFv protein relative to UreB protein experienced a
certain fluctuation trend in the first
60ns of the
simulation, and gradually converged after 60ns, vibrating around
0.5nm until the end of the
simulation. The
RMSD curve of Nb-Human protein relative to UreB protein showed a
rapid increase in the first
20ns of the
simulation, and then gradually converged, vibrating around 1.2nm
until the end of the
simulation.
Nb-6 protein is composed of
more flexible random curls at
the
interaction interface in the first 80ns of simulation, and the
interface is pulled by UreB
protein in the
process of dynamic simulation, so the protein interaction
interface maintains a certain
conformational
adjustment in the first 80 ns of simulation, which leads to a
certain fluctuation trend of RMSD
curve. After
80 ns, with the interaction between proteins gradually
stabilized, and there was no obvious
conformational
change again, the skeleton structure of the two proteins
gradually stabilized, only maintaining
a certain
change in the mutual orientation of side chains. Therefore,
after 80 ns, the RMSD curve
gradually converged
and vibrated around 0.8 nm until the end of the simulation. It
can be considered that Nb6
protein and UreB
protein formed a stable complex after 100ns molecular dynamics
simulation.
Figure 3
Schematic Diagram of Conformation
Rendering of
UreB-Nb6 Protein Complex at 0 or 100ns. Among them,
green shows UreB protein, and blue shows
Nb6
protein.
Nb-ScFv protein changed
little relative to UreB protein in
the whole
simulation process, so the RMSD value was low in the whole
kinetic simulation process. In the
first 60 ns of
the simulation, the RMSD curve kept a certain fluctuation trend,
which was due to the fact that
although the
skeleton atoms of the two proteins did not change obviously at
the interaction interface, the
orientation of
the side chains changed to some extent, so the RMSD curve still
showed a certain fluctuation
trend. After
60ns, as the interaction mode between the two proteins tends to
be stable, the RMSD curve also
tends to
converge and vibrates around 0.5nm until the end of the
simulation.
Figure 4
Schematic Diagram of Conformation
Rendering of
UreB-Nb-ScFv Protein Complex at 0 or 100ns. Among them,
green shows UreB protein, and purple
shows Nb-ScFv
protein.
In the first 20 ns of the
simulation, Nb-Human protein
appears to be
close to UreB protein, which is most obviously marked in the
above figure. In the first 20ns of
the
simulation, Nb-Human protein tends to twist and fold
spontaneously under the action of UreB
protein, and
this conformation is maintained until the end of the simulation.
After 20ns, the interaction
between the two
proteins tends to be stable, and the skeleton atoms do not
change obviously again, so the RMSD
curve tends
to converge gradually. During this period, the vibration of the
RMSD curve still comes from the
mutual
orientation change of the side chains of the two proteins, and
the whole two proteins tend to
form stable
complexes.
Figure 5
Schematic Diagram of Conformation
Rendering of
UreB-Nb-Human Protein Complex at 0 or 100ns. Among them,
green shows UreB protein, and pink
shows Nb-Human
protein.
Combined with the analysis
of the movement trajectory of the
protein
complex, it can be seen that the three proteins have no obvious
displacement relative to the
UreB protein
during the whole simulation process, so the RMSD value is
relatively stable. Similarly, in the
early stage
of the simulation, the structure of the interaction between the
three proteins and UreB protein
has changed
to some extent, so the RMSD curve shows a certain fluctuation
trend in the first certain period
of the
simulation, and then gradually converges with the dynamic
simulation. On the whole, after 100ns
molecular
dynamics simulation, all three proteins formed a stable complex
with UreB protein, which can be
combined
with subsequent analysis to further observe the binding behavior
of the three proteins.
RMSF (Root Mean Square Fluctuation)
RMSF analysis showed the
structural adaptability of each
residue in
protein, which was used in this study to analyze the flexibility
and movement intensity of amino
acid
residues in protein during the whole simulation process. It can
be seen from the figure that the
overall
RMSF of the protein remained at a low level during the whole
simulation process. The reason why
the RMSF at
the N-terminal and C-terminal of the protein is large is that
this residue is at the head end of
the protein
and is less constrained by other domains in the protein, which
leads to greater flexibility and
has a
certain impact on its own stability. In addition, the reason why
the rest of the protein RMSF is
large may
be that there is some disturbance due to the combination of
small molecules, or because the
peptide itself
is flexible and has some disturbance during the
simulation.
By comparing the RMSF curves
of three protein complexes, it
can be seen
that the RMSF value of Nb-ScFv protein and Nb-Human protein is
higher than that of Nb-6 protein
after
forming complexes with Ureb protein. Which may be due to the
stronger interaction between
Nb-ScFv protein
and Nb-Human protein and UreB protein, so the flexible residues
at the interaction interface are
more
affected, which leads to the more obvious flexible jitter of
different residues of Ureb protein.
At the same
time, it was observed that the RMSF value of N-terminal of
Nb-Human protein increased obviously.
Combined
with the above analysis, it can be seen that this is because the
Nb-Human protein moved towards
UreB protein
in the first 20ns of simulation, so the RMSF value of this
random curly segment increased, which
is
consistent with the above analysis.
On the whole, the RMSF
values of the three proteins are low
in the
simulation process, so it can be considered that the proteins
have little jitter in the
simulation process,
keep stable vibration in the solution environment, and the
sampling and analysis are
reliable.
Figure 6
RMSF curve of each protein in the
simulation
process. A. The RMSF curve of Ureb. B. The RMSF curve of
Nb-6. C. The RMSF curve of Nb-ScFv.
D. The RMSF
curve of Nb-Human.
Rg (Radius of Gyration)
The Radius of Gyration (Rg)
can be used to characterize the
compactness
of protein structure, and it can also be used to characterize
the change of protein's peptide
chain
looseness in the simulation process. Combined with the above
analysis and Rg diagram analysis,
it can be
seen that the Rg curves of the complexes formed by Nb6 protein
and Nb_ScFv protein with UreB
protein
respectively remained stable throughout the simulation process,
and changed within 2.5nm and
2.8nm
respectively until the end of the simulation. After the complex
of Nb_Human protein and UreB
protein was
formed, the Rg curve showed an upward trend in the first 20ns of
the simulation, and then it
tended to
converge, maintaining the vibration around 2.75nm until the end
of the simulation.
Combined with the above
analysis, it can be seen that after
Nb6 protein
and Nb-ScFv protein form complexes with UreB protein, the
protein complexes remain relatively
stable during
the whole dynamic simulation process. As mentioned in the
aforementioned RMSD analysis, only the
mutual
orientation of the residues at the interaction interface of the
two proteins changes during the
simulation
process, so the RMSD curve remains in a low state, showing a
certain fluctuation trend. As the
skeleton
atoms of the proteins do not change obviously as a whole, the Rg
curve does not change
significantly, but
only keeps a certain vibration until the end of the simulation.
After the complex of Nb-Human
protein and
UreB protein was formed, the random curl of Nb-Human protein was
close to UreB protein in the
first 20 ns of
simulation, and the overall compactness became larger, so the Rg
curve increased in the first 20
ns of
simulation. After 20ns, with the convergence of the two protein
complexes after exercise, the Rg
curve also
tends to be stable, vibrating around 2.75 nm until the end of
the simulation.
On the whole, the Rg curve
of protein is stable in the
simulation
process and the sampling is reliable.
Figure 7 The
Radius of Gyration
Hydrogen bond number analysis
We analyzed the changes of
the number of hydrogen bonds
formed between
proteins with the simulation time and systematically
investigated the stability of protein
binding from the
perspective of interaction. From the graph analysis, it can be
seen that the number of hydrogen
bonds formed
between Nb-6 protein and Nb-ScFv protein and UreB protein
respectively remained relatively
stable and had a
certain amplitude during the whole simulation process. The
number of hydrogen bonds formed
between Nb-Human
protein and UreB protein showed a certain fluctuation trend in
the first 40ns of the simulation,
and
gradually increased as a whole, and stabilized after 40ns until
the end of the
simulation.
Based on the above analysis
and the analysis of the protein
complex's
trajectory, it can be seen that although Nb-6 protein and
Nb-ScFv protein have no obvious
displacement
relative to UreB protein during the whole simulation process,
due to the conformational
adjustment between
the residues at the interaction interface of the two proteins,
accompanied by the residue change
and
flexible jitter of the protein flexible binding cavity, this
movement leads to the continuous
formation and
destruction of hydrogen bonds between the residues at the
protein interaction interface, which
leads to a
certain vibration of the hydrogen bond number curve with the
change of simulation time, and the
whole is
relatively stable; Nb-Human protein appeared close to UreB
protein in the first 20ns of the
simulation, and
maintained this conformation until the end of the simulation, so
the number of hydrogen bonds
formed with
UreB protein showed a certain fluctuation trend in the first
40ns of the simulation, and the
whole showed an
upward trend due to the increase of protein interaction
interface. After 40ns, with the gradual
convergence
of the residues at the interaction interface between the two
proteins, the number curve of
hydrogen bonds
tends to be stable until the end of the simulation.
In addition, it can be noted
that although the number curve
of hydrogen
bonds changes obviously with the simulation, there are hydrogen
bonds between the two proteins
at all times,
which explains that the binding between proteins is stable from
the perspective of
force.
Figure 8 The
change of the number of
hydrogen bonds
formed between proteins with the simulation
time
MM-GBSA
After eliminating the cycle
of the simulated trajectory,
combined with
the aforementioned RMSD analysis, the movement trajectory of
protein complex in the period of
90ns-100ns was
extracted for binding free energy (MM/GBSA) analysis, and the
energy term was
decomposed.
The results showed that the
binding energy between protein
and UreB was
mainly van der Waals energy (Nb6 protein was -88.53 kcal/mol,
Nb-ScFv protein was -109.01
kcal/mol, Nb-Human
protein was -139 kcal/mol), and the electrostatic energy of Nb-6
protein was -26.58 kcal/mol,
Nb-ScFv
protein is -40.13 kcal/mol and Nb-Human protein is -43.19
kcal/mol), which indicates that
electrostatic
interaction is beneficial to protein binding, and proteins
interact mainly through van der Waals
interaction. The aforementioned hydrogen bond analysis and force
analysis show that there are
certain van
der Waals interactions and hydrophobic interactions between
different amino acids of each
protein, so van
der Waals can significantly facilitate the binding between
proteins and occupy a dominant
position; At the
same time, there is a certain electrical property at the
protein-protein interaction interface,
so there is
a certain electrostatic interaction at the protein-protein
interaction interface, which leads to
the
contribution of electrostatic energy to the binding. In
addition, we noticed that the solvation
energies of
Nb6 protein, Nb-ScFv protein and Nb-Human protein were 26.53
kcal/mol, 41.46 kcal/mol and 42.60
kcal/mol
respectively in this simulation, suggesting that the addition of
solvents was not conducive to
the binding
between proteins.
The overall binding energy
of Nb-6 protein is -88.59
kcal/mol, Nb-ScFv
protein is -107.68 kcal/mol, and Nb-Human protein is -139.59
kcal/mol), which indicates that the
binding
tendency of proteins is relatively large, and the relative order
of binding tendency is Nb-Human
> Nb-ScFv >
Nb-6, which further explains the binding between proteins from
the perspective of
energy.
Based on the above computer
simulation analysis, our
nano-antibody has a
good binding effect on Ureb protein of Helicobacter
pylori in theory, among which
Nb-Human has the
best binding effect. Unfortunately, due to the limitation of
biological safety, we can't
directly verify the
inhibitory effect of antibodies on Helicobacter pylori.
Later, we may consider
cooperating with
more professional research institutions to verify the actual
inhibitory effect of our antibodies
on
Helicobacter pylori.
Figure 9
Analysis of binding free energy
(MM/GBSA)
Table 2
Energy analysis in simulation
process
Reference
[1] Wang, J., Wolf, R. M., Caldwell, J. W., Kollman, P. A. &
Case, D. A. Development and testing
of a
general amber force field. J. Comput. Chem. 25, 1157-1174
(2004).
[2] Maier, J. A. et al. ff14SB: Improving the Accuracy of
Protein Side Chain and Backbone
Parameters from
ff99SB. J. Chem. Theory Comput. 11, 3696-3713 (2015).
https://doi.org:10.1021/acs.jctc.5b00255
[3] Mark, P. & Nilsson, L. Structure and dynamics of the TIP3P,
SPC, and SPC/E water models at
298 K. J.
Phys. Chem. A 105, 9954-9960 (2001).
https://doi.org:10.1021/jp003020w
[4] Salomon‐Ferrer, R., Case, D. A. & Walker, R. C. An overview
of the Amber biomolecular
simulation
package. Wiley Interdiscip. Rev. Comput. Mol. Sci. 3, 198-210
(2013).
https://doi.org:10.1002/wcms.1121
[5] Hou, T., Wang, J., Li, Y. & Wang, W. Assessing the
performance of the MM/PBSA and MM/GBSA
methods. 1.
The accuracy of binding free energy calculations based on
molecular dynamics simulations. J.
Chem. Inf.
Model. 51, 69-82 (2011).
https://doi.org:10.1021/ci100275a
[6] Genheden, S. & Ryde, U. The MM/PBSA and MM/GBSA methods to
estimate ligand-binding
affinities. Expert
opinion on drug discovery 10, 449-461 (2015).
https://doi.org:10.1517/17460441.2015.1032936
[7] Rastelli, G., Rio, A. D., Degliesposti, G. & Sgobba, M. Fast
and accurate predictions of
binding free
energies using MM‐PBSA and MM‐GBSA. J. Comput. Chem. 31, 797-810
(2010).
[8] Nguyen, H., Roe, D. R. & Simmerling, C. Improved Generalized
Born Solvent Model Parameters
for Protein
Simulations. J. Chem. Theory Comput. 9, 2020-2034 (2013).
https://doi.org:10.1021/ct3010485
[9] Weiser, J., Shenkin, P. S. & Still, W. C. Approximate atomic
surfaces from linear
combinations of
pairwise overlaps (LCPO). J. Comput. Chem. 20, 217-230 (1999).
Appendix
1. em.mdp：
; em.mdp - used as input into grompp to generate em.tpr
; Parameters describing what to do, when to stop and what to save
integrator = steep ; Algorithm (steep = steepest descent minimization)
emtol = 1000.0 ; Stop minimization when the maximum force < 1000.0 kJ/mol/nm
emstep = 0.01 ; Minimization step size
nsteps = 50000 ; Maximum number of (minimization) steps to perform
; Parameters describing how to find the neighbors of each atom and how to calculate the interactions
nstlist = 1 ; Frequency to update the neighbor list and long range forces
cutoff-scheme = Verlet ; Buffered neighbor searching
ns_type = grid ; Method to determine neighbor list (simple, grid)
coulombtype = PME ; Treatment of long range electrostatic interactions
rcoulomb = 1.0 ; Short-range electrostatic cut-off
rvdw = 1.0 ; Short-range Van der Waals cut-off
pbc = xyz ; Periodic Boundary Conditions in all 3 dimensions
2. ions.mdp：
; ions.mdp - used as input into grompp to generate ions.tpr
; Parameters describing what to do, when to stop and what to save
integrator = steep ; Algorithm (steep = steepest descent minimization)
emtol = 1000.0 ; Stop minimization when the maximum force < 1000.0 kJ/mol/nm
emstep = 0.01 ; Minimization step size
nsteps = 50000 ; Maximum number of (minimization) steps to perform
; Parameters describing how to find the neighbors of each atom and how to calculate the interactions
nstlist = 1 ; Frequency to update the neighbor list and long range forces
cutoff-scheme = Verlet ; Buffered neighbor searching
ns_type = grid ; Method to determine neighbor list (simple, grid)
coulombtype = cutoff ; Treatment of long range electrostatic interactions
rcoulomb = 1.0 ; Short-range electrostatic cut-off
rvdw = 1.0 ; Short-range Van der Waals cut-off
pbc = xyz ; Periodic Boundary Conditions in all 3 dimensions
3. nvt.mdp：
title = Protein-ligand complex NVT equilibration
define = -DPOSRES ; position restrain the protein and ligand
; Run parameters
integrator = md ; leap-frog integrator
nsteps = 50000 ; 2 * 50000 = 100 ps
dt = 0.002 ; 2 fs
; Output control
nstenergy = 500 ; save energies every 1.0 ps
nstlog = 500 ; update log file every 1.0 ps
nstxout-compressed = 500 ; save coordinates every 1.0 ps
nstvout = 500 ; save velocities every 1.0 ps
; Bond parameters
continuation = no ; first dynamics run
constraint_algorithm = lincs ; holonomic constraints
constraints = h-bonds ; bonds to H are constrained
lincs_iter = 1 ; accuracy of LINCS
lincs_order = 4 ; also related to accuracy
; Neighbor searching and vdW
cutoff-scheme = Verlet ; Buffered neighbor searching
ns_type = grid ; search neighboring grid cells
nstlist = 20 ; largely irrelevant with Verlet
rlist = 1.2
vdwtype = cutoff
vdw-modifier = force-switch
rvdw-switch = 1.0
rvdw = 1.2 ; short-range van der Waals cutoff (in nm)
; Electrostatics
coulombtype = PME ; Particle Mesh Ewald for long-range electrostatics
rcoulomb = 1.2 ; short-range electrostatic cutoff (in nm)
pme_order = 4 ; cubic interpolation
fourierspacing = 0.16 ; grid spacing for FFT
; Temperature coupling
tcoupl = V-rescale ; modified Berendsen thermostat
tc-grps = Protein_LIG Water_and_ions ; two coupling groups - more accurate
tau_t = 0.1 0.1 ; time constant, in ps
ref_t = 310.15 310.15 ; reference temperature, one for each group, in K
comm-grps = Protein_LIG Water_and_ions
comm-mode = Angular
; Pressure coupling
pcoupl = no ; no pressure coupling in NVT
; Periodic boundary conditions
pbc = xyz ; 3-D PBC
; Dispersion correction is not used for proteins with the C36 additive FF
DispCorr = no
; Velocity generation
gen_vel = yes ; assign velocities from Maxwell distribution
gen_temp = 310.15 ; temperature for Maxwell distribution
gen_seed = -1 ; generate a random seed
disre = Ensemble
disre_weighting = conservative
disre_fc = 10000
4. npt.mdp
title = Protein-ligand complex NPT equilibration
define = -DPOSRES ; position restrain the protein and ligand
; Run parameters
integrator = md ; leap-frog integrator
nsteps = 50000 ; 2 * 50000 = 100 ps
dt = 0.002 ; 2 fs
; Output control
nstenergy = 500 ; save energies every 1.0 ps
nstlog = 500 ; update log file every 1.0 ps
nstxout-compressed = 500 ; save coordinates every 1.0 ps
; Bond parameters
continuation = yes ; continuing from NVT
constraint_algorithm = lincs ; holonomic constraints
constraints = h-bonds ; bonds to H are constrained
lincs_iter = 1 ; accuracy of LINCS
lincs_order = 4 ; also related to accuracy
; Neighbor searching and vdW
cutoff-scheme = Verlet
ns_type = grid ; search neighboring grid cells
nstlist = 20 ; largely irrelevant with Verlet
rlist = 1.2
vdwtype = cutoff
vdw-modifier = force-switch
rvdw-switch = 1.0
rvdw = 1.2 ; short-range van der Waals cutoff (in nm)
; Electrostatics
coulombtype = PME ; Particle Mesh Ewald for long-range electrostatics
rcoulomb = 1.2
pme_order = 4 ; cubic interpolation
fourierspacing = 0.16 ; grid spacing for FFT
; Temperature coupling
tcoupl = V-rescale ; modified Berendsen thermostat
tc-grps = Protein_LIG Water_and_ions ; two coupling groups - more accurate
tau_t = 0.1 0.1 ; time constant, in ps
ref_t = 310.15 310.15 ; reference temperature, one for each group, in K
; Pressure coupling
pcoupl = Berendsen ; pressure coupling is on for NPT
pcoupltype = isotropic ; uniform scaling of box vectors
tau_p = 2.0 ; time constant, in ps
ref_p = 1.0 ; reference pressure, in bar
compressibility = 4.5e-5 ; isothermal compressibility of water, bar^-1
refcoord_scaling = com
comm-grps = Protein_LIG Water_and_ions
comm-mode = Angular
; Periodic boundary conditions
pbc = xyz ; 3-D PBC
; Dispersion correction is not used for proteins with the C36 additive FF
DispCorr = no
; Velocity generation
gen_vel = no ; velocity generation off after NVT
disre = Ensemble
disre_weighting = conservative
disre_fc = 10000
5. md.mdp
title = Protein-ligand complex MD simulation
; Run parameters
integrator = md ; leap-frog integrator
nsteps = 50000000 ; 2 * 5000000 = 100000 ps (100 ns)
dt = 0.002 ; 2 fs
; Output control
nstenergy = 5000 ; save energies every 10.0 ps
nstlog = 5000 ; update log file every 10.0 ps
nstxout-compressed = 5000 ; save coordinates every 10.0 ps
compressed-x-grps = System ; save the whole system
; Bond parameters
continuation = yes ; continuing from NPT
constraint_algorithm = lincs ; holonomic constraints
constraints = h-bonds ; bonds to H are constrained
lincs_iter = 1 ; accuracy of LINCS
lincs_order = 4 ; also related to accuracy
; Neighbor searching and vdW
cutoff-scheme = Verlet
ns_type = grid ; search neighboring grid cells
nstlist = 20 ; largely irrelevant with Verlet
rlist = 1.2
vdwtype = cutoff
vdw-modifier = force-switch
rvdw-switch = 1.0
rvdw = 1.2 ; short-range van der Waals cutoff (in nm)
; Electrostatics
coulombtype = PME ; Particle Mesh Ewald for long-range electrostatics
rcoulomb = 1.2
pme_order = 4 ; cubic interpolation
fourierspacing = 0.16 ; grid spacing for FFT
; Temperature coupling
tcoupl = V-rescale ; modified Berendsen thermostat
tc-grps = Protein_LIG Water_and_ions ; two coupling groups - more accurate
tau_t = 0.1 0.1 ; time constant, in ps
ref_t = 310.15 310.15 ; reference temperature, one for each group, in K
; Pressure coupling
pcoupl = Parrinello-Rahman ; pressure coupling is on for NPT
pcoupltype = isotropic ; uniform scaling of box vectors
tau_p = 2.0 ; time constant, in ps
ref_p = 1.0 ; reference pressure, in bar
compressibility = 4.5e-5 ; isothermal compressibility of water, bar^-1
; Periodic boundary conditions
pbc = xyz ; 3-D PBC
; Dispersion correction is not used for proteins with the C36 additive FF
DispCorr = no
; Velocity generation
gen_vel = no ; continuing from NPT equilibration
comm-grps = Protein_LIG Water_and_ions
comm-mode = Angular
disre = Ensemble
disre_weighting = conservative
disre_fc = 10000