Glutathione Production Membrane System

Glutathione (GSH) is a small molecule made up of three amino acids: glutamic acid, glycine, and cysteine.[1] The sulfur in cysteine and the negatively charged carboxyl groups (when pH is above 4.5) allow GSH to bind with metals and remove them from solutions. It mainly binds with metals like gold (Au), silver (Ag), copper (Cu), mercury (Hg), and cadmium (Cd). This property makes it possible to dissolve these metals.

This study presents a dynamic model for the transport and production of Glutathione (GSH) between two containers, A and B, utilizing selective permeability through membranes with differing pore sizes. The system is designed to explore the kinetics of GSH concentration changes under the influence of production rates and selective flux influenced by membrane pore sizes relative to molecule dimensions [2, 3]. The model incorporates differential equations to describe the concentration of GSH over time, accounting for variable flux rates determined by pore size ratios and potential fouling of the membranes. This setup provides a foundation for understanding the selective transport processes in biochemical systems, which is critical in optimizing conditions for efficient separation and production of small biomolecules such as GSH.

The system consists of two interconnected containers, A and B, each with a volume of 1 litre, linked through membranes with specific pore sizes. Container A is equipped with a membrane with a pore size of 0.65 nanometers, allowing both GSH and smaller amino acids to pass through. In contrast, Container B’s membrane has a smaller pore size of 0.45 nanometers, which selectively permits only amino acids to pass, effectively blocking GSH. The differential permeability between the two membranes creates a controlled environment for studying selective molecular transport based on size exclusion.

To ensure the functionality of the system, the following components are required:

1. Containers (A and B): Each container should be of equal volume (1 litre), constructed from a material compatible with the biochemical substances involved, such as glass or chemically resistant plastic.

2. Membranes: Membranes with specific pore sizes are crucial for the selective transport of GSH and amino acids. Membrane A has a pore size of 0.65 nanometers, while Membrane B has a pore size of 0.45 nanometers. These membranes are designed to provide differential flux rates based on the molecular size of the substances.

3. Pumps: To maintain a controlled flow and pressure environment, pumps are needed to facilitate the movement of fluids through the membranes. These pumps should be calibrated to ensure the proper flow rate that matches the calculated flux rates.

4. Production System for GSH: A mechanism to continuously produce GSH in Container A at a rate of 0.00025 mol/L/min is needed. This can be achieved through an in-situ enzymatic reaction or a chemical synthesis pathway integrated into Container A.

5. Monitoring and Control Systems: Sensors to monitor GSH concentration in real-time in both containers are essential. This can include spectrophotometers or other biochemical sensors capable of detecting GSH concentrations accurately. Additionally, control systems for adjusting flow rates, monitoring fouling rates, and regulating the overall system are necessary for optimal operation.

6. Valves and Piping: Properly sized valves and piping are required to connect the containers, manage fluid flow, and maintain the desired operational pressure.

The system operates by continuously producing GSH in Container A, where it accumulates and partially diffuses into Container B through the selective membranes. The differential equations describe the concentration changes over time, driven by the production rate and the flux rates across the membranes. This setup allows the study of the transport dynamics of small biomolecules and the effect of membrane properties on the selective separation process, providing valuable insights for applications in biochemical engineering and industrial separation processes.

To model the transport and production of Glutathione (GSH) in our system, a set of differential equations that describe how the concentration of GSH changes over time in the two containers, A and B, connected by selective membranes is needed. The concentration of GSH in each container changes over time due to the production rate, flux between containers, and fouling effects of the membranes.

The model is described by the following differential equations:

Model parameters are presented below:

1. Lu, S. C. (2013). Glutathione Synthesis. Biochimica et Biophysica Acta, 1830(5), 3143–3153. Accessed: 2024-09-06. doi: 10.1016/j.bbagen.2012.09.008
2. National Institute for Physiological Sciences. (2013). Discovery of the Mechanism of Cell Death Involving Glutathione in Neurons
3. Kunkel, J.G. Amino Acid and Peptide Buffers

GSH Biosynthesis Simulation

Glutathione (GSH) is a critical tripeptide in biological systems, widely recognized for its ability to bind with metals, including those found in electronic waste (e-waste).

In this study, we model and simulate the biosynthesis of GSH. Using chemical kinetics theory, we describe the deterministic reaction dynamics of GSH biosynthesis, while Gillespie’s algorithm is employed to capture stochastic variations inherent in the process.

The deterministic model, formulated through ordinary differential equations (ODEs), represents the rates of production and consumption of chemical species involved in GSH biosynthesis. Complementing this, the stochastic approach provides deeper insights about the behavior of the chemical model, assuming that the system is in thermal equilibrium and that the volume of the spatial domain is constant.

The deterministic and stochastic approach are implemented using Python script and plots for global approach of modeling.

---

Our main goal here is to study the biosynthesis of GSH, i.e. the chemical reactions that take part to produce GSH. Below are the chemical reactions that we want to model and simulate.

These molecules take part in one or more chemical reactions, meaning that at the same time they can be reactants in one chemical reaction and products in another, as here the γ-GC. We also need to consider the initial conditions such as the concentrations of the chemical species and the reaction constants to study and simulate the rates of consumption and production.

---

- Stochastic Approach

Chemical Master Equation

One way to approach the chemical system is to use the chemical master equation (CME) to simulate the number of molecules at t . The CME is a set of ODEs, one ODE for each possible state of the system. At time t the kth equation gives the probability of the system being in the kth state. The CME provides the full probability distribution of molecular populations in the system.

For the chemical reactions (1) and (2) we have the state vector,

where each element of the above state vector is a non-negative integer that records how many molecules of species of are present. Here N=5 as we have 5 different chemical species. The state vector can change whenever one of the reactions (1) and (2) takes place. Each infinitesimal time interval [t,t+dt] one reaction may affect the state vector

and change this to

where

is the stoichiometric vector.

For the chemical reaction (1), the stoichiometric vector is:

and for the second reaction (2), the stoichiometric vector is:

The stoichiometric vectors describe the change in the number of molecules of each species due to reaction

Let

represent the probability of the system is in state X at time t. The CME provides an equation for the time evolution of the probability distribution P.

The CME expresses how the probability of the system being in state X changes due to the occurrences of various reactions. Mathematically, the CME is written as:

where:

represents the influx of probability into state X due to reaction

transforming state

to state

and the term

represents the outflux of probability from state X due to reaction

The CME is an equation that tracks the change in probability over time. It accounts for all ways the system can enter and leave each state X through individual reactions. The probability fluxes between states are determined by the propensities and the stoichiometric changes associated with each reaction. Solving the CME directly is challenging because it involves tracking the probabilities of all possible states of the system, which leads to a high-dimensional system of coupled differential equations. The CME suffers from several computational challenges, including the curse of dimensionality, nonlinearity, and stiffness, making it impractical for large or complex systems.

Given the complexity of the CME, the SSA algorithm is used as an approximate method to solve the CME.

Stochastic Simulation Algorithm

---

The stochastic simulation algorithm (SSA), also known as Gillespie’s algorithm, gives one solution to this problem and make computational effort more efficient. Gillespie’s algorithm compute realizations of the state vector {t,X} in such a way that the chance of a particular realization being computed reflects the corresponding probability given by the CME.

We consider a well-mixed system with N chemical species and M chemical reactions.

In our case we have 5 different chemical species and 2 chemical reactions. Each reaction transforms some set of reactant species into product species. The state of the system at any given time is characterized by the vector:

For each reaction we define the propensity function which denotes the probability per unit time that respective reaction will occur, given the system is in state X at time t. The propensity is proportional to the number of possible combinations of the reactant molecules and the rate constant c_j

which gives the probability of a reaction occurring per molecule.

After the reaction

occurs, the new state of the system becomes:

The key idea behind the SSA is to determine when the next reaction will occur, and which reaction will occur. For a given state X, the time τ until the next reaction event occurs follows an exponential distribution. The total propensity of the system, denoted as a₀(X)

is the sum of the propensities of all possible reactions.

In our case we have:

The probability that the next reaction occurs in the infinitesimal time interval [t+τ, t+τ+dτ] is given by:

This is the exponential distribution with parameter the sum of the propensities of all possible reactions. To sample the time τ to the next reaction event, we generate a uniform random number r₁ between 0 amd 1 and set:

This gives the time until the next reaction occurs. Once the time τ to the next reaction has been determined, we need to decide which reaction occurs. The probability that reaction

occurs given that any reaction occurs is proportional to its propensity function. Specifically, the probability that reaction is the next to occur is

To select which reaction occurs, we generate another uniform random number r₂ between 0 and 1 and select the reaction

such that:

This ensures that reactions with higher propensities are more likely to be selected.

Tau-leaping

---

The SSA algorithm is very useful at simulating chemical reactions, but it comes with high computational effort because of continuously updating the state vector and the propensity functions. This is crucial overhead that must be addressed better computational solution must be found. As said above, in the SSA, the time to the next reaction is determined by an exponential distribution, and one reaction is executed at a time. However, if many reactions occur rapidly, simulating each one individually, becomes inefficient. The solution to this problem is to use Tau-Leaping. The key idea is to approximate the number of times each reaction occurs during the time interval [t,t+τ], assuming the propensities remain approximately constant within that interval. Instead of simulating reactions one by one, tau-leaping estimates the number of firings or leaps for each reaction.

Let P_j(t) denote the number of times reaction

occurs during the time interval [t,t+τ]. Under the assumption that the propensity functions remain approximately constant during this interval, the number of firings of reaction j follows a Poisson distribution with mean a₀(X)τ.

Thus, the expected number of times the reaction

occurs in the interval [t,t+τ] is a₀(X)τ and the actual number is sampled from its Poisson distribution. To update the system, we calculate the number of times each reaction occurs in the time interval τ using the Poisson distribution. The state of the system after the time step τ is given by:

A critical part of the tau-leaping is to determine an appropriate time step τ such that the propensities a₀(X)τ do not change signigicantly over the interval [t,t+τ]. If τ is too large, the assumption that the propensities remain constant becomes invalid, leading to inaccurate results. The time step τ is usually chosen based on some heuristic that ensures that the number of molecules of any species doesn't change too drastically during the time step. One common criterion is to ensure that the relative change in any propensity function remains small. Mathematically, τ can be chosen such that:

where ε is a small tolerance parameter and v_ji is the stoichiometric coefficient of species i in reaction

In summary, the SSA simulates each reaction individually, advancing time by exponentially distributed intervals. It is computationally expensive when there are frequent reactions. Although, tau-leaping leaps forward in time, allowing multiple reactions to occur in a single time step. It is an approximation but can significantly speed up simulations, especially when many reactions occur rapidly.

---

- Deterministic Approach

Reaction Rate Equation (RRE)

---

At the section above, we explained the stochastic approach of chemical reactions. Although it is very useful to approach deterministically the system with classical chemical kinetics. Here our goal is to study the appropriate system of ODEs that model the biosynthesis of GSH. If we do not consider the stochastic part of the evolution of the system, we get a set of ODEs, according to chemical kinetics, that have the following form:

The set of ODEs for our chemical system is:

where k₁ and k₂ are the rates of reactions.

The key feature of RRE is that:
1. The RRE treat molecular concentrations as continuous variables, assuming that the number of molecules is large enough to justify ignoring discrete molecular fluctuations
2. The dynamics are deterministic, meaning that the same initial conditions will always lead to the same outcome. There is no randomness in the RRE.
3. The RRE describes the average behavior of the system, effectively averaging over many reactions events.

Investigating thermostability of gshF enzyme using computational methods

---

A common field of synthetic biology is protein engineering and especially the study of enzymes thermostability. Thermostability is a significant characteristic of enzymes, that catalyze chemical reaction, because temperatures can affect the enzyme’s structure and by extension its functionality. Here we are investigating, in vitro, the thermostability of the gshF enzyme taking into consideration appropriate factors and calculations. Also, mutations can make significant changes to enzyme’s functionality and give enhanced characteristics to increased temperatures. Especially, in industrial processes, enzymes with enhanced thermostability are more preferred than others because of their rigid structure in extreme environmental conditions. Here we study the thermostability of the wild type gshF enzyme using computational methods and aim to find ways to improve it through mutagenesis approach.

Wild type gshF enzyme

First, it is necessary to gain knowledge about gshF enzyme molecular structure. The pdb format file of gshF enzyme is used to extract useful information about the atoms and by extension the amino acids of the enzyme (residues), the coordinates of the atoms and their respective b-factors.

B-factor is one of the parameters, we took into consideration to explore the stability of the gshF. B-factor is a measure of atomic displacement or flexibility of atoms within a crystal structure. Low b-factor values indicate that the atoms (and by extension the residues) are relatively well-ordered and stable with little movement or positional uncertainty. On the other side, high b-factor values indicate that the atoms or residues are more flexible and dynamic in the crystal structure possibly due to thermal motion, disorder or less interaction with neighboring atoms.

The b-factor of each residue was calculated by averaging the b-factors of all atoms that forms the respective residue. The b-factor for individual atoms within a residue can vary. By examining the b-factors of all atoms in a residue, one can assess the overall flexibility or stability of that residue. It is worth mentioned that b-factors of loops is a significant measurement that helps to assess the flexible regions of the secondary structure. Flexible loops are the hot spots for engineering protein/ enzyme thermostability and that thermostability is greatly correlated to rigidity. The B-Factor for each loop was calculated by averaging the B-Factors of all residues within the loop.

The mathematical equation for the atom’s b-factor is 8π² ⟨ u² ⟩, where ⟨ u² ⟩ is the mean square displacement of the atom from its average position.

Another useful measure we used is the depth of loops. The depth of a loop in the context of protein secondary structure refers to how buried or exposed the loop is within the three-dimensional structure of the enzyme. Loops can either be close to the surface, where they interact with the solvent and environment, or they can be deeply embedded within the protein core, contributing to internal interactions.

The exposed loops are located on the surface of the enzyme and interact directly with the environment or other biomolecules.

They are often more flexible and dynamic. On the other hand, buried loops, are positioned within the interior of the enzyme and are shielded from the solvent. They are often more constrained, interacting with surrounding helices or sheets, and may play roles in maintaining the protein’s structural integrity.

From a computational or structural biology perspective, the depth of a loop or a residue can be quantified using the solvent accessible surface area. This measure the portion of the loop or residue that is accessible to solvent, providing insight into whether a loop or residue is on the surface or buried.

The depth of a loop or residues in gshF enzyme and the thermostability of the protein are closely interrelated. Deep loops are typically surrounded by tightly packed regions of the protein, often interacting with the hydrophobic core. These loops are usually more constrained and less flexible because of their buried position.

By reducing loop flexibility, buried loops prevent unfolding and contribute to the protein's thermostability. Surface-exposed loops are generally more flexible. This flexibility can be essential for the function of certain enzymes (e.g., allowing conformational changes for substrate binding), but it comes at the cost of thermodynamic stability.

At higher temperatures, surface loops are more likely to undergo conformational changes or unfolding, which can destabilize the entire protein structure. This is because thermal fluctuations can disrupt the weaker interactions that stabilize surface regions.

Another useful technique here that would help to investigate the thermostability and flexibility of the wild-type gshF enzyme independently from B-factors is molecular dynamics. Two main calculations interest us here.

The first one is Root Mean Square Fluctuation (RMSF) and the second one is the Root Mean Square Deviation (RMSD).

RMSF measures the average fluctuation (movement) of each atom (or residue) around its average prosition over aperiod of time. The mathematical equation that describes the fluctuation for atom i is:

where T is the time frames, r_i(t) is the position vector of atom i at time t and ⟨ uⁱ ⟩ is the average position vector of atom i over all time frames.

High RMSF values indicate flexible regions of the gshF enzyme while low RMSF values indicate rigid atoms or residues and possibly rigid structure in high thermal conditions.

RMSD quantify the structural deviation of the enzyme and provide insights into the behavior of the enzyme under different conditions. The RMSD is the average distance between the atoms of the enzyme in each conformation and those in a reference structure.

The mathematical equation of RMSD is:

where N denotes the atoms considered for the calculations, r_i the position vector of atom i in the current conformation and r_i^ref is the position vector of atom i in the reference conformation. We consider the reference conformation the conformation at the initial state t=0s.

Low RMSD values indicate that the current enzyme conformation is very similar to the reference structure, meaning that the enzyme is relatively stable and has not undergone significant structural changes, while high RMSD values indicate the opposite.

To calculate those measurements, it is required molecular simulation to assess the enzyme’s structure over time and at specific temperatures. This would help to make a chart to see the RMSD and RMSF values over time for all residues and loops at different temperatures. In this way, it would be helpful to assess the regions that could possibly unfold at a specific temperature, i.e. the most thermally sensitive ones, and determine the mutation points. Because of access limitation to suitable software for molecular simulation, we could not calculate these measurements and further analyze the thermostability of the gshF enzyme.

---

Mutation strategy

Because of access limitation to Rosetta or the relative pyrosetta python module, we could not continue the mutation process and investigation of mutated enzyme’s thermostability. Our goal, after gaining access to suitable software, is to make single point mutations to the wild-type gshF enzyme and calculate the ΔΔG between the mutated enzyme and the wild-type. Negative ΔΔG leads to a more stable enzyme while positive ΔΔG mean that the mutation destabilized the enzyme’s structure and thus the thermostability did not possibly enhanced. It would be useful to make a heatmap that indicates the changes in Gibbs free energy of the mutated enzyme related to the wild type. This can help us to decide which of the mutation can be further tested experimentally for sanity check. One visual example of a heatmap to make it easier to understand is shown below:

Heatmap showing the difference in free Gibbs energy ΔΔG after single point mutations. (source: Higham, D. J. (2008). Modeling and simulating chemical reactions.)

At the example heatmap above, we can observe the ΔΔG of a mutation in the enzyme. Red squares indicate that the Gibbs free energy reduced after the mutation and this possibly leaded to more stable enzyme.

Computational Chemistry for Bioleaching

---

The field of Computational Chemistry would assist our group to predict the interactions between molecular systems, design experiments, and study the effectiveness of bioleaching. Computational chemistry is a field of chemistry that describes the use of computer modelling and simulation – including ab initio approaches based on quantum chemistry, and empirical approaches – to study the structures and properties of molecules and materials. The theoretical chemistry methods are applied to solve complex chemical problems with the use of specific programs and software. In our case, the Gaussian 6.0.16 software and GaussView, the graphical interface used with Gaussian, developed at Carnegie Mellon University, were employed for our purpose.

In computational chemistry, the accurate description of electron orbitals and wavefunctions is important. These orbitals represent the regions in an atom or molecule where electrons are most likely to be found, and wavefunctions provide a mathematical description of the quantum state of a system. The orbitals are designed using basis sets. A basis set in theoretical and computational chemistry is a set of functions (called basis functions) that are combined in linear combinations (generally as part of a quantum chemical calculation) to generate molecular orbitals. For convenience purposes, these functions are typically atomic orbitals centered on atoms but can theoretically be any function; plane waves are frequently used in materials calculations. Our group employed the LANL2DZ basis set, an effective pseudopotential for the core electrons, while the valence electrons were correlated.

Simultaneously, the appropriate functional for the description of the exchange correlation energy (Density Functional Theory (DFT)), is also necessary. The B3LYP functional was opted for our calculations, where B3LYP stands for Becke, 3-parameter, Lee-Yang-Parr.

The B3LYP functional is :

GaussView Software

GaussView Calculation Setup

Gaussian when the calculations have begun

Gaussian was employed mainly for the optimization of the structures, and to obtain the vibrational frequencies of the optimized structure. During the optimization process, the total energy is calculated for different geometries as illustrated in the following diagram.

Optimization Chart

For each complex and interaction, Gibbs free energy was calculated. Gibbs free energy, denoted G, combines enthalpy and entropy into a single value. The change in free energy, ΔG, is equal to the sum of the enthalpy plus the product of the temperature and entropy of the system. ΔG can predict the direction of the chemical reaction. If ΔG is positive, then the reaction is nonspontaneous (i.e., the input of external energy is necessary for the reaction to occur) and if it is negative, then it is spontaneous (occurs without external energy input)

The aim is to study the different interactions between the groups of glutathione and our metals of interest (Au, Ag, Cu) when the solution pH is 11. According to the literature on cysteine leaching, the optimal pH requires the -SH groups to be ionized (-S-)

Therefore, it is essential to calculate the equilibrium constant (K) based on the Gibbs equation:

The larger the constant K, the more favored the products are compared to the reactants.

Complex Cu-1: Interaction between Cu and one sulfur, oxygen atom

Complex Cu-2 : Interaction between Cu and two oxygen atoms

Complex Cu-3: Interaction between CU and one sulfur atom

Complex Cu-4: Interaction between Cu and thiosulfate

For the reaction: Thiosulfate-Cu + GS^- -> GS-Cu + Thiosulfate

ΔG= -51,78 kcal/mol, K= 4 * 10³⁶

Complex Ag-1: Interaction between Ag and two oxygen atoms

Complex Ag-2: Interaction between Ag and one sulfur atom

Complex Ag-3: Interaction between Ag and one oxygen, one sulfur atom

Complex Ag-4: Interaction between Ag and Thiosulfate

For the reaction: Thiosulfate-Ag + GS^- -> GS-Ag + Thiosulfate

ΔG= -35,61 kcal/mol, K= 1 * 10²⁵

Complex Au-1 : Interaction between Au and one oxygen atom

Complex Au-2: Interaction between Au and one sulfur atom

Complex Au-3: Interaction between Au and one sulfur, one oxygen atom

Complex Au-5 : Interaction between Au and Thiosulfate via two sulfur and one oxygen atom

For the reaction: Thiosulfate-Au + GS^- -> GS-Au + Thiosulfate

ΔG= -48,81 kcal/mol, K= 3 * 10³⁴

In addition to the previous research, we explored the effect of the oxidation of gold.

Complex Au-7: Interaction between Au(+3) and GSH in the best geometry available

There was a discussion whether the pH of the leaching solution should be 11 or 7. Therefore, we performed the same calculations with the optimized geometry for each metal.

Glutathione at pH=7, the thiol group is not ionized while the amino group is protonated

Glutathione at pH=11, the thiol group and the amino group are deprotonated

The results and the comparisons are clear in the table below.

Thiosulfate Engineering

Glutathione may interact and stabilize metals ions, but it cannot interact directly with the metals when the oxidation state is 0. For this reason, it is importanto include a leachate that can increase the percentage of oxidized metals. However, all theoretical calculations were performed to clarify that the metal-glutathione complexes are preferred, compared to the metal-thiosulfate complexes. When performing the same calculations with the oxidation state of the Au being 0, the following results are obtained:

Results from the computational chemistry calculations

• Glutathione interacts and forms stabilized metal complexes. Therefore, it is expected to perform metal leaching as we planned.
• Gold and Copper interact better than Silver. These two metals are a big part of electronic devices.
• The optimized pH for the solution of the leaching reaction is indeed 11.
• The increase of the oxidation state of the Au will lead to better results. However, this oxidation is extremely difficult and expensive.
• The interaction between these soft metals and sulfur is stronger than the one with oxygen. A new molecule with more sulfur groups may be even better.

1. "Publisher's note: Sir John A. Pople, 1925-2004". Journal of Computational Chemistry. 25 (9): fmv–vii. 2004
2. Pacific Union College. (n.d.). Gaussian Basis Sets. Chemistry LibreTexts. Retrieved from https://chem.libretexts.org/Courses/Pacific_Union_College/Quantum_Chemistry/11%3A_Computational_Quantum_Chemistry/11.02%3A_Gaussian_Basis_Sets​:contentReference[oaicite:0]{index=0}.
3. Jochen Heyd; Gustavo E. Scuseria; Matthias Ernzerhof (2003). "Hybrid functionals based on a screened Coulomb potential". J. Chem. Phys. 118 (18): 8207
4. Konadu, K. T., & Sasaki, K. (2023). Sulfidic gold ore leaching by cysteine in the presence of Na₂SO₃. Hydrometallurgy, 216, 106018.
5. Xu, B., Kong, W., Li, Q., Yang, Y., Jiang, T., & Liu, X. (2017). A review of thiosulfate leaching of gold: Focus on thiosulfate consumption and gold recovery from pregnant solution. Metals, 7(6), 222.