Team:TU Darmstadt/Model

Modeling

Introduction

---

In synthetic biology, theoretical models are often used to gain insights, to predict and to improve experiments. In our project we are modifying Virus-like particles (VLPs) by attaching proteins to the surface of the P22 capsid through a linker. The linking is catalyzed using the enzyme Sortase A7M, which is a calcium-independent mutant of the wild type Sortase A from Staphylococcus aureus. We performed modeling to predict the unknown structure of the Sortase A7M to improve the linking between proteins and therefore to optimize the modification efficiency of our platform.
Two different modeling approaches were used to determine the structure of the Sortase A7M. We compared machine learning approaches to traditional, comparative Monte-Carlo based modeling methods. The results were evaluated using an energy-scoring function and molecular dynamics (MD) simulations. The most promising Sortase A7M structures were used to perform a docking simulation to investigate binding with linkers.

Toggle all

Structure determination

In silico modeling and simulation of proteins requires a 3D structure, which can be obtained from the RCSB Protein Data Bank. However, if no 3D structures are annotated, as is the case with Sortase A7M, the structure has to be determined by other means. The structure prediction of Sortase A7M was done using two different approaches.

Deep Learning

Background

Machine Learning is a class of algorithms that aim to determine a function between two datasets. This is commonly done by presenting the algorithm with training data, as well as a scoring function to measure its success at processing the input data. During training a feedback loop is used to allow the algorithm to automatically find a function that fits the data. In contrast, classical algorithms are often hardcoded to solve a specific problem and only allow for limited flexibility.

A neural network consists of neurons, which are commonly referred to as nodes. They process input using weights, which are adjusted during its training. Nodes in neural networks are linked together: one neuron processes the input of other neurons, loosely mimicking the structure of biological brains. While one usually has a fixed amount of input and output neurons, limited by the data one wishes to classify, adding layers of hidden neurons can improve the classification. This is often referred to as deep learning and has led to revolutions in applications like speech and image recognition.

Using Machine Learning to predict protein structures has many advantages, compared to conventional methods, especially for iGEM teams who often only have limited access to resources. After training a neural network, which is a computationally expensive process and is often done in centralized data centers, it can be used to predict the structure of a wide variety of proteins. Using pretrained models, novel protein structures can be obtained within seconds compared to conventional methods taking several hours or days.

Until earlier this year the use of Machine Learning in the prediction of protein structures has been restricted to applications within human-written algorithms. AlQuarishi demonstrated a complete deep learning approach that is able to make predictions within 1-2 Å of other approaches, while only using a fraction of the computational power. This enables accurate structural prediction with less powerful, as well as less expensive, hardware and thus significantly reduces the cost of structural modeling.

Procedure

We used AlQuarashi’s approach in combination with his pretrained model, which was trained on the ProteinNet database containing all structures released prior to the start of CASP12 (12th  Critical Assessment of Techniques for Protein Structure Prediction – 2016) our generated candidate structure CASP12 is named after. The results were tested against the CASP12 datasets and reached distance root-mean-square deviation (RMSD) values between 10 and 13 Å. The RMSD is defined as root-mean-square deviation of all atom positions compared to a template structure. It is defined as:

where x_i is a vector of the atomic coordinates of the i-th atom. All proteins in the CASP datasets were not published until after the competition, and thus represent an assessment with only little bias. We used these pretrained datasets to make structural predictions for our Sortase A7M. The predicted structure was then relaxed in a molecular dynamics simulation.

In the following, the specific steps for obtaining a tertiary structure predicted by AlQuarashi’s model are listed.

1. We used the amino acid sequence of the Sortase A7M in the FASTA format to predict the tertiary structure of the amino acid backbone using AlQuarishi’s implementation of his end-to-end differentiable learning of protein structure with the pretrained preCASP ProteinNet database. The output file was a .tertiary file which contains a sequential 3x3 Matrix with atomic coordinates from each amino acid backbone starting at the N-Terminus. The raw backbone is depicted in Animation 1.
2. As the standard format for protein structure information is the PDB (Protein Data Bank) file format, we wrote a python script to combine the structural information from the FASTA and .tertiary files into a single PDB file. For ease of use we used the Biotite Python Module.
3. Using Rosetta's fixed backbone design program 'fixbb' with the 'hpatch', the optimal position of the side-chains was determined and added to the PDB file. The fixed backbone tool adds the corresponding side-chains and optimizes their conformation. The hpatch database ensures that hydrophilic side-chains are to be preferred on the surface of the protein as our sortase is present in an aqueous environment. The resulting structure is depicted in Animation 2.

In order to evaluate the structure obtained, we constructed a Ramachandran plot by calculating the dihedral angles of the amino acid backbone, as depicted in Fig. 1. These can then be compared to the typical dihedral angles for specific secondary structures such as α-Helices and β-Sheets. The typical angles from a randomly sampled dataset are depicted in Fig. 2.

Figure 1: The dihedral angles of amino acids can be calculated to create a Ramachandran plot.

Figure 2: The dihedral angles over a range of randomly sampled proteins.

Results

Animation 1: The raw PDB File converted from the .tertiary file.

Animation 2: The PDB-File after Step 3.

For analysis the Strucure was viewed in PyMOL . As can be seen in Animation 3 below, no secondary structures could be recognized by PyMOL. Thus, a Ramachandran plot was used to evaluate the dihedral angles of the backbone. This is depicted in Fig. 3. It was found that the angles do not match with the typical angles for α-helices and β-sheets.

Animation 3: The cartoon view in PyMOL.

Figure 3: Ramachandran plot of the predicted structure.

During training the predictions in AlQuarashi’s Model were optimized for their RMSD, which is the root-mean-square deviation of the distance between the atoms of the prediction and reference structure. Thus, even though the predictions are expected to have a similar shape to the physical structure, they may not be in the energy minimum.

Rosetta Comparative Modeling

Background

In our second approach we used the RosettaCommons comparative modeling (RosettaCM), which is based on homology modeling . Homology modeling is a protein modeling method, which requires one or more template structures as a base for the protein to be modeled. The protein sequences are aligned with the sequence of the target protein. Unaligned sections are modeled using fragment or protein libraries, which leads to creating protein structures based on different sequence homologues of the protein of interest. Ab-initio or de novo modeling on the other hand attempts to find protein structures solely based on physicochemical principles applied to the primary sequence, which can be compared to the refolding of a denatured protein.

RosettaCM combines ab-initio modeling with homology modeling. The homologous structures, for which a resolved 3D structure with sufficiently similar sequence exists, are generated using homology modeling. Afterwards the unaligned sequences are modeled de novo. By combining the two methods, RosettaCM represents a precise and resource efficient tool for protein structure prediction. Rosetta applications rely on the Monte-Carlo Optimization, which is a probabilistic approach to finding a local minimum in the energy landscape of protein conformations. The underlying equation serving as the foundation of the statistical Monte-Carlo method is the Metropolis acceptance criterion:

where k_B is the Boltzmann constant, ΔE the difference in energy of the two states and T the temperature. The term k_BT can also be written as a single factor β.

During the statistical protein folding based on the Monte-Carlo method, the initial structure is changed by small random perturbations of the atom locations. Whether the structure is accepted, or not, is decided by the Metropolis acceptance criterion. If ΔE < 0, the structure is accepted, otherwise the newly proposed structure is accepted with probability p as described in the Metropolis acceptance criterion.

Procedure

The RosettaCM protocol requires evolutionary related structures and sequences, as well as fragment files of the target structure. The fragment files serve as a structure template for the proteins and consist of peptide fragments of sizes 3 and 9 amino acids. We gathered five evolutionary related structures from the PDB with the accession numbers:

• 1ija
• 1itw
• 1itp
• 1ito
• 2mlm

The five PDB entries represent different structures of sortases from Staphylococcus aureus. Fragment files can be created with the Robetta online server or with the Rosetta FragmentPicker application.

The RosettaCM procedure is best described in the following steps:

1. sequence and structural alignment of templates
2. fragment insertion in unaligned sections
3. replacement of random segment with segment from a different template structure
4. energy minimization
5. all-atom optimization

The alignment can be performed with various tools. We used MAFFT to generate the multiple sequence alignments. Prior to using the alignments as an input, they were converted to the grishin alignment format, as RosettaCM requires the alignments to be in said format. The minimization is performed using the Rosetta centroid energy function. For the centroid function to be applied, the protein is converted to the centroid representation. A protein in centroid representation consists of the backbone atoms N, C_α, O_Carbonyl and an atom of varying size representing the side chain. The advantage of using the centroid representation is that the energy landscape can be traversed easier due to the smoother nature of the centroid energy landscape. Finally, the generated structure undergoes a second minimization in an all-atom model by means of Monte-Carlo optimization. This is similar to the energy minimization, but without the amino acids being represented as centroids of their functional groups. Structures computed through all-atom optimizations can reach atomic resolutions, which is crucial for a model meant to be used to estimate atomic interactions.

Results

The run yielded 15,000 structures, which have been compared using the Rosetta scoring functions (talaris2013). From the 15,000 structures generated, we inspected the ten best scoring structures.

As can be seen in Fig. 4, the most prominent differences can be found in the regions close to the N- and C-terminus. As fluctuations in those regions are not untypical, we decided to use the best scoring structure, candidate S_14771 (Animation 4), as the input for the simulations to follow.

Figure 4: The structural alignment of the ten best scoring sortase structures displaying minor differences with the exception of the C- and N-terminal regions. N- and C-terminal regions tend to show strong fluctuations, thus it is unsurprising to find the terminal regions to be unaligned.

Animation 4: Sortase A7M candidate S_14771 created through RosettaCM.

In order to evaluate the secondary structure of the Sortase A7M candidate S_14771 a Ramachandran plot was created and compared to the five sortases used as input for the comparative modeling. Comparisons were also drawn with the sortase predicted by Deep Learning, as well as a database of randomly sampled proteins. Ramachandran plots of dihedral angles (Fig. 5) can be a first indicator whether the structures computed are valid.

Figure 5: The Ramachandran plot of randomly sampled proteins and the input structures of the comparative modeling show similar secondary structures. Secondary structure analysis of both sortase candidates reveals absence of secondary structures for candidate CASP12. This is not the case with candidate S_14771 as the Ramachandran plot shows all relevant structures.

The Ramachandran plot (Fig. 5) showing α-helices and β-sheets is a strong indicator of a successful structure determination, as those secondary structures are crucial for the functionality of sortases.

Conclusion

We used machine learning methods, as well as Monte-carlo simulations to determine the structure of the mutated transpeptidase Sortase A7M. The machine learning approach using AlQuarishi's Deep Neural Network yielded a structure which seemed to not have any secondary structures. To exclude the possibility of an error in the PyMOL visualization software by Schroedinger, a Ramachandran plot (Fig. 3) was created. The plot shows that no typical secondary structures are present, which is a strong indicator of a failed approach to determine a structure. The approach, using Rosetta Comparative Modeling, yielded 15,000 structures scored with the talaris2013 scoring function. The ten best structures were aligned and exhibited almost identical secondary structures (Fig. 4). The greatest structural differences are present in the N- and C-terminal regions. Since terminal regions tend to fluctuate more strongly than non-terminal segments of the protein, we deemed those fluctuations non-relevant for the proteins functionality.
Being the best scoring candidate, structure S_14771 was analyzed structurally using a Ramachandran plot (Fig. 5). The plot shows all the relevant and typical structures sortases exhibit and serves as an indicator for a successful structure prediction.
In the steps to follow, a molecular dynamics (MD) simulation will be performed on both structures. Even though structure CASP12 does not seem to be a valid structure, refolding processes during a MD simulation might lead to a relaxation of the protein and allow for a promising prediction of the Sortase A7M structure.

Molecular dynamics

Background

The structure predictions made so far were based on statistical methods with physical constraints. The Deep Learning algorithm uses a neural network trained to find a function associating the amino acid sequence and the final 3D positions of the atoms within the protein. On the other hand, predictions were made with Rosetta using the Monte-Carlo Method. Here random movement of individual atoms occurs, and the energy is estimated after each step.

Even though both methods use physical constraints to find plausible protein structures, neither of them actually simulates the behavior of these molecules within a physical force field. Moreover, both methods do not necessarily output fully relaxed protein structures and simulate water implicitly by preferring hydrophilic parts of the proteins to be on the outside. Thus, we conducted a molecular dynamics (MD) simulation to verify the plausibility of our protein structure and allow equilibration. The molecular dynamics simulation provides the opportunity to simulate water as discrete molecules, creating a solvated protein. This step is crucial to validate the structures, as the interaction with water is one of the primary mechanism for protein folding. Since neither candidate CASP12 nor S_14771 have been modeled with explicit water an according MD simulation is imperative, to verify the correctness of the candidates conformation. This of course is much more expensive in terms of computational ressources. As the protein has to be placed in a simulation box and said box is filled with water molecules. This is called solvation and is visualized for candidate S_14771 in Fig. 6.

Figure 6: Sortase A7M in a force field surrounded by discrete water molecules and ions.

We used GROMACS (GROningen MAchine for Chemical Simulations) as the tool for our molecular dynamic simulations. GROMACS solves Newtons equations of motion for individual atoms . While this classical simulation is much more accurate than predictions made by the other methods, approximations are used nonetheless: Forces are cut after a certain radius and the system size is quite small. Additionally, atoms are assumed to be classical particles, which is not the case, as quantum mechanics plays a role in particle-particle interactions. Still, this simulation is very computationally expensive. Therefore, only time periods less than one second could be simulated.

Methods

To perform the molecular dynamics simulations we mostly followed the GROMACS Lysozyme tutorial as it serves our purpose perfectly. We created our simulation box to be of dodecahedral shape and a 0.7 nm distance of the solute to the box borders. We used periodic boundry conditions and a Na⁺ Cl^- concentration of 0.012 mol/L. The main difference of our approach was that we used the CHARMM36 force field instead of the OPLS-AA/L force field and have adjusted our molecular dynamics parameters accordingly. The simulation was performed on a NVIDIA GTX 760 graphics card allowing us to simulate approximately 1 ns per hour.

To analyse the MD simulation we used the Python programming language and the Biotite package as well as GROMACS analysis tools as covar and anaeig. The first analyses are a root-mean-square deviation (RMSD), a root-mean-square fluctuation (RMSF) and a gyration radius analysis. RMSD calculations have been described in the structure prediction section. To compute the RMSF the movement distance of each residue is computed as a root-mean-square over time as:

where v(t)_i is the position of atom i at time t. The radius of gyration is a quantity used to describe the expansion or compression of a particle. The final analysis performed on the MD simulation is called Principle Component Analysis (PCA). By applying PCA to a protein it is possible to gain insights into the relevant vibrational motions and thereby the physical mechanism of the protein.

Results

First indicators

The first possible indicators of a stable protein structure are converging root-mean-square deviation (RMSD), small root-mean-square fluctuation (RMSF) values as well as converging radii of gyration. Using the Python software package and the module Biotite we calculated these quantities and plotted the results for both candidate S_14771 and candidate CASP12.

Figure 7: The RMSD is one of three main indicators of a stable protein structure of the MD simulation of S_14771 over the period of 200,000 ps. As time progressed the RMSD increased with a smaller slope. The value stabilizes at a time of 110,000 ps and fluctuated around the value of 6 Å.

Figure 8: At t = 40,000 ps already the RMSD has arived at a stable value, while at the same time the gyration (Fig. 10) radius decreases over time continuously. This information suggests the protein might be folding and potentially developing secondary structures not present previously.

Figure 9: The prominent fluctuations of the residues from ranges 105 to 115 might indicate a binding site or another form of functional structure. The radius of gyration, just as the RMSD (Fig. 7), stabilizes around a simulation time of of 110,000 ps and converges towards a value of 16.7 Å.

Figure 10: As from t = 40,000 ps the radius of gyration decreases constantly. At the end of the simulation the gyration radius reaches a value of 17 Å. This behavior indicates folding of the protein structure.

Figure 11: The fluctuations (RMSF) of most residues appear insignificant compared to the first, the last residues and the residues close to residue 110 . Typically the N- and C-terminus tend to fluctuate more intensively due to the lack of stabilizing structures. The prominent fluctuations in the range of residue 105 to 115 can indicate a binding site or another form of functional structure.

Figure 12: The prominent fluctuations of the residues from ranges 105 to 115 might indicate a binding site or another form of functional structure. The radius of gyration, just as the RMSD (Fig. 8), stabilizes around a simulation time of of 110,000 ps and converges towards a value of 16.7 Å.

Typical RMSDs and radii of gyration converge towards a value dependent on the size of the protein. Convergence of those quantities can be interpreted as a stable state of the protein structure. As it can be seen in Fig. 7 and Fig. 9 both the RMSD and the radius of gyration stabilize at the same time as the simulation reaches 110,000 ps (110 ns), suggesting a now stabilized structure of candidate S_14771 solvated in water. Another indicator of a functional protein is the RMSF. Instead of being averaged over all atoms, the RMSF is averaged over time with respect to each amino acid. It provides insights in both protein stability and functionality. Fig. 11 reveals the RMSF of residues 105 to 115 to be significantly higher than that of other residues. This hints at the presence of a functional unit along these residues. As commented on in the section describing our structure prediction approaches, the N- and C-terminal regions tend to fluctuate more strongly as a result of the absence of stabilizing structures.

RMSD and gyration of radius calculations of candidate CASP12 (Fig. 8 and Fig. 10) provide evidence of folding. However, the RMSF values show values significantly higher, an effect possibly caused by instability or refolding. Nevertheless, the strongest fluctuations, disregarding the terminal regions, can be seen in the region of residue 105 to 115. This insight consolidates the theory that residues 105 to 115 might be a part of a functional unit.

We were unsure whether candidate CASP12 can be considered a plausible structure and how to interpret the findings concerning the prominent fluctuations. Therefore, we decided to perform a Principle Component Analysis.

Principle Component Analysis

To analyze our system further Principle Component Analysis (PCA) was performed using GROMACS. By applying PCA to a MD simulation of a protein it is possible to extract the most relevant motions of the protein.

Animation 5: A Principle Component Analysis of a less (blue) and a more relevant (red) principal movements showing the most prominent movements of the Cα-chain of candidate S_14771. Both principal components show movement of the β6/β7 loop consisting of residues 105 to 115 towards the active site. Thus we can assume that the closing β6/β7 loop is involved in the reaction mechanism.

Animation 6: The principal movements of candidate CASP appear similar to each other and no strong single movement can be specified. This makes the most relevant (red) and less relevant (blue) principal component indistinguishable from one another. Moreover the active site amino acids do not appear to be in close proximity, which would make a reaction catalyzed by candidate CASP12 impossible.

The results from the Principle Component Analysis of candidate S_14771 (Animation 5) show a movement of the residues 105 to 115 towards the active site, supporting our theory that residues 105 to 115 are important for the reaction mechanism. Since the the most relevant vibrational movement of the sortase (red), is directed towards the active site, it is possible that the β6/β7 loop either closes the binding site of the ligand peptides or even transports one peptide towards the other.

Animation 6 shows the results of the Principle Component Analysis of candidate CASP12. As the RMSF calculations suggested (Fig. 12), the whole protein seems to be moving randomly with no directed movement. In addition the active site residues are spread across the protein confirming our assumption that the protein is not in a stable or plausible conformation.

Conclusion

We gained evidence that at least one of our Sortase A7M models is a valid and stable candidate by performing various methods to analyse the structural stability and validity of our two Sortase A7M candidates. The candidate S_14771 that was generated using RosettaCM appears to be a fitting candidate not only due to successful analyses, but also since the residues of the active site are close enough to each other to catalyze a ligation reaction. Our model created through deep learning excelled only in terms of RMSD and gyration radius calculations. Not only the RMSF and Principle Component Analysis but also the conformation of the active site have proven candidate CASP12 to be of no use for further calculations as it does not portray a valid conformation of Sortase A7M.

Docking

Now that the binding site of the Sortase had been found, the peptide ligand needed to be inserted into the binding site to create a peptide-protein complex. The procedure of choice for the introduction of a ligand into the binding site of a protein is called docking. In the following sections, we will present the protocol and methods we used as well as the results they yielded.

Background

Enzymes are one of the most relevant macromolecules in biology. Their functionality is determined through the way they interact with their ligands. Although enzymes are highly specific concerning the ligands they interact with, similar compounds can often bind to the same enzyme albeit with different affinity. To determine the best possible binding conformation of the protein-ligand complex, we use FlexPepDock, an algorithm provided by the the RosettaCommons software package.

Procedure

The ab-initio FlexPepDock protocol consists of multiple steps and is documented on the RosettaCommons online documentation. We modified the protocol as the one provided did not work with our approach. The modified protocol has the following form:

1. secondary structure determination
2. complex creation
3. FlexPepDock refinement

To determine the secondary structure of the peptide, fragment files (3-, 5- and 9-mers) had to be generated and a PSIPRED secondary structure prediction had to be performed. As the peptides had a sequence length less than 20 amino acids, we were not able to use the online services such as Robetta and the PSIPRED online service. Instead we used the Rosetta FragmentPicker application and the PSIPRED command line tool. The generated structures serve as the input for the refinement protocol.
The generation of the peptide-protein complex can be divided into three steps:

• peptide creation
• peptide relaxation
• coarse complex creation

The peptide structure was created through ab-initio modeling. Initial creation of the peptide was followed by insertion of the peptide into the sortase binding site. This lead to a coarse model of the peptide sortase complex. Here we used insight gained from the molecular dynamics simulation to place the peptide close to the binding site. This operation was performed using Biotite.
In the final step the FlexPepDock refinement protocol is executed and 50,000 complex structures are generated. We used the inputs as described in [13] , written by the authors of the FlexPepDock documentation.
To get a better overview over our data we performed a clustering in python, using the scikit-learn package. We clustered the structures with respect to:

• total score: the total score of the docking provided by the Rosetta scoring function
• interface score: the sum of the energy of the residues in the interfacing region
• reweighted score: a score calculated by double weighting the contribution of the residues in the interfacing region
• root-mean-square deviation: the root-mean-square deviation of the peptides in relation to the structure with the highest score
• peptide direction: the direction the peptide is facing

Here clustering is used to group the docking results and thereby decrease the samlple size. From the 50,000 results we picked the results with the 500 best total scores, the 500 best interface scores and the 500 best reweighted scores. As we aimed to create an unbiased set for clustering, the abscence of duplicates in the set was ensured. We decreased the sample size to 100 groups representing the best scoring structures from the three categories.

Results

For sequences MGGGGPPPPPP(M-polyG), GGGGPPPPPP(polyG) and PPPPPPLPETGG(LPETGG) 50,000 structures have been created and clustered. After the clustering the sample consisted of 100 structures of docked complexes.

Figure 13: The three best scoring structures (total score, interface score, reweighted score) of the LPETGG-tag are shown. Only two results are visible as the best reweighted score candidate is identical to the best interface score candidate. The reacting section of the LPETGG-tag namely glycine is colored yellow as is the active site. The glycin of both ligand peptides is facing the active site.

Analysis of the scores has shown a similar score for all the three dockings. The best scoring results of the LPETGG docking show a tendency of the glycines to face the active site while also being in close proximity to the active site.

Figure 14: The three best scoring structures (total score, interface score, reweighted score) of the poly-g peptide are shown. Only two results are visible as the best reweighted score candidate is identical to the best interface score candidate. Instead of facing the active site (yellow) the reacting glycines (yellow) appear to interact with the β6/β7 loop of the sortase.

Figure 15: The three best scoring structures (total score, interface score, reweighted score) of the M-polyG peptide are shown. Only two results are visible as the best reweighted score candidate is identical to the best interface score candidate. Concerning the M-polyG peptide no uniform directional orientation can be observed. The structure with the best interface score (light blue) is oriendted towards the β6/β7 loop while the structure with the best total/reweighted (dark blue) is oriented towards the β-sheets.

Fig. 13 shows the docking result of the LPETGG peptide to the sortase. The results shown are the best scoring structures of the clustering with respect to the total score, interface score and reweighted score. As the best scoring structure is the same for the total score and the reweighted score only two peptides are shown. This also applies to Fig. 14and Fig. 15. For both results the reacting glycin residues (yellow) are facing the active site. Additionally, the same residues are in close proximity to the active site.

Fig. 14 and Fig. 15 show the docking of the both polyG and M-polyG. While polyG results align well and seem to be interacting with the β6/β7 loop rather than with the active site, this does not seem to be the case for M-polyG. Instead of both structures interacting with the β6/β7 loop or active site one (best interaction score; dark blue) interacts with the β6/β7 loop and the other (best reweighted/total score; light blue-gray) appears to interact with the active site.

Figure 16: The close up of the M-polyG peptide (best total/reweighted score) indicates an interaction of methionine with arginine₁₃₉ and cysteine₁₂₆.

Figure 17: Methionine of the result with the best interface score interacted with the β6/β7 loop rather than the active site. Still the reactive glycine residues appear to be bound to the β6/β7 loop.

As can be seen in Fig. 16 visualizing the result of the the docking simulation (total/reweighted score) suggests an interaction of methionine and two of the active sites namely arginine₁₃₉ and cysteine₁₂₆. Fig. 17 shows the interaction of M-polyG with the β6/β7 loop. The glycines still interact with the β6/β7 loop. Instead of binding above the β6/β7 loop, which is the case for polyG as illustrated in Fig. 14, the interaction seems to be influenced by methionine. By interacting with the residues in the β-helix methionine could potentially hinder binding of glycine to the β6/β7 loop by partial immobilization of the peptide. Overall peptide binding and orientation is less uniform compared polyG without the leading methionine, which could be an indicator of lesser binding affinity of M-PolyG towards the β6/β7 loop.

Conclusion

To computationally investigate binding affinities of the polyG and M-polyG as well as the LPETGG tags we performed docking simulations using the Rosetta FlexPepDock application. We used a modified version of the recommended protocol as the modified version was easier to automate and served our purpose better than the standard protocol. From the calculated scores only, we could not see a difference in binding affinities. Thus, we inspected the best scoring structures regarding the total score, the interface score and the reweighted score using PyMOL. Since the best structures with respect to total score and reweighted score were the same for all simulations, only two structures have been inspected per run. A polyproline tag was appended to all the peptides to simulate the modification of the VLPs with a small peptide.

As expected, the results showed that for LPETGG, the glycines of both peptides oriented towards the active site. This is unsurprising as peptides with the sequence LPXTGG are known to be substrate of the Sortase. It was more surprising to see the polyG tag oriented away from the active site since polyG also is a known substrate of the sortase. Both polyG peptides were facing the β6/β7 loop (residues 105 to 115) uniformly and appeared to be interacting with it. The M-polyG peptides did not show a uniform orientation or interaction scheme. On one hand the visualization of the best result concerning the total and reweighted score has shown interaction of methionine with the cysteine₁₂₆ and arginine₁₃₉, two residues of the active site. On the other hand, the visualization of the best result with respect to the interface score shows the M-polyG facing the mobile β6/β7 loop. In contrast to the polyG peptide lacking the methionine, the M-polyG peptide is pulled down below the β6/β7 loop by the methionine interacting with one of the β-sheets leading to the active site. This is not the case with the polgG results, which lie aligned in one plane with the β6/β7 loop.

Conclusion

For our project it was key to understand and characterize Sortase A7M. As there is no annotated 3D structure for this specific Sortase, an in silico structure determination was performed. This problem was tackled using two different approaches. The Deep Learning approach did not yield a promising model as later analysis also confirmed. However, comparative modeling with Rosetta produced valid structures. We used the best structure, candidate S_14771, for extensive characterization. We evaluated the model with regard to its secondary structure using Ramachandran plots which suggested plausible secondary structures.

Molecular Dynamics simulations were used to investigate stability and dynamic properties of the candidate. The RMSD and radius of gyration stabilized over the course of the simulation, a first indicator of an equilibrated structure. Interestingly, RMSF analysis showed strong fluctuations of residues 105 to 115. We further investigated this by performing Principle Component Analysis. Doing so, we extracted the principle movements of the model. We could observe movement of the β6/β7 loop towards the active site, suggesting the presence of a binding site. Consequently, we performed docking simulations.

FlexPepDock was used to conduct the docking simulations with target peptides. Each run yielded 50,000 structures. In multiple steps we reduced the amount of complexes to 100 clusters with respect to total, reweighted and interface score. We extracted the best scoring complexes and investigated interactions.

For LPETGG we observed a uniform binding to the active site, fullfilling our expectation. Strikingly, polyG appeared to bind to the β6/β7 loop in a uniform manner. As it is known from literature polyG is a functioning ligand of sortase. Supported by literature and our data, we postulate the following mechanism: the β6/β7 loop transports bound polyG towards the active site of Sortase A7M, thereby lowering the activation energy of the linking reaction.

As the theory is neither backed up by nor contradicts experimental data, further research is required.

Acknowledgements

We would like to thank the working group of Prof. Dr. Kay Hamacher. Especially Benjamin Mayer, Maximilian Dombrowsky and Patrick Kunzmann for their generous advice and support.

Furthermore, we would like to thank the LAB3 for providing us with the computing power necessary to execute our Modeling.

References

---

---