Team:UNIPV-Pavia/Project/Modelling

From 2011.igem.org

Revision as of 16:04, 17 September 2011 by Fedesampi (Talk | contribs)

UNIPV TEAM 2011

Modelling

Contents



Mathematical modelling: introduction

Mathematical modelling plays a central role in Synthetic Biology, due to its ability to serve as a crucial link between the concept and realization of a biological circuit: what we propose in this page is a mathematical modelling approach to the entire project, which has proven extremely useful before and after the "wet lab" activities.

Thus, immediately at the beginning, when there was little knowledge, a mathematical model based on a system of differential equations was derived and implemented using a set of reasonable values of model parameters, to validate the feasibility of the project. Once this became clear, starting from the characterization of each simple subpart created in the wet lab, some of the parameters of the mathematical model were estimated thanks to several ad-hoc experiments we performed within the iGEM project (others were derived from literature) and they were used 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.

After a brief overview on the importance of the mathematical modelling approach, we deeply analyze the system of equations, underlining the role and function of the parameters involved.
Experimental procedures for parameter estimation are discussed and, finally, a different type of circuit is presented. Simulations were performed, using ODEs with MATLAB and used to explain the difference between a closed-loop control system model and an open one.


The importance of mathematical modelling

Mathematical modeling reveals fundamental in the challenge for understanding and engineering complex biological systems. Indeed, these are characterized by a high degree of interconnection among the single constituent parts, requiring a comprehensive analysis of their behavior through mathematical formalisms and computational tools

Synthetically, we can identify two major roles concerning mathematical models:
  • Simulation: mathematical models allow to analyse complex system dynamics and to reveal the relationships between the involved variables, starting from the knowledge of the single subparts behavior and from simple hypotheses of their interconnection. (Endler et al, 2009)
  • Knowledge elicitation: mathematical models summarize into a small set of parameters the results of several experiments (parameter identification), allowing a robust comparison among different experimental conditions and providing an efficient way to synthesize knowledge about biological processes. Then, through the simulation process, they make possible the re-usability of the knowledge coming from different experiments, engineering complex systems from the composition of its constituent subparts under appropriate experimental/environmental conditions (Braun et al, 2005; Canton et al, 2008).

Equations for gene networks

Below is provided the system of equations of our mathematical model.




Hypotheses of the model

HP1: in equation (2) only HSL is considered as inducer, instead of the complex LuxR-HSL. This is motivated by the fact that our final device offers a constitutive LuxR production due to the upstream constitutive promoter Pλ. Assuming LuxR is abundant in the cytoplasm, we can understand this simplification of attributing pLux promoter induction only by HSL.

HP2: in system equation, LuxI and AiiA amounts are expressed per cell. For this reason, the whole equation (3), 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 freely inside/outside bacteria.

HP3: as regards promoters pTet and pLux, we assume their strengths (measured in PoPs), due to a given concentration of inducer (aTc and HSL for Ptet and Plux respectively), to be independent from the gene downstream. In other words, in our hypothesis, if the mRFP coding region is substituted with a region coding for another gene (in our case, AiiA or LuxI), we would obtain the same synthesis rate: this is the reason why the strength of the complex promoter-RBS is expressed in Arbitrary Units [AUr].

HP4: considering the exponential growth, the enzymes AiiA and LuxI concentration is supposed to be constant, because their production is equally compensated by dilution.


Equations (1) and (2)

Equations (1) and (2) have identical structure, differing only in the parameters involved. They represent the synthesis, degradation and diluition of both the enzymes in the circuit, LuxI and AiiA, respectively in the first and second equation: in each of them both transcription and translation processes have been condensed. The mathematical formulation is analogous to the one used by Pasotti et al 2011, Suppl. Inf., even if we do not take LuxR-HSL complex formation into account, as explained below.

These equations are composed of 2 parts:

  1. The first term describes, through Hill's equation, the synthesis rate of the protein of interest (either LuxI or AiiA) depending on the concentration of the inducer (anhydrotetracicline -aTc- or HSL respectively), responsible for the activation of the regulatory element composed of promoter and RBS. In the parameter table (see below), α refers to the maximum activation of the promoter, while δ stands for its leakage activity (this means that the promoter is slightly active even if there is no induction). In particular, in equation (1), the almost entire inhibition of pTet promoter is given by the constitutive production of TetR by our MGZ1 strain. In equation (2), pLux is almost inactive in the absence of the complex LuxR-HSL.

    Furthermore, in both equations k stands for the dissociation constant of the promoter from the inducer (respectively aTc and HSL in eq. 1 and 2), while η is the cooperativity constant.


  2. The second term in equations (1) and (2) is in turn composed of 2 parts. The former one (γ*LuxI or γ*AiiA respectively) describes, with an exponential decay, the degradation rate per cell of the protein. The latter (μ*(Nmax-N)/Nmax)*LuxI or μ*(Nmax-N)/Nmax)*AiiA, respectively) takes into account the dilution factor against cell growth which is related to the cell replication process.

Equation (3)

Here the kinetics of HSL is modeled, through enzymatic reactions either related to the production or the degradation of HSL. This equation is composed of 3 parts:

  1. The first term represents the production of HSL due to LuxI expression. We modeled this process with a saturation curve in which Vmax is the HSL maximum transcription rate, while kM,LuxI is LuxI dependent half-saturation constant.

  2. The second term represents the degradation of HSL due to the AiiA expression. Similarly to LuxI, kcat represents the maximum degradation per unit of HSL concentration, while kM,AiiA is the concentration at which AiiA dependent HSL degradation rate is (kcat*HSL)/2. The formalism is similar to that found in the Supplementary Information of Danino et al, 2010.

  3. The third term (γHSL*HSL) is similar to the corresponding ones present in the first two equations and describes the intrinsic protein degradation.


Equation (4)

This is the typical cells growth equation, depending on the rate μ and the maximum number Nmax of cells per well reachable (Pasotti et al, 2009).


Table of parameters and species


Parameter & Species Description Unit of Measurement Value
αpTet maximum transcription rate of pTet (dependent on RBSx efficiency) [(AUr/min)/cell] -
δpTet leakage factor of promoter pTet basic activity [-] -
ηpTet Hill coefficient of pTet [-] -
kpTet dissociation constant of aTc from pTet [nM] -
αpLux maximum transcription rate of pLux ((dependent on RBSx efficiency) [(AUr/min)/cell] -
δpLux leakage factor of promoter pLux basic activity [-] -
ηpLux Hill coefficient of pLux [-] -
kpLux dissociation constant of HSL from pLux [nM] -
γpLux LuxI constant degradation [1/min] -
γAiiA AiiA constant degradation [1/min] -
γHSL HSL constant degradation [1/min] -
Vmax maximum transcription rate of LuxI [nM/(min*cell)] -
kM,LuxI half-saturation constant of LuxI from HSL [AUr/cell] -
kcat maximum number of enzymatic reactions catalyzed per minute [1/(min*cell)] -
kM,AiiA half-saturation constant of AiiA from HSL [AUr/cell] -
Nmax maximum number of bacteria per well [cell] -
μ rate of bacteria growth [1/min] -
LuxI kinetics of LuxI enzyme [AUrcell] -
AiiA kinetics of AiiA enzyme [AUrcell] -
HSL kinetics of HSL [nM(min)] -
N number of cells cell -


NOTE: In order to better investigate the range of dynamics of each subpart, every promoter has been studied with 4 different RBSs, so as to develop more knowledge about the state variables in several configurations of RBS' efficiency (Salis et al, 2009). Hereafter, referring to the notation "RBSx" we mean, respectively, RBS30, RBS31, RBS32, RBS34.

Parameter CV

Parameter & Species BBa_B0030 BBa_B0031 BBa_B0032 BBa_B0034
αpTet 2.13 imm 10.9 12.2
δpTet 888.5 1.47*10^11 411.3 61.74
ηpTet 12 imm 76671474 2705172
kpTet 2.54 imm 393180 55889
αpLux 10.14 7.13 2.78 5.8
δpLux 179.7 57.04 1317.7 187.2
ηpLux 47.73 29.13 9.75 19.3
kpLux 27.5 25.81 8.46 17.86


Parameter estimation

The aim of the model is to predict the behavior of the final closed loop circuit starting from the characterization of single BioBrick parts through a set of well-designed ad hoc experiments. This section presents the experiments performed. As explained before in NOTE, considering a set of 4 RBSs for each subpart expands the range of dynamics and helps us to better understand the interactions between state variables.

Promoter (PTet & pLux)

These are the first parts tested, with the target of learning more about pTet and pLux promoters. In particular, as previously explained in NOTE, for each promoter, we tested four different combinations of promoter-RBS, providing us a set of fundamental building blocks for the subsequent assebly of the closed-loop circuit.
As shown in the figure below, we considered a range of inductions and we monitored, in time, absorbance (O.D. stands for "optical density") and fluorescence; the two vertical segments for each graph highlight the exponential phase of bacterial growth. Scell (namely, synthesis rate per cell) 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:
However, also Relative Promoter Unit (RPU, Kelly 2009) has been calculated as a ratio of Scell of promoter of interest and the Scell of BBa_J23101 (reference to HP3).
As shown in the figure, α, as already mentioned, represents the protein maximum synthesis rate, which is reached, in accordance with Hill equation, when the inducer concentration tends to infinite, and, more practically, when the inducer concentration is sufficiently higher than the dissociation constant. Meanwhile the product α*δ stands for the leakage activity (at no induction), liable for protein production (LuxI and AiiA respectively) even in the absence of inducer. The paramenter η is the Hill's cooperativity constant and it affects the ripidity of transition from the lower and upper boundary of the curve relating Scell to the inducer concentration. 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.

T9002 introduction

LuxI and AiiA tests have been always performed exploiting the well-characterized BioBrick BBa_T9002, by which it's possible to quantify exactly the concentration of HSL.
This is a biosensor which receives HSL concentration as input and returns GFP intensity (more precisely Scell) as output. (Canton et al, 2008). According to this, it is necessary to understand the input-output relationship: so, a T9002 "calibration" curve is plotted for each test performed.


AiiA & LuxI

This paragraph explains how parameters of equation (3) are estimated. The target is to learn the AiiA and LuxI degradation and production mechanisms in addiction to HSL intrinsic degradation, in order to estimate Vmax, KM,LuxI, kcat and KM,AiiA parameters. These tests have been performed using the following BioBrick parts:
By now, parameter identification about promoters has already been performed. Furthermore, as explained before, the HP3 is also valid in this case. So, our idea is to control the degradation of HSL in time. ATc activates pTet and, later, a certain concentration of HSL is introduced. Then, at fixed times, O.D.600 and HSL concentration are monitored using Tecan and T9002 biosensor.
Referring to HP4, in exponential growth enzymes equilibrium is conserved. Due to a known induction of aTc, the steady-state level per cell can be calculated:
Considering, for a determined promoter-RBSx couple, several induction of aTc and, for each of them, several samples of HSL concentration during time, parameters Vmax, kM,LuxI, kcat and kM,AiiA can be estimated, through numerous iterations of an algorithm implemented in MATLAB.

N

The parameters Nmax and μ can be calculated from the analysis of the OD600 produced by our MGZ1 culture. In particular, μ is derived as the slope of the log(O.D.600) growth curve. Counting the number of cells of a saturated culture would be considerably complicated, so Nmax is determined with a proper procedure. The aim here is to derive the linear proportional coefficient Θ between O.D'.600 and N: this constant can be estimated as the ratio between absorbance (read from TECAN) and the respective number of CFU on a petri plate. Finally, Nmax is calcultated as Θ*O.D'.600 (Pasotti et al, 2010).


Degradation rates

The parameters γLuxI and γAiiA are taken from literature since they contain LVA tag for rapid degradation. Instead, approximating HSL kinetics as a decaying exponential, γHSL can be derived as the slope of the log(concentration), which can be monitored through BBa_T9002.


Simulations

On a biological level, the ability to control the concentration of a given molecule reveals itself as fundamental in limiting the metabolic burden of the cell; moreover, in the particular case of HSL signalling molecules, this would give the possibility to regulate quorum sensing-based population behaviours. In this section we present some simulations of another circuit, which could validate the concept of the closed-loop model we have discussed so far.
In order to see that, we implemented and simulated in Matlab an open loop circuit, similar to 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 hypotheses, 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.

Steady state of enzyme expression


References

  1. Braun D, Basu S, Weiss R (2005) Parameter estimation for two synthetic gene networks: a case study ICASSP '05 5:v/769-v/772.

  2. Canton B, Labno A, Endy D. (2008) Refinement and standardization of synthetic biological parts and devices. Nat Biotechnol. 26(7):787-93.

  3. Danino T, Mondragón-Palomino O, Tsimring L et al. (2010) A synchronized quorum of genetic clocks. Nature. 463(7279):326-30.

  4. Endler L, Rodriguez N, Juty N et al. (2009) Designing and encoding models for synthetic biology. J. R. Soc. Interface 6:S405-S417.

  5. Kelly JR, Rubin AJ, Davis JH et al. (2009) Measuring the activity of BioBrick promoters using an in vivo reference standard. J. Biol. Eng. 3:4.

  6. Pasotti L, Quattrocelli M, Galli D et al. (2011) Multiplexing and demultiplexing logic functions for computing signal processing tasks in synthetic biology. Biotechnol. J. 6(7):784-95.

  7. Salis H M, Mirsky E A, Voight C A (2009) Automated design of synthetic ribosome binding sites to control protein expression. Nat. Biotechnol.27:946-950.

Retrieved from "http://2011.igem.org/Team:UNIPV-Pavia/Project/Modelling"