Team:Cambridge/Modelling
From 2008.igem.org
x
We are interested in using peptidebased quorumsensing systems from grampositive bacteria to create selforganising biological systems. We would like to demonstrate how to engineer multicellular biological systems that are capable of expressing spatial patterns of gene expression, such as with GFP, on a bacterial lawn. This would be a first step in the direction of engineering multicellular behaviour that will become immensely useful when thinking about further applications: Instead of using GFP in our patterns, we could use genes that increase cellular electric conductivity to arrive at selfassembling systems responding to electric input and output.
In the following article, we would like to introduce a model of the agrQS of S.aureus to illustrate how a typical QS system works. Furthermore, the model will predict theoretical values for biological parameters (such as the threshold cell density) that will be verified and we will be also be able to predict how the system behaves if we change a number crucial parameters. This can be extremely useful in informing design decisions when building a synthetic device. Finally, we would like to use the model in a hypothetical setup with a second parallel agrsystem to investigate how to engineer a biological patterning systems. A .pdf download of this article can be found by clicking here: [1]
Modelling Bacterial Quorum Sensing Systems  ParametersModelling biological systems involves finding relations between a great number of variables and parameters. Usually, these systems cannot be solved in an analytical way as might be accomplished with simpler systems. The problem is compounded by the fact that we only have limited knowledge of the specific parameters involved (production rates, decay rates, etc.). In this article, we shall make a number of assumptions: 1) we work with enzyme concentrations [A] etc. and assume enzymes are wellmixed 2) we use standard enzyme kinetics to infer relations between enzyme concentrations 3) and we scale rate parameters relative to a fixed value. To expand on the last point, we note that all parameters we will be working with denote rates and thus have the same dimensions. Therefore, we can work with scalar values only. Secondly, we can scale these parameters w.r.t one of them: We shall fix the maximum rate of promoter expression at 1 and determine the other parameters in relation to this rate. For example, basal expression rates must be significantly lower than the maximum rate of promoter expression. In fact, in the agrQS system of Staphylococcus aureus that we are going to model the p2 promoter seems to exhibit a significant increase^{[2]} upon activation in wildtypes. We shall take basal rates of 0.05, implying a 20fold increase of promoter activity translating into even a greater increase in AIP production (which is 50fold according to Gustaffson^{[3]} who has also modelled part of the agrsystem). Furthermore, we assume that all enzymes decay at the same rate of 0.2. In addition, phosphorylation and enzyme complex formation rates must be significantly higher than protein expression rates. We shall take phosphorylation rates at 50 (and dephosphorylation rates at 5) and enzyme complex formation rates at 20 (deformation rates at 2). The important thing to remember is that these rates are much higher than the previous ones and will thus allow us to use quasisteady state assumptions in our model. An example of this method can be found here^{[4]}. To some degree, the quasisteady state assumption allows us to ignore actual parameter values. This comes in quite handy when modelling biological systems. To summarise, the parameter values we shall be using are: maximum promoter expression rate  1, basal expression rate  β=0.02, enzyme decay rate  δ=0.2, enzyme complex formation and deformation rates  c_{+}=20 and c_{}=2, enzyme phosphorylation and dephosphorylation rates  p_{+}=50 and p_{}=5. The agrQuorum Sensing System of S.aureus
Synthesis EquationLet us at first look at protein synthesis. The RNA Polymerase binds to the promoter and transcribes mRNA, [m], at a specific rate (assuming MichaelisMenten kinetics, which we can deduce from simple chemical kinetics theory) depending on the promoter strength and the inducer concentration, [A]. It is then translated at a rate depending on the strength of the ribosome binding site and folds into enzyme [E]. mRNA decay rates are much greater than enzyme decay rates so we can assume that [m] is in a quasisteady state. The rate equations (after scaling) can be simplified into a single equation (3) as follows: and further scaling of [E] will leave us with the rate equation, where p=1 following previous assumptions of maximal promoter expression. By assuming quasisteady state behaviour, we have gotten rid of the mRNA dynamics. We shall use this synthesis equation in the following sections. Recasting the agrreceiver part as a biobrickWe would like to have the receiver part of the agrQS system in form of a biobrick. It becomes useful to think about inputs and outputs to and from such a biobrick. We shall use PoPS as a universal signal to "connect" biobricks to each other. A constitutive promoter biobrick for instance has no inputs but only PoPS output. Our agrreceiver device can be thought of as follows: In designing the biobrick we use two principles introduced in a recent paper and wellknown in iGEM circles^{[5]}: Standardisation and characterisation. Here, our biobrick part is modular in the sense that the "interface" between biobricks consists of PoPS in and output. In essence, this demonstrates a higherlevel analogue to the biobrick restriction site standard that is used in assembling the biobrick parts. We will also look at the output, or transfer function, of the biobrick. Modelling the agrquorum sensing response/receiverTo model and test the agrreceiver, we need to put a p2 promoter (and rbs) upstream of the device we just introduced. Let [P] denote the concentration of AIP, [A] denote the concentration of AgrA and similarly for the other proteins, [CP] denote the obvious enzyme complex and [A^{+}] be AgrA phosphorylated. We also need to consider concentrations in extra and intracellular space carefully. AIP exists (or is only a chemical player) in the extracellular environment whereas cell enzymes (such as AgrA and AgrC) exist only within a cell. Assuming that the local density (i.e. the volume fraction) of cells is ρ then 1ρ will be occupied by AIP with a certain concentration. We assume that concentrations have values relative to extra or intracellular space i.e. if [A] = 0.2 when ρ=0.1 then doubling intracelullar space to ρ=0.2 while holding the total amount of AgrA constant, will result in [A]=0.1. The equations of this model read: Assuming that concentration of [CP] is in a quasisteady state we arrive at [CP] = c_{+}/(c_{} + δ} [C] [P] and can substitute it back into the system to simplify it into a 4component system: Since A and C are surface enzymes, we can also assume that their levels change much slower than the levels AIP or A^{+} implying that A and C are in a quasisteady state. We would then end up with a 2component system in P and A^{+} and could easily determine the dynamics of this system (and show that the systems has two stable steady states (low and highdensity behaviour) and an unstable steady state (the 'switch' density). In the following sections we shall, however, only use minimal assumptions and use the 4component system in modelling only. Characterising the agrreceiver deviceLet us now look at how the receiver behaves under a number of conditions. In all cases, we assume that [P] is held constant (last rate equation would thus fall away).  In the following, we can observe PoPS output at a fixed point in time while varying AIP levels. This is also called a transfer function (between input and output). Incidentally, PoPS output is identical to PoPS input in this case (both are using the same promoter). The steplike shape of the curve agrees with previously published results^{[6]}.  We can also look at a surface plot of transfer functions, adding an additional time dimension. In this case, we can also clearly see static dynamic performance  PoPS output plotted against time while holding AIP levels constant. This is important in determining the timing of the system output, which can be crucial when using this part in more complex settings. Throughout the plots we observe steplike behaviour w.r.t AIP concentrations  this is typical for quorumsensing systems that react sensitively to local cell densities. Our model seems to work. It is interesting to note that even when we change phosphorylation and complex formation rates to different scalar values (we still require them to be large), the transfer function will only change negligibly. No qualitative change is observed. Modelling the agrsender deviceTo complete the agrsystem, let us look at the sender that works as follows: AgrB and AgrD lie downstream of a promoter, AgrB is a membrane protein that forms a complex with AgrD and produces AIP. We assume the reaction is B + D <=> BD > B + P. A corresponding biobrick would consist of agrB and agrD with a terminator and would take PoPS input, but would not generate PoPS ouput since the biobrick would only change AIP levels. Again, we assume quasisteady state behaviour for BD. The model reads: where p_{AIP} is the production rate of AIP (we take it to be 1). After substituting BD into the rate equations, we end up with a threecomponent system again. In the wildtype system, AIP synthesis lies downstream of the p2 promoter (1st figure). There is again a term involving ρ to account for extra and intracellular medium. Modelling the complete agrquorum sensing systemLet us now consider the complete agrQS system as illustrated in the 1st figure. Clearly, this will be described by solving the equations for the agrreceiver and sender simultaneously (where PoPS in the sender part will now be the p2 synthesis term from before). The system of equations to be solved will consist of the rate equations of [A][D] and [P]. For the rate equation involving [P] we need to add both contributions from the sender and receiver systems and only take a single decay and basal expression term so that the resulting equation reads: In the lab one can "simply" (Biologists clearly agree with this statement) put the sender and receiver device downstream of a p2 promoter to recreate the agrQS system with biobricks. For completeness, the entire model reads: and we substitute in the appropriate terms that we arrived at using quasisteady state assumptions. Remember we took δ=0.05. We solve these equations and plot PoPS output w.r.t. the volumefraction of intracellular space in the following figure. In the case of S. aureus, PoPS output drives a virulence gene. While the curve is surprisingly robust when changing initial conditions and most parameters, it is extremely sensitive when varying basal rates of expression (assuming maximum expression is fixed at 1). The shape of the curve and particularly the threshold density at which the system switches behaviour, is therefore largely dependent on the basal rate of expression. The QSsystem switches from low into highdensity behaviour at a local density of roughly 0.005 per unit volume. While we were not able to find data that would verify this estimate, we can consider the gramnegative E.coli. In ideal conditions, E.coli grows to roughly (25) x 10^{9} cells/ml in LB until it reaches stationary phase. Assuming that a single E.coli cell has a volume on the order of 1μm^{3}, the cell density at stationary phase will be on the order of 10^{3} per unit volume. Assuming that growth in ideal conditions implies higher resulting cell densities, this number will be an upper limit to cell densities, and the E.coliQS system must be activated at roughly those densities (or slightly below). This validates that our guess of 5x10^{3} for the threshold density in agrsystems is in the correct order of magnitude. Finally, it is interesting to note that we have arrived at this result by purely theoretical considerations. We have seen that the basal rate of expression plays a critical role in the system response. This in turn will be useful in informing possible design choices for engineered applications: To change the threshold density one can simply increase basal rates of expression, e.g. by constitutively expressing the agrsystem in parallel to the agrgenes that are regulated by p2. As an example, if one was to increase basal expression rates (i.e. β) the threshold value would decrease. This makes it extremely interesting to engineered applications in Synthetic Biology. Pattern formationIs it possible to engineer bacteria that exhibit multicellular behaviour and can express spatial patterns (e.g. on a bacterial lawn)? We would like to use AIP quorumsensing systems to demonstrate how this might be accomplished. AIPbased quorumsensing systems have the advantage that there exist numerous derivatives in grampositive bacteria that are mutually exclusive. This stands in constrast with OHLbased systems (in gramnegative bacteria) that exhibit a high occurance of interference between OHL derivatives. Before we try to model such behaviour, let us look at a simple mathematical model for pattern formation. Reactiondiffusion systemsThe theory of patterns is vast and it would (and does) extend beyond the scope of this project. However, we will look at a very simple mathematical system that is capable of generating patterns. This might give us some intuition in how patterns are formed even though in our case there might be no direct application. Turing Patterns Turing's idea was to look at systems in the absence of diffusion. He required such a system to have a stable steadystate and upon addition of diffusion, required this to become unstable. A simple example is the Schnakenberg system that is analysed in Mathematical Biology by J. Murray. This is a simple system that is capable of creating diffusiondriven patterns. The mathematical analysis goes as follows: One considers the system in the absence of diffusion and looks for a stable steadystate. Next, one considers the same system with the diffusion term added and substitutes in a trial solution that is composed of a Fourier series while requiring the steadystate to be unstable. This will result in a number of further conditions (inequalities) and it is easy to find the (Turing) space of parameters that will satisfy all conditions for pattern generation. A gallery of patterns produced by Schnakenberg system can be found on the next page. We can see that the type of patterns created is highly dependent on our parameter values, particularly α and β. Can we arrive at similiar patterns with our AIPQS system?
Pattern formation in Synthetic BiologyHow can we enginner bacteria to express spatial patterns? (It is interesting to note that many bacteria do form patterns, particularly when exposed to different nutrient levels) We would like to express patterns on a bacterial lawn with constant bacterial density everywhere and we can think of patterns as spatial distributions of PoPS output (downstream of which we can put GFP) or cell differentiation. Since we have modelled the agrsystem and know that noninterfering AIPbased systems exist, we would like to construct a pattern generating system as follows: Both the AIPA and AIPB systems inhibit production of the "competitor" molecule. Furthermore, both systems are selfactivating. This would be a more complex analogue to the simple reactiondiffusion system introduced previously. In the following part, we assume that nutrients are abundant (we do not model nutrient consumption) and that the bacterial density throughout the lawn is constant and high by setting ρ=0.01. Note that we are disregarding the fact that these bacteria now live on a 2d surface. We would arrive at a more detailed estimate for bacterial densities in a bacterial lawn if we introduced a diffusion term to account for signalling molecules diffusing across the nutrient layer. Furthermore, the synthesis terms in this model are bounded in our system. Meinhardt mentioned^{[7]} stripelike patterns can be generated by such systems where activator production is bounded. While he was analysing a simple twocomponent system, we would like to generate patterns by using two agrsystems in parallel, both selfactivating and inhibiting the parallel system. The rate equations for one of the parallel agrsystems (with peptide P_{1}) read as follows: where f(x)=40/(0.1+x). The equations for the the second agrderivative system with peptide P_{2} read analogously with P_{1} and P_{2} interchanged (and independent variables for intracellular enzymes). Note that we have added a diffusion term to account for spatial diffusion across the plane. Solving those equations with random initial conditions in levels of P_{i} gives us:
Some regularity in the patches of differentiation can be observed although these "patterns" depend to a greater degree on initial conditions compared to the patterns seen previously in the Schnakenberg system. However, this is a promising indicator that demonstrates how coordinated multicellular behaviour might originate from an engineered system. Further QuestionsA concise modelWe have just tried to go from intra to intercellular behaviour by naively applying the primarily intracellular agrQS model to the intercellular case. We cannot say how accurate this approach is since we did not run experiments. However, we can consider a more concise and simpler version of our patterngenerating system. The most important feature of our envisioned systems is that both "populations" inhibit or compete against each other and that the autosynthesis term is bounded. Assuming that competion can be modelled by simple massaction, a simplified version of our model would read: Of course, one can also try other systems that incorporate the main features of our system. Further work can be done on the dynamics that these systems exhibit. Other patterngenerating systemsThere is a vast number of patterngenerating systems. While we trying to model two competitively inhibiting populations, patterns might be created by a completely different system. In particular, we are interested in the question how pattern generation might occur using a single signalling molecule only. Bounded patternsPatterning is everywhere in biology  one can think of the whole development process as "pattern formation". Here, patterns generation is switched on by higherlevel mechanisms and patterns are expressed in bounded domains (e.g. if we regard organs as patterns, these are clearly seperated from other organs). How does this work? How can we generate patterns that are constraint to a domain? Can we engineer a system that might exhibit this kind of behaviour?
ReferencesAll equations have been solved numerically with a simple difference method  this method is easy to implement and gives accurate results, even though one must be careful that stability does not become an issue.
x
