Team:ETH Zurich/Model/Environment Sensing/parameter space

Model-based design guidelines used to choose our parts rationally

Why do we need to search for functional parameters?

A biological circuit often has different functioning regimes and can only achieve a particular and interesting behavior for a few given combinations of parameters (among which protein expression levels, or promoter sensitivity and leakiness for example). This is why we had to get some insights into the sets of parameters of our circuit and determine a target combination that would make our circuit work, to give out design guidelines for the choice of the parts we would use in the lab.

Please check the full detailed model if you are interested in knowing how our whole model works.

Parameter search

Even before we got our first parts cloned and characterized, we attempted to predict the requirements that they should meet to achieve the criteria previously established from literature data. For this, we extensively explored the parameter space controlling our model, and simulated the response of potential systems to find subsets of parameter combinations satisfying the performance specifications needed to get a sensitive as well as specific tumor sensing circuit.

Modeling process - Parameter search
Figure 1. Goal of the parameter search: deduce genetic design guidelines

The parameters satisfying the specifications are called functional space (green ellipse on Figure 1). We selected the biological parts that were the most likely to fall in the functional parameter space.

Different categories of parameters

Our model is relying on a dozen of parameters, some of which we can have a leverage on (typically maximal expression of the proteins, via RBS tuning), and others not (binding constants, promoter leakiness...). Some of these latter parameters have been precisely characterized and others are not very well known. This is why we have chosen to set some parameters to a certain value when we could find a reasonably reliable source in the literature, or when their influence would be redundant with other parameters (such as protein degradation rates and maximal expression rates that can alleviate each other's influence when co-varied), and leave other parameters free to vary to check their influence on our system.

Fixed parameters, because well known

Symbol Description Value Reference
aAHL AHL synthesis rate by LuxI 0.01 min-1 [2]
dAHL AHL degradation rate 5x10-4 min-1 [4]
D AHL diffusion coefficient in water 3x10-8 m2min-1 [4]
KLuxR-AHL LuxR-AHL quadrimer binding constant 5x10-10 nM-3 [5]
w Width of the colonized shell area 5x10-10 nM-3 [5]

Fixed parameters, not very well known but redundant with other parameters

Symbol Description Value Reference
dLuxI LuxI degradation rate 0.017 min-1 [4]
dLuxR LuxR degradation rate 0.023 min-1 [4]
dAzu Azurin degradation rate 0.1 min-1 estimated

Parameters allowed to vary because not very well known and which may have a significant effect on our circuit

Symbol Description Typical value (initial value in the parameter search) Reference Lower bound Higher bound
aLuxR Maximum expression of luxR 5 nM min-1 iGEM ETH 2014 1x10-2 nM min-1 1x104 nM min-1
aLuxI Maximum expression of luxI 1x103 nM min-1 [5] 1x10-2 nM min-1 1x104 nM min-1
KLac Half-activation lactate concentration of the hybrid promoter 2x106 nM Characterized lactate sensing part on which our AND-gate is based 1x104 nM 1x108 nM
kLuxI Leakiness of the hybrid promoter 0.01 Characterized lactate sensing part on which our AND-gate is based 0.0001 0.1
KLuxR Half-activation LuxR-AHL concentration of the hybrid promoter 5 nM iGEM ETH 2013 1 nM 100 nM
nLuxR Hill coefficient of the hybrid promoter regarding LuxR-AHL concentration 1.7 iGEM ETH 2015 1.1 1.9
nLac Hill coefficient of the hybrid promoter regarding lactate concentration 1.7 iGEM ETH 2015 1.1 1.9
kAzu-LuxI Relative expression of azurin compared to LuxR 10 times the luxI expression estimated 10-5 105
How do we quantitatively define the functional space?

Cost function

To be able to distinguish systems satisfying the criteria about specificity and azurin production and the ones that do not, we need to use a numerically evaluable condition quantifying how well the criteria are met. Based on this, the script will either accept or discard the parameter set. For this, we will use the following cost function, taking a parameter vector as argument:

\[\text{cost}(p) = \max \left( \frac{10\cdot {[Azu]_{\mathrm{(low\; lac,HIGH\; d_{cell})}}}}{[Azu]_{\mathrm{(HIGH\; lac,HIGH\; d_{cell})}}}, \frac{10\cdot {[Azu]_{\mathrm{(HIGH\; lac,low\; d_{cell})}}}}{[Azu]_{\mathrm{(HIGH\; lac,HIGH\; d_{cell})}}}, \frac{1\times10^{6}}{[Azu]_{\mathrm{(HIGH\; lac,HIGH\; d_{cell})}}}\right)\]

Each argument of the max function represents in the same order the following criteria:

  1. Specificity of the sensing for lactate
  2. Specificity of the sensing for bacterial cell density
  3. Achieving a large amount of produced azurin

      4. Interpretation of the result of the cost function goes as follows: the smaller the value the better better the criterium is met. If the highest value (e.g. the value for the criterium is met the worst) is below 1, the paramter set is accepted. Over 1, a ratio is not good enough. This monotonicity enables us to rely on optimization algorithms to reach the best combination of parameters available. Also, we can say that every system that has a cost function value below 1 is good enough for us, while "the smaller the better" still applies.

Intermediate modeling result: Analysis of the functional parameter space

Using an optimization toolbox developed for biological systems, MEIGO [6], followed by a package exploring parameter spaces, HYPERSPACE [7], we could obtain the following graphs describing, in the high-dimension space of all possible circuits, a subset of systems satisfying our performance criteria:

Parameter search

On this figure are shown the systems suitable for our application. All the axis are logarithmic, except for nLac and nLuxR. The yellow points are good systems, the blue ones are even better and surpass the specifications that we demand. From this figure, we can draw the following interpretations (see corresponding sub-graphs referred to on the figure).

  1. Only some given combination of expression of LuxI and LuxR are suitable for our needs. This is expected as the tuning of the bacterial cell density at which the quorum sensing is triggered is mainly done with these two proteins
  2. High amounts of azurin are more easily achieved when LuxI maximal expression is high: then the expression of azurin does not need to be that much more compared to luxI to reach the desired level.
  3. The tipping point of the lactate sensing must be either around or above the lactate levels to be distinguished (1 mM in healthy tissues and 5mM in tumors). The first possibility makes sense as the promoter should ideally be unactivated at low lactate level and activated above. However, the combination of this lactate sensing and quorum sensing into the hybrid promoter seems to allow for a second possibility: that the full activation of the promoter happens at much higher concentrations. In both cases, the differential expression at 1 mM and 5 mM plays the role of "increasing the leakiness" of the promoter in regard to luxR so that the quorum sensing is more easily activated in presence of lactate.
  4. The leakiness is a very important parameter to be able to achieve a good performance for our system. The smaller the leakiness, the more probable it is to find a good system.
  5. The Hill coefficient of our hybrid promoter in regard to lactate will allow more or less possibilities of systems: when over 1.5, a population of systems is present (more on the yellow side) that allows for a larger set of aLuxR/aLuxI combinations (see also nLac vs aLuxR and nLac vs aLuxI graphs). As we won't be able to tune it, we should prepare for the worst and try to aim for the best systems (the blue ones) on graph 1 to keep a security margin
  6. KLuxR and nLuxR don't have a significant influence on our system, we can stop studying them

Final modeling result: Experimental guidelines used for circuit design

From these observations, we can deduce guidelines regarding the parameters on which we can exert an active control, that is to say the expression level of the genes LuxI and LuxR (aLuxR and aLuxI here) as well as a judicious choice of a previously characterized lactate sensor circuit (comprising lldR and lldP genes) among the iGEM ETH 2015 part collection .

Target parameters and restrictions

To translate these insights into experimental results in the lab, we need to chose a target in the range of parameters that works for our application. With the help of the previously characterized initial values for aLuxR (5 nMmin-1) and aLuxI (1x103 nMmin-1), we can hope to tune our system and reach our target in the parameter space via simple RBS tuning which is supported by the Salis Lab RBS Calculator [7].

As it turned out, the regulatory sequence in front of the luxR gene on the part at our disposal induced already a relatively high expression level. It was hard to get more than 10 times more expression for this gene on the Salis calculator, this is why the range aLuxR > 1x102 nMmin-1 is inaccessible to us (grey area), and that we have to choose LuxI in consequence. We also get to chose KLac among the ones available in the promoter collection of parts ranging from BBa_K1847002 to BBa_K1847009: between 0.3 mM and 2.4mM.

Taking into account the experimental constraints (forbidden grey area), the targeted parameters (red squares) were chosen on the following plot, with an extensive compatibility for different potential leakiness of our hybrid promoter (red frame):

Parameter search iteration 2
Parameter space of vectors of parameters satisfying our criteria (cost below 1). Yellow points are good enough systems, blue ones are even better and surpass the specifications that we demand. Red features highlight our choice for the operating choice that we will implement in the lab through genetic design guidelines.

With aLuxR = 1x102 nMmin-1, aLuxI = 1x104 nMmin-1 and KLac = 1x106 nM, we should be at a suitable operating point for our system and still have some security margin in case the genetic design does not yield the exact expression levels that we would expect from it. To achieve these parameters, we gave the following directions for the design of our parts:

  1. Use a 10 times stronger RBS than on the piG0047 sequence of iGEM ETH 2014 team for the expression of LuxR
  2. Use a 10 times stronger RBS than on the piG0050 sequence of iGEM ETH 2014 team for the expression of LuxI
  3. Use the BBa_K1847008 part with J23118-B0034 regulatory sequences, giving Klac = 1.8 mM

These value were the basis for the design of our parts and the subsequent experimentations.

Sanity check: in silico behavior for the chosen target parameters We can validate on our model that they would work well to distinguish the specific levels dictated by our application:

System response after optimization
Expected behavior of the circuit for the chosen parameter set target: intracellular azurin concentration (in nM) depending on both lactate level and bacterial cell density. The white lines correspond to low levels of lactate and bacterial cell density (in healthy tissues) and the black lines represent the high levels (in tumor tissues). The output level of azurin is only significant if both inputs are high, which is the behavior that suits our application.

We can confirm that the obtained parameter set target would lead to a functioning circuit.

Going further: a 3D tumor model

To assess the validity of this analytical model, and especially validate the assumptions regarding the AHL diffusion model, we developed a comprehensive 3D model of a colonized tumor on the software COMSOL, on which we implemented the same equations governing our bacterial circuit.

References

  1. Kaplan HB, Greenberg EP. Diffusion of autoinducer is involved in regulation of the Vibrio fischeri luminescence system. Journal of Bacteriology. 1985;163(3):1210-1214.
  2. Jordi Garcia-Ojalvo, Michael B. Elowitz, and Steven H. Strogatz Modeling a synthetic multicellular clock: Repressilators coupled by quorum sensing PNAS 2004 101 (30) 10955-10960
  3. Fekete, A., Kuttler, C., Rothballer, M., Hense, B. A., Fischer, D., Buddrus-Schiemann, K., Lucio, M., Müller, J., Schmitt-Kopplin, P. and Hartmann, A. (2010), Dynamic regulation of N-acyl-homoserine lactone production and degradation in Pseudomonas putida IsoF. FEMS Microbiology Ecology, 72: 22–34. doi:10.1111/j.1574-6941.2009.00828.x
  4. A.B. Goryachev, D.J. Toh T.Lee, Systems analysis of a quorum sensing network: Design constraints imposed by the functional requirements, network topology and kinetic constants Biosystems, Volume 83, Issues 2–3, February–March 2006, Pages 178-187
  5. A synthetic multicellular system for programmed pattern formation Subhayu Basu, Yoram Gerchman, Cynthia H. Collins, Frances H. Arnold & Ron WeissNature 434, 1130-1134 (28 April 2005) | doi:10.1038/nature03461
  6. Egea JA, Henriques D, Cokelaer T, Villaverde AF, MacNamara A, Danciu DP, Banga JR and Saez-Rodriguez J. (2014) MEIGO: an open-source software suite based on metaheuristics for global optimization in systems biology and bioinformatics. BMC Bioinformatics 15:136.
  7. Zamora-Sillero E, Hafner M, Ibig A, Stelling J, Wagner A (2011) Efficient characterization of high-dimensional parameter spaces for systems biology. BMC Syst Biol 5: 142.http://doi.org/10.1186/1752-0509-5-142
  8. Salis, Howard M., Ethan A. Mirsky, and Christopher A. Voigt. "Automated design of synthetic ribosome binding sites to control protein expression."Nature biotechnology 27.10 (2009): 946-950. http://dx.doi.org/10.1038/nbt.1568