Abstract

Our project focuses on developing a biofilter system capable of sequestering heavy metal ions from water using engineered biological components. The dry lab team focused on two key proteins: phytochelatin synthase (PCS) and metallothionein (MT), which are known for their high metal-binding potential. We employed computational tools to screen, model, and evaluate candidate sequences based on docking with relevant metal ions. Our pipeline includes sequence selection, homology modeling, docking simulations, and protein engineering to improve selectivity and affinity.

Goal

The goal of our dry lab work was threefold:

1. Engineering: Rationally redesign Metallothionein (MT) using the HSAB principle to achieve selective and strong binding of hard toxic metal ions (Cr, Fe, Al).
2. Screening: Identify and select the optimal Phytochelatin Synthase (PCS) candidate exhibiting the best molecular docking affinity to its substrate.
3. Validation: Use Molecular Dynamics (MD) simulations to structurally and dynamically confirm the success of the engineered MT and check stability of Phytochelatin synthase and Glutathione complex.

Phytochelatin Synthase

Phytochelatin synthase (PCS) is an inducible enzyme that catalyzes the biosynthesis of phytochelatins (PCs) from glutathione (GSH) in response to heavy metal exposure. These oligomeric peptides, with the general structure (γ-Glu-Cys)_n-Gly (n = 2–11), chelate metal ions such as cadmium, lead, and mercury via their thiol groups, forming stable complexes that are sequestered into vacuoles.
Unlike constitutively expressed metallothioneins, PCS is activated upon metal stress, providing a dynamic detoxification response. PCS is found in plants, fungi, and nematodes, and has potential in bioengineering for enhancing metal tolerance or designing biosorbent systems. Synthetic incorporation of phytochelatins into polymer matrices offers a sustainable strategy for removing toxic metals from contaminated environments.

The objective is to reduce the initial dataset of 6,858 phytochelatin syntase (PCS) sequences to a curated set of 5 candidates for downstream analyses or experimentation. The filtering process involved multiple stages of computational screening and biological relevance checks.
This protocol outlines the computational workflow used to identify and select phytochelatin synthase (PCS) sequences based on docking analysis with glutathione, the natural substrate. Sequences were prioritised according to their binding affinity with glutathione, which serves as the substrate for phytochelatin synthesis.

Methodology

1. Database Construction: A database was created containing the following fields for each PCS candidate:
   • Protein name
   • Source organism
   • Accession number
   • URL to protein entry
   • Titles of associated publications or studies
2. Sequence Selection Based on Title Relevance: Sequences were included in the study only if the titles of the associated studies were relevant to the project scope, ensuring biological and functional relevance.
3. Representative Sequence Identification Per Organism: For organisms with multiple PCS sequences, a representative sequence was selected by:
   • Performing multiple sequence alignment in KEGG CLUSTALW of all sequences per organism.
   • Calculating percentage identity.
   • Choosing the sequence with the highest identity representative of that organism’s PCS variants.
4. 3D Structure Prediction: The representative sequences were used to generate three-dimensional protein structures using AlphaFold, enabling structural analysis.
5. First Round Docking (Blind Docking):
   • Glutathione was docked onto each predicted PCS structure without specifying binding sites (blind docking).
   • Binding affinities were recorded and stored in the representative sequence database.
6. Threshold-based Sequence Filtering: Sequences with binding affinities better than or equal to −7.5 kcal/mol from the first round were selected for further analysis.
7. Second Round Docking (Focused Docking with Genetic Algorithm):
   • All PCS sequences from organisms meeting the affinity threshold were dockedagainst glutathione.
   • Docking parameters included rounds of genetic algorithm optimization increased to 50 to refine binding predictions.
8. Final Sequence Selection: Binding affinities from the second round were recorded. The top-performing sequence, based on binding affinity, was selected for experimental validation or downstream analyses.

Flow chart of PCS Selection pipeline Drylab Documentation

Docking protocol for Phytochelatin Synthase

Molecular docking simulations were conducted using AutoDock4 to assess the interaction between glutathione (GSH) and selected phytochelatin synthase (PCS) candidates. The docking procedure consisted of two stages: an initial blind docking screen to filter low- affinity binders, followed by a refined focused docking round

• Ligand Preparation: The structure of glutathione was downloaded from RCSB PDB. Gasteiger charges were assigned, and rotatable bonds were defined via AutoDock- Tools.
• Protein Preparation: For each PCS candidate, the predicted 3D structure from AlphaFold was processed by adding polar hydrogens and assigning kollman charges.
• First-Round Docking (Blind Docking):
  • Objective: Broad screening to eliminate low-affinity PCS candidates.
  • Docking was performed without predefined active sites.
  • Grid box: Large enough to encompass the entire protein structure (blind docking).
  • Genetic Algorithm (GA) runs: 10 .
  • Top-ranked binding affinity (kcal/mol) for each PCS was recorded.
  • Candidates with affinity better than or equal to -7.5 kcal/mol were selected for refinement.
• Second-Round Docking (N-terminal Focused Docking):
  • Objective: N-Terminal focused docking for top candidates.
  • Grid box: Reduced and centered at predicted binding site.
  • Genetic Algorithm runs: Increased to 50.
• The best scoring conformations were selected for each PCS for final ranking. This two-tiered docking protocol enabled high-throughput sequence filtering followed by precise binding affinity estimation for the most promising PCS candidates.

Structural validation of top hits

The 3D structures of the top five hits from our pipeline were generated and analyzed using the DALI server against the full PDB database. All five hits showed significant similarity to experimentally determined phytochelatin synthase structures available in the PDB, confirming the reliability and accuracy of our pipeline. Furthermore, RMSD was calculated by aligning our top hit with the top DALI hit, as reported in the table below, which further validates the structural similarity and demonstrates the robustness of our approach.

Interaction of GSH with with PCS of P.S

Pre MD

Pre MD P.S. PCS - GSH Strucutre Drylab Documentation

Pre MD P.S. PCS - GSH Interaction Drylab Documentation

Pre MD P.S. PCS - GSH Legend Drylab Documentation

Intermediate MD Simulation

Intermediate MD P.S. PCS - GSH Strucutre Drylab Documentation

Intermediate MD P.S. PCS - GSH Interaction Drylab Documentation

Intermediate MD P.S. PCS - GSH Legend Drylab Documentation

Post MD Simulation

Post MD P.S. PCS - GSH Structure Drylab Documentation

Post MD P.S. PCS - GSH Interaction Drylab Documentation

Post MD P.S. PCS - GSH Legend Drylab Documentation

MD Simulation

To assess the stability of our PCS selected from our pipeline, we performed Molecular Dynamics (MD) simulations using YASARA Structure with the AMBER14 force field. The objective was to evaluate whether pipeline selected PCS can interact with Glutathione and maintain structural stability over time.
The PCS-GSH complexe were subjected to MD simulations for 10 nanoseconds under physiological conditions. This allowed us to observe atomic-level dynamics and interactions in a solvated environment, providing insights into the behaviour of the protein–ligand systems beyond static docking models.

Analysis

We performed several structural analyses to characterize the simulations:

• Root Mean Square Deviation (RMSD): to monitor overall protein stability and conformational changes over time.
• Radius of Gyration (Rg): to examine structural compactness during the simulation.
• Inter- and Intramolecular Hydrogen Bonds: to study the strength and stability of protein–metal interactions as well as internal stabilization.
• Approximated Binding Energy (ΔEpot​): Estimated binding affinity qualitatively by analyzing the difference in total potential energy between the final complex and the apo-receptor, calculated as:
  ΔE pot ≈ EComplex ​− EReceptor​

RMSD of P.S. PCS - GSH and P.S. PCS- GSH Drylab Documentation

ROG of P.S. PCS - GSH and P.S. PCS- GSH Drylab Documentation

Intramolecular H bonds P.S. PCS and P.S. PCS- GSH Drylab Documentation

Intermolecular H bonds P.S. PCS and P.S. PCS- GSH Drylab Documentation

Approximate Binding Energy of P.S. PCS with GSH Drylab Documentation

Conclusion

By conducting a DALI search on our predicted PCS, we demonstrated that our top hits are identical to experimentally verified 3D structures of PCS available in the RCSB Protein Data Bank (PDB).
In our molecular dynamics (MD) simulation, the PCS-GSH complex exhibited a rapid decrease in RMSD, followed by a brief increase, and then a decline. In contrast, the RMSD for PCS alone showed significant fluctuations.
The Radius of Gyration for the PCS-GSH complex remained stable and was lower compared to that of PCS alone. Additionally, both intermolecular and intramolecular hydrogen bonds were found to be more numerous and more stable in the PCS-GSH complex than in PCS alone.
We estimated the binding affinity by subtracting the total potential energy of the PCS-GSH complex from that of PCS, which resulted in a negative value. This indicates that glutathione is binding to phytochelatin synthase.

Abstract

Our project focuses on developing a biofilter system capable of sequestering heavy metal ions from water using engineered biological components. The dry lab team focused on two key proteins: phytochelatin synthase (PCS) and metallothionein (MT), which are known for their high metal-binding potential. We employed computational tools to screen, model, and evaluate candidate sequences based on docking with relevant metal ions. Our pipeline includes sequence selection, structural modeling, docking simulations, and protein engineering to improve selectivity and affinity.

Goal

The goal of our dry lab work was threefold:

1. Engineering: Rationally redesign Metallothionein (MT) using the HSAB principle to achieve selective and strong binding of hard toxic metal ions (Cr, Fe, Al).
2. Screening: Identify and select the optimal Phytochelatin Synthase (PCS) candidate exhibiting the best molecular docking affinity to its substrate.
3. Validation: Use Molecular Dynamics (MD) simulations to structurally and dynamically confirm the success of the engineered MT and check stability of Phytochelatin synthase and Glutathione complex.

metallothionein

Metallothioneins (MTs) are a class of low-molecular-weight, cysteine-rich proteins that play a key role in metal homeostasis and detoxification. Their high cysteine content, typically around 20–30 percent, allows them to bind various heavy metal ions such as cadmium (Cd²⁺), mercury (Hg²⁺), zinc (Zn²⁺), and copper (Cu²⁺) through thiolate bonds. These proteins form metal-thiolate clusters, effectively sequestering toxic metals and reducing their bioavailability within the cell.

Methodology

To ensure biological relevance and feasibility of downstream application, a single metallothionein sequence was selected based on literature evidence of its metal-responsive expression profile. Specifically, the metallothionein-like protein 1 from Triticum aestivum (NCBI Accession: NP 001414877.1) was chosen for further computational analysis. A peer-reviewed [7] study associated with this protein reported its upregulation in the presence of aluminum ions, and one of the proteins shared homology with metallothionein suggesting its native ability to respond to and potentially sequester aluminum under abiotic stress conditions. This biologically informed selection enabled a focused and realistic docking study, allowing assessment of its binding affinity toward multiple target metal ions of interest, including Al³⁺, Cr⁶⁺, Fe²⁺, and Fe³⁺. Structural modeling and docking simulations were subsequently performed using this candidate to evaluate its potential for engineering enhanced metal-binding capabilities.

Rationale for Metal Ion Charge Assignment

In molecular docking studies involving metal ions, accurate representation of charge is critical to model electrostatic interactions. However, assigning partial charges requires quantum mechanical (QM) or hybrid QM/MM calculations. These methods, while more accurate, are computationally intensive and not directly supported in standard docking pipelines such as AutoDock4. To maintain computational feasibility and ensure consistency across all docking simulations, formal oxidation states (valency) of metal ions were used as proxies for partial charges. Although this approximation does not capture full electronic behavior, it provides a reasonable and reproducible way to estimate the relative strength of electrostatic interactions with protein residues, especially in a comparative docking framework.
This approach enabled relative comparison of docking affinities while acknowledging the limitations of classical force-field-based docking in modelling metal–protein interactions.

Rationale for Using Individual Metal Ions Instead of Complexes or Hydrated Ions

In our study, individual metal ions were used for docking rather than metal complexes or hydrated ion forms. This choice aligns with the practical goal of designing a biofilter system that incorporates metallothionein and phytochelatin fused within a biopolymer matrix, aimed at remediating contaminated water in mining-affected communities.
People living near mining and industrial sites are often exposed to elevated levels of free metal ions such as Al³⁺, Cr⁶⁺, Fe²⁺, and Fe³⁺ in their water sources. To ensure that our engineered metallothionein is relevant to real environmental conditions and effectively captures the forms of metals most commonly encountered, we performed docking simulations using the individual ionic species. This approach allows us to design and optimize metal-binding proteins with direct applicability to contaminated site conditions.

Docking Protocol for Metallothionein

• Molecular docking simulations were carried out using AutoDock4 to evaluate the interaction between the selected metallothionein and target metal ions (Al³⁺, Cr⁶⁺, Fe²⁺, Fe³⁺). The overall docking workflow consisted of protein structure preparation, metal ion parametrization, grid generation, and docking execution using genetic algorithm (GA)- based sampling.
• Protein Preparation: The amino acid sequence of metallothionein-like protein 1 from Triticum aestivum (Accession: NP 001414877.1) was used to predict the 3D structure via AlphaFold. Polar hydrogens were added, and Kollman charges were assigned using AutoDockTools.
• Ligand (Metal Ion) Preparation: Their formal charges were retained to approximate electrostatic interactions, as partial charge assignment for metal ions would require quantum mechanical methods, which were beyond the scope of this study.
• Grid Generation: A blind docking approach was used. A grid box large enough to encompass the entire protein was generated using AutoGrid.
• Docking Parameters: The Lamarckian Genetic Algorithm was employed with the following settings:
  • Number of GA runs: 50
  • Population size: 150
  • Maximum number of energy evaluations: 2,500,000
  • Maximum generations: 27,000
• Scoring and Ranking: Binding affinities (G, kcal/mol) from the top-ranked pose for each metal ion were recorded. The lowest energy conformation was considered for further analysis.
• This docking protocol was designed to provide a consistent and reproducible framework for evaluating relative metal-binding affinities of the wild-type and engineered metallothionein variants

Molecular Docking of Metallothionein

>NP_001414877.1 metallothionein-like protein 1 [Triticum aestivum]

The initial molecular docking of the Wild-Type Metallothionein (WT-MT) with the target hard acids confirmed the necessity of protein engineering, as the Al³⁺ ion exhibited the least favourable binding affinity at −4.0 kcal/mol. This Al³⁺ score represented the lowest binding affinity and was strategically selected as the primary optimisation target to ensure the robustness of the new binding site.
To maximize the affinity for this hard acid, five custom-designed, carboxylate-rich peptide sequences were generated and subsequently docked with Al³⁺. Only the most successful engineered sequence was then advanced for testing against the remaining hard acid ions (Cr⁶⁺, Fe³⁺, and Fe²⁺).

Design of Engineered Metallothionein

Repurpose metallothionein into a flexible protein that effectively binds Al³⁺, Cr⁶⁺, Fe²⁺, and Fe³⁺

N-Terminal ←− [Cys→Ala MT]+(GGGGS)2+Custom Metal Binding Peptide −→ C-Terminal

Cysteine replacement

Metallothioneins and cysteine rich proteins and cysteine are involved in sequestring metals like lead, mercury, zinc, arsenic. Since our project is focused on sequestering metal ions like Al ³⁺, Cr⁶⁺, Fe²⁺, Fe³⁺, We have mutated Cysteine ( C ) with Alanine ( A ). Alanine is chosen because it is smaller than cysteine and nonpolar amino acids and distorts the thiol based metal binding. ⁵
> Mutated MT [Triticum aestivum]

Chemical Principles Behind the Design

HSAB Theory

According to Pearson’s HSAB theory, binding preference is governed by the matching of acid/base hardness:

• Hard acids (e.g., Al³⁺, Fe³⁺, Cr³⁺): small, highly charged, weakly polarizable → prefer hard bases (O^-, carboxylates).
• Soft acids (e.g., Cd²⁺, Cu²⁺): large, low charge, highly polarizable → prefer soft bases (S^-).

Metallothionein is inherently a soft-base protein, optimized for soft metal ions. To redirect its selectivity toward hard metals, its binding site chemistry must be reprogrammed to incorporate hard-base donors.

Designs of 5 peptides

Peptide A: DGEGDEGDGDEE
Peptide B: EGGEDGEDGGEG
Peptide C: DPGEEPGDDEEG
Peptide D: GEGDEEGGDDEEDG
Peptide E: EDEGEDEGEDEG

Common Design Features Across All Peptides

Glu/Asp : Strong O-donors for hard acid
Gly/Pro : Induce flexibility, turns, and prevent folding
No Cys : Avoids soft ligands that don’t favour hard acid binding

> Engineered MT [Triticum aestivum]

Binding affinities of Engineered Mt with metal of interest

Wildtype Metallothionein vs Engineered Metallothionein

Before Engineering Drylab Documentation

After Engineering Drylab Documentation

Ramachandran Plot

The quality of the protein structure is judged by how many residues fall into the favored, allowed, and outlier regions.

When we compare both the plots, Engineered Metallothionein
Engineered MT has less outliers compared to Wildtype Metallothionein suggesting that engineering improved the secondary structure of metallothionein.

Engineered Metallothionein Drylab Documentation

Ramachandran Plot of Engineered Metallothionein Drylab Documentation

Wildtype Metallothionein Drylab Documentation

Ramachandran Plot of Wildtype Metallothionein Drylab Documentation

Pre and Post MD Simulation Metal and Engineered Metallothionein Interaction

Aluminium (Al³⁺)

Pre MD binding site of Al on Engineered Metallothionein Drylab Documentation

Post MD binding site of Al on Engineered Metallothionein Metallothionein Drylab Documentation

Chromium (Cr⁶⁺)

Pre MD binding site of Cr on Engineered Metallothionein Drylab Documentation

Post MD binding site of Cr on Engineered Metallothionein Metallothionein Drylab Documentation

Iron (Fe^2/3+)

Pre MD binding site of Fe on Engineered Metallothionein Drylab Documentation

Post MD binding site of Fe on Engineered Metallothionein Metallothionein Drylab Documentation

The comparative analysis of pre- and post-MD simulation structures confirms the effective sequestration of the hard metal ions (Al³⁺, Cr⁶⁺, and Fe³⁺) by the Engineered Metallothionein. The metal ions remained tightly bound to the designed binding pocket, forming persistent coordination spheres with the carboxylate oxygen atoms of the engineered peptide. This dynamic stability, demonstrated by the lack of metal dissociation from the site throughout the trajectory, validates the successful shift in binding preference predicted by the HSAB rationale.

MD Simulation

To assess the stability and metal-binding efficiency of our engineered Metallothionein, we performed Molecular Dynamics (MD) simulations using YASARA Structure with the AMBER14 force field. The objective was to evaluate whether the engineered protein can effectively sequester target metal ions — Al³⁺, Cr⁶⁺, and Fe²⁺/Fe³⁺ — and maintain structural stability over time. The protein–metal complexes were subjected to MD simulations for 10 nanoseconds under physiological conditions. This allowed us to observe atomic-level dynamics and interactions in a solvated environment, providing insights into the behaviour of the protein–metal systems beyond static docking models.

Analysis

We performed several structural analyses to characterize the simulations:

• Root Mean Square Deviation (RMSD): to monitor overall protein stability and conformational changes over time.
• Root Mean Square Fluctuation (RMSF): to identify flexible regions and assess local mobility of residues.
• Radius of Gyration (Rg): to examine structural compactness during the simulation.
• Inter- and Intramolecular Hydrogen Bonds: to study the strength and stability of protein–metal interactions as well as internal stabilization.
• Approximated Binding Energy (ΔEpot​): Estimated binding affinity qualitatively by analyzing the difference in total potential energy between the final complex and the apo-receptor, calculated as:
  ΔE pot ≈ EComplex ​− EReceptor​

Eng Metallothionein for Al

PRMSD of WT MT, Al-WT MT, Eng MT, Al-Eng MT Drylab Documentation

ROG of WT MT, Al-WT MT, Eng MT, Al-Eng MT Drylab Documentation

Solute H Bonds of WT MT, Al-WT MT, Eng MT, Al-Eng MT Drylab Documentation

Solute-Solvent H bonds of WT MT, Al-WT MT, Eng MT, Al-Eng MT Drylab Documentation

Approximate Binding Affinity of Al with WT MT and Eng MT Drylab Documentation

Eng Metallothionein for Cr

RMSD of WT MT, Cr-WT MT, Eng MT, Cr-Eng MT Drylab Documentation

ROG of WT MT, Cr-WT MT, Eng MT, Cr-Eng MT Drylab Documentation

Solute H bonds of WT MT, Cr-WT MT, Eng MT, Cr-Eng MT Drylab Documentation

Solute-Solvent H bonds of WT MT, Cr-WT MT, Eng MT, Cr-Eng MT Drylab Documentation

Approximate Binding Affinity of Crwith WT MT and Eng MT Drylab Documentation

Eng Metallothionein for Fe

RMSD of WT MT, Fe-WT MT, Eng MT, Fe-Eng MT Drylab Documentation

ROG of WT MT, Fe-WT MT, Eng MT, Fe-Eng MT Drylab Documentation

Solute H bonds of of WT MT, Fe-WT MT, Eng MT, Fe-Eng MT Drylab Documentation

Solute-Solvent H bonds of of WT MT, Fe-WT MT, Eng MT, Fe-Eng MT Drylab Documentation

Approximate Binding Affinity of Fe with WT MT and Eng MT Drylab Documentation

all mentioned in kJ/mol

Binding affinity was approximated by calculating the change in the system's potential energy (ΔE pot​) upon complex formation, defined as the difference between the final potential energy of the solvated metal-Metallothionein complex and the potential energy of the solvated apo-receptor:
ΔEpot​≈EComplex​−EReceptor​
While this metric serves as a robust qualitative indicator of favorable electrostatic and van der Waals interactions, it is understood to be an approximation, as it neglects entropic and solvation free energy terms (ΔS and ΔGsolv​). However, the massive disparity observed between the Wild-Type and Engineered ΔEpot​ values provides clear evidence of successful Metallothionein engineering.

Conclusion

When comparing the RMSD and Radius of Gyration plots for Wildtype and Engineered Metallothionein, we observe an increase in both parameters for Engineered Metallothionein compared to Wildtype. This can likely be attributed to the substitution of Cysteine with Alanine, which results in the loss of disulfide bonds. Additionally, the introduction of a flexible (GGGS)2 linker between the custom-designed metal-binding peptide and the Cysteine-to-Alanine mutated Wildtype Metallothionein may also contribute to this increase.
Furthermore, an analysis of intramolecular and intermolecular hydrogen bonds reveals an increase in hydrogen bonds in Engineered Metallothionein compared to Wildtype Metallothionein. This increase likely indicates enhanced interactions between metal ions and Metallothionein, as well as within the Metallothionein structure itself.
The MD stability analysis (comparing pre- and post-simulation structures) confirms the effective sequestration of the target hard metal ions (Al³⁺, Cr⁶⁺, and Fe³⁺). The metal ions remained tightly bound to the designed pocket, forming persistent coordination spheres with the carboxylate oxygen atoms of the engineered peptide. This dynamic retention, demonstrated by the lack of metal dissociation from the binding site throughout the trajectory, substantiates the successful shift in binding preference predicted by the Hard-Soft Acid-Base (HSAB) rationale.
Finally, when estimating binding affinity by subtracting the total potential energy of the metal-Metallothionein complex from that of Metallothionein alone, a huge difference is observed when comparing the two in a bar plot, highlighting the disparity between Wildtype Metallothionein and Engineered Metallothionein. Therefore, we conclude Metallothionein Engineering as Successful.