Team:USTC/Model
Modelling
Theory meets experiment
Leaders
Our modeling mainly consists of three parts. One of them models the chemical reaction in the system to show the effect of GdnHCl and heat shock on the expression level of final product, which is replaced and thus indicated by GFP level in our project. The second part focuses on binding reactions intermediated by prion, which is just sup35 in our experiment. Each part includes two gene circuits, binding of AD (activate domain) and BD (binding domain) of promoters (circuit 1), and binding of sfGFP1-10 and sfGFP11, two GFP splits (circuit 2). The last part models the application of our products in hardware.
In the first part, we vary the heat shock temperature and concentration of GdnHCl, together with the time point of heat shock and introduction of GdnHCl. It turns out that those time points almost have no effect. The concentration of GdnHCl has the boldest effect, and the temperature has an auxiliary effect. Therefore, we can quantitatively control the outcome of GFP (or other protein products) in both great and tiny extent. The effects in gene circuit 1 are more significant and robust than in gene circuit 2.
In the second part, the binding reactions intermediated by prions are discussed in details. We set the maximum possible number of prions in prion aggregation. Under this upper limit, all sizes of aggregation are allowed to exist. Based on the configuration and size of different molecules compared to prion, we calculate or set the minimal number of prions which can cover all the molecules bound on it originally, such as AD and BD. Varying the maximal size of aggregation and the rate of binding and decomposing reactions, we get the final GFP level in different conditions. Thus we verify the effect of heatshock and GdnHCl on expression level of final product from another aspect. In addition, there may be some doubts on whether AD and BD will be covered by prions and fail to function well. Here we have disproved this suspicion using our modeling result.
In the last part, our product bacteria solution is bound outside a tube as a sensor of the inner temperature of the tube. The tube can be fixed outside a water bath kettle, and its inner temperature can reflect the temperature in the kettle. With some signal processors, we can design an automatically controllable water bath kettle, whose inner temperature is kept steady by it own.
Introduction
In this part, we involve the process of heat shock in our simulation of this reaction system. Together with GdnHCl, we hope to use heat shock to express product proteins(GFP is used in our simulation) quantitatively. By varying the time of heat shock and introduction of GdnHCl, the temperature of heat shock and the concentration of GdnHCl, we change the final balanced state and thus the final level of GFP. The relationship between the final state and the process is also studied quantitatively. In order to get a good effect of controlling the final level of GFP, we conduct all the process at the period when the number of the yeast (S.cerevisiae) cells is almost a constant.
Modified modeling of the mechanism of Hsp
In all of our previous modeling, we didn't take Hsp into consideration when modeling the aggregation of PRs. This was reasonable because we didn't need to concern about the detail of the aggregation, for example, the influence of temperature. However, we must involve Hsp in aggregation reaction because we are looking into the heat shock, to which temperature, thus the detail of the reaction is very significant.
The mechanism of Hsp involves two aspects now. First, before introduction of GdnHCl, the Hsp acts as the molecular chaperone of aggregation of PRs. After introducing GdnHCl, a part of Hsp is converted to another kind of enzyme that catalyzes the deccomposition of PR aggregates.
Gene circuit 1(AD and BD)
1.Reactions
The gene circuit is shown below(Figure 1), which is the same as the one before. So description of it is omitted here.
Figure 1. Gene circuit 1
Because we make the simulation when the colony hardly grows, we made a little adjustment to the parameters used before, and they are listed here.
Figure 2. Parameters
Because we modify the mechanism of Hsp, so there are some changes about the reaction. Now the reactions in the system are listed below:
\[mRNA(PR-AD) \rightarrow PR-AD\] \[mRNA(PR-BD) \rightarrow PR-BD\] \[PR-AD+PR-BD+E(pre) \leftrightarrow CPX3+E(pre)\] \[CPX1+pGal(off) \rightarrow CPX1+pGal1(on)\] \[pGal1(on)\rightarrow pGal1(off)\] \[mRNA(GFP)\rightarrow GFP\] \[mRNA(Epre)\rightarrow E+GdnHCl\] \[E+CPX1\leftrightarrow CPX2\rightarrow PR-AD+PR-BD+E\]
There are 2 catalytic reactions in this system. One is aggregation of PRs and the other is the decomposition of the aggregates. This two catalytic reactions are regarded as Michaelis-Menton reactions, and the catalytic rate are set to 1000 folds and 10000 folds, respectively. E(pre) is Hsp before adding GdnHCl, and E is Hsp after introducing GdnHCl. Then the differential equations are:
\begin{align*} \frac{d[mRNA(PR-AD)]}{dt}=&N*k_{m}-d_{m}[mRNA(PR-AD)]\\\\ \frac{d[mRNA(PR-BD)]}{dt}=&N*k_{m}-d_{m}[mRNA(PR-BD)]\\\\ \frac{d[PR-AD]}{dt}=&k[mRNA(PR-AD)]-d(PR-AD)+M1* k_{f}[CPX2]\\&-M*k_{f}[E(pre)]* [PR-AD]*[PR-BD]\\\\ \frac{d[PR-BD]}{dt}=&k[mRNA(PR-BD)]-d(PR-BD)+M1* k_{f}[CPX2]\\&-M* k_{f}[E(pre)]*[PR-AD]*[PR-BD]\\\\ \frac{d[CPX3]}{dt}=&-d(CPX3)-M1* k_{f}[CPX1]*[E]+k_{b}*M1*[CPX2]\\&+k_{f}*M[CPX3]\\\\ \frac{d[CPX1]}{dt}=&-d(CPX1)-M1* k_{f}[CPX1]*[E]+k_{b}*M1*[CPX2]\\&+k_{f}*M[CPX3]\\\\ \frac{dp}{dt}=&k_{p1}*[CPX1]*(1-p)-k_{p2}*p\\\\ \frac{d[mRNA(GFP)]}{dt}=&N*k_{m}-d_{m}[mRNA(GFP)]\\\\ \frac{d[GFP]}{dt}=&k[mRNA(GFP)]-d_{GFP}[GFP]\\\\ \frac{mRNA(Epre)}{dt}=&N*K1(T)-K2(T)[mRNA(GFP)]\\\\ \frac{d[E(pre)]}{dt}=&k[mRNA(Epre)]-d[E(pre)]+M1*k_{f}[E(pre)]\\&*[GdnHCl]-M* k_{f}[E(pre)]*[PR-AD]*[PR-BD]+\\&M*k_{f}*[CPX3]\\\\ \frac{d[E]}{dt}=&-d[E]+M1*k_{f}[E(pre)]*[GdnHCl]-M*k_{f}[E]*[CPX1]\\&+M1* (k_{f}+k_{b}*[CPX2])\\\\ \frac{d[GdnHCl]}{dt}=&-d[GdnHCl]\\\\ \frac{d[CPX2]}{dt}=&-d[CPX2]+M1*k_{f}[E(pre)]*[GdnHCl]\\&-M1*(k_{f}+k_{b})* [CPX2]\\\\ \end{align*}
where M is the catalytic rate of aggregation, M1 is the catalytic rate of decomposition. M=1000 and M1=10000. K1(T) and K2(T) are functions of temperature.2.Introduction of GdnHcl
After introduced to the system, the GdnHCl reacts with Hsp, so we control the concentration of the GdnHCl according to the amount of Hsp at the moment of introduction. We first try to introduce the same amount of GdnHCl as Hsp. Then we increase the amount of GdnHcl to research its relationship with the final GFP level.
The time of introduction is set to 60min at the very beginning. Similarly, we will modify it to research its relationship with the final GFP level.
3.Heat shock
We first set the heat shock temperature to 42℃, and the time is set to 10min after the reaction begins. We will also modify the time and the temperature in order to study their relationship with the final GFP level.
According to Guisbert et al.(1), when the temperature rises from 30℃ to 42℃, the expression rate increases about 10 folds, but the degradation rate also increases about 3 folds.
According to Arrhenius formula, the rate constant of a reaction is k=Aexp(-E/RT). From the information in Guisbert et al.(1), we calculate the value of E/R of expression and degradation of Hsp, which are 18314K and 8738K respectively. So
\[K1(T)=k_{f}\frac{exp(-\frac{18314}{T})}{exp(-\frac{18314}{273+30})}\] \[K2(T)=d_{m}\frac{exp(-\frac{8738}{T})}{exp(-\frac{8738}{273+30})}\]
4.Result
(1)The variation trend of GFP
We simulate the level of GFP in a scope of 1000 minutes, which is shown in Figure 3 below.
Every parameter uses the value mentioned in previous parts.
Figure 3
(2)The influence of concentration of GdnHCl
As mentioned above, when introducing GdnHCl, the concentration is set according to the amount of Hsp in the yeast. We vary the concentration of GdnHCl from 1 to 1000 fold of the amount of Hsp.
For every amount of GdnHCl, there is a final GFP level. We normalize the GdnHCl concentration to the amount of Hsp(express the GdnHCl concentration by the fold), then we plot the final GFP to the fold, which is shown in Figure 4 below.
Figure 4
The slope is -0.956, very approximate to -1, which means the final GFP level is about inversely proportional to the concentration of GdnHCl.
We also plot two reactions, the concentration of GdnHCl in one of them is twice of the other(Figure 5).
Figure 5
(3)The influence of time point of introducing GdnHCl
We vary the time point of introducing of GdnHCl. The result is that the final GFP level is not relevant to the time point of introducing GdnHCl(Figure 6).
Figure 6
But the trend of GFP value may be different when we take a different introducing time point. For example, if we introduce GdnHCl as early as 20min, the GFP level increases monotonously and reach the equilibrium state(Figure 7).
Figure 7
(4)The influence of heat shock temperature
The influence of heat shock temperature on final GFP level is plotted below(Figure 8). We also fit it, and it is a exponential curve. However, it seems that the influence is not very significant.
Figure 8
(5)The influence of heat shock time point
We research the relationship between the final GFP level and the time point of heat shock.
It is a linear relationship, with a very small slope. So the time point when the heat shock takes place hardly affects the final GFP level(Figure 9 and Figure 10).
Figure 9
Figure 10
We also tried to heat shock the yeast after introduction of GdnHCl(300 min after the reaction begins).The trend is quite different, but the final GFP level does not change much(Figure 11).
Figure 11
From this figure, we can find that there is an abrupt increase at 300min, which suggests that heat shock can increase the GFP level(or more essentially, the amount of PR aggregation) in a short time.
5.Conclusion
We can get the conclusion that the final level of GFP(or other product proteins) is controlled by heat shock and introduction of GdnHCl. However, the time point and even the order of these operations will not affect the final outcome. Only the volumes(or temperature) count. While the concentration of GdnHCl plays a major role in controlling the final outcome, the heat shock temperature plays an auxiliary role. Combining these two operations, we can control the final outcome of product proteins in both large and small order of volume.
The source code can be found here:Matlab Code 1 Matlab Code 2
Gene circuit 2: sfGFP1-10 and GFP-11
1.Reactions
The second gene circuit we use is shown in Figure 12 below.
Figure 12
The principle of this gene circuit is similar to gene circuit 1, except for lacking a transcriptional activator. The relationship between GFP level and aggregation is contrary to gene circuit 1, that is to say, the more GdnHCl introduced, the higher GFP level will come out.
The reactions related to this gene circuit are:
\[mRNA(sfGFP1-10) \rightarrow sfGFP1-10\] \[mRNA(PRsfGFP11) \rightarrow PRsfGFP11\] \[sfGFP1-10+PRsfGFP11 \rightarrow CPX1\] \[PRsfGFP11+E(pre) \leftrightarrow CPX3\rightarrow aggregates+E(pre)\] \[mRNA(Epre)\rightarrow E(pre)\] \[E(pre)+GdnHCl\rightarrow E+GdnHCl\] \[E+aggregates\leftrightarrow CPX2\rightarrow PRsfGFP11+E\]
And the differential equations are:
\begin{align*} \frac{d[mRNA(sfGFP1-10)]}{dt}=&N*k_{mRNA}*[mRNA(sfGFP1-10)]\\\\ \frac{d[mRNA(PRsfGFP11)]}{dt}=&N*k_{mRNA}*[mRNA(PRsfGFP11)]\\\\ \frac{d[sfGFP1-10]}{dt}=&k[mRNA(sfGFP1-10)]-d[GFP-10]\\&-k_{f}[PRsfGFP11]*[sfGFP1-10]\\\\ \frac{d[PRsfGFP11]}{dt}=&k*[mRNA(PRsfGFP11)]-d*[PRsfGFP11]\\&-k_{f}[PRsfGFP11]* [sfGFP1-10]-k_{f}*M*[PRsfGFP11]\\&*[E(pre)]+k_{b}*M*[CPX3]+k_{f}*M1*[CPX2]\\\\ \frac{d[aggregate]}{dt}=&k_{agg}*[PRsfGFP11]-d*[aggregate]+k_{f}*M1* [CPX2]\\&-k_{f}*M1*[aggregate]*[E]\\\\ \frac{d[CPX1]}{dt}=&k_{f}*[PRsfGFP11]*[sfGFP1-10]-d*[complex]\\\\ \frac{d[CPX2]}{dt}=&-d[CPX2]+M1*k_{f}[E]*[aggregate]-M1*(k_{f}+k_{b})* [CPX2]\\\\ \frac{d[CPX3]}{dt}=&-(d+M*(k_{f}+k_{b}))[CPX3]+M*k_{f}*[E(pre)]*[PRsfGFP11]\\\\ \frac{d[mRNA(Epre)]}{dt}=&N*K1(T)-K2(T)[mRNA(GFP)]\\\\ \frac{E(pre)}{dt}=&k[mRNA(Epre)]-d[E(pre)]-M1*k_{f}[E(pre)]*[GdnHCl]\\ &-M*k_[f][E(pre)]*[PRsfGFP11]+[PRsfGFP11]+M*(k_{f}+k_{b})*[CPX3]\\\\ \frac{d[E]}{dt}=&-d[E]+M1*k_{f}[E(pre)]*[GdnHCl]-M*k_{f}[E]*[aggregate]\\& +M1*(k_{f}+k_{b})*[CPX2]\\\\ \frac{d[GdnHCl]}{dt}=&-d[GdnHCl]\\\\ \end{align*}
Some functions and parameters, such as K1(T), K2(T), M, M1, are the same as those mentioned before.Introduction of GdnHCl and heat shock are discussed before, and are omitted here.
2.Result
(1)The variation trend of GFP
We simulate the level of GFP in a scope of 2000 minutes, which is shown in Figure 13 below.
Every parameter uses the value mentioned in previous parts.
Figure 13
(2)The influence of concentration of GdnHCl
Result of this part is not the same as gene circuit 1. In this gene circuit, the GFP level is only sensitive to concentration of GdnHCl in a particular range(Figure 14).And the sensitivity is not as obvious as in gene circuit 1.
Figure14
The concentration of GdnHCl is normalized to the concentration of Hsp at the moment of introduction.
(3)The influence of time point of introducing GdnHCl
We vary the time point of introducing of GdnHCl. And the concentration of GdnHCl is set to be 1% of the Hsp. This is because according to Figure 14, at this concentration, the GFP level is sensitive to changes and not saturated. The influence of the time point when introducing GdnHCl is not significant, as shown in Figure 15.
Figure 15
Observing the figure that the GdnHCl is added at 1000min, we can see that GdnHCl can motivate the expression of GFP(Figure 16).
Figure 16
To sufficiently show the influence of GdnHCl, we introduce it at 1000min in subsequent modeling.
(4)The influence of heat shock temperature
The influence of heat shock temperature on final GFP level is plotted below(Figure 17).
It seems to be more obvious than that in gene circuit 1. It is an linear relationship. Heat shock activates the aggregation of PR and thus represses the GFP level.
Figure 17
(5)The influence of heat shock time point
This influence is not very important, as shown in Figure 18.
Figure 18
We show the result of heat shock at 500min and 1500min in Figure 19 and 20, respectively.
Figure 19
Figure 20
From these two figures, we can conclude that heat shock activates the aggregation of PR and thus represses the GFP level.
5.Conclusion
We can get the conclusion that the final level of GFP(or other product proteins) is controlled by heat shock and introduction of GdnHCl. However, the time point and even the order of these operations will not affect the final outcome. Only the volumes(or temperature) count. The volume of GdnHCl only affects the final output in a range(about 2 orders), and saturates both above and below this range. The heat shock temperature also affects the final GFP level, but not very significant(only about 1.5 folds less from 30℃ to 42℃).
The source code can be found here:Matlab Code 3 Matlab Code 4
Reference 1.Eric Guisbert, Takashi Yura, Virgil A. Rhodius and Carol A. Gross. 2008. Convergence of Molecular, Modeling, and Systems Approaches for an Understanding of the Escherichia coli Heat Shock Response. MICROBIOL. MOL. BIOL. REV. VOL. 72, 545-554.
Introduction
A PR-sfGFP11 complex can bind with another same complex, repressing binding between sfGFP11 and sfGFP1-10. This indicates that we can control the final GFP level by modulating the binding reaction among PR-sfGFP11 complexes. We must take multistep binding into consideration. That is to say that three or even more PR-sfGFP11 complexes can bind together and form large aggregation. This factor increases the complexity of this system and thus our modeling.
Calculation
1.Aggregation of PRs
In fact, the binding mechanism of PRs are not very well understood. We only know that upon heat shock, PRs will aggregate. And there are some inhibitors of this process, such as GdnHCl. In this part, we use detailed calculation to prove that the GFP level is sure to decrease after aggregation of PRs.
We are not aware that how many PRs form an aggregation. And there is a large possibility that this number is not definite. We assume that the upper limit of the number of PRs in an aggregation is N. Then, we denote the mass of sfGFP1-10, sfGFP11-PR complex (including a linking peptide between sfGFP11 and PR) as M1 and M2, respectively. We also denote the number of aggregations containing i complexes as N(i). Each of them has a possibility of p1(i) to bind a complex again, and a possibility of p2(i) to drop a complex, in a unit of time. The binding and dissociation of one single complex and one single aggregation, are regarded as Poisson processes, which will not happen more than once in a time unit. And the possibilities are independent with time. These assumptions are shown in the diagrams below (Figure 1).
Figure 1
In the final equilibrium state, there must be N(i)p2(i)+N(i)p1(i)=N(i-1)p1(i-1)+N(i+1)p2(i+1) The formula above is not suitable for i=1. When i=1, it should be replaced by N(1)p1(1)=N(2)p2(2). It is noteworthy that if N(i-1)p1(i-1)= N(i)p2(i), there will be N(i)p1(i)= N(i+1)p2(i+1), i>1. Since we already have N(1)p1(1)=N(2)p2(2), we can use the method of mathematical induction to draw the conclusion that for every positive integer i, N(i)p1(i)= N(i+1)p2(i+1).
Binding of sfGFP1-10 and sfGFP11 is a diffusion-controlled reaction, whose velocity is directly determined by the speed of diffusion of the reactants in the solutions. According to statistical physics, the average velocity of a molecule is inversely proportional to the square root of its mass. And it is also easy to prove, the relative velocity of two molecules is inversely proportional to the square root of their reduced mass. The frequency of collision between two molecules are proportional to the product, vS. v is the relative velocity, and the S denotes the effective collisional section. The relative velocity of a single complex and an aggregates containing i complexes is proportional to (i/(i+1))^(-1/2). And the effective collisional section, is proportional to (i^(1/3)+1)^2, because the scale is proportional to i^(1/3). As for dissociation, the more complexes an aggregation has, the larger likelihood that a single complex drops off there will be. So p2(i) is proportional to i. So we have
\[p1(i)\propto(\frac{i+1}{i})^\frac{1}{2} (i^\frac{1}{3}+1)^2\] \[p2(i)\propto i\]
p1(1) and p2(2) are just possibility of binding and dissociation of two complexes. They should be constants. We denote them as K1 and K2. So\[p1(i)=K_{1}(\frac{i+1}{i})^\frac{1}{2} (i^\frac{1}{3}+1)^2\] \[p2(i)=K_{2}i/2\]
Recalling the relationship proved above, N(i)p1(i)= N(i+1)p2(i+1), we can derive\[\frac{N(i+1)}{N(i)}=\frac{p1(i)}{p2(i+1)}=\frac{2K_{1}}{K_{2}} (\frac{1}{i(i+1)})^\frac{1}{2}(i^\frac{1}{3}+1)^2\]
In steady state, the total amount of sfGFP11-PR complexes is constant, that is\[\sum_{i=1}^N iN(i)=N_{0}\]
Based on the two equations above, we can calculate the proportion of every kind of aggregates.2.Binding of sfGFP1-10 and sfGFP11
In order to confine sfGFP11 in the aggregation and prevent the reaction between sfGFP11 and sfGFP1-10, there must be at least 4 complexes in the aggregation forming a structure of regular tetrahedron. In this case, sfGFP11 will be restricted in the ‘central hole’ of the PR tetrahedron, since sfGFP11 is much smaller than PR (as shown in Figure 2 below).
Figure 2
For aggregates containing fewer than 4 complexes, sfGFP11 molecules are sure to be bare and can react with sfGFP1-10. However, the aggregates lower down the diffusional velocity, thus reduce the level of GFP. As mentioned before, the relative velocity of two molecules is inversely proportional to the square root of their reduced mass. So we have
\begin{align*} [GFP]'(i)=&\sqrt{\frac{\frac{M_{1}M_{2}}{M_{1}+M_{2}}}{\frac{M_{1}*iM_{2}}{M_{1}+i*M_{2}}}} [GFP]\\=&\sqrt{\frac{M_{1}+i*M_{2}}{i(M_{1}+M_{2})}}[GFP] &(i=1,2,3.) \end{align*}
Where [GFP] is the level of GFP in the absence of aggregation. Combing this two mechanism, the total GFP level is\[[GFP]'=\frac{\sum_{i=1}^3 iN(i)\sqrt{\frac{M_{1}+i*M_{2}}{i(M_{1}+M_{2})}}}{\sum_{i=1}^N iN(i)}[GFP]\]
3.Quantitative result
We set the constants above to some values, and to see the quantitative repressive effect of PR aggregation on GFP level. We first assume that the mass of sfGFP1-10 (200-300Aa) is about identical to that of PR (200~300Aa), and 5 folds larger than that of sfGFP11. So M2=1.2*M1.
We first assume that the binding and dissociation of two sfGFP11-PR complexes are comparable, namely K1=K2, and we vary N from 1 to 10. The result is shown in Figure 3 below. We can see an abrupt drop at N=4, which means big aggregation that confine sfGFP11 can greatly reduce GFP level, while the repression of diffusion is only a subordinate factor.
Figure 3
Then we fix N=10. And we adjust K1 to change the dominance of binding or dissociation. The result is shown in Figure 4.
Figure 4
Figure 4 shows that even if when N is large, if binding of sfGFP11-PR complexes is not dominant, the GFP level can still be high. The situation is that when switched off (introducing GdnHCl), dissociation of complexes is dominant and yield a high GFP level. When switched on (heat shock), the binding the complexes is dominant, leading to a low GFP level. This allows us to switch the GFP between two widely apart levels and detect the difference easily.
MATLAB code is available at: Matlab Code 6
Discussion
We made two hypothesis in our calculation above.
First, we hypothesize that for large aggregation (i>=4), the sfGFP11s are totally restricted. This is not ensured, but it does not affect the final conclusion. At least this is the case in a part of large aggregation, and thus can also cause significance decrease of GFP level.
Second, we neglect the fusion of two aggregation. Since one single complex can bind with an aggregation, two aggregation can also ‘bind’. We think this is not very likely to happen. Aggregation is in a steady structure. At this time, introduction of another aggregation rather than a single complex, will destroy the steady structure, which requires a lot of energy and is almost impossible.
Conclusion
Taking advantage of the mechanism that sfGFP11 will be constrained in PR aggregates, we can control the GFP level by switching on or switching off the aggregation of PRs.
Introduction
In previous part, we discussed binding of GFP splits. A PR-sfGFP11 complex can bind with another same complex, repressing binding between sfGFP11 and sfGFP1-10. The mechanism is also right for binding of AD and BD, which must be assisted by aggregation of PRs. However, it is much more complicated than the model for GFP splits before. For GFP splits, the PR aggregates consist only one kind of monomers, PR-sfGFP11. On the contrary, there are two different kinds of monomers in the PR aggregates, PR-AD and PR-BD.
In fact, the problem has evolved into a non-linear problem, we will avoid non-linear equations and develop a linear method to give a solution.
Calculation
1.Reaction
First, we list all the parameters and functions we used here: K11: possibility of reaction between PR-AD and PR-AD in a time unit. K12: possibility for a PR-AD to bind a PR-BD in a time unit. K2: possibility for a dimer to drop a PR-AD or PR-BD complex in a time unit. M1: mass of PR-AD. M2: mass of PR-BD. M3: mass of PR. n: the largest number of PR in aggregation. N(i,j): the relative amount of aggregation containing i PRs and j ADs. p11(i,j): the possibility for aggregation containing i PRs and j ADs to bind an AD in a time unit. (p11(1,1)=K11) p12(i,j): the possibility for aggregation containing i PRs and j ADs to bind an BD in a time unit. (p12(1,1)=K12) p2(i): the possibility for aggregation containing i PRs to drop a PR-AD or PR-BD complex in a time unit. (p2(2)=K2) row(i,j): the position of N(i,j) in the matrix we used in calculation.
We show all the possible reactions, together with some functions listed above, in Figure 5 below. The colors on each aggregate are not nonsense. They characterize the type of reactions and different methods we take in the calculation. Note that in the figure, there are some arrows marked 'p2(i)j/i' or some others in this form. This means that in an aggregate containing i PRs and j ADs, the possibility of dropping a PR-AD complex is 'p2(i)j/i', and the possibility of dropping a PR-BD complex is 'p2(i)(i-j)/i'. In other words, it is to say that the possibility for a PR-AD or a PR-BD to drop from an aggregate is the same.
Figure 5
We have discussed that the possibility of reaction is proportional to the relative diffusional velocity and the effective section, in the previous part about GFP. And we also mentioned that the relative diffusional velocity is inversely proportional to the square root of reduced mass.
Based on these points, we derive the general formula for p11(i,j), p12(i,j) and p2(i).
\begin{align*} p_{11}(i,j)=&K_{11}\left(\frac{\frac{M_{1}(M_{1}j+M_{2}(i-j))}{M_{1}(j+1)+M_{2}(i-j)}}{M_{1}/2}\right)^{-\frac{1}{2}} \left(\frac{M_{1}^\frac{1}{3}+(M_{1}j+M_{2}(i-j))^\frac{1}{3}}{2M_{1}^\frac{1}{3}}\right)^2\\ =&K_{11}\left(\frac{2(M_{1}j+M_{2}(i-j))}{M_{1}(j+1)+M_{2}(i-j)}\right)^{-\frac{1}{2}}\left(\frac{M_{1}^\frac{1}{3}+(M_{1}j+M_{2}(i-j))^\frac{1}{3}}{2M_{1}^\frac{1}{3}}\right)^2\\\\ p_{12}(i,j)=&K_{12}\left(\frac{\frac{M_{2}(M_{1}j+M_{2}(i-j))}{M_{1}j+M_{2}(i-j+1)}}{\frac{M_{1}M_{2}}{M_{1}+M_{2}}}\right)^{-\frac{1}{2}} \left(\frac{M_{2}^\frac{1}{3}+(M_{1}j+M_{2}(i-j))^\frac{1}{3}}{M_{1}^\frac{1}{3}+M_{2}^\frac{1}{3}}\right)^2\\ =&K_{12}\left(\frac{(M_{1}+M_{2})(M_{1}j+M_{2}(i-j))}{M_{1}(M_{1}j+M_{2}(i-j+1))}\right)^{-\frac{1}{2}} \left(\frac{M_{2}^\frac{1}{3}+(M_{1}j+M_{2}(i-j))^\frac{1}{3}}{M_{1}^\frac{1}{3}+M_{2}^\frac{1}{3}}\right)^2\\\\ \end{align*}
2.Set up equations
Till now, we have make every reaction in this system clear. The task now is to calculate relative amount of each type of aggregation (i.e. N(i,j)). We use the requirements of steady state to derive a series equations. Unlike the PR-sfGFP11 aggregation discussed before, the 'reaction web' forms a triangle here, which means whichever angle we start from, we will encounter greater and greater number of unknown variables, and the number of equations are permanently insufficient. For example, if we are given the value of N(n,0), using steady-state condition for N(n,0), we can calculate N(n-1,0), but fail to calculate N(n,1). Even if we are further given N(n,1), we cannot calculate N(n-2,0), N(n-1,1), N(n,2), as shown in Figure 6.
Figure 6
It is also not feasible to start from one whole side. This is because when move inward from a side, the number of equations are always larger than variables. And the equation sets are always contradictory which derive no solution (Figure 7).
Figure 7
The only way to calculate all the N(i,j) is to calculate them as a whole. There all totally (n+3)n/2 N(i,j). And there are the same number of conditions of steady state (one condition for one N(i,j)). However, there are one repeated condition. In fact, for example, if we know that all N(i,j) except for N(n,0) are already in steady state. If N(n,0) changes, the N(n-1,0), which is relevant to N(n,0) through a reaction, is sure to change. So when ((n+3)n/2-1) N(i,j) are steady, the rest one is also steady. So there are only ((n+3)n/2-1) valid equations. But it is enough for us to complete the calculation. What we only want is relative amount of each N(i,j). So we can simply set N(1,0) to 1, and normalize other N(i,j). Then we only need to determine ((n+3)*n/2-1) unknown variables, whose amount is identical to the number of equations and can be surely solved out.
We list all the steady-state conditions for all N(i,j) below:
At the end of every equation, there is a color, which is identical to the colors in Figure 5. Colors here and those in Figure 5 have the same meaning. We use color here to help distinguish various conditions.
It is necessary to calculate using matrix. We put all the unknown N(i,j) into a column vector in the following order: N(1,1), N(2,0), N(2,1), N(2,2),N(3,0),...,N(3,3),...,N(n,0),...,N(n,n). We use a function row(i,j) to determine the position of N(i,j) in the vector. Here is its pseudo-code:
*if i==1
y=1;
else
y=(i+1)i/2+j-1;
end*
Then we need to determine the coefficient matrix a(p,q) and non-linear terms b(p). This is based on the equations listed above. Here we give the pseudo-code.
a(row(n,n),row(n,n))=0;
b(row(n,n))=0;
if n>2
for i=1:n
for j=0:i
if i==1&&j==1 %light green
a(row(i,j),row(1,1))=p12(1,1+1)+p11(1,1+1);
a(row(i,j),row(i+1,j))=-p2(i+1)*(i+1-j)/(i+1);
a(row(i,j),row(i+1,j+1))=-p2(i+1)*(j+1)/(i+1);
end
if i==n&&j==0 %purple
a(row(i,j),row(n,0))=p2(n);
a(row(i,j),row(n-1,0))=-p12(n-1,0+1);
end
if i==n&&j==n %red
a(row(i,j),row(n,n))=p2(n);
a(row(i,j),row(n-1,n-1))=-p11(n-1,(n-1)+1);
end
if i==n&&j~=0&&j~=n %light blue
a(row(i,j),row(n,j))=p2(n);
a(row(i,j),row(n-1,j))=-p12(n-1,j+1);
a(row(i,j),row(n-1,j-1))=-p11(n-1,(j-1)+1);
end
if i==j&&i~=1&&i~=n %dark blue
a(row(i,j),row(i,i))=p2(i)+p11(i,i+1)+p12(i,i+1);
a(row(i,j),row(i+1,i))=-p2(i+1)/(i+1);
a(row(i,j),row(i+1,i+1))=-p2(i+1);
a(row(i,j),row(i-1,i-1))=-p11(i-1,(i-1)+1);
end
if i~=1&&i~=n&&j==0 %orange
a(row(i,j),row(i,0))=p11(i,0+1)+p12(i,0+1)+p2(i);
a(row(i,j),row(i+1,1))=-p2(i+1)/(i+1);
a(row(i,j),row(i+1,0))=-p2(i+1);
if i==2
b(row(i,j))=p12(i-1,0+1);
else
a(row(i,j),row(i-1,0))=-p12(i-1,0+1);
end
end
if i~=1&&i~=n&&j~=0&&j~=i %white
a(row(i,j),row(i,j))=p11(i,j+1)+p12(i,j+1)+p2(i);
if i==2&&j==1
b(row(i,j))=p11(i-1,(j-1)+1);
else
a(row(i,j),row(i-1,j-1))=-p11(i-1,(j-1)+1);
end
a(row(i,j),row(i+1,j+1))=-p2(i+1)*(j+1)/(i+1);
a(row(i,j),row(i+1,j))=-p2(i+1)*(i+1-j)/(i+1);
a(row(i,j),row(i-1,j))=-p12(i-1,j+1);
end
end
end
else
if n==2
a(1,1)=p12(1,1+1)+p11(1,1+1);
a(1,3)=-p2(2)/2;
a(1,4)=-p2(2);
a(2,2)=p2(2);
a(3,1)=-p12(1,1+1);
a(3,3)=p2(2);
a(4,1)=-p11(1,1+1);
a(4,4)=p2(2);
b(2)=p12(1,0+1);
b(3)=p11(1,0+1);
end
Then the matrix equation is
\begin{align*} (a_{pq})_{row(n,n)\times row(n,n)}(N(1,1),\cdots,N(i,0),\cdots,N(i,j),\cdots, \\ N(i,i),\cdots,N(n,0),\cdots,N(n,j),\cdots,N(n,n))^T\\\\ =(b_{p})_{row(n,n)\times 1} \end{align*}
For example, the matrix equation for n=2 is:\[ \left( \begin{array}{cccc} p_{12}(1,1)+p_{11}(1,1) & 0 & -\frac{p_{2}}{2} & -p_{2}(2) \\ 0 & p_{2}(2) & 0 & 0 \\ -p_{12}(1,1) & 0 & p_{2}(2) & 0 \\ -p_{11}(1,1) & 0 & 0 & p_{2}(2) \end{array} \right) \left( \begin{array}{c} N(1,1)\\ N(2,0)\\ N(2,1)\\ N(2,2) \end{array} \right) = \left( \begin{array}{c} 0\\ p_{12}(1,1)\\ p_{11}(1,1)\\ 0\\ \end{array} \right) \]
Then we can simply calculate all the N(i,j) out.
It is time to recall what we mentioned in the beginning of this chapter. We mentioned that this problem is non-linear. All we have completed is only the linear method. But to use this method, we must do some preparation to transfer a non-linear problem to a linear one.
It is noteworthy that 'a PR-AD binds a PR-BD' and 'a PR-BD binds a PR-AD' is the same reaction! Due to the uniqueness of this reaction, there should be p11(1,0)/p12(1,1)=N(1,1)/N(1,0)=N(1,1). In other words, K1 and K2 are cannot directly set by ourselves, they have relationship with N(1,1) and N(1,0), which is actually intuitive. In fact, we can only set K1, and then express K2 using N(1,1) and N(1,0). However, that will lead to a non-linear equation set.
In order to avoid this problem, we set K1 first, then we change K2 in a proper interval and calculate p11(1,0)/p12(1,1)-N(1,1)/N(1,0). When the result is 0 (or smaller than a limit of accuracy), we think we have find the actually K2. And then, we use K1 and K2 to finish the linear calculation and get the result of N(i,j). This is the algorithm we actually used. Though it is troublesome and time-consuming for computer, it at least make this problem solvable.
3.Data Analysis
We define a function called 'efficiency' here. It means the ratio of effective AD-PR-BD complex divided by the largest possible amount of effective AD-PR-BD.
The most efficient way to generate AD-PR-BD complex is to let every PR-AD only bind with one PR-BD. So the largest possible amount of effective AD-PR-BD is a half of the amount of PRs. 'Effective AD-PR-BD' means the aggregation at least has an AD and a BD out of 'PR envelop', thus has the ability to bind and activate the promoters on another gene circuit.
Now we discuss when all ADs or BDs will be enveloped by PRs. Unlike sfGFP11 discussed before, AD and BD have comparable size with PR. This means we must perform a more detailed calculation to determine the number of PRs needed. We postulate that N PRs form a sphere whose radius is R. And we denote the radius of a single PR is r. Then, we have
\[N\pi r^2=4\pi R^2,\mbox{and then} R=\frac{\sqrt{N}r}{2}\]
The maximal efficiency of occupying space in a pile of stacking sphere is 74.05% (face centered cubic or hexagonal close packed). So the valid space is\begin{align*} V=&0.7405\times\frac{4}{3}*\pi R^3\\=&0.7405\times\frac{4}{3}*\frac{N^\frac{3}{2}}{8}\\\\ \end{align*}
Take AD as an example. Denote the radius of AD as r1. When all the AD are just covered, there should be\begin{align*} V=&N*\frac{4}{3}*\pi r_{1}^3\\=&N*\frac{4}{3}*\pi r_{1}^3\frac{M_{1}-M_{3}}{M_{3}}\\\\ \end{align*}
Then we have\[N_{1}=\frac{64(M_{1}-M_{3})^2}{0.7405^2M_{3}^2}\]
If it is for BD, we just change M1 into M2. For the smaller one of M1 and M2, the corresponding N is also smaller. We think the smaller molecule are more prone to be covered, which means once the amount of PR exceed the smaller one of N1 and N2, the aggregation loses its function.Then we can give the formula of the efficiency function:
\[\mbox{efficiency}=\frac{2\sum_{i=2}^{min(N_{1},N_{2})}\sum_{j=i}^(i-1)} {\sum_{i=1}^N\sum_{j=0}^i i*N(i,j)}\]
We discuss the efficiency in three different conditions in following part.
1) The forming of aggregation is dominant
We set K11=10, K2=1. Then we determine K12 and calculate efficiency. We only increase I from 2 to 8. The result is shown below in Figure 8 and Figure 9.
Figure 8
Figure 9
We can conclude that efficiency decreases and K12 increases as n increases.
2)The binding and dropping of PR-AD or PR-BD complex are comparable.
We first set K11=K12=K2=1, when i ranges from 2 to 100, the maximal absolute value of p11(1,0)/p12(1,1)-N(1,1)/N(1,0) is only 0.0324. So we think K12 just equals to 1. Then we plot efficiency to n, the result is shown in Figure 10 below.
Figure 10
We find that in this case, the efficiency is even higher than the condition that binding reaction is dominant. This may be due to two reasons: (a)The amount of PR-AD and PR-BD are almost the same, thus the dimerization is facilitated; (b)Because binding reaction is not dominant, the proportion of big aggregation are relatively small, thus enhancing the efficiency.
We also calculated N1=27.5354, N2=46.1147. So min(N1,N2)=27.5354. However, we observe no abrupt change in efficiency at n=28. On the contrary, the efficiency seems to remain on a level. This may be because big aggregation (for example, i=28) already lead to a low efficiency. The emerging of 0-efficiency aggregation will not influence the total efficiency much. What's more, we calculated the amount of aggregation of different sizes when n=40 (Figure 11).The result shows that there are nearly no aggregation containing more than 20 complexes. And this is another reason for which there is no abrupt change at n=27.
Figure 11
3)The decomposition is dominant
We first set K11=K12=1, K2=10. When i ranges from 2 to 100, the maximal absolute value of p11(1,0)/p12(1,1)-N(1,1)/N(1,0) is 0.0031. So we think the K12 is just 1.
The we plot the efficiency with regard to n (Figure 12). The efficiency is much lower than before.
Figure 12
Theoretically, we have proven the switch between binding-dominant state and decomposition-dominant state can obviously affect the downstream promoters. However, since when K2=10,K11=K12=1, efficiency is about 0.1 and about a quarter of the efficiency when K11=10 and K2=1. This extent of reduction may be not satisfying. But if the switch has a weaker leaky binding rate of PRs, the efficiency will continue to decrease. For example, if K2=100,K11=K12=1, for all i from 2 to 100, the efficiency is always around 0.0133, which is much lower than those when forming of aggregation is dominant.
Conclusion
We have developed a valid method to describe the aggregation of PR-AD and PR-BD. The result shows that the switch is effective and can actually have effect on downstream promoters.
This model and computing method, can also promoted to other binding and decomposing reactions mediated by PR.
MATLAB code for this part is available at: Matlab Code 7 Matlab Code 8 Matlab Code 9 Matlab Code 10
We designed a water bath kettle as shown in Figure 1. The main body of the kettle is the same as normal water bath kettles. However, there is a long and thin tube, extending out of the kettle from 1/3 of the total height. There are 3 sensors on the tube, which measure the local temperature respectively. To adjust the temperature automatically and quickly, we also designed a tube for cold water at the top of the kettle. We denote the length of sides of the bottom as a (The bottom is a square), and the total height as H, the height of the water as h.
Figure 1
For simplicity, we think the diameter of the tube is too small to affect the temperature field in the main body of the kettle. Similarly, we also neglect the effect on sensors on the tube. The tube and the sensors are all covered with heat-insulated material. In our model, we just consider this material ideally and consider no leakage of heat through it.
We first give the solution to the steady state of this system, and then we discuss the common time-varying solution. Based on these discussions, we give the formula of how the three sensors indicate the temperature in the kettle. In the kettle, what we would like to consider the most is the surface of the water, where the test tubes we heat float.
Main Body of the Kettle
We denote the temperature filed in the water as φ1, and that in the air inside the kettle as φ2. The coefficient of heat conducting through the crust of the kettle is denoted as α. The heat that the heater provides on a unit area of the bottom in a time unit is denoted as Q. And the coefficients of heat conducting inside uniform water and uniform air are denoted as k1 and k2 respectively.
Then we can derive the functions of temperature field: In the water,
\[\frac{\partial{\phi_{1}}}{\partial{t}}=a_{1}^2 \Delta \phi_{1}\] \[k_{1}\frac{\partial{\phi_{1}}}{\partial{n}}+\alpha\phi_{1}=Q.\quad (z=0) \]
Above the water,\[\frac{\partial{\phi_{2}}}{\partial{t}}=\alpha_{2}^2\Delta\phi_{2}\] \[k_{2}\frac{\partial{\phi_{2}}}{\partial{n}}+\alpha\phi_{2}=0.\quad (z=H) \] \[\phi_{1}=\phi{2}.\quad (z=h)\] \[k_{1}\frac{\partial{\phi_{1}}}{\partial{z}}= k_{2}\frac{\partial{\phi_{2}}}{\partial{z}}.\quad(z=h)\]
Here a^2=k/(cρ). c is the specific heat capacity, and ρ is the density of the matter. x and y are the dimensions of the two perpendicular sides of the bottom, and z is the vertical direction. In addition, we cover the periphery of the kettle with heat-insulating materials to achieve better results of thermal control. Then we have
\[\frac{\partial{\phi_{1}}}{\partial{n}}=0, \qquad \frac{\partial{\phi_{2}}}{\partial{n}}=0. \quad (0
\[\Delta\phi_{1}=0,\quad\Delta\phi_{2}=0\]
\[k_{1}\frac{\partial{\phi_{1}}}{\partial{n}}+\alpha\phi_{1}=Q.\quad(z=0)\]
\[k_{1}\frac{\partial{\phi_{1}}}{\partial{z}}= k_{2}\frac{\partial{\phi_{2}}}{\partial{z}}.\quad(z=h)\]
\[\phi_{1}=\phi{2}.\quad (z=h)\]
\[\phi_{1}=Az+B,\]
\[\phi_{2}=Cz+D,\]
\[\alpha B-k_{1}A-Q=0,\]
\[\alpha(CH+D)+k_{2}C=0\]
\[Ah+B=Ch+D,\]
\[k_{1}A=k_{2}C.\]
\[A=-\frac{k_{2}Q}{k_{2}(\alpha h+2k_{1})+k_{1}\alpha (H-h)}\]
\[B=\frac{k_{1}A+Q}{\alpha}\]
\[C=-\frac{k_{1}Q}{k_{2}(\alpha h+2k_{1})+k_{1}\alpha (H-h)}\]
\[D=-\frac{\alpha H+k_{2}}{\alpha} C\]
\[T_{h}=Ah+B=\frac{k_{1}Q(k_{2}+(H-h)\alpha)}{\alpha(k_{2}(\alpha h+2k_{1})+k_{1}\alpha(H-h))}\]
\[T_{0}=\frac{AH}{3}+B=\frac{Q}{3\alpha}\left(3-\frac{k_{2}(3k_{1}+H\alpha)}{k_{2}(\alpha h+2k_{1})+k_{1}\alpha(H-h)} \right)\]
We now discuss the temperature field in the tube. In the last step, we have already got the temperature T0 at the mouth. Due to the pressure at the mouth, the water in the tube has a velocity, which makes the temperature field different from that in stationary water. According to the equation of energy in fluid mechanics,
\[\frac{d\epsilon}{dt}=\frac{1}{\rho}\sigma_{ji}\frac{\partial{v_{i}}}{\partial{x_{j}}} -\frac{1}{\rho}\frac{\partial{q_{i}}}{\partial{x_{i}}}\]
The long tube
According to Poiseuille's formula, in the tube, the velocity of fluid is parallel with the axis, and it has a relationship with r(radial position):
\[v(r)=\frac{G}{4\mu}(r_{0}^2-r^2)\]
where r_0 is the radius of the tube, G is the gradient of pressure, and μ is the coefficient of viscosity. We use polar coordinates(r, θ, z) to describe the position in the tube. Since the velocity is in z direction and independent on z, we have\begin{align*} \sigma_{ji}\frac{\partial{v_{i}}}{\partial{x_{j}}}=& \sigma_{j3}\frac{\partial{v_{3}}}{\partial{x_{j}}}\\=& \sigma_{13}\frac{\partial{v_{3}}}{\partial{x_{1}}}+ \sigma_{23}\frac{\partial{v_{3}}}{\partial{x_{2}}}\\=&\mu \left( (\frac{\partial{v_{3}}}{\partial{x_{1}}})^2+ (\frac{\partial{v_{3}}}{\partial{x_{2}}})^2 \right)\\=&\mu \left( (\frac{\partial{v_{3}}}{\partial{r}})^2+\frac{1}{r^2} (\frac{\partial{v_{3}}}{\partial{\theta}})^2 \right)\\ =&\mu\left(\frac{\partial{v_{3}}}{\partial{r}} \right)^2\\ =&\frac{G^2}{4\mu}r^2\\\\ \end{align*}
In the fluid,\[q_{i}=-k_{1}\frac{\partial{T}}{\partial{x_{i}}}\]
So\[\frac{\partial{q_{i}}}{\partial{x_{i}}}=-k_{1}\Delta T\]
Substitute these results into the ‘equation of energy’, we get\[c\rho\frac{\partial{T}}{\partial{t}}=\frac{G^2}{4\mu}r^2+k_{1}\Delta T\]
In the steady state,\[\Delta T+\frac{G^2}{4\mu k_{1}}r^2=0\]
To solve this equation, we introduce\[T_{1}=T+\frac{G^2}{64\mu k_{1}}r^4\]
which transforms the equation to\[\Delta T_{1}=0\]
Let T1=R(r)Z(z), separate the variables, we have\[\frac{(rR')'}{rR}+\frac{Z''}{Z}=0\]
Because the tube is heat-insulated, R' (r0 )=0, where r0 is the radius of the tube. At the end of the tube, the water is directly exposed in the outside environment, whose temperature is set to 0. All sets of possible solutions are:\begin{equation*} \left\{ \begin{aligned} Z_{0}=&l-z\\ R_{0}=&1 \end{aligned} \right. \qquad \left\{ \begin{aligned} Z_{n}=&sinh\omega_{n}(l-z)\\ R_{n}=&J_{0}(\omega_{n}r) \end{aligned} \right. \end{equation*} \[\omega_{n}\mbox{is the nth root of}J'_{0}(\omega_{n}r)\]
So\[T_{1}=A_{0}(l-z)+\sum_{n=1}^{+\infty} A_{n}sinh\omega_{n}(l-z)J_{0}(\omega_{n}r)\]
Note that\[T_{1}(z=0)=A_{0}l+\sum_{n=1}^{+\infty}A_{n}sinh(\omega_{n}l) J_{0}(\omega_{n}r)=T_{0}+\frac{G^2}{64\mu k_{1}}r^4\]
Thus we can calculate all the An.\begin{align*} A_{0}=&\frac{\int_{0}^{r_{0}} r(T_{0}+\frac{G^2}{64\mu k_{1}}r^4)dr}{l\int_{0}^{r_{0}} rdr}\\ =&\frac{T_{0}+\frac{G^2}{192\mu k_{1}}r_{0}^4}{l}\\\\ A_{n}=&\frac{\int_{0}^{r_{0}}r(T_{0}+\frac{G^2}{64\mu k_{1}}r^4)J_{0}(\omega_{n}r)dr}{sinh(\omega_{n}l) \int_{0}^{r_{0}}rJ_{0}^2(\omega_{n}r)dr}\\ =&\frac{r_{0}J_{3}(\omega_{n}r_{0})(8-\omega_{n}^2 r_{0}^2)G^2}{32J_{0}^2(\omega_{n}r_{0})\omega_{n}^3 \mu k_{1}sinh(\omega_{n}l)}\\\\ \end{align*}
We give all the parameters above some proper values, and we plot the temperature in z direction and r direction (Figure 2 and Figure 3). We find that the temperature is almost independent with r, and simply change linearly with z, which indicates that the effect of velocity of fluid is negligible. The data we use is listed here, all the unit is omitted because they are in International System. l=0.3, k1=0.58, k2=0.023, Q=2e4, a=1, H=0.6, h=0.5, α=800, r0=2.5e-3 μ=8e-4.
Figure 2
Figure 3In Figure 3, the red line corresponds to the temperature at the mouth of the tube, and the blue line corresponds to the middle of the tube. The independence of temperature on the radial position is because of the extremely small radius of the tube, which allows heat to diffuse evenly. We also study the relationship of temperature at the water surface with Q and h, which is shown in Figure 4.
Figure 4We can conclude from Figure 4 that temperature difference is proportional to Q. And it decreases as h increases. So we should not fill the kettle with too much water.
Sensors
The sensors are all covered with heat-insulating materials, except the face connected to the tube. Actually, a sensor is just a short column filled with bacterial solution fixed around the tube (as shown in Figure 1). Because the small volume of the sensors and the heat-insulating material, we can know that the temperature inside the sensors is constant, which equals to the temperature in the tube at that point. This can also be rigorously proved.
It is noteworthy that we used three sensors and distribute them at different locations on the tube. This design aims to improve the accuracy of temperature measurement, which we will prove below.
The distances from sensors to the mouth of the tube are respectively denoted as x1, x2 and x3. The temperature they measure are respectively T1, T2 and T3. So the temperature at the mouth of the tube, T0, from which we can know the temperature at the water surface as long as h is given, can be expressed as
\[T_{0}=t_{1}=\frac{l}{l-x_{1}}T_{1}\] \[T_{0}=t_{2}=\frac{l}{l-x_{2}}T_{2}\] \[T_{0}=t_{3}=\frac{l}{l-x_{3}}T_{3}\]
Considering these three results from the three sensors at the same time, we have\begin{align*} T_{0}=&a_{1}t_{1}+a_{2}t_{2}+(1-a_{1}-a_{2})t_{3}\\ =&\frac{a_{1}l}{l-x_{1}}T_{1}+\frac{a_{2}l}{l-x_{2}}T_{2}+ \frac{(1-a_{1}-a_{2})l}{l-x_{3}}T_{3}\\\\ \end{align*}
where a1 and a2 are undetermined coefficient, we now determine them to make the final result T_0 as accurate as possible.We assume that the errors of T1, T2 and T3 are subjected to one same normal distribution, with the average value 0 and variance σ^2. Because of additivity of normal distribution, the average value of error of T0 is 0, and the variance of the error is
\[\sigma_{T}^2=l^2 \sigma^2\left(\frac{a_{1}^2}{(l-x_{1})^2}+ \frac{a_{2}^2}{(l-x_{2})^2}+\frac{(1-a_{1}-a_{2})^2}{(l-x_{3})^2} \right)\]
When σT2 reaches its minimal value, there should be\[\frac{\partial{\sigma_{T}^2}}{\partial{a_{1}}}=0, \qquad \frac{\partial{\sigma_{T}^2}}{\partial{a_{2}}}=0\]
Substitute σT2, we have\[\frac{a_{1}}{(l-x_{1})^2}=\frac{1-a_{1}-a_{2}}{(l-x_{3})^2} =\frac{a_{2}}{(l-x_{2})^2}\]
Solve this equation, we get
So
\begin{align*} T_{0}=&\frac{a_{1}l}{l-x_{1}}T_{1}+\frac{a_{2}l}{l-x_{2}}T_{2} +\frac{(l-a_{1}-a_{2})l}{l-x_{3}}T_{3}\\= &\frac{l}{\sum_{i=1}^3(l-x_{i})^2}((l-x_{1})T_{1}+(l-x_{2})T_{2} +(l-x_{3})T_{3})\\\\ \end{align*}
And the variance of error\begin{align*} \sigma_{T}^2=&l^2\sigma^2\left(\frac{a_{1}^2}{(l-x_{1})^2} +\frac{a_{2}^2}{(l-x_{2})^2}+\frac{(l-a_{1}-a_{2})^2}{(l-x_{3})^2}\right)\\ =&\frac{l^2\sigma^2}{\sum_{i=1}^3(l-x_{i})^2}\\\\ \end{align*}
If we set x1=0, which means place a sensor st the mouth, we can ensure
Thus we enhanced the accuracy of temperature measurement.
However, the mouth can have a high temperature, which is bad for the bacteria's survival. So we can place the sensors far from the mouth of the tube, for instance, at the middle position of the tube, which can control the variance of error less than 4 folds of that of a single sensor.
Automatic control
If we process the information of temperature from the three sensors, we calculate the temperature at the mouth. If we also know the height of water surface, we can calculate the temperature there, which we want most. What’s more, if we connect the processor with the valve of cold water at the top of the kettle and the heater at the bottom, we can automatically control the temperature in the kettle. The relationship of temperature with Q and h has been shown in Figure 4. We also plot the temperature of water surface to h in Figure 5, which shows the surface temperature will be obviously lower than that at the mouth of the tube. However, if we control the surface below 0.52m (total height is 0.6m), 80% temperature difference with the outside environment can be achieved, which is efficient enough. So we should also have a sensor that control h.
If we want to increase the temperature, we can only turn up the heater, and let the steady temperature reach the value we want, based on the formula we derived before. If we want to cool it down, we can turn down the heater and turn on the valve of cold water. However, these two measures are linked together, because both Q and h will affect the steady temperature. Since cold water can not only balance the temperature, but also raise the water surface, which lead to the decrease of temperature, adding cold water is more efficient. However, the raised water surface can confine the increase of temperature in the future, which is quite troublesome. So we think the best strategy is that if we are in emergency or we want to cool the kettle down quickly, for instance, we found the heater out of control and the temperature inside the kettle is excessively high, we add a large amount of water. On the other side, if we do not hurry, we only turn down the heater to let the kettle naturally cool down to the set temperature. An outlet is also necessary to expel unwanted water in the kettle.
Based on previous discussion, we can safely draw the conclusion that the temperature in the kettle can be well controlled automatically.
MATLAB code is available at: Matlab Code 5
