Full text

Turn on search term navigation

1. Introduction

The lac operon in Escherichia coli is a paradigmatic example of a genetic regulation system involving the interaction of positive and negative regulatory molecules. The system encodes a pathway for lactose catabolism that is hierarchically controlled by glucose availability, and it has been reported to exhibit a bistable behavior, since the catabolic genes are either uninduced (OFF) or induced (ON) in a single cell, depending on previous activation history and specific extracellular lactose and/or glucose concentrations [1,2]. The lac operon has been used as a model system for gene regulation since its initial description in [3], where the concept of an operon was first introduced. Since the molecular components of the regulatory system have been well characterized, it is an excellent candidate for global analysis and modeling. The operon includes three genes involved in the uptake and catabolism of lactose and/or structurally similar sugars. The first gene, lacZ, encodes for the enzymeβ-galactosidase, which converts lactose to allolactose and to subsequent catabolic intermediates; lacY, the second structural gene of the pathway, encodes for lactose permease, a membrane transporter for lactose uptake. Finally, the product of gene lacA is an acetyltransferase that takes part in the degradation and excretion of non-lactose sugars that may be misrouted through the lactose degradation pathway.

The fundamental states of the lac operon are described schematically in Figure 1. In brief, when a low level of extracellular lactose occurs, a strong binding of the LacI protein to the operator sequences strongly represses lac gene expression, thus allowing only a very low amount of mRNA to be transcribed (Figure 1a), which can be considered as an OFF state. However, even at this low expression level, a few LacY permease molecules can reach the cytoplasmic membrane and take up any eventual lactose that appears outside the cells. The presence of low amounts of LacZ in the cytoplasm leads to transformation of this lactose into allolactose, the true natural inducer of the lac operon. When this happens, allolactose interacts allosterically with the LacI repressor protein, decreasing its affinity for operator sequences and releasing the negative control of the lacZ promoter, rendering the system into an ON state (Figure 1b). Positive feedback loops, like the one described here for the lac system, have been proposed as a typical property of bistable systems [4].

On the other hand, an essential feature of the regulatory system is that high extracellular glucose concentrations can also inhibit the expression of lac genes, a phenomenon traditionally considered as a type of catabolite repression that is directly mediated by the glucose molecule. Such a repression has usually been assumed to involve the transcriptional activator protein Crp (cAMP receptor protein), which also takes part in lac gene expression as a cAMP-dependent positive regulator (Figure 1). In the repression model, glucose uptake by the cells produces a drastic reduction in the intracellular cAMP levels, negatively affecting Crp activity and abolishing lac gene expression, which would explain the metabolic preference of glucose over lactose in E. coli [5]. An additional explanation of the catabolic hierarchy of glucose is given by the inducer exclusion effect, which involves glucose-mediated inhibition of lactose uptake, by means of a direct reduction of LacY permease activity [6]. This would greatly contribute to the regulatory effect of glucose by preventing the inducer molecule to reach the LacI repressor.

Bistability is normally defined as a condition where a system can respond to the same external signal or input in two different manners, depending on the internal state and/or the previous history of the system. However, true bistability requires the system to switch in an all-or-none way between alternative fixed points, without the appearance of intermediate states or a “quantitative” cellular response. Although a quantitative continuous response to the presence of lactose, or the gratuitous inducer thiomethyl-β -D-galactoside (TMG), was initially shown in experiments measuring the average response of whole culture populations, containing several millions of cells, observations with individual bacteria showed unequivocally that each cell responds in a discrete manner, switching alternatively among the OFF and ON states [1,2], with virtually no evidence of intermediate responses. This shows that bistability can be a characteristic of the lac system. Using genetically engineered E. coli reporter cells, gratuitous inducers and sophisticated experiments, the authors in [7] showed a discrete hysteretic response of lac gene expression in individual bacteria when pre-incubated in the presence or absence of extracellular lactose and then transferred to fresh medium with different lactose concentrations. Their results showed that, depending on the cell’s previous history, the response profile to lactose in the fresh medium varied significantly within a specific concentration range, while outside its limits, cells would respond independently of their previous incubation conditions: Any concentration lower than 3μM would not lead to lacZ induction, whereas any concentration higher than 30μM would induce the reporter gene. However, most values within the 3–30μ M window were found to induce lacZ expression only for lactose preincubated cells. This window of bistability was very clearly defined in the absence of extracellular glucose but when this sugar was present, the concentration window moved towards higher lactose concentrations (see Figure 2c in [7]).

Thus, the above all-or-none and bistable behavior that characterizes the lac system makes it suitable for discrete modeling. In particular, Boolean models have been used successfully to analyze genetic regulatory systems on numerous occasions [8,9,10]. They are useful for the analysis of large regulatory networks when mechanistic details underlying are scarce, insufficiently defined, or even controversial, as it has been discussed that, in such models, the interaction type (inhibition or activation) and network topology are enough for capturing dynamics of gene networks [11,12], without the need of estimating parameters, thereby reducing model complexity. For example, they are being used extensively to represent interactions revealed by high throughput sequencing technologies in the context of systems biology [13,14]. Although some of these assumptions make them less applicable for certain biological processes, they have proved very valuable for predicting the behavior of genes during bacterial metabolism [15] or in the context of interspecies interactions [16]. Even more, other researchers have been able to reproduce bistability for eukaryotic cell differentiation [17], surfactant production in bacteria [18] and microbial signaling [19] using Boolean models.

In this context, the authors of [20] modeled the lac regulation system using a Boolean network that was able to reproduce the bistability of the system under certain but not all of the conditions described in [7]. Figure 2 displays a scheme of this model, whereGeand(Le,Lem) are parameters that account for different concentration levels of glucose and lactose, respectively. This allows to assign specific values as inputs in order to predict the final outcome of the network in terms of fixed points (steady states) that are consistent with the operon being ON, OFF, or in a bistable condition. When the nodes are updated in a parallel regime, the model can reproduce the experimental results in [7] when extracellular glucose is low or absent. However, the predicted outcome of the model differs from the actual behavior of the lac system when lactose and glucose concentrations are high, a feature that limits the usefulness of this initial model for representing lactose fermentation by E. coli.

In this work, we first analyze the dynamics of the Boolean models considered but under all different updating schemes in order to know how robust they are in the sense of whether or not new attractors appear. In general, this is a difficult task because the number of different updates schemes associated to a network grow exponentially with its size [21]. Fortunately, there are currently mathematical tools that can significantly reduce this analysis [22]. Secondly, we propose improvements for the performance of the model considered by a rational redesign of the specific Boolean functions in order to predict the operon behavior in each of the sugar combinations tested experimentally in [7], including combined scenarios of low, medium and high concentrations of glucose and lactose. Furthermore, we have taken into account current biological data, regarding catabolic repression and inducer exclusion mechanisms, to bring function definitions up to date.

The paper is organized as follows; Section 2 summarizes the basic mathematical concepts and definitions used along the manuscript, Section 3 contains the most relevant aspects taking into account of [20], in Section 4 our main mathematical results are presented, in Section 5 alternative improvements for the studied models are proposed. Finally, we discuss on our findings in the conclusions section.

2. Mathematical Background

A Boolean Network (BN) is a couple(G,F), whereG=(V,E) is a finite directed graph named interaction digraph (e.g., Figure 2(left)), V is a set of nnodes (n also called the network size),E⊆V×Vis the set of edges andF=(f1,...,fn):{0,1}n→{0,1}nis a Boolean function composed by n local functionsfi:{0,1}n→{0,1}such thatx∈{0,1}n→xi=fi(x)∈{0,1}(x is called a configuration while the value of the variablexkassociated with node k is known as an state). Besides, the local functionsfidepends only on variablesxjsuch that(j,i)∈E (e.g., Figure 2(right)).

Among the different ways that the states of a BN can be updated is the family of (deterministic) update schedules [22] defined as functionss:V→{1,...,n}such thats(V)={1,...,m}for somem≤n. In particular, the well-known parallel (or synchronous) and sequential (or asynchronous) update schedules are obtained whenm=1(i.e., all states are updated at the same time) andm=n(i.e., all states are updated at different times but following a predetermined order), respectively. We denote bySnthe set of deterministic update schedules that exist for a BN of size n. If the states of a configurationx∈{0,1}nare updated according to a given update schedule s, a new configurationx′=F(x)∈{0,1}nwill be obtained. The change from x tox′can be represented by the arcx→x′and is known as transition from x tox′. Thus, the2nconfigurations and their respective transitions give rise to the dynamics of the system which can be represented in a state transition graph Furthermore, because it is finite, the dynamics has limit behaviors called attractors, which are of 2 types:

  • Fixed point (or steady state): configurationx∈{0,1}nsuch thatF(x)=x.

  • Limit cycle: sequence of configurationsx0,...,xl−1∈{0,1}npairwise distinct such thatF(xj)=xj+1, for allj=0,...,l−2,F(xl−1)=x0andl>1is a integer named the length of the limit cycle.

The set of configurations that converge to a specific attractor is called the attraction basin. 3. The lac Operon Boolean Models to Be Considered: Main Aspects and Choice Justification

In [20] the authors proposed different Boolean models explaining the bistability of the lac operon in Escherichia coli. These are BNs that have the following characteristics in common:

1.

All of them have qualitative behaviors that match very well with the experiments performed in [7].

2. The interaction digraph is composed by nodes that represent mRNA, proteins, and sugars. Its edges represent the type of interaction between the nodes (activation/inhibition).

3. The dynamic is obtained considering the parallel update schedule.

Next, we will specify which are the models of [20] that we choose. Hence, we will summarize the main points that their authors studied and with which we will work from Section 4 onwards. Finally, we justify our choices.

3.1. The Chosen Models of Veliz-Cuba and Stigler 2011: Those without Catabolite Repression

We will consider the alternative model without catabolite repression proposed in [20] which differs from the initial one presented in that work only in the definition of the local functionfC; in the initial onefC=¬Gewhile in the alternative model without catabolite repressionfC=1 (see our justification in Section 3.4). In addition, we will also consider its reduced version (without catabolite repression) and, from now on, we will refer to this alternative model without catabolite repression and its reduced version as the original and reduced one, respectively. We show the details of these models in Figure 2 and Figure 3, where we point out that the networks have been drawn “disarming” those of [20] which consider some pairs of nodes as a single node, but the reader can easily check that our networks and the corresponding of [20] are equivalents.

3.2. The Dynamics Produced by the Original and Reduced Models

According to [20], the original model (Figure 2) assumes two concentration of glucose, low (Ge=0) or high (Ge=1), and three concentrations of extracellular lactose, low, medium or high, represented by the pairs of parameters(Le,Lem)=(0,0),(0,1)and(1,1), respectively (the pair(Le,Lem)=(1,0)has no meaning, therefore is not considered). This gives a total of six different combinations of parameters and, consequently, each one will have its corresponding dynamic. Here, a configuration corresponds to the vector(M,P,B,C,R,Rm,A,Am,L,Lm)∈{0,1}10 and, for purposes of relating the dynamics to the results of [7] (see details in Section 3.3), only its steady states matter, i.e., limit cycles (if they exist) are ignored. Besides, when(M,P,B)=(1,1,1)(when(M,P,B)=(0,0,0)) it will be said that the operon is ON (OFF). Abusing this notation, henceforth, we will call ON (OFF) to a steady state that verifies the above condition.

Therefore, if the dynamics for each of the six combinations of parameters are obtained and only are observed its steady states, then, there exist the following four cases:

Case 1:

Ge=1. Any configuration eventually converges to the unique fixed point OFF =(0,0,0,1,1,1,0,0,0,0).

Case 2:

Ge=Le=Lem=0. Any configuration eventually converges to the unique fixed point OFF=(0,0,0,1,1,1,0,0,0,0).

Case 3:

Ge=0∧Le=Lem=1. Any configuration eventually converges to the unique fixed point ON=(1,1,1,1,0,0,1,1,1,1).

Case 4:

Ge=Le=0∧Lem=1. Any configuration eventually converges to one of the two fixed points; OFF=(0,0,0,1,1,1,0,0,0,0)or ON=(1,1,1,1,0,0,0,1,0,1), i.e., bistability is obtained.

Notice that case 1 implicitly covers three of the six combinations of parameters (those obtained when high concentration of glucose is combined with low, medium and high concentration of extracellular lactose, respectively), all three having the same dynamic. The remaining three combinations of parameters are those obtained in cases 2, 3 and 4, each with a different dynamic.

We summarize the four cases in Table 1 and their respective dynamics in Figure 4.

Regarding the reduced model (Figure 3), four cases similar to the original one are obtained, which can be summarized using the same Table 1. It should be mentioned that, in this reduced model, a configuration is of the form(M,L,Lm)∈{0,1}3, and the operon is ON (OFF) whenM=1(M=0 ). In the Figure 5, we show the dynamics of its four cases.

Notice that although only 3 different dynamics are generated (because its cases 1 and 2, involving 4 out of 6 combinations of parameters, generate a single dynamic), they are still in accordance with the summary of Table 1.

3.3. Stochastic Simulations in the Original Model Which Largely Coincide with the Biological Experiments of Ozbudak et al., 2004

The bistability of the lac Operon, according to [23], means that it must be a region of bistability, i.e., for a range of parameters, a set of “cells” may have both, OFF and ON, fixed points. In [20], this is accomplished by considering stochasticity in the uptake of the inducer in the following way, considering that glucose is absent (Ge=0) while that extracellular lactose is given as the discretization of a random variableL∼N(μ,σ)taken from a normal distribution withμequal to the mean of{0.6,...,2.8}incremented by 0.2 andσ=0.3 (these values are chosen in order to be consistent with those measured in the biological experiments of [7]). Then, the extracellular lactose(Le,Lem)is given as a function ofLdefined by

(Le,Lem)=(0,0)=low,ifL≤1(0,1)=medium,if1≤L<2(1,1)=high,if2≤L

Below, we summarize the stochastic experiment carried out in [20]:

(1)

Starting with the normal distributionN(0.6,0.3)forL, to generate randomly a set of values forLand calculate for each of them the corresponding value for the pair(Le,Lem).

(2)

AssumingGe=0, to find which are the steady states of the dynamic obtained for each value(Le,Lem) of (1), it is simply one of the three possibilities showed in the second column of Table 1; OFF, bistable (i.e., OFF and ON) or ON.

(3)

To repeat (1) and (2) but forN(μ,0.3)withμ∈{0.8,1.0,1.2,...,2.8}.

The results of this experiment were those that matched with the biological experiments of [7] in 5 out of 6 combinations of parameters described in Section 3.2. In the remaining one, where glucose and lactose concentrations are high (see Table 1), the original model (as well as the reduced one) predicts the operon being OFF while the biological experiments of [7] show a bistable behavior.

3.4. Our Justification for the Choice of Models without Catabolite Repression

One of the most relevant outcomes of the Boolean models proposed in [20] is that the presence of glucose at a medium or high concentration would be able to completely repress the expression of the lac mRNA, which is a direct consequence of the inhibition functionfC=¬Ge , which aims to take into account one of the proposed aspects of catabolite repression. However, this assumption fails to consider the well-described fact that glucose is not able to shut down LacZ expression in bacteria when the operon has been pre-induced experimentally [24]. Furthermore, the whole concept of glucose repression mediated by a drop in cAMP and subsequent lack of the CRP-cAMP activator complex has been under critical revision [25,26,27,28] based on the following compelling challenges: experimental data indicate that intracellular cAMP concentrations are similar when lactose, glucose, or a mixture of both substrates is present at concentrations lower than 300μM and that cAMP abundance is only reduced by a factor of 5–8 when external glucose concentrations are above 300μ M [29,30], which is higher than the bistability window explored in [7]. Exogenous addition of cAMP did not abolish glucose-mediated inhibition of lac gene expression, even though the metabolite was shown to enter the cells [31]. On the other hand, when the repressor protein LacI is inactivated by deletion of the lacI gene, no catabolite repression of LacZ is observed, indicating that this process relies solely on the activity of the repressor protein, shutting down transcriptional activity by means of inducer exclusion [32]. Finally, the CRP-cAMP activator complex is also required for the expression of the PtsG glucose uptake transporter [33], which would imply that, if this pathway of catabolite repression were valid, glucose would eventually inhibit its own transport, which has never been observed. This inconsistency with experimental data and the relevant objections discussed above justify our choice of the models without catabolite repression of [20], wherefC=1, to account for the current view of glucose-induced lac repression.

4. Results: Dynamical Robustness of the Original and Reduced Models

Our goal is to study the dynamics of the original and reduced models not only for the parallel update schedule but for the whole family of deterministic update schedules defined in Section 2. As evidenced in [21], this set grows exponentially with the network size. In particular, its sizes are|S3|=13and|S10|=102,247,563for the reduced and original models, respectively.

Observe that it is easy to deduce that the steady states OFF and ON obtained in the dynamics of the original, and reduced models will also continue to appear with any other deterministic update schedule because their are known to be invariant under update schemes; however, limit cycles could also appear making the desired modeling less robust. We will analyze this aspect below. 4.1. Dynamical Robustness of the Original Model

We begin by giving a short known Lemma that allows us to prove that the dynamics of the first three cases mentioned in Section 3.2 are highly robust for any deterministic update schedule.

Lemma 1.

LetN=(G,F)be a Boolean network updated by the update schedule s. If G is acyclic and it has no loop [i.e., an edge of the form (i,i)], then the attractors are only fixed points.

Proof.

It is known that the nodes of an acyclic digraph can be ordered in m layers, being the first and last those that have the root and leaves nodes, respectively, such that the arcs belongs between the i-layer and the j-layer, for1≤i≤m−1andi<j. Since there are no loops, the dynamic evolves fixing the states of the nodes, in the worst case, from the first to the last layer. Therefore, the only possible attractor is a fixed point. □

4.1.1. Cases 1, 2 and 3 for the Original Model

Proposition 1.

The dynamics associated with the cases 1, 2 and 3 of the original model of Section 3.2 have no limit cycle, whatever the update schedule considered.

Proof.

Let s be an arbitrary update schedule and(M,P,B,C,R,Rm,A,Am,L,Lm)∈{0,1}10an initial configuration at timet=0 . We will show that in every case, we can obtain one of the acyclic digraphs without loops showed in Figure 6:

Cases 1 and 2.

It is easy to check thatC=1andL=Lm=0,∀t≥1. This implies that,A=Am=0,∀t≥2 . Therefore, we have the left digraph of Figure 6.

Case 3.

Notice thatC=Lm=1,∀t≥1. Therefore, we have the following sequence of implications;Am=1,∀t≥2⇒R=0,∀t≥3⇒Rm=0,∀t≥4 . Therefore, we have the right digraph of Figure 6.

Thus, in every case, Lemma 1 assures that the attractors are only fixed points, i.e., there is no limit cycle. □

Proposition 1 implies that for the first 3 cases, and whatever be the update schedule, the attraction basin of its dynamics will be the whole space of configurations{0,1}10 , as summarized in Table 2.

4.1.2. Case 4 for the Original Model

It is easy to see that for an arbitrary update schedule and a initial configuration(M,P,B,C,R,Rm,A,Am,L,Lm)∈{0,1}10att=0, we will have thatC=1andL=0,∀t≥1. This implies thatA=0,∀t≥2 , and we obtain the digraph of Figure 7 that has three cycles; one of length 6 and two of length 5. Therefore, we cannot apply the previous Lemma.

We will analyze its dynamical behavior through a detailed study of all its different dynamics by using the algorithms developed in [22] that, roughly speaking, list all the sets of schemes that generate exactly the same dynamics, significantly reducing the number of dynamics to analyze. These algorithms have been shown to be effective in obtaining new and relevant information for the study of concrete biological networks [34,35] as well as in the field of the cellular automata theory [36,37,38,39,40,41].

Simply, we generate one dynamic for each of the above representative update schedules and summarize its results in Table 3.

We can observe that less than 30% of all the dynamics present limit cycles and with a balanced proportion between ON/OFF basins, but they are smaller than those associated to the limit cycles, which are on average 5 times larger than that of ON and OFF, respectively. On the other hand, more than 70% of all the dynamics did not have any limit cycle, i.e., they only had the fixed points OFF and ON (being, on average, about 8 times greater the OFF basin than the ON basin).

In Figure 8 we exhibit one example of the dynamical behavior that may have the original model when the update schedule considered is slightly different from the parallel one. There are, in addition to the fixed points ON an OFF, three limit cycles of length 4 and one of length 2.

4.2. Dynamical Robustness of the Reduced Model Clearly, the dynamical robustness of the reduced model is quick and straightforward to analyze, so we only summarize its results.

4.2.1. Cases 1, 2 and 3 for the Reduced Model

Proposition 2.

The dynamics associated with the cases 1, 2 and 3 of the reduced model of Section 3.2 have no limit cycle, whatever the update schedule considered.

The dynamics for the above three cases also are quite robust for any deterministic update schedule, and this is summarized in Table 4.

4.2.2. Case 4 for the Reduced Model

Now the exhaustive analysis boils down to just 13 possible update schedules and the results are summarized in Table 5.

Observe that about 77% of the dynamics have only the fixed points OFF and ON, but, in this case the average size of the ON attractor basin is almost 2 times bigger than that of OFF. The same occurs in the other 23%. Furthermore, the OFF and ON basins add up, on average, exactly the same as that of the limit cycles.

Examples of the dynamical behavior of the reduced model exhibiting limit cycles are those of Figure 5(left).

5. Alternative Improvements for the Studied Models

As mentioned at the end of Section 3.3, of the 6 parameter combinations,(Ge,Le,Lem)=(1,1,1) –meaning glucose and lactose concentrations in high levels (see Table 1)–it is the only one that does not matches with the experimental data (see Figure 2c in [7]) where bistability should be observed again, instead of the operon being OFF as obtained from the models. Therefore, the first improvement naturally consists in solving that in both models. The second one will consist of increasing the parameters of the original model from 6 to 9, matching in all these cases with the biological experiments of [7].

5.1. Improvement 1: The Original and Reduced Models Match in All 6 Parameter Combinations with Ozbudak et al., 2004

Observe that according to what was discussed at the beginning of Section 4, some local functions of both models must necessarily be modified to make the steady states OFF and ON appear at the same time (bistability). On the other hand, by considering carefully the arguments we gave in Section 3.4, the local functionfC=1is not sufficient to adjust the Boolean model to represent the actual behavior of the system when both glucose and lactose are present in the extracellular environment of E. coli. To further refine the response of the Boolean network to the balance of each different sugar and to highlight the relevance of the previous history of the cell, a new definition has been proposed forfLandfLm , considering the inhibitory effect of glucose on lactose transport will be significant in the short term as long as a low concentration of lactose was initially present, while the opposite will be true if a high lactose concentration was already inside the cell. These changes are shown in Figure 9 and Figure 10.

Besides, we show in Figure 11 and Figure 12 the interaction digraph and the dynamics of the improved original and reduced models, respectively.

Thus, the improved original and reduced models predicts:

(1)

Bistability: when(Ge,Le,Lem)∈{(0,0,1),(1,1,1)}.

(2)

OFF: when(Ge,Le,Lem)∈{(0,0,0),(1,0,0),(1,0,1)}.

(3)

ON: when(Ge,Le,Lem)=(0,1,1).

The above is summarized in Table 6, which match for all the 6 combinations of parameters with the biological experiments of [7].

5.2. Improvement 2: (Improved) Original Model Extended to 9 Parameters

As explained in the introduction, there is a window of bistability (see Figure 2c in [7]) which justifies that three parameters can also be considered for the glucose (low, medium and high). Taking into account again the arguments given in Section 3.4 as well as those discussed for the first improvement, the following changes are proposed over the improved original model (see Figure 13 and Figure 14).

At this point, we repeat the stochastic experiment described in Section 3.3, but for different glucose values. Its results are shown in Figure 15 and summarized in Table 7.

Observe how the results of this extended model match the entire window of bistability (see Figure 2c in [7]) in each of its 9 combinations of parameters.

6. Conclusions

The lac Operon Boolean model without catabolite repression and its respective reduced version proposed in [20] were analyzed. These models have the particularity of being simple but capable of reproducing the operon being OFF, ON and bistable for different levels of lactose and glucose, matching very well with biological experiments such as those published in [7]. Unlike the models of [20], the Boolean networks proposed here predict bistability even at high glucose concentrations. This feature has been observed experimentally when inducer (a non-metabolizable lactose analogue) concentrations are also high [7]. Furthermore, our models take into account intermediate glucose concentrations, thus increasing sensibility to changes within the “window of bistability” of the lac system (see more details in the second part below).

In a first part of this paper, we have studied the dynamical robustness of these models where we made the following contributions:

  • For the first 3 cases described in Section 3.2 and that included 5 of the 6 combinations of parameters allowed in the original and reduced models, we establish two Propositions proving the non-existence of any limit cycle, whatever the update schedule used.

  • For the case 4, where bistability appears, we made an exhaustive analysis of all its possible dynamics generated with any update schedule. Here we detail for both models the average sizes of their attraction basins, the number of dynamics without limit cycles (i.e., only with fixed points) and the number of dynamics with limit cycles (being less than 30% in both models).
  • Again in the case 4, its predominant attractor (those that have the bigger attraction basin), changes dramatically; OFF attraction basin being, in average, 8 times bigger than that of ON for the original model while that in the reduced one, the ON basin is almost 2 times bigger than that of OFF.

The above findings allow us to conclude that the effect of the glucose and lactose parameters in the interaction digraph associated to each network breaks its cycles for cases 1, 2 and 3, transformed into an acyclic network that will have only fixed points as attractors and, consequently, being completely robust in these situations. However, in the case of bistability, the role of the update schedule can change this property significantly because limit cycles with large attraction basins can appear. Moreover, this work supports the hypothesis of [20] that network topology is a key factor for qualitative dynamical properties but not for quantitative ones.

As a second part, we have proposed two alternative improvements for the Boolean models studied, with biological support and that involve small modifications in their local functions. In the first improvement, the prediction is corrected for the case in which the parameters of the models represent high glucose and lactose levels, achieving the bistability observed in the biological experiments of [7]. In this way, with such an improvement, the original and reduced models match perfectly with the above experiments of [7] for all the 6 combinations of parameters. The second improvement increases the possible combinations of parameters for the original model, going from 6 to 9, enriching the dynamics of the models and matching the bistability window observed in [7] with each of the 9 possible combinations of parameters. By keeping in mind that in our models, inducer exclusion can effectively explain operon regulation, which takes into account the current challenges against the glucose-mediated repression model, our results can also be compared with some continuous models such as those of [42], based on differential equations, where bistability windows are displayed in a wide range of glucose concentrations. It is worth mentioning that the main softwares used in this manuscript were Matlab (for the exhaustive analysis of all the different dynamics of each model and which allowed us to build Table 3 and Table 5) and RStudio (for the visualization of most of the state transition graphs presented in this work as well as for the stochastic experiment; Figure 15).

Although the lac operon is a classical model and many of its key molecular players were identified a long time ago, our understanding of how these players interact with each other has evolved continuously through the years [43]. Full comprehension of the system requires robust and thorough knowledge of key regulatory features, like catabolite repression, where the view of a direct inhibitory effect of glucose on the cAMP-CRP regulatory system has been challenged during the last decades [25,26,27,28,29,30]. Furthermore, lac operon bistability is a feature that is hard to reproduce by models in general, and the applicability of a Boolean network including glucose-mediated regulatory systems has only been tried previously in [20]. Our present work shows how the translation of updated mechanistic information about the actual role of glucose into the network allows taking better advantage of a Boolean model to test and reproduce bistability in a way that is faithfully representative of experimental data. This has multiple implications as, in general, bistability has been considered a ubiquitous feature in bacteria, involving several processes [44], and understanding the dynamics of global gene regulatory systems revealed by high throughput technologies has become increasingly complex, thus demanding simpler, robust modelling algorithms, such as Boolean networks. Future models could include additional molecular features of the operon for the lac system, like repressor oligomerization, inducer (lactose) degradation by LacZ activity, and consideration of additional binding sites within the lac promoter region to provide a more detailed description of biological outcomes. Another natural future extension consists of giving a more quantitative approach to the models here studied, similar to what has been done in works such as [45,46,47].

Figure 1.(a) The lac operon OFF state is characterized by transcription of low mRNA levels, due to the strong binding of the LacI repressor to operator sequences. (b) The lac operon ON state involves transcription of high mRNA levels and increased production of the different lac proteins.

Figure 2. The original model of [20] (left) and its local Boolean functions (right). A solid (dashed) edge represents an activation (inhibition). We use the red box to differentiate the parameters of the model (outside) from its variables (inside).

Figure 3. The reduced model of [20] (left) and its local Boolean functions (right). Dashed/solid edges represent inhibitions/activations.

Figure 4. State transition graphs of the original model under parallel update schedule for: (a) Case 1, (b) Case 2, (c) Case 3 and (d) Case 4 (bistability appears). The circles represent all the210=1024configurations(M,P,B,C,R,Rm,A,Am,L,Lm)∈{0,1}10, in particular, the red one correspond to the steady state OFF while the green one correspond to the steady state ON. The sizes of its attraction basins are 1024 for cases 1, 2 and 3 and 1006 and 18 for the OFF and ON steady states of case 4, respectively.

Figure 5. State transition graphs of the reduced model under parallel update schedule. Case 4 (left): low glucose and medium lactose. Case 3 (center): low glucose and high lactose. Cases 1 and 2 (right): the other four combination of parameters.

Figure 6. Acyclic digraphs without loops obtained when considering the parameters of the cases 1 (left), 2 (left also) and 3 (right) of Section 3.2.

Figure 7. Digraph obtained when considering the parameters of case 4 (left) of Section 3.2 and its local functions after consideringGe=Le=L=A=0andLem=C=1(right).

Figure 8. State transition graph of the original model under the parameters of case 4 and considering that the first nodesM,P,B,C,R,A,L,Lmare updated (in parallel) and then the nodesRm,Am(also in parallel). The size of the attraction basins are 18 and 98 for the steady states ON and OFF, respectively (both in red), and 908 considering the four limit cycles (in blue).

Figure 9. The original model (left) and the improved original model (right), where changes are marked in red.

Figure 10. The reduced model (left) and the improved reduced model (right), where changes are marked in red.

Figure 11. Interaction digraph associated to the improved original model.

Figure 12. Dynamic of the improved reduced model associated with; low glucose and low lactose (top left), low glucose and medium lactose (top center), low glucose and high lactose (top right), high glucose and low lactose (bottom left), high glucose and medium lactose (bottom center) and high glucose and high lactose (bottom right).

Figure 13. The improved model shown in Figure 9 (right) but with the additional parameterGem. Modified local functions appear in red.

Figure 14. Interaction digraph of the improved original model with additional parameterGem and the local functions of Figure 13.

Figure 15. Bifurcation diagram for the stochastic experiment performed over the improved original model of Figure 13 and Figure 14. Blue, red and green marks represent dynamics with steady states ON, bistability (OFF and ON at the same time) and OFF, respectively.

Glucose LowGe=0 HighGe=1
Lactose
Low:(Le,Lem)=(0,0)Case2(OFF)Case1(OFF)
Medium:(Le,Lem)=(0,1)Case4(Bistability)Case1(OFF)
High:(Le,Lem)=(1,1)Case3(ON)Case1(OFF)

Attractor OFF ON
Case
1:Ge=110240
2:Ge=Le=Lem=010240
3:Ge=0∧Le=Lem=101024

Case 4:Ge=Le=0∧Lem=1 (Bistability)
AttractorsSFPLC
OFF684.7908.9153.7
ON124.4115.1146.5
Limit cycles214.90723.8
Total102410241024

Attractor OFF ON
Case
1:Ge=180
2:Ge=Le=Lem=080
3:Ge=0∧Le=Lem=108

Case 4:Ge=Le=0∧Lem=1 (Bistability)
AttractorsSFPLC
OFF2.52.81.3
ON4.65.22.7
Limit cycles0.904

Glucose LowGe=0 HighGe=1
Lactose
Low:(Le,Lem)=(0,0)OFFOFF
Medium:(Le,Lem)=(0,1)BistabilityOFF
High:(Le,Lem)=(1,1)ONBistability

Glucose Low(Ge,Gem)=(0,0) Medium(Ge,Gem)=(0,1) High(Ge,Gem)=(1,1)
Lactose
Low(Le,Lem)=(0,0)OFFOFFOFF
Medium(Le,Lem)=(0,0)BistabilityBistabilityOFF
High(Le,Lem)=(0,0)ONONBistability

Author Contributions

Conceptualization, M.M.-M. and T.L.; methodology, M.M.-M. and T.L.; software, M.M.-M. and G.A.R.; validation, M.M.-M., T.L., G.A.R. and E.G.; formal analysis, M.M.-M. and T.L.; investigation, M.M.-M., T.L., G.A.R. and E.G.; writing-original draft preparation, M.M.-M., T.L., G.A.R. and E.G.; writing-review and editing, M.M.-M., T.L., G.A.R. and E.G.; visualization, M.M.-M., T.L. and G.A.R.; funding acquisition, M.M.-M., G.A.R. and E.G. All authors have read and agreed to the published version of the manuscript.

Funding

This research was partially funded by STIC-AmSud CoDANet project 19-STIC-03 (M.M.-M., E.G.), ANID FONDECYT grant number 1200006 (E.G.), ANID FONDECYT grant number 1180706 (G.A.R), and ANID PIA/BASAL FB0002 (T.L.).

Institutional Review Board Statement

Not applicable.

Informed Consent Statement

Not applicable.

Data Availability Statement

The data used to support the findings of this study are included within the article.

Conflicts of Interest

The authors declare no conflict of interest.

1. Novick, A.; Weiner, M. Enzyme induction as an all-or-none phenomenon. Proc. Natl. Acad. Sci. USA 1957, 43, 553.

2. Cohn, M.; Horibata, K. Inhibition by glucose of the induced synthesis of the β-galactoside-enzyme system of Escherichia coli. Analysis of maintenance. J. Bacteriol. 1959, 78, 601.

3. Jacob, F.; Monod, J. Genetic regulatory mechanisms in the synthesis of proteins. J. Mol. Biol. 1961, 3, 318-356.

4. Griffith, J. Mathematics of cellular control processes II. Positive feedback to one gene. J. Theor. Biol. 1968, 20, 209-216.

5. Monod, J.; Pappenheimer, A., Jr.; Cohen-Bazire, G. La cinétique de la biosynthèse de la β-galactosidáse chez E. coli considérée comme fonction de la croissance. Biochim. Biophys. Acta 1952, 9, 648-660.

6. Deutscher, J. The mechanisms of carbon catabolite repression in bacteria. Curr. Opin. Microbiol. 2008, 11, 87-93.

7. Ozbudak, E.M.; Thattai, M.; Lim, H.N.; Shraiman, B.I.; Van Oudenaarden, A. Multistability in the lactose utilization network of Escherichia coli. Nature 2004, 427, 737-740.

8. Jenkins, A.; Macauley, M. Bistability and Asynchrony in a Boolean Model of the L-arabinose Operon in Escherichia coli. Bull. Math. Biol. 2017, 79, 1778-1795.

9. De Jong, H. Modeling and simulation of genetic regulatory systems: A literature review. J. Comput. Biol. 2002, 9, 67-103.

10. Albert, R. Boolean modeling of genetic regulatory networks. In Complex Networks; Springer: Berlin/Heidelberg, Germay, 2004; pp. 459-481.

11. Albert, R.; Othmer, H.G. The topology of the regulatory interactions predicts the expression pattern of the segment polarity genes in Drosophila melanogaster. J. Theor. Biol. 2003, 223, 1-18.

12. Mendoza, L.; Xenarios, I. A method for the generation of standardized qualitative dynamical systems of regulatory networks. Theor. Biol. Med. Model. 2006, 3, 1-18.

13. Wang, R.S.; Saadatpour, A.; Albert, R. Boolean modeling in systems biology: An overview of methodology and applications. Phys. Biol. 2012, 9, 055001.

14. Chaves, M.; de Jong, H. Qualitative Modeling, Analysis and Control of Synthetic Regulatory Circuits. In Synthetic Gene Circuits; Humana: New York, NY, USA, 2021; pp. 1-40.

15. Silva-Rocha, R.; de Lorenzo, V. The TOL network of P seudomonas putida mt-2 processes multiple environmental inputs into a narrow response space. Environ. Microbiol. 2013, 15, 271-286.

16. Timmermann, T.; González, B.; Ruz, G.A. Reconstruction of a gene regulatory network of the induced systemic resistance defense response in Arabidopsis using boolean networks. BMC Bioinform. 2020, 21, 1-16.

17. Offermann, B.; Knauer, S.; Singh, A.; Fernández-Cachón, M.L.; Klose, M.; Kowar, S.; Busch, H.; Boerries, M. Boolean modeling reveals the necessity of transcriptional regulation for bistability in PC12 cell differentiation. Front. Genet. 2016, 7, 44.

18. Richard, A.; Rossignol, G.; Comet, J.P.; Bernot, G.; Guespin-Michel, J.; Merieau, A. Boolean models of biosurfactants production in Pseudomonas fluorescens. PLoS ONE 2012, 7, e24651.

19. Kothamachu, V.B.; Feliu, E.; Cardelli, L.; Soyer, O.S. Unlimited multistability and Boolean logic in microbial signalling. J. R. Soc. Interface 2015, 12, 20150234.

20. Veliz-Cuba, A.; Stigler, B. Boolean models can explain bistability in the lac operon. J. Comput. Biol. 2011, 18, 783-794.

21. Aracena, J.; Demongeot, J.; Fanchon, E.; Montalva, M. On the number of update digraphs and its relation with the feedback arc sets and tournaments. Discret. Appl. Math. 2013, 161, 1345-1355.

22. Aracena, J.; Demongeot, J.; Fanchon, E.; Montalva, M. On the number of different dynamics in Boolean networks with deterministic update schedules. Math. Biosci. 2013, 242, 188-194.

23. Santillán, M.; Mackey, M.C. Influence of catabolite repression and inducer exclusion on the bistable behavior of the lac operon. Biophys. J. 2004, 86, 1282-1292.

24. Cohn, M. Contributions of studies on the β-galactosidase of Escherichia coli to our understanding of enzyme synthesis. Bacteriol. Rev. 1957, 21, 140.

25. Görke, B.; Stülke, J. Carbon catabolite repression in bacteria: Many ways to make the most out of nutrients. Nat. Rev. Microbiol. 2008, 6, 613-624.

26. Narang, A. cAMP does not have an important role in carbon catabolite repression of the Escherichia coli lac operon. Nat. Rev. Microbiol. 2009, 7, 250.

27. Narang, A. Quantitative effect and regulatory function of cyclic adenosine 5'-phosphate in Escherichia coli. J. Biosci. 2009, 34, 445-463.

28. Crasnier-Mednansky, M. Is there any role for cAMP-CRP in carbon catabolite repression of the Escherichia coli lac operon? Nat. Rev. Microbiol. 2008, 6, 954.

29. Notley-McRobb, L.; Death, A.; Ferenci, T. The relationship between external glucose concentration and cAMP levels inside Escherichia coli: Implications for models of phosphotransferase-mediated regulation of adenylate cyclase. Microbiology 1997, 143, 1909-1918.

30. Bettenbrock, K.; Sauter, T.; Jahreis, K.; Kremling, A.; Lengeler, J.W.; Gilles, E.D. Correlation between growth rates, EIIACrr phosphorylation, and intracellular cyclic AMP levels in Escherichia coli K-12. J. Bacteriol. 2007, 189, 6891-6900.

31. Wanner, B.L.; Kodaira, R.; Neidhardt, F.C. Regulation of lac operon expression: Reappraisal of the theory of catabolite repression. J. Bacteriol. 1978, 136, 947-954.

32. Inada, T.; Kimata, K.; Aiba, H. Mechanism responsible for glucose-lactose diauxie in Escherichia coli: Challenge to the cAMP model. Genes Cells 1996, 1, 293-301.

33. Kimata, K.; Takahashi, H.; Inada, T.; Postma, P.; Aiba, H. cAMP receptor protein-cAMP plays a crucial role in glucose-lactose diauxie by activating the major glucose transporter gene in Escherichia coli. Proc. Natl. Acad. Sci. USA 1997, 94, 12914-12919.

34. Goles, E.; Montalva, M.; Ruz, G.A. Deconstruction and dynamical robustness of regulatory networks: Application to the yeast cell cycle networks. Bull. Math. Biol. 2013, 75, 939-966.

35. Ruz, G.A.; Goles, E.; Montalva, M.; Fogel, G.B. Dynamical and topological robustness of the mammalian cell cycle network: A reverse engineering approach. Biosystems 2014, 115, 23-32.

36. Ruivo, E.L.; Montalva-Medel, M.; de Oliveira, P.P.; Perrot, K. Characterisation of the elementary cellular automata in terms of their maximum sensitivity to all possible asynchronous updates. Chaos Solitons Fractals 2018, 113, 209-220.

37. Perrot, K.; Montalva-Medel, M.; De Oliveira, P.P.; Ruivo, E.L. Maximum sensitivity to update schedules of elementary cellular automata over periodic configurations. Nat. Comput. 2020, 19, 51-90.

38. Ruivo, E.L.; Balbi, P.P.; Montalva-Medel, M.; Perrot, K. Maximum sensitivity to update schedules of elementary cellular automata over infinite configurations. Inf. Comput. 2020, 274, 104538.

39. Goles, E.; Montalva-Medel, M.; Maclean, S.; Mortveit, H. Block Invariance in a Family of Elementary Cellular Automata. J. Cell. Autom. 2018, 13, 15-32.

40. Goles, E.; Montalva-Medel, M.; Mortveit, H.; Ramirez-Flandes, S. Block invariance in elementary cellular automata. J. Cell. Autom. 2015, 10, 119-135.

41. MacLean, S.; Montalva-Medel, M.; Goles, E. Block invariance and reversibility of one dimensional linear cellular automata. Adv. Appl. Math. 2019, 105, 83-101.

42. Santillán, M. Bistable behavior in a model of the lac operon in Escherichia coli with variable growth rate. Biophys. J. 2008, 94, 2065-2081.

43. Santillán, M.; Mackey, M.C. Quantitative approaches to the study of bistability in the lac operon of Escherichia coli. J. R. Soc. Interface 2008, 5, S29-S39.

44. Veening, J.W.; Smits, W.K.; Kuipers, O.P. Bistability, epigenetics, and bet-hedging in bacteria. Annu. Rev. Microbiol. 2008, 62, 193-210.

45. Jack, J.; Wambaugh, J.F.; Shah, I. Simulating quantitative cellular responses using asynchronous threshold Boolean network ensembles. BMC Syst. Biol. 2011, 5, 1-13.

46. Samaga, R.; Klamt, S. Modeling approaches for qualitative and semi-quantitative analysis of cellular signaling networks. Cell Commun. Signal. 2013, 11, 1-19.

47. Wittmann, D.M.; Krumsiek, J.; Saez-Rodriguez, J.; Lauffenburger, D.A.; Klamt, S.; Theis, F.J. Transforming Boolean models to continuous models: Methodology and application to T-cell receptor signaling. BMC Syst. Biol. 2009, 3, 1-21.

Marco Montalva-Medel

1,*,

Thomas Ledger

1,2,

Gonzalo A. Ruz

1,2 and

Eric Goles

1

1Facultad de Ingeniería y Ciencias, Universidad Adolfo Ibáñez, Av. Diagonal Las Torres 2700, Peñalolén, Santiago 7941169, Chile

2Center of Applied Ecology and Sustainability (CAPES), Santiago 8331150, Chile

*Author to whom correspondence should be addressed.

Word count: 7284

Show less

© 2021. This work is licensed under http://creativecommons.org/licenses/by/3.0/ (the “License”). Notwithstanding the ProQuest Terms and Conditions, you may use this content in accordance with the terms of the License.

Abstract

In Veliz-Cuba and Stigler 2011, Boolean models were proposed for the lac operon in Escherichia coli capable of reproducing the operon being OFF, ON and bistable for three (low, medium and high) and two (low and high) parameters, representing the concentration ranges of lactose and glucose, respectively. Of these 6 possible combinations of parameters, 5 produce results that match with the biological experiments of Ozbudak et al., 2004. In the remaining one, the models predict the operon being OFF while biological experiments show a bistable behavior. In this paper, we first explore the robustness of two such models in the sense of how much its attractors change against any deterministic update schedule. We prove mathematically that, in cases where there is no bistability, all the dynamics in both models lack limit cycles while, when bistability appears, one model presents 30% of its dynamics with limit cycles while the other only 23%. Secondly, we propose two alternative improvements consisting of biologically supported modifications; one in which both models match with Ozbudak et al., 2004 in all 6 combinations of parameters and, the other one, where we increase the number of parameters to 9, matching in all these cases with the biological experiments of Ozbudak et al., 2004.

Details

Title
Lac Operon Boolean Models: Dynamical Robustness and Alternative Improvements
First page
600
Publication year
2021
Publication date
2021
Publisher
MDPI AG
e-ISSN
22277390
Source type
Scholarly Journal
Language of publication
English
DOI
https://doi.org/10.3390/math9060600
ProQuest document ID
2501873915
Copyright
© 2021. This work is licensed under http://creativecommons.org/licenses/by/3.0/ (the “License”). Notwithstanding the ProQuest Terms and Conditions, you may use this content in accordance with the terms of the License.