Team:CAU China/Model
Ordinary Differential Equations: Deeper Understanding of Our Project
1. Overview
Before the experiment, we need to prove that the system we designed is feasible through the model, and we hope that the mathematical model we established can let us have a deeper understanding of the behavior and mechanism of our system and project. We designed gene circle to enable Corynebacterium glutamicum to produce γ-PGA, and consider the adjustment of the expression of key metabolic enzymes. In addition, in order to ensure biosafety, we also consider designing a kill switch sensitive to both salt and alkali. We did substantial jobs to describe the behavior of the whole system by establishing ODEs and analyze the sensitivity of a large amount of parameters.2. Background
According to our design, capA, capB and capC from Bacillus subtilis will be introduced into Corynebacterium glutamicum for guidance and catalysis the synthesis of γ-PGA. P0864 promoter was connected upstream of the three genes, and the expression value of the promoter in Corynebacterium glutamicum has been determined. In addition, we also designed to change the metabolic pathway in Corynebacterium glutamicum so that it can flow more carbon to glutamate synthesis. For biosafety, we consider designing a kill switch to make Corynebacterium glutamicum actively die after they complete the task. In order to verify the feasibility of the whole system, that is, whether it can produce sufficient polyglutamic acid and verify the effectiveness of kill switch, we simulated the whole system through modeling. At the same time, according to this model, we can be given a reasonable direction for the further transformation of our bacteria.
3. Problem restatement
Verify the feasibility of the modified Corynebacterium glutamicum secretion polyglutamic acid system, and look for the parameters that have a great influence on the secretion.
Verify the effectiveness of kill switch and point out the most influential parameters.
4. Assumption & symbols
- In the environment, except salinity (relative value) and alkalinity (relative value), the levels of other factors (temperature, humidity, light, etc.) remain unchanged. In order to make our model more responsive to the factors we want to focus on, it is essential to ignore the impact of other secondary factors.
- Corynebacterium glutamicum does not divide and is not affected by the survival threat and other effects that may be caused by other bacteria. In order to better reflect the performance of our design in chassis bacteria and simplify our model, we need to ignore the part of energy that bacteria need to resist external environmental stress and meet their own growth needs.
- Corynebacterium glutamicum does not synthesize flagella. The values of external salinity and alkalinity do not change. The purpose of this model is to verify the feasibility of design and make us have a deeper understanding of the system we designed. Therefore, it is not necessary to consider the changes of external factors caused by bacteria movement.
5. Model Establishment
5.1 Overview
Our model includes two independent parts, including secretion model and kill switch model. Different models perform different functions.
The secretion model corresponds to the γ-PGA synthesis and high yield of glutamate in the experiment. Through this model, we can verify the usability of the system we designed, and let us learn and build a more perfect system.
The other model is the kill switch model, which corresponds to the construction and verification of kill switch in the experiment. The model mainly describes the expression of ndoA protein under different salinity (relative value) and alkalinity (relative value). According to our design, the repressor protein expressed by Patp2 and Pgsib promoters will affect the promoter of ndoA, and these two promoters are regulated by salinity (relative value) and alkalinity (relative value) respectively. This complex regulation line needs a model to verify its effectiveness.
6. Secretory model
Corynebacterium glutamicum absorbs nutrients from the outside by a certain mechanism, whether in the fermentation environment or in the natural state. In the design of this part, we highlight the reaction process of PEP to oxaloacetic acid and α-ketoglutarate to glutamate. Closely coupled with these metabolic pathways are the partial reactions of TCA cycle and glycolysis pathway. For glycolysis reaction, we ignore the reaction before PEP and replace it with a comprehensive hypothetical enzyme (absorption enzyme), which can absorb external nutrients and convert into PEP. Accordingly, PEP and oxaloacetic acid need to go through the TCA before they can reach through a series of reaction pathways α- ketoglutarate. In the model, we simplified this point and also designed a hypothetical enzyme, which can catalyze citric acid into α-ketoglutarate in one step. The same is true for α-ketoglutarate to oxaloacetic acid. The model is very important for us to understand the operation of the whole system and the transformation relationship between substances.
The reaction pathway diagram we designed is as follows:
Fig.1 The metabolic pathway map based on the secretion model includes a variety of hypothetical enzymes
After simplifying the pathway of TCA and glycolysis, six related substances were retained, including PEP, oxaloacetic acid (OAA), citric acid, α-Ketoglutarate (KGA), glutamic acid (Glu) and γ-PGA. We establish the ordinary differential equation model of the mutual conversion relationship between them.
Our hypothetical absorption enzyme can absorb external influences and catalyze the formation of PEP. The nutrient kinetic absorption equation of osmotic nutrition organisms has been described by Michaelis equation1. In addition, PEP generates citric acid through our hypothetical carboxylase. Both are described by Michaelis equation. Also, PEP generates OAA by the enzyme encoded by ppc gene, and the kinetic equation of the enzyme is determined by reversible Michaelis equation. $$\frac{d[PEP]}{dt} = \frac{v_{m,a}[Nutrient]}{K_{M,a}+[Nutrient]}-\frac{v_{m,PEP2Citrate}[PEP]}{K_{M,PEP2Citrate}+[PEP]}- \frac{\frac{v_{m,f,ppc}[PEP]}{K_{s,ppc}}-\frac{v_{m,r,ppc}[OAA]}{K_{P,ppc}}}{1+\frac{[PEP]}{K_{s,ppc}}+\frac{[OAA]}{K_{p,ppc}}}$$
The hypothetical enzyme does not actually exist, so it is difficult to determine its maximum reaction rate and Michaelis constant. In determining the parameters of this hypothetical enzyme, we rely on two methods. One is to determine the reaction rate and reaction constant range of some hypothetical enzymes through experimental data. The second is to determine the influence of these parameters on the experimental results through the sensitivity test of parameters.
Oxaloacetic acid is also involved in three reaction pathways. KGA is catalysed to OAA by the hypothetical enzyme of downstream TCA. This hypothetical enzyme is actually the synthesis of all enzymes from KGA to OAA. OAA generates citric acid through citrate synthase, which is a definite enzyme. In addition, OAA participates in the reversible reaction with PEP, which is catalyzed by the enzyme encoded by ppc gene. $$\frac{d[OAA]}{dt} = \frac{v_{m,d}[KGA]^n}{K_{P,d}+[KGA]^n}-\frac{v_{m,CS}[OAA]}{K_{M,CS}+[OAA]}+ \frac{\frac{v_{m,f,ppc}[PEP]}{K_{s,ppc}}-\frac{v_{m,r,ppc}[OAA]}{K_{P,ppc}}}{1+\frac{[PEP]}{K_{s,ppc}}+\frac{[OAA]}{K_{p,ppc}}}$$
Citric acid is produced by PEP and OAA respectively, which is catalyzed to KGA by a hypothetical enzyme of downstream TCA. $$\frac{d[Citrate]}{dt} = -v_{m,u}*\frac{[Citrate]^n}{K_{P,u}+[Citrate]^n}+\frac{v_{m,PEP2Citrate}[PEP]}{K_{M,PEP2Citrate}+[PEP]} +\frac{v_{m,CS}[OAA]}{K_{m,CS}+[OAA]}$$
KGA is catalyzed by citric acid through hypothetical enzyme, and then transformed into OAA through another hypothetical enzyme. At the same time, the enzyme encoded by gdh gene can reversibly produce glutamate, which conforms to the reversible Michaelis equation. $$\frac{d[KGA]}{dt} = v_{m,u}\frac{[Citrate]^n}{K_{P,u}+[Citrate]^n}-\frac{v_{m,d}[KGA]^n}{K_{P,d}+[KGA]^n}-\frac{\frac{v_{m,f,gdh}[KGA]}{K_{s,gdh}}\frac{v_{m,r,gdh}[Glu]}{K_{p,gdh}}}{\frac{1+[KGA]}{K_{s,gdh}}+\frac{[Glu]}{K_{P,gdh}}}$$
Glutamate is reversibly generated from KGA. At the same time, polyglutamic acid is generated under the action of the enzyme encoded by capABC gene cluster and secreted to the outside of the cell. $$\frac{d[Glu]}{dt} = -\frac{v_{m,capABC}[Glu]}{K_{M,capABC}+[Glu]}+\frac{\frac{v_{m,f,gdh} [KGA]}{K_{s,gdh}}\frac{v_{m,r,gdh}[Glu]}{K_{p,gdh}}}{\frac{1+[KGA]}{K_{s,gdh}}+\frac{[Glu]}{K_{P,gdh}}}$$
Polyglutamic acid is catalyzed by an enzyme encoded by the capABC gene cluster. Since γ-PGA is a high polymer, the concentration value here does not represent the actual PGA concentration value, but represents the concentration value after PGA is degraded into monomer. In addition, PGA is decomposed in vitro, which conforms to the law of mass action. $$\frac{d[\gamma-PGA]}{dt} = \frac{v_{m,capABC}[Glu]}{K_{M,capABC}+[Glu]}-k_{f,deg}[\gamma-PGA]$$
In general, these equations are simplified. We imagined some enzymes to omit the intermediate complex steps, because if they are not omitted, many parameters will be involved, which greatly increases the complexity of the model and the difficulty of sensitivity analysis. These equations describe the mechanism of PGA synthesis by genetically modified Corynebacterium glutamicum. These parameters jointly determine the synthetic efficiency of engineering bacteria.
5.3 Kill switch model
In order to ensure biosafety, we also set up kill switch in the cell. This kill switch can respond to the changes of external salt ion concentration and pH value, so as to regulate the transcription of toxin protein ndoA gene.
Our model describes the transcriptional regulation of the repressor protein, the transcriptional inhibition of ndoA by repressor protein, and the translation of each protein. At the same time, we considered the self degradation of each protein and RNAs. Relative diagram is shown below:
Fig.2 Relationship between components in kill switch model
In order to express the regulatory effect of the protein product of gene X as a transcription factor on the transcription rate of gene Y, we can use input function: $$Y=f(X)$$
Here, Y represents the amount of transcript of gene y, and X represents the concentration of X protein (in the active state). Therefore, the input function represents the transcription rate of gene y regulated by active transcription factors. Generally, the function is represented by Hill function, which can effectively describe the real input function form. Transcription factors can activate or inhibit the transcription process of target genes, so there are also two forms of Hill function.2
When transcription factor X activates the expression of Y: $$Y=f(X)=a_0+\frac{a_{tr}X^n}{K^n+X^n}$$
When transcription factor X inhibits the expression of Y: $$Y=f(X)=a_0+\frac{a_{tr}}{K^n+X^n}$$
Pgsib and Patp2 regulate the expression of repressor protein mRNA downstream of these two promoters in response to changes in the level of external factors. At the same time, mRNA has its own degradation effect, which is expressed in the form of mass action law. $$\frac{[Pgsib-mRNA]}{dt}=a_0+a_{tr}\frac{[S]^n}{K_M^n+[S]^n}-k_{f,deg}[Patp2-mRNA]$$ $$\frac{[Patp2-mRNA]}{dt}=a_0+a_{tr}\frac{[A]^n}{K_M^n+[A]^n}-k_{f,deg}[Patp2-mRNA]$$
The equation describes the translation of repressor proteins. Downstream of pgsib and patp2 promoter are CDS of repressor protein. Both expressed mRNAs are translated into repressor proteins. The role of translation is expressed in the form of mass action law because there is no post transcriptional regulatory element on the mRNA. The repressor protein itself has degradation effect. $$\frac{[Repressor]}{dt}=k_f[Pgsib-mRNA]+k_f[Patp2-mRNA]-k_{f,deg}[Repressor]$$
The repressor protein inhibited the transcriptional initiation of Ptac. Therefore, it is reflected in the suppression form of Hill equation. $$\frac{d[Ptac-mRNA]}{dt}=a_0+a_{tr}\frac{K_M^n}{K_M^n+r_e[Repressor]^n}-k_{f,deg}[Ptac-mRNA]$$
The CDS region of ndoA protein is connected downstream of Ptac. This equation describes the translation of the transcriptional sequence initiated by the Ptac promoter. $$\frac{d[NdoA]}{dt}=k_f[Ptac-mRNA]-k_{f,deg}[NdoA]$$
In general, the above equation describes the mechanism of kill switch, including the interaction between each regulatory element, protein, mRNA and environmental factor level. In order to simplify the processing, we assume that the background expression rate and the maximum expression rate are constant in all transcription processes, and are the same in each equation. In this way, we can focus more on the interaction between these components.
6. Simulation & feasibility analysis
6.1 Secretory model
In order to test the feasibility of the whole system and the feasibility and effect of kill switch, we simulated the above dynamic equation network and the relationship between the product and the level of influencing factors. The action relationship of each component is described by kinetic equation, and the rationality is analyzed.
The above ordinary differential equations (ODEs) are solved by ode15s solvertype in MATLAB. The setting of initial values and parameters is based on the literature and experimental results. For the parameter values that are difficult to determine, the sensitivity analysis will be given later.
Fig.3 The fitting results of secretion model showed the change trend of each component in 300s.
The prediction results of the above model reveal the dynamic behavior of the system. When the external nutritional conditions are given, appropriate and unchanged, the genetically modified Corynebacterium glutamicum can finally produce γ-PGA, and will not cause the lack of TCA cycle intermediates. The simulation results of the model show the feasibility of the system, so that we can make further progress on the basis of our design.
6.2 Kill switch model
The efficiency of kill switch determines the security of our project. Our kill switch is affected by salinity and alkalinity in the external environment (the two factors are expressed in relative values). The toxin protein was not expressed when one factor was at a high level, but a large number of toxin proteins were expressed when both factors were at a low level. To verify our results, we verified our kill switch efficiency under four conditions.
Fig.4 Changes of expression of various substances in 10 s under different initial conditions of salinity and alkalinity
We set up a combination of two levels of two factors separately, a total of four sets of parameter settings. When both A and S are set to 1, the expression of repressor is very high and the amount of toxin protein is low. When either of them was 0.4 and the other was 1, the expression of repressor decreased slightly and the expression of toxin protein increased slightly. When both values were 0.4, the expression of repressor decreased significantly and the expression of toxin protein increased significantly. The result of the simulation can prove that our Kill Switch can work effectively to ensure biosafety.
7. Sensitivity Analysis
7.1 Secretory model
In the secret model, some parameters are difficult to determine and may have a great impact on the model fitting results. We take all the relative parameters for analysis. In addition, our design also involves the control of ppc and odhA (TCA downstream hypothetical enzyme) expression. In order to observe the effect of controlling the expression of key metabolic enzymes, we also modified the model.
Fig.5 Variation of components under initial conditions of different nutrient values (from 0.8 to 1.2). Y axis represents Relative Value(au)
Scan the value of nutrition concentration from 0.8 to 1.2. It is found that when the amount of nutrient increases, the maximum amount of γ-PGA is also gradually increasing.
Because the model contains a large number of parameters, we analyzed the γ-PGA sensitivity of each parameter. It is found that the most influential factor is about the expression of the protein encoded by capABC gene cluster. In addition, the decomposition rate of γ-PGA and the expression of the protein encoded by gdh gene will also affect γ- PGA has a great impact on the yield. The influence degree of other factors is relatively low.
Fig.6 Effect of changes in various parameters on γ-PGA yield (sensitivity)
In order to describe the transformation of ppc gene and gdh gene expression value, we multiplied the reaction rate of PEP to OAA, KGA to Glu and KGA to OAA by a variable to simulate the expression. $$PEP↔OAA,v=a_1\frac{\frac{v_{m,f,ppc}[PEP]}{K_{s,ppc}}-\frac{v_{m,r,ppc}[OAA]}{K_{P,ppc}}}{1+\frac{[PEP]}{K_{s,ppc}}+\frac{[OAA]}{K_{p,ppc}}}$$ $$KGA↔Glu,v=a_2\frac{\frac{v_{m,f,gdh}[KGA]}{K_{s,gdh}}\frac{v_{m,r,gdh}[Glu]}{K_{p,gdh}}}{\frac{1+[KGA]}{K_{s,gdh}}+\frac{[Glu]}{K_{P,gdh}}}$$ $$KGA→OAA,v=a_3\frac{ v_{m,d} [KGA]^n}{K_{P,d}+[KGA]^n }$$
Fig.7 Effects of changes in expression parameters of three enzymes on PGA yield (sensitivity)
As can be seen from the figure, \(a_2\), \(a_3\) has a great impact on the yield of γ-PGA, of which \(a_2\) has the greatest impact. And the sensitivity contribution of \(a_1\) to PGA was almost 0. This may be because this model only considers the reactions of glycolysis and TCA cycle, so the influence of the enzyme amount of PPC enzyme is not shown. This result also plays a guiding role in our experiment. We will give priority to adjusting gdh gene expression and α- Expression of ketoglutarate dehydrogenase related genes.
7.2 Kill Switch model
There are also many parameters to be determined in the kill switch model. Therefore, we analyzed the sensitivity of each parameter. The sensitivity of each parameter to the amount of toxin protein and repressor was observed.
Fig.8 Effects of different parameter values on the expression of ndoA and repressor protein (sensitivity)
It can be seen from the figure that the most decisive parameters for the change of the amount of toxin protein are the rate constant during the translation of toxin protein mRNA and the maximum rate of transcription of toxin protein mRNA itself. The most decisive parameters for the change of the amount of repressor are the expression of repressor mRNA under the regulation of salt promoter and the expression of repressor mRNA under alkali promoter. The sensitivity of other parameters is low. It shows that the overall expression of toxin protein is most affected by the change of repressor, which depends on the expression promotion of salt promoter and alkali promoter. It shows that the construction of the whole system is successful. We can focus on the influence of salt promoter and alkali promoter itself on the changes of external salinity and alkalinity.
8.Conclusions
8.1 Secretory Model
With a certain initial value and reasonable parameter setting, our secretion model shows that it has a certain amount of γ-PGA production, and other TCA intermediate metabolites will not be reduced to 0, so it will not affect the normal metabolic pathway of bacteria too much. The feasibility of the system is proved. The changes of various parameters have more or less affected the yield of γ-PGA. The parameters related to capABC enzyme and γ-PGA degradation rate had the highest influence on the results. The related parameters of the enzyme downstream of ketoglutarate and the parameters of the enzyme encoded by gdh gene are also very high. We can use these enzymes as a starting point to improve the production of γ-PGA as much as possible.
8.2 Kill Switch Model
When setting the initial values of different salinity and alkalinity, our model results show that the kill switch system has different response to the changes of external conditions. When the salinity and alkalinity were at least one high, the expression of Repressor was maintained at a high level, and the expression of ndoA was maintained at a low level. When the salinity and alkalinity were low, the expression of Repressor decreased rapidly and the expression of ndoA increased rapidly. This result is in line with our design. In the sensitivity test of the parameters, we found that the inhibition ability of the Reperssor and the response ability of the salt and alkali inducible promoter to the changes of environmental factors are the key to determine the efficiency of kill switch.
9.How the model results affected the project design and development
By modeling the system in which Corynebacterium glutamicum produces and secretes γ-PGA, and Corynebacterium glutamicum responds to changes in the external environment and starts the kill switch mechanism, we can have a deeper understanding of the behavior, reliability and feasibility of the system we designed, and get further guidance to our experiment. Therefore, we have done substantial and excellent work to achieve this goal.
The results of secretory system model showed the secretion of γ-PGA. Due to the secretion of γ-PGA, the content of metabolic intermediates of bacteria itself decreased. Therefore, the viability of this kind of bacteria must be lower than that of wild bacteria. This makes us consider that in the actual design of bacterial agents, we need to provide some nutrients that can be absorbed by engineering bacteria, and take into account the cost. Through the sensitivity analysis of parameters, we found that many parameters have an impact on the secretion of γ-PGA. The parameters related to the enzymes encoded by capABC gene cluster, the enzymes related to glutamate production, and the parameters related to the decomposition of γ-PGA in vitro all play a considerable role in the synthesis of polyglutamic acid. Therefore, we need to seriously consider the expression of capA, capB and capC genes in Corynebacterium glutamicum.
Considering the results of sensitivity analysis, in addition to the parameters mentioned above, we also note that the relevant parameters of hypothetical TCA downstream enzyme also has a significant impact on the results. Therefore, we consider modifying the expression of the enzyme encoded by the odhA gene.
The equations of kill switch also play a role in our design of gene circuit. We find that the whole circuit is feasible under appropriate parameters. However, more attention should be paid to the sensitivity of the two inducible promoters to the external environment and the inhibitory effect of Repressor on downstream promoters. In the sensitivity analysis, the corresponding parameters are also outstanding. Therefore, before constructing switch, we need to fully determine the effects of the two promoters and Repressor.
10.References
[1] Ding.,et al.(2007)The 7th general meeting and 14th Symposium of algal branch of China Society of marine Limnology.
[2] Yang, H.-T., et al.(2007)An Analytical Rate Expression for the Kinetics of Gene Transcription Mediated by Dimeric Transcription Factors. The Journal of Biochemistry 142, 135-144.
[3] Zhao. (2016)Research of the Synthetic Gene Oscillator Based on Simulation. Zhejiang University.
Cellular Automation Model
1.Overview
Saline alkaline soil has increasingly become an important and harmful agricultural and land problem. The main purpose of the project is to use Corynebacterium glutamicum as the chassis bacteria and construct engineering bacteria for the restoration of saline alkali soil. For biosafety, the bacteria we designed will express suicide protein after treating saline alkali soil. In order to guide the design of our project and the implementation of bacterial agents, and enable us to further understand the dynamic change process and mechanism of bacterial agent treatment of soil, we use cellular automata model to simulate the performance of bacterial agents after release, predict the results of different release strategies, and evaluate the impact of different factors (suicide efficiency, treatment efficiency, miscellaneous bacteria growth efficiency, etc.) on the work efficiency of the system. In addition, this model can also be extended to simulate the release and leakage of various microorganisms.
2.Background
2.1 Background of the Project
Saline alkali soil contains a large amount of soluble saline alkali, which endangers environmental safety and agricultural production. Saline alkali soil has become a global problem, causing great agricultural and environmental hazards (Wu, Li et al. 2013). In the context of increasing population, mankind urgently needs more arable land to maintain food supply. Therefore, our team is committed to the treatment and restoration of saline alkali soil. It is expected that this goal can be achieved by engineering bacterial repair. We take Corynebacterium glutamicum as the chassis bacterium. Relying on its strong glutamic acid production capacity, we can produce it through synthetic biology γ-Polyglutamic acid has the ability to commit suicide according to the change of environmental salinity. As a result, bacterial agent that is able to treat the saline alkaline soil is our project application form.
2.2 Background of Corynabacterium glutamicum
Corynebacterium glutamicum has become one of the most important bacteria in bioengineering today, and has produced more than 2 million tons of amino acids every year, of which more than 1 million is sodium glutamate, which is an important flavor enhancer in Food Engineering (Eggeling and Bott 2005). Corynebacterium glutamicum was once classified as Micrococcus because of its morphological characteristics (Udaka 1960). It is a Gram-positive bacterium with a culture cycle of 1 ~ 2 days at 28 ℃ (Baumgart, Unthan et al. 2013). Corynebacterium glutamicum is not pathogenic. According to the German TRBA standard, the risk group of the strain is 1, which has sufficient biosafety. The culture and transformation methods of the strain and the expression vectors related to the strain have been reported(Jakoby, Ngouoto-Nkili et al. 1999, Jang and Britz 2000).
3.Problem statement
Establish the model to observe the performance and distribution of engineering bacteria with the ability of secreting polyglutamic acid and suicidal ability in saline alkali soil.
The effect of suicidal capacity on the distribution of engineering bacteria.
Possible impact of colonization of other flora on the distribution of engineering bacteria after soil remediation.
The possible effects of bacterial agents with different repair efficiency on the repair results.
4.Assumptions & Symbols
- There is few bacteria in newly restored soil. Through our literature research, the biodiversity of saline alkali land itself is worse than that of normal soil. Therefore, we assume that when bacteria just restore saline alkali soil, the restoration of biodiversity lags behind.
- The genes of engineering bacteria will not lost, and natural bacteria do not mutate or obtain edited genes. In the actual natural environment, horizontal gene transfer and mutation are common. However, this phenomenon adds many uncontrollable factors, so we choose to ignore this phenomenon.
- Saline alkali soil environment is not affected by other factors, such as weather, temperature, etc. And it is the most suitable environment for engineering bacteria. Saline alkali land itself is affected by many other factors, but it is not the focus of our model.
- Ignore the temperature difference due to diurnal variation. This is an environmental factor. In order to simplify the model, we need to assume that the temperature of the environment is constant, and this assumption can also save the trouble of considering the influence of temperature on various parameters.
- It is assumed that the whole saline alkali soil system is a two-dimensional system, and all indexes are evenly distributed except water potential (representing salinity) and pH (representing alkalinity). In this system, we only consider the possible effects of salinity and alkalinity and the improvement of saline alkali land by bacteria. In order to simplify the model and avoid other negative effects, we adopt a two-dimensional model for simulation.
- The impact of the surrounding environment of saline alkali land is not considered. Saline alkali land is locally distributed, and saline alkali land is surrounded by various other soil types. The model only considers a certain saline alkali land and there are saline alkali land all around. Because saline alkali land may have boundaries with various types of soil, considering the impact of all types of soil on saline alkali land is undoubtedly a complex and rarely useful work.
The specific value of the parameter can be seen in the source code (at the bottom of the web page)
5.Model Establishment
5.1 What is Cellular Automation?
Cellular automata is a discrete dynamic model. It is discrete in time, space and state. The model is not determined by strictly defined physical equations or functions, but consists of a series of rules constructed by the model. This concept was first proposed by Stanislaw Ulam and John von Neumann in 1940s to simulate the self replication function of life systems. This kind of model can be used to study many phenomena and is widely used in sociology, ecology and computer science. The basic unit in the model is called cell. A cell has a variety of discrete states, and the change of its state with time depends on the cell and nearby cell states. Therefore, cellular automata has significant advantages in simulating the change of dynamic system, which is helpful for people to understand and deeply understand the change mechanism of dynamic system.
5.2 Cellular states & Rules
- Environmental Cellular In this project, we consider the process of saline alkali soil improvement. Due to its low water potential value and high pH value, many ordinary microorganisms and plants are difficult to survive in this environment, and the biodiversity is low (Zhang, Liu et al. 2021). Shortly after the soil is repaired, even if its water potential and pH value return to normal, the biodiversity has not been restored. At this time, the soil is called nascent soil. This kind of soil will lead to the suicide of engineering bacteria, but it does not affect the repair process of engineering bacteria on saline alkali soil. After a period of time, due to the colonization of living bacteria or spores floating in the air in the nascent soil, biodiversity will recover soon. This kind of soil for restoring biodiversity is called improved soil. Because there may be bacteria in the improved soil that have a certain impact on the survival of engineering bacteria, the improved soil will not only lead to the suicide of engineering bacteria, but also affect the efficiency of engineering bacteria in repairing saline alkali soil.
- Bacterial Cellular In a micro environment, we consider the distribution of bacteria in saline alkali soil. Due to the action of engineering bacteria, the saline alkali soil around engineering bacteria will be slowly restored to nascent soil. The engineering bacteria existing in the nascent soil died quickly. In addition, the colonies colonized in the nascent soil will have a certain impact on the survival of the engineering bacteria. We consider setting the cell of engineering bacteria and set certain rules for the cell.
- Rules
γ- Polyglutamic acid can well inhibit the occurrence of secondary salinization because of its water retention. However, an important factor of soil salinization is climatic conditions. even if γ- Polyglutamic acid can maintain the water in the soil to a certain extent, but can not change the climatic conditions. Therefore, both the nascent soil and the improved soil will become saline alkali soil again to a certain extent.
Based on the above discussion, this model sets up three soil States, which have a certain transformation relationship. We set the size of cellular automata to 500 × 500 cells, each cell represents a soil unit, regardless of the internal differences of each unit. In the conversation with Mr. Li Huixin (see Human_Practice), we learned that soil colloids have changes in indicators on the micro domain scale. And considering that when the bacterial agent is actually applied, the bacterial agent is covered on the soil and permeated, and covered on the soil surface on the macro scale, rather than concentrated application of the bacterial agent at a certain point. Therefore, when this 500 × 500 system of CA corresponds to the actual space, it only represents a very small area, such as 1mm×1mm.
Fig.9 The 500 × 500 cells system of this model represents only a small part of the soil under the actual size.
Conclusion:
Three cellular states: saline alkali soil, nascent soil and improved soil
The three soils have a certain initial proportion.
The size of cellular automata is 500 × 500 cells, corresponding to only a small area in the actual space.
Conclusion: There are two cellular states: engineering bacteria and natural bacteria (equivalent to improved soil)
There is a certain relationship between the two.
Based on the above cell and project background, we set the conversion rules between cell states and give a reasonable explanation.
Rule 1. SAS has no less than one EB cell in its molar neighborhood and is transformed into EB under a certain probability
Corynebacterium glutamicum has good salt and alkali tolerance. It still has a good growth rate when pH = 9.0 and sodium ion concentration is 0.6M (Xu, Zheng et al. 2018). However, its growth efficiency is lower than that in normal medium. In the model, it shows diffusion under a certain probability, which is lower. Since the type, size, surface charge and nutritional status of microorganisms will affect the diffusion of bacteria in soil (Han 2016), this probability is difficult to determine accurately. Therefore, based on the literature data, we set \(p_G= 0.003\) as an example, sensitivity analysis will be carried out later.
Rule 2. SAS has no less than one EB cell in its molar neighborhood and is transformed into NS under a certain probability
For saline alkali soil, polyglutamic acid released by engineering bacteria around it will improve saline alkali soil into nascent soil. And the improvement efficiency of polyglutamic acid is limited, so it is reflected in a certain transformation probability in this model. The transformation probability here reflects the remediation ability of the engineering bacteria to the soil. If the repair ability is too weak, the repaired small area will soon be occupied by other colonies, which will quickly endanger the survival of bacteria. According to the experimental data, the ability of bacteria to secrete polyglutamic acid in ODE model and the performance of polyglutamic acid in soil remediation, we set the probability as 0.004.
Rule 3. The molar neighborhood cells of NS are all NSs and are transformed into SAS under a certain probability
This rule indicates the secondary salinization of nascent soil. Due to the existence of polyglutamic acid and the effect of climate takes a certain time, it is reflected in a certain transformation probability in the model, and the probability is very low. In order to estimate the probability of land secondary salinization, we obtained the area data of saline alkali land in Songnen Plain of Northeast China and Western Jilin (Li et al. 1998). As shown in the figure below. The saline alkali land diffusion area of the two places in the past 40 years is equivalent to 41% and 54% of the area in the 1950s. This is the change of saline alkali land area under the condition of human activities. With the water retention of polyglutamic acid, its diffusion rate will be further reduced. We assume that polyglutamic acid will reduce the ratio of diffusion area to initial saline alkali land area in 40 years to 40%. We use ergodic algorithm to search the probability value that can best meet this hypothesis. Here, we set \(p= 0.0000018\)
Fig.10 Changes of saline alkali soil area with time in Songnen Plain and Western Jilin.
Rule 4. EB becomes SAS under certain probability
This rule describes the natural death of bacteria in saline alkali soil environment.
Rule 5. NS includes EB in its molar neighborhood and is converted to EB under a certain probability
This rule describes the growth of engineering bacteria around the soil and the diffusion of bacteria at the interface between saline alkali soil and nascent soil to nascent soil.
Rule 6. EB is converted to NS under certain probability when its molar neighborhood does not include SAS
This rule describes that the bacteria start kill switch and commit suicide when the water potential and pH of the surrounding environment reach normal. This probability reflects the efficiency of kill switch. According to our experimental data and the results of ODE model, we set this probability to 0.8 to reflect the role of the actual kill switch.
Rule 7. NS includes SAS in its molar neighborhood and is converted to SAS under certain probability
This rule describes the secondary salinization of fresh soil. Unlike rule 3, this rule emphasizes the diffusion of SAS. Due to the saline alkali stress, the biodiversity of organisms in the saline alkali soil newborn soil interface will be reduced, so as to accelerate the process of salinization. The probability of this rule is higher than that of rule 3.
Rule 8. The molar neighborhood cells of NS are all NSs and transformed into IS under a certain probability
This rule describes the colonization of airborne spores or living bacteria in nascent soil and the restoration of biodiversity. Considering the vast area of saline alkali land, even if the water potential and pH value are restored rapidly in the core area of saline alkali land (i.e. far away from the normal soil area), its biodiversity will not return to the level of normal soil quickly.
Rule 9. The molar neighborhood cell of NS does not contain SAS and contains IS, which is transformed into IS under a certain probability
This rule is different from rule 8. This rule reflects the diffusion of bacteria in IS in NS. Natural bacteria grow faster and do not face salt and alkali stress. Here, we set \(p = 0.1\) as an example.
Rule 10. EB contains IS in its molar neighborhood and is transformed into IS under a certain probability (the probability increases with the increase of the number of is cells in the molar neighborhood)
This rule describes the existence of harmful bacteria unfavorable to the growth of engineering bacteria in soil after Biodiversity Restoration. In addition, harmful bacteria and engineering bacteria rob natural resources, resulting in the death of engineering bacteria. The probability is closely related to the competition between natural bacteria and engineering bacteria, and the probability is set to 0.01 for the time being. We will conduct sensitivity analysis on it later and look for a threshold. When the probability is less than the threshold, the engineering bacteria can complete the task well; When it is above the threshold, the engineering bacteria will be greatly disturbed by natural bacteria and can not work normally.
Rule 11. IS is transformed into SAS under a certain probability (which increases with the increase of the number of SAS cells in the molar neighborhood)
This rule also describes the secondary salinization. Because the climate has not been changed, there is still a certain probability that the improved soil will be transformed into saline alkali soil, and the higher the salt alkali stress around it, the faster the salt alkalization rate. This probability must be less than the probability of rule 3. Because NS in Rule 3 has poor and fragile soil biodiversity, is in this rule has restored biodiversity to a certain extent. Here we set \(p = 0.0000009\).
The above 11 rules illustrate the model's description of objective reality. Summarized as the following figure. Since the occurrence conditions of these events are not decisive, they all need to meet a certain probability to trigger. By changing these probability values, we can know the impact of the probability change of each rule on the result when the probability of other rules is certain.
Fig.11 Schematic diagram of intercellular transformation relationship
5.3 Simulation & Analysis
When adding bacterial agent, the end user can control the amount of bacterial agent. Therefore, we set the delivery probability \(p_d\) to reflect the actual investment.
When the model is initialized, $$p_d=\frac{n_0 (EB)}{n_0(SAS)}$$
Take liquid fungicide as an example. When the fungicide is actually applied, every part of the target soil will be sprayed. But on the scale of 1mm × 1mm, not every place is covered with liquid bacterial agent. In the small part not covered by bacterial agent, it needs to be realized by the growth and diffusion of adjacent engineering bacteria. On this scale, we simulate the distribution and quantity relationship of various cells with the increase of iteration times under different delivery probabilities. The simulation results of each delivery probability are displayed in charts and videos.The information in the video can tell us the expression and migration forms of four kinds of cells in soil. We can start with the information in the video to further understand the working mechanism of our system. The meanings of four cell colors in the video are shown in following picture. It should be noted that in the early stage of video playback, there will be color dislocation because IS cells do not exist yet. The spread of bacteria in the soil can be seen in the video. The faster the bacteria spread, the smaller the proportion of white cells (saline alkali land), the better the effect.
\(p_d=0.0002\)：
\(p_d=0.0004\)：
\(p_d=0.0006\)：
\(p_d=0.0008\)：
\(p_d=0.0010\)：
\(p_d=0.0012\)：
\(p_d=0.0016\)：
\(p_d=0.0020\)：
\(p_d=0.0025\)：
Through the above simulation results, the following common phenomena and conclusions can be obtained:
- Due to the treatment ability of engineering bacteria and the role of kill switch, the living range of engineering bacteria is concentrated at the boundary between newborn soil and saline alkali soil.
- The flora harmful to engineering bacteria in the improved soil has an impact on the growth and diffusion of engineering bacteria, but the impact is not as serious as expected.
- Nascent soil is intermediate. After being improved into nascent soil, saline alkali soil is quickly replaced by improved soil.
- The larger the amount of bacteria, the better and faster the improvement and treatment effect.
- When the amount of bacteria is small (\(p_d<0.0010\)), finally the ratio of improved soil to saline alkali soil increases with the increase of bacterial dosage, and the time to reach this ratio decreases rapidly with the increase of bacterial dosage.
- When the amount of bacteria is large (\(p_d≥0.0010\)), the ratio of improved soil to saline alkali soil further decreases with the increase of bacterial dosage, but the time to reach this ratio still decreases with the increase of bacterial dosage.
- When the amount of bacteria is small, the cells of nascent soil and engineering bacteria exist for a long time, but the peak value is low. The larger the amount of bacteria, the shorter the survival time of nascent soil and engineering bacteria, but the higher the peak value.
When the amount of bacteria is small, the increase of the amount of bacteria can significantly improve the improvement efficiency and speed. However, after reaching a certain value (0.0010), the increase of bacterial dosage mainly plays a role in the improvement speed. This plays a guiding role in the actual delivery strategy, that is, when considering the actual cost and demand time at the same time, users can choose the most appropriate delivery rate for delivery.
When p was 0.0002, 0.0007, 0.0012, 0.0017 and 0.0022, after all the engineering bacteria died, we collected the cell number ratio of improved soil and saline alkali soil at this time.
Fig.12 Change curve of IS/SAS with bacterial feeding rate
It can be seen the final IS/SAS has an obvious logarithmic relationship with the bacterial rate. This result of the model further confirms that when using the bacterial agent, we should not blindly pursue more, but should seek a certain balance between the amount and the cost.
Of course, if the model iteration continues, the proportion of saline alkali soil will eventually be back to 100% due to rule 3, rule 7 and rule 11.
5.4 Development：Factors that have a significant impact on the effect
In addition to the bacterial feeding rate, there are two important factors that may have a great impact on the treatment efficiency including kill switch efficiency and effectiveness of bacterial agents. The two indexes are expressed by \(p_k\) and \(p_e\). Set the bacterial feeding rate \(p_e=0.001\), then observe the change of model result output by changing these two probability values.
5.4.1. Kill Switch Efficiency
The probability of kill switch in the model is \(p_k\). When other parameters are constant, we observe the change of the number of iterations required when the engineering bacteria is finally 0. Set \(p_k\) increases from 0.1 to 0.9, effective probability is set to 0.004, and the model is used for simulation. Due to the randomness of cellular automata, we set up two repetitions.
Fig.13 The number of iterations at the time of death of all engineering bacteria varied with the suicide rate. A total of 2 repetitions were performed
It can be seen that the effect of suicide probability on the number of iterations when all engineering bacteria die is not prominent. This may be because the different probability of suicide also affects the migration of bacteria. Although the results show that the suicide efficiency has no significant impact on the treatment results, the efficiency reflects the important issue of biosafety. Therefore, the higher the parameter, the better.
5.4.2. Effectiveness of bacterial agents
The effectiveness of bacterial agents (i.e. the sum of the efficiency of engineering bacteria secreting polyglutamic acid and the treatment efficiency of polyglutamic acid itself) can directly affect the final treatment efficiency. We set \(p_k\) to 0.8, and \(p_e\) increases from 0.001 to 0.01, Record the IS/SAS and the number of iterations.
Fig.14 Variation curve of IS/SAS and the number of iterations at the time of death of all engineering bacteria with bacterial agent effectiveness.
The treatment efficiency of engineering bacteria has a great impact on the results. With the increase of \(p_e\), IS/SAS has a decreasing trend. After the engineering bacteria have treated the surrounding land, they will be killed by kill switch and lose the opportunity to spread to the surrounding areas. Therefore, the best bacterial agent needs to find a certain balance between the limited diffusion rate of bacteria and the treatment efficiency.
6.Sensitive Analysis
It is difficult to determine the bacterial migration growth rate \(p_G\) and bacterial colonization rate \(p_H\) , and may have a great impact on the experimental results. In order to measure the impact of these two parameters on the results of our model and optimize the model, we analyze the results by changing these two parameters.
Fig.15 Relationship between iteration times and \(p_G\) and \(p_H\) when all engineering bacteria die
Fig.16 Relationship between \(\frac{IS}{SAS}\) and \(p_G, p_H\)
We set \(p_G∈[0.001,0.011],p_H∈[0.00001,0.00061]\), then count the number of iterations and \(\frac{IS}{SAS}\). It is observed that \(p_G\) has a significant effect on the number and ratio of final iterations. With the increase of \(p_G\), the number of iterations when all engineering bacteria die gradually decreases and the ratio gradually increases. It shows that the growth of engineering bacteria has an obvious effect on the treatment efficiency of bacterial agents. However, the impact of \(p_H\) is not as great as expected, indicating that it is not necessary to consider too much the impact of local flora on the treatment efficiency when actually adding bacterial agents. However, the survival of engineering bacteria in the target saline alkali land needs to be considered.
7.Conclusions
The final treatment effect changed significantly when different bacterial dosage was applied. And there is an obvious logarithmic relationship between them. It shows that after the amount of bacteria feeding rate is increased to a certain value, the treatment effect will not be significantly improved with the increase of the amount of rate. Therefore, when selecting the delivery strategy, we need to take into account the cost to find an optimal method. The number of iterations when all engineering bacteria die has no obvious relationship with the change of suicide efficiency, indicating that we do not need to over consider the possible side effects caused by too high kill switch efficiency.
It should be noted that the treatment effect does not increase monotonicly with the bacterial treatment capacity（ γ-PGA secretion capacity) increased. This may be because too good treatment results incur faster growth of harmful bacteria, then squeeze the diffusion space of engineering bacteria.
Under a certain treatment efficiency, the impact of harmful bacteria may not be as great as we think. Because the growth of harmful bacteria depends on the soil treated by engineering bacteria, the growth of harmful bacteria itself is subject to certain self inhibition. One parameter closely related to the treatment ability is the diffusion rate of bacteria, which can be achieved by increasing the amount of bacteria.
8.How the model results affected the project design and development
Through the cellular automata model, we understand the behavior of the engineering bacteria we designed in practical work. This model helps us understand the mechanism of the system.
The cellular automata model evaluates the behavior of our engineered bacterial system in a certain range, including how it spread and how it treat. Based on assumptions and our design, we designed 11 rules to describe and define the behavior of the system. Through the simulation of the model, we found that there was a logarithmic relationship between the amount of bacteria and the treatment effect (expressed in \(\frac{IS}{SAS}\) here), and the fitting effect was very good. The results remind us that we need to design according to the actual situation, cost budget and other factors. In addition, the number of model iterations when all engineering bacteria die has little to do with the efficiency of kill switch. It shows that we do not need to worry about the possible impact on the effect when the suicide efficiency is too high. We only need to fully consider improving the suicide efficiency of switch, so as to ensure the biosafety as much as possible.
9.References
[1] Baumgart, M., et al. (2013). "Construction of a prophage-free variant of Corynebacterium glutamicum ATCC 13032 for use as a platform strain for basic research and industrial biotechnology." 79(19): 6006-6015.
[2] Eggeling, L. and M. Bott (2005). Handbook of Corynebacterium glutamicum, CRC press.
[3] Jakoby, M., et al. (1999). "Construction and application of new Corynebacterium glutamicum vectors." 13(6): 437-441.
[4] Jang, K.-H. and M. L. J. B. L. Britz (2000). "Improved electrotransformation frequencies of Corynebacterium glutamicum using cell-surface mutants." 22(7): 539-545.
[5] Udaka, S. J. J. o. b. (1960). "Screening method for microorganisms accumulating metabolites and its use in the isolation of Micrococcus glutamicus." 79(5): 754-755.
[6] Wu, Y., et al. (2013). "Organic amendment application influence soil organism abundance in saline alkali soil." 54: 32-40.
[7] Xu, N., et al. (2018). "The lysine 299 residue endows the multisubunit Mrp1 antiporter with dominant roles in Na+ resistance and pH homeostasis in Corynebacterium glutamicum." 84(10): e00110-00118.
[8] Zhang, Z., et al. (2021). "Organic fertilizer enhances rice growth in severe saline–alkali soil by increasing soil bacterial diversity."
[9] Han. (2016).Transport Mechanisms and Mathematical Model of E.coli in Soil. Tianjin University of Technology.
[10] Li., et al. (1998). Study on land secondary salinization in Songnen Plain.".Geographical Science. (03): 268-272.
[11] DSMZ.https://www.dsmz.de/collection/catalogue/details/culture/DSM-102070.