Team:IISER Kolkata/Molecular Dynamics

Loading . . .

Introduction
Protein dynamics in water
Membrane protein interactions
Interaction with Nisin Resistance Protein
Key Takeaways
A word of caution

Introduction

Choice of systems

Considering that one of the key objectives of our modelling approach was to establish the efficacy of Nisin PV as compared to Nisin A, we went ahead with creating replicate systems with both antimicrobial peptides. The biological relevance of for this choice of systems is outlined in our design page.

Model setup

All systems were prepared using CHARMM GUI and were simulated using the GROMACS molecular dynamics engine. The PDB structure 1WCO was used as the reference structure for all simulations. The lipid-II moiety was stripped off during system preparation. For NisinA, the native 1WCO structure was used and to generate NisinPV the following mutations were performed within CHARMM GUI itself.:

Serine to Proline in 29th position
Isoleucine to Valine in 30th position

The below animation shows the system preparation of Nisin A and Nisin PV peptides in PyMol by invoking the aforementioned mutations.

Considering that we have non-traditional (or engineered amino acid residues in our protein), we tried to gather previously reported topologies and parameters for some of the non-standard amino acids. Amino acids for which we could not find compatible topologies and parameters, we used AMBER LEaP and antechamber program to generate parameters. The topology and parameter files that were used in our systems are present in our GitHub repository. For some non-traditional amino acids, we could not generate the topology and parameters using AMBER's LEaP program because of conflicting dihedral angle values. Thus, we had to resort to converting these amino acids to their closest traditional amino acid counterparts while taking care to maintain the connectivity.

All systems were equilibrated as per the following equilibration scheme:

Suggested Equilibrium Scheme [Reducing Force Constants] (5 cycles, 1 cycle = 50 - 100 ps)
	1 cycle	2 cycle	3 cycle	4 cycle	5 cycle	6 cycle
BB	10.0	5.0	2.5	1.0	0.5	0.1
SC	5.0	2.5	1.0	0.5	0.1	0.0
wforce	2.5	2.5	10	0.5	0.1	0.0
tforce	2.5	2.5	1.0	0.5	0.1	00
mforce	2.5	2.5	1.0	0.5	0.1	0.0
ion	10.0	0.0	0.0	0.0	0.0	0.0

Protein dynamics in water

We performed molecular dynamics simulations of NisinA and NisinPV in water to determine the stability of our Pro-Val mutant protein as compared to the wild type NisinA. Both systems were simulated using GROMACS MD engine and the CHARMM 36m force field and explicit water(TIP3P water model). The simulation temperature was 310K and 0.15 Molar NaCl(with neutralising ions) was used. Trajectory length is 30 nanoseconds.

Root-mean-square deviation (RMSD) Analysis

RMSD is a measure of how the protein residues change structurally at each frame of the simulation. Each residue is compared to its position to a reference frame, and its RMSD is plotted for every fame. The RMSD is given by:

$$ RMSD(t) = \left[ \frac{1}{M} \sum^{N}_{i=1} m_{i} | r_{i}(t) - r_{i}^{mof} | ^{2} \right] ^{\frac{1}{2}} $$

Where $M = \sum_{i} m_{i}$ and $ r_{i} (t) $ is the position of atom i at time t after least square fitting the structure to the reference structure.

Usually fitting is done with respect to the starting structure. It is not necessary to use the same set of atoms for the RMSD calculation and fitting, e.g. a protein is usually fitted on the backbone atoms but the RMSD can be computed for the backbone or for the whole protein.

We performed RMSD analysis using the GROMACS “gmx rms” tool. We are using RMSD to indicate any major structural fluctuations that would have occurred during the progression of the simulation. The RMSD was calculated for the entire protein molecule for both systems and the least square was fitted to the protein backbone. Additionally RMSD analysis can be used to check the similarity between both protein models.

Initially, both models quickly deviate from their reference structure by approximately 1 nm (this happens over the course of the first few nanoseconds during the heating, equilibration and the start of the production phase of the simulation). From then on in the simulation, the Nisin A model continues to deviate until it reaches approximately 2.5 nm and exhibits temporary stability. The Nisin PV model exhibits similar behaviour but shows a stark contrast beyond 15ns. Although it too initially deviates further to almost 3 nm early in the simulation it then returns to just under 1nm deviation from its reference structure and remains virtually constant for the rest of the simulation. There is some major fluctuation around the end of the trajectory for both models.

Solvent Accessible Surface Area (SASA) Analysis

We performed SASA calculations along with free energy of solvation calculation for both protein-water systems. The SASA was calculated with the GROMACS “sasa” tool which uses

It is interesting to note that major fluctuations in RMSD correspond to fluctuations in SASA. We further performed residue-wise SASA calculation to understand the dynamic of the system in greater detail. The residue-wise SASA was calculated using the VMD timeline analysis extension with the distance cutoff set as 0.14 nm.

Residue-wise SASA for NIsinA (left) and Nisin PV (right)

Connolly surfaces

The SASA is closely related to the concept of the solvent-excluded surface (also known as Connolly surface), which is imagined as a cavity in bulk solvent. The Connolly surfaces for both systems are as follows.

Secondary structure analysis

Residue-wise secondary structure for NIsinA (left) and Nisin PV (right). The calculation was performed using DSSP and visualised in VMD

Free energy surfaces

Biomolecular processes can be described in terms of the molecule's free energy:

$$ \Delta G (R) = - k_{B}T \left[\ln P(R) - \ln P_{max} \right] $$

where kB is the Boltzmann constant, P is the probability distribution of the molecular system along with some coordinate R, and Pmax denotes its maximum, which is subtracted to ensure that δG = 0 for the lowest free energy minimum. Popular choices for R (so-called order parameters) are the RMSD, Rgyr, number of hydrogen bonds or native contacts, or principal components.

The free energy is typically plotted along with two such order parameters, giving rise to a (reduced) free energy surface (FES). Here we use a custom python script to plot the FES as a function of dihedrals in the first case and as a function of RMSD and radius of gyration (RGyr)in the second case. The radius of gyration is a measure for the compactness of a structure and is calculated as:

$$ R_{g}=\left(\frac{\sum_{i}\left|r_{i}\right|^{2} m_{i}}{\sum_{i} m_{i}}\right)^{1 / 2} $$

We performed the radius of gyration calculation with the GROMACS “gmx gyrate” tool which computes the radius of gyration of a group of atoms as a function of time. The atoms are explicitly mass-weighted.

The free energy surface with Phi-Psi as reaction coordinates is similar for both NisinA and NisinPV. The position of the free energy minima for both systems are also in the same vicinity.

For the FES with RMSD-RGyr as reaction coordinates, we, unfortunately, do not have sufficient data points to draw any inferences. The distribution is very sparse. Performing longer and higher frequency simulation would have enabled us to draw inferences in this regard.

Diffusion behaviour (MSD)

The mean-squared displacement (MSD) is a measure of the deviation of the position of a particle with respect to a reference position over time. The MSD at time ${\displaystyle t}$ is defined as an ensemble average.

$$ \mathrm{MSD} \equiv\left\langle\left|\mathbf{x}(t)-\mathbf{x}_{\mathbf{0}}\right|^{2}\right\rangle=\frac{1}{N} \sum_{i=1}^{N}\left|\mathbf{x}^{(\mathbf{i})}(t)-\mathbf{x}^{(\mathbf{i})}(0)\right|^{2} $$

where $N$ is the number of particles to be averaged, vector $\mathbf{x}^{(\mathbf{i})}(0)=\mathbf{x}_{0}^{(\mathbf{i})}$ is the reference position of the $i$-th particle, and vector $\mathbf{x}^{(\mathbf{i})}(t)$ is the position of the $i$-th particle at time $t .{ }^{[3]}$

We calculated the MSD for both Nisin A and Nisin PV using the GROMACS “gmx msd” tool. The MSD was fitted to a line (excluding the first and last 10% of the trajectory to avoid the initial ballistic region and the noise at the end) to obtain the diffusion coefficient for the system. The stope of the MSD v/s time plot is indicative of the diffusion coefficient (with dimensions)

For n dimensions, MSD is related to diffusion coefficient (D) as :

$$ \mathrm{MSD}=2 n D t $$

The observed diffusion coefficients for both models are tabulated below:

Model	Diffusion Coefficient (1e-5 cm^2/s)
Nisin A	0.0470 (+/- 0.0051)
Nisin PV	0.3208 (+/- 0.0997)

We shall be using these values in our biofilm degradation model.

Our simulations showed that MSD for Nisin PV is significantly higher than MSD of Nisin A, thereby the diffusion coefficient of NisinPV would also be higher likewise. This aspect of Nisin PV might confer it an advantage against the pathogen as it can theoretically diffuse faster in the system and therefore cause faster action against the pathogenic bacteria.

Membrane protein interactions

We started out by setting up simulations of NisinA and NisinPV with a biological membrane to be able to delineate the dynamics and the interaction of these two bacteriocides with bacterial membrane. We went ahead with a membrane lipid composition of POPE: POPG = 3:1. This particular membrane composition is representative of a typical gram-positive bacterial cell membrane and has been used extensively in the literature for similar simulations involving bacterial membrane.

The following aspects of the nature of lipids are important to understand the nature of membrane protein interactions.

Lipid	Structure	Net charge (at pH 7)	Tm (K)
POPE	16:0 /18:1 PE	0	298
POPG	16:0 /18:1 PG	-1	269

POPE has 16 carbon atoms in one tail and 18 carbon atoms in the other with one unsaturation and the head group is Phosphatidylethanolamine. POPG acyl tail structure is the same as POPE and the head group, in this case, is Phosphatidylglycerol. The important point to note here is the head group charge and melting temperate (Tm). POPG is much more negatively charged as compared to POPE and has a significantly lower melting temperature. We started off with an inkling that the higher negative charge will play a crucial role in membrane protein interactions.

We set up two identical membrane systems and added NisinA to one and NisinPV to the other. All systems were created using CHARMM GUI. The thermodynamics and composition parameters were the same as previously outlined. Although membrane-protein systems are huge in terms of the number of particles to be simulated, we did not go for a coarse-grained approach as we intended to delineate the intricacies of the membrane-protein interactions. Therefore, we went ahead with fully atomistic systems and the simulations were run for 150 ns each. A summary of the setup is given as follows:

Force field	CHARMM36m
Water Model	TIP3P
Ions	0.15M NaCl (with neutralising ions)
Temperature	310.15K
Trajectory length	150 ns
Lipid composition on each leaflet	POPE:POPG = 3:1

Membrane bilayer thickness

Lipid membrane thickness is a standard parameter to check the stability of a membrane simulation. Additionally, if there is the generation of any membrane curvature, it should be indicated in the membrane thickness. We calculated the membrane thickness over time by defining a selection group consisting of the Phosphate head groups of the lipids in the bilayer. The local normal to the nearest phosphate neighbour would produce the local membrane thickness. This when averaged out over the entire membrane gives one value of membrane thickness as a function of time. FATSLiM is a nifty python-based tool that can be used to calculate membrane thickness for courage-gained simulation. However, since we performed atomistic simulation we had to play around with the distance cutoff for the local normal and leaflet identification in the FATSLiM source code.

The membrane bilayer thickness remains more or less constant and within the margins of thermal fluctuations. Certain significant dips in the bilayer thickness might be indicative of the generation of local membrane fluctuations.

Area per lipid (APL) calculation

Membrane APL is another parameter to describe the dynamic of membrane systems. APL is indicative of not only membrane stability but also indicates probate regions of interaction between the lipid bilayer and the membrane. APL calculation is similar to the thickness and used a section group consisting of the phosphates. The Area of the lipid boundary is approximately defined by the coordinates of the phosphates in the head groups. This area divided by the total number of lipids in the bilayer gives the area per lipid. We again used a modified implementation of FATSLiM for this calculation

The area per lipid for both membrane-protein systems is virtually constant within the margins of thermal fluctuations. However, the APL for NisinPV is marginally higher than the APL for NisinA. The more significant fluctuations of APL could be indicative of membrane-protein interactions.

Hydrogen bonds between protein and lipid head-groups

Hydrogen bonds provide most of the directional interactions that underpin protein folding, protein structure and molecular recognition. Hydrogen bonding confers rigidity to the protein structure and specificity to intermolecular interactions. An H-bond is formed by the interaction of a hydrogen atom that is covalently bonded to an electronegative atom (donor D) with another electronegative atom (acceptor A). The accepted geometry for an H-bond is a distance of less than 3.5 $Å$ between hydrogen D and A and an D-H-A angle of 180° ± 30°.

Having established that the membrane system is stable and there are instances of membrane-protein interaction. We went ahead to try to explain the nature of these interactions. We started plotting the number of Hydrogen bonds between the protein and the lipids over time. The analysis was performed using the “gmx hbond” tool in GROMACS. This analyses all H-bonds existing between two groups of atoms (which must be either identical or non-overlapping) or in specified D-H-A triplets using specified D-A distance and D-H-A angle criteria. The distance cutoff or hydrogen bonds was set up to be 0.35 nm so that even slightly weaker hydrogen bonds could be taken into account.

The time evolution of the total number of H-bonds between the selected groups shows that the number of hydrogen bonds between the protein and the POPG is similar for both NisinA and NisinPV. However, the same is drastically different for POPE wherein Nisin PV forms significantly lower H-bond contacts with the lipid. It must be noted here that POPE is three times more abundant in the system. Still, the number of H-bond contacts between NisinPV and the POPE is significantly lower than NisinA. In fact, even the total membrane H-bond contacts for NisinPV as lower than the NisinA-POPE contacts.

To have a deeper look at this behavior, we plotted the residue-wise H-bond contacts for Nisin and realised that the majority of contacts formed between the protein and the POPG lipids involve the first few residues in the N-terminal region of Nisin (which in our case in unmodified between NisinA and NisinPV) specifically in residue 1 to 11 region and the contacts are formed with the oxygen atoms in the POPG head group.. Here, we would like to invoke the fact that the net charge of POPG at pH7 id -1 whereas for POPE is 0. We know from literature that the Nisin 1-11 region is responsible for interactions with lipid-II and is in fact responsible for the capture of membrane-bound lipid-II, causing membrane degradation. Nisin interacts with the pyrophosphate group of lipid-II (as reported by Panina et al.). We suspect that the interaction with POPG is along the same lines as the interaction with Lipid-II given the charged nature of the groups involved. We wanted to perform Potentials of Mean Force (PMF) calculations to further explore this aspect, but due to a shortage of computational resources and cluster time, we could not perform any enhanced sampling simulations in this regard.

Given here is an animation of the residues involved in the interaction of Nisin with Lipid-II. The animation has been made with the data published by Panina et al.

Free energy surfaces

We calculated the free energy for NisinA and NisinPV in the membrane environment just the way we did in the section where we looked at the protein dynamics in water. We took two cases of reaction coordinates, the first being the Phi-Psi angles (Ramachandran plot) and the second being RMSD-RGyr.

The free energy surface with Phi-Psi as reaction coordinates is similar for both NisinA and NisinPV. The position of the free energy minima for both systems are also in the same vicinity.

The free energy surface with RMSD-RGyr as the reaction coordinates is more dispersed (Spread out) for NisinA as compared to NisinPV. The free energy minimums are also positioned almost linearly in both systems, indicating that RMSD and RGyr are positively correlated in terms of the free energy associated with it. The maximum RMSD and maximum RGyr are also significantly higher for the NisinPV system (almost two times that of the NisinA system). Although, these extremes are sparsely populated.

Interaction with Nisin Resistance Protein

The data in this section is adapted from the simulations by Field, Des et al. They carried out MD simulations to predict the difference in molecular dynamics of the NSR protein in its interaction with the nisin C-terminus (residues 22–34) for nisin-A and nisin PV models.

Molecular docking

We tried to replicate the setup used by Des et al. For the NSR (Nisin Resistance protein) molecule, the PDB file 4Y68 was used and for the Nisin the 1WCO PDB file was used. NisinA structure was mutated using PyMol and a PRO-VAL mutation was added to the 29th and 30th residue respectively. To determine a starting configuration for the NSR-nisin complex a docking program Autodock for ligand and protein binding was used where the Nisin (A/PV) was set as ligand and NSR was set as receptor. We first performed the docking of NisinA with the tunnel region of NSR and selected the lowest energy stage as the starting configuration. The docs location was observed to be consistent with the ligand-binding site of Nisin as discussed in literature. The ligand-binding site for NisinA was noted and NisinPV was docked using the same coordinates.

The free energies of the configurations obtained from docking are as follows:

position	Affinity (kcal/mol)
	Nisin A	Nisin PV
1	-9.3	-9.5
2	-8.9	-9.1
3	-8.7	-9
4	-8.7	-8.9
5	-8.7	-8.8
6	-8.6	-8.8
7	-8.6	-8.7
8	-8.5	-8.4
9	-8.5	-8.4

The 1st position is the best docking position for both molecules. Subsequent docking positions were also found by this study. This shows both Nisin A and Nisin PV docks similarly to NSR protein, which indicates that binding is not the mechanism by which NSR cleaves Nisin A.

The above figure elucidates the docking of the 22-34 N-terminus end of NisinA (cyan) and Nisin PV (magenta) to the tunnel region of NSR. Only the lowest free energy structures are shown.

RMSD Analysis

Initially, both models quickly deviate from their reference structure by approximately 3 Å (this happens over the course of the first 0.3 ns during the heating, equilibration and the start of the production phase of the simulation). From then on in the simulation, the nisin A model continues to deviate until it reaches approximately 5 angstroms. In contrast, the nisin PV model exhibits different behaviour. Although it too initially deviates further to almost 5 angstroms early in the simulation it then returns to just under 4 angstroms deviation from its reference structure and remains virtually constant for the rest of the simulation.

This evidence of greater structural change in the nisin A model that is less apparent in the nisin PV model supports the findings of the binding energy analyses for both bond and total energies. Indeed, it is clear that the nisin PV model is more durably bound to the NSR protein and is certainly more robust than the nisin A model.

Hydrogen bond occupancy

These simulations revealed the hydrogen bond occupancy (defined as whether a hydrogen bond is present or absent at every frame over the duration of simulation) for all nisin residues greater than 20%. It is evident that the nisin PV model exhibits higher occupancy values than the nisin-A model for residues corresponding to the region 27-29. In particular, there is a hydrogen bond with 47.65% occupancy between the sidechain oxygen of proline in the nisin-PV model and the hydrogen from the nitrogen sidechain on residue 265 (asparagine) in NSR that may play a role in the inhibition of proteolytic cleavage for nisin PV.

Free energy per residue

Analysis of the total binding energies for each residue in both nisin A and nisin PV models reveals that at the critical residue 29, the total energy for proline is seven times greater than serine (4.6292 kcal/mol for PRO and 0.8775 kcal/mol for SER).

Key Takeaways

Considering the binding of Nisin(A/PV) to NSR, there is compelling evidence that NSR has a much more difficult time cleaving the bond at the 29th in the case of NisinPV.
The total energy at the 29th residue is approximately 7 times larger for NisinPV as compared to NisA. This is consistent with the results for the biofilm assay where the minimum inhibitory concentration (MIC) for NisinA was measured to be approximately 7 times that of NisinPV. We shall use this information in our mathematical model for the degradation of biofilm where we would make the assumption that the rate of degradation of the biocide is higher for NisinPV as compared to NisinA. biofilm assay here.
Both NisinA and PV majorly interact with the highly negatively charged lipids in the biological membrane. This would mean that NisinPV has at least equal membrane degradation potential (in the absence of any resistance proteins of course) to NisinA.
Our simulations show that NisinA interacts with lipids with a neutral charge (like POPE) more preferentially than NisinPV. It is yet to be determined whether this differential interaction with POPE would affect the membrane destabilisation and degradation.
Our simulations showed that MSD for Nisin PV is significantly higher than MSD of Nisin A, thereby the diffusion coefficient of NisinPV would also be higher likewise. This aspect of Nisin PV might confer it an advantage against the pathogen as it can theoretically diffuse faster in the system and therefore cause faster action against the pathogenic bacteria.

A word of caution

Even though we tried to model the topologies and parameters of all the non-standard amino acids in the structure, we had to settle for replacing some of the non-standard amino-acid with standard amino acids. This resulted in higher fluctuations in our model as compared to other modes in the literature.

References

GROAMCS

Berendsen, et al. (1995) Comp. Phys. Comm. 91: 43-56. (DOI, Citations of this paper)
Lindahl, et al. (2001) J. Mol. Model. 7: 306-317. (DOI, Citations of this paper)
van der Spoel, et al. (2005) J. Comput. Chem. 26: 1701-1718. (DOI, Citations of this paper)
Hess, et al. (2008) J. Chem. Theory Comput. 4: 435-447. (DOI, Citations of this paper)
Pronk, et al. (2013) Bioinformatics 29 845-854. (DOI, Citations of this paper)
Páll, et al. (2015) Proc. of EASC 2015 LNCS, 8759 3-27. (DOI, arxiv, Citations)
Abraham, et al. (2015) SoftwareX 1-2 19-25 (DOI, Citations)

VMD - Visual Molecular Dynamics

http://www.ks.uiuc.edu/Research/vmd/

CHARMM- GUI

S. Jo, T. Kim, V.G. Iyer, and W. Im (2008) CHARMM-GUI: A Web-based Graphical User Interface for CHARMM. J. Comput. Chem. 29:1859-1865

B.R. Brooks, C.L. Brooks III, A.D. MacKerell, Jr., L. Nilsson, R.J. Petrella, B. Roux, Y. Won, G. Archontis, C. Bartels, S. Boresch, A. Caflisch, L. Caves, Q. Cui, A.R. Dinner, M. Feig, S. Fischer, J. Gao, M. Hodoscek, W. Im, K. Kuczera, T. Lazaridis, J. Ma, V. Ovchinnikov, E. Paci, R.W. Pastor, C.B. Post, J.Z. Pu, M. Schaefer, B. Tidor, R. M. Venable, H. L. Woodcock, X. Wu, W. Yang, D.M. York, and M. Karplus (2009)

CHARMM: The Biomolecular Simulation Program. J. Comput. Chem. 30:1545-1614 J. Lee, X. Cheng, J.M. Swails, M.S. Yeom, P.K. Eastman, J.A. Lemkul, S. Wei, J. Buckner, J.C. Jeong, Y. Qi, S. Jo, V.S. Pande, D.A. Case, C.L. Brooks III, A.D. MacKerell Jr, J.B. Klauda, and W. Im (2016)

CHARMM-GUI Input Generator for NAMD, GROMACS, AMBER, OpenMM, and CHARMM/OpenMM Simulations using the CHARMM36 Additive Force Field. J. Chem. Theory Comput. 12:405-413

FATSLiMhttps://fatslim.github.io/
Secondary structure using DSSP https://swift.cmbi.umcn.nl/gv/dssp/
Nisin-NSR interactions

Field, Des et al. “Bioengineering nisin to overcome the nisin resistance protein.” Molecular microbiology vol. 111,3 (2019): 717-731. doi:10.1111/mmi.14183

CHARMM Parameters and topologies for engineered aminoacids

Van de Ven, F J et al. “NMR studies of lantibiotics. The structure of nisin in aqueous solution.” European journal of biochemistry vol. 202,3 (1991): 1181-8. doi:10.1111/j.1432-1033.1991.tb16488.x

Nisin-LipidII interactions

Mulholland, Sam et al. “Docking and molecular dynamics simulations of the ternary complex nisin2:lipid II.” Scientific reports vol. 6 21185. 18 Feb. 2016, doi:10.1038/srep21185

AutoDock

http://autodock.scripps.edu/resources/references
http://autodock.scripps.edu/resources/references