Team:SJTU-BioX-Shanghai/modeling

Modeling

Our whole project consists mainly of two approaches to OA: the eukaryotic one and the prokaryotic one. Our mathematicians have built two models respectively for them.

Click to jump to:

• Eukaryotic part
• Prokaryotic part
• References

The eukaryotic approach

Click to jump to:

• Abstract
• Introduction
• Hypotheses
• Development
• Simulation
• Conclusion

Abstract

Back to top ↑

Our eukaryote circuit is designed to bring about the expression of the functional protein, both Oct4 and Col2, in response of the light-controlled calcium influx. It is undoubtedly challenging to model and predict the reliable operational outcomes of a dynamic system in the context of differentiating and morphologically complex eukaryote cells, however, promising and exciting at the same time. As is shown below, ordinary nonlinear differentiation equation is applied to illustrate the complex gene regulatory network, in which we determine the vital parameters and parameter spaces in which the system will function as intended and verify the robustness of the system in more realistic situations. Though given the largely uncharacterized infrastructure and potential perturbation, our model has demonstrated the feasibility in manipulating the eukaryote preexist pathways and expressing gene networks for the use of synthetic biology.

Introduction

Back to top ↑
First of all, our model is based on the cell signaling pathways which have been rigorously researched and are with known mechanisms. The one we use here as an exemplary is the Ca2+ /CAM/CAMK dependent signaling pathway, of which the complex network is well illustrated below. The gene regulatory network can be characterized via three pathways’ synergic effect on the gene expression with a cascade of interactions between protein modules involved. The key nodes of the network are CAM, MEF2-Cabin1 complex, and MEF2-NFAT complex.

Furthermore, the two main type of chemical reaction involved are the simple forward reaction and phosphorylation and dephosphorylation enzymatic reaction.



Reactions are as follows:

Assumptions and hypotheses

Back to top ↑
Complex and uncharacterized as the eukaryote gene network is, the use of assumptions and hypothesis for simplification has become a necessity for mathematical modeling. We assume that:
• the molecules involved in the reaction are sufficient to ensure the reaction to take place continuously and deterministically
• the activating of the protein and kinase, which is described in the form of complex formation, is treated as the same principle as the simple forward reaction
• the necessary components of the gene regulatory network that are not tissue specific are sufficient constitutively provided by the cell at a constant rate
• the components of the gene regulatory network that are tissue specific such as Cabin1 and that are not adequately provided such as MEF2 are produced under the genetic control of our self-designed promoter mainly
• either protein or mRNA are decreased though degradation, diffusion and growth dilution at the rate proportional to its concentration
• the phosphates are sufficient in the environment

Model development

Back to top ↑
To simplify our model, we eliminate the intermediary node – CAMKIV, which expects no influence on how our system works. The equations composed according to the assumptions above are as follows:



where:

• r_i (i=1,2,…,7) represents the production rates of the protein, in combination of the promoter strength of either the constitutive promoters of the cell or the self-designed promoters as well as the ribosome strength in the translation process
• k_i (i=1,2,…,6) represents the reaction coefficient of the simple forward reaction
• b_i (i=1,2,…,15) represents the decreased (degradation) coefficient of either the protein or the mRNA
• K_mi (i=1,2) represents the Michaelis constant of the enzyme-catalyzed reaction
• K_cati (i=1,2) represents the number of substrate molecules turned over per enzyme molecule per second
• e represents the signal-dependent gene expression coefficient
• H represents the Hill constant value of input that gives 50% response
• n represents the Hill coefficient slope of signal-response curve at this input signal

Simulation and sensitivity analysis

Back to top ↑
Values of parameters:
• We should admit that there are quite some parameters that are unavailable to us after searching for a variety of large data, because our lack of quantitative knowledge about eukaryote cells. We set value of some of the parameters to ‘1’ for normalization purpose, and at the same time others adjusted for reasonable order of magnitude in order to demonstrate the relative tendency of the reaction. The Km values are chosen according to the usual magnitude of enzyme-catalyzed reaction compared with the normalization result.

Parameter values
b1	1	b11	1	k5	1	r2	1
b2	1	b12	1	k6	1	r3	1
b3	1	b13	1	kcat1	40	r4	1
b4	1	b14	1	kcat2	40	r5	1
b5	1	b15	1	km1	1	r6	1
b6	1	b16	1	km2	1	r7	0.01
b7	1	k1	1	e	1	r8	1
b8	1	k2	1	n	10
b9	1	k3	10	H	unknown
b10	1	k4	1	r1	0

• We realized the simulation aided by MATLAB SimBiology package.

The most crucial curve in these reactions is the responsive curve of the transcription, which is demonstrated by the Hill curve, shown below.



In order to show the entire curve, we figured out the appropriate value of the unknown H and modified the b13, which is the degradation of the mRNA, into 2.5 after observing and characterizing our system using simulation. After parameters are fully defined, we determine the initial amount of the components in our system by fitting it into the first stable status it arrives at under these pre-specified parameters. This step is aimed at mimicking the usual resting status of the eukaryote cell in preparation of the disturbance of various amount of calcium. Let us take a look at how the system changes after calcium influxes at the rate of 1.









Now let us see the critical points in the amount of mRNA of the designated protein in response of the different calcium influx rate.

When there is almost no calcium (the initial state)



After calcium intriguing:

• When the influx rate of calcium reaches 0.6:



• When the influx rate of calcium reaches 0.9:



• When the influx rate of calcium reaches 2:



• When the influx rate of calcium reaches 5:



We can almost delineate the Hill curve of the mRNA with above five critical values.

Conclusion

Back to top ↑
As we can see from the above, the mRNA of the designated protein Col2/OCT4 is almost zero without intriguing. When we gradually increase the influx rate of the calcium, the transcription is not induced until the calcium reaches nearly 0.6. The transcription hits almost saturated point --- 0.4 when the amount reaches 2, and the midpoint arrives at the amount of 0.9. When the concentration of the complex NFAT-MEF2 formed by the induction of calcium influx rate is over the saturated point of the mRNA, there will be no change in the final production of mRNA of the protein Col2/OCT4.

---

The prokaryotic approach

Back to top ↑

Click to jump to:

• Abstract
• Introduction
• Hypotheses
• Development
• Simulation
• Conclusion

Abstract

Back to top ↑
Our prokaryote circuit is designed in the content of Ecoli, which consists of several subsystems. We have constructed the mathematical models to see how well the main components of the whole system interact and cooperate, proving the feasibility of our project, and predicting the concentration of our targeted protein.

Introduction

Back to top ↑
First of all, our model is based on a system induced by extracellular signaling and controlled by negative feedback of the quorum sensing. The interaction between the subsystems is depicted below. The main reactions involved are the complex formation of SoxR-NO and the induction of the transcription of the mRNA-Col2-OCT4. And the cell-cell communication integrates the entire population as an essential component of the population-control circuit.



Assumptions and hypothesis

Back to top ↑
• the molecules involved in the reaction are sufficient to ensure the reaction to take place continuously and deterministically
• any protein and its related mRNA are generated at a constant rate
• either protein or mRNA are decreased though degradation, diffusion and growth dilution at the rate proportional to its concentration

Model development

Back to top ↑
The equations composed according to the assumptions above are as follows:

Subsystem:

(latex)

where

• r_i (i=1,2) represents the strength of promoters
• t_i (i=1,2,3) represents the strength of the ribosome in translation of protein
• k_i (i=1) represents the reaction coefficient of the complex formation reaction
• b_i (i=1,2,…,6) represents the decreased (degradation) coefficient of either the protein or the mRNA
• e represents the signal-dependent gene expression coefficient
• H represents the Hill constant value of input that gives 50% response
• n represents the Hill coefficient slope of signal-response curve at this input signal

Considering that our interest lies in only the amount of the proteins rather than various mRNA, we decided to only focus on the protein interaction. Thus, we proposed the simplified version of our model in the following equations:

(latex)

where

the difference is r_i (_i=1,2) represents the combination of the strength of promoters and ribosome.

Supervisor:

We model the population-control circuit in a simplified way, demonstrated and experimentally validated by Lingchong, Y^[4].

Further assumptions:

• changes in viable cell density follow logistics kinetics
• the cell death rate is in direct proportional to the concentration of the killer protein - Barnase
• the production rate of the killer protein is proportional to AHL concentration (to be same inside and outside the cells)
• AHL production rate is proportional to number of cell



where:

• K represents the growth rate
• N_m represents the carrying capacity
• K_B represents the production rate of Barnase related to the concentration of AHL
• K_A represents the production rate of AHL to the number of cell
• d represents the death rate of the cell related to the amount of Barnase
• b_i (i=4,5) represents the decreased (degradation) coefficient of the protein

In order to make the protein concentration compromised with the negative feedback of the population supervisor, we revise the equation in the subsystem as follows:



Simulation

Back to top ↑

Parameters' Values
Param.	Value	Unit	Annotation
b1	1	h-1	degradation rate of the various proteins at their base values[4]‘1’s are adjusted for keeping the same order of magnitude as above
b2	1	h-1
b3	1	h-1
b4	2	h-1
b5	0.639	h-1
d	0.004	nM-1·h-1	death rate of the cell related to the amount of Barnase
e	5.00E-7	/	unknown, adjusted signal-dependent gene expression coefficient
h	1	/	Hill constant value of input that gives 50% response
n	10	/	Hill coefficient slope of signal-response curve for simulation
k	0.97	h-1	growth rate chosen at PH=7
K_A	4.80E-7	nM·ml·h-1	production rate of various protein at their base values[4]
K_B	5	h-1
k1	1	h-1	unknown, adjusted reaction coefficient of the complex formation reaction
Nm	1.24E+9	ml-1	carrying capacity at PH=7 under certain media composition
r1	5.00E-7	h-1	unknown, adjusted strength of promoters which are available to be chosen
r2	5.00E-7	h-1

As you can see, some parameters are unavailable though we have searched from a variety of databases. However, the approachable parameters can help us to keep a reasonable order of magnitude for the parameters for sake of system realization. We may pursue the estimation of these parameters depending on the external experiment data.

We realized the simulation aided by MATLAB SimBiology package.

First, the initialized state of the system is determined by the first steady state, in which most of the value is almost zero, except that the amount of NO is under control and remains constant throughout the simulation process.

Now take a glance at how the population-controlled system works.





As we expected, there are two steady states in the system, and the second state after some moderate oscillation is reached when b4+b5>b5>K, which is feasible in biological sense. Then we change the amount of the NO to see where the minimum and saturated points lie in the response of the protein regulated by the SoxR-NO complex, shown by five critical points contained in the Hill curve.

• When there is no NO present:



• When the amount of NO reaches 0.1:



• When the amount of NO reaches 0.2:



• When the amount of NO reaches 0.3:



• When the amount of NO reaches 1:



---

Conclusion

Back to top ↑
As we can see from the above, the designated protein Col2/OCT4 is expressed at a very low level by the constitutive promoter without being induced by the SoxR-NO complex. When we gradually add slight amount of NO, the transcription is not induced until the NO reaches nearly 0.1. The transcription hits almost saturated point when the amount reaches 0.3, and at the same time the complex concentration of SoxR-NO is quite small, which to some extent reflects that the transcription can be induced without consuming too much of the intermediate reagents. When the concentration of SoxR-NO complex formed after induction of NO is over the saturated point, though the the amount of SoxR-NO tends to be raised, there will be no change in the final production of protein Col2/OCT4.

It is really delightful to successfully validate the feasibility of the mathematical models we have constructed.

---

References:

Back to top ↑
[1] Hidde D.J., Modeling and simulation of genetic regulatory systems: a literature review, Journal of Computational Biology Volume 9, Number 1, 2002

[2] A. Polynikis, S.J Hogan & M.di Bernardo, Comparing different ODE modeling approaches for gene regulatory networks, Journal of Theoretical Biology 261, 2009

[3] Didier Gonze, Stochastic simulations application to molecular networks, March 27, 2007

[4] Lingchong Y, Robert S.C., Ron W & Frances H.A, Programmed population control by cell-cell communication and regulated killing, April 4, 2004