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:

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

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.

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

Modeling of the expression of Gene(RCu)

---

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)

---

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).

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.

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

This part is shown in website of human practice