Team:Wageningen UR/Model/QStest

Quorum sensing and cell lysis

We developed a mathematical model of the quorum sensing and cell lysis system. The goals of the mathematical model are the following. To characterize the theoretical properties of the system. To identify which kinetic parameters dictate the behavior of the system. To determine which parameter, that can be changed in the lab, can be varied to obtain a system with the correct behavior. To achieve this goal, the Simple model was first developed. The Simple model was used to sample the global parameter space. The second model that we used is the Spatial Model. This model was used to test hypothesis and to guide wet-lab experiments.
With the modeling of the quorum sensing and cell lysis system we had the following goals.

  • Characterize system behavior
  • Sample global parameter space
  • Identify parameters with high sensitivity
  • Find lab-tunable parameters
  • Test hypothesis

Summary

Our mathematical model shows that the QS system can spontaneously activate, function properly, or be insensitive to AHL. Changing the parameters of LuxR degradation and AHL cellular degradation allow a system to transition from one behavior to the other. In the diffusion model of our QS and cell lysis system, 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.

For the system to function properly, it needs the following behaviors.
Bistability: Cells must be able to switch from an OFF-state to an ON-state. This requires that the system of ODE equations has more than one equilibrium state solutions. 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 means unclear if the native lux box 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 quorum sensing system.[4]
Delayed cell lysis: Cells in the ON-state should lyse to release their split fluorescent protein. But they should do so only after signaling cells still in the OFF-state.

Simple Model

The main goal of the simple model is to allow for a sampling of the parameter space of the kinetic rate constants of the model. 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 assumptions have been made.

System

The system consists out of three compartments; the external medium, antigen-activated bacterial cells, and signal-receiving bacterial cells. 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 halfway during the simulation. This simulates the addition of antigen-containing blood to the bacteria and detecting 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 volume unit in 1000 total volume units. This approximates a cell density of 1.0 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 1000 cells in this system. Out of these cells, 900 cells are receiver cells, 100 are activator cells. The two cellular compartments in the model actually represent only a single cell. When cellular compartments exchange species with the external medium compartment, the total amount of exchange is calculated by using the total for each cell type. This means that all receiver cells are identical copies of each other. The same is true for the and activator cells.

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

Each cell has concentration values for the following species: AHL, LuxI, LuxR, [AHL-LuxR] complex, split fluorescent protein and lysis protein. In the external medium, there are concentration values for AHL, C-terminal split Fluorescent protein, N-terminal split fluorescent protein, and the fused and active fluorescent protein. 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. There is no spatial component to any reaction in this model.

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 much slower than transcription. Therefore, there are no concentration values for any mRNA species in the model. Similarly, transcription factors binding to DNA is also assumed to quasi steady state. The assumption here is that 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 [AHL-LuxR] is also in a quasi steady state. While there has been some debate on whether [AHL-LuxR] complex formation or [AHL-LuxR] dimer formation occurs faster. We decided to go with the assumption that dimer formation is faster and in quasi steady state. Therefore, there is no concentration value for the [AHL-LuxR] dimer in the model. To summarize, transcription factor-DNA complexes, mRNA, and [AHL-LuxR] 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, giving a Hill coefficient of 2.

The lysis mechanism is not implicitly modeled. Lysis is assumed to occur instantaneously the moment the lysis protein value exceeds a value of 1.0. Before lysis protein reaches this concentration, the lysis protein 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 out ouf of the cell. As stated earlier, the compartments are considered to be well-mixed, which equates to assuming that both the diffusion of AHL inside each compartment is instantaneous. Similarly, the split fluorescent proteins released upon cell lysis is assumed to be instantaneous.

Except for the fluorescent proteins in the external medium, all species have degradation rate. The split and fused fluorescent proteins are considered to be sufficiently stable so that degradation outside the cell plays no role on simulation length timescales. There are no values of LuxR, LuxI and [AHL-luxR] inside the external medium. Upon lysis of the cell, only AHl and the split fluorescent proteins are released into the external medium. AHL bond to LuxR and LuxI inside lysing cells simply disappear from the model.
Some of these assumptions may result in a model that fails to incorporate important dynamics. To accurately sample the 26-dimensional parameter space, the computational time of a single simulation can not be longer than several seconds. The spatial model, described in the second part of this page, alleviates some of the assumptions that are made here.


Ordinary differential equations

${\mathrm{d}AHL_{cell} \over \mathrm{d}t} = \alpha_{1}LuxI + D_{AHL}(AHL_{ext}-AHL_{cell})+k_{-1}RA-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}RA \over \mathrm{d}t} = k_{-1}LuxR \cdot AHL - \beta_{3}RA^{2} - \beta_{4}RA$

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


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


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

${\mathrm{d}Lysis \over \mathrm{d}t} = \alpha_{10} + \alpha_{11} {{LuxR \cdot RA^2} \over {km_2 + RA^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

$RA$ The LuxR-AHL complex

$RA_2$ The dimer of the LuxR-AHL complex and active transcription factor

$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

$\alpha_1$ AHL production rate of LuxI

$\alpha_2$ Leaky production rate of LuxI from pLuxA

$\alpha_3$ Maximum production rate of LuxI from pLuxA

$\alpha_4$ Leaky production rate of LuxI from pLuxB

$\alpha_5$ Maximum production rate of LuxI from pLuxB

$\alpha_6$ Leaky production rate of LuxR from pLuxB

$\alpha_7$ Maximum production rate of LuxR from pLuxA

$\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 pLuxB

$\alpha_{11}$ Maximum production rate of Lysis from pLuxB

$\alpha_{12}$ Constitutive production of $SplitFP_{cell}$

$k_1$ formation rate of [LuxR-AHL]

$k_{-1}$ dissociation rate of [LuxR-AHL]

$k_{-2}$ dissociation rate of RA

$k_{3}$ formation & maturation rate of the full fluorescent protein

$km_{1}$ Michaelis-Menten constant of pLuxA for $RA$

$km_{2}$ Michaelis-Menten constant of pLuxB for $RA$

$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 the $RA$ dimer

$\beta_4$ degradation rate of $RA$ complex

$\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 system of ODE equations 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. These parameter sets were simulated for a total of 240 time units, where between t=120 and t=180, the external signal S was set from 0 to 1. A value for 1.0 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.0, the cell would be removed from the simulation and the AHL and split fluorescent protein contents would be released into the external medium. All released fluorescent protein was always 50-50 of both C- and N-terminal versions. This was done to avoid having to introduce random numbers.

Scoring function

To properly score how well the parameter set resulted in the desired behavior, the following scoring scheme was used. The area under the curve of the fluorescence timeseries was approximated using the numpy.trapezoid() function. As an ideal parameter set generates no fluorescence before t=120 and very quickly generates a lot of fluorescence after t=120, the area under the curve before t=120 should be minimal, while after t=120, the area should be as large as possible. The area under the curve represents both how much, as well as how quickly, fluorescence is being generated. Secondly, in absence of an external input, so in case S = 0 during the entire simulation, there should be no fluorescence generated. Fluorescence in 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=120 and T=180. The final score of a parameter set is calculated by combining both simulation runs. As the objective function has to be minimalized. This leads to the following equations for the scoring function.

$Score_{simulation} = \displaystyle\int_{t=0}^{t=120} F \mathrm{d}t + {1 \over {\displaystyle\int_{t=120}^{t=240} F \mathrm{d}t + 0.01}}$

$Score_{set} = {900 \over { |Score_{noAntigen} - Score_{Antigen}|+ 0.9} } + 0.1 \cdot Score_{Antigen}$

Results

Global parameter sampling

A total of 1 million parameter sets were sampled for the system with positively regulated LuxR.

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 778326 77.8%
Activator lysis - -
Spontaneous lysis 10986 1.1%
Induced full lysis 21 0.002%
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.

Figure 2: Comparing parameter sets that lyse and those that do not lyse, low values for the degradation rates of the [LuxR-AHL] dimer, the [LuxR-AHL] complex, and LuxR are overrepresented in sets that lyse. Low degradation rates of species containing LuxR are associated with systems are more likely to lyse. The same is true for the LuxI protein, responsible for AHL production.

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

Figure 3: Comparing parameter sets that lyse and those that do not lyse, low values for the degradation rates of the [LuxR-AHL] dimer, the [LuxR-AHL] complex, and LuxR are overrepresented in sets that lyse. Low degradation rates of species containing LuxR are associated with systems are more likely to lyse. The LuxI rate of parameter sets that lyse are uniformly distributed, suggesting the LuxI rate has no significant effect.

As the number of parameter sets that quorum sense and lyse in response to antigen were very low in both cases, it was impossible to say something about the difference between the parameters with the desired behavior and those that lysed spontaneosly. Therefore, the local parameter space around these parameter sets was explored. By scoring parameter sets in the local vicinity of the parameter sets discovered by the global search, parameter sets with the desired behavior could be compared to those that lysed spontaneously. While the two types of parameter sets were still quite similar, a different behavior could be seen for the extracellular degradation of AHL and the degradation rate of LuxR (figure 4 and 5).

Figure 4: PosReg ExtAHLdeg
Figure 5: PosReg LuxRdeg
Figure 6: NegReg ExtAHLdeg
Figure 7: NegReg LuxRdeg

Spatial Model

As the assumption that AHL diffuses instantaneous breaks down when the system becomes larger, or when the number of activator cells becomes lower, we developed a spatial model that incorporates both diffusion of AHL, but also of the mixing of the split fluorescent proteins.

This spatial 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. Diffusion occurs through the face that each square has with it's top, left, right, and bottom neighbor.

The finite volume method was implemented in the following manner. For squares with four neighbors, the sum over all of the concentration values of the four adjacent squares mimus four times the concentration value in that square.

$U_{i,j} = p_{i-1,j} + p_{i,j-1} + p_{i+1,j} + p_{i,j+1} - 4 \cdot p_{i,j}$
$U$ is the scalar value of the vector field describing the flux in square with coordinates i and j
$p_{i,j}$ is the concentration in square with coordinates i and j

$U_{i,j}$ is positive when the concentration in a square is increasing. In that case, more flux enters the square than leaves the square. As there are four faces through which flux can flow, 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 that describes the d

The change in concentration is calculated by multipliying U with the diffusion constant.

$f = D \cdot U \Delta t$
$f$ is the change concentration in a square
D is the diffusion constant.
$U$ is vector describing the flux

The divergence describes how large the gradients between squares are. THe diffusion coefficient describes how quickly flux flows along the gradient field.

Anim pariatur cliche reprehenderit, enim eiusmod high life accusamus terry richardson ad squid. 3 wolf moon officia aute, non cupidatat skateboard dolor brunch. Food truck quinoa nesciunt laborum eiusmod. Brunch 3 wolf moon tempor, sunt aliqua put a bird on it squid single-origin coffee nulla assumenda shoreditch et. Nihil anim keffiyeh helvetica, craft beer labore wes anderson cred nesciunt sapiente ea proident. Ad vegan excepteur butcher vice lomo. Leggings occaecat craft beer farm-to-table, raw denim aesthetic synth nesciunt you probably haven't heard of them accusamus labore sustainable VHS.

Results

Spatial model and aiiA

We took a spontaneously lysing parameter set from the ODE model and transfered it to the spatial model. Production and degradation terms for aiiA were added to the system of ODE's, 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).

Table 3: A system with positive regulation of LuxR that self-activated spontaneously was modified by adding aiiA, resulting in the desired behavior.
Original Set Original Set & aiiA
Score 99.98 0.22
Animation no antigen - -
Animation antigen 10986 1.1%


Conclusion

Conclusion/Discussion on the affibody experiment.

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.