<h1 class="headline">Model-based design guidelines used to choose our parts rationally</h1>
−
<section class="first">
+
<section class="abstract">
−
<h1>Model of our circuit</h1>
+
<h1 class="goal">Why do we need to search for functional parameters?</h1>
−
<p>The tumor sensing circuit is composed of several proteins that interact with small molecules (AHL and lactate) and DNA (at the transcription factors binding sites). To establish a model describing the behavior of our circuit, we first had to understand the way these interactions are happening inside of the cell. Building on the description of the <a href="https://2017.igem.org/Team:ETH_Zurich/Circuit/Fa_Tumor_Sensor">Tumor Sensor circuit</a>, here is a more detailed overview:</p>
+
<p>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.
<p><strong>Please check the <a href="https://2017.igem.org/Team:ETH_Zurich/Model/Environment_Sensing/model">full detailed model</a> if you are interested in knowing how our whole model works.</strong><p>
−
<p>Let us first focus on the lactate sensing part of the circuit. In the cell, two proteins are produced:</p>
+
−
<ul>
+
−
<li>LldP: a transmembrane protein enabling the transport of extracellular lactate into E. coli.</li>
+
−
<li>LldR: a transcription factor, repressing the activity of the hybrid promoter when not bound to lactate. Lldr releases repression once it binds to lactate.</li>
+
−
</ul>
+
−
<p>To model precisely the regulation of the hybrid promoter by lactate, it would be necessary to take into account all the following points:</p>
−
<ul>
−
<li>How the intracellular lactate concentration behaves in regard to the expression level of LldP and the extracellular lactate concentration</li>
−
<li>What is the binding constant between LldR and the lactate</li>
−
<li>What is the binding dynamics of Lldr to the operon, and how it affects the transcription rate downstream</li>
−
</ul>
−
−
<p>In an effort to simplify our model to reduce it to the most meaningful parameters, and because it has already extensively been studied and characterized by previous iGEM teams, we have chosen not to take into account the complexity of the lactate sensing pathway and rather use a phenomenological model to describe its influence. We rely on the characterization of the lactate sensor using several expression regulation sequences <a href="https://2015.igem.org/Team:ETH_Zurich/Results#Characterization_of_synthetic_promoter_library">done by the ETH 2015 iGEM team</a>. We consider therefore that lactate sensing follows a Hill function as following:</p>
<p>The sensing of the bacterial cell density is done via a quorum sensing circuit. The principles behind quorum sensing is that, via the expression of the enzyme LuxI, each bacteria produces a basal amount of a small molecule (here N-acyl homoserine lactone, AHL) that diffuses in the environment and into neighboring cells. When AHL is present in sufficient quantity, it binds to the intracellular LuxR and induces the production of more LuxI, which in turn results in the production of more AHL. This positive feedback loop results in the activation of the operon containing the LuxI gene when the cell density reaches a critical threshold.
<p>Even before we got our first parts cloned and characterized, we attempted to predict the requirements that they should meet to achieve the criteria <a href="https://2017.igem.org/Team:ETH_Zurich/Model/Environment_Sensing/system_specifications">previously</a> 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.</p>
−
</figure>
+
−
<p>LuxR is under a constitutive promoter of strength <span class="math">\( a_{\text{luxR}} \)</span> and its degradation rate is <span class="math">\( d_{\text{luxR}} \)</span>. <span style="font-variant: small-caps;">AHL</span> binds and stabilizes LuxR; LuxR-<span style="font-variant: small-caps;">AHL</span> molecules can only act as transcription factors when they form a tetramer (2*<span style="font-variant: small-caps;">AHL</span>+2*LuxR). Since we are modeling the steady state, the following simplifications apply:</p>
<li><p>We consider that the total amount of LuxR present in the cell is constant, and only depends on its constitutive expression and degradation rate.</p></li>
+
/> <figcaption>Figure 1. Goal of the parameter search: deduce genetic design guidelines</figcaption>
−
<li><p>We consider the global binding equilibrium between LuxR and <span style="font-variant: small-caps;">AHL</span> without taking into account the intermediary dimers.</p></li>
+
−
</ul>
+
−
+
−
<p>We can therefore write the following equations:</p>
+
−
+
−
<p><span class="math">\[\begin{aligned}
+
−
[\text{LuxR}]_0 &= \frac{a_{\text{LuxR}}}{d_{\text{LuxR}}} & \text{steady state concentration} \\
+
−
[\text{LuxR-AHL}] &= K_{LuxRAHL} [\text{LuxR}]^2 [\text{AHL}]^2 & \text{rapid binding equilibrium} \\
+
−
[\text{LuxR}] &= [\text{LuxR}]_0 - 2 [\text{LuxR-AHL}] & \text{mass conservation}\end{aligned}\]</span>
+
−
</p>
+
−
+
−
<h2>Hybrid Lux-Lac promoter</h2>
+
−
<p>The expression of the main operon containing the LuxI, Bfr and Azurin genes is regulated by a hybrid promoter activated by the quorum sensing and repressed by the lactate sensing (the repression being released in presence of lactate). This hybrid promoter should behave as a AND-gate: mathematically, this corresponds to multiplying the Hill functions describing their behavior.</p>
+
−
+
−
<p>Along with the lactate concentration, the intracellular levels of LuxR-<span style="font-variant: small-caps;">AHL</span> complexes affect LuxI expression. With <span class="math">\( a_{\text{LuxI}} \)</span> being the maximal production rate of LuxI, <span class="math">\( d_{\text{LuxI}} \)</span> the degradation rate, <span class="math">\( k_{\text{LuxI}} \)</span> the leakiness of the promoter and <span class="math">\( P_{\text{Lux-Lac}} \)</span> the combined effect of the <span class="math">\( P_{\text{Lux}} \)</span> and <span class="math">\( P_{\text{Lac}} \)</span> regulating sequences behavior, the ODE governing the production of <span class="math">\( [\text{LuxI}] \)</span> can be written as following:</p>
<p>AHL is produced intracellularly by LuxI and diffuses then freely through the membrane <a href="#bib1" class="forward-ref">[1]</a>. Modeling the production of AHL is quite straightforward: it is proportional to the amount of LuxI present intracellularly. To describe the production per unit of volume though, we have to take into account the bacterial cell density present locally and take it as a dilution coefficient (for instance, if the cells occupy locally half of the volume, then the intracellularly produced AHL would be instantly diluted two times as it diffuses into the surrounding environment).</p>
<p>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.</p>
<summary>Different categories of parameters</summary>
−
+
<p>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.</p>
−
<section>
+
−
<h1>The missing link between AHL production and its concentration: DIFFUSION </h1>
+
−
+
−
<p>To be able to close mathematically the feedback loop, we still miss an equation: we need to know how the intracellular production of AHL translates into the AHL concentration in the environment. For that, we have to consider the diffusion of AHL around the colonized area of the tumor. By solving the equation governing AHL transport, we can, under certain assumptions which will be detailed below, get the relationship between AHL production and its local concentration around bacterial cells.</p>
+
−
+
−
<h2 id="assumption-negligible-ahl-degradation-in-the-layer">Assumption: negligible AHL degradation in the layer</h2>
+
−
+
−
<p>The colony produces <span style="font-variant: small-caps;">AHL</span> which diffuses out of the layer. Because of the symmetry of the problem, we only consider diffusion in the radial direction, being characterized by the diffusion constant <span class="math">\( D = \SI{4.9e-6}{cm^2 s^{-1}} \)</span> (diffusion constant of <span style="font-variant: small-caps;">AHL</span> in water) <a href="#bib4" class="forward-ref">[4]</a>. We also take into account extracellular degradation of <span style="font-variant: small-caps;">AHL</span>, described by the degradation constant <span class="math">\( k_{\text{deg}} = \SI{5e-4}{min^{-1}} \)</span> <a href="#bib4" class="forward-ref">[4]</a>. Using the relationship between the mean square distance and time in a Brownian movement, we can estimate the average time needed for <span style="font-variant: small-caps;">AHL</span> to diffuse out of the layer:</p>
<p>During this time <span class="math">\( \Delta t \)</span>, we can estimate the magnitude of the degradation of <span style="font-variant: small-caps;">AHL</span>. The proportion of AHL that gets degraded before diffusing out of the layer can be estimated as:</p>
<p>Therefore, we can consider for further work that the degradation of AHL happening in the layer is negligible.</p>
+
−
+
−
<h2 id="assumption-only-care-about-local-ahl">Assumption: AHL doesn't diffuse far from where it is produced</h2>
+
−
+
−
<p>To assess whether we would have to consider in a given point of the tumor the AHL coming from every part of the tumor or only the closest area, we have to estimate how far <span style="font-variant: small-caps;">AHL</span> diffuses before being degraded. Considering the half-life of AHL and the mean distance covered in this given time by diffusion:</p>
d &= \sqrt{D t_{1/2}} \simeq \SI{6}{mm}\end{aligned}\]</span></p>
+
−
+
−
<p>Since <span class="math">\( d < r_1 = \SI{10}{mm} \)</span>, <span style="font-variant: small-caps;">AHL</span> will be substantially degraded before it reaches a colonized area far from where it is produced. We will thus assume that we can restrict the study of the diffusion of AHL to a local one, without considering the AHL that could come from the other side of the tumor.</p>
+
−
</section>
+
−
+
−
<section>
+
−
<h1 id="simplified-model-of-ahl-diffusion">Model of AHL diffusion</h1>
+
−
+
−
<p>Given the latter assumption, and given <span class="math">\( w \ll r_1 \)</span>, we assume that the colonized area appears locally as an infinite sheet of width <span class="math">\( w \)</span>. According to the first assumption stated above, we consider only <em>production and diffusion</em> to be significant within the sheet, and <em>diffusion and degradation</em> to be significant outside the sheet. We model the sheet as perpendicular to the <span class="math">\( x \)</span> axis, centered at <span class="math">\( 0 \)</span>, covering the interval <span class="math">\( [-w/2, w/2 ] \)</span>.</p>
+
−
+
−
<p>For a small parallelepiped of surface <span class="math">\( \mathrm{d} S \)</span> and width <span class="math">\( \mathrm{d} x \)</span> perpendicular to <span class="math">\( x \)</span> axis (see figure below), we can use <em>Fick’s law</em> to model diffusion of <span style="font-variant: small-caps;">AHL</span>.</p>
<p><span class="math">\[\Phi(x) = - D \left. \frac{\mathrm{d} [\text{AHL}]}{\mathrm{d} x} \right|_{x}\]</span></p>
+
−
<p>When inside (<span class="math">\( x < |w/2| \)</span>) of the layer, the change of protein amount within time length <span class="math">\( \mathrm{d} t \)</span> is equal to diffusive transports plus production:</p>
+
−
+
−
<p><span class="math">\[\begin{aligned}
+
−
\mathrm{d} n &= (( \Phi(x) - \Phi(x + \mathrm{d} x) ) \mathrm{d} S + P_{\text{AHL}} \mathrm{d} V) \, \mathrm{d} t \\
+
−
\frac{\mathrm{d} n}{\mathrm{d} V \mathrm{d} t} &= D \frac{\mathrm{d}^2 [\text{AHL}]}{\mathrm{d} x^2} + P_{\text{AHL}} \\ \frac{\mathrm{d} [\text{AHL}]}{\mathrm{d} t} &= D \frac{\mathrm{d}^2 [\text{AHL}]}{\mathrm{d} x^2} + P_{\text{AHL}}\end{aligned}\]</span></p>
+
−
+
−
<p>where <span class="math">\( P_{\text{AHL}} \)</span> is the <em>volumic production</em> of <span style="font-variant: small-caps;">AHL</span> in <span class="math">\( \si{mol L^{-1} s^{-1}} \)</span>.</p>
+
−
+
−
<p>Outside (<span class="math">\( x > |w/2| \)</span>) of the layer, where there is diffusion and degradation, we get:</p>
+
−
+
−
<p><span class="math">\[\begin{aligned}
+
−
\mathrm{d} n &= (( \Phi(x) - \Phi(x + \mathrm{d} x) ) \mathrm{d} S - k_{\text{deg}} \mathrm{d} n) \, \mathrm{d} t \\
+
−
\frac{\mathrm{d} n}{\mathrm{d} V \mathrm{d} t} &= D \frac{\mathrm{d}^{2} [\text{AHL}]}{\mathrm{d} x^2} - k_{\text{deg}} \frac{\mathrm{d} n}{\mathrm{d} V} \\
<h2 id="diffusion-inside-the-layer">Solving for diffusion inside the layer</h2>
+
−
+
−
<p>We are interested in the concentration profile of <span style="font-variant: small-caps;">AHL</span> at steady state. Indeed, we assume that diffusion happens faster than the colonization of the bacteria (happening over 2 days <a href="#bib1" class="forward-ref">[1]</a>), so <span class="math">\( P_{\text{AHL}} \)</span> is considered constant in this quasi steady state assumption (QSSA), and we proceed as follows for <span class="math">\( x \in (-w/2, w/2) \)</span>:</p>
[\text{AHL}](x) &= - \frac{P_{\text{AHL}}}{2D} x^2 + c_1 & \text{\( c_0 = 0\) because of symmetry}\end{aligned}\]</span></p>
+
−
+
−
<p>The concentration profile has a parabolic shape.</p>
+
−
</section>
+
−
+
−
<section>
+
−
<h2 id="diffusion-outside-the-layer">Solving for diffusion outside the layer</h2>
+
−
+
−
<p>Applying a similar QSSA, we solve the differential equation as follows for <span class="math">\( x \in (-\infty, -w/2) \cup (w/2, +\infty) \)</span>:</p>
<p>If we assess the order of magnitude of each term, we notice that one can simplify this result since <span class="math">\( \kappa \simeq \SI{7.7}{\milli\metre} \)</span> and <span class="math">\( w/4 \simeq \SI{0.125}{\milli\metre} \)</span>:</p>
<p>The final concentration of <span style="font-variant: small-caps;">AHL</span> is:</p>
+
−
+
−
<p><span class="math">\[[\text{AHL}] =
+
−
\frac{P_{\text{AHL}} w \kappa}{2 D} \cdot
+
−
\begin{cases}
+
−
\exp{(\frac{1}{\kappa} (x + w/2))} & for \space x < -w/2 \\
+
−
\ 1+\frac{1}{w \kappa}\left(\frac{w^2}{4} - x^2 \right) & for \space -w/2 < x < w/2 \\
+
−
\exp{(\frac{1}{\kappa} (w/2 - x))} & for \space x > w/2
+
−
\end{cases}\]</span></p>
+
−
+
−
<p>A dimensional analysis can confirm that [AHL] is a concentration in M, as <span class="math">\(P_{\text{AHL}}\)</span> is a production rate in M.min<sup>-1</sup>, w and <span class="math">\(\kappa\)</span> two lengths in m, and D a diffusion coefficient in m<sup>2</sup>.min<sup>-1</sup> (the remaining expressions are dimensionless).
+
−
</p>
+
−
</section>
+
−
+
−
<section>
+
−
<h2>AHL concentration inside of the layer</h2>
+
−
+
−
<p>To complete our model at the bacterial circuit level, we only need to know the AHL concentration inside the colonization layer. This is the AHL concentration the bacteria will be exposed to and to which they will react. If we compute the concentration of AHL at x=0, which is where the concentration is maximal, we get:</p>
\frac{P_{\text{AHL}} w \kappa}{2 D}(1+\frac{w}{4 \kappa} )
+
−
\end{aligned}\]</span></p>
+
−
+
−
<p>With the numerical values of w and <span class="math">\(\kappa\)</span> applying to our problem, since we have <span class="math">\(\frac{w}{4 \kappa} \simeq 0.02\)</span>, we can neglect this component of the equations, which amounts to ignore the intra-layer variation. The intuition is that this simplification is allowed to us because of the very fast diffusion compared to the width of the colonization layer (high <span class="math">\(\kappa\)</span> and small w), which results in a high homogenization of AHL concentration into the layer.</p>
+
−
+
−
<p>We therefore get the following equation relating AHL concentration to its volumetric production:</p>
+
−
+
−
<p><span class="math">\[[\text{AHL}] =
+
−
\frac{P_{\text{AHL}} w \kappa}{2 D}
+
−
\]</span></p>
+
−
</section>
+
−
+
−
<section>
+
−
<h1>Initial test of our model</h1>
+
−
+
−
<p>Our final in-vivo model comprises the following equations:</p>
+
−
+
−
<p><span class="math">\[\begin{aligned}
+
−
\text{[LuxR]}_0 &= \frac{a_{\text{LuxR}}}{d_{\text{LuxR}}} & \text{LuxR steady state concentration} \\
+
−
[\text{LuxR-AHL}] &= K_{LuxRAHL} [\text{LuxR}]^2 [\text{AHL}]^2 & \text{rapid binding equilibrium} \\
+
−
[\text{LuxR}] &= \frac{a_{\text{LuxR}}}{d_{\text{LuxR}}} - 2 [\text{LuxR-AHL}] & \text{mass conservation}\end{aligned}\]
&= d_{\text{cell}} a_{\text{AHL}} [\text{luxI}] & \text{AHL concentration} \\
+
−
\frac{d_{\text{cell}} a_{\text{AHL}} [\text{luxI}]} w \kappa}{2 D} & \text{AHL concentration} \\
+
−
\]</span></p>
+
−
+
−
<p>To simulate responses of our system under different values of lactate and bacterial cell density input, we solve this system using the fzero function of Matlab. Here is an example of the response of our system with typical values (see tables in the following part) to see how it behaves:</p>
<p>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). We can see that our system behaves well like an AND-gate as expected, but that the levels at which the transitions happen are not the right ones for our application. To find under which circumstances the system behaves as we need it to, we have to proceed to a parameter search.</p>
+
−
</section>
+
−
+
−
<section>
+
−
<h1>Parameter search</h1>
+
−
<p>Even before we got our first parts cloned and characterized, we attempted to predict the requirements that they should meet to achieve the criteria <a href="https://2017.igem.org/Team:ETH_Zurich/Model/Environment_Sensing/system_specifications">previously</a> established from literature data. For this, we extensively explored the space of the parameters controlling our model, and simulated the response of potential systems to find the subsets of parameters combinations satisfying the performances needed to get a sensitive and specific and tumor sensing circuit.</p>
+
−
+
−
<h2>Different categories of parameters</h2>
+
−
<p>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 parameter (typically protein degradation rates, which have an influence opposed to maximal expression rates of the respective protein), and leave other parameters free to vary to check what would be their influence on our system.</p>
<td>Half-activation lactate concentration of the hybrid promoter</td>
<td>Half-activation lactate concentration of the hybrid promoter</td>
−
<td>2.10<sup>6</sup> nM</td>
+
<td>2x10<sup>6</sup> nM</td>
<td><a href="http://parts.igem.org/Part:BBa_K1847008">Characterized lactate sensing part on which our AND-gate is based</a></td>
<td><a href="http://parts.igem.org/Part:BBa_K1847008">Characterized lactate sensing part on which our AND-gate is based</a></td>
−
<td>1.10<sup>4</sup> nM</td>
+
<td>1x10<sup>4</sup> nM</td>
−
<td>1.10<sup>8</sup> nM</td>
+
<td>1x10<sup>8</sup> nM</td>
</tr>
</tr>
<tr>
<tr>
−
<td>k_luxI</td>
+
<td>k<sub>LuxI</sub></td>
<td>Leakiness of the hybrid promoter</td>
<td>Leakiness of the hybrid promoter</td>
<td>0.01</td>
<td>0.01</td>
Line 414:
Line 147:
</tr>
</tr>
<tr>
<tr>
−
<td>K_luxR</td>
+
<td>K<sub>LuxR</sub></td>
<td>Half-activation LuxR-AHL concentration of the hybrid promoter</td>
<td>Half-activation LuxR-AHL concentration of the hybrid promoter</td>
<td>5 nM</td>
<td>5 nM</td>
Line 422:
Line 155:
</tr>
</tr>
<tr>
<tr>
−
<td>n_luxR</td>
+
<td>n<sub>LuxR</sub></td>
<td>Hill coefficient of the hybrid promoter regarding LuxR-AHL concentration</td>
<td>Hill coefficient of the hybrid promoter regarding LuxR-AHL concentration</td>
<td>1.7</td>
<td>1.7</td>
Line 430:
Line 163:
</tr>
</tr>
<tr>
<tr>
−
<td>n_lac</td>
+
<td>n<sub>Lac</sub></td>
<td>Hill coefficient of the hybrid promoter regarding lactate concentration</td>
<td>Hill coefficient of the hybrid promoter regarding lactate concentration</td>
<td>1.7</td>
<td>1.7</td>
Line 438:
Line 171:
</tr>
</tr>
<tr>
<tr>
−
<td>ar_azu</td>
+
<td>k<sub>Azu-LuxI</sub></td>
−
<td>Relative expression of azurin compared to luxR</td>
+
<td>Relative expression of azurin compared to LuxR</td>
<td>10 times the luxI expression</td>
<td>10 times the luxI expression</td>
<td>estimated</td>
<td>estimated</td>
Line 446:
Line 179:
</tr>
</tr>
</table>
</table>
+
</details>
+
<details>
+
<summary>How do we quantitatively define the functional space?</summary>
+
<h2>Cost function</h2>
+
<p>To be able to distinguish systems satisfying the <a href="https://2017.igem.org/Team:ETH_Zurich/Model/Environment_Sensing/system_specifications">criteria</a> 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:</p>
<p>Each argument of the max function represents in the same order the following criteria:</p>
−
<p>To be able to distinguish between systems satisfying or not the <a href="https://2017.igem.org/Team:ETH_Zurich/Model/Environment_Sensing/system_specifications">criteria</a> about specificity and azurin production that we have set for our circuit, we need to use a numerically evaluable condition that will enable our optimization script whether a system is good or not, and "how much good or bad". For this, we will use the following cost function:</p>
+
<ol>
+
<li>Specificity of the sensing for lactate</li>
+
<li>Specificity of the sensing for bacterial cell density</li>
+
<li>Achieving a large amount of produced azurin</li>
<p>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.
−
<p>This function takes the value of the most badly respected criterium, a criterium being respected when the corresponding ratio is below 1. Over 1, a ratio is not good enough, but get better and better as it gets smaller. This monotonicity enables us to rely on optimization algorithms to reach the best combination of parameters available. Also, we can say that every system having a cost function value below 1 is good enough for us, with still the smaller the better.
+
</p>
</p>
+
</details>
+
<details>
+
<summary>Intermediate modeling result: Analysis of the functional parameter space</summary>
<p>Using an optimization toolbox developed for biological systems, MEIGO <a href="#bib6" class="forward-ref">[6]</a>, followed by a package exploring parameter spaces, HYPERSPACE <a href="#bib7" class="forward-ref">[7]</a>, we could obtain the following graphs describing, in the high-dimension space of all possible circuits, a subset of systems satisfying our performance criteria:
<p>Using an optimization toolbox developed for biological systems, MEIGO <a href="#bib6" class="forward-ref">[6]</a>, followed by a package exploring parameter spaces, HYPERSPACE <a href="#bib7" class="forward-ref">[7]</a>, we could obtain the following graphs describing, in the high-dimension space of all possible circuits, a subset of systems satisfying our performance criteria:
</p>
</p>
Line 464:
Line 208:
</figure>
</figure>
</div>
</div>
−
<p>On this figure are drawn the systems suitable for our application. All the axis are logarithmic, except for n_lac and n_lux. 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).
+
<p>On this figure are shown the systems suitable for our application. All the axis are logarithmic, except for <em>n</em><sub>Lac</sub> and <em>n</em><sub>LuxR</sub>. 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).
<ol>
<ol>
−
<li>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</li>
+
<li>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</li>
−
<li>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.</li>
+
<li>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.</li>
<li>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.</li>
<li>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.</li>
<li>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.</li>
<li>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.</li>
−
<li>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 a_luxR/a_luxI combinations (see also n_lac vs a_luxR and n_lac vs a_luxI 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</li>
+
<li>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 <em>a</em><sub>LuxR</sub>/<em>a</em><sub>LuxI</sub> combinations (see also <em>n</em><sub>Lac</sub> vs <em>a</em><sub>LuxR</sub> and <em>n</em><sub>Lac</sub> vs <em>a</em><sub>LuxI</sub> 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</li>
−
<li>K_luxR and n_luxR don't have a significative influence on our system, we can stop studying them</li>
+
<li><em>K</em><sub>LuxR</sub> and <em>n</em><sub>LuxR</sub> don't have a significant influence on our system, we can stop studying them</li>
</ol>
</ol>
</p>
</p>
−
<p> 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 (a_luxR and a_luxI here) as well as a judicious choice of a previously characterized lactate sensor circuit (comprising lldR and lldP genes) among the <a href="https://2015.igem.org/Team:ETH_Zurich/Part_Collection">iGEM ETH 2015 part collection </a>.
+
</details>
+
</div>
+
<h1>Final modeling result: Experimental guidelines used for circuit design</h1>
+
<p> 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 (<em>a</em><sub>LuxR</sub> and <em>a</em><sub>LuxI</sub> here) as well as a judicious choice of a previously characterized lactate sensor circuit (comprising lldR and lldP genes) among the <a href="https://2015.igem.org/Team:ETH_Zurich/Part_Collection">iGEM ETH 2015 part collection </a>.
<p>To translate these insights into experimental results in the lab, we need to chose a target in the range of parameters that work for our application. With the help of the previously characterized initial values for a_luxR (5 nM.min<sup>-1</sup>) and a_luxI (1.10<sup>3</sup> nM.min<sup>-1</sup>), we can hope to tune our system and reach our target in the parameter space via simple RBS tuning, thanks to the Salis calculator.
+
<details>
+
<summary>Target parameters and restrictions</summary>
+
<p>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 <em>a</em><sub>LuxR</sub> (5 nMmin<sup>-1</sup>) and <em>a</em><sub>LuxI</sub> (1x10<sup>3</sup> nMmin<sup>-1</sup>), we can hope to tune our system and reach our target in the parameter space via simple RBS tuning which is supported by the <a href="https://salislab.net/software/forward">Salis Lab RBS Calculator</a> <a href="#bib8" class="forward-ref">[7]</a>.
</p>
</p>
−
<p>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 expresssion for this gene on the Salis calculator, this is why the range a_luxR > 1.10<sup>2</sup> nM.min<sup>-1</sup> is inaccessible to us (grey area), and that we have to chose luxI in consequence. We also get to chose K_lac among the ones available in the promoter collection of parts ranging from BBa_K1847002 to BBa_K1847009: between 0.3 mM and 2.4mM.
+
<p>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 <em>a</em><sub>LuxR</sub> > 1x10<sup>2</sup> nMmin<sup>-1</sup> is inaccessible to us (grey area), and that we have to choose LuxI in consequence. We also get to chose <em>K</em><sub>Lac</sub> among the ones available in the promoter collection of parts ranging from <a href="http://parts.igem.org/Part:BBa_K1847002">BBa_K1847002</a> to <a href="http://parts.igem.org/Part:BBa_K1847009">BBa_K1847009</a>: between 0.3 mM and 2.4mM.
</p>
</p>
−
<p>Taking into account these experimental constraints, 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):
+
</details>
+
</div>
+
+
<p>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):
</p>
</p>
<div>
<div>
Line 488:
Line 240:
alt="Parameter search iteration 2"
alt="Parameter search iteration 2"
/>
/>
+
<figcaption>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.</figcaption>
</figure>
</figure>
</div>
</div>
−
<p>With a_luxR = 1.10<sup>2</sup> nM.min<sup>-1</sup>, a_luxI = 1.10<sup>4</sup> nM.min and K_lac = 1.10<sup>6</sup> 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:
+
−
<ol>
+
<section class="emphasize" id="GuidelineToParts">
−
<li>Use a 10 times stronger RBS than on the <a href="https://2014.igem.org/Team:ETH_Zurich/lab/sequences">piG0047 sequence of iGEM ETH 2014 team</a> for the expression of luxR</li>
+
−
<li>Use a 10 times stronger RBS than on the <a href="https://2014.igem.org/Team:ETH_Zurich/lab/sequences">piG0050 sequence of iGEM ETH 2014 team</a> for the expression of luxI</li>
+
<p>With <em>a</em><sub>LuxR</sub> = 1x10<sup>2</sup> nMmin<sup>-1</sup>, <em>a</em><sub>LuxI</sub> = 1x10<sup>4</sup> nMmin<sup>-1</sup> and <em>K</em><sub>Lac</sub> = 1x10<sup>6</sup> 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:
−
<li>Use the <a href="http://parts.igem.org/Part:BBa_K1847008">BBa_K1847008 part</a> with J23118-B0034 regulatory sequences, giving K_lac = 1.8 mM</li>
+
+
<ol>
+
<li>Use a 10 times stronger RBS than on the <a href="https://2014.igem.org/Team:ETH_Zurich/lab/sequences">piG0047 sequence of iGEM ETH 2014 team</a> for the expression of LuxR</li>
+
<li>Use a 10 times stronger RBS than on the <a href="https://2014.igem.org/Team:ETH_Zurich/lab/sequences">piG0050 sequence of iGEM ETH 2014 team</a> for the expression of LuxI</li>
+
<li>Use the <a href="http://parts.igem.org/Part:BBa_K1847008">BBa_K1847008 part</a> with J23118-B0034 regulatory sequences, giving <em>K</em><sub>lac</sub> = 1.8 mM</li>
</ol>
</ol>
</p>
</p>
−
<p>These value were the basis for the design of our parts and the <a href="https://2017.igem.org/Team:ETH_Zurich/Experiments/Tumor_Sensor">subsequent experimentations</a>. We can validate on our model that they would work well to distinguish the specific levels dictated by our application:
+
<p>These value were the basis for the design of our parts and the <a href="https://2017.igem.org/Team:ETH_Zurich/Experiments/Tumor_Sensor">subsequent experimentations</a>.
+
</section>
+
+
<div class="multi-summary">
+
<details>
+
<summary>Sanity check: in silico behavior for the chosen target parameters</summary>
+
+
We can validate on our model that they would work well to distinguish the specific levels dictated by our application:
</p>
</p>
<div>
<div>
Line 504:
Line 268:
alt="System response after optimization"
alt="System response after optimization"
/>
/>
+
<figcaption>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.</figcaption>
</figure>
</figure>
</div>
</div>
−
<p>We can see on this figure that compared to the initial system using existing parameters, the optimization of crucial and tunable parameters we could perform thanks to the parameter search enables us to tune the circuit so that it detects the right levels of inputs.
+
<p>We can confirm that the obtained parameter set target would lead to a functioning circuit.</p>
+
</details>
+
</div>
+
+
<section>
+
<h1>Going further: a 3D tumor model</h1>
+
<p>To assess the validity of this analytical model, and especially validate the assumptions regarding <a href="https://2017.igem.org/Team:ETH_Zurich/Model/Environment_Sensing/model#diff_model">the AHL diffusion model</a>, we developed a <a href="https://2017.igem.org/wiki/index.php?title=Team:ETH_Zurich/Model/In_Vivo">comprehensive 3D model</a> of a colonized tumor on the software COMSOL, on which we implemented the same equations governing our bacterial circuit.</p>
</section>
</section>
+
<section class="references">
<section class="references">
Line 519:
Line 291:
<li id="bib6"><a href="http://gingproc.iim.csic.es/meigo.html">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. <source>BMC Bioinformatics</source> 15:136.</a></li>
<li id="bib6"><a href="http://gingproc.iim.csic.es/meigo.html">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. <source>BMC Bioinformatics</source> 15:136.</a></li>
<li id="bib7"><a href="https://gitlab.com/csb.ethz/HYPERSPACE">Zamora-Sillero E, Hafner M, Ibig A, Stelling J, Wagner A (2011) Efficient characterization of high-dimensional parameter spaces for systems biology. <source>BMC Syst Biol</source> 5: 142.http://doi.org/10.1186/1752-0509-5-142</a></li>
<li id="bib7"><a href="https://gitlab.com/csb.ethz/HYPERSPACE">Zamora-Sillero E, Hafner M, Ibig A, Stelling J, Wagner A (2011) Efficient characterization of high-dimensional parameter spaces for systems biology. <source>BMC Syst Biol</source> 5: 142.http://doi.org/10.1186/1752-0509-5-142</a></li>
+
<li id="bib8"><a href="https://www.nature.com/nbt/journal/v27/n10/full/nbt.1568.html">Salis, Howard M., Ethan A. Mirsky, and Christopher A. Voigt. "Automated design of synthetic ribosome binding sites to control protein expression."<source>Nature biotechnology</source> 27.10 (2009): 946-950. http://dx.doi.org/10.1038/nbt.1568</a></li>
</ol>
</ol>
</section>
</section>
Latest revision as of 15:07, 23 November 2017
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.
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.
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:
Each argument of the max function represents in the same order the following criteria:
Specificity of the sensing for lactate
Specificity of the sensing for bacterial cell density
Achieving a large amount of produced azurin
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:
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).
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
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.
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.
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.
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
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 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:
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:
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.
