Overview

In our experiment we use engineered yeast cells to absorb and enrich heavy metals such as cooper and cadmium. At first heavy metal ions diffuse into the cell surface from the liquid phase body, and then heavy metal ions are combined with those heavy metal-treated proteins inside the yeast cells.

Summary

Treating heavy metal pollution by means of biosorption is a complicated process. First, it is very meaningful to study the growth of yeast in heavy metal ions solution. Considering that the toxic effects of heavy metal ions on yeast can’t be ignored, we use the matrix inhibition growth model to simulate the growth kinetics of yeast in heavy metal ions solution. Next, we decide to study the process of biological adsorption from the thermodynamic and kinetic point of view. In terms of thermodynamics, we use the basic thermodynamic function to explain the adsorption process, and the conclusions can guide the further optimization of the biosorption. In addition, different static adsorption models are used to simulate the adsorption process, and the conclusions are able to explain the mechanism of the part of the biosorption process. Then we discuss the change of heavy metal ions with time in the process of biosorption from the point of view of dynamics, and compare with the actual measured data.

Yeast growth model

Heavy metal ions inhibits the growth of yeast. In order to describe the kinetics of cell growth accurately，these crucial factors should be taken into account。Unlike the traditional Monod equation, Andrew equation takes the presence of matrix anticompetitive inhibition into consideration.

\[\mu = {\mu _{\max }}\frac{S}{{{K_s} + S + \frac{{{S^2}}}{{{{\rm{K}}_{l,s}}}}}}\]

\(\mu \)—— Specific growth rate，\({{\rm{s}}^{ - 1}}\);
\({\mu _{\max }}\)—— Maximum specific growth rate，\({{\rm{s}}^{ - 1}}\);
S ——The concentration of limiting substrate，\(g/L\);
\({K_s}\)——Saturation constants, its value is equal to the concentration of limiting substrate when the specific growth rate is exactly half the maximum specific growth rate，\(g/L\).

Fig.1 Comparison between Andrew equation and Monod equation

Taking the presence of matrix inhibition into consideration, when the concentration of heavy metal ions is low, the cell growth rate increases with the increase of heavy metal ions concentration and could reach the maximum value. When the heavy metal concentration continues to increase, the cell growth rate decreases. But when there is no matrix inhibition (Monod equation), the cell growth rate increases with the concentration of the matrix until it approaches the maximum value \({\mu _{\max }}\)

Thermodynamics of Adsorption of Heavy Metals

In order to study the ability of yeast to treat the pollution of heavy metal, we quantify the process of biosorption from a thermodynamic point of view. Therefore the conception of separation constant Kc is introduced to measure the equilibrium concentration ratio of intracellular and liquid heavy metal ions.

\[{K_c} = \frac{{{C_i}}}{{{C_l}}}\]

\({K_c}\): equilibrium concentration;
\[{{C_i}}\]：the concentration of intracellular heavy metal ions，mol/L ; \[{{C_l}}\]:the concentration of liquid heavy metal ions，mol/L ;

Bistability is a common phenomenon in single-cell microbes, that two types of cell phenotype coexist. Bistability is very important for many single-cell microbes adapting to environmental changes. Single-cell microbes can choose the appropriate form according to changes of the environment, and bistability is the basis for achieving this change.

the reason for the existence of bistability in single-cell microbes is complex, and is generally thought to be related to the positive feedback of the gene network. We simulate the bistability in single-celled microbes by establishing a simplified gene regulation model.

\(y\) represents the concentration of the bistable substance, \(x\) is the activator of \(y\), which promotes the expression of the \(y\) gene to increase the concentration of \(y\), and \(R\) is the intracellular inhibitor of \(y\), whose concentration is in accordance with Hill Function. Establish the following equation(1)：

\[ \begin{cases} \frac{R}{R_{T}}=\frac{1}{1+{(x/x_{0})}^{n}}; \\ \tau_{x} \frac{dx}{dt}=\beta y-x; \\ \tau_{y}\frac{dy}{dt}=\alpha \frac{1}{1+R/R_{0}}-y; \end{cases} \]

Here \(R\) represents the concentration of the active inhibitory factor, \(R_{T}\) representing the total concentration, \(n\) is the Hill coefficient, \(x_0\) is the concentration at which the activation rate reaches the half, and the generation of \(y\) is described by the Michaelis-Menten equation, \(\tau_x\) 、\(\tau_y\) is the average survive time.

Under normal circumstances, the cells will be in a steady state, then the derivative of time on the equation (1) is zero, we can get equation (2):

\[ y=\frac{1+{(\beta y)}^{n}}{1+RT/R_{0}+{(\beta y)}^{n}} \]

Note that all of these formulations of derivation of Hill function from mass action kinetics assume that the protein has n sites to which ligands can bind. In practice, however, the Hill Coefficient n rarely provides an accurate approximation of the number of ligand binding sites on a protein.[2] We assume that the Hill Coefficient is 2 according to experience.

We can get the following cubic equation form equation(2), we can get equation(3)

\[ \begin{cases} y^3-ay^2+(\rho /{\beta}^2)y-(\alpha /{\beta}^2)=0 \\ \rho=1+R_{T}/R_{0} \\ \end{cases} \]

For the cubic equation, there can be 1, 2, 3 positive real solutions, to achieve the bistability in this model, the function should have two solutions, we assume that the equation form is equation(4):

\[ (y-a)^{2}(y-ka)=y^3-(2+k)ay^2+(1+2k)a^2y-ka^3=0 \]

Comparing equations 3 and 4, we can get the following parametric equation(5):

\[ \begin{cases} \rho=(1+2k)(1+2/k) \\ \alpha \beta=(2+k)^{1.5}k^{-0.5} \\ \end{cases} \]

Draw the parameters equation, can be obtained under the range of parameters in bistability:

Fig.X. parameter equation of \(\alpha \beta \) and \(\rho \)

The two curves are the bistable curves we want, and between the two curves the cells can achieve bistability.

Fig.X.parameter equation of \(\frac{\alpha \beta}{\rho} \) and \(\frac{1}{\rho} \)

The above result is obtained by the equation (1) assuming that the cell is in a steady state, and the stability analysis is given below for equation (1)

Let \(x^{*}\) and \(y^{*}\) represent stability：

\[ x=x^{'}+x^{*};y=y^{'}+y^{*} \]

Combine equation (1), we can get equation(5):

\[ \begin{cases} \tau_{x}\frac{dx^{'}}{dt}=\beta y^{'}-x^{'} \\ \tau_{y}\frac{dy^{'}}{dt}=a x^{'}-y^{'} \\ \end{cases} \]

In this equation:

\[ a=\frac{\partial}{\partial x}(\alpha\frac{1}{1+R/R_{0}})|_{x=x^{*}}=\frac{\partial}{\partial x}(\alpha\frac{1}{{1+(x/x_{0})}^n})|_{x=x^{*}} \]

According to the steady state theory of ordinary differential equations, the equilibrium point is stable if the eigenvalues of the coefficient matrix of equation (5) contain negative real parts, we can find the eigenvalues of equation (5)

\[ \lambda_{1,2}=-1\pm\sqrt{\alpha \beta} \]

When \(\alpha \beta \)<1, the equations reach the steady state. Combine the equation with equation(2):

\[ g(x):=a\frac{\beta(1+x^n)}{\rho+x^n}-x=0 \]

When \(\alpha \beta \)<1, \(g(x)^{'}＜0\), thus we can judge the stability of the resulting roots of equation(2).

Bistability is a common phenomenon in single-cell microbes, that two types of cell phenotype coexist. Bistability is very important for many single-cell microbes adapting to environmental changes. Single-cell microbes can choose the appropriate form according to changes of the environment, and bistability is the basis for achieving this change.

the reason for the existence of bistability in single-cell microbes is complex, and is generally thought to be related to the positive feedback of the gene network. We simulate the bistability in single-celled microbes by establishing a simplified gene regulation model.

\(y\) represents the concentration of the bistable substance, \(x\) is the activator of \(y\), which promotes the expression of the \(y\) gene to increase the concentration of \(y\), and \(R\) is the intracellular inhibitor of \(y\), whose concentration is in accordance with Hill Function. Establish the following equation(1)：

\[ \begin{cases} \frac{R}{R_{T}}=\frac{1}{1+{(x/x_{0})}^{n}}; \\ \tau_{x} \frac{dx}{dt}=\beta y-x; \\ \tau_{y}\frac{dy}{dt}=\alpha \frac{1}{1+R/R_{0}}-y; \end{cases} \]

Here \(R\) represents the concentration of the active inhibitory factor, \(R_{T}\) representing the total concentration, \(n\) is the Hill coefficient, \(x_0\) is the concentration at which the activation rate reaches the half, and the generation of \(y\) is described by the Michaelis-Menten equation, \(\tau_x\) 、\(\tau_y\) is the average survive time.

Under normal circumstances, the cells will be in a steady state, then the derivative of time on the equation (1) is zero, we can get equation (2):

\[ y=\frac{1+{(\beta y)}^{n}}{1+RT/R_{0}+{(\beta y)}^{n}} \]

Note that all of these formulations of derivation of Hill function from mass action kinetics assume that the protein has n sites to which ligands can bind. In practice, however, the Hill Coefficient n rarely provides an accurate approximation of the number of ligand binding sites on a protein.[2] We assume that the Hill Coefficient is 2 according to experience.

We can get the following cubic equation form equation(2), we can get equation(3)

\[ \begin{cases} y^3-ay^2+(\rho /{\beta}^2)y-(\alpha /{\beta}^2)=0 \\ \rho=1+R_{T}/R_{0} \\ \end{cases} \]

For the cubic equation, there can be 1, 2, 3 positive real solutions, to achieve the bistability in this model, the function should have two solutions, we assume that the equation form is equation(4):

\[ (y-a)^{2}(y-ka)=y^3-(2+k)ay^2+(1+2k)a^2y-ka^3=0 \]

Comparing equations 3 and 4, we can get the following parametric equation(5):

\[ \begin{cases} \rho=(1+2k)(1+2/k) \\ \alpha \beta=(2+k)^{1.5}k^{-0.5} \\ \end{cases} \]

Draw the parameters equation, can be obtained under the range of parameters in bistability:

Fig.X. parameter equation of \(\alpha \beta \) and \(\rho \)

The two curves are the bistable curves we want, and between the two curves the cells can achieve bistability.

Fig.X.parameter equation of \(\frac{\alpha \beta}{\rho} \) and \(\frac{1}{\rho} \)

The above result is obtained by the equation (1) assuming that the cell is in a steady state, and the stability analysis is given below for equation (1)

Let \(x^{*}\) and \(y^{*}\) represent stability：

\[ x=x^{'}+x^{*};y=y^{'}+y^{*} \]

Combine equation (1), we can get equation(5):

\[ \begin{cases} \tau_{x}\frac{dx^{'}}{dt}=\beta y^{'}-x^{'} \\ \tau_{y}\frac{dy^{'}}{dt}=a x^{'}-y^{'} \\ \end{cases} \]

In this equation:

\[ a=\frac{\partial}{\partial x}(\alpha\frac{1}{1+R/R_{0}})|_{x=x^{*}}=\frac{\partial}{\partial x}(\alpha\frac{1}{{1+(x/x_{0})}^n})|_{x=x^{*}} \]

According to the steady state theory of ordinary differential equations, the equilibrium point is stable if the eigenvalues of the coefficient matrix of equation (5) contain negative real parts, we can find the eigenvalues of equation (5)

\[ \lambda_{1,2}=-1\pm\sqrt{\alpha \beta} \]

When \(\alpha \beta \)<1, the equations reach the steady state. Combine the equation with equation(2):

\[ g(x):=a\frac{\beta(1+x^n)}{\rho+x^n}-x=0 \]

When \(\alpha \beta \)<1, \(g(x)^{'}＜0\), thus we can judge the stability of the resulting roots of equation(2).

Modeling the functional switch function of yeast cells after mating. After the yeast cells are mated, the expression of thefunctional gene(RCu) is activated, thereby inhibiting the expression of the gene(Cu), and the expression of the downstream gene(Cr) of the inhibited gene(Cu) is activated. As the process involves many genes, and the regulatory process is complex, so we use Cell Illustrator Online to simulate the process.

We build up the model shown in Figure 1, the establishing process shows here.

Fig.X. Simulation model of functional switch

Modeling of the expression of Gene(Cu)

The rate of transcription and translation depends on the concentration of DNA and mRNA, but the transcription and translation process does not directly consume DNA and mRNA, so in these two processes use the auxiliary connections. In addition, we establish the self-degradation process of mRNA and Cu protein

Fig.X. Simulation model of gene(Cu) experssion

Modeling of the expression of Gene(RCu)

Fig.X. Simulation model of gene(RCu) experssion

We add the substance youdao to our model to simulate the expression of the gene RCu. This has two advantages: youdao's self-degradation process makes the inhibitor gradually self-degraded over time, so you can observe whether the effect of RCu protein on Cu DNA is irreversible; youdao can simulate the Influence of the addition of foreign substance to yeast, thus increasing the scope of using of our model.

The rate of transcription and translation depends on the concentration of RCu DNA and RCu mRNA, but the transcription and translation process does not directly consume DNA and mRNA, so in these two processes use the auxiliary connections. In addition, we establish the self-degradation process of RCu mRNA and RCu protein.

Modeling of the binding of protein(RCu) and DNA(Cu)

Fig.X. Simulation model of binding of DNA(Cu) and protein(RCu)

The rate of binding depends on the concentration of DNA（Cu） and protein（RCu）, and this process also consumes DNA（Cu） and protein(RCu). After binding, the product will not express the Cu gene and express the Cr gene. THus achieve functional switch.

Modeling of the expression of Gene(Cd)

The process is similar to the process of the expression of gene(Cu).

Fig.X. Simulation model of Cd expression

Simulation results

The simulation results of the model are shown in the figure above, and the concentration of Cu protein in the initial stage of simulation is still increased due to the presence of mRNA. After the reduction of DNA (Cu) and degradation of mRNA (Cu), Cu protein will gradually disappear. Cr protein is very slow in the beginning due to lack of initial DNA (Cr), and then the Cr protein gradually increases due to the accumulation of DNA (Cr), and then the number gradually becomes stable.

Fig.X. Result of simulation

This part is not common in modeling, because this part is not a model of our biological system to simulate and predict it behavior. What we do in this part is an assessment of its value and scope of use of our biological systems. In addition, we get a preliminary database of copper and cadmium pollution in the world through the data collection, which is also within the scope of modeling.In this part, we work with teammates who responsible for human practice to carry out the main work of data collection and map processing. Because different countries and regions have different monitoring and treatment methods for copper and cadmium pollution, and many countries and regions even lack of monitoring and treatment methods of copper and cadmium pollution, we spent a lot of time to collect data, even so, we did not collect enough data to conduct detailed data analysis of copper and cadmium pollution worldwide. Nevertheless, our existing work can still lead to some preliminary and regional results for the pollution of copper and cadmium.

Data collection

As the United States for the detection and processing of pollution is more perfect, we first from the United States to start the data collection of data. We find the National Aquatic Resource Surveys in the EPA (United States Environmental Protection Agency) website. There are four national assessments in NARS: National Coastal Condition Assessment (NCCA), National Lakes Assessment (NLA), National Rivers and Streams Assessment (NRSA), National Wetland Condition Assessment (NWCA). Only the NWCA have the data of copper and cadmium, so we use it as one of data sources of our maps. The methods and standards used in NWCA are shown in EPA websites.

The methods we used to collect European data is similar to that of the United States. We use the FOREGS-EuroGeoSurveys Geochemical Baseline Database as one of data sources of our maps. The methods and standards used in NWCA are shown in EPA websites.

The data collection process of China is more difficult than that of Europe and the USA, although China also has a comprehensive monitoring and treatment system for copper and cadmium pollution, these data have not published in detail, so we cannot get these authoritative national census data.

So we turn to other sources to collect data of China. The team of Jianchao Li from Shaanxi Normal University retrieved 2450 papers on six kinds of heavy metal pollution (Cd, Pb, Zn, Zn, Cu, Cu, Cr) in China's soil for nearly ten years. After finishing the analysis, the data of 850 papers is plotted the distribution of heavy metal pollution, and provided a detailed distribution of the key pollution areas. The work was published in June 2016 at Bulletin of Environmental Contamination and Toxicology.

Duan Q., Lee J., Liu Y., Chen H., Hu H., (2016) Distribution of Heavy Metal Pollution in Surface Soil Samples in China: A Graphical Review. Bull Environ Contam Toxicol. Jun 24.

For the data of other regions, we did not find similar data sources, so we refer to the practice of the Jiaochao Li team, through the web of science core collection on the literature in recent years, we retrieve 309 sets of data, although these data is not enough Detailed global analysis, we draw some preliminary conclusions.

Map Drawing

We stored data we collected as .mat file in MATLAB, and then we use the worldmap and our map drawing command to draw our World Copper and Cadmium Pollution Map.

Analysis of World Copper and Cadmium Pollution Map

This part is shown in website of human practice