Logicalcontinuous modelling of posttranslationally regulated bistability of curli fiber expression in Escherichia coli
 Kaveh Pouran Yousef^{1}Email author,
 Adam Streck^{1},
 Christof Schütte^{1},
 Heike Siebert^{1},
 Regine Hengge^{2} and
 Max von Kleist^{1}
Received: 9 December 2014
Accepted: 29 June 2015
Published: 23 July 2015
Abstract
Background
Bacteria have developed a repertoire of signalling mechanisms that enable adaptive responses to fluctuating environmental conditions. The formation of biofilm, for example, allows persisting in times of external stresses, e.g. induced by antibiotics or a lack of nutrients. Adhesive curli fibers, the major extracellular matrix components in Escherichia coli biofilms, exhibit heterogeneous expression in isogenic cells exposed to identical external conditions. The dynamical mechanisms underlying this heterogeneity remain poorly understood. In this work, we elucidate the potential role of posttranslational bistability as a source for this heterogeneity.
Results
We introduce a structured modelling workflow combining logical network topology analysis with timecontinuous deterministic and stochastic modelling. The aim is to evaluate the topological structure of the underlying signalling network and to identify and analyse model parameterisations that satisfy observations from a set of genetic knockout experiments. Our work supports the hypothesis that the phenotypic heterogeneity of curli expression in biofilm cells is induced by bistable regulation at the posttranslational level. Stochastic modelling suggests diverse noiseinduced switching behaviours between the stable states, depending on the expression levels of the cdiGMPproducing (diguanylate cyclases, DGCs) and degrading (phosphodiesterases, PDEs) enzymes and reveals the quantitative difference in stable cdiGMP levels between distinct phenotypes. The most dominant type of behaviour is characterised by a fast switching from curlioff to curlion with a slow switching in the reverse direction and the second most dominant type is a longterm differentiation into curlion or curlioff cells. This behaviour may implicate an intrinsic feature of the system allowing for a fast adaptive response (curlion) versus a slow transition to the curlioff state, in line with experimental observations.
Conclusion
The combination of logical and continuous modelling enables a thorough analysis of different determinants of bistable regulation, i.e. network topology and biochemical kinetics, and allows for an incorporation of experimental data from heterogeneous sources. Our approach yields a mechanistic explanation for the phenotypic heterogeneity of curli fiber expression. Furthermore, the presented work provides a detailed insight into the interactions between the multiple DGC and PDEtype enzymes and the role of cdiGMP in dynamical regulation of cellular decisions.
Keywords
Logical modelling Stochastic modelling Bistability Phenotypic heterogeneity Biofilm CdiGMP Escherichia coliBackground
The ability to adapt to external inputs is a crucial property of cellular systems [1]. Bacteria, for example, can outlast antibiotic stress by forming biofilm colonies [2], or they may produce bacteriocins when sensing nutrient competition [3]. This form of acute adaptation is often encoded in the bacterial signal transduction network, which enables multiple stable states and may thus give rise to phenotypic heterogeneity [4–6].
In Escherichia coli, the regulatory system of amyloid curli fibers is a central component of the stressinduced biofilm formation under the control of the master regulator σ ^{s} [7–9]. While heterogeneous allornothing expression of the curli protein CsgB in isogenic wild type cells under identical environmental conditions was shown in [10, 11], subsequent work [7] allowed to infer interactions of the molecular components of the underlying signalling network. However, the ability of the suggested network to produce bistable behaviour (curli on/off) and the influence of individual molecular components on singlecell dynamics still remains to be elucidated.
Rigorous mathematical modelling may further our mechanistic understanding of this system. Two different methodological views for analysing bistability in signalling networks have been established in previous studies. Logical modelling approaches assign discrete states to different activity levels of network components and represent reactions by logical functions. The advantage of this methodology is that it allows to test a large amount of alternative model topologies using a limited amount of data, or semiquantitative data such as promoter activities in genetic knockout strains or immunoblotting data [12, 13].
Particularly regarding asymptotic behaviour, logical analysis has been shown to be a useful tool for uncovering fundamental characteristics and functionalities of biological systems, (see e.g. [14–16]). Consequently, it is a wellsuited approach for investigating bistable systems, where there is evidence that qualitative properties, in particular the network topology, are crucial determinants [17, 18]. However, the low detail resolution of such models only allows for a preliminary understanding, and the biological interpretation of the results is not always straightforward.
Timecontinuous models (ODEbased or stochastic) entail more biological detail and model parameters are readily interpretable. In this approach, besides topology, the signalling network is defined by wellestablished reaction rate kinetics. Among others, the approaches for analysing bistability in continuous models are based on a feasibility analysis of unstable states [19] and the chemical reaction network theory [20], as e.g. applied to study bistable regulation of split histidine kinases in twocomponent signalling networks [21]. In analytically less tractable models, heuristic parameter search may be used, such as Monte Carlo [22] or the genetic algorithm [23]. However, an exhaustive analysis of a large space of alternative model topologies and parameterisations is often infeasible, requiring a confinement to certain regions of the model and parameter space. Thus, a drawback of the continuous modelling approach lies in the containment to a network topology and the requirement of sufficient information about the underlying biochemical reactions.
In addition to the deterministic modelling approach, timecontinuous stochastic modelling [24–26] may deliver further insights into the system behaviour such as the probability of the cells to be in either of the two stable states and the noiseinduced switching dynamics between the stable states [27, 28]. The high computational effort for analysing stochastic models precludes, with only few exceptions, e.g. [29], from addressing inverse problems (model and parameter estimation) at all.
Methods
Regulatory network of curli expression
The second messenger molecule bis(3’5’)cyclic dimeric guanosine monophosphate (cdiGMP) is the central component within the regulatory signalling network of curli fimbriae during the stationary phase of the bacterial population growth cycle and upon induction of various stress conditions. CdiGMP positively contributes to the expression of CsgD, the key biofilm regulator, activating the expression of the curli regulon (csgBAC), with CsgB and CsgA proteins constituting the curli fiber [9, 35]. As shown in Fig. 1 b, the synthesis and degradation of cdiGMP within this network is maintained by two diguanylate cyclase (DGC)/phosphodiesterase (PDE) pairs: YegE/YhjH (module I) and YdaM/YciR (module II).
While DGCtype enzymes synthesise cdiGMP from two GTP substrates, PDEtype enzymes are responsible for the degradation of cdiGMP. As shown in Fig. 1 b, besides the degradation function of cdiGMP, YciR exhibits a second activity: the inhibition of YdaM and MlrA. Importantly, YdaM is the key activator of the transcription factor MlrA, which directly activates transcription at the csgD promoter region. A sufficiently high amount of cdiGMP prevents YciR from inhibiting YdaM, since YciR binds cdiGMP and starts to degrade it, acting as a trigger enzyme [7]. In turn, unbound (active) YdaM activates MlrA and induces the signalling cascade leading to the expression of the curli operon csgBAC.
Previously, the activity of the csgB promoter was measured in all possible single, double, triple and quadruple knockout strains of the set of genes yegE, yhjH, ydaM and yciR [7]. As depicted in Fig. 1 a, despite 15 different genetic backgrounds, only four main expression levels of the csgB gene were observed. E. coli strains containing a ydaMknockout exhibit a completely absent expression of the curli gene csgB (hyperrepressed) or a very low (basal) level expression, as compared to wild type (Fig. 1 a). Furthermore, strains containing a yciRknockout reveal an increased (hyperactivated) expression of csgB if ydaM is present. This knockout background is “blind” with respect to knockouts of module I genes (yegE and yhjH). Genetic and biochemical analyses allowed Lindenberg et al. [7] to derive a model of the regulatory network, which is depicted in Fig. 1 b. In this model different functional states are assigned to YciR (I/II) and YdaM (active/inactive).
In order to find evidence for the bistability of the curli signalling system we assigned each of the four different expression levels observed in the genetic knockout experiments at population level (Fig. 1 a) a certain singlecell phenotype. To this end, we used the complementary singlecell measurements of the wild type expression of the CsgB protein in Serra et al. [11]. These singlecell data suggest that the intermediate expression level of the csgB gene at population level corresponds to a heterogeneous mix of curlion and curlioff cells in the bacterial population. This indicates that the basal and the hyperrepressed level of csgB reflects a mixture of cellular phenotypes where the fraction of curlion cells is significantly lower as compared to the wild type, or it becomes undetectable. In contrast, the hyperactivated phenotype may be due to two different types of singlecell expression levels. Either all cells in the population are in the curlion mode, or there is a subset of cells in the curlioff mode but their frequency is significantly lower than the frequency of curlioff cells in the wild type population.
Derivation of the logical model
In order to set up the logical model of the curli regulation system, we described the involved components in terms of logical variables whose values were interpreted as functional states of the components, e.g. active vs. inactive states. For deriving the model we applied a timescale separation. Thus we assumed that the total levels of involved proteins do not significantly change in contrast to their activity states. We based our assumption on results indicating that protein transcription and translation are significantly slower than posttranslational interaction dynamics [36].
Having described the components by means of logical states, logical rules were added to characterise the principles of mutual regulation by the components. In particular, all activating or inhibitory effects indicated by experimental data were incorporated, resulting in the logical model depicted in Fig. 1 c. Most of the components were assigned two states (on or off) corresponding to the situation where the component has either an effect on the system or not. An exception was made for the protein YciR, which was assigned a knockout state (offstate) and two different observable activities that are mutually exclusive (PDE activity and YdaM/MlrAinhibition activity) [7]. Note that YdaM has also two different activity levels. However, since one of them does not have any observable effect, we identified it with the offstate.
In addition to the logical variables, we described the semantics of the model via the intertwined regulatory effects. Thus, a regulation is effective if the respective regulator is on. In the case of YciR we distinguished between alternative regulations depending on their corresponding effect. A regulation is effective only when YciR occurs in the respective state, including the positive and negative effects (Fig. 1 c). Most of the included regulations correspond to the experimentally derived network topology in Fig. 1 b, with a few exceptions. The first difference is given by the negative feedback induced by the allosteric product inhibition of cdiGMP on the catalytic activity of YegE. When molecular numbers are considered, product inhibition (PI) gives rise to an upper limit on the synthesis rate and thus contributes to setting up a homeostatic steadystate level of cdiGMP [37]. However the semantics of this negative feedback loop do not transfer to the discrete logical model i.e. it is not possible for cdiGMP to completely inhibit YegE. Since in the logical model, only one discrete value is used as an abstraction for the concentrations of cdiGMP that are considered as the onstate, this negative feedback inhibition has no observable logical effect. Therefore we eliminated the corresponding edge from the logical interaction graph (see Fig. 1 c).
Furthermore, the DGC/PDEpair YhjH and YegE constitute the inputs of the network, which means that they maintain the values that they were initially set to. Finally, we used the system property that the expression of csgB is induced if its transcription factor CsgD is expressed. This is the case if MlrA is on, since it is in turn the functional activator of csgD transcription (Fig. 1 b). Therefore, in order to reduce the model, we removed csgD and used the boolean component MlrA as the output component of the network. Thus, the onstate of MlrA represents the induction of curli expression in the logicalmodel.
At this point we incorporated a sufficient amount of information into the model for deriving all regulatory functions except the functions of cdiGMP and MlrA. These two components are influenced by multiple regulators, giving rise to a set of alternative logical functions describing their effect on cdiGMP or MlrA (see Additional file 1 for details). After resolving the constraints derived from the interaction effects, we identified 114 possible regulatory functions for cdiGMP and 2 for MlrA. Finally, by generating all possible combinations of the functions of the two components, we obtained 228 alternative models for the whole network.
Formalisation of the experimental data

The intermediate (bistable) phenotype was related to the ability to reach a stable state with MlrA on and a stable state with MlrA off in the logical modelling framework.

In the hyperrepressed phenotype no curli is being expressed, therefore we required that the stable state with MlrA off is reachable, while the onstate is not.

The hyperactivated phenotype exhibits strong expression of curli and therefore we stated that the stable state with MlrA on must be reachable. However, unlike the hyperrepressed phenotype, we did not prohibit the inactive stable state, since we could not rule out that this phenotype is generated by a mix of curlion and curlioff cells.

Lastly, the basal phenotype exhibits only little, but still measurable expression of curli. However, very weak expression of curli was also observed in the hyperrepressed phenotype (see Fig. 1 a). It would therefore be possible to identify the basal with the hyperrepressed phenotype, but this constraint could be spurious. We have therefore decided to not include the basal phenotype in the main analysis.
However, we also tested the validity of the resulting constraint in a second, less conservative analysis step. Here, we assumed that the basal phenotype is equivalent to the hyperrepressed phenotype, see Additional file 1 for details. Since the measurements were conducted in the early stationary phase of the growth cycle (see Fig. 1 a), we expected the input components YegE and YhjH to be on unless knocked out. Using the stability constraints resulting from experimental data (the knockout data), we formulated 15 properties (bold indices in Table S1 in Additional file 1) that a valid model should fulfil (using the conservative analysis step) and 32 properties using the less conservative analysis step (all rows in Table S1 in Additional file 1). We used these properties to test the entire set of alternative logical models. Subsequently, we identified a valid reduced model topology, as well as parameter constraints that were passed on to the continuous modelling step, as exemplified in Additional file 1 and shown in the “Results” section.
Derivation of the continuous reactionrate model
Based on the results from the logical modelling step (see Additional file 1), we included three system components as dynamical variables into the continuous model: the amount of cdiGMP (x _{1}), the amount of YciR molecules in its role as a YdaM inhibitor (x _{2}) (i.e. YciR II in Fig. 1 b) and the amount of active YdaM molecules (x _{3}). Due to the timescale separation assumption, transcription and translation were assumed to be negligible and thus the total amount of YciR and YdaM was assumed to be constant (parameters YciRtot and YdaMtot, respectively). Two further system components, YegE and YhjH, were included as parts of the maximum catalytic velocity parameters (V _{max1} and V _{max2}, respectively). We did not include MlrA into the model since its activity reflects the activity of YdaM. This was supported by the results of the logical modelling step indicating that MlrA does not influence the interactions between module I and module II proteins (we assumed that YdaM and MlrA do not compete for YciR, as discussed later). Details of the derivation of reaction rates of the timecontinuous model are given in the Additional file 2.
According to the network wiring validated in the logical model checking step (see “Results”), the final continuous model was composed of interactions given either by the production or degradation reactions of cdiGMP. Based on the elementary catalytic reactions, we derived MichaelisMenten rates by applying the quasisteadystate assumption (see Additional file 2 for details and derivations).
Biochemical rates of molecular interactions involved in the curli regulation network. Reaction rates were derived according to catalytic properties of cdiGMP regulation and proteinprotein interactions identified in [7]. Since the focus was on posttranslational dynamics, the protein expression levels of all DGCs and PDEs were held constant using the parameters V _{max1} (YegE), V _{max2} (YhjH), YdaMtot (YdaM) and YciRtot (YciR)
Rate function  Description of the reaction 

\(V_{1} = \frac {V_{\max 1}}{1 + x_{1}/K_{\mathrm {i}}^{\text {YegE}}}\)  Synthesis of cdiGMP by YegE 
\(V_{2} = \frac {V_{\max 2} x_{1}}{x_{1} + K_{\mathrm {m}}^{\text {YhjH}}}\)  Degradation of cdiGMP by YhjH 
\(V_{3} =(\textit {YciRtot} x_{2})\frac {k_{\text {YciRact}} x_{1}}{x_{1} + K_{\mathrm {m}}^{\text {YciR}}}\)  Degradation of cdiGMP by YciR 
\(V_{4} = \frac {k_{\text {YdaMact}} {x_{3}^{n}}}{\left (K_{\mathrm {d}_{\text {polymer}}}^{\text {YdaM}} \right)^{n} + {x_{3}^{n}}}\)  Synthesis of cdiGMP by YdaM (no PI) 
\(V_{5} = k_{\text {YciRde}} x_{2} \frac {x_{1}}{x_{1} + K_{\mathrm {d}}^{\text {YciR}}}\)  Transition from YciR II (YdaM inhibiting) 
to YciR I (PDE activity) due to cdiGMP  
binding  
V _{6}=c _{6}(Y c i R t o t−x _{2})  Transition YciR I → YciR II 
\(V_{7} = k_{\text {YdaMde}} x_{3} \frac {x_{2}}{x_{2} + K_{\mathrm {d}}^{\text {YdaM}}}\)  Inhibition of YdaM due to binding of YciR II 
V _{8}=c _{8}(YdaMtot−x _{3})  Reactivation of YdaM 
where the system variables x _{1}, x _{2}, and x _{3} denote the molecular concentrations of cdiGMP, YciR in the YdaMinhibition state and active YdaM, respectively. The reaction rate functions are stated in Table 1.
Substituting x _{2} and x _{3} with \(x_{2}^{\text {ss}}\) and \(x_{3}^{\text {ss}}\) in the first equation of the ODE system (1) yields a onedimensional ODE. Its roots indicate the levels of cdiGMP \(x_{1}^{\text {ss}}\), which, in conjunction with the steady state Eqs. 2 and (3), give rise to the fixed points \(\mathbf {x}^{\text {ss}} = \left \{x_{1}^{\text {ss}},x_{2}^{\text {ss}}, x_{3}^{\text {ss}}\right \}\) of the entire ODE system. A fixed point x ^{ss} is a stable steady state if all eigenvalues of the Jacobian matrix of the ODE system (1) at x ^{ss} are real and negative.
Parameter identification
Ranges of the uniform proposal distribution used for parameter sampling. Kinetic rate parameters were sampled using a uniform proposal distribution with ranges as depicted in the table. Physiologically meaningful boundaries were chosen if no specific information was available from previous studies (last three rows). In the case where the parameter distribution was posthoc trimmed, the effective range is shown and the original sampling range is given in brackets
Parameter names  Unit  Lower bound  Upper bound  References 

\(K_{\text {m}}^{\text {YhjH}}\)  molecules  100  3000  
\(K_{\text {m}}^{\text {YciR}}\)  molecules  100  1000 (3000)  
\(K_{\text {d}}^{\text {YciR}}\)  molecules  100  3000  
\(K_{\text {i}}^{\text {YegE}}\)  molecules  0  2000  [37] 
V _{max1}, V _{max2}, k _{YdaMact}, k _{YciRact}  molecules /s  0  2000  
\(K_{\text {d}}^{\text {YdaM}}\)  molecules  100  3000  physiol. range 
YciRtot, YdaMtot, \(K_{\text {d}_{\text {polymer}}}^{\text {YdaM}} \)  molecules  0  2000  physiol. range 
k _{YciRde}, k _{YdaMde}  molecules /s  0  2000  physiol. range 
In order to reduce the parameter search space, we applied a set of constraints. Firstly, we required that the proposal values of the parameters YciRtot and YdaMtot are equal, as suggested by protein assays performed by Lindenberg et al. [7]. Furthermore, we included a constraint resulting from the logical modelling step requiring that the catalytic activity of YciR is weaker than the catalytic activity of YdaM, i.e. k _{YciRact}≤k _{YdaMact} (see Additional file 1 for derivation).

the wild type model and the double knockout yegE/yhjH exhibit two stable states,

the singlegene knockout mutant lacking yhjH either exhibits one or two stable states,

the singlegene knockout mutant lacking yegE exhibits one stable state.
Note, that the remaining knockout data in Fig. 1 a were not included for kinetic parameter sampling since they contain genetic knockouts of ydaM or yciR. Setting the amount of each of these two components to zero disrupts the structure of the ODE (1) since these two components are described by system variables (as opposed to the activities of YegE and YhjH, described by model parameters). The parameters were estimated using a fixed Hill constant n=4. This is supported by the tendency of YdaM to build tetramers invitro [7], in line with previous work suggesting multimerisation as a possible source for ultrasensitivity [40].
In the second step we conducted a rejection sampling with a multivariate Gaussian proposal distribution by accepting the sampled parameter sets if all stability constraints from the first step were fulfilled. The mean μ of the proposal distribution was given by the parameter sets identified in the first (uniform) sampling step and the standard deviation σ was chosen sufficiently large to cover a large sampling space but minimise the probability of negative sampling proposals (coefficient of variation σ/μ=0.25). In the case of the parameter \(K_{\text {m}}^{\text {YciR}}\) all sampled parameter sets, where the range of 1000 molecules was exceeded, were posthoc rejected due to experimental results suggesting this upper bound [41]. Additionally, we posthoc selected only those parameter sets which fulfilled the validity criterion for the MichaelisMenten approximation in deterministic and stochastic models (see Additional file 2, Eq. (S12)).
Chemical Master Equation of the curli regulation system
The stochastic model of curli regulation and the validity of the translation of continuous rates into stochastic propensities are discussed in Additional file 2. In order to solve the corresponding Chemical Master Equation and to compute switching probabilities between the stable states, we generated a sufficiently large amount of sample trajectories using the Stochastic Simulation Algorithm [42]. For computing the probability of stable states (stationary distribution) we used sufficiently large simulation times ensuring the equilibration of the system. For instance, a simulation time of 100s for the parameter set used for generating Fig. 4 was sufficient, as shown in Additional file 5: Figure S3.
where ε is a neighbourhood parameter that we set to 0.1 and N denotes the number of system components (3 components in the case of our model). The stable states X ^{ss1} and X ^{ss2} were computed as the steady states of the deterministic reactionrate model (1) and confirmed by identifying the modes of the solution of the CME after a sufficiently long stochastic simulation time.
Results
We applied logical model checking in order to assess the potential of the experimentally derived wild type interaction network (Fig. 1 b) to exhibit bistable dynamics and at the same time to fulfil the stability constraints from the various genetic knockout strains (Fig. 1 a).
Following the logical model checking step, we then incorporated currently available information on biochemical reactions in order to derive a more detailed reactionrate model of curli regulation.
Logical analysis: topology and parameters
As a first analysis step, we tested the consistency of the candidate models, constructed in the “Methods” section, with the 15 constraints (Table S1 in Additional file 1, rows indicated with a bold index) derived from the experimental observations (Fig. 1 a). As a result, we found that 10 out of 228 possible models were in agreement with all constraints (see Table S2 in Additional file 1). Using a comparative analysis of this set of 10 feasible models we obtained a number of results yielding meaningful insights into signalling mechanisms of curli regulation and supporting subsequent continuous modelling. In the following, we first describe the essential features of the system resulting from common properties of the feasible models and afterwards we present the consequences for the continuous modelling step.
Apart from yielding insights into the modelled system, the logical approach could also be used to validate crucial assumptions underlying the continuous reactionrate model: Firstly and most importantly, we validated the network structure underlying the continuous reactionrate model. To this end, as described above, we identified a nonempty set of logical models fulfilling the experimental constraints and exhibiting bistability. Furthermore, by inspecting the set of 10 feasible models, we observed that despite of two different choices for the regulation of MlrA, the realisation of the core circuit YdaM—YciR—cdiGMP is independent of the choice of the regulatory function of MlrA. Furthermore, MlrA is an output component and thus does not influence any other components in the network, allowing for model reduction by completely removing MlrA in the continuous modelling step. This gave rise to a core model that generates bistability, containing only cdiGMP, YdaM and YciR as dynamic variables.
In a second step, we asked if, besides the already incorporated (conservative) constraints, the network topology also fulfils the constraints related to the less conservative interpretation of the basal expression phenotype, as e.g. observed in the mutant strain with a yegEknockout (Fig. 1 a, see “Methods” for a description of model constraints). To this end, we made the assumption that the basal phenotype is equivalent with the hyperrepressed phenotype. This assumption is partly supported by the knockout data, which suggest that these two phenotypes are indeed equivalent w.r.t. to the core network (YdaM, YciR and cdiGMP) and only differ by the basal YdaMindependent activity of MlrA inhibited by the activity of YciR (e.g. compare csgB expression in the mutant strain lacking ydaM and the doubleknockout mutant strain lacking ydaM and yciR). Overall we generated two constraints from each knockout experiment in Fig. 1 a, giving rise to 32 constraints (Additional file 1: Table S1). As a result, we identified one valid model fulfilling all constraints (including the less conservative constraints obtained by identifying the basal and hyperrepressed phenotype, see first row in Table S2, Additional file 1).
Finally, we were able to translate our findings on the regulation of cdiGMP into a parameter constraint for the ODE model. Our result stating that the activating effect of YdaM overcomes the inhibitory influence of YciR under certain conditions yields the kinetic parameter inequality k _{YciRact}<k _{YdaMact} in the continuous setting (see Additional file 1 for details).
In summary, we were able to transfer four results into the continuous modelling step. Firstly, the network structure postulated in Fig. 1 c, is a valid basis for a continuous reactionrate model of curli regulation and supports the experimentally derived network topology. Secondly, the core model giving rise to bistable behaviour is based on the interaction between YdaM, YciR and cdiGMP (with YegE and YhjH as model inputs). Thirdly, we have shown that the network topology also fulfils the constraints resulting from the equivalence assumption of the basal and hyperrepressed phenotypes. Thus, this assumption can be used as a constraint in the continuous modelling step. Finally, we obtained a parameter constraint for the catalytic activity of YciR in comparison with the activity of DGC enzyme YdaM.
Reactionrate model of interaction between the DGC/PDE modules
In the subsequent step, we analysed the dynamics of the system after translating the logical network into a continuous reaction rate model. We then determined the stable states of the system (see “Methods”). Since in the logical modelling step the activity of MlrA was shown to almost completely reflect the activity of YdaM, we only focused on interactions of the module I and module II proteins, upstream of MlrA, as the regulatory core of the network (Fig. 1 d and “Methods”). Note, that this model reduction is valid under the assumption that MlrA and YdaM can simultaneously bind YciR i.e. there is no competition (see “Discussion”). We used the bistability of the wild type strain and the stability properties of the singlegene knockout strains lacking either yegE or yhjH (Fig. 1 a) as acceptance criteria for parameter inference (see “Methods”).
Bifurcation analysis and comparison to genetic knockout experiments
Previously, it was shown that proteinprotein interactions between the module II proteins YciR and YdaM induce a switching point in this system requiring a sufficient amount of cdiGMP in order to relieve the inhibition of YdaM by YciR and enable downstream signal transduction [7]. In the following, we analysed the potential role of the regulatory parameters of cdiGMP in generating bistability. Since experimental results suggest that YdaM is the main component transmitting the signal further downstream to MlrA, we used the amount of active YdaM molecules as the readout parameter of our analysis. Note that the stable states with a high amount of YdaM molecules are characterised by a high level of cdiGMP and a low level of YciR II and vice versa, in line with the stable states of the logical model (Fig. 2). By construction of the rejection sampling, all identified parameter sets shared equivalent stability properties of the wild type and genetic knockout strains. Therefore, in order to analyse the dependence of the stable states on model parameters we selected a representative parameter set from the sampled distribution (see Additional file 3, parameter set 1 listed in the first data row of the table).
In the third experiment we varied the complex parameter c _{6}/k _{YciRde}. This parameter is inversely proportional to the binding affinity of cdiGMP to YciR and it has an opposite bifurcation behaviour as compared to V _{max1}. Values of c _{6}/k _{YciRde} below 7.92·10^{−3}, i.e. a high affinity of cdiGMP to YciR, give rise to a single steady state with a high level of active YdaM molecules. Between 7.92·10^{−3} and 8.39·10^{−3} the system becomes bistable with two different levels of active YdaM molecules, which corresponds to the wild type dynamics. The corresponding stochastic simulations indicate that within this parameter region the probability of the system is divided between the two stable states, which is in agreement with the bistable dynamics of csgB expression. Finally, at values higher than 8.39·10^{−3}, i.e. a further decrease of the binding affinity of cdiGMP to YciR, the system becomes monostable again, exhibiting low level of active YdaM molecules (curlioff). This is supported by experimental data indicating a basal expression of csgB in mutant strains where the EAL domain was mutated to an AAL domain, thus diminishing the binding affinity of cdiGMP to YciR [7]. The model suggests that the observed basal expression of curli might be a result of a single steady state with a low level of active YdaM proteins. Note that despite individual differences in the exact parameter values V _{max1}, V _{max2} and c _{6}/k _{YciRde} between distinct feasible parameter sets at bifurcation points, the qualitative dynamics described above can be generalised to all identified parameter sets shown in Fig. 3.
Dynamics of stochastic switching between different activity levels of YdaM
Although the identified parameter sets are qualitatively equivalent with respect to their stability properties (i.e. number of stable states in the wild type, yegE and yhjHknockout mutant strains and bifurcation behaviour), the equilibrium probabilities of each stable state and the dynamics of switching between the stable states might differ significantly. Potential spontaneous switching between the stable states is an important property of stochastic multistable systems. By generating sample trajectories of the stochastic model one can show that on very short time scales the probability of each stable state is fully determined by the initial level of the key components in the system, while after sufficiently long simulation times(∼100s) the switching dynamics between the stable states marginalises the effect of initial molecular levels by inducing equilibration (see Additional file 5: Figure S3).
In addition, we computed the time to switch between the two stable states. A summary plot for all simulations (all parameter sets) is shown in Fig. 5 b, indicating that the time required to switch from the low to the high level of active YdaM (switch from curlioff to curlion) was significantly shorter than for the reverse direction. When starting from the curlioff (low YdaM activity) state, a switch was observed in 93 % of all model parameterisations within 100 seconds. However, only slightly over 41 % exhibited a switch from the high YdaM activity state to the low YdaM activity state after this period (Fig. 5 b). This difference in the times for switching in different directions is in line with the overrepresentation of parameter sets with a high probability of switching from the low to the high YdaM activity state and low probability of switching in the reverse direction, as discussed above (Fig. 5 a).
This suggests that changes of expression and activity levels of the DGC and PDE enzymes due to e.g. transcriptional regulation might modulate the responsiveness of the system.
Overall, besides the potential of the curli regulation system to induce bistable dynamics, our stochastic modelling results indicate alternative feasible scenarios of switching between the stable states. At the one extreme, the cells may be able to rapidly switch their phenotype and at the other extreme they may differentiate into two stable phenotypes. A difference between these two scenarios may be given by distinct expression levels of the involved proteins which regulate the energetic barrier separating the two stable states. It should be noted that certain reactions of the continuous model do not allow to distinguish between the impact of individual parameters as, for instance, the catalytic velocity of cdiGMP degradation is determined by the product of two parameters: YciRtot and k _{YciRact}. Therefore, additional kinetic data may help in future to further dissect the impact of reaction rate parameters on bistability and switching dynamics of curli expression.
Discussion
Multistability is a recurring mechanism in biological systems, enabling cellular heterogeneity and differentiation at the phenotypic level [6]. In order to understand the principles of its regulation, a thorough combination of wetlab experimentation and mathematical modelling is required. Many computational studies of bistable systems focused either on the network topology [17, 18] or on reaction kinetics [19, 21]. Here, we combined these two approaches. This enabled us to deal with heterogeneous data sources, such as expression levels in various knockout strains and kinetic rate parameters from invitro studies.
While multistability in signalling networks may be attributed to transcriptional and translational regulation as well as posttranslational dynamics, in this study, we focused on molecular interactions at the posttranslational level.
In our modelling pipeline we utilised three modelling frameworks in sequence to fully exploit the available data and investigate different aspects of the system on appropriate levels of abstraction. Insights from each modelling and analysis step were passed on in the pipeline and constituted, together with additional biological information accessible within the next higher detailresolving formalism, the foundation for more elaborate models.
In a first step we used a constraintbased logical modelling approach enabling us to extract essential model characteristics needed to reproduce all available experimental observations. More precisely, we identified a set of logical models consistent with the suggested network topology that generates bistability in the wild type strain and fulfils all constraints given by the genetic data. Furthermore, the logical modelling step revealed that the bistable distribution of the molecular pool of active YdaM molecules (not bound by YciR) is the key determinant of phenotypic heterogeneity of curli expression observed in singlecell experiments. This insight allowed for a model reduction (compare Fig 1 b and d). The implicit assumption made here was that MlrA and YdaM do not compete for binding to YciR, making MlrA a downstream component of the core regulatory network. Structurally, this is justified by the possibility of independent and simultaneous binding of YdaM and MlrA to YciR due to its large interaction domain (EALdomain) [43] and different binding modes of YdaM and MlrA on that domain. In addition, genetic data from Lindenberg et al. [7] yielded functional arguments against competition between YdaM and MlrA. Within our pipeline, the main result from the logical analysis passed on to the higher resolution modelling is the validation of the underlying network structure. We verified that the reduced topology shown in Fig 1 d can carry a dynamical model in agreement with all available data.
Beyond purely topological insights, a closer look at the regulatory mechanisms encoded by the logical functions can also be exploited to derive parameter constraints for the continuous reactionrate model. Analysing the impact of different regulators on a target in different system states may yield information on relations between production and decay rates. Theoretical results in this direction have been shown for specific classes of ODEs, see e.g. [38], but they are not yet widely applicable. We suggested how to exploit this idea in a wellsupported case for the curli regulation network and identified its potential for the sampling of kinetic parameters. Still, for a systematic application of this strategy the theoretical groundwork needs to be further extended.
An even closer intertwining of the formalisms would also have been possible. Automatic conversion procedures could have been used to lift a more abstract model into a more resolved formalism, e.g., to derive a generic ODE system from a given logical model as suggested in [44]. However, this approach yields models that cannot incorporate additional mechanistic knowledge (e.g. reaction parameters) in contrast to our approach, where we aimed at constructing a biologically realistic model by utilising all available kinetic information. As a consequence, the continuous model in our approach does not necessarily mimic the dynamics of the coarser models in all system states. It rather provides a complementary view based on the much higher detail resolution that allows to reevaluate and broaden the results of the logical analysis. Thus, a crucial advantage of our approach is the ability to incorporate experimental data from heterogeneous sources, carrying qualitative, semiquantitative or fully quantitative information. A further possibility would have been to take smaller steps on the modelling scale and add hybrid formalisms to the pipeline, e.g. discretecontinuous systems modelled as hybrid automata or Petri nets based on the boolean models that we derived in the first pipeline step [41, 45]. However, we found that in our context this does not provide any clear advantages, since the computational cost of parameter and dynamic analysis that go beyond qualitative aspects is already close to that of an ODE model while still incorporating many abstractions of the underlying mechanisms [46, 47].
The ODE model, derived in the second step of the pipeline, allowed us to identify regions in the multidimensional kinetic parameter space inducing bistable behaviour. Using a Monte Carlo sampling procedure, we identified a sufficiently large amount of parameter sets inducing bistability in the wild type strain and monostable behaviour in particular knockout strains and strains with point mutations (see “Methods”). Note that sampling procedures are frequently used for identifying parameter regions inducing bistability [23, 48]. We posthoc compared the identified parameter distributions with experimentally measured parameters in related systems. Thus, we could show that the network topology validated in the logical modelling step is capable of generating bistable dynamics within a biologically meaningful range of kinetic parameters, as shown in Fig. 3.
By conducting a bifurcation of key model parameters we exemplarily analysed the ranges where these parameters induce bistable regimes. This is of particular relevance since it was previously shown that during entry into the stationary phase of the growth cycle the σ ^{S}controlled DGC enzyme YegE is induced, while the expression of the PDE enzyme YhjH (controlled by the flagellar master regulators FlhDC and FliA) is turned down. The remaining levels of the YhjH protein are diluted and degraded within the following cell divisions [8]. Our results suggest that by finetuning the levels of YegE and YhjH the cells maintain parametric configurations, which generate bistable behaviour. Stochastic simulations indicate that a further increase of the protein levels of YegE during the stationary phase might reduce the probability of the curlioff state and increase the probability of the curlion state (see Fig. 4 a, inlay plots). This is in line with experimental measurements of the distribution of curli cells in bacterial colonies. Serra et al. [11] showed a spatial variation of this distribution, where the relative amount of curlion cells increases with an increasing vertical distance to the nutritional source at the bottom of the colony. In addition, Grantcharova et al. [10] observed a variation of this distribution within different stages of the stationary phase of the bacterial growth cycle. Furthermore, the bifurcation results suggested a narrow parameter region where the system is bistable. Changing these parameters in a way that eliminates the catalytic activity of the cdiGMPregulating enzymes YegE and YhjH, or decreases the binding affinity of cdiGMP to YciR yielded monostable systems which could explain experimental results of corresponding singlegene knockout strains lacking yegE and yhjH or carrying a y c i R ^{AAL} point mutation (see Fig. 4).
Finally, we addressed a possible mechanism inducing spontaneous switching between different stable states. Stochastic simulations of the system indicated that depending on kinetic parameters, qualitatively different types of noiseinduced switching dynamics may be observed, where most of the identified parameter sets indicated a higher probability of switching from the curlioff state to the curlion state than vice versa (72 %). In line with these results, small chains of curli cells were observed within macrocolonies of E. coli, suggesting that once the cells switch to curli expression, they maintain this state and inherit it to successive cell generations [11]. We identified a significant subset of parameters (12 %) exhibiting a low switching probability in either direction (lower left corner in Fig. 5 a). This suggests that certain parameter configurations might enable a longterm differentiation of the cells into curli producers and nonproducers, retaining the statusquo until system parameters change again e.g. due to a variation of environmental conditions. A possible regulatory mechanism for controlling the switching behaviour might be given by the level of separation of steady states within distinct phenotypes. In support of this, we observed that the deviation of stable cdiGMP levels between curlioff and curlion phenotypes was much more pronounced in parameter sets that favour a differentiation (rare switchers) than those parameter sets that allow frequent switching between the two phenotypes. A greater separation of these levels may induce a larger (energetic) barrier and robustness against stochastic fluctuations, which could be a potential mechanism controlling the phenotypic plasticity of E. coli. A comparison of the steady state levels in Fig. 7 with experimentally measured numbers of system components, in particular cdiGMP, in curlion and curlioff cells may shed further light on this topic.
Note that in this study we only focused on the noise resulting from posttranslational interactions. The implicit timescale separation assumption can be justified by the longer time scales of fluctuations in protein copy numbers as compared to the longest simulations in this study (≤100 s). In our simulations, we observed that switching from curlioff to curlion states was generally fast (≤100 s, see Fig. 5 b) and may thus not be affected by slower noise processes, that are related to protein copy number variations. However, switching in the opposite direction (curlion to curlioff) may occur at a slower timescale (≥100s, see Fig. 5 b) in some parameter sets and may thus be influenced by the noise resulting from transcriptional and translational events. Further kinetic data may enable a refinement of this model by explicitly including gene expression noise.
Conclusion
In this study we have presented a framework for the analysis of multistable signalling networks with an application to the stationary phaseinduced curli regulation system of E. coli. Using a hybrid sequential logicalcontinuous approach we first verified the potential of the curli regulation system for inducing bistability by incorporating genetic knockout experiments as discrete model constraints. Based on the validated network topology, we derived a reactionrate model of the system and identified several parameter sets that are capable of inducing bistable dynamics. Notably, most identified parameter sets are located in a biologically meaningful region. In addition, we analysed the dependence of the probability of the stable curlion and curlioff states on system parameters. We showed that certain parameter variations, e.g. the amount and activities of the DGCs and PDEs, give rise to different switching dynamics between curlion and curlioff states (i.e. frequent switching or stable differentiation).
The present study introduces the first mathematical model of bistable regulation of curli fibers in E. coli based on posttranslational interactions of multiple DGC and PDE proteins with the signalling molecule cdiGMP as the central player. Our results deliver new potential targets for further experimental investigations, which may help to test the modelling hypotheses and to successively deepen the understanding of biofilm control in bacteria.
Declarations
Acknowledgements
KPY and MvK acknowledge funding through the BMBFfunded JRG “Meth4SysPharm”, grant number 031A307.
Authors’ Affiliations
References
 Guo MS, Gross CA. Stressinduced remodeling of the bacterial proteome. Curr. Biol. 2014; 24(10):424–34.View ArticleGoogle Scholar
 Martínez JL, Rojo F. Metabolic regulation of antibiotic resistance. FEMS Microbiol Rev. 2011; 35(5):768–89.PubMedView ArticleGoogle Scholar
 Cornforth DM, Foster KR. Competition sensing: the social side of bacterial stress responses. Nat Rev Microbiol. 2013; 11(4):285–93.PubMedView ArticleGoogle Scholar
 Arnoldini M, Vizcarra IA, PeñaMiller R, Stocker N, Diard M, Vogel V, Beardmore RE, Hardt WD, Ackermann M. Bistable expression of virulence genes in Salmonella leads to the formation of an antibiotictolerant subpopulation. PLoS Biol. 2014; 12(8):1001928.View ArticleGoogle Scholar
 Deris JB, Kim M, Zhang Z, Okano H, Hermsen R, Groisman A, Hwa T. The innate growth bistability and fitness landscapes of antibioticresistant bacteria. Science. 2013; 342(6162):1237435.PubMed CentralPubMedView ArticleGoogle Scholar
 Smits WKK, Kuipers OP, Veening JWW. Phenotypic variation in bacteria: the role of feedback regulation. Nat Rev Microbiol. 2006; 4(4):259–71.PubMedView ArticleGoogle Scholar
 Lindenberg S, Klauck G, Pesavento C, Klauck E, Hengge R. The EAL domain protein YciR acts as a trigger enzyme in a cdiGMP signalling cascade in E. coli biofilm control. The EMBO journal. 2013; 32(14):2001–014.PubMed CentralPubMedView ArticleGoogle Scholar
 Pesavento P, Becker G, Sommerfeldt N, Possling A, Tschowri N, Mehlis A, Hengge R. Inverse regulatory coordination of motility and curlimediated adhesion in Escherichia coli. Genes & Development. 2008; 22:2434–446.View ArticleGoogle Scholar
 Weber H, Pesavento C, Possling A, Tischendorf G, Hengge R. CyclicdiGMPmediated signalling within the σ ^{ S } network of Escherichia coli. Mol Microbiol. 2006; 62(4):1014–1034.PubMedView ArticleGoogle Scholar
 Grantcharova N, Peters V, Monteiro C, Zakikhany K, Römling U. Bistable expression of CsgD in biofilm development of Salmonella enterica Serovar Typhimurium. J Bacteriol. 2010; 192(2):456–66.PubMed CentralPubMedView ArticleGoogle Scholar
 Serra DO, Richter AM, Klauck G, Mika F, Hengge R. Microanatomy at cellular resolution and spatial order of physiological differentiation in a bacterial biofilm. mBio. 2013; 4(2):e10010313.View ArticleGoogle Scholar
 Klarner H, Siebert H, Bockmayr A. Time series dependent analysis of unparametrized Thomas networks. IEEE/ACM Trans Comput Biol Bioinformatics. 2012; 9(5):1338–1351.View ArticleGoogle Scholar
 SaezRodriguez J, Simeoni L, Lindquist JA, Hemenway R, Bommhardt U, Arndt B, Haus UUU, Weismantel R, Gilles ED, Klamt S, Schraven B. A logical model provides insights into T cell receptor signaling. PLoS Comp Biol. 2007; 3(8):163.View ArticleGoogle Scholar
 Morris MK, SaezRodriguez J, Sorger PK, Lauffenburger DA. Logicbased models for the analysis of cell signaling networks. Biochemistry. 2010; 49(15):3216–224. http://pubs.acs.org/doi/pdf/10.1021/bi902202q.PubMed CentralPubMedView ArticleGoogle Scholar
 Calzone L, Tournier L, Fourquet S, Thieffry D, Zhivotovsky B, Barillot E, Zinovyev A. Mathematical modelling of cellfate decision in response to death receptor engagement. PLoS Comput Biol; 6(3):1000702.Google Scholar
 Mbodj A, Junion G, Brun C, Furlong EEM, Thieffry D. Logical modelling of drosophila signalling pathways. Mol BioSyst. 2013; 9:2248–258.PubMedView ArticleGoogle Scholar
 Stigler B, VelizCuba A. Boolean models can explain bistability in the lac operon. J. Comput. Biol. 2011; 18(6):783–94.PubMedView ArticleGoogle Scholar
 Thattai M. Using topology to tame the complex biochemistry of genetic networks. Philos Trans R Soc A: Math Phys Eng Sci. 2013; 371(1984):20110548.View ArticleGoogle Scholar
 Tiwari A, Ray, Narula J, Igoshin OA. Bistable responses in bacterial genetic networks: Designs and dynamical consequences. Math Biosci. 2011; 231(1):76–89.PubMed CentralPubMedView ArticleGoogle Scholar
 Shinar G, Feinberg M. Concordant chemical reaction networks. Math Biosci. 2012; 240(2):92–113.PubMedView ArticleGoogle Scholar
 Amin M, Porter SL, Soyer OS. Split histidine kinases enable ultrasensitivity and bistability in twocomponent signaling networks. PLoS Comput Biol. 2013; 9(3):1002949.View ArticleGoogle Scholar
 Eissing T, Waldherr S, Allgöwer F, Scheurich P, Bullinger E. Response to bistability in apoptosis: Roles of Bax, Bcl2, and mitochondrial permeability transition pores. Biophys J. 2007; 92(9):3332–334.PubMed CentralPubMedView ArticleGoogle Scholar
 Chickarmane V, Paladugu SR, Bergmann F, Sauro HM. Bifurcation discovery tool. Bioinformatics. 2005; 21(18):3688–690.PubMedView ArticleGoogle Scholar
 Ferm L, Lötstedt P. Adaptive solution of the master equation in low dimensions. Appl Numerical Math. 2009; 59(1):187–204.View ArticleGoogle Scholar
 Munsky B, Khammash M. The finite state projection algorithm for the solution of the Chemical Master Equation. J Chem Phys. 2006; 124(4):044104.PubMedView ArticleGoogle Scholar
 Menz S, Latorre JC, Schütte C, Husinga W. Hybrid stochasticdeterministic solution of the Chemical Master Equation. Multiscale Model Simul. 2012; 10(4):1232–1262.View ArticleGoogle Scholar
 Dandach SH, Khammash M. Analysis of stochastic strategies in bacterial competence: A Master Equation approach. PLoS Comput Biol. 2010; 6(11):1000985.View ArticleGoogle Scholar
 Gérard C, Gonze D, Lemaigre F, Novák B. A model for the epigenetic switch linking inflammation to cell transformation: Deterministic and stochastic approaches. PLoS Comput Biol. 2014; 10(1):1003455.View ArticleGoogle Scholar
 Rath BA, Yousef KP, Katzenstein DK, Shafer RW, Schütte C, von Kleist M, Merigan TC. In vitro HIV1 evolution in response to triple reverse transcriptase inhibitors & in silico phenotypic analysis. PLoS One. 2013; 8(4):61102.View ArticleGoogle Scholar
 AbouJaoudé W, Ouattara DA, Kaufman M. From structure to dynamics: frequency tuning in the p53Mdm2 network I. Logical approach. J Theoret Biol. 2009; 258(4):561–77.View ArticleGoogle Scholar
 Ouattara DA, AbouJaoudé W, Kaufman M. From structure to dynamics: frequency tuning in the p53Mdm2 network. II Differential and stochastic approaches. J Theoret Biol. 2010; 264(4):1177–1189.View ArticleGoogle Scholar
 Albert R, Othmer HG. The topology of the regulatory interactions predicts the expression pattern of the segment polarity genes in drosophila melanogaster. J Theoret Biol. 2003; 223(1):1–18.View ArticleGoogle Scholar
 VelizCuba A, Arthur J, Hochstetler L, Klomps V, Korpi E. On the relationship of steady states of continuous and discrete models arising from biology. Bull Math Biol. 2012; 74(12):2779–792.PubMedView ArticleGoogle Scholar
 Wittmann DM, Krumsiek J, SaezRodriguez J, Lauffenburger S. D. A, Klamt, Theis FJ. Transforming boolean models to continuous models: methodology and application to Tcell receptor signaling. BMC Syst Biol. 2009; 3(1):1–21.View ArticleGoogle Scholar
 Hengge R. Principles of cdiGMP signalling in bacteria. Nat Rev Microbiol. 2009; 7(4):263–73.PubMedView ArticleGoogle Scholar
 Alon U. An Introduction to Systems Biology: Design Principles of Biological Circuits (Chapman & Hall/CRC Mathematical and Computational Biology). Boca Raton, FL: Chapman and Hall/CRC; 2006.Google Scholar
 Christen B, Christen M, Paul R, Schmid F, Folcher M, Jenoe P, Meuwly M, Jenal U. Allosteric control of cyclic diGMP signaling. J Biol Chem. 2006; 281(42):32015–2024.PubMedView ArticleGoogle Scholar
 Thomas R, d’Ari R. Biological Feedback. Boca Raton, FL: CRC Press; 1990.Google Scholar
 Kubitschek HE, Friske JA. Determination of bacterial cell volume with the Coulter Counter. J Bacteriol. 1986; 168(3):1466–1467.PubMed CentralPubMedGoogle Scholar
 Zhang Q, Bhattacharya S, Andersen ME. Ultrasensitive response motifs: basic amplifiers in molecular signalling networks. Open Biology. 2013; 3(4):130031.PubMed CentralPubMedView ArticleGoogle Scholar
 Chaouiya C. Petri net modelling of biological networks. Brief Bioinform. 2007; 8(4):210–9.PubMedView ArticleGoogle Scholar
 Gillespie DT. Exact stochastic simulation of coupled chemical reactions. J Phys Chem. 1977; 81:2340–381.View ArticleGoogle Scholar
 Schirmer T, Jenal U. Structural and mechanistic determinants of cdiGMP signalling,. Nat Rev Microbiol. 2009; 7(10):724–35.PubMedView ArticleGoogle Scholar
 Krumsiek J, Pösterl S, Wittmann DM, Theis FJ. Odefy – from discrete to continuous models. BMC Bioinformatics. 2010; 11:233.PubMed CentralPubMedView ArticleGoogle Scholar
 Siebert H, Bockmayr A. Temporal constraints in the logical analysis of regulatory networks. Theoret Comput Sci. 2008; 391(3):258–75.View ArticleGoogle Scholar
 Alur R, Henzinger TA, Lafferriere G, Pappas GJ. Discrete abstractions of hybrid systems. In: Proc. IEEE. IEEE: 2000. p. 971–84. www.ieee.org.
 Lygeros J, Johansson KH, Simic SN, Zhang J, Sastry SS. Dynamical properties of hybrid automata. IEEE Trans Automat Control. 2003; 48:2–17.View ArticleGoogle Scholar
 Qiao L, Nachbar RB, Kevrekidis IG, Shvartsman SY. Bistability and oscillations in the HuangFerrell model of MAPK signaling. PLoS Comput Biol. 2007; 3(9):184.View ArticleGoogle Scholar
 Christen M, Christen B, Folcher M, Schauerte A, Jenal U. Identification and characterization of a cyclic diGMPspecific phosphodiesterase and its allosteric control by GTP. J Biol Chem. 2005; 280(35):30829–0837.PubMedView ArticleGoogle Scholar
 Lindenberg S. CdiGMP Signaltransduktion in der Regulation der Expression des Biofilmregulators CsgD in Escherichia coli. PhD thesis. 2013. http://www.diss.fuberlin.de/diss/receive/FUDISS\_thesis\_000000094976.
 Spangler C, Kaever V, Seifert R. Interaction of the diguanylate cyclase YdeH of Escherichia coli with 2’,(3’)substituted purine and pyrimidine nucleotides. J Pharmacol Exp Ther. 2011; 336(1):234–41.PubMedView ArticleGoogle Scholar
 Milo R, Jorgensen P, Moran U, Weber G, Springer M. BioNumbers  the database of key numbers in molecular and cell biology. Nucleic Acids Res. 2010; 38(suppl 1):750–3.View ArticleGoogle Scholar
Copyright
This is an Open Access article distributed under the terms of the Creative Commons Attribution License (http://creativecommons.org/licenses/by/4.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly credited. 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.
Comments
View archived comments (1)