<!DOCTYPE html>
MODEL
Overview
The models we built included four parts. First, we established a fluorescent protein model to screen out the most suitable fluorescent protein, the main modeling method here is grayscale analysis. Then, for the large amount of measured OD values, we drew the growth curve of yeasts and it fitted logistic model. It described the growth situation of the yeasts after plasmid introduction, and we compare it with yeasts without any foreign plasmid. The growth curve also offers the best measuring point and the best measuring interval. What’s more, we drew the degradation curve of the fluorescent protein, which helps us know different characteristics of the two chosen fluorescent proteins better. Finally, we constructed a model to illustrate the oscillation of KaiA, KaiB and KaiC protein called Mars Model, it explained the reason why the cycle reduced in yeasts nicely. Modeling work integrated with experiments tightly made our project complete and convincing.
As a lot of fluorescent proteins and luciferase existing, an evaluation model to evaluate the performance of every fluorescent protein was established to find the most suitable one.
relationship  a_{i}=a_{j}  a_{i}>a_{j}(a little)  a_{i}>a_{j}  a_{i}>a_{j}(obviously)  a_{i}>>a_{j} 
score  1  3  5  7  9 
There were five factors we should consider. To determine the weight of them, we discussed with Prof. Zhao Kun, school of chemical engineering and technology, Tianjin University and professors in Tianjin Institute of Industrial Biotechnology, Chinese Academy of Sciences, they gave us a suggested order as shown below.
Between every two standards, select numerical value 1, 3, 5, 7 to measure the relative relation between index.
Stokes  QY*Brightness  Bleaching time  ph  Maturing time  
Stokes  1  3  5  1/3  5 
QY*Brightness  1/3  1  3  1/5  3 
Bleaching time  1/5  1/3  1  1/7  1 
ph  3  5  7  1  7 
Maturing time  1/5  1/3  1  1/7  1 
Construct standard matrix shown as below.
$$A = \begin{bmatrix} 1 & 3 & 5 & 1/3 & 5 \\ 1/3 & 1 & 3 & 1/5 & 3 \\ 1/5 & 1/3 & 1 & 1/7 & 1 \\ 3 & 5 & 7 & 1 & 7 \\ 1/5 & 1/3 & 1 & 1/7 & 1 \end{bmatrix}$$
Next, calculate the concordance index (CI) to judge whether the matrix is reasonable.
$$CI = {\lambda_{max}n \over n1} = 0.0340$$
$$\lambda_{max} = 5.1361$$
λ_{max}  the largest eigenvalue of the evaluation matrix and n is the dimension of the matrix.
According to the dimension of the matrix, RI value could be found from the form below.
n 
1 
2 
3 
4 
5 
6 
7 
8 
9 
RI 
0 
0 
0.58 
0.90 
1.12 
1.24 
1.32 
1.41 
1.45 
In our model, because n=5, RI=1.12.
$$CR = {CI \over RI} = 0.0304$$
Because CI<0.1, the matrix passes the consistency check. Therefore, the method we used to find the weight above is suitable.
Calculate the corresponding eigenvector of λ_{max} and normalize it, W=[0.0762 0.1607 0.3621 0.0389 0.3621].
At last, using the method mentioned in the [1], the score table is shown as below.
Score table  

score  0  25  50  75  100 
Stokes  0  00.02  0.020.05  0.050.1  0.1 
QY*Brightness  0  00.01  0.010.04  0.040.06  0.06 
Bleaching time  0  00.01  0.010.10  0.100.20  0.20 
pH  0  1.52.5 8.59.5  2.53.5 7.58.5  3.54.5 6.57.5  4.56.5 
Maturing time  0  00.01  0.010.05  0.050.30  0.30 
The score of every fluorescent protein can be obtained shown as Figure1.
Figure1 Selection of report genes
Figure2 Fluorescent proteins could be chose from
It shows that mCherry, mOrange and EYFP ranks top 3. However, after our preexperiment, it was found that the signal of mOrange was too weak. The reason maybe was that we sequenced the part of mOrange and found that there were some mutations. Therefore, we finally chose mCherry(66.0336 ranking 12/478) and EYFP(34.9807 ranking 82/478) as our report genes.
The degradation of the fluorescent protein is also of great importance to our experiment. Therefore, the change of fluorescence intensity with time was measured, and here are the results.
Figure3 EYFP Degradation Curve
Figure4 mCherry Degradation Curve
To better explain the degradation of the fluorescent protein, we consulted a large number of documents^{[2,3,4]}. It shows that the degradation of the fluorescent protein is exponential.
The EYFP degradation curve can be described by function (1):
$$f(x) = a·e^{bx} + c·e^{dx} (1) $$
Coefficients (with 95% confidence bounds):
a = 133, b =0.005066, c =44.38, d =0.02168
Figure5 Fitted EYFP Degradation Curve
Goodness of fit:
SSE: 76.77
Rsquare: 0.9443
Adjusted Rsquare: 0.9363
RMSE: 1.912
The mCherry degradation curve also can be described by function (1), but the coefficients are different.
Coefficients (with 95% confidence bounds):
a = 613.8, b = 0.0003886, c = 0.0003207, d =0.06852
Figure6 Fitted mCherry Degradation Curve
Goodness of fit:
SSE: 82.27
Rsquare: 0.9482
Adjusted Rsquare: 0.9404
RMSE: 2.028
The fitted degradation curve of mCherry approximates a straight line, so it can be simplified as a linear function:
$$ f(x) = p_1·x + p_2 (2)$$
Coefficients (with 95% confidence bounds):
$$ p_1 = 0.2287 , p_2 = 613.6$$
Figure7 Linear mCherry Degradation Curve
Goodness of fit:
SSE: 82.67
Rsquare: 0.9479
Adjusted Rsquare: 0.9455
RMSE: 1.938
It can be seen that the degradation rate of different fluorescent protein is very different, which made it clearer for us to know about the characteristics of the two fluorescent protein and offered important information to our experiment.
For the OD_{600} values we got, we did some processing and modeling work. And here are our steps and results.
There were three groups in our experiment. They were blank control group, partial control group and experimental group. After getting all the data, first, we drew a histogram and a scattergram of time and maximum OD_{600} values (Figure8, 9). These results were very instructive to experiments that these results told us the best measuring point and the best measuring interval.
Figure8 histogram of TimeMaximum OD Value
Figure9 Scatter gram of TimeMaximum OD Value
From the beginning to the maximum OD_{600} value, it fits the logistic model. The block effect of resource and environment for the growth of yeasts is reflected in the growth rate r, which makes r decrease with the increase in the number of yeasts x. Express r as a function r(x) of x, and take a simple and convenient linear reduction function r(x)=a+bx. In order to give a real meaning to the coefficients a and b in the growth rate function, we introduced two parameters:
(1)Intrinsic growth rate r : r is the growth rate when x=0 (in theory);
(2)Population capacity x_{m} : x_{m} is the largest yeast amount that can be accommodated by resources and the When x=x_{m}, the quantity of yeasts is no longer increasing, that is r(x_{m})=r+bx_{m}=0, then b=r/x_{m}.
r and x_{m} values in our experiments are shown in the chart below.
 YPD  SC  BY4741  dthree 
r 

 0.0164  0.0172 
x_{m} 

 0.8523  0.8034 
 pABaC＋p1m  pABaC＋p1E  pABaC＋p2N  pABaC＋p1F 
r  0.002364402  0.001746617  0.002826764  0.001905785 
x_{m}  0.402523944  0.508816901  0.424323944  0.542298592 
 pCiRbS＋p1m  pCiRbS＋p2N  pCiRbS＋p1F  pbCiRS＋p1m 
r  0.006367923  0.007098618  0.007176452  0.007853975 
x_{m}  0.410507042  0.254873239  0.446169014  0.315098592 
 pbCiRS＋p1E  pbCiRS＋p2N  pbCiRS＋p1F  pABaC+pCiRbS＋p1m 
r  0.024143608  0.012145451  0.002428334  0.006280764 
x_{m}  0.413985915  0.458239437  0.270442254  0.337278873 
 pABaC+pCiRbS+p1E  pABaC+pCiRbS＋p2N  pABaC+pCiRbS＋p1F  pABaC+pbCiRS＋p1m 
r  0.002305512  0.00217225  0.002272595  0.002039534 
x_{m}  0.33171831  0.293661972  0.303701408  0.289346479 
 pABaC+pbCiRS+P1e  pABaC+pbCiRS＋p2N  pABaC+pbCiRS＋p1F 

r  0.001894111  0.003848457  0.007151104 

x_{m}  0.301574648  0.345819718  0.329769014 
The resulting growth rate function is $$r(x) = r(1  {x \over x_m}) $$ Replacing intrinsic growth rate with r(x), get $${dx \over dt} = rx(1  {x \over x_m}) , x(0) = x_0 (3)$$
Factor rx in the function shows the growth trend of yeast amount itself, while factor (1  x/x_{m}) reflects the block effects of resources and environment to the yeast quantity growth. Obviously, the bigger x is, the former factor is bigger and the latter factor is smaller. The growth of yeast amount is the result of the two factors.
Take x as the horizontal axis and dx/dt as the vertical axis, we obtained a parabola (Figure 10), when x = x_{m}/2, dx/dt reaches the maximum. As shown in Figure 10, dx/dt changes with the increasing x, and we can do the following analysis to the curve x(t).
Setting when t = 0 x_{0} < x_{m}/2, with the increase of t, dx/dt increases, so x is growing faster and faster and the curve x(t) is raised downward; dx/dt decreases when x_{0} > x_{m}/2, x grows slower and slower and the curve x(t) is raised upward. x = x_{m}/2 is the inflection point of the curve. When x → x_{m} , dx/dt → 0, so x = x_{m} is the asymptote of x(t). From the above analysis, we can draw the figure as shown in Figure 11.
Actually, equation (1) can be solved by separation of a variable method as
$$x(t) = {x_m \over {1+({x_m \over x_0}1)e^{rt}}} (4) $$
Figure10 example xdx/dt curve
Figure11 example tx curve
By Analyzing our data, following results were obtained similar to the above. Figure12,13 shows the tx curve and xdx/dt curve of our blank control group, partial control group and experimental group. To make the results clearer, we enlarged two of them.
Figure12 tx curve
Figure13 xdx/dt curve
Model Construction
Oscillation in KaiC phosphorylation is the bestobserved parameter in this system and represents a key state variable for the clock in vivo. Thus we have sought to closely mimic this output in our project. Nakajima et al. ^{[6]} suggest, given the dual function of KaiC and ‘‘cooperation between KaiA and KaiB,’’ that autonomous oscillation of KaiC phosphorylation might be achieved. We established a model based on known biological and biochemical observations and our experiments that did not involve transcription or translation. In Figure14, we summarized the key steps of three Kai proteins oscillation when ATP is provided in excess. It was well established that we used three circles to represent all possible combinations of three Kai proteins, just like Mars and its two satellites. This was also why we call it Mars Model.
The model (Figure 14) contained twelve processes (R1R12) describing all the proteinprotein interactions and phosphorylationdephosphorylation between the Kai proteins. KaiXY represents KaiX and KaiY compound and KaiC* represents fully phosphorylated KaiC. Process R1, R2 and R3 are six aggregations of KaiC protein, two aggregations of KaiA protein and four aggregations of KaiB protein respectively. In process R4, KaiC_{6} binds KaiA_{2}, forming KaiA_{2}C_{6} compound. Since KaiA_{2} facilitates the autokinase activity of KaiC_{6}, KaiA_{2}C_{6} first converts to partial phosphorylated form, KaiA_{2}C_{6}C_{6}*, by process R5, and then rapidly converts to fully phosphorylated form, KaiA_{2}C_{6}*, by process R10. Then, fully phosphorylated protein KaiA_{2}C_{6}* combines with KaiB_{4}, forming KaiA_{2}B_{4}C_{6}*, by process R6. In process R7, KaiA_{2} is displaced from KaiA_{2}B_{4}C_{6}*. When KaiA_{2} no longer exists in KaiA_{2}B_{4}C_{6}*, KaiB_{4} dissociates from KaiB_{4}C_{6}*, by process R8. Process R9, R11, and R12 are depolymerization of KaiC_{6}, KaiA_{2} and KaiB_{4} protein, respectively^{[7]}.
Figure14 A dynamic model of KaiABC proteins oscillation.See text for description
Figure15 Rate equations of every reaction
Figure16 Oscillatory curve of every protein
The effect of temperature change
(Note: All the blue curves represent the initial data and red curves represent the revised data.)
(a)Period shortens with temperature rising temperature falling
(b) Period prolongs with
Figure17 Period changes with temperature changing(k changes proportionally)
(a) Period shortens with temperature rising temperature falling
(b) Period prolongs with
Figure18 Period changes with temperature change(k changes slightly)
When the reaction rate constants change slightly with temperature changing, the period shortens while amplitude shortens too. Therefore, if k changes disproportionately, when the temperature increases, the cycle is shortened and the oscillation is unsteady. The curve tends to be gentle with time, which means the oscillation disappears shown asFigure 19.
Figure19 The disappearance of oscillation with temperature changing
The effect of phosphorylation rate
The involved process of phosphorylation is R4 and the relating reaction rate constant isk_{4}.
Phosphorylase in yeasts may have promoting effect to the phosphorylation of protein and yeasts offer enough ATP/ADP in vivo, which increase the rate of phosphorylation. Therefore, k_{4}may increases in yeasts, which makes oscillation cycle shortens shown as Figure20.
Figure20 Period changes with phosphorylation rate changing
(a)Period change with the supply rate of KaiA increasing
(b)Period change with the supply rate of KaiB increasing
(c)Period change with the supply rate of KaiC increasing
Figure21 Period changes with oscillation environment changing
Analysis of other process
The above analysis involves some reaction rate constants that has real biological meaning. Besides, other reaction rate constants were also analyzed and we found them had an impact to the oscillation.
In the following table, we list the impact of k_{5}k_{12} changing to the oscillation. It can be seen in Figure22 that k_{5} and k_{6} have great influence on the disappearance of the oscillation. The impact of k_{5}k_{8} may relates to phosphorylation and temperature changing as well as other factors. Therefore, figuring out biological factors related to these reaction rate constants is one of our future work.
The impact of k_{5}k_{12} changing  

k_{5}  Oscillation reduces largely with it increasing slightly; oscillation disappears with it increasing by order of magnitude. 
k_{6}  The similar effect as k_{5}. 
k_{7}  Oscillation reduces slightly with it increasing. 
k_{8}  The similar effect as k_{7} 
k_{9}  Moderate impact on the oscillation. Oscillation disappears with it increasing by one order of magnitude. 
k_{10}  Has little impact on the oscillation. Phase changes with it increasing. 
k_{11}  The oscillation is fine when k_{11} <0.1; KaiA and KaiB curve oscillates when the order of magnitude is 10^{2}; KaiA and KaiB oscillation disappear gradually with the order of magnitude decreasing 
k_{12}  The similar effect as k_{11} 
Figure22 Period changes with k_{5} k_{10} changing
Figure23 Period changes with k_{11} and k_{12} changing
Summary
We propose some assumptions in combination with the model based on our experimental results.
The following are the three main factors:
(1)Temperature: There are two possible patterns that temperature may effects. If the reaction rate constants k changes proportionately, the period will shorten. If the reaction rate constants k changes slightly and equally, the period will shorten too while the oscillation disappears gradually.
(2)Phosphorylation: We have demonstrated through model tests that the oscillations decay rapidly with accelerated phosphorylation. Therefore, we have two conjectures: one is that phosphorylase in yeast plays a better role in promoting phosphorylation, and the other is that yeast provides sufficient or excess ATP/ADP to accelerate the rate of phosphorylation. This direction is also the focus of our further research.
(3)Concentrations of KaiA, KaiB and KaiC: Unlike the envisaged results, the concentrations of KaiA, KaiB and KaiC did not have much effect on our model in the testing of mathematical models.
Besides the above three factors, other reaction rate constants were also analyzed and their biological significance needs to be figured out in the future.