Introduction

This year HBUT built engineered Saccharomyces cerevisiae to remove nickel ions in the effluent for the goal of addressing nickel ion pollution problem. In order to make sure that our system could work perfect as expected, we built 4 mathematical models to simulate and predict it with help from HUST-China, HZAU-China, OUC-China.

The first one is pathway’s expression model. By the law of mass action, we first obtained a basic model of the gene expression. We thoroughly predicted the movement process of nickel ions in water during reaction in a three-dimensional homogenous space occupied by Ni²⁺ by making a dispersion model. Thanks for HUST-China, who built these two models for us.

What’s more, we made a mathematical Yeast growth model to demonstrate that TgMTP1t2 protein can alleviate the cytotoxicity of nickel ions with the help of HZAU-China.

OUC-China also helped us made a Molecular dynamic simulation model of the process that nickel ion entering the cell through the NixA transporter protein.

Equilibrium enrichment model. By calculating the equilibrium enrichment and adsorption constant, we can know how long our engineered yeast can be adsorbed in the device and what is the equilibrium enrichment.

With these models, we could improve the experiment significantly and find out the best strategy to make full use of Gluttonous Yeast.

Gene pathway’s expression model

We use mass-action model to describe the process of DNA to Protein. First, there are many copies of DNA in the cell the number of copies is n, we assume that some gene are repressed while others are activated, then the activated DNA can be transcribed into mRNA with a coefficient, then be translated to protein with another coefficient. For mRNA and Protein, they can both degrade with different coefficient. By the law of mass action, we can obtain a basic model of the gene expression.^[1]

We can use ordinary differential equation the describe this.

Since the copy of gene and the activation and repression of gene is constant compared to the transcription process, we can assume that $λ[ X_1]= C_r$, which is the constitutive transcription rate.^[2]

Symbols

	Description	Value	Unit	Source
C_r	Constitutive transcription rate	Differ with different promotor (5.58?)	molecules.s^-1	[1]
X₂	Number of mRNA		molecules
X₃	Number of proteins		molecules
$λ_2$	Transcription rate of gene	0.07	s^-1	[2],[3],[4]
$λ_3$	Translation rate of mRNA	0.2	s^-1	[2],[3],[4]
δ₂	Degradation rate of mRNA	0.005	s^-1	[2],[3],[4]
δ₃	Degradation rate of Protein	0.004	s^-1	[2],[3],[4]

The results is as below. At last it will reach a homeostasis which X_3 reaches to 55804.

Fig 3.Printing

Convection dispersion model

We thoroughly predicted the absorption of Ni²⁺ in a three-dimensional homogenous space occupied by Ni²⁺.

To simplify the model of the real world, we can take only one well-grown cell into our consideration. At first the space was full of Ni²⁺, but in the cell, there is none, with channel protein and hexa-His, Ni²⁺ diffuse into the cell and are combines by the hexa-His protein. Then the Ni²⁺ concentration in the environment decrease gradually.

Using the results of Gene pathway’s expression model, we regard that there are 55804 active molecules on the plasmid.^[3]

For hexa-His gene, we can ignore the time for chelate reaction, and only consider the combination, for the ion channel, we can ignore how Ni²⁺ is absorbed.

Then, according to Stokes-Einstein equation:

$D_AB={\frac{kT}{(6πμ_B R_A )}}$

By diffusion equation:

${\frac{∂N}{∂t}}=D_AB ∇^2 N$

$∇^2$ is Laplace operator, N is a function of x, y, z, t, which we defined to be the number density of Ni²⁺.

We define a cell occupied a globular space V, so the initial state of the system must be:

$$N(x,y,z,0)= \begin{cases} n_0& \text{(x,y,z)∉V}\\\\ 0& \text{(x,y,z)∈V} \end{cases}$$

For the absorption, it depends on the diffusion and we can consider the Ni²⁺ number density in plasmid is always 0 to simplify the model.

$$N(x,y,z,0)= \begin{cases} {\frac{∂N}{∂t}}=D_AB ∇^2 N& \text{(x,y,z)∉V}\\\\ N(x,y,z,t)=0& \text{(x,y,z)∈V}\\\\ N(x,y,z,0)=\begin{cases} n_0& \text{(x,y,z)∉V}\\\\ 0& \text{(x,y,z)∈V} \end{cases}& \text{} \end{cases}$$

These two model parts are made by Dongfang Zhang from HUST-China

Yeast growth model

We hope that the mathematical model can demonstrate that TgMTP1t2 protein plays a role in alleviating the cytotoxicity of nickel ions. The three growth curves of the yeast were then plotted to illustrate the problem.

The growth curve is determined by u_max [maximum growth rate] and x_max [maximum biomass]^[4]

The curve of the change is obtained by determining how Nickel and tgmtp1t2 affect these two parameters.

Logistic curve:^[5]

$μ_max=μ_max^o-K_1 e^{K_2 [N_i^{2+}]}$

$μ={\frac{dX}{dt}}{\frac{1}{X}}$

$ln⁡{\frac{X}{X_0}}=μt$

$μ=μ_{max} (1-{\frac{X}{X_max}})$

$X={\frac{X_0 e^{μ_{max} t}}{(1-{\frac{X_0}{X_{max}}} )(1-e^{μ_{max} t}))}}$

Result analysis: ^[6]

Abscissa: culture time, unit h (hours)

Vertical coordinate: biomass concentration (g/l dry matter)

1. blue line: growth curve at a nickel ion concentration of 0

2. red line: the effect of nickel ions on the growth of yeast without the addition of TgMTP1t2

3. yellow line: the addition of TgMTP1t2 cell biomass rebound, alleviating the toxicity of nickel.

This model part is made by Yun Liu from HZAU-China.

Molecular dynamic simulation model

This year, we are doing the job about the NiCoT transporter, OUC-China are going to make a molecular dynamic simulation of the process that nickel ion entering the cell through the NixA transporter protein.

We all know that NixA is the high-affinity cytoplasmic membrane nickel transport protein of Helicobacter pylori which is capable of importing Ni²⁺into the cell for insertion into the active site of the urease metalloenzyme.^[7]But we also want to use the molecular dynamic method to explain and simulate the process how could the nickel ion go through the cytomembrane with the help of NixA. So OUC-China helped us made a model to explore that question. The result is finally visualized using videos.

At the very beginning, we need to find the protein NixA coordinate file (knowns as PDB files). We search the biggest protein tertiary structure base: RCSB PDB (Protein Data Bank). Unfortunately, we can’t find anything about NixA, even any protein from the NiCoT family. We continue to search the Modbase and obtain the Model-predict structure by the sail lab (Primary Database Link is Q48262.2). However, there is a big hole in the center, which is different from common transporter protein, indicating that this structure is probably not so precise about its true structure.What’s more, the model coverage sketch is from 23 to 275, which is much shorter than the 350 residues we need.

Fig1: NixA model-predict structure from Modbase.

Under that circumstance, we try another way to get the PDB file. That is Swiss Model . We provided the protein sequence to them.

In order to make the simulation more accurate, we use a transporter protein as a template (PDB ID: 2XUT, which demonstrate the crystal structure of a proton dependent oligopeptide (POT) family transporter.) and obtain a well-predict model structure of NixA.

Fig2 SWISS model-predict structure of NIxA.

This model has an eight-transmembrane-domain, which is well-qualified by the characteristics of Nickel/cobalt transporters (NiCoTs)(Fig3). So we take it as our NixA model.

Fig3: Topology of nickel/cobalt transporters (NiCoT)^[8].

However, there is no paper demonstrating where Ni²⁺ connects to the NixA and how to go through it, we have to do it with another simulaion in silico. We making docking between Ni²⁺ and NixA protein using autodock. Due to the fact that autodock couldn’t accept single ion docks with the protein (we also try the FoldX, neither can accompolish satisfying result), so we use the NiCl2 which is very common in the acidic waste water to dock. We get the sdf files of NiCl2 from PubChem and transfer it into pdb format using Openbabel.

Fig4 NiCl2 3Dstructure.

Docking proves to be very difficult because Ni²⁺ is not so common as other metal ion like Cu²⁺, Ca²⁺, so autodock couldn’t recognize its atom type. Finally, we get the parameters from the links below and solve the problem. (http://mgldev.scripps.edu/pipermail/autodock/2009-March/005439.html). Fully unaware of the connecting location, we set the grid box large enough to contain the whole protein, hoping to get a propriate docking conformation. The certain parameter when setting the box are shown in Fig 5. We get 50 docking results of the NiCl2 and NixA and choose the one which NiCl2 is above the NixA protein and close to the transmembrane domains II (Fig6), since requisite motifs for Ni²⁺ transport locates entirely within transmembrane domains II and III.^[7]

Fig5: Parameter we use when setting the docking box.

Fig6: Docking result we choose. The red triangle represents NiCl2 while the green section represents transmembrane domains II of NixA.

For some PDB files, there may be some missing heavy atoms in flexible regions that could not be clearly resolved from the electron density. This may include anything from a few atoms at the end of a sidechain to entire loops or there may be multiple locations listed for some atoms^[9]. So it’s necessary to fix the PDB files before next move. We use PDBfixer to fix the docking result file.

After fixing, we use Membrane Builder tools of CHARMM-GUI to build a protein/membrane complex for molecular dynamics simulations. We use insertion method to insert NixA protein into DPPC^[11] membrane since it’s ubiquitous and a classic membrane model. The result is shown as Fig7(a). We compare it to the topoiogical model of NixA in H. pylori Fig7(b) and find some similarities between each other, further proving the accuracy of the model prediction.

(a)

(b)

Fig 7. Comparison between model and previous paper work. (a)The result of CHARMM-GUI, the membrane is DPPC. (b)Topoiogical model ol NixA in H. pylori.^[10]

Then we put the protein\membrane complex into a water box and add 0.15M K+ and Cl- to make charge balance and systems stable. The specific parameters of the box are as follows.

(a)

(b)

Fig 8. The final result of CHARMM-GUI. (a)The certain parameters of the water box we built in tthis simulation. (b)Visual result of the whole system. Pink dots represents solvent (that is, TIP3), green dots represnts DPPC while red dots in the center of the box represents NixA protein.

We use output of CHARMM-GUI as the input of Further molecule dynamic simulations are all carried out in GROMACS, and we refer to GROMACS Tutorial--Membrane Protein: KALP15 in DPPC to guide us for the rest work. We run the energy minimization, nvt equilibration, npt equilibration and production MD in GROMACS. We use a temperature that is higher than the phase transition temperature of the lipid molecule, which is 323 K. As for tc-grps, we use four groups: protein, DPPC, SOL_CL and NiCl2 to make simulation more accurate.

We use server to run the GROMACS simulations. Due to lake of time, the simulation only lasts for 1.34ns. The group consists of DPPC membrane, protein NixA and NiCl2. The simulation generates about 6G pdb files. Constrained by computing power of PC, such file is far too large for PyMOL to accept. We have to extract first 1000 frames to generate pdb file and use PyMOL to visualize it. The result is shown in Fig 9.

Fig 9. Visualization of the simulation using PyMOL. The green part represents NixA protein, red part represents NiCl2, and yellow part represents DPPC membrane.

We make an 8s video to display the process of Ni²⁺ going through the membrane. As we extract more frames, the video becomes longer, so it can show more parts of the process. We believe we can see finally the Ni²⁺ goes from periplasm side to the cytoplasm side in this way.

This model part is made by Chen Donglin (Lyon, from OUC-China, server supported by Aziya.

Enrichment rate model

The bioenrichment of nickel ions by original yeast is a rapid reaction process. The Lagergren quasi-second-order dynamic model is used to describe:^[12]

$t/q_t=1/(K×q_e^2)+t/q_e$

formula derivation: $dq_t/dt=K(q_e-q_t)^2$

integral：$\int_0^{q_t} {\frac{1}{(q_e-q_t)^2}} {\rm d}{q_t}$ = $K * \int_0^t {\rm d}t$

$q_t=(K×q_e^2×t)/(1+K×q_e×t)$

Take the countdown on both sides at the same time：

$t/q_t=1/(K×q_e^2)+t/q_e$

q_e---Enrichment of Ni²⁺ by yeast in adsorption equilibrium (mg/g)

q_t---Enrichment of Ni²⁺ by yeast at t time (mg/g)

K---Adsorption constant (mg/L)

The adsorption equilibrium qe and adsorption equilibrium constant K were obtained by fitting the enrichment rate model:

Original S. cerevisiae:

K=0.3897 q_e=2.44 (mg/g)

S. cerevisiae/BBa_k3126018(hexa-His):

K=0.1696 q_e=3.663 (mg/g)

S. cerevisiae/BBa_k3126019(nixA):

K=0.2817 q_e=3.3245 (mg/g)

S. cerevisiae/BBa_k3126020(TgMTP1t2):

K=0.4694 q_e=3.201 (mg/g)

S.cerevisiae/BBa_k3126021(hexa-His+nixA)

K=0.298 q_e=3.701 (mg/g)

S. cerevisiae/BBa_k3126022(nixA+TgMTP1t2):

K=0.135 q_e=5.136 (mg/g)

S. cerevisiae/BBa_k3126023（hexa-His+nixA+TgMTP1t2):

K=0.17 q_e=6.349 (mg/g)

By calculating the equilibrium enrichment and adsorption constant, we can know how long our engineered yeast can be adsorbed in the device and what is the equilibrium enrichment.

References:

[1]U. Alon. An Introduction to Systems Biology. Desing Principles of Biological Circuits. Champan and Hall/CRC, Edition, 2007.

[2]Golding I, Paulsson J, Zawilski S M, el al. Real-Time Kinetics of Gene Activity in Individual Bacteria. Cell, 2005, 123:1025-1036.

[3]Shahrezeai V, Ollivier J, Swain P. Colored extrinsic fluctuations and stochastic gene expression. Molecular Systems Biology, 2008, 4:1-9

[4]Oliveira R P D S , Basso L C , Junior A P , et al. Response ofSaccharomyces cerevisiaeto Cadmium and Nickel Stress: The Use of the Sugar Cane Vinasse as a Potential Mitigator[J]. Biological Trace Element Research, 2012, 145(1):71-80.

[5]Najafpour G D . Growth Kinetics and Ethanol Productivity of Saccharomyces cerevisiae PTCC 24860 on Various Carbon Sources[J]. World Applied Sciences Journal, 2009:140-144.

[6]Persans M W , Nieman K , Salt D E . Functional activity and role of cation-efflux family members in Ni hyperaccumulation in Thlaspi goesingense[J]. Proceedings of the National Academy of Sciences, 2001, 98(17):9995-10000.

[7]Fulkerson J F , Mobley H L T . Membrane Topology of the NixA Nickel Transporter of Helicobacter pylori: Two Nickel Transport-Specific Motifs within Transmembrane Helices II and III[J]. Journal of Bacteriology, 2000, 182(6):1722-1730.

[8]Eitinger T , Suhr J , Moore L , et al. Secondary Transporters for Nickel and Cobalt Ions: Theme and Variations[J]. Biometals, 2005, 18(4):399-405.

[9]http://htmlpreview.github.io/?https://raw.github.com/pandegroup/pdbfixer/master/Manual.html PDBfixer manual

[10]Mobley H L T , Garner R M , Bauerfeind P . Helicobacter pylori nickel‐transport gene nixA: synthesis of catalytically active urease in Escherichia coli independent of growth conditions[J]. Molecular Microbiology, 1995, 16(1).

[11]http://www.mdtutorials.com/gmx/membrane_protein/index.html

[12]Zhao Shao. Study on Adsorption of Heavy Metal Ionsby Yeast [D].Shenyang Living University,2016