Team:Wageningen UR/Model/QS

Quorum sensing and cell lysis

We developed a mathematical model of the signal amplification module of the Mantis system. Using quorum sensing (QS) and a lytic mechanism, the goal of this module is to take a weak antigen signal detected only by a minority of the biosensor cells, and to turn this into a population-wide response. Thus, a signal can be detected despite a low concentration of antigens in the tested sample. We explored the role of various reaction rates on the behavior of the system. Our mathematical model shows that the QS system can spontaneously activate, be insensitive to N-Acyl homoserine lactone (AHL), or function properly, depending on the reaction rates of the system. In our most complete model, we can change a spontaneously activating system into a system with proper signaling behavior by either increasing LuxR degradation, or by adding the enzyme aiiA to the genetic circuit.

Alongside the development of this model, the amplification module described here was built in the lab. We described how the modeling and wet-lab projects were integrated on a separate page. The QS system depends on a detection module that is able to detect disease antigens.

Introduction

The Mantis biosensor can use an optional amplification module. A detection module can be directly linked to a reporter module. When using the QS system, detection of an antigen will trigger the production of the QS signaling molecule AHL. Only when AHL levels are sufficiently high, cells will activate their reporter module. A general overview of the QS and lysis module is given in Figure 1.

Figure 1: A general overview of the QS module of Mantis. Some cells will detect antigen and start producing AHL. Other cells will detect AHL and the entire population will start QS, triggering cell lysis. The release of cellular components upon lysis allows the two halves of the split fluorescent protein to generate a fluorescent signal. In green dots an entire fluorescent molecule (GFP) is depicted.

Adding an amplification module to Mantis can make it more sensitive to very low levels of antigen, making it able to diagnose diseases that are difficult to detect. As the QS system depends on transcription and translation, the QS system will make the diagnostic test slower. But different circumstances require diagnostics that are either very rapid, but not very sensitive, or are very sensitive, but take more time until a result is produced. Therefore, the amplification module can be very useful to whole-cell biosensors.

The signaling molecule AHL is central to QS. Together with the two proteins LuxR and LuxI, the system allows bacteria to produce a signaling molecule and activate genes when the bacteria reach a high enough cell density. AHL can freely diffuse across the cell membrane. Inside the cell, AHL can bind with LuxR to form the the LuxR-AHL complex. This complex can dimerize to become an active transcription factor (TF). The LuxR-AHL dimer can result in positive and negative regulation of genes. The pLuxR promoter is positively regulated while pLuxL is negatively regulated by the LuxR-AHL dimer. In our Mantis module (Figure 2), QS is coupled to a lytic mechanism. Initially, there will be two bacterial populations. Each population contains one part of split fluorescent protein (FP). Together, proteins from each population can generate a fluorescent signal. The cellular membranes will prevent the proteins from interacting, so initially there will be no fluorescence.

Figure 2: The molecular components of the Mantis QS and lytic system. Antigen will activate QS (top left). QS will initiate cell lysis (top right). Cell lysis will result in a fluorescent signal (bottom).

Once some cells detect antigen, they will start producing AHL. This will spread to other cells within the vicinity of those cells that detect antigen, eventually causing all cells in Mantis to start QS. The active form of LuxR will activate the lytic mechanism, resulting in lysis of all cells engaging in QS. Once cell contents are lysed, the two complementary parts of the split FP will start to interact. This will result in a fluorescent signal, which will give a read-out on the Mantis device.

In the QS system, LuxI positively regulates its own expression. This can lead to positive feedback, where more and more AHL is produced. Cells with this positive feedback motif can therefore become activated. This is important, as the lytic mechanism should either be active or completely inactive.

To achieve this goal, a temporal model was first developed. The temporal model was used to sample the global parameter space. The second model that we used is the spatiotemporal model. This model was used to test hypotheses and to guide wet-lab experiments.

For the system to function properly, it needs the following properties.

  • Bistability: Upon activation by an external signal, cells must be able to switch from an OFF-state to an ON-state. This requires that the time-dependent system of ordinary differential equations (ODEs) has more than one equilibrium state solution. The bistability of switches based on the lux operon have been explored in the literature. A motif with two feedback loops has been shown to have bistability [1,2]. However, it remains unclear if the native lux box (BBa_K546000) has the ability to switch from a clear OFF-state to an ON-state [3].
  • Cell signaling: Cells in the ON-state should be able to signal cells in the OFF-state to also transition to the ON-state. This is a core property of the QS system [4].
  • Delayed cell lysis: Cells in the ON-state should lyse to release their split FP. But they should do so only after signaling cells still in the OFF-state.

Activator and receiver cells

In the Mantis system, antigen concentrations will be low and therefore only a small number of cells come into contact with antigens. For the temporal model, we assumed that 10% of all cells would come into direct contact with the antigen, allowing them to detect the antigen. These cells are the activator cells. In the model, detection of antigen is simulated by switching on LuxI and LuxR production in 10% of the cells. The receiver cells are identical to the activator cells in all respects, except for this activation of LuxR and LuxI production.

Biological components

In the activator cells, both LuxI and LuxR will be produced under control of a promoter pInd induced by antigen. LuxI can also be produced from a promoter pLuxA, which is positively regulated by LuxR-AHL. LuxR is either positively or negatively regulated by LuxR-AHL and produced from pLuxB. The lysis protein is produced from pLuxC and is always positively regulated by LuxR-AHL.


Temporal model

The main goal of the temporal model is to assess the dynamics of the system and what kind of responses it has to the antigen input. To effectively do this, the biological picture has to be translated to mathematical equations. Furthermore, this model is more manageable and more straightforward to interpret when the number of kinetic parameters are low. To achieve this, the following simplifying assumptions have been made.

System

The system consists of three compartments; the external medium, antigen-activated bacterial cells, and signal-receiving bacterial cells (Figure A). The equations of both cells are identical, except that the activator cells receive an external signal intended to put the cells in the ON-state at t = 240. This simulates the addition of antigen-containing blood to the bacteria and detection of this low concentration of antigen by a minority of the cells. The total size of the system is 1 000 000 volume units. The ratio of external volume to cells is 1 cell in 1000 volume units. This approximates a cell density of 1 OD600. Therefore, the total cell volume is 1000 volume units. As the volume unit is defined as the volume of a single cell, there are a total of 1000 cells in the simulated system. Out of these cells, 900 cells are receiver cells, 100 are activator cells. Both of them are represented by single cell compartments. When the compartment representing each cell exchanges species with the external medium compartment, the total amount of exchange is calculated by using the total for each cell type. This exchange of AHL allows the cells to signal across the medium, communicating the presence of antigen. This leads to cell lysis and activation of fluorescence.

Figure A: The temporal model has three different compartments, the external medium, activator cells, and receiver cells. Diffusion is assumed to be instantaneous inside each compartment allowing us to neglect spatiotemporal dynamics. Diffusion across the cell membrane is computed using a linear approximation, as indicated by the 'flux' arrows. ODEs describe the rate of change of all species inside the three compartments.

Each cell has concentration values for the following species: AHL, LuxI, LuxR, LuxR-AHL, split FP and lysis protein. In the external medium, there are concentration values for AHL, C-terminal split FP, N-terminal split FP, and the fused and active FP. All compartments are assumed to be well-mixed, having a homogeneous concentration. All cells are identical in size and their size is constant throughout the simulation. There is assumed to be no cell growth and therefore no dilution term.

Kinetics

The levels of mRNA of all genes expressed is assumed to be in a quasi-steady state. The basis for this assumption is that translation is slower than transcription. Therefore, the concentration of mRNA species in the model is assumed to be constant. Similarly, transcription factors binding to DNA are also assumed to occur quickly since binding and unbinding of a transcription factor to a promoter occurs at timescales much faster than translation. The third assumption is that the dimerization of LuxR-AHL takes place quicker than protein dynamics. Therefore, the concentration of the LuxR-AHL dimer is assumed constant in the model. To summarize, transcription factor-DNA complexes, mRNA, and LuxR-AHL dimer concentrations are implicitly modeled. All promoters are assumed to have a single transcription binding site and the transcription factor is assumed to form a dimer only, thus we use a Hill coefficient of 2.

The lysis mechanism is not explicitly modeled. Lysis is assumed to occur instantaneously the moment the lysis protein concentration exceeds a value of 1. Before the lysis protein reaches this concentration, the lysis mechanism is inert. Diffusion of AHL across the cell membrane is modeled by a linear diffusion equation depending on the concentration difference inside and outside of the cell, and a diffusion constant for diffusion of AHL across the cell membrane. This diffusion constant is assumed to be the same whether AHL diffuses into or out of the cell. As stated earlier, the cells as well as the external medium are considered to be well-mixed, which equates to assuming that both the diffusion of AHL inside each compartment is instantaneous. Similarly, the release of split FPs upon cell lysis is assumed to be instantaneous.

Protein production in the model is dependent on the expression rate of the gene. Production of AHL is dependent on LuxI. Except for the FPs in the external medium, all species also have degradation rates. The split and fused FPs are considered to be sufficiently stable so that degradation outside the cell plays no role on the timescales of our simulations. Upon lysis of the cell, only AHL and the split FPs are released into the external medium. All other model components are assumed to be destroyed.


The change in time of every species in the system is described by the following ordinary differential equations.

Ordinary differential equations

${\mathrm{d}AHL_{cell} \over \mathrm{d}t} = \alpha_{1}LuxI + D_{AHL}(AHL_{ext}-AHL_{cell}) + k_{-1}[LuxR-AHL] - k_1 LuxR \cdot AHL -\beta_{1}AHL_{cell}$

${\mathrm{d}AHL_{ext} \over \mathrm{d}t} = V_{ratio} \cdot D_{AHL} (AHL_{ext}-AHL_{cell}) - \beta_{2}AHL_{cell}$

${\mathrm{d}[LuxR-AHL] \over \mathrm{d}t} = k_{1}LuxR \cdot AHL - \beta_{3}[LuxR-AHL]^{2} - \beta_{4}[LuxR-AHL]$

${\mathrm{d}LuxI \over \mathrm{d}t} = \alpha_2 + \alpha_3 {{LuxR \cdot [LuxR-AHL]^2} \over {km_1 + [LuxR-AHL]^2}} + \alpha_4 + S \cdot \alpha_5 - \beta_5 LuxI$


Positive regulation of LuxR:
${\mathrm{d}LuxR \over \mathrm{d}t} = k_{-1}[LuxR-AHL] + \alpha_6 + \alpha_7 {{LuxR \cdot [LuxR-AHL]^2} \over {km_1 + [LuxR-AHL]^2}} + \alpha_8 + S \cdot \alpha_9 - \beta_6 LuxR$


Negative regulation of LuxR:
${\mathrm{d}LuxR \over \mathrm{d}t} = k_{-1}[LuxR-AHL] + \alpha_6 + \alpha_7 {{LuxR \cdot km_1} \over {km_1 + [LuxR-AHL]^2}} + \alpha_8 + S \cdot \alpha_9 - \beta_6 LuxR$

${\mathrm{d}Lysis \over \mathrm{d}t} = \alpha_{10} + \alpha_{11} {{LuxR \cdot [LuxR-AHL]^2} \over {km_2 + [LuxR-AHL]^2}} - \beta_7 Lysis$

${\mathrm{d}SplitFP_{cell} \over \mathrm{d}t} = \alpha_{12} - \beta_8 SplitFP$

${\mathrm{d}SplitFP_{N} \over \mathrm{d}t} = -k_3 SplitFP_{N} \cdot dSplitFP_{C}$

${\mathrm{d}SplitFP_{C} \over \mathrm{d}t} = -k_3 SplitFP_{N} \cdot dSplitFP_{C}$

${\mathrm{d}ActiveFP \over \mathrm{d}t} = k_3 SplitFP_{N} \cdot dSplitFP_{C}$


Species
$AHL_{cell}$ Signaling molecule AHL inside the cell
$AHL_{ext}$ Signaling molecule AHL inside the external medium
$LuxR$ Free inactive LuxR protein
$[LuxR-AHL]$ The LuxR-AHL complex
$LuxI$ The AHL-producing enzyme LuxI
$LuxR$ The AHL-producing enzyme LuxI
$Lysis$ The lysis protein responsible for the lytic mechanism
$SplitFP_{cell}$ The split fluorescent protein in the cell
$SplitFP_{N}$ The N-terminal part of the split fluorescent protein in the external medium
$SplitFP_{C}$ The C-terminal part of the split fluorescent protein in the external medium
$ActiveFP$ The active full form of the split fluorescent protein in the external medium
Parameters
$S$ is 0 or 1 depending on if the cell detected antigen
$\alpha_1$ AHL production rate of LuxI
$\alpha_2$ Leaky production rate of LuxI from pIndnst
$\alpha_3$ Maximum production rate of LuxI from pInd
$\alpha_4$ Leaky production rate of LuxI from pLuxA
$\alpha_5$ Maximum production rate of LuxI from pLuxA
$\alpha_6$ Leaky production rate of LuxR from pInd
$\alpha_7$ Maximum production rate of LuxR from pInd
$\alpha_8$ Leaky production rate of LuxR from pLuxB
$\alpha_9$ Maximum production rate of LuxR from pLuxB
$\alpha_{10}$ Leaky production rate of Lysis from pLuxC
$\alpha_{11}$ Maximum production rate of Lysis from pLuxC
$\alpha_{12}$ Constitutive production of $SplitFP_{cell}$
$k_1$ formation rate of LuxR-AHL
$k_{-1}$ dissociation rate of LuxR-AHL
$k_{2}$ formation & maturation rate of the full fluorescent protein
$km_{1}$ Michaelis-Menten constant of pLuxA and pLuxB for $LuxR-AHL$
$km_{2}$ Michaelis-Menten constant of pLuxC for $LuxR-AHL$
$D_{AHL}$ diffusion rate of AHL across the cellular membrane
$\beta_1$ degradation rate of $AHL_{cell}$
$\beta_2$ degradation rate of $AHL_{ext}$
$\beta_3$ degradation rate of $LuxR-AHL$ complex
$\beta_4$ degradation rate of the $LuxR-AHL$ dimer
$\beta_5$ degradation rate of $LuxI$
$\beta_6$ degradation rate of $LuxR$
$\beta_7$ degradation rate of $Lysis$
$\beta_8$ degradation rate of $SplitFP_{cell}$

The initial concentration for all species was zero. The system of ODEs was solved using the scipy.integrate.odeint module available in Python. To sample the global parameter space, 1 000 000 parameter sets were generated using the lhs() function from pyDOE. The parameters were uniformly distributed between 0 and 10. These parameter sets were simulated for a total of 480 time units, where between t=240 and t=300, the external signal S was set from 0 to 1. A value for 1 was chosen to represent enough lysis protein for cell lysis to occur. Every 0.01 timesteps, the value of lysis protein in each cell was checked. In the case that a cell has a lysis protein value above 1, the cell would be removed from the simulation and the AHL and split FP contents would be released into the external medium. Half of the cells contained the C-terminal FP, whilst the other half of the cells contained the N-terminal version.

Scoring function

To assess how well a parameter set of reaction rates performed, we needed to generate scoring functions that described our ideal system behavior. We want our system to produce a lot of fluorescence quickly, but only after detection of antigen. This can be quantified by calculating the area under the curve, which we approximate using the numpy.trapz() function. As an ideal parameter set generates no fluorescence before t=240 and very quickly generates a lot of fluorescence after t=240, the area under the curve before t=240 should be minimal, while after t=240, the area should be as large as possible. The area under the curve represents both how much, as well as how quickly, fluorescence is generated. Secondly, in absence of an external input, (when S = 0 during the entire simulation), no fluorescence should be generated. Fluorescence in the absence of an antigen would amount to a false positive. Therefore, each parameter set is simulated twice. Once with no external signal, and once with an external signal between t=240 and t=300. The final score of a parameter set is calculated by combining both simulation runs. We define this as

$Score_{set} = {w_1 \over { |Score_{noAntigen} - Score_{Antigen}|+ w_2} } + Score_{Antigen}$

$Score_{noAntigen}$ is the $Score_{simulation}$ computed in the absence of antigen

$Score_{Antigen}$ is the $Score_{simulation}$ computed in the presence of antigen

$w_1$ sets an upper bound to the score in case the simulations produce identical levels of fluorescence and was set to $w_1 = 1000$

$w_2$ determines the weight of the difference in fluorescent production between the negative and positive simulation and was set to $w_2 = 1$

$Score_{noAntigen}$ is the $Score_{set}$ for the simulation in absence of antigen

$Score_{Antigen}$ is the $Score_{set}$ for the simulation in the presence of antigen

$Score_{simulation} = \displaystyle\int_{t=0}^{t=240} ActiveFP \mathrm{d}t + {1 \over {w_3 + \displaystyle\int_{t=240}^{t=300} ActiveFP \mathrm{d}t}}$

$ActiveFP$ is the fluorescence at time t

$w_3$ determines the weighting of fluorescent production before and after t=240, and was set to $w_3 = 0.01$

Results

After initial investigation of the model, the following types of system behaviors can be categorized:

  • No lysis: None of the cells in the system lyse. Even the activator cells, that are directly induced by the external signal, are unable to produce enough lysis protein within the duration of the simulation.
  • Activator lysis only: Activator cells are able to lyse, release their FPs and generate some fluorescence. But the majority of the cell population remain unaffected. There is no cell-to-cell signaling through QS.
  • Spontaneous lysis: All cells lyse within the duration of the simulation, but the lysis is not caused by induction of the activator cells by the external signal. Spontaneous auto-activation and cell lysis are caused by the build-up of AHL over time. This build-up is caused by leaky expression of LuxR and LuxI and, in the case of a system overly sensitive to AHL, this leaky expression will lead to auto-activation. Receiver cells will thus auto-activate using AHL they themselves produce, in the absence of antigen.
  • Induced full lysis: Activator cells produce AHL and this triggers the receiver cells to auto-activate as well. At the end of the simulation, all cells will have lysed and released their split FPs, resulting in a significant fluorescent signal. The fluorescent signal is caused by the presence of antigen and not auto-activation. The increase of fluorescence over time of a simulation with the desired behavior is displayed in Figure 3.
Figure 3: The increase of fluorescence for a simulation with induced full lysis. At t=240, 10% of the cells detect antigen and they lyse at t=242, releasing split FP. At t=247, the receiver cells lyse, releasing much more split FP. Click for parameters.


Global parameter sampling

A total of 1 million parameter sets were sampled for the system with positively regulated LuxR. The parameter sets were scored with the scoring functions and sorted by score. In the case of negative regulation of LuxR by the LuxR-AHL dimer, the percentage of parameter sets for each behavior is displayed in Table 1.

Table 1: Results of simulating 1 million parameter sets using Latin hypercube sampling. LuxR was positively regulated by the dimer of the LuxR-AHL complex
Behavior # parameter sets Percentage
No lysis 789430 78.9%
Activator lysis 199560 19.9%
Spontaneous lysis 10986 1.1%
Induced full lysis 21 0.002%

The parameter space for the system with negative LuxR regulation was explored in an identical manner, and the results are displayed in Table 2.

Table 2: Results of simulating 1 million parameter sets using Latin hypercube sampling. LuxR was negatively regulated by the dimer of the LuxR-AHL complex
Behavior # parameter sets Percentage
No lysis 813045 81.3%
Activator lysis 146105 14.6%
Spontaneous lysis 40842 4.1%
Induced full lysis 8 0.001%

The variation in parameters for each behavior was visualized by generating a histogram showing how many parameter sets fell within a certain parameter range for each behavior (positive regulation of LuxR shown in Figure 4, negative regulation regulation of LuxR in Figure 5). This indicates that parameter sets that have low degradation values for LuxR, and all species that contain LuxR, are more likely to display spontaneous lysis. This trend is stronger in the case of positive regulation of LuxR, where the active TF is needed for high levels of LuxR expression. For negative LuxR regulation, basal levels of LuxR are high and the system can spontaneously activate more easily with higher LuxR degradation factors. The parameter distributions for the degradation of LuxI differ significantly between the cases of positive and negative regulation of LuxR by the LuxR-AHL dimer. In the case of positive LuxR regulation, levels of LuxR limit the amount of active TF. When LuxR is negatively regulated, LuxR levels are higher and the levels of AHL are stronger correlated with the total amount of active TF.

Figure 4: The parameter distribution of the LuxR degradation rate for spontaneous lysis (blue) and no lysis (green) is compared. When there is no relationship between a parameter and the behavior of the system, the distribution is uniform. LuxR was positively regulated.

A similar trend is seen in the data of the negative regulation of LuxR.

Figure 5: The parameter distribution of the LuxR degradation rate for spontaneous lysis (blue) and no lysis (green) is compared. When there is no relationship between a parameter and the behavior of the system, the distribution is uniform. LuxR was negatively regulated.


Local parameter space

To test the robustness of parameter sets that produced the system responses we desired, we explored the local parameter space around these optimal parameter sets. By scoring the local parameter sets, set with induced full lysis could be compared to those that lysed spontaneously. While the parameter distributions for many parameters when comparing spontaneous lysis and induced full lysis were still quite similar, a different behavior could be seen for the extracellular degradation of AHL and the degradation rate of LuxR. The results for the external degradation rate of AHL and the degradation rate of LuxR in the situation that LuxR is positive regulated are displayed in Figures 6 and 7.

Figure 6: Parameter sets with spontaneous self-activation (green) have higher external AHL degradation rates. Parameter sets with the desired behavior (blue) have lower external AHL degradation rates.
Figure 7: Parameter sets with induced full lysis (blue) and compared with spontaneous lysis (green). The parameter sets with high values of LuxR degradation are very likely to have the desired behavior. The parameter distribution of induced full lysis appears shifted to the right compared to the distribution of spontaneous lysis.

The same analysis was carried out on the parameters obtained through Latin hypercube sampling for the situation with negative regulation of LuxR (Figure 8 and 9).

Figure 8: Induced full lysis (blue) compared with spontaneous lysis (green). Several high blue peaks are visible for low values of external AHL degradation rates. Of the parameter sets simulated that had high external AHL degradation rates, very few of them had the desired behavior.
Figure 9:Induced full lysis (blue) compared with spontaneous lysis (green). The parameter sets induced full lysis have a clear shift towards higher rates of LuxR degradation when compared to the spontaneous lysis parameter sets.

Both the external AHL degradation rates and LuxR degradation rates seemed to play a large role in determining if a system would spontaneously lyse, or have induced full lysis, the effect of changing these parameters was explored further. We took one spontaneously activating parameter set and in this set varied both these two parameter values from 0 to 10. The scores for these 10 000 (100 by 100) parameter sets can be visualized in a heatmap (Figure 10). This heatmap shows the previously described types of behavior and how the desired behavior is the transition area between spontaneous lysis and lysis only of the activator cells and that by reducing the sensitivity of the system to AHL, for example by increasing the degradation rate of LuxR. This heatmap also shows that the external AHL degradation controls the thickness of the transition area and that for high values the transition area of desired system behavior disappears.

Figure 10: Heatmap of the scoring function when varying both the external AHL degradation rate and the LuxR degradation rate. This shows the transition area of desired behavior in between spontaneous lysis and lysis of only the activator cells.

Summary

The global and local parameter sensitivity of our system was explored. When parameters allow for QS to occur, spontaneous lysis is the most likely outcome. But parameters were identified that distinguish between induced full lysis and spontaneous lysis. The external AHL degradation rate and the LuxR degradation rate have a large effect on system behavior. By decreasing LuxR stability or increasing external AHL stability, a simulation outcome can change from spontaneous lysis to induced full lysis. LuxR degradation can be changed experimentally by adding a degradation tag to a protein. The external AHL degradation cannot be changed. However, it's relative value can be decreased by increasing the cellular degradation rate. Therefore, these results provide sufficient justification for exploring these strategies in the lab.



Spatiotemporal model

As the assumption that AHL diffuses instantaneously breaks down when the system becomes larger, or when the number of activator cells decreases, we developed a spatiotemporal model that incorporates both diffusion of AHL, but also the mixing of the split FPs. To more realistically model the biological system, we used the spatiotemporal model with the following features:

  • Diffusion of the signaling molecule, AHL
  • Diffusion of fluorescent proteins
  • Cell movement
  • Asymmetric cell division
  • Variable growth rates

In the temporal model the cell density had remained constant throughout the simulation. Now, in the spatiotemporal model, cell density will increases exponentially.

This spatiotemporal model is comprised out of a 2D grid where each square is the size of a single bacterial cell. Squares are either external medium, or completely occupied by a single cell. The grid system discretizes space while at the same time gives every cell a definite position.

The finite volume method approximates component concentrations over a spatial domain by dividing the space into smaller volume elements. At every coordinate in space, the concentration value there is set equal to the average concentration over that small volume element. By computing the flux through the faces of each volume piece, the concentration change over time can then be approximated. Diffusion rates from a cell to a medium square or from one medium square to another medium square have different diffusion coefficients, where the diffusion coefficient of AHL inside the external medium is much larger. When two cells are in adjacent squares, diffusion is still possible, but at a much slower rate.

Cell movement in this system is modeled by a random movement where cells have a certain probability to move to the four adjacent squares as well as the four diagonal squares. The model therefore does not realistically model chemotaxis or other forms of bacterial movement [5]. When a cell tries to move into a square that is already occupied, movement is prevented and the cell simply stays in their current square. When a cell moves into an unoccupied square, the AHL concentration of the original square and the destination square get switched around, as squares occupied by cells have no external medium AHL values. Cells don't change in size in this model. A bacteria is always the size of exactly one square. All cells do have their own growth rate. Once the cell has become big enough to divide, the simulation simply creates a new cell in an adjacent square and distributes all cellular species among both daughter cells. Every species is independently and asymmetrically distributed by drawing a number from a normal distribution (µ=0.5, σ=0.2). Once the cells are divided, the new cell is awarded a growth rate independent of that of the parent cell and both cells have their growth rate clock reset to zero. Cells that lyse still release their AHL and FP contents the same way as in the temporal model, except in the spatiotemporal method the square the lysing cell occupies becomes a external medium square with the same concentration of extracellular AHL as had been the cellular AHL concentration of that cell upon lysis.

Most equations for the spatiotemporal model are identical to those in the temporal model. The exceptions are described here. In the spatial model, production of both N- and C-terminal split FP is now explicitly modeled and the following two equations replace the one from the temporal model.

${\mathrm{d}SplitFPC_{cell} \over \mathrm{d}t} = \alpha_{12} - \beta_8 SplitFPC$

${\mathrm{d}SplitFPN_{cell} \over \mathrm{d}t} = \alpha_{12} - \beta_8 SplitFPN$

The change in AHL and split FP in a square as a result of diffusion is calculated using the following equation:

$G_{i,j} = p_{i-1,j} + p_{i,j-1} + p_{i+1,j} + p_{i,j+1} - 4 \cdot p_{i,j}$
$G$ is a scalar value indicating the size and the sign of the flux gradient in the square with coordinates $i$ and $j$
$p_{i,j}$ is the concentration in square the with coordinates $i$ and $j$

$G_{i,j}$ is positive when the concentration in a square is increasing. In that case, more flux enters the square than leaves. The flux out of a square is proportional to the concentration inside that square. As there are four faces through which flux can leave, this term is present four times. Flux can also flow into a square from either of the four neighbors. These represent the four positive terms. Together, this provides a scalar quantity that describes the intensity of the flux.

The change in concentration is calculated by multiplying $U$ with the diffusion constant $D$.

$f = D \cdot G \Delta t$
$f$ is the change concentration in a square
$D$ is the diffusion constant.
$G$ is the gradient

With the addition of $aiiA$, the equation for $AHL_{cell}$ gains a degradation term:

${\mathrm{d}AHL_{cell} \over \mathrm{d}t} = \alpha_{1}LuxI + D_{AHL}(AHL_{ext}-AHL_{cell}) + k_{-1}LuxR-AHL - k_1 LuxR \cdot AHL - k_3 aiiA -\beta_{1}AHL_{cell}$

And the equation describing $aiiA$ itself:

${\mathrm{d}aiiA \over \mathrm{d}t} = (1-S) \cdot \alpha_{13} - \beta_9 aiiA$

Simulations start with an initial value of zero for all species. Again, odeint is used to solve the system of equations. The algorithm checks for cell lysis, movement and division every 0.01 time steps. The duration of a simulation is now 120 time units, with activators detecting antigen at t=60. The diffusion constant of AHL in the external medium is 100, unless otherwise stated. The diffusion constant of split FP inside the external medium was 10.


Results

From temporal to spatiotemporal model

After developing the spatiotemporal mode, we first explored the differences between the spatiotemporal model and the original temporal model by taking parameter sets with the desired behavior from the temporal model and simulating them in the spatiotemporal model. The main difference between the spatiotemporal and the temporal model is neighbor-to-neighbor cell signaling. In the spatiotemporal model, AHL is generated by an activator cell, generating a local cloud of higher AHL concentrations in the medium (Figure 11).

Figure 11:After three cells were activated by antigen, AHL (black) is being produced and diffusing out of the cell, generating a local cloud of AHL.

This changes the sensitivity of the model to the ratio between activator cells and receiver cells. Under the right parameters, an activator cell can activate its nearest neighbors regardless of the total system size. In some simulations (Video 1), a single cell activated by antigen can signal to and activate the entire population, as a wave of activation travels through the system. This suggests that for specific parameter values, the QS system has a great capacity for signal amplification, where a single cell can start a chain reaction, activating the entire cell population.

Video 1: Fluorescence produced by the original system (top) versus the system with aiiA included (bottom). On the right side is shown the system induced with antigen. On the left, the system in the abscence of antigen. This system was simulated with positive regulation of LuxR. Click for parameters.

Inclusion of aiiA

We took a spontaneously lysing parameter set from the temporal model and transferred it to the spatiotemporal model. Production and degradation terms for aiiA were added to the system of ODEs, as well as degradation of AHL by aiiA. Simulations were run both with and without antigen induction. With the right parameters associated with aiiA, a parameter set with spontaneous self-activation and lysis can be adjusted into a parameter set where cells fully activate and lyse only in response to antigen (Table 3 and 4).

Table 3: A system with positive regulation of LuxR that self-activated spontaneously was modified by adding aiiA or increasing the degradation rate of LuxR, resulting in the desired behavior.
Score Animation, no antigen Animation, antigen Parameters
Original set 99.98 link link link
Original Set & aiiA 0.22 link link link
Original Set & increased LuxR deg. 0.24 link link link
Table 4: A system with negatively regulation of LuxR that self-activated spontaneously was modified by adding aiiA or increasing the degradation rate of LuxR, resulting in the desired behavior.
Score Animation, no antigen Animation, antigen Parameters
Original set 99.97 link link link
Original Set & aiiA 0.059 link link link
Original Set & increased LuxR deg. 0.084 link link link

By adding aiiA to our model, we can simulate a model where cells are able to communicate the presence of antigen, resulting in a population-wide signal (Video 2). But in the absence of antigen, cells remain inactive and no fluorescence is generated. This suggests that aiiA can remove enough AHL to prevent spontaneous lysis, but does not remove so much that it prevents antigen-induced cell lysis.

Video 2: Fluorescence produced by the original system (top) versus the system with aiiA included (bottom). On the right side is shown the system induced with antigen. On the left, the system in the abscence of antigen. This system was simulated with negative regulation of LuxR.


Influence of LuxR stability

Analogous to adding aiiA, simulations with increased LuxR degradation rates were also carried out. All parameters are identical, except for the LuxR degradation rate, and no aiiA is present. Here, spontaneous lysis was also prevented without the system losing the ability to lyse upon induction (Video 3). As with our temporal model, this suggests another system parameter that could be altered experimentally to improve performance.

Video 3: Fluorescence produced by the original system (top) versus the system with decreased LuxR stability (bottom). On the right side is shown the system induced with antigen. On the left, the system in the abscence of antigen. This system was simulated with negative regulation of LuxR.

Summary

Two strategies for tuning the system into having the desired behavior using experimentally tunable parameters have been successfully tested in silico with both strategies giving very comparable results.

In the spatiotemporal model, the system behaved similarly to the temporal model. One major difference was that the number of activator cells can be lower while still being able to activate the entire population after detecting antigen. This is because in the temporal model, activator cells have to produce enough AHL to saturate the entire system, as AHL is instantaneously diluted homogeneously across the entire system. In the spatiotemporal model, activator cells can trigger lysis of all the cells in the population by producing enough AHL to activate their direct neighbors.

Cell density is very important to QS. However, our proposed QS system was able to function even though the cell density was increasing during the simulation. The system can therefore be robust to changes in cell density, which is important for the application of Mantis.

Furthermore, two straightforward ways to change an overly sensitive system into one with the desired behavior were demonstrated. By only increasing the rate of LuxR degradation, or by only adding the aiiA enzyme, the system could be tuned to only generate fluorescence in the presence of antigen. This is important because in our experimental subproject we encountered the same problem of spontaneous auto-activation as is predicted in many of our simulations.


Conclusion

In the global parameter space, only a very small number of parameter sets have the desired behavior. This can simply be explained by the size of the parameter space. From the lab we know that the lux operon is prone to self-activation, even at low cell densities. When visualizing the entire parameter space. the desired behavior lies in between the area representing spontaneous activation and the area representing no QS (Figure 10). By lowering the sensitivity of the QS system, a system will move towards that area of the parameter space where the desired behavior is found. We found that by changing the degradation rate of LuxR or with the addition of aiiA the system could transition from either of these two types of behaviour into the desired behavior.

We hypothesized that negatively regulated LuxR may have a faster response time and that they would be able to switch from a state where they would take up and sense AHL to a state where they would produce and release AHL. While we did find some evidence to suggest negatively regulated systems are faster, and indeed the parameter set found that had both the desired behavior and the quickest generation of fluorescent signal was one with negative LuxR regulation.

To solve the problem of spontaneous auto-activation observed in the lab, we considered introducing the aiiA enzyme that degrades AHL. Therefore, we tested this hypothesis with our mathematical model and found that by introducing aiiA, the sensitivity of a system to AHL can be reduced while maintaining the ability of the system switch from an OFF-state to an ON-state, self-activate, and cell-signal to its neighboring cells. Both for positive and negative regulation of LuxR, we were able to show that only by adding aiiA, the desired functionality can be obtained. This hypothesis was tested and confirmed in the lab.

Similarly, we found that the degradation rate of LuxR is a parameter that plays a large role in determining the system properties. As the LuxR degradation is tunable in the lab, by adding a degradation tag to the protein, we also explored increasing the LuxR degradation rate. Here we find similar results, where for some spontaneously self-activating parameters, we can obtain the desired behavior by only adding aiiA to the system. This was also successful for both positive and negative regulation of LuxR.

References

  1. Williams, Joshua W., et al. "Robust and sensitive control of a quorum‐sensing circuit by two interlocked feedback loops." Molecular systems biology 4.1 (2008): 234.
  2. Haseltine, Eric L., and Frances H. Arnold. "Implications of rewiring bacterial quorum sensing." Applied and environmental microbiology 74.2 (2008): 437-445.
  3. Miyashiro, Tim, and Edward G. Ruby. "Shedding light on bioluminescence regulation in Vibrio fischeri." Molecular microbiology 84.5 (2012): 795-806.
  4. Bassler, Bonnie L., and Richard Losick. "Bacterially speaking." Cell 125.2 (2006): 237-246.
  5. Sourjik, Victor, and Ned S. Wingreen. "Responding to chemical gradients: bacterial chemotaxis." Current opinion in cell biology 24.2 (2012): 262-268.
  6. Moukalled, F., L. Mangani, and M. Darwish. "The finite volume method in computational fluid dynamics." (2016).