Team:Freiburg/Model
Modelling of mirror-image D-proteins
Fig. 1: Schematic representation of the L-to-D Converter's function. In cyan the NMR structure of L-PSMα3 (5kgy) obtained from the Protein Data Bank. By inversion of the X-coordinate of each atoms' coordinates a model of D-PSMα3 (red) is generated.
In our project, we want to promote the power of D-peptides and to accelerate research of D-peptides. Computational modelling is an important tool for experimental design and system analysis, therefore we wanted to make computational modelling accessible for D-proteins. For this purpose we wrote the program L-to-D Converter and made it available for the research community on github. The L-to-D Converter generates the mirror-image D-protein structure to any given L-protein structure. Protein structure models are written as .pdb files, that list the coordinates of every atom in the protein. The L-to-D Converter works by inverting the algebraic sign of the X-coordinate of every atom present in the .pdb file corresponding to mirroring the whole molecule along the Y-Z plane of a coordinate system.
We verified the accuracy of this method by comparing the crystal structure of D-Monellin, which was chemically synthesized from D-amino acids (PDB id:2q33)1, with a D-Monellin structure which was generated by the L-to-D converter from a crystal structure of natural L-Monellin (PDB id:1mol)2. Except for minor differences due to the diverging crystallization conditions of the two structures the experimentally obtained structure and the modeled one are almost identical, with an RMSD of 0.051 between the two structures.
Fig. 2: Structures of L-Monellin (blue PDBid: 1mol) and D-Monellin (red PDBid: 2q33), The structure of L-Monellin was converted to a model for D-Monellin and aligned with the crystal structure for D-Monellin using PyMOL.
Obtaining realistic models of proteins and peptides through molecular dynamics simulations with GROMACS3
In our wet lab MIPD, we identified novel D-peptides which bind to our target PSMα3. Because these peptides were derived from a random phage library we knew nothing more about them than their sequence. To gain more insights on those ligands and how they might interact with PSMα3 we decided to model them in silico.
Fig. 3: NMR structure of PSMα3 (PDBid: 5kgy), visualized with PyMOL
One of the biggest constraints of computational modelling is the transferability of simulated results and predictions to the reality under wet lab conditions. This also applies for the structural modelling of proteins and peptides. The starting point for such a model is usually a 3D structure which can be obtained from the Protein Data Bank, in the case of PSMα3 we chose the structure PDB id: 5kgy.
However, proteins, especially short ones are highly dynamic molecules and they can exist in various conformations constantly changing over time influenced by their surroundings. Thus they cannot be accurately represented by a static crystal structure. Furthermore these crystal structures show a highly artificial state of the protein, densely packed with many copies of the same molecule allowing barely any movement or interaction with a molecule other than itself.
To avoid working with a structure that is biased by these highly artificial, unnatural conditions and to ensure a high comparability of our simulations with wet lab reality, we performed a molecular dynamics simulation of our target PSMα3 using the program GROMACS4.
GROMACS simulates the movements of atoms in a given protein at every time point by solving Newton’s equation of motion to predict its movement through a solvent and conformational shifts over time. The parameters that determine the kinetics of different atoms depend on their chemical identity but also on type and angle of the bonds connecting them. All this information is provided by so called force fields.
Attempting a simulation of the toxin PSMα3 of the multiresistant Staphylococcus aureus we encountered the problem that the N-terminal formylation of this protein couldn’t be recognized by GROMACS because it was lacking an AMBER99 force field containing the properties of N-formylmethionine. We managed to solve this problem by implementing a formylmethionine that we parameterized by adapting a force field for formylvaline (which was provided by Dr. Steffen Wolf, originally from Na, S. et. al.5). For further detail on the process visit our methods page. With this we facilitated the faithful modelling of bacterial proteins which are important drug targets in antibiotics research.
Hit the button to see the adapted charge and bond parameters for N-formymethionine.
Fig. 4: Molecular Dynamics simulation of the movements of PSMα3 over a timespan of 10 ns and the most occuring state of the molecule derived from a 50 ns molecular dynamics simulation through clustering.
In our molecular dynamics simulation with GROMACS we simulated the movement of PSMα3 over a timespan of 50 nanoseconds, generating 5000 sequential structure files depicting its conformation at different timeframes. This was repeated 4 times for greater statistics. Finally, to obtain the most common conformation of PSMα3 out of all these structures, we clustered them according to root mean square deviation (RMSD) and extracted the most-occurring states of the molecule for docking.
Prediction of 3D structure models of the peptide ligands from mirror-image phage display
Obtaining a computational model for the ligands which were identified in mirror-image phage display (MIPD) required another approach, since there is no crystal structures available for these ligands. Therefore we constructed 3D structures from scratch using PyMOL. To verify these structures and to extract the most more realistic state we performed molecular dynamics simulations with GROMACS.
Fig. 5: To generate a faithful model of peptide ligands based on an amino acid sequence in the one letter amino acid code we use the program PyMOL to create an artificial peptide model as a first step. To obtain a faithful model of the peptide structures a short molecular dynamics simulation is performed and the most ocurring state is extracted.
Modelling of the mechanisms behind concentration dependent behaviour of PMSα3 cytotoxicity
Having obtained a faithful modelling strategy for both D-PSMα3 and its L-ligands we aimed to analyze how they bind to each other and if the ligands we obtained from wet lab MIPD would interact with the functionally important residues of PSMα3. Cheung et. al. found that especially the lysine residues in positions 6, 12 and 19 were important for the cytotoxic function of PSMα3. Additionally they showed that the C-terminal asparagine residues have an impact on the FPR2 receptor dependent neutrophil activation6.
Using the program AutoDock Vina we performed molecular docking of the peptide ligands from MIPD against PSMα3 and consequently analyzed the predicted intermolecular contacts. AutoDock Vina works by searching for a global minimum of the binding free energy of two molecules using an iterative algorithm that starts independently from thousands of different initial positions. The following models of ligand-target interaction are the most likely as they yielded the highest absolute value of binding free energy and were identified multiple times in independent simulations.
Our model of the ligands’ interaction with PSMα3 showed that all three ligands interact with several of the functionally important residues. Therefore, they might have an impact on toxicity by masking these residues. Furthermore all ligands interact with either the formylated N-terminus or the C-terminus which are both important for FPR2 receptor activation.
Based on these quite promising in silico predictions, we tested the effect of the ligands on the lytic and inflammatory properties of PSMα3 in wet lab and validated our predictions on the ligands’ effect on toxicity.
Evaluation of the mechanisms behind concentration dependent behaviour of PMSα3 cytotoxicity in wetlab
Fig. 6: Association of L- to D-PSMα3. L-PSMα3 (blue) was docked to D-PSMα3 (red) using AutoDock Vina after molecular dynamics simulation (rigid docking, exhaustiveness 10). Residues involved in the binding interface are highlighted in darker colors.
In our wet lab characterization of PSMα37. we observed that its cytotoxicity is only concentration dependent in the low micromolar range but not in the millimolar range. We suspected that this paradox drop of cytotoxicity at higher concentrations might be due to the formation of amyloid structures as it is described in the literature for PSMα3 (Tayeb-Fligelman, E Science 2017).
To verify the assumption that L-PSMα3 molecules bind to each other we docked them. The binding free energy of the 2 molecules was -12.5 kcal/mol, confirming that the association of L-PSMα3 molecules with each other is favorable and reflecting our wet lab data.
This gave us the idea to test the binding between the mirror-image enantiomers, L- and D-PSMα3 by docking. Interestingly, our simulations predicted that these two associate to each other with a quite high binding free energy of -19.2 kcal/mol. For small peptides like PSMα3 this is a quite high value, indicating that L- and D-PSMα3 might also form some kind of fibrils
We decided to explore the association of L- and D-PSMα3 by conducting further experiments. Firstly, we aimed to verify the predictions we made in silico by wet lab binding analyses. The binding of L-PSMα3 to D-PSMα3 was measured via Bio Layer Interferometry (BLI) using the Octet RED96 system (FortéBio, Molecular Devices). BLI is based on optical measurements of the refractive index of a target-coated sensor. Depending on ligand binding to this sensor the biological layer becomes thicker causing a change in the refractive index which can be measured by the Octet RED system.
For the BLI assay, biotinylated D-PSMα3 was bound to streptavidin sensors (Fig. 7) and the association and dissociation of L-PSMα3 to its D-form was measured. a binding constant KD of 13.1 +/- 2.5 µM was calculated using the software ForteBIO.
We could show that the heterochiral association of L- and D-PSMα3 also has an effect on cytotoxicity, even at a concentration of 8 µM where there was amyloid formation between homochiral PSMα3 molecules.
Upon further research we found also a cutting-edge paper mentioning this phenomenon in the context of protein crystallization8.
Fig. 7: Binding of L- to D-PSMα3. Real time BLI measurement of the association of 50 µM L-PSMα3 to a D-PSMα3 coated sensor.
Fig. 8: Influence of the association of L- and D-PSMα3 on cytotoxicity. The cytotoxicity of D- and L-PSMα3 on Jurkat T-cells with or without 24h prior preincubation was tested using lactate dehydrogenase assay. Depicted is the relative cytotoxicity as compared to Triton X-100 lysis buffer.
Establishing a model for our MIPD system in silico
We also aimed to model our MIPD system itself in silico. Having established representative model structures for the target D-PSMα3 we still needed models for a peptide library and a method to determine high-affinity binders from non binding peptides to complete our model.
Generation of two random libraries of potential ligands
Fig. 9: Schematical depiction of the helix extraction from proteins.
A diverse library of potential ligands is key to every kind of phage display, wet lab or in silico. Since there are no in silico peptide libraries openly available we decided to generate one on our own which we also provide for the iGEM community.
The largest data set of 3D crystal structures is the Protein Data Bank (PDB), however it contains proteins not peptides. To render this huge amount of information useful for our purposes we extracted only the helical segments of all proteins in the database.
We focused on the helical parts since most binding interfaces between proteins are located in helices9 and following the advice of Dr. Steffen Wolf who pointed out the relative stability of this secondary structure and the fact that helices can form more independently of their surrounding, starting from peptide lengths exceeding 8 amino acids. The folding of the extracted regions is an important matter because the library depends on the assumption that the extracted peptides would fold in a similar manner as they do within the protein complex. To address this issue we generated two libraries applying different approaches. Since these libraries are incredibly large in terms of data volume we decided not to provide the library itself to the iGEM community but rather the tools to extract these libraries.You can find the PDB Library tools and more useful programs on our software page.
Unsliced helix peptide library
For our first library we exclusively extracted helices between the sizes of 8 to 12 amino acids which are flanked by unstructured regions in proteins. This conformation suggests that the helicity is a property of the extracted sequence itself and that it would fold similarly if isolated. Therefore we hypothesize that the potential ligands in this unsliced library have a high chance to fold into a conformation which is similar to the structure used for modelling. This makes them decent candidates for computational screening with the aim to transfer in silico results to the wet lab. The total number of unsliced peptides from the PDB is 4 million, duplicates excluded. For reasons of computational capacity we performed all further mentioned docking simulations with 1324 randomly selected helices of this pool.
Sliced helix peptide library
Expanding the pool of potential ligands we created a second library. It is also based on helix extraction from the proteins of the PDB, but applying less exclusionary criteria. All helical elements of any length were extracted and then chopped to smaller pieces of desired length. The extraction is a process taking up very much computational capacity to complete , so we couldn’t even calculate the number of resulting peptides, but we estimated that it exceeds 100 million structures by far. For all simulations mentioned later on we only use a random 30,034 of these ligands. This sliced peptide library offers a higher diversity at the cost of increasing the risk that the ligands’ real folding is not predicted accurately by the extracted structure files.
However, the question of folding can easily be approached after docking by performing a longer (50 ns) MD simulation to analyse any changes in structure or alpha helicity pointing to lacking transferability of the modeled .pdb structure to the real folding. Therefore we provide the sliced peptide library with the prerequisite that post-docking analysis of the best potential ligands should be performed to ensure transferability of the results. However, as nature shows, a higher diversity combined with an additional selection step often yields better results.
Molecular docking to select high affinity peptide binders
The central part of our in silico equivalent to mirror-image phage display is the molecular docking of an L-ligand to a D-target using AutoDock Vina for evaluation of the binding affinity. This program is based on a scoring algorithm, aiming to find the global minimum of the binding energy of two given molecules, considering their distance in different poses and the resulting intra- and intermolecular interactions. AutoDock Vina allows to manipulate different parameters such as the exhaustiveness of the search which corresponds to the number of random starting positions of the scoring algorithms search. Additionally one can model the ligand’s properties as either rigid or flexible. Flexible docking allows movement of the side chains making the simulation more representative at the cost of increasing the computation time needed for a veritable result. Several experts on the field whom we contacted on the matter of the accuracy of AutoDock Vina emphasized that most modelling programs come with a considerable error. Therefore, we saw it fit to evaluate our docking protocol using a set of controls from the literature. 10
We modeled the well-characterized interactions between the cancer-related protein MDM2 and several of its small peptide inhibitors, comparing our results with the crystal structures and then experimentally determined binding affinities. As negative control we chose 12mer poly-serine which should not show a significant binding towards MDM2. Since the relevant binding domain of MDM2 is only 86 amino acids long and the ligands are small 12mer peptides this system is ideal for comparison with the PSMα3 system which consists also of small peptides.
Overall the binding free energies ranged from -8 kcal/mol to approximately -12 kcal/mol.
Fig. 10: Evaluation of docking protocols for AutoDock Vina on the example of known ligands for MDM2. S12 = 12mer serine as negative control. rigid = simulation of the ligand as rigid structure; flex = simulation of the ligand allowing flexibility for the side chains.
Fig. 11: Absolute values of the free binding affinities of unsliced helixes docked to MDM2. Docking was performed with AutoDock Vina rigid at exhaustiveness 10.
As one can see the binding energy of our negative control S12 (polyserine) corresponds to the left edge of the histogram of possible binding energies. The positive controls pDI, pMI and 6W however show higher binding affinities in our model. The mediocre binding energy which was predicted for the less-affine p53 shows that our model detects only high-affinity binders accurately.
The evaluation of different combinations of exhaustiveness and ligand states led us to the conclusion that dockings at exhaustiveness 10 with a flexible ligand model gives the best separation, however yields lower absolute values for the free binding energy. The reason for this is not yet completely clear but our simulations point to an intrinsic offset within the algorithm of Vina. Since flexible docking requires significantly more computational capacity we chose rigid docking at exhaustiveness 10 for our system and purposes. Our evaluation shows that also this docking protocol yields a good separation between high affinity binders and the negative control.
In our modelling approach we encountered the problem that it takes up lots of computational capacity to model the ligand in a flexible way for docking, even if this seems to be the more accurate model. This phenomenon is size-dependent, meaning that it might not occur smaller ligands like chemical compounds. Yet for peptide ligands the number of rotatable bonds and consequently the number of different possible conformations seems to be too large. Also for larger target molecules and hence larger search spaces the computational time needed to gain the same accuracy would increase massively.
Fig. 12: Visualization of favorable binding pockets on PSMα3. The unsliced library of 1324 helical peptides was docked to D-PSMα3 and visualized in PyMOL, representing peptides as sticks according to the orientation of their backbone, colored according to affinity (the darker the higher).
Identification of gridboxes for optimizing docking efficiency
Responding to this challenge and to keep our modelling approach universally applicable even to larger systems we created an analysis for the visualization of binding pockets of a target protein, meaning the surface patches of a protein that are accessible to ligand binding.
To determine binding pockets of a protein one run of docking with a random library of ligands is required, for example the “unsliced library”. Based on the binding of those random ligands and their binding the most likely spaces of binding can be visualized. As an example we did this for PSMα3.
With this model one can create custom tailored grid boxes, covering the most probable sites of binding in order to minimize the computational requirements while losing as little potential binding candidates as possible. To draw design gridboxes see our tool on the methods page.
We applied this analysis on PSMα3 but due to the small size of this protein there was only a minor advantage in using the new gridboxes instead of a global one around the whole protein. Therefore we decided not to apply the results of this analysis to our docking.
In silico modelling of MIPD
Fig. 14: Absolute values of the binding free energy of 30,034 sliced and unsliced helices docked against D-PSMα3. Docking was performed with AutoDock Vina rigid at exhaustiveness 10.
Combining our in silico models of all the components needed for MIPD we proceeded to simulate the whole MIPD system. We performed molecular docking of the whole 30,034 L-peptides from the unsliced helix peptide library against our model of D-PSMα3 using AutoDock Vina.
The best binders we identified are listed in the Table below.
Transferring our model back to wet lab reality we synthesized the best candidates of both libraries by SPPS and analyzed their binding affinity to PSMα3 by BLI. We could only detect a minor and quite unstable binding to PSMα3, which is much lower than to the binding affinities which were achieved by wet lab MIPD.
Fig. 15: Binding of the ligands derived from our in silico model of MIPD to D-PSMα3. Real time BLI measurement of the association of each ligand at 50 µM concentration to a D-PSMα3 coated sensor. First a baseline of the D-PSMα3 coated sensor is measured, then the association of the ligands and the resulting shift in refractive index. The last phase shows the dissociation of the ligands. As reference the wetlab MIPD ligand 8 is also depicted, water control in orange.
A possible reason for this could be that the library we used for docking was not large and diverse enough.
In silico optimization of MIPD results using an genetic algorithm
To optimize the binders from the in silico MIPD model and also the wet lab MIPD we implemented the power of evolution at molecular level. Using our software finDr GA we subjected an initial population of the wet lab and in silico MIPD ligands to several generations of selection according to their binding affinity to PSMα3 as determined by AutoDock Vina, recombination with each other and random point mutations of the amino acid sequence. With this protocol, which is described to more detail on our software page, we achieved an improvement of binding affinities over 14 generations of directed evolution in silico.
Fig. 16: Improvement of peptide ligands to PSMα3 in the genetic algorithm. Depicted is the binding energy of the best binder to D-PSMα3 in each generation of the genetic algorithm.
Again we transferred this iteration of in silico modelling of our ligands to the wet lab. We synthesized the best two peptides that were predicted to be good binders by the GA using SPPS and analyzed their binding affinity to PSMα3 with BLI (see our methods page page for the more detailed analysis).
Fig. 17: Binding of the ligands from in silico optimization of MIPD ligands to D-PSMα3. Real time BLI measurement of the association of each ligand at 50 µM concentration to a D-PSMα3 coated sensor. First a baseline of the D-PSMα3 coated sensor is measured, then the association of the ligands and the resulting shift in refractive index. The last phase shows the dissociation of the ligands. As reference the wetlab MIPD ligands are also depicted.
With this final transfer of our in silico modelling predictions we demonstrated that the basis of these predictions, our computational model of D-PSMα3 depicts reality in a truthful manner.
Attributions
We gratefully acknowledge the support by the High Performance and Cloud Computing Group at the Zentrum für Datenverarbeitung of the University of Tübingen, the state of Baden-Württemberg through bwHPC and the German Research Foundation (DFG) through grant no INST 37/935-1 FUGG. Having access to the BinAC high performance computing cluster was of essential help for our modelling, especially for performing the simulations with GROMACS and the docking of large libraries.
References
[1] Hung, L., et al. Structure of an Enantiomeric Protein,D-Monellin at 1.8 Å Resolution (1998). Acta Crystallographica Section D Biological Crystallography 54, 494-500.
[2] Somoza, JR., et al. Two crystal structures of a potently sweet protein. Natural monellin at 2.75 A resolution and single-chain monellin at 1.7 A (1993). J Mol Biol. 234(2):390-404.
[3] Berendsen, H.J.C. et al. GROMACS: A message-passing parallel molecular dynamics implementation (1995). Comp. Phys. Comm. 91: 43-56.
[4] Rockey, W., et al. Rapid Computational Identification of the Targets of Protein Kinase Inhibitors (2005). Journal of Medicinal Chemistry 48, 4138-4152.
[5] Na, S., et al. Thermodynamic integration network approach to ion transport through protein channels: Perspectives and limits. (2018). Journal of Computational Chemistry. 39(30):2539-2550
[6] Cheung, G., et al. Insight into structure-function relationship in phenol-soluble modulins using an alanine screen of the phenol-soluble modulin (PSM) α3 peptide (2014). The FASEB Journal 28, 153-161.
[7] Tayeb-Fligelman, E., et al. The cytotoxic Staphylococcus aureus PSMα3 reveals a cross-α amyloid-like fibril (2017). Science 355, 831-833.
[8] Yao, Z., et al. Use of a Stereochemical Strategy to Probe the Mechanism of Phenol-Soluble Modulin A3 Toxicity (2019). Biophysical Journal 116, 87a.
