Finite state machine and Markovian equivalents of the lac Operon in E. coli bacterium

1.   Introduction

The main source of nutrition for E. coli is glucose. However, bacteria exist and grow even when it is absent, because they use lactose in its place. The genetic region of bacterial DNA that manages lactose metabolism is known as the lac operon [22], [47]. An operon is a group of regulatory elements that affect the expression of more than one genes, producing a polycistronic mRNA strand [6]. Figure 1 depicts the genetic regions that constitute the lac operon in E. coli. Three genes, lacZ, lacY and lacA, are included in the lac operon. lacZ gene produces β–galactosidase. lacY gene produces Beta–galactoside permease. lacA gene produces Galactoside acetyltransferase (GAT). Beta–galactoside permease transports lactose into the bacterium [7]. β–galactosidase cleaves lactose to give glucose and galactose [33]. The role of GAT in lactose metabolism is yet unclear. Studies of the lac operon generally exclude lacA and its product from discussion. Upstream of the lac operon, the lacI gene produces the lac repressor [19]. The lacI gene is constitutively expressed [6].

The lac operon is inactive when glucose is present for the bacterium to survive. Additionally, it remains inactive as long as the cellular environment does not have lactose to metabolize. The lac repressor (R), binds tightly to the operator O₁ (in Figure 1) and hinders RNA polymerase (RNAP) from transcribing the genes of the operon [6]. R also binds to the other two operators O₂ and O₃. However, binding of R to only O₂ or only O₃ does not deter the transcription of the lac genes [43].

The lac repressor normally occurs as a tetramer [6] and is an allosteric protein. When it binds with allolactose, it undergoes a conformational change that hinders its DNA binding site from binding to bacterial DNA [11]. Allolactose, a combination of galactose and glucose, is the inducer of the lac operon. It is produced at low levels by the few molecules of β–galactosidase that are present before induction.

The lac operon is functional when glucose is absent and lactose is present in the environment. As long as lactose remains present, β–galactosidase cleaves it into glucose and galactose to give allolactose, which keeps R from turning the operon off [6].

When glucose is absent from the system, the concentration of cyclic Adenosine MonoPhosphate (cAMP), goes up in the bacterium [32]. cAMP binds with the Catabolite Activator Protein (CAP) to create the cAMP-CAP complex. It is this complex that binds to the CAP binding sites, C₁ and C₂. One bound, the cAMP-CAP complex activates and promotes transcription through interactions with RNA polymerase [4], [25]. The repression of the lac operon due to low cAMP levels in the presence of glucose is known as catabolite repression [52]. For brevity, the cAMP-CAP complex will be referred to as CAP only.

The lac operon acts as a bi-stable switch, by existing in one of two possible states, the operon being on and off [42]. Figure 2 depicts the regulatory elements that constitute the lac operon in E. coli.

O₃ lies at the end of lacI gene and O₂ lies in the early part of lacZ gene [35]. P₂ is considered a secondary promoter, as it does not play a significant role in transcription initiation [12], [39]. Mathematical models including that presented by [43] only consider the role of P₁ in the activation of the operon. O₃ and C₁ are adjacent to each other. R bound to O₃ bends DNA in the same direction as CAP bound to C₁ does [50]. If R binds to O₃, it is not possible for CAP to bind to C₁, and vice versa. O₁ and C₂ are at the same position but on opposite sides of the DNA chain. Hence, if R binds to O₁, it is not possible for CAP to bind to C₂, and vice versa [43].

The operators O₂ and O₃ play a major role in bringing about DNA folds. This happens as a result of a lac repressor tetramer binding to more than one operators [34], [35], [43]. R may bind simultaneously to O₁ and O₂. No transcription can take place because of the attachment of the repressor to O₁. R may bind simultaneously to O₁ and O₃. This leaves no space for RNAP to bind to the promoter. R may bind simultaneously to O₂ and O₃. Although the DNA strand is folded, RNAP can bind to the promoters. As O₁ is unoccupied, transcription can take place too.

A list of possible configurations of regulatory elements of the lac operon is drawn in Santillán and Mackey [43]. The authors developed a set of states L={l0,l1,...,l18} from the list given by [43] which is shown in Table 1. States are classified on the basis of the occupancy of O₁, O₃, C₁, C₂ and P₁. Each state name is composed of three characters. An unbound entity (operator, CAP binding site or promoter) is represented by e. A bound operator is represented by r. A bound CAP binding site is represented by c. A bound promoter is represented by p. In each state name, the first character represents O₃ or C₁, the second character represents P₁, and the last character represents O₁ or C₂. fef signifies that O₁ and O₃ are bound to the same tetramer of R, f representing a fold. It must be noted that R_F signifies the repressor tetramer bound simultaneously to two operators, bringing about a DNA fold. States in which RNAP is bound to the promoter are transcriptionally active states.

2.   Existing models of the lac Operon

The lac operon has been the focus of mathematical modeling ever since the conception of computational biology. Wong et al. developed a mathematical model of the lac operon which included catabolite repression, inducer exclusion, lactose hydrolysis to glucose and galactose, and synthesis and degradation of allolactose [52]. Mahaffy and Savev proposed a mathematical model for induction of the lac operon using biochemical kinetics [27]. Suen and Jacob presented a symbolic, grammar-based model for the operon in Mathematica [47]. Yildirim and Mackey proposed a mathematical model for the regulation of induction in the lac operon [55]. Yildirim et al. derived a reduced model from an existing model to explore the bi-stability of the operon [56]. Santillán and Mackey used a mathematical model to investigate the influence of catabolite repression and inducer exclusion on the bi-stable behavior of the operon [43]. Jacob and Burleigh developed a swarm-based, 3-dimensional model of the lactose operon gene regulatory system [21]. Santillán investigated bi-stability of the operon using a mathematical model [42]. Veliz-Cuba and Stigler presented a Boolean network as a discrete model for the lac operon [49]. Angelova and Ben-Halim proposed a deterministic model of the lac operon with a noise term, representing the stochastic nature of regulation [3]. Yildirim and Kazanci revisited the Yildirim-Mackey model [55] to show how deterministic and stochastic methods can be used to investigate various aspects of the operon [54]. Esmaeili et al. built a 3-dimensional, interactive computer model of the lactose operon using PROKARYO, an agent-based cell simulator [14]. Choudhary and Narang solved a detailed model of lac regulation using singular perturbation theory in [10]. Li et al. investigated the controllability of logical control networks and applied the obtained results to the reachability of the lac operon in E. coli [26].

3.   Digital and Markov models of genetic circuits

FSMs are machines with memory [48]. A finite state machine is composed of two finite sets, S and E. S is the set of all possible states of the machine. E is the set of all possible events for the machine. At any time instant, the machine persists in one of the states of S, s_i, which is considered the current state. When an event occurs, the machine transits from s_i to the next state, s_j. The transition from s_i to s_j depends on both s_i and the event that triggered the transition [29]. FSMs are immensely useful in modeling a closed environment which accepts external stimuli and reacts accordingly, changing itself in alignment to a set of rules. They can be programmed on a field-programmable gate array (FPGA) using hardware description language (HDL).

A stochastic process has value Xn(n=0,1,2...) for each time instant n. X_n = i represents that the process is in state i at time instant n. There is a fixed probability P_ij, which signifies transition of the process from state i to j. This process is called a Markov chain if the next state depends only on the current state, and never on states attained in the past [40].

Eq 3.1 signifies that the probability of any future state depends only on the present state. Markov chains are used to simulate randomly changing events.

Muri et al. used hidden Markov models to study bacterial genomes [30]. Krogh et al. described a membrane protein topology prediction method based on a hidden Markov model [24]. Markov chains were used to ponder the effect of individual genes on global dynamical network behavior [45]. Calder et al. analyzed signal transduction networks using continuous time Markov chains [8]. Julius et al. used Markov chains to model the lactose regulation system of E. coli [23]. Wang et al. studied the dynamic properties of the regulatory network governing lytic and lysogenic growths of coliphage lambda using a Markov chain stochastic model [51]. Vergne presented a drifting Markov model for the genetic content of lambda phage [31]. Yang et al. modeled proteins with the help of molecular finite automata (MFA) [53]. A Markov chain was presented for the process of target finding by lambda phage on the surface of E. coli [9]. Gao and Hu used FSMs to study gene expression [16]. Gao et al. defined a finite state machine and a hidden Markov model for genetic mutation [17]. FSMs were created using DNA [13]. Synthetic biologists created a framework to create FSMs using gene regulatory networks [36]. Shenker and Lin proposed Markov chains for cooperative binding in E. coli's infection with lambda phage [44]. Fang et al. created a hidden Markov model for the λ switch [15]. A finite state machine and a Markov chain were created by Ainuddin et al. [1] for the lambda switch. This paper extends the same modeling technique to the lac operon.

4.   Finite state machine design

The set L in Table 1 is used to create a finite state machine which transits in response to presence and absence of glucose and lactose. Six transitions have been defined for the finite state machine, shown in Table 2. Figure 3 depicts the binding transitions of the equivalent FSM. Figure 4 depicts the unbinding transitions of the equivalent FSM. The FSM is translated to a field-programmable gate array (FPGA) with the help of Verilog hardware description language (HDL). The implementation is named Digital lac Operon (DLO). The input signals of DLO are listed in Table 3. The output signals of DLO are listed in Table 4.

When the initial state is either l₃ or l₉, and the machine transits with RNAP unbinding from the operon, out = 01. When the initial state is any one of l₅, l₁₁ and l₁₅, and the machine transits with RNAP unbinding from the operon, out = 10. When the initial state is l₁₇ and the machine transits with RNAP unbinding from the operon, out = 11. In all remaining transitions, out = 00. A module called Lac, and a permease counter are brought together in DLO. Figure 5 depicts the internal circuitry of DLO.

An 8-bit linear feedback shift register (LFSR) is used for pseudo-random number generation [2], [37]. The feedback polynomial of this LFSR is f(x)=x8+x4+1. The output byte is called rand. sel=rand[0], loc1=rand[1], loc0=rand[2], f1=rand[3] and f0=rand[4].

All outputs of DLO are delivered by the Lac module. All inputs to DLO are fed on to the Lac module. Additional inputs to the Lac module are listed in Table 5.

Consider states l₀ and l₃, when R = 1, C = 0 and O₂ is unbound. If the f bit pair is 00, R binds to any one of the available repressors. If the f bit pair is 01, R binds to both O₂ and O₃. If the f bit pair is 10, R binds to both O₁ and O₂. If the f bit pair is 11, R binds to both O₁ and O₃. Consider states l₁, l₂, l₄ – l₆, l₉, l₁₂ and l₁₅, when R = 1, C = 0 and O₂ is unbound. If f₀ = 0, R binds to any one of the available operators. Otherwise, R binds to two operators, one of which is O₂, creating a DNA fold.

Consider states l₀ and l₃, when R = 1, C = 0, O₂ is unbound and f = 00. If the loc bit pair is 00, R binds to O₃. If the loc bit pair is 01, R binds to O₁. If loc₁ = 1, R binds to O₂. Consider states l₇ and l₁₀, when R = 0, O₂ is bound and no DNA fold exists. If the loc bit pair is 00, R unbinds from O₃. If the loc bit pair is 01, R unbinds from O₁. If loc₁ = 1, R unbinds from O₂. loc₀ is used in all states of L, when two operators are available for binding, or two repressor tetramers are bound to the operon, or two CAP binding sites are available for binding, or two cAMP-CAP complexes are bound to the operon. If loc₀ = 0, the available operator closest to the lacI gene is selected for binding or unbinding. Otherwise, the available operator farthest from the lacI gene is selected. If loc₀ = 0, C₁ is selected, otherwise C₂.

The sel bit is used to select R or CAP in states l₂, l₅, l₈, and l₁₁ – l₁₇, when at least one instance of both is attached to the operon. If sel = 0, R is removed, otherwise CAP is removed.

A 3–bit memory register called mem is defined inside the Lac module. O₂ is represented by the most significant bit. If O₂ is unbound, its relevant bit is 0, otherwise, 1. If DNA is not folded, the rightmost bits are 00. If O₂ and O₃ are bound to the same R tetramer, the rightmost bits are 01. If O₁ and O₂ are bound to the same R tetramer, the rightmost bits are 10. If O₁ and O₃ are bound to the same R tetramer, this is not represented by mem, but by the state itself (l₁₈).

DLO employs an up-down counter to keep track of the permease transcribed. It is a hexadecimal counter of 4 bits. Table 6 depicts the inputs of the permease counter. Table 7 depicts the outputs of the permease counter.

The bits of out are ORed to create a signal which advances the counter when HIGH. The inlac bit decreases the counter when HIGH. If the 4–bit output counter value is greater than 0H, the output bit permease is HIGH.

The memory mem is reset at initialization. The Lac module is initialized to state l₀. The permease counter is initialized to 0H. The FSM was synthesized on Xilinx Virtex-4 FPGA using ISE Design Suite from Xilinx Design Tools. The machine runs for a simulation time of 1µs. Time resolution is 1ps. Simulation results for seed = FEH are summarized in Table 8.

The complete truth table of the DLO appears at the end of this paper as Appendix A.

5.   Markov chain for the lac operon

We define α as the probability of R binding to the operon, β as the probability of CAP binding to the operon and γ as the probability of RNAP binding to the operon. We also define the following probabilities:

Figure 6 depicts the transition matrix M for the lac operon. This is a left stochastic square matrix.

5.1.   Steady state analysis of M

A Markov model attains steady state when the probability of the system to exist in any state does not change with t or n, as n → ∞, as given by Eqs 5.4,5.5 [18]. For a scalar value ψ, eigenvector e for a matrix T satisfies Eq 5.6 [41]. We can deduce that v_ss = e when M = T and ψ = 1. The normalized eigenvector corresponding to unity eigenvalue is shown in Eq 5.7.

Table 9 details configurations of regulatory elements that lead to operon induction [43].

We refer to the steady state probability of activation as p_acss.

6.   Discussion

The inverter at input C of the Lac module is the digital manifestation of catabolite repression.

In the first 200 epochs (0–190 ps) of simulation (as shown in Table 8), the DLO is functional even in the absence of both glucose and lactose. The cell is starving, and some β–galactosidase, still present from the operon's previous activity, produces a small amount of allolactose, which inactivates the lac repressor. This cause the operon to be induced.

The top-left subplot of Figure 7 depicts the condition where γ = 0. This implies that RNAP does not bind to the promoter. When γ ≠ 0, the probability of activation always remains less than or equal to γ, as can be seen from the Figure 7. Figure 7 clearly displays that bacterial cultures with higher growth rates have their operons transitioning swiftly to full induction. This is a positive feedback.

In Figure 8, p_acss decreases with α. This means that as the concentration of the lac repressor increases, the operon transitions to the off state.

In Figure 9, p_acss increases with β. This means that as the concentration of the cAMP-CAP complex increases, the operon transitions to full induction. Figure 9 shows a highly activated operon for lower values of α, or lower concentration of R.

6.1.   Comparable results from existing literature

Wong et al. have reported similar findings regarding concentrations of the lac repressor, β–galactosidase and Beta–galactoside permease [52]. Yildirim and Mackey have reported similar findings regarding concentrations of lactose and β–galactosidase [55]. Santillán and Mackey furnish graphs conveying the same information as Figures 8,9 [43]. Esmaeili et al. produce similar results from their model of the operon [14].

7.   Conclusion

This work focuses on digital and stochastic counterparts to the lac operon and translates the process from the cytoplasm to silicon. FPGA implementation of the developed FSM is simulated and documented, creating a silicon mimetic [20], [38]. The FSM is found to bear integrity with the modeled genetic switch. A Markov chain with the same set of states as the FSM, undergoing the same set of transitions, furnishes stochastic results that are in line with experimental data.

Although the lac operon has been the subject of mathematical modeling in numerous other works, this paper sets itself apart because it converts it to a digital machine. This is an emerging approach to modeling biological circuits. Other ventures in digital modeling of the lac operon have been listed in Section 2, but the techniques used are all together different from this research. It should be noted that in [28], the circuit diagram proposed is genetic (and not electronic) in nature. In [5], the modeling approach is not Boolean; however it is logical and uses discrete variables and functions. Table 10 clearly shows the novelty of modeling presented in this work. This work helps to bring the genetic switch to electronic hardware, where it can be tested in the same ways as an electronic circuit. Although this work restricts itself to a genetic switch, the modeling strategy can be adapted to genetic processes with more than two stable behaviors.