Team:USTC/Model

---

Overview

Why do we model? What do we want to achieve?
Our project consists of two parts: gene circuit and electron transfer. The goal of our modeling is to optimize the electron transfer rate, which is related to the amounts of 4 kinds of proteins: CymA, MtrA, MtrB and MtrC. Therefore, we construct the gene circuit model to calculate the amounts of 4 proteins in each Shewanella oneidensis MR-1 cell. Then, we build the electron transfer model to simulate the process of electron transfer in our system and calculate the electron transfer rate. We can thus analyze the relationship between electron transfer rate and the amounts of different proteins. This will guide our experiment design to optimize the electron transfer rate.

Gene circuit

Gene circuit model is composed of 8 nonlinear ordinary differential equations (ODEs). We use SimBiology Toolbox in MATLAB to solve these equations and then obtain the amount of each protein. By modifying the parameters in the equations, we can also know how to control the amount of each protein.

Electron transfer

As for electron transfer, our model is an approximate simulation to compute the possibility of electrons' transfer. We consider the proteins' motions as Brownian motion. The electron transfer happens when different proteins get close enough to each other. The movement of proteins can be described as a function of possibilities according to Markov pathway. By this method, we successfully simulate and predict the electron transfer rate. Besides, we analyze how the transfer rate would change when the amount of certain proteins double and explore the reasons. This is of guiding significance when WetLab considering gene circuit modifying to maximize electron transfer rate.

Gene circuit

Introduction

Since the WetLab have designed the circuits, it is natural for us to think of these questions: how the circuit works and how much protein can be produced from the gene circuit. Thus, we apply nonlinear ordinary differential equations (ODEs) to build gene circuit model. Gene circuit model can give out the concentration of each type of proteins in the model and predict their amounts in a single S.oneidensis MR-1 cell.
Additionally, this model helps us to understand the logical relationship between the promoter, terminator, protein, inducer, and so on of the genetic circuits. Furthermore, this model will guide us on controlling the amount of proteins.

Procedure

We use Simbiology toolbox in MATLAB to help us build this model. We establish four genetic circuits, which are responsible for producing MtrA, MtrB, MtrC and CymA proteins. (Fig 1)

Fig 1. The relating gene circuit

The change of mRNA concentration can be divided into two components: transcription from DNA and self-degradation. The change of protein concentration can also be divided into two components: translation from mRNA and self-degradation. (Fig 2.a) To simulate this process, we use the ordinary differential equations as showed in Fig 2.b to build our model. Most of the kinetic parameters in the equations are from literatures, but some without experimental data are estimated. The model is simulated with MATLAB solvers to get the change of every component.

(a)

(b)

Fig 2. (a) A diagram showing the process of how the amounts of mRNAs and proteins change. (b) The ordinary differential equations describing the process of transcription, translation and degradation

Results

Since transcription and translation processes proceed in a low speed and the degradation of mRNA and protein are also going on simultaneously, it takes a long time for all macromolecule quantities to stabilize. Therefore, we set the termination time to be 40000 seconds. By solving the ODEs, the figures we obtained of the concentration changes over time is shown below.

(a)

(b)

(c)

(d)

Fig 3. Outcomes of our simulations showing the change of concentrations over time, (a) shows the change of MtrA, (b) of MtrB, (c) of MtrC and (d) of CymA. As is shown, all concentrations have a steep increase at the beginning and reach a rather stable stage over a long period of time. In all cases, the concentrations of DNA and RNA are so small compared with those of proteins as their curves are on the very bottom of the figures.

Conclusions

It can be seen from the above curves that all the gene circuits work well and all the proteins can be expressed successfully. We can also get the final stable quantities of all the proteins, which are shown in the table below.

protein	MtrA	MtrB	MtrC	CymA
amount	9.702E4	2.206E4	2.298E4	8.215E4

Table 1. the amount of proteins when the system reaches a stable state. It fits well with the fact that MtrB and MtrC formed a complex as their amounts are roughly the same.

These data will be applied to our second model, the electron transfer model. The modeling results are in good agreement with the actual situation because the amounts of MtrB and MtrC are almost the same for they are bound together on the cell membrane.

Electron transfer

Introduction

To quantitatively analyze the electron transfer rate, we apply a model based on simple symmetric random walk process, in which the membrane of the cell is approximately seen as an infinite flat plane and MtrA, MtrB, MtrC and CymA are seen as particles on it with certain densities. The density can be calculated using the outcome of the former model. As MtrB and MtrC form a binding, we only take the motions of MtrA, MtrB and CymA into consideration, while the motion of bound MtrC is considered exactly the same as MtrB while the motion of unbound MtrC have little influence on the results. Electron transfer occurs in cells only when MtrA collides with a combination of MtrC and MtrB. Additionally, CymA is indispensable to the process of electron transfer. Therefore, we build our model in the following way. The membrane plane is divided into cubic lattices, then one electron can be transferred from one to another protein only when they move into a vertical line, or into the lattices with the same x and y coordinates. To run our model, the interactions of proteins like physically collision and electrostatic force are ignored, and each protein can move into any of its eight immediate neighbor lattices of the same 2-D plane in different directions in each simulation step with equal possibility. Under these assumptions, we use MATLAB to find the optimized transfer rate with five different concentrations of the proteins. The WetLab could try to increase the concentration of the protein as our optimal solution suggests.

Procedure and results

The local structure of the membrane is shown in Fig 4. We can see that CymA is binding on the intracellular membrane, and MtrC is binding on MtrB which is binding on the extracellular membrane, MtrA in the periplasm between the two layers. For the proteins are on a much smaller scale than the cell membrane, the membranes can naturally be seen as flat planes. As the thickness of periplasm is tiny in contrast with the membrane, the periplasm can be considered as another flat plane in this model. Thus the proteins are assumed to move freely on three 2-D planes respectively, their motions seen as Brownian Motion, which can be described by Markov pathways.

Fig 4. This figure shows the complex process of electron transfer. Fig 4. This figure shows the complex process of electron transfer. We are only interested in the process related to CymA, MtrA, MtrB and MtrC, and their positions described above can be seen in the middle part of this picture.

The two transfer processes, from CymA to MtrA and from MtrA to MtrB/MtrC, are counted respectively in our model and the actual number of transferred electrons from intracellular to extracellular takes the smaller of the two results.
A single S. oneidensis MR-1 cell can be seen as a column with the height of 2μm and the diameter of 0.4μm, its surface area then comes as 0.88 μm^2. In order to reduce the amount of our calculation, we only consider a small area with the size of 20nm*20nm, the number of MtrA, MtrB and CymA is then given as 15, 5 and 13 by calculating the densities. The area is meshed into 200*200 grinds. To imitate an infinite plane, when a particle reached the boundary and is about to move out of this area in the next step, it is posted on the symmetric edge. Then we write a MATLAB code to describe the random walk process and add up the total time of MtrA, MtrB and CymA overlapping each other in 8000 steps. The process is repeated 80 times, each time the number of electrons been transferred and the average current are given out and plotted. Our MATLAB script can be seen here. Then we wonder how protein concentrations can affect the average current. We run 4 other simulations by doubling the amount of MtrA, doubling the amount of MtrB, doubling the amount of CymA and doubling the amount of all four proteins. The result of our 5 simulations are plotted below. The exact data are shown in Table 2.

(a)

(b)

(c)

(d)

(e)

Fig 5. The output of our simulation under the 5 conditions: the origin situation(a), doubling the amount of MtrB(b), doubling the amount of MtrA(c), doubling the amount of CymA(d) and doubling the amount of all four proteins(e). We can see that the overlapping time doubles as the related protein doubles.

Then we wonder how protein concentrations can affect the average current. We run 4 other simulations by doubling the amount of MtrA, doubling the amount of MtrB, doubling the amount of CymA and doubling the amount of all four proteins. The result of our 5 simulations are plotted below. The exact data are shown in Table 2.

	(a)	(b)	(c)	(d)	(e)
electrons from CymA to MtrA	79.1	147.9	152.8	75.8	294.9
electrons from MtrA to MtrB	185.1	183.5	378.2	348.1	708.4
averaging current	79.1	147.9	152.8	75.8	294.9

Table 2. the average overlapping times under 5 conditions. The number of electrons transferred outside the cell is no more than the smaller one of the two processes, so we take the smaller rate as the approximate current rate.

From our outcomes, we can draw that rise in the amount of proteins can help increase the rate of electron transfer process associated with it as that effectively increase the possibility of two proteins overlapping each other. The doubling of all proteins offers the biggest average current, almost 4 times of the premier rate. The doubling of MtrB and the doubling of MtrA gives out roughly the same average current, about 2 times of the premier rate, and the doubling of CymA brings about almost no difference.
Although the doubling of all proteins gives out the maximal transfer rate, the doubling of MtrB might be a better choice, because only increasing MtrB offers the largest increase in transfer rate per extra-produced protein, so the same increase in current consuming less energy for the cell’s vital movements.
However, if we want to further increase the electron transfer rate, only increasing MtrB is not a wise choice. Notice that the doubling of MtrB in this case just fills up the huge density gap between proteins in initial conditions. From the table we can conclude that the smaller the difference of the concentration of proteins, the less the overlapping times would be wasted. So the biggest bargain for WetLab might be trying to modify the gene circuit for an approximate equal yield for the four proteins in practice.

Parameters

parameter	description	value	unit	reference
t1	transcription rate of cymA gene	0.0886	𝑠-1	[1]af
t2	transcription rate of mtrA gene	0.05	𝑠-1	[1]af
t3	transcription rate of mtrB gene	0.0238	𝑠-1	[1]af
t4	transcription rate of mtrC gene	0.0248	𝑠-1	[1]af
p1	translation rate of cymA protein	0.0802	𝑠-1	[1]bg
p2	translation rate of mtrA protein	0.045	𝑠-1	[1]bg
p3	translation rate of mtrB protein	0.0215	𝑠-1	[1]bg
p4	translation rate of mtrC protein	0.0224	𝑠-1	[1]bg
f1	degradation rate of mRNA	0.0012	𝑠-1	[2]c
f2	degradation rate of mRNA	0.0012	𝑠-1	[2]c
f3	degradation rate of mRNA	0.0012	𝑠-1	[2]c
f4	degradation rate of mRNA	0.0012	𝑠-1	[2]c
d1	degradation rate of cymA protein	0.000193	𝑠-1	[2]d
d2	degradation rate of mtrA protein	0.000193	𝑠-1	[2]d
d2	degradation rate of mtrB protein	0.000193	𝑠-1	[2]d
d2	degradation rate of mtrC protein	0.000193	𝑠-1	[2]d
N	plasmid number	10	per cell	estimated

Description

a. The average transcription rate is about 50 nucleotide per second in Shewanella. The transcription rate of each transcript is calculated according to their lengths.
b. The average translation rate is about 15 amino acids per second in Shewanella. The translation rate of each protein is calculated according to their amino acid numbers.
c. The average half-life of mRNAs is set to 10 min. The degradation rate is calculated according to the half-life value.
d. The average half-life of proteins is set to 2 h. The degradation rate is calculated according to the half-life value.
f．the gene’s length is from Genbank database
g. the length of protein is from PDB database

Script

Mark FinalVersion.m

Reference

[1] 2018 iGEM USTC wiki
[2] Stamatakis, Michail, and Nikos V. Mantzaris. "Comparison of deterministic and stochastic models of the lac operon genetic network." Biophysical journal 96.3 (2009): 887-906.
[3] https://2017.igem.org/Team:USTC