Team:Lethbridge/Model

CyaNoMore

Modelling
  1. Table of Contents
  2. Background
  3. Modelling Active Site
  4. Materials and Methods
  5. Results
  6. Acknowledgements
  7. References

Background

The modeling team began our focus on predicting mechanisms for different cyanobacterial-specific cargos, such as the cyanophage lysozyme (BBA_K2888000). Initial ideas were also focused on the interaction of antimicrobial peptides and the cyanobacterial membrane of phycobilisome. The lack of available crystal structures, or even sequences of interacting molecules, made set-up for modeling these structures a difficult task.

The modeling team presented some of the initial ideas to Dr. Stacey Wetmore, who is Board of Governors Research Chair in the Department of Chemistry and Biochemistry of the University of Lethbridge whose research is in the field of computational chemistry. Dr. Wetmore discussed ideas with the team for using molecular modeling to improve understanding and characterization of our system components. Dr. Wetmore highlighted the difficulties we may have with simulations that do not have hypothesized active sites or crystal structures.

The discussion with Dr. Wetmore led us to the other primary aspect of our project: the degradation of cyanotoxins. The full sequence of mlrA was available and therefore modeling the structural behaviour at an atomic level would be useful for characterization of this important part and build on previous work in iGEM. Our toxin degradation system involves the microcystinase-A (mlrA) enzyme, which vastly reduces toxicity of microcystin–LR (MC­–LR). Thus, we looked at modeling the enzyme-ligand interaction through molecular dynamics (MD) simulations and docking to improve characterization.

As indicated by Dexter et al. (2020), one of the primary directions of mlrA research is genetic and biochemical characterization. The structural dynamic characterization will be invaluable to understanding of the enzyme for both our team and future teams. Through our computational approach, we can more evenly balance the well-characterized mechanisms of CRISPR-crRNA and its interaction with mRNA with the relatively limited understanding of the mlrA-catalyzed degradation of MC–LRs by building on it with structural dynamic prediction tools. As a consequence, wet lab experimentation can be complemented and informed through our model results. Furthermore, development and implementation of our technology will be improved when we are able to predict the most efficient application of our mlrA system. The model also provides an efficient method to predict behaviour with different structural modifications which will aid in developing the most efficient system.

MlrA: The microcystin degrading enzyme

MlrA is a member of a gene cluster (mlrABCD) required for the full degradation and inactivation of the microcystin cyclic peptide. The critical step in toxin inactivation is carried out by mlrA; and further degradation by the other proteins coded for in the gene cluster. Due to its specificity to toxins and known function, it is a promising candidate for combatting water toxicity due to cyanoblooms via cellular uptake (Maseda et al., 2012).

Fig. 1. Homology model of mlrA based on template 4cad.1.C. Constructed using the SWISS-MODEL engine (Bienert et al., 2017).

MC—LR: The harmful cyanotoxin

MC–LR is what is known as a hepatotoxin: toxins that have been shown to be incredibly harmful by damage the liver tissue of animals. They are heptapeptides, and MC–LR adopts a cyclic structure with non-standard amino acids (Greer et al., 2018). The mechanism of MC–LR interacting with its target protein in epithelial cells, phosphatase 2A, leading to functional inhibition as studied by Liang et al. (2011). Xing et al. (2006), obtained the bound MC–LR to Phosphatase 2A crystal structure, with PDB ID 2IE3, and the averaged minimized NMR structure of MC–LR by itself (PDB ID 1LCM) has been obtained by Trogen et al. (1996).

Fig. 2. MC–LR cyclic hepatotoxin structures, extracted from phosphatase 2A–bound PDB ID 2IE3 (left), and NMR–spectroscopy minimized PDB ID 1LCM (right) MC–LR.

The difficulties associated with cyclic peptides were noted in an interview with James McFarlane, a postdoctoral researcher of biomolecular simulations in Dr. Wetmore’s lab, who indicated that the high variability in cyclic torsion angles of MC-LR would cause challenges in finding the most accurate docking poses. This can be seen in the bound and unbound structures (Fig. 2), where torsions of the cyclic region deviate significantly whether MC–LR is bound to Phosphatase 2A or unbound. Furthermore, minimum–energy conformers of MC–LR may not be the most optimal structure for docking as they do not necessarily represent the bound conformation. This was addressed in our computational methodology in choosing to dock many host–ligand combinations and determine the most accurate binding complexes for subsequent MD simulation. Dr. McFarlane suggested using a conformational search software designed for predicting the unusual conformational behaviours of macrocyclic peptides. The details of these softwares and their methods of calculation relevant to our project are explained more thoroughly in our computational methodology section.

The interaction between host (MlrA) and ligand (MC—LR): previous experimental and computational work

Early studies of the interaction between these two molecules have concluded that the mlrA–B gene cluster enzymes reduce toxicity in a degradation pathway (Bourne et al., 2001).

Fig. 3. Degradation pathway of MC–LR by the mlrA peptidase (Bourne et al, 2001).

The cyclic heptapeptide with the highest degree of toxicity is degraded by mlrA, and further by the mlrBC enzymes. Firstly, mlrA cleaves at the ADDA–Arg peptide bond of cyclic MC–LR (MW = 994 u), resulting in MC–LR linearization (MW = 1012 u) and 160–fold loss of MC–LR reactivity with target phosphatase 2A. Secondly, mlrB converts the linearized MC–LR into the tetra–peptide NH2–ADDA–isoGlu–Ala–OH (MW = 614 u). Finally, mlrC cleavage results in smaller peptides and amino acids.

MlrA contains a metal–binding motif characteristic of metalloproteases (Fig. 4.), HEXXH, thought to interact with MC–LR. Xu et al. (2019) previously studied this interaction, and based on the H260AIH263NE265 motif, thought to be a variant of the HEXXH motif. Their study examined this as the active site for the linearization of the cyclic peptide, but in their studies found that this active site was not possible due to the inability of the Glu265 to act in catalysis with the mlrA H260AIH263NE265 sequence.

Fig. 4. MC–LR previously hypothesized metalloendopeptidase active site. Residues previously thought to be a part of the zinc–binding domain shown, with H263 & H260 as the zinc–binding residues, and N264 as the nucleophilic water hydrogen–bonding residue.

Fig. 5. Mechanism of MC-LR ADDA-Arg cleavage by mlrA proposed by Xu et al.(2019).

Alternatively, Xu et al. (2019) suggests that mlrA is instead a glutamate intramembrane protease of class II CAAX prenyl endopeptidases. The template homology model belonged to a glutamate protease from Methanococcus maripaludis (PDB ID 4CAD_F) categorized as a class II CAAX prenyl endopeptidase. The proposed mechanism involves Glu172 and His205 first activating a water molecules which permits nucleophilic attack of the MC-LR ADDA-Arg peptide bond, with Trp176 and Trp201 potentially raising the pKa for accelerating reaction kinetics, and finally the His260 and Asn264 from the metalloprotease-hypothesized active site serving to stabilize the transition state as an oxyanion hole. Furthermore, their experimental studies which involved comparative EDTA enzymatic inhibition assays also suggest mlrA is likely not a metalloprotease, but a glutamate protease.

The objective statement in the “Modelling the Active Site Interactions with Microcystin-LR Ligand and Microcystinase-A Enzyme” section guided the modeling team in the motivation for choosing different methods to study the enzyme-ligand interaction. We posed different questions and hypotheses that reflect this objective, and informed our choices for the methods used.

Modelling the Active Site Interactions with Microcystin-LR (MC-LR) Ligand and Microcystinase-A (mlrA) Enzyme





The system modelling has been broken up into three parts (A,B, and C), which improve on the dynamic structural understanding of each component for our cyanotoxin-degrading system. Firstly, part A examined the mlrA enzyme in isolation to achieve suitable structures that more accurately represent native folding using Assisted Model Building with Energy Refinement (AMBER) MD simulations (Case et al., 2021). Part B is focussed on both (1) improvements to the MC-LR structure characterization and (2) obtaining the minimum free energy structures for docking by using macrocyclic conformational search tools. Lastly, part C uses the structures obtained in A and B to be used in docking and MD simulations. The general objective of our study is given below in bold, then the research questions, hypothesis, and plan for each part. The Materials and Methods, Results & Discussion follow after.

Objective: To improve the characterization of the mlrA and MC—LR interaction using computational modeling.

With this new approach in mind, we met once again with Dr. Wetmore recommended that we construct the homology model based on the available sequence first and conduct literature research to find the current hypothesized active sites. She suggested the software AutoDock Vina with a grid docking method which would allow us to determine the most stable binding complexes. This way we could improve the understanding of how the mlrA enzyme and MC-LR ligand bind and provide valuable predictions for our experimental work and future iGEM Teams. We then reached out to other computational chemists in her lab, and with advice from Rebecca Jeong, James McFarlane, and Makay Murray we decided to first obtain MD production trajectories based on the homology model. This way, we built on the limitations of the static structure and provide a significantly more accurate look into the conformational dynamics of the enzyme in simulated physiological conditions.

Part A. How can the 3D structural prediction of the mlrA enzyme be improved?

Previous work provides the foundation for understanding enzymatic behaviour, however there are improvements that can be made. Namely, the homology model is static and accuracy of structural prediction as well as dynamic behaviour is greatly improved by using MD simulations. We hypothesize that through using AMBER force field-based MD simulations, the structural behaviour of mlrA can be better characterized. To ensure reliability in prediction and reproducibility of convergence in the conformational ensemble, we chose a biologically relevant time-scale of 500 ns for the MD production runs with a consistent starting minimized, heated, and equilibrated structure ran in triplicate. This allowed us to observe general trends and minimize the effect of statistically irrelevant fluctuations in behaviour.

Our Plan was to:


The SWISS-MODEL Homology model provides an initial crude predicted structure based on the glutamate protease; thus obtaining these approximate atomic coordinates of the enzyme was the first stage. MD allowed us to build on the biases of structure prediction solely using homology modeling by firstly introducing a water box which gives a much better prediction of the different interactions affecting the main and side-chain folding of the protein. The equilibrated model was first obtained, which is a standard first stage in order to obtain input coordinates for MD production simulations of the protein that has been acclimated to solvent conditions by restraining . Next, the production MD is carried out allowing the protein to move without restraints to obtain simulated behavioural data.

Different technical aspects during trajectory setup are highly dependent on the objective. Our rationale for different decisions in methodological setup are described in the Part A computational method section, however the criteria was based on: what the most up-to-date water box is and which ion parameters it is compatible with that allow us to best predict different interactions leading to active site and structural behaviour. Secondly, molecular dynamics is concerned with atomic motion and therefore also simulates the convergence of the enzyme as well as how it arrives to its preferred structure from the homology model starting point. Without crystallographic information of the protein being available, our model provides an ensemble of minimized structures as well as motion of atoms over time calculated with the best computational tools available to capture the structural bases of biochemical behaviour.

Following the MD production simulations, analysis was completed using CPPTRAJ (Roe and Cheatham., 2013). This program reads our production trajectories and outputs analysis data based on analysis scripts that are customizable using CPPTRAJ keywords. Using CPPTRAJ, we completed both visual analysis of enzyme structure as well as different quantitative analyses. Firstly, inter-replicate structure validation was performed due to the highly chaotic nature of simulations containing the motion of thousands of atoms with consistent initial coordinates but randomized initial velocities.

Part B. What is the best structure for MC-LR to be used in docking and MD simulations?

Obtaining the ideal pre-binding conformations of the ligand is difficult due to the unusual structure of MC-LR. First we looked at previous research done computationally, and found that Xu et al. (2019) first obtained the crystal structure and computed restrained electrostatic potential (RESP) by using Gaussian 09 (Frisch et al., 2016) at the HF/6–31G* level of theory to fit electrostatic potentials of MC-LR for non-standard residues. Following this, 100 ns simulations using the general AMBER force field (GAFF) were carried out and clustered for use as ligands in their following docking simulations. We propose that using the QM charge calculating method combined with a conformational search software designed for macrocycles as suggested by James McFarlane will produce the most accurate minimum-energy MC-LR structures for use in docking and MD simulations.

The limitations of Autodock Vina were also stressed to us, as the ligand is constantly rotated and translated and this affects the torsional search space due to the large number of torsion angles in MC-LR. The macrocyclic conformational search softwares are typically used for drug design, but since MC-LR is a macrocycle and a major issue associated with obtaining reliable structures is the cyclic torsions, we determined that this would likely be a better approach than solely using Gaussian QM and MD production simulations.

Since there is discrepancy between the two MC-LR crystal structures (Fig. 2.), we wanted a tool that searches independently of starting coordinates. A study by Poongavanam et al. (2018) assessed different macrocyclic drug conformational search tools and found that OMEGA uses an algorithm that assesses candidate macrocycle structures and the calculations are independent of starting conformation (Hawkins et al., 2010).

The method OMEGA uses is built from the distance geometry (DG) method developed by Spellmeyer et al., 1997, which utilizes an embedding algorithm to generate randomized 4D coordinate generation which are refined into 3D coordinates via a distance geometry error function; this refinement may also include an MD-calculation step. In the ‘Macrocycle’ mode of OMEGAs conformational search tool, the DG error function is replaced with a direct error function, where atomic coordinates are randomly chosen, but comparatively rigid fragments are randomly placed while intact. MMFF94 force field parameters are used in optimization of the direct error function for pairs of atoms in bond angles, torsion angles, Van der Waals radii and tetrahedral constraints. Another advantage offered is a continuum solvation model (Grant et al., 2007) that can be incorporated during refinement to calculate solvent forces, which would allow for better prediction of optimal binding conformations adopted by the free ligand. Thus we have decided on OMEGA as a candidate software, however due to the large cost we were unable to reach this point as of yet.

Part C. How can we better understand the interactions of mlrA and MC-LR at the active site?

Our trajectories of the mlrA enzyme already provide a look into the hypothesis from Xu et al. (2019), as they propose a water coordinates between Glu172 and His205. By examining the interactions between these residues and waters first through radial distribution and followed by AMBER Grid Inhomogeneous Solvation theory (GIST) analysis, not only can we test this hypothesis, but we can also obtain the optimal water geometry for the following docking simulations. From our docking simulations of the best predicted possible structures that are improved from previous homology models, and our understanding of the active site location, the chemical interactions, structural behaviour, and energetics of binding can be better understood. When the coordinates of initial active site orientations are modelled, MD trajectories can be run with the bound ligand-enzyme structure.

Our Plan for Part C would then be:

Materials and Methods

Part A. How can the 3D structural prediction of the mlrA enzyme be improved?

The SWISS-MODEL homology model was first smoothed with the pdb4AMBER tool and tleap was used to generate the coordinate and parameter files for the solvated system. The ff19SB force field was used in conjunction with the OPC water box with 10.0 Angstroms distance from the box edge, which has been shown to improve the accuracy of explicit water simulations. General AMBER parameters were used for the amino acids and TIP4PEW ion parameters compatible with OPC. The TIP4PEW 4-point water box was chosen to emulate the lone pair electrons better and improve approximation of the long-range electrostatics The OPC water model was chosen as it works with many different analysis tools such as GIST and has been shown to improve the accuracy of bulk water properties and works the best with ff19SB (Case et al., 2021). The model was neutralized with 29 sodium ions and 29 chlorine ions as calculated from the Screening Layer Tally by Container Average Potential (SLTCAP) tool (Schmit et al., 2018).

Fig. 5. Mechanism of MC-LR ADDA-Arg cleavage by mlrA proposed by Xu et al. (2019).

Part B. What is the best structure for MC-LR to be used in docking and MD simulations?

Though we haven’t had the opportunity to begin with part B and fully establish the computational method in detail, we have a general scheme for how we will approach the ligand model. Once the receptor model with appropriate water geometry is obtained, as an initial test, we will do a first docking simulation to see if the NMR structure is capable of binding well. Then to improve the ligand structure accuracy, the conformational search engine OMEGA will be used to generate the lowest energy ligand poses to be used in Autodock Vina. Once docking simulations are completed, the following MD simulations on the bound complex require that MC-LR is properly parameterized due to its non-standard residues. Following conformational sampling, computing charges via Gaussian 09 single-point energy calculation would be performed on the candidate lowest-energy structures. The level of theory chosen (HF/6-31G**(2d)) is reflected by the conclusions from Heidari et al. (2019) who found that each of these methods were reliable in replicating FT-IR and NMR data on MC-LR studies. Restrained electrostatic potential (RESP) fitting for MD trajectories would be carried out using antechamber to parameterize modified residues based on the General Amber Force Field 2 (GAFF2).

Results and Discussion

Part A. MlrA: Equilibrated Structure

Fig. 6. Overlayed structures of equilibrated mlrA (blue) and initial homology model (red). From the 10 stages of minimization, heating, and equilibration, we produced the final equilibrated model. This structure is used as the input for subsequent MD trajectory simulations.


Trajectory Motion: Movie, Root-Mean-Squared-Deviation (RMSD) and Radius of Gyration (RoG)

Fig. 7. MlrA backbone RMSD from homology model.

Fig. 8. MlrA backbone RMSD from equilibrated structure.

MlrA backbone RMSD from converged structures. Frames taken from 50ns for replicates 1–3, and 300ns for replicate 2.


The backbone root-mean-squared deviation (RMSD) calculation allows for an initial approximation into the motion of a molecule during a trajectory. The following equation describes the deviation of the set of target coordinates from the reference (e.g. an RMSD of 0.0 indicates that the target coordinates are perfectly overlayed with the reference). N is the number of atoms, the mass of the ith atom is mi, and X and Y are coordinate vectors for the target (X) and reference (Y) ith atoms. The total mass is described by M. Our calculation method utilized a best-fit RMSD which minimized each frame of the RMSD by rotation and translation of the structure to best align with each reference. Multiple references were chosen to have greater confidence with concluding structural convergence. Deviations from the homology model and equilibrated structure were similar among replicates, and the mlrA backbone appears to plateau at similar deviations between replicates. The 50ns frames were chosen as references for the final RMSD calculations, and indicate similar behaviour between replicates. Replicate 2 converged very similarly to 1 and 3 towards the end of the simulation with both 50ns and 300ns coordinates chosen as reference.

Fig. 10. MlrA backbone RoG from center of mass.


The next analysis used was radius of gyration (RoG), which estimates the density/compactness of the structure over the course of the simulation.

Rg is the radius of gyration, N is the number of atoms, r is the radius from the root mean squared center of the kth atom and rmean is the average radius from the center of mass.

The RoG is defined as the root mean squared distance of all the atoms from their center of mass. The folding of the protein can thus be monitored to obtain information about the conformational flexibility of the protein, as well as determining whether or not behaviour between replicates is consistent (Sneha and Doss., 2016). The RoG analysis from the unbound protein can then be compared to the bound to assess the differences in gyration for both models.

Cluster Analysis

The information contained in the MD trajectory includes the movement of the atoms based on 3-dimensional configurations. Clustering is a method used to correlate these configurations via pairwise comparison, where the distance between a set of points is measured. The points are then organized into separate clusters based on their similarities determined by the distances. There are multiple different algorithms available to be used with cluster analysis, and the choice of algorithm depends on selected atoms for dataset pairwise comparisons.



Linkage algorithms are distinct from others such as hierarchical, and tend to categorize clusters into larger groups with small outliers. It is important to know beforehand how many clusters you want to produce with the linkage algorithms. In MD simulations, a cluster centroid is calculated, which requires that the motion of atoms is relative to the same reference frame. These linkage algorithms are categorized as “Bottom-up” clustering algorithms, which utilize an iterative approach. First, the points are assigned to a particular cluster and merged into the larger clusters. Parameters for these can be defined in the input file, such as the epsilon value in CPPTRAJ which defines how large the minimum distance is between clusters.

Fig. 11. Xu et al. (2019) postulated active site where centroids representing clusters (a) and (b) had occupancies of 50% during 3 replicates of 500 ns MD.

Production trajectories were clustered using CPPTRAJ and the above described method with an đœș value of 4.0 Å, the minimum average distance to define a new cluster. The primary catalytic residues appear to not maintain an optimal structure to bridge a water for ~50% of the simulation, but visual inspection does not provide a conclusive result. Otherwise, the active site residues have a somewhat similar geometry as observed by Xu et al. (2019) in both dominant clusters.

Radial Distribution of Water

Fig 12. MlrA model in a water box.

The radial distribution of water serves as a preliminary test of the glutamate-active site. Any two atom masks can be chosen and the distribution function of the occupied distances is generated. This is useful as it allows for the active residue atom distances to any water molecule to be determined. The hydrogen-bonding distance should be maintained for a significant portion of the simulation time, which is ~1.9 Angstroms and can be seen from the output normalized distribution data through peaks on a graph. The radial distribution from the NE2 of His205 to any hydrogens of water molecules, and from OE1 or OE2 to any hydrogens of water molecules was calculated for triplicate production trajectories of mlrA.

Fig. 13. Radial distribution of water from His205 and Glu172

The radial distribution plot reveals that a hydrogen or hydrogens of surrounding water molecules is occupied within the predicted hydrogen-bonding distance for a significant number of frames during the simulation for all interacting residues. The largest peak is seen with OE2 and OE2 of Glu which maintain a near-perfect overlap, whereas the peak is less pronounced with NE2 of His205. As seen in Fig. 11 from the cluster analysis, the rotation of its imidazoline ring and direction of the hydrogen bond acceptor NE2 may be the source of its inability to interact as strongly at Glu172 if a water is coordinated between the residues. Furthermore, the inherent relative hydrophilicity of the residues may result in less hydrogen bonding occurring between a water and His205.

Fig. 14. Radial distribution comparison between catalytic residues (His205 and Glu172)

Following this, the same analysis was carried out on other various Glu and His residues to see if the behaviour was significant. Though it may appear that the behaviour was generally consistent, the His and Glu residues that are placed on the exterior of mlrA and always exposed to solvent (E217, E44, E93, E52 and E265) were comparable in radial distribution to the proposed catalytic residue. Furthermore, E173, which is also within the enzyme cavity like E172, showed a relative decrease in the interaction and the largest deviation from the proposed catalytic Glu173. This provides additional evidence that the water molecule is held in place by Glu172 and His205 considering that their hydrogen bonding interactions were both on-par with residues entirely exposed to solvent, despite having additional steric hindrance preventing water access in the protein core.

Acknowledgements

We would like to acknowledge Dr. Stacey Wetmore, Dr. James, McFarlane, Rebecca Jeong, and Makay Murray from the Wetmore Computational Team for their guidance to the modeling Lethbridge Collegiate iGEM team throughout the 2021 competition.

We would also like to acknowledge Dr. Angeliki Pantazi for sponsoring the modeling team and providing the computational allocation with resources to run our simulations.



References

Bienert S, Waterhouse A, de Beer Tjaart AP, Tauriello G, Studer G, Bordoli L, et al. The SWISS-MODEL Repository—new features and functionality. Nucleic Acids Research. 2017;45(D1):D313-D9. doi: 10.1093/nar/gkw1132.

Bourne D, Riddles P, Jones G, Smith W, Blakeley R. Characterisation of a gene cluster involved in bacterial degradation of the cyanobacterial toxin microcystin LR. Environmental toxicology. 2001;16:523-34. doi: 10.1002/tox.10013.abs.

D.A. Case, H.M. Aktulga, K. Belfon, I.Y. Ben-Shalom, S.R. Brozell, D.S. Cerutti, T.E. Cheatham, III, G.A. Cisneros, V.W.D. Cruzeiro, T.A. Darden, R.E. Duke, G. Giambasu, M.K. Gilson, H. Gohlke, A.W. Goetz, R. Harris, S. Izadi, S.A. Izmailov, C. Jin, K. Kasavajhala, M.C. Kaymak, E. King, A. Kovalenko, T. Kurtzman, T.S. Lee, S. LeGrand, P. Li, C. Lin, J. Liu, T. Luchko, R. Luo, M. Machado, V. Man, M. Manathunga, K.M. Merz, Y. Miao, O. Mikhailovskii, G. Monard, H. Nguyen, K.A. O’Hearn, A. Onufriev, F. Pan, S. Pantano, R. Qi, A. Rahnamoun, D.R. Roe, A. Roitberg, C. Sagui, S. Schott-Verdugo, J. Shen, C.L. Simmerling, N.R. Skrynnikov, J. Smith, J. Swails, R.C. Walker, J. Wang, H. Wei, R.M. Wolf, X. Wu, Y. Xue, D.M. York, S. Zhao, and P.A. Kollman (2021), Amber 2021, University of California, San Francisco.

Dexter J, McCormick A, fu P, Dziga D. Microcystinase - a review of the natural occurrence, heterologous expression, and biotechnological application of MlrA. Water Research. 2020;189:116646. doi: 10.1016/j.watres.2020.116646.

Gaussian 09, M. J. Frisch, G. W. Trucks, H. B. Schlegel, G. E. Scuseria, M. A. Robb, J. R. Cheeseman, G. Scalmani, V. Barone, G. A. Petersson, H. Nakatsuji, X. Li, M. Caricato, A. Marenich, J. Bloino, B. G. Janesko, R. Gomperts, B. Mennucci, H. P. Hratchian, J. V. Ortiz, A. F. Izmaylov, J. L. Sonnenberg, D. Williams-Young, F. Ding, F. Lipparini, F. Egidi, J. Goings, B. Peng, A. Petrone, T. Henderson, D. Ranasinghe, V. G. Zakrzewski, J. Gao, N. Rega, G. Zheng, W. Liang, M. Hada, M. Ehara, K. Toyota, R. Fukuda, J. Hasegawa, M. Ishida, T. Nakajima, Y. Honda, O. Kitao, H. Nakai, T. Vreven, K. Throssell, J. A. Montgomery, Jr., J. E. Peralta, F. Ogliaro, M. Bearpark, J. J. Heyd, E. Brothers, K. N. Kudin, V. N. Staroverov, T. Keith, R. Kobayashi, J. Normand, K. Raghavachari, A. Rendell, J. C. Burant, S. S. Iyengar, J. Tomasi, M. Cossi, J. M. Millam, M. Klene, C. Adamo, R. Cammi, J. W. Ochterski, R. L. Martin, K. Morokuma, O. Farkas, J. B. Foresman, and D. J. Fox, Gaussian, Inc., Wallingford CT, 2016.

Grant JA, Pickup BT, Sykes MJ, Kitchen CA, Nicholls A. A simple formula for dielectric polarisation energies: The Sheffield Solvation Model. Chemical Physics Letters. 2007;441(1):163-6. doi: https://doi.org/10.1016/j.cplett.2007.05.008.

Greer B, Meneely JP, Elliott CT. Uptake and accumulation of Microcystin-LR based on exposure through drinking water: An animal model assessing the human health risk. Scientific Reports. 2018;8(1):4913. doi: 10.1038/s41598-018-23312-7.

Hawkins, P.C.D.; Skillman, A.G.; Warren, G.L.; Ellingson, B.A.; Stahl, M.T. Conformer Generation with OMEGA: Algorithm and Validation Using High Quality Structures from the Protein Databank and the Cambridge Structural Database J. Chem. Inf. Model. 2010, 50, 572-584.

Heidari A, Esposito J, Caissutti A. “Microcystin–LR Time–Resolved Absorption and Resonance FT–IR and Raman Biospectroscopy and Density Functional Theory (DFT) Investigation of Vibronic–Mode Coupling Structure in Vibrational Spectra Analysis”, Malaysian Journal of Chemistry, Vol. 21 (1), 70-95, 2019.

Liang J, Li T, Zhang Y-L, Guo Z-L, Xu L-H. Effect of microcystin-LR on protein phosphatase 2A and its function in human amniotic epithelial cells. J Zhejiang Univ Sci B. 2011;12(12):951-60. doi: 10.1631/jzus.B1100121. PubMed PMID: 22135143.

Maseda H, Shimizu K, Doi Y, Inamori Y, Utsumi M, Sugiura N, et al. MlrA Located in the Inner Membrane Is Essential for Initial Degradation of Microcystin in Sphingopyxis sp. C1. Japanese Journal of Water Treatment Biology. 2012;48(3):99-107. doi: 10.2521/jswtb.48.99.

Morris GM, Huey R, Lindstrom W, Sanner MF, Belew RK, Goodsell DS, et al. AutoDock4 and AutoDockTools4: Automated docking with selective receptor flexibility. Journal of computational chemistry. 2009;30(16):2785-91. doi: 10.1002/jcc.21256. PubMed PMID: 19399780.

Poongavanam V, Danelius E, Peintner S, Alcaraz L, Caron G, Cummings MD, et al. Conformational Sampling of Macrocyclic Drugs in Different Environments: Can We Find the Relevant Conformations? ACS Omega. 2018;3(9):11742-57. doi: 10.1021/acsomega.8b01379.
Schmit JD, Kariyawasam NL, Needham V, Smith PE. SLTCAP: A Simple Method for Calculating the Number of Ions Needed for MD Simulation. Journal of Chemical Theory and Computation. 2018;14(4):1823-7. doi: 10.1021/acs.jctc.7b01254.

Sneha P, Doss CGP. Molecular Dynamics: New Frontier in Personalized Medicine. Advances in protein chemistry and structural biology. 2016;102:181.

Spellmeyer DC, Wong AK, Bower MJ, Blaney JM. Conformational analysis using distance geometry methods. J Mol Graph Model. 1997;15(1):18-36. Epub 1997/02/01. doi: 10.1016/s1093-3263(97)00014-4. PubMed PMID: 9346820.

Trogen G-B, Annila A, Eriksson J, Kontteli M, Meriluoto J, Sethson I, et al. Conformational Studies of Microcystin-LR Using NMR Spectroscopy and Molecular Dynamics Calculations. Biochemistry. 1996;35(10):3197-205. doi: 10.1021/bi952368s.

Trott O, Olson AJ. AutoDock Vina: improving the speed and accuracy of docking with a new scoring function, efficient optimization, and multithreading. Journal of computational chemistry. 2010;31(2):455-61. doi: 10.1002/jcc.21334. PubMed PMID: 19499576.

Xing Y, Xu Y, Chen Y, Jeffrey PD, Chao Y, Lin Z, Li Z, Strack S, Stock JB, Shi Y. Structure of protein phosphatase 2A core enzyme bound to tumor-inducing toxins. Cell. 2006 Oct 20;127(2):341-53. doi: 10.1016/j.cell.2006.09.025. PMID: 17055435.

Xu Q, Fan J, Yan H, Ahmad S, Zhao Z, Yin C, et al. Structural basis of microcystinase activity for biodegrading microcystin-LR. Chemosphere. 2019;236:124281. Epub 2019/07/17. doi: 10.1016/j.chemosphere.2019.07.012. PubMed PMID: 31310980.