 Methodology
 Open Access
 Published:
Spatial quorum sensing modelling using coloured hybrid Petri nets and simulative model checking
BMC Bioinformaticsvolume 20, Article number: 173 (2019)
Abstract
Background
Quorum sensing drives biofilm formation in bacteria in order to ensure that biofilm formation only occurs when colonies are of a sufficient size and density. This spatial behaviour is achieved by the broadcast communication of an autoinducer in a diffusion scenario. This is of interest, for example, when considering the role of gut microbiota in gut health. This behaviour occurs within the context of the four phases of bacterial growth, specifically in the exponential stage (phase 2) for autoinducer production and the stationary stage (phase 3) for biofilm formation.
Results
We have used coloured hybrid Petri nets to stepwise develop a flexible computational model for E.coli biofilm formation driven by Autoinducer 2 (AI2) which is easy to configure for different notions of space. The model describes the essential components of gene transcription, signal transduction, extra and intra cellular transport, as well as the twophase nature of the system. We build on a previously published nonspatial stochastic Petri net model of AI2 production, keeping the assumptions of a limited nutritional environment, and our spatial hybrid Petri net model of biofilm formation, first presented at the NETTAB 2017 workshop. First we consider the two models separately without space, and then combined, and finally we add space. We describe in detail our stepwise model development and validation.
Our simulation results support the expected behaviour that biofilm formation is increased in areas of higher bacterial colony size and density. Our analysis techniques include behaviour checking based on linear time temporal logic.
Conclusions
The advantages of our modelling and analysis approach are the description of quorum sensing and associated biofilm formation over two phases of bacterial growth, taking into account bacterial spatial distribution using a flexible and easy to maintain computational model. All computational results are reproducible.
Background
Human gut microbiota. The human body hosts a large number of prokaryotic and unicellular eukaryotic cells, which are found in most tissues. This huge population can be up to 100 trillion in the human body, almost the same number as human cells, with the total mass of around 0.2 kg [1].
Human gut microbiota play an important role in human health and are necessary for a human to survive. For instance, human gut microbiota help in vitamin synthesis, dairy digestion, gut development, nutrient processing, resistance to pathogens, brain development and function, immune cell development and immune responses [2]. It is notable that some of these bacteria in human intestine can prevent certain diseases such as diabetes [3] or liver diseases [2]. In addition, gut microbiota have a direct effect on the diseases related to intestines such as Irritable Bowel Syndrome (IBS), Inflammatory Bowel Disease (IBD) [2], metabolic diseases and obesity [4].
The composition of gut microbiota is unique for each individual host. However, family members can have more similar gut microbiota compared to unrelated people. The composition of gut microbiota is influenced by diet, drugs, lifestyle, stress and genetics [2, 4].
Quorum sensing. Bacteria need to act in a concerted manner in order to achieve large scale effects in, for example, a multicellular host such as human. This involves communication, and the detection of the presence of other bacteria of the same species or strain whose activities can be coordinated and aggregated. This process of stimulus and response correlated to population density is called quorum sensing which can happen among microbes, particularly among the bacteria of the same strain, with other strains and with human cells [5]. This communication mechanism between cells is based on chemical broadcasting, i.e. oneway communication whereby molecules acting as chemical signals are transmitted via diffusion without requirement for an acknowledgement of receipt. The concentration of the molecules involved stands as a proxy for population density. The receiving cells detect these molecules via their membrane receptors which employ a threshold mechanism over the concentration of the molecules. This results in the activation of pathways which induce changes in gene expression.
Bacteria use quorum sensing to initiate some changes in their behaviour or phenotype, such as biofilm formation, pathogenicity and virulence factor production, motility and toxin production. These changes play an important role in the relationship between the host and the microorganism [6].
This process is achieved by the production, release and detection of a signalling molecule called autoinducer (AI) [7]. According to Elias and Banin, distance has an impact on the ability of the cells to intercommunicate. Their study shows that the distance between the cells and hence the compactness of the society of the bacteria is a more important factor than the size of the population itself [8]. When the population density is low, the total amount of AI produced is too low for the quorum sensing process to be enabled. As soon as AI is produced and exported, it diffuses into the environment. Thus it can not be detected by the individual cell to activate biofilm formation. An increase in the population increases the amount of AI observable in the intercellular environment, facilitating the ability of the cell to detect and import it. AI accumulates in the cell and can activate transcriptional signals by binding to a regulator protein which is a gene transcriptional suppressor. This is when quorum sensing is switched on [6].
Bacterial population growth occurs in four different phases, see Fig. 1. In the first stage, Lag Phase, the bacteria adapt to the environment, and growth is not fast. In the second stage, Exponential Phase, the bacteria start to grow fast such that their number doubles logarithmically. In the third stage, Stationary Phase, the growth rate and death rate of bacteria equalise due to limited availability of growth factors and the accumulation of waste products. The final Death Phase takes place when waste products dominate the environment and the growth rate dramatically decreases while death rate is high. Autoinducer production occurs in the Exponential Phase (Phase 2) and the response of the bacteria to it occurs in the Stationary Phase (Phase 3) [9].
Looking at the phases of bacterial growth in detail, AI is produced and exported to the intracellular environment during the midExponential Phase and early Stationary Phase. In the case of a sufficiently large population of bacteria, the level of AI exceeds its threshold, which happens during the late Phase 2 and the start of the Phase 3, when the bacteria start to import this chemical to their cytoplasm [10]. At this point (Phase 3), AI production is slow while the import is so rapid that no AI in the environment is detectable [11].
In Gramnegative bacteria, there are three classes of AI according to the way they are produced: AI1, AI2, AI3 [12]. The main focus of this study is on Autoinducer2 (AI2) which is the only AI molecule that is produced and detected by both Gramnegative and Grampositive bacteria [10]. AI2 is synthesized by LuxS or its homologues from methionine and is actively transported to the extracellular environment throughout cell growth [12, 13].
When nonpathogenic E.coli produce AI2 in the intestine, an inflammatory response is initiated which is quickly arrested. Not only can AI2 act as an interkingdom signal between bacteria and their hosts [12], but it is also nonspecific between bacterial species. Thus when the AI2 system is activated, the signal molecule is quickly removed from the environment by absorbtion by the receiving cells in order to avoid crosstalk with other strains of bacteria. This mechanism prevents these strains from using AI2 to change their behaviour — for example, interacting inappropriately with the host. In E.coli, AI2 affects virulence factors, motility, pathogenicity and biofilm formation [12]; in this paper we concentrate on the latter.
Biofilm formation. Quorum sensing can cause changes in behaviour and morphology of a bacterium. One of the consequences is the formation of biofilm. Biofilm formation is a structural process where a polymer layer of carbohydrates or small amounts of protein and DNA, are produced and exported to the microenvironment [14]. This process results in phenotype changes while enabling the bacteria to adhere to each other and the surface.
The bacteria in this community are structurally organised by the way in which they communicate with each other and respond to one another’s signals. The formation of biofilm can have several beneficial effects for bacterial colonies: it (i) can act as a physical protection against antibiotics, thus facilitating bacterial antibiotic resistance, (ii) helps them in colonization of the host by facilitating adherence of the bacteria to host tissues, (iii) is helpful in collecting nutrients from the environment [15], and (iv) facilitates DNA exchange between bacteria in the biofilm which results in changes in DNA recombination or DNA repair [14].
Biofilm production leading to intercell matrix formation does not form in sparse or scattered colonies because the distances between the bacteria are too high to permit high enough autoinducer concentration to pass the threshold required to trigger the relevant internal pathways [16].
The study of biofilm is an important field for microbiologists and immunologists since it can be a threat for human health. Biofilms are involved in 65% of hospital infections, causing serious problems due to their high resistance to antibiotics [15].
In a microbiology lab, one of the constraining resources is time. In order to obtain results a researcher should wait at least for 24 h for a sufficient population of bacteria to have developed. By simulating biological systems insilico it is possible to predict and analyse the behaviour of biological systems faster. This method, aside from being less time and money consuming, gives us a clear detailed result to study and to guide lab experiments.
Related work
The first description of quorum sensing was given in the early 1970s for Vibrio fischeri [17, 18] and many scientific descriptions have subsequently been published; see the review by Waters and Bassler [19] for an overview of this large area.
Many mathematical models of the quorum sensing mechanism have been proposed, and are reviewed in depth in [16, 20]. Models of the intracellular molecular mechanism at the single cell level have mostly been continuous, using ordinary differential equations (ODEs). Some models are stochastic, for example Weber and Buceta [21] give both a deterministic model to describe AI2 production in E.coli which treats the cell as one entity, and a stochastic model which divides the cell into individual compartments, including a noise term on the luxR gene expression, which depends on the cell density and may influence phenotypic changes stochastically. Another set of models describe selfcontrolling mechanisms in the quorum sensing process, where the general purpose of AI systems is held to be (i) the homeostatic control of costly cooperative behaviours [22] and (ii) that bacteria have evolved mechanisms to repress certain components of quorum sensing if needed, see e.g. [23].
One of the most complete works in this field is the stochastic Petri net model of E.coli AI2 production, constructed by Li and his colleagues [24]; it has quantified values for the reaction rate constants and the initial concentrations of species. This model describes the AI2 production pathway as outlined in “Phase 2 model – AI2 production” section, which forms one component of our overall model.
The modelling of biofilms has a long history over the last 30 years, and has been extensively reviewed in [25]. The representation of space is a critical component of such models, and initially models comprised partial differential equations (PDEs) modelling a biofilm as a flat layer; subsequently multidimensional models have been developed to describe nonuniformities.
The mechanism of biofilm production itself has been described in papers by Li et al. [26] and Novak et al. [27]. However the descriptions are textual only without any quantified values for concentrations or rate constants, and no formal models are given.
Janowski et al. [28] have modelled both the production of autoinducer signals as well as their diffusion and cellular import in a quorum sensing scenario using standard Petri nets. They achieve a notion of space by constructing a large model explicitly comprising several copies of the Petri net descriptions of individual cells, connected by the diffusion of the AI molecules outside the cell. One drawback of this modelling approach is that each copy of the cellular mechanism needs to be updated if the intracellular description is changed during model development, and also the extracellular connections are explicitly represented individually and likewise would require individual updating.
Other approaches to describing quorum sensing at the population scale include agentbased models which permit the allocation of a range of conditions to different individuals. For example, Müller et al. [29] have proposed a spatial singlecell model of quorum sensing using two layers: intracellular and intercellular, and a population of such single cells, where the quorum sensing signalling system is modelled by nonlinear ODEs and PDEs are used to give spatial structure. P systems exploit a hierarchical membrane approach for modelling and have been used by PérezJiménez and RomeroCampero [30] to describe quorum sensing in terms of an agentbased model where an individual agent or cell is described by a P system exploiting Gillespiestyle [31] stochastics.
Our contributions
We developed a spatial model describing E.coli biofilm formation driven by the Autoinducer 2 (AI2) covering Phase 2 and Phase 3 of bacterial growth and comprising different abstraction levels. The model is represented as coloured Petri net, which can be equally read as stochastic or hybrid Petri net. We investigate different model configurations comprising the single bacterium, dense colony and sparse colonies of different sparseness, by simulating them with the approximative deltaleaping stochastic simulation algorithm [32].
We describe in detail our stepwise model development and validation strategy. All results are reproducible. Our approach of encoding space is flexible, it can easily be configured for the 1D, 2D or 3D scenario, and conveniently adjusted to different notions of space, such as different boundary conditions or neighbourhood relations [33]. Our modelling strategy can be equally applied to other problems evolving in time and space.
Outline. This paper is organised as follows. In the next section we describe our methods, comprising standard Petri nets in the stochastic, continuous and hybrid paradigms, as well as their coloured counterparts, and simulative model checking of linear time temporal logic properties. Next, we present our model and its stepwise development and validation; first we consider the model components separately without space, and then combined, and finally we add space. We conclude our paper with a discussion and conclusions suggesting further work. We also provide Additional files with supplementary material illustrating in more details our methodology and results.
Methods
Petri nets
To simulate a model in silico, a powerful tool is needed, which is able to show the complex interactions within a cell [34] in the simplest way possible. At the same time, this tool should be easy to work with and provide the required data such as graphs or data tables to be analysed. We chose to use Petri nets as they support a graphical and mathematical environment in silico; they ease communication among professionals with a diverse background and permit to computationally simulate a biological system.
Standard Petri nets are inherently discrete and free of any notion of time, which means that they consider behaviour under any timing constraints, supported by a substantial body of Petri net theory [35]. They have been proven to be useful for a wide range of applications, among them biochemical networks, such as metabolic networks, gene regulatory networks, or a combination of them [36]. If required for clarity, we call them qualitative Petri nets (QPN).
Over the years, Petri nets have been extended in various ways, for example by assigning firing rates to transitions, which can either be read as stochastic or deterministic rates. Stochastic rates yield stochastic Petri nets (SPN), which can be seen as highlevel descriptions of continuoustime Markov Chains [37], while deterministic rates yield continuous Petri nets (CPN), which in turn can be seen as a structural descriptions of ordinary differential equations (ODEs) [38, 39]. Combining stochastic and continuous rates within one model brings us to hybrid Petri nets (HPN), which contain SPN and CPN as proper subsets [40]. Firing rates are often statedependent and follow specific kinetic laws. In the case study considered in this paper, all firing rates follow massaction kinetics [41].
We use the Petri net tool Snoopy [42] which supports Petri nets in different paradigms – as required for the case study discussed in this paper, comprising QPN, SPN, CPN, and HPN, as well as their coloured counterparts, see next section. Importantly, it is possible to conveniently convert the individual net classes into each other according to the user’s need.
In the following we assume basic knowledge of Petri nets; see [43] for a gentle introduction in the context of Systems Biology, and [44] for formal definitions.
Coloured Petri nets
Enriching Petri nets by a colouring concept yields a kind of highlevel Petri nets. Colours are organised in colour sets, which can be seen as a synonym for (finite) discrete data types as known from traditional programming languages.
Coloured Petri nets allow for concise descriptions of similar network structures, which in turn permits, among others, to conveniently encode space, as we have first shown in [45] and later applied for various notions of space, for example to deal with planar cell polarity requiring the nesting of two spatial concepts [46], or to investigate phase variation in bacterial colonies, where we applied Cartesian and polar coordinates [47]. The tutoriallike paper [33] showcases the use of coloured continuous Petri nets to encode partial differential equations.
A special strengths of this approach of colouring space is its flexibility: models can be conveniently adjusted to different notions of space, such as different boundary conditions or neighbourhood relations [33]. As we will see, models can also be designed in such a way that they subsume as special cases the onedimensional, twodimensional or threedimensional settings. Thus, the basic modelling strategy can be applied to a large variety of problems evolving in time and space. See [48] for a recent review paper illustrating the wide use of coloured Petri nets for multilevel, multiscale, and multidimensional modelling of biological systems.
Coloured Petri nets consist, as standard Petri nets, of places, transitions and arcs. Additionally, a coloured Petri net is characterised by a set of discrete data types, the colour sets, and related net inscriptions.

Places get assigned a colour set and may contain a multiset of distinguishable tokens coloured with a colour of this colour set.

Transitions get assigned a guard, which is a Boolean expression. The guard must be evaluated to true for the enabling of the transition.

Arcs get assigned an expression; the result type of this expression is a multiset over the colour set of the connected place.
We consider diffusion in 3D space to recall the basic concepts of coloured Petri nets, as supported in our toolkit in different paradigms, i.e. as coloured Petri nets (no notion of time), coloured stochastic Petri nets (coloured SPN), coloured continuous Petri nets (coloured CPN), and coloured hybrid Petri nets (coloured HPN).
In order to encode spatial information we introduce the following.

Constants D1, D2 and D3: integers which define the length of the X, Y and Z axes of the environment space.

Constant D: an integer to be used to define uniform lengths of the nonzero axes. For example, in order to model a 3D cube, they are all set to the same value D. For a square 2D grid, D1 and D2 are set to D while D3 is set to 1, and for a 1D linear array, D1 is set to D while both D2 and D3 are set to 1.
Having these constants we can now define the following colour sets specifying a 3D grid, and a colour function specifying the neighbourhood relation over the grid positions, which are triples of the X, Y, Z coordinates.
The advantage of this approach relying on colour is to concisely define a potentially large model which has many repeated elements, one example of which is a spatial grid whose elements can be addressed by colour tuples specifying the grid positions, and movement in space translates into recolouring of tokens. See Fig. 2 for a Petri net showing diffusion in 3D space which can be easily configured for 1D or 2D by merely changing the constants D1…D3. This coloured representation can be automatically unfolded generating the corresponding uncoloured Petri net in standard notation amenable to simulation and analysis techniques available for uncoloured Petri nets. See Additional file 1: Figures S2–S4 for some examples of unfolded nets.
We can further define a subspace, called ‘region’ for brief, within the overall environment, using the following constants:

Constants M1, M2 and M3: integers representing the central position of the X, Y, Z dimensions of the environment,

Constant R: an integer controlling the radius of a region in each dimension (assuming a circular shape which maps to a rectangle in our representation). Note that when R=1, the side of the region is of length 3 and hence a 2D region comprises 3×3=9 grid positions and a 3D region comprises 27 positions. When R=2, the side is of length 5 giving 25 positions in the 2D case and 125 in the 3D case, and so on.

Constants S1, S2 and S3: integers controlling the sparseness (density) of the bacterial colony in the X, Y and Z dimensions, respectively, where S _{i}=1 results in no gaps in the region, S _{i}=2 generates a region spread s.t. there is one empty grid position between each occupied position, S _{i}=3 gives 2 empty grid positions, etc.
This allows us to introduce the following colour set and colour function:
which defines a region around the centre (x0,y0,z0) with a max distance from the centre of ‘radius’, and a space of xd, yd, zd inbetween in the X, Y, Z axes, respectively. In our current scenario, the main grid represents an environment within which a bacterial colony is located using the region function.
The entire diffusion mechanism with its associated colour definitions will be one of the three components of our final model, discussed when adding space – see “Adding space” section. A complete selfcontained description is given in the Additional file 1.
Temporal logics
Model checking can be used to determine whether the behaviour of a model conforms to some desired properties specified in temporal logic. We use here Probabilistic Lineartime Temporal Logic with numerical constraints (PLTLc) [49], based on Lineartime Temporal Logic (LTL) [50], extended with probabilities [51] and numerical constraints over real value variables [52]. Several features of PLTLc facilitate the expression of the behaviour of biochemical networks, including the ability to express properties relative to an absolute time value or range, the use of functions which compare the concentration of a protein to its peak value, and the derivative function enabling the description of transient, sustained or oscillatory behaviour.
The operators of PLTLc include the usual first order logic connectives: & (And), ∨ (Or), ¬ (Not), and → (Implication). In addition the language includes the following temporal operators over formulae: F(ϕ) (ϕ holds eventually (finally)), G(ϕ) (ϕ holds forever (globally)), ϕ_{1} U ϕ_{2} (ϕ_{1} holds until ϕ_{2} holds).
In this research we have used simulative model checking over time series traces of the species in the model which include the small molecules, metabolites, proteins, RNAs, genes and complexes. This can be done for single output of a continuous model, or for the several runs of a stochastic model – either separately or averaged. In our case, the behaviours that we check can also be generated by a hybrid model yielding both deterministic traces (for the continuous places) as well as stochastic traces (for the discrete places). We give the properties in PLTL format following the usage in the simulationbased MC2 model checker [49] which we have employed, with the results belonging to the set {0,1} rather than in being in the range [0..1], when applied to continuous or averaged traces.
MC2 has two builtin functions, the first of which, max, operates over all the values of a species to return the maximum of the species’ value in the simulation run, thus the peak of a species can be expressed using max(Protein). The other function d returns the derivative of the concentration of the species at each time point, thus increasing and decreasing species values can be expressed as d(Protein)>0 and d(Protein)<0, respectively.
Commonly, the property of interest refers to known values, often motivated by observations in the wet lab; e.g. The concentration of metabolite A is always below a certain threshold, for example 100, and always decreasing. This can be expressed by the temporal logic property:
P_{≥1}[ G(A<100 & d(A)<0) ].
In the research reported in this paper we often are presented with a large number of species, for example all those in the particular network of interest in one bacterium, and would like to know which of these fulfil a particular property, then we can express the previous property by:
P_{≥1}[ G($x <100 & d($x)<0) ], where $x is a meta variable ranging over all species in the model. This feature can also be applied to check all the species in the often very large models generated by uncolouring a coloured model by unfolding. As result we obtain the set of all variables fulfilling the given property, with the set being empty, when no variable fulfils the property.
Results  models
In this section we show how our model has been developed using a stepwise approach. First we consider AI2 production (Phase 2 model) and biofilm formation (Phase 3 model) separately as two models without space, then we combine them, and finally we add space to the combined model. Each step involves model validation.
Phase 2 model – AI2 production
We have used the nonspatial stochastic Petri net model of AI2 production given in Li [24] to describe the production of AI2 in Phase 2, keeping the assumptions of a limited nutritional environment. The model does have elements of AI2 activity in Phase 3, which we separate out into our description of biofilm formation (see “Phase 3 model – biofilm formation” section). The model is available in SBML format from the website of the PMC journal where the paper was published. By importing this SBML specification, we obtained the Petri net model shown in Fig. 3, which can be read as a stochastic Petri net or a continuous Petri net alike.
The detailed mechanism of the AI2 production pathway is as follows: AI2 is derived from Sadenosylmethionine (SAM) which is a methyl donor for the methyl transferase enzyme, which produces a methylated product from SAM called Sadenosylehomocysteine (SAH) which can be toxic if accumulated. Therefore, the cell rapidly acts on SAH to produce adenine and Sribosyl homocysteine (SRH) using the nucleosidase Pfs. While transforming SAH to SRH, Pfs uses a water molecule to produce Adenine. In the next step the LuxS protein catalyses reactions which use SRH to produce both homocysteine and 4,5dihydroxy2,3pentamidine (DPD). The effect of this is to create a feedback loop in the pathway by transforming homocysteine to methionine, and also to rearrange DPD to eventually produce more AI2 [24].
On another branch of the pathway, SAM decarboxylase converts SAM to Decarboxylated SAM (DeSAM) by releasing CO_{2}. After that Spermidine synthase transforms DeSAM to 5 ^{′}methylthioadenosine (MTA) using putrescine and producing spermidine. The Pfs enzyme also converts MTA to 5 ^{′}methylthioribose (MTR) while using a water molecule and producing polyamines; MTR is later excreted from E.coli [53]. Note that ubiquitous molecules such as CO_{2} and water are not represented in our model.
The model has 21 places, and 23 transitions. The Li model has a specific finite number of tokens in each place and specific rate constants, all being derived from experiments [24]. There are two source places, i.e. with only outgoing arcs and thus an always limited number of tokens, reflecting limited environmental resources: Nutrients and Putrescine. Two other places maintain a constant value: the Pfs and LuxS genes, modelling the fact that genes are not used up by transcription. In summary, the model describes the pathway starting from the sources places to AI2_Out which represents the AI2 produced in Phase 2 excreted to the environment. See Tables 1 and 2 respectively for an explanation of the names of the places and transitions in the Phase 2 model.
In the following, place and reaction names are the short form given in the actual model and indicated in italic.
Phase 3 model – biofilm formation
Our Phase 3 biofilm formation model follows the wording in the papers by Li et al. [26] and Novak et al. [27], neither of which presents a formal model, but rather provide textual descriptions of the various submechanisms. In the following we consider the production of the biochemical components of the biofilm matrix rather than the formation and structural properties of the matrix itself.
In order for the process of biofilm formation to be activated, bacteria actively transport AI2 into the cell if they detect that the external concentration of the molecule is above some threshold, implying that the population of bacteria is high enough to initiate this pathway as well as other pathways [26]. Then the next stage of quorum sensing starts; AI2 enters the bacteria in the stationary phase through a transporter complex called LuxS Regulator ABCD (LsrABCD). Thus, when there is insufficient AI2, the cell will not detect it and biofilm formation will not be initiated [6, 54].
As soon as AI2 enters the cell, it is phosphorylated by the kinase LsrK. The phosphorylated AI2 counteracts the repression of the lsr genes by LsrR, by binding to and thus separating the repressor LsrR from the lsr genes [26, 53]. This allows the lsr genes to actively produce the proteins LsrABCD, LsrK and LsrR, all of which are degraded.
The biochemical components of biofilm (carbohydrates etc.) – hereafter called ‘biofilm’ for short – can be formed in two ways; firstly by the interaction between AI2_In and QSeBC via a complex process represented overall by the transition BiofilmFormationAI2. QSeBC (Quorum sensing E.coli regulator B and C) is a coupling mechanism that senses and responds to environmental changes in order to control biofilm formation [27]. The alternative route for the formation of biofilm is via LsrR which requires the presence of AI2_In as well [26], again via a complex process represented overall by the transition BiofilmFormationLsrR.
See Fig. 4 for a graphical representation of the model, which comprises 11 places and 15 transitions; in addition we have introduced a control mechanism (“Go system”) to ensure that AI2 will enter the cell only when a certain threshold is reached [54], and also that AI2 can build up in the environment and consequently diffuse as a signal to other cells in the environment. This models the fact that under low population densities AI2 will keep diffusing in the environment and never enter the cell, with the consequence that biofilm will never be formed [12].
In order to impose a threshold on the activity of the Transport_in transition, we exploit two special arcs available in Extended Petri nets. Special arcs always connect a place to a transition. Read arcs (represented by black circles as arrow head) test if the marking of the preplace is larger or equal than the arc weight, which we specify by the constant THR. Inhibitory arcs (represented by hollow circles as arrow head) test if the marking of the preplace is smaller than the arc weight; in our case we use (so far) the weight 1, which is usually not shown. Both special arcs influence the enabledness of a transition, but do not change the marking of the tested places upon firing. Combining both as shown in the Fig. 4 (yellow subnet) brings the required control mechanism. As we want to obtain a sharp on/off mechanism, we implement the Go control mechanism as a discrete component. The stochastic transition open is enabled if m(AI2_Out)≥THR∧m(go)<1, and the stochastic transition close if m(AI2_Out)<1. Obviously, both conditions preclude each other for any THR≥1; thus there is always at most one token on the discrete place go, which appears as a sidecondition for the transition Transport_in. Consequently, the rate (following the massaction pattern) is zero, if m(go)=0, and else the rate depends only on the other two preplaces (AI2_Out, LsrABCD).
Due to the stochastic nature of the go mechanism, the Phase 3 model can either be read as an entirely stochastic model, or as a hybrid model keeping the go mechanism stochastic and the remainder being continuous. Moreover, the model implements the repressionderepression transcription cycle of lsr genes which also should be modelled stochastically because gene expression is a fundamentally stochastic process [55]. This suggests that the two places and transitions involved in this cycle can be kept as stochastic in the hybrid model as well.
See Tables 3 and 4 respectively for an explanation of the names of the places and transitions in the Phase 3 model.
Please note that this is a new model and no precise rates, rate parameters or initial concentrations are available in the literature. We assume massaction kinetics, and the rate parameters have been uniformly set in three categories, see Table 5 for details. The initial marking comprises four categories, see Table 6 for details.
Combining phase 2 & 3 models
We created an overall model, see Fig. 5, describing AI2 production and biofilm formation by combining the Phase 2 model based on Li [24] described in “Phase 2 model – AI2 production” section with the Phase 3 model for production of biofilm described in “Phase 3 model – biofilm formation” section. This was achieved by linking these two models using the place AI2_Out, which we represented in each model as a logical place (a shared place appearing in both models). To keep all place names unique, we named AI2_In from Phase 2 as AI2_In_Phase2, and likewise for AI2_In in Phase 3. In preparation for the spatial interpretation of the combined system, we provide the basis of the AI2 broadcast mechanism. For this purpose we add a sink transition Diffusion summarising diffusion to the neighbourhood in any direction; the kinetic constant for diffusion is multiplied by the number of neighbours (2, 8 and 26 in 1D, 2D and 3D, respectively). In this way, the Diffusion transition will enable transmission of the broadcast signal via dispersion thus facilitating longrange communication in the spatial model below.
We also adjusted the estimated rates in the Phase 3 model to concord with the general ranges of the specific rates in the Phase 2 model, low being set to 10^{−4}, medium to 10^{−3}, and high to 1. Note that there are two thresholds, one for open and one for close – for the time being we keep the threshold for close at its default value of 1, while varying that for open.
Adding space
In this section we reuse the coloured Petri net definitions introduced in “Coloured Petri nets” section to describe diffusion in space, which are summarised in the Additional file 1.
In order to add spatial information regarding both the environment of the bacteria, and the size, position and density (“sparseness”) of the bacterial colony we need to define the environment size over a 1, 2 or 3D grid, and the colony as a smaller grid within this environment. The topology of both the environment and the bacterial colony (or “cluster”) are by default linear in the 1D case, a square in 2D and a cube in 3D. We first export the combined model obtained in the previous section to a coloured Petri net and define all colourrelated definition as discussed in “Coloured Petri nets” section. Next we add the following annotations to the coloured model.

We uniformly assign the colour set Grid3D to all places. Thus, every (coloured) place may hold coloured tokens which are colour triples (x,y,z), the x, y, z components of which are defined over a range from 1 to D1, D2 or D3, respectively.

All (coloured) places are initialised with markings using the region function. This will allocate bacteria to all (x,y,z) positions fulfilling the region criterion as specified by its parameterised call region(x,y,z,M1,M2,M3,R,S1,S2,S3).

We enhance the definition of the Diffusion transition, making it two way in order to model the longrange reception of broadcast AI2 signals.

All arcs get assigned the arc expression (x,y,z), with one exception: one arc adjacent to Diffusion has to get the triple (a,b,c) to provide the parameters for the neighbourhood function serving as guard, see next item.

All transitions get assigned as guard the function region, with one exception: the transition Diffusion gets the guard neighbour3D26. This ensures that the combined biofilm model for Phase 2 and Phase 3 will only be generated upon unfolding) for those grid positions as specified by the parameters of the function call region, and only AI2_Out can diffuse over the whole of the environment grid.
See Fig. 6 for the resulting model and Additional file 2 for its complete CANDL specification.
Discussion  model validation
Sound model engineering requires a stepwise model development strategy. Here, we present selected aspects of our stepwise model validation.
Phase 2 model validation
Note that we consider the Li model as being validated by its publication, but take this opportunity to demonstrate the power of our analysis techniques. We start with standard Petri analysis techniques as known from Petri net theory; see [44] for a summary and related formal definitions in the context of systems biology.
When adding source and sink transitions to model an environment ensuring infinite in/out flow, which makes the Petri net transitionbordered (no sink/source places anymore), the net becomes covered with Tinvariants, which is known as a generally criterion for “all transitions may be involved forever in some basic elementary behaviour”.
There are two trivial 1Pinvariants (made of a single place, having always exactly one token), which correspond to the two genes. Due to the added source transitions, the net is structurally unbounded, thus, without timing constraints (i.e. forgetting the rates), all other places (besides those involved in Pinvariants) are unbounded.
There are only the two trivial siphons, established by the two 1Pinvariants; consequently, the Siphon Trap Property (STP) holds (because a Pinvariant is a siphon and a trap as well). As the net also belongs to the net class Extended Simple (ES), we can conclude that the net is live, i.e. all transitions will forever contribute to the system behaviour. Please note, this conclusion for a structurally unbounded net is possible thanks a structural criterion (STP & ES), and this is true for any timing constraints because it is known that an ES Petri net is timeindependently live [56]. See Additional file 3 for more details of these structural analysis results.
We continue with simulative model checking. It is known that the results may depend on the simulation trace, specifically on the step size of the ODE solver and the granularity of the recorded trace, see, e.g., [49, 57]. We chose to simulate the continuous Phase 2 model using the classic RungeKutta solver in Snoopy over 10,000 time points and recorded 10,000 observations and truncating the accuracy up to 2 digits after the decimal place. There are 21 species (places), which we categories using the simulative Model Checker MC2 by means of a property library, we developed previously [58]; compare Fig. 7.

1
always steady state zero:
P_{≥1}[ G([$x]=0) ]
There are no corresponding variables, i.e. there are no entities which are never active.

2
always steady state above zero:
P_{≥1}[ G(d[$x]=0 & [$x]>0) ]
holds for two entities, which are genes – in this model these are side conditions, they are always active, but never consumed: luxS_gene and pfs_gene.

3
decreasing weakly until steady state zero:
P_{≥1}[ G(d[$x]≤0) & ¬G(d[$x]=0) & ((G(d[ $x] ≤0) & ¬G(d[ $x]=0)) U (G(d[ $x]=0 & [$x]=0))) ]
holds for the two source places, i.e. places without pretransitions, which reflect our assumption in the model of finite supply: Nutrients and Putrescine.

4
increasing weakly, until steady state:
P_{≥1}[ G(d[$x]≥0) & ¬G(d[$x]=0) & ((G(d[$x]≥0) & ¬G(d[$x]=0)) U (G(d[$x]=0))) ]
holds for six entities comprising the sink place AI_2_Out, and two pairs of protein and mRNA which are produced and degraded continuously until production and degradation are balanced LuxS_mRNA, LuxS_prot, Pfs_mRNA, Pfs_prot, and also for Decarb_SAM where the inflow and outflow eventually become balanced.

5
peaks once and falls weakly until steady state zero:
P_{≥1}[ d[$x]>0) & (d[$x]>0 U (G(d[$x])≤0 U (G(d[$x]=0 & [$x]=0)))) ]
holds for seven entities: Adenine, DPD, Homocys, MTA, MTR, SAH, and SRH.

6
falls and rises to peak before falling weakly until steady state zero:
P_{≥1}[ ((d[$x]<0) & (d[$x]<0 U F((d[$x]>0) & (d[$x]>0 U (F((d[$x]<0) & (d[$x]<0 U (G(d[$x]=0) & [$x]=0)))))))) ]
holds for four entities: AI2_intra, Met, SAM, and Spermidine.
Phase 3 model validation
Standard Petri net analysis techniques do not cover inhibitory arcs; so we preclude the go mechanism from the structural analysis. Applying the same approach as we did for the first component relating to Phase 2 confirms again that the transitionbordered version of the Petri net is covered with Tinvariants, and reveals two Pinvariants  a trivial one, comprising QSeBC, and another one comprising the two states (repressed, derepressed) of the lsr genes, i.e., the places lsrGenes_LsrR and lsrGenes. Moreover, the STP holds; but because the net structure is beyond ES we can only structurally deduce that the Petri net is free of dead states; see Additional file 3 for details.
Next, we considered our model under two conditions: (i) unlimited amount of AI_2_Out and (ii) limited amount of AI2_Out. Each condition was simulated in the stochastic setting using the Gillespie [31] algorithm. The simulation was also performed in the hybrid setting using a simulator comprising three components: (1) an ODE solver, we specifically use the CVODE library [59], for the continuous part, (2) Gillespie simulation for the stochastic part, and (3) the synchronisation between the continuous and stochastic net components is done via the improved Hybrid Rejectionbased Stochastic Simulation Algorithm (HRSSA) [60], which combines the accelerated method introduced in [60] with the hybrid rejectionbased stochastic algorithm from [61].
Under the first condition we observed a repeatedly peaking behaviour for most of the metabolites because the go place is repeatedly opening and closing. These metabolites were AI_2_In, AI_2_Out, AI_2_P, LsrABCD, LsrK, LsrR_AI_2_P, lsrGenes, lsrGenes_LsrR and go. Because of its accumulatory nature, Biofilm increases stepwise. See Fig. 8 for the hybrid results The stochastic results are not reported as they are similar for a sufficiently high number of runs.
In contrast, with a limited amount of AI_2_Out no repeated peaks were observed for most of the places because the go place is immediately occupied once and then becomes empty later on.
To explore the repression–depression machinery causing changes between the states of lsrGenes bound or unbound to LsrR, we performed stochastic runs in the limited AI2_out scenario. The onoff nature of the binding is clearly observable in single runs, whereas the nature of the overall state switch driven by the availability of AI2_P is apparant in the average of many stochastic runs, see Fig. 9. The average behaviour can equally be observed in the hybrid setting. The initial repressed state will be reached again as AI2_P approaches zero steady state; this is not observed in the unlimited AI_2_Out setting because AI2_P never approaches a zero steady state.
Validation of the nonspatial combined model
We start with the following observations of some specific aspects of the model behaviour. Increasing the threshold THR of the transition open (in the following briefly called THR_open) postpones the inflow into the biofilm component and allows higher values of AI2_Out to accumulate. When THR of the transition close (in the following briefly called THR_close) is close to 1, i.e. the gate rarely closes, then it takes longer until the AI2_In_Phase3 is processed by the biofilm component, because there is an higher inflow.
When THR_close is close to THR_open (THR_close≤THR_open), then the thresholds strongly influence how much of the AI2 is used by the biofilm component and how much is left to diffuse in space by the transition Diffusion. The experimental evidence [62] is that the biological threshold mechanism is symmetrical in that the opening and closing of the AI2 inward transport channel occur at the same threshold value, which we have reflected by making the values for open and close identical.
We found that the lsr genes eventually reach the repressed state in the steady state using the following values: diffusion parameter kd=0.1,THR_open=THR_close=40, with a simulation time = 15000, see top row in Fig. 10, and compare with Fig. 9. While performing these experiments we noticed that the runtimes to simulate the hybrid model were considerably greater than those for the stochastic simulation (using Gillespie: 95s, hybrid: 3609s), caused by increasing stiffness and discontinuities which are problematic for continuous solvers. Thus in the following we use stochastic simulations. We further reduced the runtime required for stochastic simulation to 13 s, i.e. by a factor of over 7, by employing the approximative deltaleaping algorithm [32, 63]. Although deltaleaping is an approximative approach, the results are sufficiently close to the exact Gillespie algorithm, see second row in Fig. 10.
Furthermore the behaviour of the combined model can be characterised by a combination of the validation results reported in “Phase 2 model validation” and “Phase 3 model – biofilm formation” sections, which reflects the fact that the two components have a simple interface comprising just one species (place).
Validation of the spatial combined model
The validation is carried out by simulating an automatically unfolded stochastic Petri net. The size of the unfolded net obviously depends on the size of the grid, i.e., its resolution, and on the number of bacteria we position on the grid (the exact position does not matter). Due to the regular structure of our coloured model, we are in the favourite position to be able to provide closed formulae specifying the size of the underlying uncoloured model. The coloured Petri net comprises 34 places and 40 transitions. Thus, with N being the number of bacteria and D the size of the square grid, we obtain by unfolding (33N+D^{2}) places and (39N+8D^{2}−12D+4) transitions. For example, with D=21 and N=25, we obtain an unfolded Petri net comprising 1266 places and 4255 transitions, and with D=101 and N=25 we obtain 11,026 places and 81,375 transitions.
We explored the behaviour of the model for different configurations: a single bacterium, and for colonies with 25 bacteria under different sparseness conditions, in two different environment grids (21x21 and 101x101). We simulated each configuration with the stochastic approximative Delta leaping solver for the same time period (100 runs for 20,000 time steps and recording 1000 time points).
In all cases the AI2_Out peaks before the Biofilm reaches a steady state due to the exhaustion of Nutrients and Putrescine in the Phase 2 part of the model.
We found that the single bacterium produced an extremely small amount of biofilm. In all the other configurations, although the production of AI2_Out is the same in all these cases because the number of bacteria is the same, the total amount of AI2_Out remaining in the environment increased as the sparseness increased, because less is absorbed by the colony under sparse conditions.
Further, we observed that the total amount of Biofilm noticeably decreased with sparseness, because AI2_Out diffuses more rapidly away from the bacteria in sparse situations, thus reducing the signal required to trigger the threshold and maintain the production of biofilm; see Tables 7 and 8. We also found that that the values of biofilm at the central positions of a compact colony are higher than those at the edges, see Figs. 11 and 12.
These observations were in concord with our expectations; so we consider the model as validated. To illustrate the spatial behaviour, we provide movies as supplementary material showing the behaviour of AI2_Out and Biofilm over time and space, see Additional file 1 for download instructions.
We conclude our model validation with questioning our stochastic simulation strategy by comparing the results produced by the exact Gillespie algorithm with those generated by the Deltaleaping algorithm, for 10 and 100 runs each. For doing this, we generated four traces of simulation length 20,000 recorded at 1,000 time points for a grid square of size 9 in the 2D scenario, hosting a sparse colony of 9 bacteria distributed over a region with a radius of 2. We use the following formula to calculate the absolute difference between two simulation traces over time:
The result is represented as a heat map, showing the propagation of the absolute difference encoded by colour over the given time period, where time and the trace variables are represented by x and y axes, respectively; see Fig. 13.
The stripes in the middle relate to the biofilm produced by the 25 bacteria positioned in the middle of the grid, the two thinner stripes in the lower half relate to the two genes in the Phase 2 model, which both have lower token numbers, thus stochasticity has an higher effect. Both stripes can be explained by the stochastic effect not levelled out by the small number of runs. We take the close resemblance of the two heat maps in the first row as a confirmation that the delta leaping simulation method sufficiently well approximates the exact stochastic method, while substantially reducing the runtime, as previously reported [32, 58]. All heat maps together illustrate that there is not much difference between 10 and 100 runs, which we take as a justification not to increase the number of runs over which we average.
Conclusions
The spatial model which we have developed and presented in this paper describes E.coli biofilm formation driven by the Autoinducer 2 (AI2). We have incorporated the nonspatial model of AI2 production in Phase 2 of bacterial growth previously published by Li et al. [24], and linked it with our description of biofilm formation in Phase 3. We have described the method by which we have rigorously developed and tested our model in a stepwise manner covering continuous, stochastic, hybrid, and spatial aspects at different levels of abstraction. The spatial representation is achieved using coloured stochastic and coloured hybrid Petri nets, and the behaviour of this version reveals emergent properties not evident in the nonspatial one. For example, we have investigated the model under different configurations ranging from a single bacterium, to colonies under different conditions of sparseness, all contained within environments of different sizes. The model can also be used to investigate the behaviour of models under different host environments, for example finite resources versus those renewed under steady or even cyclical regimes – as occurs in the gut or in oral biofilm. Our results confirm that our model behaves as expected, i.e. that biofilm formation is increased in areas of higher bacterial density. Our research into related work has shown that there is a lack of quantitative experimental data on the effect of population density on quorum sensing driven biofilm formation, and we suggest that this is an open area for exploration by experimentalists.
Although we have explored configurations in 2D in order to keep computations within reasonable time bounds, our approach to encode space is highly flexible, and easily configured for the 1D, 2D or 3D scenario. Our spatial modelling strategy can be equally applied to other problems evolving in time and space.
Reproducibility
We provide all models in their source format, and use only publicly available tools. Modelling and simulation was done with Snoopy [42]. Additionally, we used Marcie [64] and MC2 [49] for simulative model checking. The 2D visualisation was done by help of Octave [65]; all our scripts are available on request.
Abbreviations
 1D/2D/3D space:

One/two/threedimensional space
 AI:

Autoinducer
 AI2:

Autoinducer 2
 CANDL:

Coloured abstract net description language
 CPN:

Continuous Petri nets
 ES net:

Extended simple net
 HPN:

Hybrid Petri nets
 IBD:

Inflammatory bowel disease
 IBS:

Irritable bowel syndrome
 LTL:

Lineartime temporal logic
 ODEs:

Ordinary differential equations
 PDEs:

Partial differential equations
 PLTLc:

Probabilistic lineartime temporal logic with numerical constraints
 QPN:

Qualitative Petri nets
 SBML:

Systems biology markup language
 SPN:

Stochastic Petri nets
 STP:

Siphon trap property
References
 1
Sender R, Fuchs S, Milo R. Revised estimates for the number of human and bacteria cells in the body. PLoS Biol. 2016; 14(8):1002533.
 2
Distrutti E, Monaldi L, Ricci P, Fiorucci S. Gut microbiota role in irritable bowel syndrome: New therapeutic strategies. World J Gastroenterol. 2016; 22(7):2219.
 3
Maslowski KM, Mackay CR. Diet, gut microbiota and immune responses. Nat Immunol. 2011; 12(1):5–9.
 4
Goodrich JK, Waters JL, Poole AC, Sutter JL, Koren O, Blekhman R, Beaumont M, Van Treuren W, Knight R, Bell JT, et al. Human genetics shape the gut microbiome. Cell. 2014; 159(4):789–99.
 5
Bäckhed F, Ley RE, Sonnenburg JL, Peterson DA, Gordon JI. Hostbacterial mutualism in the human intestine. Science. 2005; 307(5717):1915–20.
 6
González JE, Keshavan ND. Messing with bacterial quorum sensing. Microbiol Mol Biol Rev. 2006; 70(4):859–75.
 7
Papenfort KA, Bassler BL. Quorum sensing signalresponse systems in Gramnegative bacteria. Nat Rev Microbiol. 2016; 14(9):576–88.
 8
Elias S, Banin E. Multispecies biofilms: living with friendly neighbors. FEMS Microbiol Rev. 2012; 36(5):990–1004.
 9
Surette MG, Bassler BL. Quorum sensing in escherichia coli and salmonella typhimurium. Proc Natl Acad Sci. 1998; 95(12):7046–50.
 10
De Keersmaecker SC, Sonck K, Vanderleyden J. Let LuxS speak up in AI2 signaling. Trends Microbiol. 2006; 14(3):114–9.
 11
Xavier KB, Bassler BL. Regulation of uptake and processing of the quorumsensing autoinducer AI2 in Escherichia coli. J Bacteriol. 2005; 187(1):238–48.
 12
Rossi E, Cimdins A, Lüthje P, Brauner A, Sjöling Å, Landini P, Römling U. “It’s a gut feeling”–Escherichia coli biofilm formation in the gastrointestinal tract environment. Crit Rev Microbiol. 2018; 44(1):1–30.
 13
Pereira CS, Thompson JA, Xavier KB. Ai2mediated signalling in bacteria. FEMS Microbiol Rev. 2013; 37(2):156–81.
 14
Nadell CD, Xavier JB, Foster KR. The sociobiology of biofilms. FEMS microbiology reviews. 2008; 33(1):206–24.
 15
Li YH, Tian X. Quorum sensing and bacterial social interactions in biofilms. Sensors. 2012; 12(3):2519–38.
 16
PérezVelázquez J, Gölgeli M, GarcíaContreras R. Mathematical modelling of bacterial quorum sensing: a review. Bull Math Biol. 2016; 78(8):1585–639.
 17
Nealson KH, Platt T, Hastings JW. Cellular control of the synthesis and activity of the bacterial luminescent system. J Bacteriol. 1970; 104(1):313–22.
 18
Eberhard A. Inhibition and activation of bacterial luciferase synthesis. J Bacteriol. 1972; 109(3):1101–5.
 19
Waters CM, Bassler BL. Quorum sensing: celltocell communication in bacteria. Annu Rev Cell Dev Biol. 2005; 21:319–46.
 20
Goryachev AB. Understanding bacterial cell cell communication with computational modeling. Chem Rev. 2011; 111(1):238–50.
 21
Weber M, Buceta J. Dynamics of the quorum sensing switch: stochastic and nonstationary effects. BMC Syst Biol. 2013; 7(1):6.
 22
Hense BA, Schuster M. Core principles of bacterial autoinducer systems. Microbiol Mol Biol Rev. 2015; 79(1):153–69.
 23
Ward JP, King J, Koerber A, Croft J, Sockett R, Williams P. Cellsignalling repression in bacterial quorum sensing. Math Med Biol. 2004; 21(3):169–204.
 24
Li J, Wang L, Hashimoto Y, Tsao CY, Wood TK, Valdes JJ, Zafiriou E, Bentley WE. A stochastic model of Escherichia coli AI2 quorum signal circuit reveals alternative synthesis pathways. Mol Syst Biol. 2006; 2(1):67.
 25
Eberl H, Morgenroth E, Noguera D, Picioreanu C, Rittmann B, van Loosdrecht M, Wanner O, Vol. 18. Mathematical Modeling of Biofilms: IWA Publishing; 2006.
 26
Li J, Attila C, Wang L, Wood TK, Valdes JJ, Bentley WE. Quorum sensing in escherichia coli is signaled by ai2/lsrr: effects on small rna and biofilm architecture. J Bacteriol. 2007; 189(16):6011–20.
 27
Novak EA, Shao H, Daep CA, Demuth DR. Autoinducer2 and QseC control biofilm formation and in vivo virulence of Aggregatibacter actinomycetemcomitans. Infect Immun. 2010; 78(7):2919–26.
 28
Janowski S, Kormeier B, Töpel T, Hippe K, Hofestädt R, Willassen N, Friesen R, Rubert S, Borck D, Haugen P, et al. Modeling of celltocell communication processes with petri nets using the example of quorum sensing. In silico Biol. 2010; 10(1, 2):27–48.
 29
Müller J, Kuttler C, Hense BA, Rothballer M, Hartmann A. Cell–cell communication by quorum sensing and dimensionreduction. J Math Biol. 2006; 53(4):672–702.
 30
PérezJiménez MJ, RomeroCampero FJ. P systems, a new computational modelling tool for systems biology In: Priami C, Plotkin G, editors. Transactions on Computational Systems Biology VI. Lecture Notes in Computer Science, vol 4220. Berlin: Springer: 2006. p. 176–97. https://doi.org/10.1007/11880646_8.
 31
Gillespie DT. Exact stochastic simulation of coupled chemical reactions. J Phys Chem. 1977; 81(25):2340–61.
 32
Rohr C. DiscreteTime Leap Method For Stochastic Simulation. Fundamenta Informaticae, Special Issue ’Advances in Biological processes and Petri nets (BioPPN)’. 2018; 160(12):181–98. https://doi.org/10.3233/FI20181680.
 33
Liu F, Blätke M, Heiner M, Yang M. Modelling and simulating reaction – diffusion systems using coloured Petri nets. Comput Biol Med. 2014; 53:297–308. https://doi.org/10.1016/j.compbiomed.2014.07.004.
 34
Balagaddé FK, Song H, Ozaki J, Collins CH, Barnet M, Arnold FH, Quake SR, You L. A synthetic escherichia coli predator–prey ecosystem. Mol Syst Biol. 2008; 4(1):187.
 35
Heiner M, Gilbert D. In: Kristensen L, Petrucci L, (eds).How Might Petri Nets Enhance Your Systems Biology Toolkit. Lecture Notes in Computer Science, vol. 6709: Springer; 2011, pp. 17–37. https://doi.org/10.1007/9783642218347_2.
 36
Chaouiya C. Petri net modelling of biological networks. Brief Bioinform. 2007; 8(4):210–9.
 37
Ajmone Marsan M, Balbo G, Conte G, Donatelli S, Franceschinis G. Modelling with Generalized Stochastic Petri Nets. New York: Wiley; 1995.
 38
Gilbert D, Heiner M. From Petri nets to differential equations  an integrative approach for biochemical network analysis. In: Proc. ICATPN 2006. LNCS 4024, Springer: 2006. p. 181–200. https://doi.org/10.1007/11767589_11.
 39
Breitling R, Gilbert D, Heiner M, Orton R. A structured approach for the engineering of biochemical network models, illustrated for signalling pathways. Brief Bioinform. 2008; 9(5):404–21. https://doi.org/10.1093/bib/bbn026.
 40
Herajy M, Heiner M. Hybrid Representation and Simulation of Stiff Biochemical Networks. J Nonlinear Anal Hybrid Syst. 2012; 6(4):942–59. https://doi.org/10.1016/j.nahs.2012.05.004.
 41
Herajy M, Heiner M. Adaptive and Biosemantics of Continuous Petri Nets: Choosing the Appropriate Interpretation. Fundamenta Informaticae. 2018; 160(1–2):53–80. https://doi.org/10.3233/FI20181674.
 42
Heiner M, Herajy M, Liu F, Rohr C, Schwarick M. Snoopy – a unifying Petri net tool. In: Proc. PETRI NETS 2012. LNCS, vol. 7347. Berlin: Springer: 2012. p. 398–407. https://doi.org/10.1007/9783642311314_22.
 43
Blätke M, Heiner M, Marwan W. Biomodel engineering with Petri nets In: Robeva R, editor. Algebraic and Discrete Mathematical Methods for Modern Biology. Academic Press: 2015. p. 141–92. https://doi.org/10.1016/B9780128012130.000071.
 44
Heiner M, Gilbert D, Donaldson R. In: Bernardo M, Degano P, Zavattaro G, (eds).Petri Nets for Systems and Synthetic Biology. LNCS, vol. 5016.Berlin: Springer; 2008, pp. 215–64. https://doi.org/10.1007/9783540688945_7.
 45
Gilbert D, Heiner M, Liu F, Saunders N. Colouring spacea coloured framework for spatial modelling in systems biology. In: International Conference on Applications and Theory of Petri Nets and Concurrency. Berlin: Springer: 2013. p. 230–49. https://doi.org/10.1007/9783642386978_13.
 46
Gao Q, Gilbert D, Heiner M, Liu F, Maccagnola D, Tree D. Multiscale Modelling and Analysis of Planar Cell Polarity in the Drosophila Wing. IEEE/ACM Trans Comput Biol Bioinform. 2013; 10(2):337–51. https://doi.org/10.1109/TCBB.2012.101.
 47
Pârvu O, Gilbert D, Heiner M, Liu F, Saunders N, Shaw S. Spatialtemporal modelling and analysis of bacterial colonies with phase variable genes. ACM Trans Model Comput Simul (TOMACS). 2015; 25(2):25. https://doi.org/10.1145/2742546.
 48
Liu F, Heiner M, Gilbert D. Coloured Petri nets for multilevel, multiscale, and multidimensional modelling of biological systems. Brief Bioinform. 2017. https://doi.org/10.1093/bib/bbx150. published 03 November 2017.
 49
Donaldson R, Gilbert D. A model checking approach to the parameter estimation of biochemical pathways. In: International Conference on Computational Methods in Systems Biology. Berlin: Springer: 2008. p. 269–87. https://doi.org/10.1007/9783540885627_20.
 50
Pnueli A. The Temporal Semantics of Concurrent Programs. Theor Comput Sci. 1981; 13:45–60.
 51
Baier C. On Algorithmic Verification Methods for Probabilistic Systems. Habilitation thesis, University of Mannheim. 1998.
 52
Fages F, Rizk A. On the analysis of numerical data time series in temporal logic. In: International Conference on Computational Methods in Systems Biology. Berlin: Springer: 2007. p. 48–63. https://doi.org/10.1007/9783540751403_4.
 53
Wang L, Hashimoto Y, Tsao CY, Valdes JJ, Bentley WE. Cyclic AMP (cAMP) and cAMP receptor protein influence both synthesis and uptake of extracellular autoinducer 2 in Escherichia coli. J bacteriol. 2005; 187(6):2066–76.
 54
Quan DN, Tsao CY, Wu HC, Bentley WE. Quorum sensing desynchronization leads to bimodality and patterned behaviors. PLoS Comput Biol. 2016; 12(4):1004781.
 55
Raj A, van Oudenaarden A. Nature, nurture, or chance: stochastic gene expression and its consequences. Cell. 2008; 135(2):216–26.
 56
Starke P. Some properties of timed nets under the earliest firing rule In: G R, editor. European Workshop on Applications and Theory in Petri Nets. Berlin: Springer: 1988. p. 418–32. https://doi.org/10.1007/3540524940_39.
 57
Heiner M, Lehrack S, Gilbert D, Marwan W. Extended stochastic petri nets for modelbased design of wetlab experiments. Trans Comput Syst Biol XI. 2009; 5750:138–63. https://doi.org/10.1007/9783642041860_7.
 58
Gilbert D, Heiner M, Jayaweera Y, Rohr C. Towards dynamic genomescale models. Brief Bioinform. 2017. https://doi.org/10.1093/bib/bbx096.
 59
Sundials. Sundials CVODE – Suite of Nonlinear and Differential/algebraic Equation Solvers. http://computation.llnl.gov/projects/sundials/sundialssoftware. Accessed 03 Mar 2019.
 60
Herajy M, Heiner M. Accelerated simulation of hybrid biological models with quasidisjoint deterministic and stochastic subnets. In: International Workshop on Hybrid Systems Biology. Cham: Springer: 2016. p. 20–38. https://doi.org/10.1007/9783319471518_2.
 61
Marchetti L, Priami C, Thanh V. HRSSA–Efficient hybrid stochastic simulation for spatially homogeneous biochemical reaction networks. J Comput Phys. 2016; 317:301–17.
 62
Goryachev A, Toh D, Lee T. Systems analysis of a quorum sensing network: design constraints imposed by the functional requirements, network topology and kinetic constants. Biosystems. 2006; 83(23):178–87.
 63
Rohr C. Simulative analysis of coloured extended stochastic Petri nets. 2017. PhD thesis, BTU Cottbus, Dep. of CS.
 64
Heiner M, Rohr C, Schwarick M. MARCIEmodel checking and reachability analysis done efficiently In: Colom J, Desel J, editors. International Conference on Applications and Theory of Petri Nets and Concurrency. Berlin: Springer: 2013. p. 389–99. https://doi.org/10.1007/9783642386978_21.
 65
Eaton J, Bateman D, Hauberg S, Wehbring R. GNU Octave Version 4.0.0 Manual: a Highlevel Interactive Language for Numerical Computations; 2015. http://www.gnu.org/software/octave/doc/interpreter. Accessed 03 Mar 2019.
 66
Correia PJBC. Regulation of the mechanism of interference with quorum sensing in escherichia coli. PhD thesis, Universidade De Lisboa, Faculdade De Ciencias, Departamento De Biologia Vegetal. 2011.
Acknowledgements
We would like to thank Christian Rohr for his assistance during model construction.
Funding
The open access fee has been covered by Brunel University London.
Availability of data and materials
All models in their source format and some video clips visualising analysis results can be found at https://wwwdssz.informatik.tucottbus.de/DSSZ/Software/Examples. Additional files (pdf) provide a Primer for the use of coloured Petri nets in Snoopy, a complete specification of the coloured SPN developed, and some further model validation results.
About this supplement
This article has been published as part of BMC Bioinformatics Volume 20 Supplement 4, 2019: Methods, tools and platforms for Personalized Medicine in the Big Data Era (NETTAB 2017). The full contents of the supplement are available online at https://bmcbioinformatics.biomedcentral.com/articles/supplements/volume20supplement4.
Author information
Author notes
Affiliations
Contributions
The Petri net models have been jointly developed by DG, MH and LG. Model validation by use of model checking has been done by DG and MH, JC took care of 2D visualisations. All authors contributed in writing the paper; all authors have read and approved the final manuscript.
Corresponding author
Correspondence to David Gilbert.
Ethics declarations
Ethics approval and consent to participate
Not applicable.
Consent for publication
Not applicable.
Competing interests
The authors declare that they have no competing interests.
Publisher’s Note
Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.
Additional information
Authors’ information
David Gilbert is a Professor in the Department of Computer Science, Brunel University London. His research interests include Bioinformatics, Systems Biology, Synthetic Biology, multilevel/multiscale modelling, model checking and computational methods for design of biological systems.
Monika heiner is a Professor in the Department of Computer Science, Brandenburg University of Technology CottbusSenftenberg. Her research interests include multilevel/multiscale modelling and analysis of technical as well as biochemical networks using qualitative and quantitative Petri nets, model checking and simulation techniques.
Leila Ghanbar is a PhD student in the Department of Computer Science, Brunel University London, conducting her PhD project in multiscale BioModel Engineering, supervised by DG.
Jacek Chodak is a research assistant in the Department of Computer Science, Brandenburg University of Technology CottbusSenftenberg, conducting his PhD project in efficient simulation techniques, supervised by MH.
Additional files
Additional file 1
Colouring space. A brief primer illustrating how to deal with coloured Petri nets in Snoopy by means of diffusion in 3D. (PDF 567 KB)
Additional file 2
Complete CANDL specification of the case study. The complete listing of the CANDL specification for the coloured SPN given in Fig. 6. (PDF 184 KB)
Additional file 3
Supplementary material for model validation. Additional explanations and figures illustrating various aspects of the model validation performed. (PDF 3614 KB)
Rights and permissions
Open Access This article is distributed under the terms of the Creative Commons Attribution 4.0 International License (http://creativecommons.org/licenses/by/4.0/), which permits unrestricted use, distribution, and reproduction in any medium, provided you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license, and indicate if changes were made. The Creative Commons Public Domain Dedication waiver (http://creativecommons.org/publicdomain/zero/1.0/) applies to the data made available in this article, unless otherwise stated.
About this article
Published
DOI
Keywords
 Quorum sensing
 Biofilm formation
 Coloured hybrid Petri nets
 Coloured stochastic Petri nets
 Diffusion in 3D space
 Simulative model checking