CTRL + E
Signalling is nothing without control...
Contents |
Mathematical modelling: introduction
Thus, immediately at the beginning, when there was little knowledge, a mathematical model based on a system of differential equations was derived and performed using a set of reasonable parameters, so as to validate the feasibility of the project. Once this became clear, starting from the characterization of the single subparts created in the wet lab, some of the parameters of the mathematical model were estimated (the others are known from literature) and they have been implemented in the same model, in order to predict the final behaviour of the whole engineered closed-loop circuit. This approach is consistent with the typical one adopted for the analysis and synthesis of a biological circuit, as exemplified by Pasotti et al 2011.
Therefore here, after a brief overview about the advantages that modelling engineered circuits can bring, we deeply analyze the system of equation formulas, underlining the role and the function of the parameters involved.
Experimental procedures for parameters estimation are discussed and, finally, a different type of circuit is presented and simulations performed, using ODEs with MATLAB and explaining the difference between a closed-loop model and an open one.
The importance of mathematical modelling
Equations for gene networks
Equations (1) and (2)
The equations are composed of 2 parts:
- The first term describes, through Hill's equation formalism, the synthesis rate of the protein of interest (either LuxI or AiiA) depending on the concentration of the inducible protein (anhydrotetracicline -aTc- or HSL respectively). As can be seen in the parameters table (see below), α refers to the maximum activation of the promoter, δ stands for its leakage activity (this means that the promoter is quite active even if there is no induction). In particular, in equation (1), the quite total inhibition of pTet promoter is due to the constitutive production of TetR by our MGZ1 strain, while, in equation (2), pLux is almost repressed in the absence of the complex LuxR-HSL.
In equation (2) only HSL seems to be the inducer, instead of the complex LuxR-HSL. This is motivated by the fact that our final device offers a constitutive production of LuxR (due to the upstream constitutive promoter pLac), so that, assuming it abundant in the cytoplasm, we can derive the semplification of attributing pLux promoter induction only by HSL: this is the reason why we didn' t consider LuxR in the equations system as well as LuxI and AiiA. Furthermore, in both equations k stands for the dissaciation constant of the promoter from the inducer (respectively aTc and HSL in eq. 1 and 2), while η is the cooperativity constant.
The second term in equations (1) and (2) is in turn composed of 2 parts. The first one (γ*LuxI or γ*AiiA respectively) describes, with a linear relation, the degradation rate per cell of the protein. The second one (μ*(Nmax-N)/Nmax)*LuxI or μ*(Nmax-N)/Nmax)*AiiA, respectively) takes into account the dilution term due to cell growth and is related to the cell replication process.
Equation (3)
3 parts have been identified in this equation:
- The first term represents the production of HSL due to LuxI expression. We model this process with saturation curve in which Vmax is HSL maximum transcription rate, while KM is the dissociation constant of LuxI from the substrate HSL and it represents the concentration of LuxI at which HSL synthesis rate is Vmax/2.
- The second term represents the degradation of HSL due to the AiiA expression. Similarly to LuxI, Kcat represents maximum degradation per unit of HSL concentration, while KM1 is the concentration at which AiiA dependent HSL concentration rate is (Kcat*HSL)/2. The formalism is similar to that found in the Supplementary Information of Danino et al, 2010.
- The third term (γHSL*HSL) is similar to the corresponding ones present in the first two equations and describes the intrinsic protein degradation.
NOTE2: the whole equation, except for the term of intrinsic degradation of HSL, is multiplied by the number of cells N, due to the property of the lactone to diffuse free inside/outside bacteria. Notice that, in system equation, LuxI and AiiA amounts are expressed per cell.
Equation (4)
Table of parameters and species
Parameter & Species | Description | Unit of Measurement | Value |
αpTet | maximum transcription rate of pTet (related with RBSx efficiency) | [(AUr/min)/cell] | - |
δpTet | leakage factor of promoter pTet basic activity | [-] | - |
ηpTet | Hill coefficient of pTet | [-] | - |
kpTet | dissociation costant of aTc from pTet | [nM] | - |
αpLux | maximum transcription rate of pLux (related with RBSx efficiency) | [(AUr/min)/cell] | - |
δpLux | leakage factor of promoter pLux basic activity | [-] | - |
ηpLux | Hill coefficient of pLux | [-] | - |
kpLux | dissociation costant of HSL from pLux | [nM] | - |
γpLux | LuxI costant degradation | [1/min] | - |
γAiiA | AiiA costant degradation | [1/min] | - |
γHSL | HSL costant degradation | [1/min] | - |
Vmax | maximum transcription rate of LuxI | [nM/(min*cell)] | - |
KM | dissociation costant of LuxI from HSL | [AUr/cell] | - |
Kcat | maximum number of enzymatic reactions catalysed per minute | [1/(min*cell)] | - |
KM1 | dissociation costant of AiiA from HSL | [AUr/cell] | - |
Nmax | maximum number of bacteria per well | [cell] | - |
μ | rate of bacteria growth | [1/min] | - |
LuxI | kinetics of enzyme LuxI | [AUr⁄cell] | - |
AiiA | kinetics of enzyme AiiA | [AUr⁄cell] | - |
HSL | kinetics of HSL | [nM⁄(min)] | - |
N | number of cells | cell | - |
Parameter estimation
Promoter (PTet & pLux)
As shown in the box below, we consider a range of induction and we monitor, during the time, absorbance (line1, line2) and fluorescence (line3); the two vertical segments for each figure highlight the exponential phase of bacteria' s groth. Scell (explained few lines below) can be derived as a function of inducer concentration, thereby providing the desired input-output relation (inducer concentration versus promoter+RBS activity), which was modelled as a Hill curve.
After that, we can calculate the Scell as: In the end, plotting Scell VS induction, we obtain the activation Hill curve of the considered promoter. As shown in the box above, α as already mentioned, represents the protein maximum synthesis rate, which is reached, in accordance with Hill's formalism, when the inducer concentration tends to infinite, and, more practically, for sufficently high concentrations of inducer, meanwhile the product α*δ stands for the leakage activity (at no induction), liable for protein production (LuxI and AiiA respectively) even in the absence of autoinducer. The paramenter η is the Hill's cooperativity constant and it affects the rapidity and ripidity of the switch like curve relating Scell with the concentration of inducer. Lastly, k stands for the semi-saturation constant and, in case of η=1, it indicates the concentration of substrate at which half the synthesis rate is achieved. The unities of the various parameters can be easily derived considering the Hill equation(for more details see the Table of parameters above).
AiiA & LuxI
So, our idea is to control the degradation of HSL in time. aTc activates pTet and, after having waited enough for the enzyme to become in stationary phase, a certain concentration of HSL is given. Then, in precise time samples absorbance is controlled and HSL concentration monitored, reading the fluorescence of T9002. Now, considering the exponential growth, the concentration of the AiiA and LuxI is supposed to be constant: after the cell division, fewer enzyme is present into the single bacteria but, on the other hand, the number of cells has increased, and so the enzyme equilibrium is conserved. Due to a well-known induction of aTc, the steady-state level per cell can be calculated: Then considering, for the same match of promoter and RBS, several induction of aTc and, for each of it, several samples of HSL concentration during the time, parameters Vmax, KM, Kcat and KM1 can be estimated, through numerous iterations of an algorithm that includes the functions lsqnonlin and ODE of MATLAB.
N
Degradation rates
Simulations
In order to see that, we implemented and simulated in Matlab an open loop circuit, totally like CTRL+E, except for the constitutive production of AiiA.
Sensitivity Analysis of the steady state of enzyme expression in exponential phase
In this paragraph we want to theoretically investigate our circuit behaviour in cell's culture exponential growth phase. According to this, we first derive, under feasible hypoteses, the steady state condition for our enzymes concentration. Then we perform a sensitivity analysis relating the output of our system (HSL) to input (atc) and system parameters.
Referring to the system of equations, it can be easily noted that it provides two equilibrium points, associated with the values 0 and Nmax for the state variable N. Moreover, as can be seen in simulations paragraph, there is also a region of functioning in which the system appears to be subject to small changes over time. Indeed, during the cellular population exponential growth phase, the enzymes (and also HSL) undergo only minor changes in their values. This can be mathematically demonstrated, if we make the feasible hypothesis to have N<<Nmax, surely verified since in the exponential growth phase N is of the order of 10^3, while Nmax=2*10^9. In addiction, in equation 3, we can disregard the third term (-γ*HSL) compared to the others, since that it is the only one not multiplied by N, and γ is lower than 10^-2 1/min. Under this assumptions, we can approximate the exponential growth phase as another steady state for the first three variables system:
Now, the information brought by the fourth equation is simply that we are in exponential growth phase and we can neglect it afterwards. Second equation is independent from the third and fourth ones, enabling us to directly determine LuxI steady state expression during the exponential growth phase. On the contrary, third and fourth equations depend each other in defining the value of AiiA and HSL respectively, because the former is a function of HSL, while the latter is a function of AiiA. So we could resolve a system of two equations, first by expliciting one of the two variables with respect to the other, and then substituting its expression in order to determine the other variable possible values. This would bring a complex mathematical formulation, which is not helpful in understanding the influence of the various model parameters on the output HSL. On the other hand, AiiA and HSL values can also be graphically determined from the intersection of the curves derived from these two equations, if we explicit HSL as a function of AiiA (or, alternatively, AiiA as a function of HSL). It is easy to discover that these two curves represent rectangular hyperbolae (the first one only under a simple approximation) whose tails intersect each other at a particular point, corresponding to the searched values for AiiA and HSL.
centered at O(-d/c;a/c), with the vertical asymptote x=-d/c and the horizontal asymptote y=a/c
From equation 2, we have:
We can introduce the simplification to remove the ηplux exponent to the entire expression in the right hand side of the equation, obtaining a rectangular hyperbola:
Parameter | Expression |
a | |
b | |
c | |
d |
The table below provides the corresponding asymptotes
Horizontal asymptote | Vertical asymptote |
As pertains to equation 3, its steady state approximation during the exponential growth is more immediately identifiable as a rectangular hyperbola. In particular, the table below clarifies it has the ordinate axis as its vertical asymptote.
a | b | c | d | Ver. asym. | Hor. asym. |
0 | 0 |
Sensitivity Analysis
Now, it is interesting to conduct some qualitative and quantitative considerations about our system sensitivity to its parameters and atc input signal.
From the above figures, it is also clear that HSL is not very sensitive to αpLux, at least when this parameter presents values greater than unity, because this brings the two curves to intersect in their low slope regions.
Referring to the second RH, the only adjustable asymptote is the horizontal one, that we can move upward or downward by altering VLuxI, which in turn depends on LuxI.
A further deepening in the exponential phase analysis involves determining the characteristics of the input-output relation. Our closed loop system can be realized with two alternative purposes in mind:
- realizing a circuit able to adapt HSL output depending on atc input concentration. This requires a good sensitivity between input and output.
- designing a robust HSL concentration controller, which is immune to the input noise and offers a constant and defined amount of HSL. In this case HSL level should be appropriately tuned during the design stage, by choosing the correct strength of promoter-RBSx complexes.
Now,again considering the system of equations, it is easy to observe that HSL dependence on atc input passes through two Hill equations. The former describes atc driven LuxI synthesis, while the latter models LuxI dependent HSL synthesis rate. Therefore, in order to achieve a high atc sensitivity, it is advisable to tune atc and LuxI levels so that they place outside the saturation regions of their Hill curves. In this regard, it is possible to determine a closed form expression relating HSL to atc, if we hypothesize that both atc and LuxI are far lower than their respective half-saturation constant (k_ptet and Km), thus simplifying the Hills with a first order relation: