We developed a molecular dynamics (MD) simulation of our main protein PETase binding to PET to
understand the reaction it catalyzes.
This page gives details on the MD simulation we developed and the things we have learned from it.
What is MD?
A 200ps simulation of water molecules
Molecular dynamics is a method to computationally study the physical movements of atoms and
molecules. It starts with an ensemble of atoms and molecules that all have well-defined
Next, the respective physical laws are evaluated between these atoms and the forces
van-der-waals, covalent bonds,) are calculated. Newton’s second law of motion teaches us that
acting on masses will result in acceleration. In MD, we now allow for the atoms to move
the forces that act upon them for a really short time (usually 1–2 femtoseconds). This timestep
simulation results in new coordinates for the atoms. These new coordinates are now evaluated to
provide new atomic forces. With these updated forces we can run the ensemble for another
This procedure has to be repeated for 1 million times to make for an overall simulation of one
nanosecond. The nice thing about MD is that, given that the correct forces are calculated, the
system will behave exactly like a real pendant. But now this version of reality can be looked at
an incredible zoom and in super slow-motion. If we want to, we can observe a single molecule of
water slowed down 1000000000000 times. With this powerful tool, we can run virtual experiments
would be impossible to conduct in real life. We can simulate extreme pressures and temperatures,
toxic, expensive or hard to synthesize compounds. In the past, MD was a niche method for small
systems, but as computing power develops exponentially, MD is making its way into mainstream
and anyone with a modern PC can now do their own simulations.
Building the Ensemble
The ensemble consists of the enzyme PETase, the PET strand as well as the ions and water molecules surrounding it.
We used the softwares NAMD (Nanoscale Molecular Dynamics) and VMD (Visual Molecular Dynamics) –
developed and run by the University of Illinois at Urbana–Champaign².
When modeling any kind of protein, the first thing one needs is the tertiary structure. This can be determined by X-Ray Cristallography, NMR or even cryo-electron microscopy. In our case, we found an X-ray protein structure of PETase with a resolution of 1.54 Å in the protein data bank (pdb)¹. We downloaded the .pdb file, which gives information about all the proteins’ atoms as found in a protein crystal. Since modeling our protein by itself would have been pretty boring, we also needed a structure file for the PET polymer. Luckily, the University of Illinois at Urbana–Champaign had our backs once again. They had already published a paper on PET nanopores⁴ and with this, uploaded .pdb files of PET monomers, dimers, nonamers and bulk structures. Since monomers and dimers were too small to create a realistic model of PETase activity and because large bulk structures can result in long calculation runtimes, we decided to use the nonamer structure for our interaction studies.
Since .pdb files only contain atom coordinates and no information about the interatomic forces, we had to find matching topology files first. These topology files tell our program about what types of connectivities and interatomic forces to expect between the atoms, given their relative coordinates. We used the CHARMM force field topology files, because that is what is recommended by the programmers of VMD. In VMD, using the .pdb files and the CHARMM topology files, we created the protein structure files (.psf). We merged the structures of the enzyme and the PET-strand and brought them into the right orientation and proximity. Next, it was important to solvate the ensemble, as we did not want to model PETase and PET in space. We decided to create a water box with 20 Å padding around our structures, to ensure that they do not leave their compartment. This way, we ended up with a roughly cubic simulation volume with dimensions of around 100 Å. To have our virtual experiment as close to reality as possible, we also wanted to add the right amount of ions, which, in our wet lab experiments, were 80 mM NaCl. Luckily, in VMD we can easily add ions to our simulation volume.
Running the Simulation
A 20 ns simulation of PET bound to PETase
Now that we had our atom ensemble set up, it was time to start running some simulations. All
simulations were done in NAMD, while the set-up and visualization happened in VMD. We thought
20 ns should be enough time for the binding to happen. Specifically, we ran an NPT simulation,
meaning we kept the number of particles, the pressure (1 atm), and the temperature (310 K)
We used Langevin dynamics for both the pressure and temperature control. Furthermore, we ran the
cubic simulation with periodic boundary conditions and used particle mesh ewald method (PME) for
evaluating the long range electrostatics.
We were lucky enough to observe our enzyme binding to its substrate. We had started the simulation from a state where the PET strand is already aligned with a binding cleft that can be found on the surface of the PETase enzyme. After x ns of simulation and indecisive random movement, the strand finally slipped into the binding pocket and got into the right conformation and proximity to the active site. Here, we speculate, the strand would now usually be cut by the enzyme. However, in MD it is not possible to break bonds, so we just had to be content with observing the binding.
In this gif we can see the PET strand moving into the cleft on the surface of PETase, where it is now accessible for cleavage by the catalytic triad.
We were able to verify that the strand really was in the correct position by comparing it to the crystal structure of the substrate bound protein¹.
Crystal structure of PET bound to PETase¹
The strand fully immersed in the hydrophobic cleft.
It is noteworthy that the PETs aromatic ring is perfectly sandwiched by Y87 and W185, where it is stabilized by π-π interactions¹.
π-π interactions stabilize the PET ring structure between the two aromatic amino acids Y87 and W185.
We also learned a lot from
having a close up look at the active site and from comprehending it’s catalytic mechanism,
works analogously to the standard textbook example mechanism of the serine proteases. We can
see the catalytic triad consisting of the Ser160 Asp206 His237 residues in close proximity
ester bond it is about to cleave.
The catalytic triad is very close to the substrate.
The Sequence Asn205 Asp206 Ser207 is a signal for the protein to be N-glycosylated in Chlamy's endoplasmatic reticulum.
The Glycosylation site is in immediate proximity to the PET binding site.
Something that also came to our minds when looking at the enzyme was how possible glycosylations
would impact substrate affinity and cleaving activity. We noticed quite a lot of Asn residues in
proximity of the binding site/active site. In eukaryotes like Chlamy, these residues will be
glycosylated in the ER if found in a sequence of Asn, X, Thr/Ser, where X can be any amino acid
proline. With Asn205, just this was the case. We can imagine that a glycosylation in a place
close to the active site would severely hinder the binding of the substrate. Unfortunately, we
noticed this only very late in our project and did not have the time to check these implications
the wet lab. It is something that we will have to look into in the future.
Comparing the Isozymes
The arginine residue is bulky and positively charged, repulsing the PET strand from the binding site.
The uncharged alanine residue blends in with the hydrophobic cleft.
In our wetlab experiments, we used the wild type PETase as well as one with three mutations:
for increased activity and improved binding of longer substrates¹ as well as S238F and W159H for
narrowing the active site and thus increasing activity³. In the protein data bank we also found
crystal structure of a PETase containing the R280A mutation. So we wanted to compare these two
variants side by side to see whether the point mutation really makes a difference. In the side
side view we can see that while both isozymes bind to the strand, the binding in the R280A
tighter, and that the alanine is less bulky than the arginine. Furthermore, the neutral alanine
less repulses the strand than the positively charged Arginine. After all, the binding is driven
mostly by hydrophobic interactions. These results are in agreement with the paper that states
the R280A mutant has a 32% higher activity on PET than the wild type¹.
Our modeling and analysis was focused on achieving a better theoretical understanding of the
catalyzed by PETase. Our model as well as the insights we have gained from it will benefit
iGEM teams who choose to work with PETase. Having the opportunity to see, turn, move, color,
visible/invisible any part of the protein substrate complex is incredibly valuable for all work
on the protein. We believe that anyone planning to intensively study a protein and the reaction
catalyzes should get a 3D look on their subject of interest, if not a real MD simulation.
-  Structural insight into molecular mechanism of poly(ethylene terephthalate) degradation; Seongjoon Joo, Sang Yup Lee & Kyung-Jin Kim; Nature Communications volume 9, Article number: 382 (2018)
-  James C. Phillips, Rosemary Braun, Wei Wang, James Gumbart, Emad Tajkhorshid, Elizabeth Villa, Christophe Chipot, Robert D. Skeel, Laxmikant Kale, and Klaus Schulten. Scalable molecular dynamics with NAMD. Journal of Computational Chemistry, 26:1781-1802, 2005.
-  Characterization and engineering of a plastic-degrading aromatic polyesterase; Austin HP, Beckham GT; Proc Natl Acad Sci U S A. 2018 May 8;115(19):E4350-E4357. doi: 10.1073/pnas.1718804115. Epub 2018 Apr 17.
-  Eduardo R. Cruz-Chu, Thorsten Ritz, Zuzanna S. Siwy, and Klaus Schulten. Molecular control of ionic conduction in polymer nanopores. Faraday Discussion, 143:47-62, 2009. (PMC: 2907245)