Team:HUST-China/Model

iGEM_2019
model-abstract
Model
banana-tree
Abstract
With the benefit of Synthetic Biology, we built BanaMax which is an optimized degumming kit for banana fiber. In order to make sure that our kit could fulfill its duty, we built two mathematical models to simulate it.

One model is the ODEs model, this model could simulate the gene circuit of Pichia pastoris, evaluate the properties of pH control system and predict the amount of secreted enzymes.

What’s more , for the purification of the enzymes, we needed to put a His-Tag, either at the N- or C-termini of the protein sequence. So we realized a de novo modeling using the webserver I-TASSER to modelize the 3D structures of the enzymes with the His-Tag at both ends.

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

banana-tree
Gene pathway model

Preface
1. Why do we need a model

Developing a mathematical model of the systems we established (or plan to establish) is a vital step in our project. A model, namely, is a mathematical representation of the system. The model could help us in the following aspects:
- Confirming that our system has the potential to work properly.
- Determining optimal conditions for the system.
- Estimating biological parameters.
- Understanding the mechanism of the system better.
- Optimizing the design of our final project.

2. What can we get from modeling

Based on logistic function, Monod kinetics,etc, we can describe a cell or the group of cells mathematically. By quantifying cells, a number of parameters can be predicted, including the colony growth rate of Pichia pastoris , the expressions of proteins and the changes of the PH during the growth of cells. More significantly, the model of the project gives us a visualize and intuitive result, which could guide the experiment before running it.

3. General assumptions

In modeling the growth of the pastoris, in order to analyze the growth of cells efficiently and mathematically, we use unsegregated model, in which average cellular properties and balanced growth are assumed. Also, each cell is deemed identical, meaning that population is regarded as one-component solute and cells are produced of the same quality with the common components — at various stage of the co-culture, new cells are assumed to be equal to old cells.
The model is also unstructured, in the sense that the intermediate reactions are omitted.To simply the model, we assume that the environment remain unchanged during the growth of Pichia pastoris and the toxic effects of the metabolites from these reactions are ignored, shows the K is a non-time-value parameter. Likewise, the popluation of Pichia pastoris should not exceed the carrying capacity of the bioreactor.
And in modeling the gene path, we ignoring the fluctuation and the external noise in expression process, so that the expression speed of the protein can be regard as a Time-dependent function and later be described as a partial differential equation.

Results and Conclusions
Through the analysis of the mathematical model of the system, we use MATLAB to simulate, through parameter fitting and data query, and finally get the following graphics.


Figure1 pH and the Concentration of Pichia pastoris

Our project is based on pH control system. Under the analysis of the picture and mathematical model, we can get the state of the system under given assumptions. Assuming that the initial condition is that the concentration of Pichia pastoris is 10g/L and pH is 7, then under these conditions, we can see from the above figure that Pichia pastoris, after the retarded period, enters the logarithmic phase and begins to grow in large quantities due to suitable conditions. The concentration of bacteria increases. In this process, Pichia pastoris secretes PelA, and PelA degrades pectin to produce galacturonic acid, which leads to the decrease of pH in the solution. Through the inquiry data, we learned that Pichia pastoris had good growth at pH 4.5 and 6.5, so the initial decrease of pH did not slow down the growth rate of Pichia pastoris. However, when pH decreased to 3 or less, the growth of Pichia pastoris was seriously inhibited, resulting in an increase in the mortality of Pichia pastoris. In addition, considering that the cell concentration will reach the maximum environmental capacity of Pichia pastoris, the concentration of Pichia pastoris around t=4h began to decrease. After that, when pH was near 2, the expression of basic protein gene in Pichia pastoris is initiated in acidic environment, the expression of alkaline protein gene in Pichia pastoris increases, gradually stabilizes, and finally the environment is slightly acidic.Finally, the cell concentration is also maintained at a certain value. Through this model, on the one hand, we deepen our understanding of the system, on the other hand, it is also verified that our system is stable and feasible, so that the change of pH will not lead to the collapse of the system.


Figure2 the Concentration of Enzymes

From the picture, we can see the expression of three enzymes: PelA, SALC and VP secreted by Pichia pastoris. When the pH value was about 6, the first large amount of PelA is expressed, which lead to the further decrease of pH. When pH is about 5, SLAC and VP begin to express a lot. However, when pH is reduced to around 2, the basic protein gene is highly expressed, which leads to a rebound of pH, thus ensuring that our system will not collapse.


Figure3 the Concentration of Pectin,Lignin and Hemicellulose

In the third picture, we can see the degradation of the three main substances that need to be degraded. There is the first large amount of PelA expression, corresponding to it, pectin is also the first to be degraded. Subsequently, due to the high expression of SLAC and VP, hemicellulose and lignin began to degrade.



As a mathematical representation of the system, the model sets some parameters and confirms that our system has the potential to work properly. Therefore, We can provide the experimental group with the data and methods we obtained in the process of consulting data and analysis, as well as the parameters obtained in the process of simulation, so that the experimental group can obtain the best results in a shorter time and lower cost.
Through this model, we deepen our understanding of the system, from the most basic chemical reactions, enzymatic reactions, to dynamic changes in substance concentration and the growth of Pichia pastoris colonies.
From the first picture, we can see that the initial growth rate of Pichia pastoris is very slow. In order to speed up the expression of enzyme gene and reduce the degradation time of non-cellulosic substances in banana, we suggested to the experimental group that measures should be taken to shorten the retard time of Pichia pastoris in the new environment. For example, in the process of the experiment, we took the following measures: select the yeast in the logarithmic phase of growth for inoculation, increase the initial amount of inoculation, try to match the composition of the fermentation medium with the seed medium.

Modeling Analysis
I.Modeling of colony growth of Pichia pastoris
Using logistic function to predict the growth of Pichia pastoris

By the hypothesis above, we can assume that Pichia pastoris and yeast follow logistic growth.

In the equation, the concentration of Pichia pastoris should not exceed the carrying capacity of the bioreactor, K, which is assumed to be non-time-varying. μP-S is specific growth rate.dP represents the specific cell death rate of Pichia pastoris.

2. Deriving specific growth rate based on Monod kinetics

Specific growth rate, μ, affected by the limited substrate and the fluctuation of the pH, is expressed using Monod kinetics.
a. Specific growth rate with the concentration of substrate
By assuming the total substrate uptake follows Michaelis-Menten kinetics and the substrate is used for non-growth associated metabolic maintenance during stationary phase and microbial growth, the substrate uptake rate rs is expressed as:
In which KmaxP-S is the maximumsubstrate uptake rate by cells. S represents substrate concentration, P represents the concentration of Pichia pastoris,YP-S is the yield coefficient of the biomass to the substrate consumed, and KP-S is the substrate affinity. With constant maintenance requirement assumed, m could be regarded as a maintenance constant.
Therefore, the specific growth rate can be expressed as:

And assuming S >> KP-S the maximum growth rate should be:

And the substrate affinity can be expressed as:

Thus, the specific growth rate can be reformulated into:

The substrate consumption rate becomes:

b.The influence of the fluctuation of the pH

The pH of a culture may influence both the maximum growth rate and the inhibitory potentials of substances, particularly those of organic acids. To account for the effect of pH on μP-Smax the following expression is used:


where H+ is the hydrogen ion concentration;μP-Smax, KH, and KOH are constants.

3. The concentration of substrate

In order to keep the system working properly, we should keep adding substrates. In our system, we use formaldehyde as the carbon source of Pichia pastoris. For the sake of simple treatment, we assume that the addition of the substrate changes with time is a constant. That is:


kS-add is a constant.

Therefore, the concentration of the substrate can be expressed as



Combined with the equations above, the differential equation of substrate concentration according to time is:


II.Modeling of single Cell Gene Pathway
In Pichia pastoris, three heterologous genes are expressed:PelA, SLAC and VP. First we model their expressions separately. In addition, the expression of alkaline protein in Pichia pastoris is beneficial to regulate the low pH of the system and maintain the stability of the system.
During the growth of Pichia pastoris, the number of Pichia pastoris can be calculated by the following equation:


Where NumP is the number of Pichia pastoris, P is the concentration of Pichia pastori, V is the volume of the bioreactor, ρis the mass of a Pichia pastori.

  1. Kinetic Model of PelA Gene expression

The kinetic equation of PelA gene transcription is:


where trcPelA is the transcription rate of PelA,degPelA-mRNA is the degradation rate of mRNAPelA.

The last part of this equation takes into account the effect of pH on gene expression. KH-PelA , KOH-PelA and nPelA are constants.In view of the fact that the expression of different genes varies from the change of pH, add parameters(ngene) to adjust the equation.



  2. Kinetic Model of SLAC Gene expression

The kinetic equation of SLAC gene transcription is:


where trcSLAC is transcription rate of SLAC,degSLAC-mRNA is the degradation rate of mRNASLAC. As discussed above, KH-PelA , KOH-SLAC and nSLAC are constants .

The kinetic equation of SLAC gene transcription is:


where trlSLAC is translation rate of mRNASLAC,degSLAC-Pro is the degradation rate of ProSLAC.

  3. Kinetic Model of VP Gene expression

The kinetic equation of VP gene transcription is:


where trcVP is transcription rate of VP,degVP-mRNA is the degradation rate of mRNAVP. As discussed above, KH-VP , KOH-VP and nVP are constants .

The kinetic equation of VP gene transcription is:


  4. Kinetic Model of alkaline protein Gene expression

The kinetic equation of alkaline protein(abbreviated as a-Pro) gene transcription is:


where trca-Pro is transcription rate of a-Pro,dega-Pro-mRNA is the degradation rate of mRNAa-Pro. As discussed above, KH-a-Pro , KOH-a-Pro and na-Pro are constants .

The kinetic equation of a-Pro gene transcription is:


where trla-Pro is translation rate of mRNAa-Pro,dega-Pro is the degradation rate of Proa-Pro.

III. Modeling of the degradation of Banana components

The components of banana stem include cellulose, hemicellulose, lignin and pectin.The purpose of our project is to use the three enzymes heterogeneously expressed in Pichia pastoris to degrade lignin, pectin and hemicellulose to obtain higher quality fiber.
Simple Mass Action doesn't fit enzymatic reactions well. The reason for it is the fact that the enzyme and the substrate form a complex, which is then converted to a product and the original enzyme.In the process of substrate degradation, the degradation rate of enzyme satisfies Michaelis-Menten equation:



It has several parameters including the maximal reaction rate vmax , the Michaelis constant Km and the turnover number kcat.

  1. Kinetic Model of degradation of pectin and galacturonic acid

In the process of pectin degradation, the vast majority of the degradation product is galacturonic acid(abbreviated as GA), which lead to the decrease of pH. The equation of GA is as follows:


  2. Kinetic Model of degradation of lignin



  3. Kinetic Model of degradation of hemicellulose



IV.Regulation of pH system

In our project, pH, as an important environmental variable, affects many processes, such as gene expression, enzyme activity, yeast growth, and so on.

  1. Effect of pH on enzyme activity

The effect of acid H+ ions or basic OH− ions on the activity of an enzyme is probably caused by a change in stereo configuration at or in the neighborhood of the active sites (Fersht, 1984; Whitaker, 1994). As in almost all protonation reactions, these reactions occur very rapidly. The different configurations are instantaneously in equilibrium. The protonated and hydroxylated enzymes are assumed to be completely inactive or at least less active. This can be represented by the following mechanisms:

where KEH and KEOH are the equilibrium constants of the reactions.

The water dissociation is defined, as usual, as:


The amount of EnH and EnOH can now be expressed in terms of actual amount of active enzyme and pH by:


The total amount of enzyme in any configuration has to remain constant:


` According to the equation above, solving for En, an expression is obtained for the active enzyme at any H+ concentration (or pH):


Again, this model, for the amount of available active enzyme configuration, can be converted into an apparent activity (Act), as for the temperature model. This results in:


Again, in this equation, k s and only appear in combination with each other. It is therefore impossible to estimate both variables at the same time. Both parameters are therefore combined in a new parameter called Act 0 . This results in the final equation:



The equation describes how the activity of the total pool of enzyme changes with H + concentration over the entire range from pH0 up to pH14.

  2. Effect of pH on the growth of Pichia pastoris

The pH effect on bacterial growth may be directly attributed to the H+ in the medium. The hydrogen ion has been described as a noncompetitive inhibitor to the growth of some microorganisms at acidic pH. The hydrogen ions may inhibit one metabolic reaction that is the rate-limiting step and thus limit the growth rate of the culture. In general, the effect of H+ on the specific growth rate may be represented as follows:


where H+ is the hydrogen ion concentration;μP-Smax*, KH, and KOH are constants.

  3. Effect of pH on gene expression

The kinetic equation of gene transcription is:


where trcgene is transcription rate of gene,deggene-mRNA is the degradation rate of mRNAgene.
The last part of this equation takes into account the effect of pH on gene expression. KH-gene , KOH-gene and ngene are constants.
In view of the fact that the expression of different genes varies with the change of pH, add parameters(ngene) to adjust the equation.
  4.Calculation of hydrogen ion concentration

In this system, there are several reactions that are directly related to the concentration of H+.They are:


Therefore, the differential equation can be written as :


References:
[1]CHEN Letian, WANG Huiting, HAN Jingluan, et al. Research progress and perspective of plant pectin lysase[J]. Journal of South China Agricultural University, 2019, 40(5): 71-77.

[2]A.‐P. Zeng, A. Ross, H. Biebl, C. Tag, B. Günzel, W.‐D. Deckwer. Multiple product inhibition and growth modeling of clostridium butyricum and klebsiella pneumoniae in glycerol fermentation.Biotechnology and Bioengineering, Vol.44, Pp.902-911(1994).

[3]Bailey J E, Ollis D F. Biochemical Engineering Fundamentals.[J]. Chemical Engineering Education, 1976, 134(3):N/A.

[4]Costa J M . Growth Inhibition Kinetics for the Acetone-Butanol Fermentation[C]// Acs Symposium Series. 1983.

[5]Widder S , Schicho J , Schuster P . Dynamic patterns of gene regulation I: Simple two-gene systems[J]. Journal of Theoretical Biology, 2007, 246(3):395-419.

[6]Tijskens L M M , Greiner R , Biekman E S A , et al. Modeling the effect of temperature and pH on activity of enzymes: The case of phytases[J]. Biotechnology and Bioengineering, 2001, 72(3):323-330.

[7]John G T , Heinzle E . Quantitative screening method for hydrolases in microplates using pH indicators: Determination of kinetic parameters by dynamic pH monitoring[J]. Biotechnology and Bioengineering, 2001, 72(6):620-627.

[8]Ryhiner G B , Heinzle E , Dunn I J . Modeling and Simulation of Anaerobic Wastewater Treatment and Its Application to Control Design: Case Whey[J]. Biotechnology Progress, 1993, 9(3):332-343.

banana-tree
3-D Structural Model
We constructed the three-dimensional model of three kinds of enzymes, and calculated the Connolly surface of the enzyme and His-tag.Using the 3-D structures and the coordinate system of VP,PelA,sLAC gotton from webserver I-TASSER, the software Materials Studios can present the surface of the three enzyme.

1.Why do we construct the protein 3-D structures

His-tag can be added to the C-terminal or the N-terminal of the peptide, so that our model can help us with the following aspects:

-to predict where the His-tag should be attached to
-learn more details about these enzymes.
-predict the the possibility of the chelation
2. I-TASSER

I-TASSER (Iterative Threading ASSEmbly Refinement) is a hierarchical approach to protein structure and function prediction.

The server first tries to retrieve template proteins of similar folds (or super-secondary structures) from the PDB library by LOMETS, a locally installed meta-threading approach.

In the second step, by using replica-exchange Monte Carlo simulations with the threading unaligned regions, it constructed a thorough model. The low free-energy states are identified by SPICKER through clustering the simulation decoys.

In the third step, the fragment assembly simulation is performed again starting from the SPICKER cluster centroids, where the spatial restrains collected from both the LOMETS templates and the PDB structures by TM-align are used to guide the simulations. The purpose of the second iteration is to remove the steric clash as well as to refine the global topology of the cluster centroids. The decoys generated in the second simulations are then clustered and the lowest energy structures are selected. The final full-atomic models are obtained by REMO which builds the atomic details from the selected I-TASSER decoys through the optimization of the hydrogen-bonding network

For predicting the biological function of the protein, the I-TASSER server matches the predicted 3D models to the proteins in 3 independent libraries which consist of proteins of known enzyme classification (EC) number, gene ontology (GO) vocabulary, and ligand-binding sites. The final results of function predictions are deduced from the consensus of top structural matches with the function scores calculated based on the confidence score of the I-TASSER structural models, the structural similarity between model and templates as evaluated by TM-score, and the sequence identity in the structurally aligned regions.

3.The structures of the three enzymes with different positions of His-tag

4.Analysis the possibility level of the chelation mathematical

5.Conclusion
Form the analysis above, we can draw the conclusion as below:

For VP with His-tag in the C-terminal, the ratio of surface of His-tag to the whole protein is 0.01479, the other is 0.01062.

For PelA with His-tag in the C-terminal, the ratio of surface of His-tag to the whole protein is 0.01975 the other is 0.01872.

For SLAC with His-tag in the C-terminal, the ratio of surface of His-tag to the whole protein is 0.03589 the other is 0.03568.

For all three enzymes, with His-tag attached to the C-terminal the ratio is higher so we can regard that put the His-tag in the end of the enzyme will be more helpful. Therefor, we suggest that His-tag should be added to the C-terminal of the proteins.

6.References

1.J Yang, Y Zhang. I-TASSER server: new development for protein structure and function predictions, Nucleic Acids Research, 43: W174-W181, 2015.

2.C Zhang, PL Freddolino, Y Zhang. COFACTOR: improved protein function prediction by combining structure, sequence and protein–protein interaction information. Nucleic Acids Research, 45: W291-W299, 2017.

3. Jmol: an open-source Java viewer for chemical structures in 3D.