The missing link between AHL production and its concentration: DIFFUSION

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.

Assumption: negligible AHL degradation in the layer

The colony produces AHL 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 \( D = \SI{4.9e-6}{cm^2 s^{-1}} \) (diffusion constant of AHL in water) [4]. We also take into account extracellular degradation of AHL, described by the degradation constant \( k_{\text{deg}} = \SI{5e-4}{min^{-1}} \) [4]. Using the relationship between the mean square distance and time in a Brownian movement, we can estimate the average time needed for AHL to diffuse out of the layer:

\[\Delta t = \frac{w^2}{D} = \frac{0.05^2}{4.9 \cdot 10^{-6}} \simeq \SI{510}{s}\]

During this time \( \Delta t \), we can estimate the magnitude of the degradation of AHL. The proportion of AHL that gets degraded before diffusing out of the layer can be estimated as:

\[\begin{aligned} \Delta t \, k_{\text{deg}} &= \frac{510}{60}\cdot 5 \cdot 10^{-14} \simeq 0.4 \%\end{aligned}\]

Therefore, we can consider for further work that the degradation of AHL happening in the layer is negligible.

Assumption: AHL doesn't diffuse far from where it is produced

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 AHL diffuses before being degraded. Considering the half-life of AHL and the mean distance covered in this given time by diffusion:

\[\begin{aligned} t_{1/2} &= \frac{\ln(2)}{k_{\text{deg}}} \simeq \SI{1400}{min} \\ \\ d &= \sqrt{D t_{1/2}} \simeq \SI{6}{mm}\end{aligned}\]

Since \( d < r_1 = \SI{10}{mm} \), AHL 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.

Model of AHL diffusion

Given the latter assumption, and given \( w \ll r_1 \), we assume that the colonized area appears locally as an infinite sheet of width \( w \). According to the first assumption stated above, we consider only production and diffusion to be significant within the sheet, and diffusion and degradation to be significant outside the sheet. We model the sheet as perpendicular to the \( x \) axis, centered at \( 0 \), covering the interval \( [-w/2, w/2 ] \).

For a small parallelepiped of surface \( \mathrm{d} S \) and width \( \mathrm{d} x \) perpendicular to \( x \) axis (see figure below), we can use Fick’s law to model diffusion of AHL.

The flux of the protein is:

\[\Phi(x) = - D \left. \frac{\mathrm{d} [\text{AHL}]}{\mathrm{d} x} \right|_{x}\]

When inside (\( x < |w/2| \)) of the layer, the change of protein amount within time length \( \mathrm{d} t \) is equal to diffusive transports plus production:

\[\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}\]

where \( P_{\text{AHL}} \) is the volumic production of AHL in \( \si{mol L^{-1} s^{-1}} \).

Outside (\( x > |w/2| \)) of the layer, where there is diffusion and degradation, we get:

\[\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} \\ \frac{\mathrm{d} [\text{AHL}]}{\mathrm{d} t} &= D \frac{\mathrm{d}^{2} [\text{AHL}]}{\mathrm{d} x^2} - k_{\text{deg}} [\text{AHL}]\end{aligned}\]

Solving for diffusion inside the layer

We are interested in the concentration profile of AHL at steady state. Indeed, we assume that diffusion happens faster than the colonization of the bacteria (happening over 2 days [1]), so \( P_{\text{AHL}} \) is considered constant in this quasi steady state assumption (QSSA), and we proceed as follows for \( x \in (-w/2, w/2) \):

\[\begin{aligned} \frac{\mathrm{d} [\text{AHL}]}{\mathrm{d} t} &= 0 \\ \frac{\mathrm{d}^2 [\text{AHL}]}{\mathrm{d} x^2} &= - \frac{P_{\text{AHL}}}{D} \\ [\text{AHL}](x) &= - \frac{P_{\text{AHL}}}{2D} x^2 + c_0 x + c_1 \\ [\text{AHL}](x) &= - \frac{P_{\text{AHL}}}{2D} x^2 + c_1 & \text{\( c_0 = 0\) because of symmetry}\end{aligned}\]

The concentration profile has a parabolic shape.

Solving for diffusion outside the layer

Applying a similar QSSA, we solve the differential equation as follows for \( x \in (-\infty, -w/2) \cup (w/2, +\infty) \):

\[\begin{aligned} \frac{\mathrm{d} [\text{AHL}]}{\mathrm{d} t} &= 0 \\ \frac{\mathrm{d}^2 [\text{AHL}]}{\mathrm{d} x^2} - \frac{k_{\text{deg}}}{D} [\text{AHL}] &= 0 \\ \frac{\mathrm{d}^2 [\text{AHL}]}{\mathrm{d} x^2} - 1/\kappa^2 [\text{AHL}] &= 0 & \kappa = \sqrt{D/k_{\text{deg}}} \simeq \SI{7.7}{\milli\metre} \space \text{is the characteristic length of diffusion}\\ [\text{AHL}](x) &= c_2 \exp{(-x/\kappa)} + c_3 \exp{(x/\kappa)} + c_4 \end{aligned}\]

We have \( c_2 = c_3 := c \) due to symmetry of the problem.

Boundary conditions

AHL concentration is \( 0 \) at \( \pm \infty \). This implies \( c_4 = 0 \) and:

\[\begin{aligned} [\text{AHL}] = \begin{cases} c \exp{(x/\kappa)} & for \space x < -w/2 \\ c \exp{(-x/\kappa)} & for \space x > w/2 \\ \end{cases}\end{aligned}\]

The flux of AHL should be continuous at the boundary of the colonized area:

\[\begin{aligned} \frac{\mathrm{d} [\text{AHL}]}{\mathrm{d} x}(w^{+}/2) &= \frac{\mathrm{d} [\text{AHL}]}{\mathrm{d} x}(w^{-}/2) \\ -(1/\kappa) c \exp{(- w/2\kappa)} &= -\frac{P_{\text{AHL}}}{D} \frac{w}{2} \\ c &= \frac{P_{\text{AHL}}w}{2D\kappa} \exp{(w/2\kappa)}\end{aligned}\]

Plus, due to continuity of the concentration of AHL:

\[\begin{aligned} [\text{AHL}](w^{+}/2) &= [\text{AHL}](w^{-}/2) \\ c \exp{(-\kappa w/2)} &= -\frac{P_{\text{AHL}}}{2D} \frac{w^2}{4} + c_1 \\ c_1 &= c \exp{(-w/2\kappa)} + \frac{P_{\text{AHL}}}{2D} \frac{w^2}{4} \\ c_1 &= \frac{P_{\text{AHL}}w\kappa}{2D} + \frac{P_{\text{AHL}}w^2}{8D} \\ c_1 &= \frac{P_{\text{AHL}}w}{2D}(\kappa + w/4)\end{aligned}\]

If we assess the order of magnitude of each term, we notice that one can simplify this result since \( \kappa \simeq \SI{7.7}{\milli\metre} \) and \( w/4 \simeq \SI{0.125}{\milli\metre} \):

\[c_1 \simeq \frac{P_{\text{AHL}}w\kappa}{2D}\]

Full solution

The final concentration of AHL is:

\[[\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}\]

A dimensional analysis can confirm that [AHL] is a concentration in M, as \(P_{\text{AHL}}\) is a production rate in M.min^-1, w and \(\kappa\) two lengths in m, and D a diffusion coefficient in m².min^-1 (the remaining expressions are dimensionless).

AHL concentration inside of the layer

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:

\[\begin{aligned}\text{AHL}(x=0) = \frac{P_{\text{AHL}} w \kappa}{2 D}(1+\frac{w}{4 \kappa} ) \end{aligned}\]

With the numerical values of w and \(\kappa\) applying to our problem, since we have \(\frac{w}{4 \kappa} \simeq 0.02\), 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 \(\kappa\) and small w), which results in a high homogenization of AHL concentration into the layer.

We therefore get the following equation relating AHL concentration to its volumetric production:

\[[\text{AHL}] = \frac{P_{\text{AHL}} w \kappa}{2 D} \]

Initial test of our model

Our final in-vivo model comprises the following equations:

\[\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}\] [\text{luxI}] &= \frac{a_{\text{luxI}}}{d_{\text{luxI}}} (k_{\text{luxI}} + (1 - k_{\text{luxI}}) P_{\text{Lux-Lac}})\end{aligned}\] & \text{LuxI concentration} \\ P_{\text{AHL}} &= 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} \\ \]

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:

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.

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.

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

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

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

Parameter search

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:

\[\max\left(\frac{10\times azu(low\space lac,HIGH\space d_{cell})}{azu(HIGH\space lac,HIGH\space d_{cell})},\frac{10\times azu(HIGH\space lac,low\space d_{cell})}{azu(HIGH\space lac,HIGH\space d_{cell})},\frac{1.10^{6}}{azu(HIGH\space lac,HIGH\space d_{cell})}\right)\]

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.

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 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).

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 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
6. K_luxR and n_luxR don't have a significant influence on our system, we can stop studying them

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 iGEM ETH 2015 part collection .

Our target for parameters

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 a_luxR (5 nM.min^-1) and a_luxI (1.10³ nM.min^-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 expresssion for this gene on the Salis calculator, this is why the range a_luxR > 1.10² nM.min^-1 is inaccessible to us (grey area), and that we have to choose 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.

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):

With a_luxR = 1.10² nM.min^-1, a_luxI = 1.10⁴ nM.min and K_lac = 1.10⁶ 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 K_lac = 1.8 mM

These value were the basis for the design of our parts and the subsequent experimentations. We can validate on our model that they would work well to distinguish the specific levels dictated by our application:

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.