One of the first mathematical models of the tryptophan operon was introduced by Bliss et al. in 1982 [14]. The independent variables of this model are the trp mRNA concentration (M), the concentration of anthranilate synthase (E) which, according to Bliss et al., of all the enzymes formed with the polypeptides of the tryptophan operon is the most important from a regulatory point of view, and the tryptophan concentration (T). The equations governing the dynamic evolution of these three quantities are

Bliss et al. assumed that R can be modeled by a Hill function

This model clearly considers two different regulatory mechanisms, repression and enzymatic feedback inhibition by tryptophan. R(T) accounts for repression, while I(T) stands for feedback inhibition. Bliss et al. also assert that I(T) can be modeled by a Hill function

The authors made a careful estimation of the parameters in this model based on the available experimental data. They solved the model equations numerically and also performed an analytical stability analysis of the steady states. They were able to reproduce the results of derepression experiments with cultures of wild and mutant E. coli strains reported in [15]. From these experiments, the tryptophan operon loses stability in mutant E. coli cell cultures that have a partial loss of of feedback inhibition. The previously stable steady state is replaced by an oscillatory production of tryptophan [15]. This mutation is modeled by increasing K_i to ten times its wild type value. However, to our knowledge, these experiments have never been successfully repeated. Moreover, more recent experimental evidence has demonstrated that the function R(T), which Bliss et al. used to model repression, is in disagreement with the experimental facts about the interaction between the trp operon and repressor molecules. The issue of oscillations in tryptophan control was revisited by Xiu et al. [16] in their investigation of the role of growth and dilution rates on stability.

In 1988, Sinha [17] introduced a different model for the tryptophan operon regulatory system in which the DNA-repressor interaction is modelled in a more detailed way. The independent variables of the Sinha model are the same as in the Bliss et al. model, i.e. the trp mRNA concentration (M), the anthranilate synthase enzyme concentration (E), and the tryptophan concentration (T). The dynamic equations for these variables are

When produced by the trpR operon, repressor molecules are unable to repress trp operons. For this repression to take place, they need to be activated by two tryptophan molecules which sequentially bind non-cooperatively in two independent places. From this, and assuming that this reaction takes place in a quasi-steady state, Sinha obtained the concentration of active repressor given by R_A=nTR/(T+K_d), where n=2, R is the total repressor concentration, and K_d is the repressor activation dissociation constant.

Compared with the Bliss et al. [14] model, the Sinha model considers the DNA-repressor and repressor-tryptophan interactions in a more detailed way. It also ignores the fact that anthranilate synthase is feedback inhibited by tryptophan as well as ignoring the time delays inherent to the system. These are both important features of the tryptophan operon regulatory system.

It is also important to note that Sinha assumes a constant tryptophan consumption rate G_max. Observing this, Sen and Liu [18] modified the Sinha model to study the case of a non-linear tryptophan consumption rate given by a Michaelis–Menten-like function G(T)=G_maxT/(T+K_g). Sinha, as well as Sen and Liu, studied the stability of their models for various values of parameters K_o and K₂ and found the corresponding stability regions in the (K_o,K₂) parameter space. However, they neither compare their models with experimental data nor do they discuss whether the values for K_o and K₂ where stability is lost are physiologically attainable in wild type or mutant bacterial strains.

Recently Santillán and Mackey [19] introduced a more detailed mathematical model of the tryptophan operon regulatory system, which is shown schematically in Fig. 1. In this model, all of the known regulatory mechanisms of the tryptophan operon (repression, transcriptional attenuation, and enzyme feedback inhibition) are taken into account. The system inherent time delays due to transcription and translation are also included, as well as the mRNA dynamics. The model considers four independent variables: The concentration of trp operons with operators free to be bound by mRNA polymerase molecules (O_F), the concentration of trp mRNA molecules with free TrpE-related ribosome binding sites (M_F), the concentration of the enzyme anthranilate synthase (E), and the tryptophan concentration (T).

Although five different polypeptides are synthesized during the expression of the trp operon and they combine to form the enzymes that catalyze the reaction pathway which synthesizes tryptophan from chorismic acid, this model concentrates on anthranilate synthase. The reason for this focus is that this enzyme is the most important from a regulatory point of view since it is subject to feedback inhibition by tryptophan and is the first to play its catalytic role in the reaction pathway [13].

Anthranilate synthase is a complex of two TrpE and two TrpD polypeptides. From this, Santillán and Mackey assume that the anthranilate synthase production rate is one half that of the TrpE polypeptide and look only at the TrpE-related ribosome binding sites of the mRNA chain. The equations governing the dynamics of these variables are:

The factor

μ is the growth rate of the bacterial culture. The term μO stands for the operon production rate due to DNA replication, which was assumed to be such that it keeps the total operon concentration at a constant value O, balancing dilution by growth.

The term

In Eq. (8), τ_m is the time it takes for a polymerase to produce a mRNA chain long enough to have an available TrpE-related ribosome binding site. A(T) is the probability for transcription to be prematurely terminated due to transcriptional attenuation. Santillán and Mackey [19] take A(T)=b[1−exp(−T/c)]. The term

The first term on the right-hand side of Eq. (9), is the anthranilate synthase production rate. The factor 1/2 is incorporated because, as mentioned earlier, the enzyme production rate is assumed to be one half that of TrpE polypeptide. The TrpE production rate equals the rate of ribosome binding free TrpE-related ribosome binding sites a time τ_e ago. τ_e is the time it takes for a ribosome to synthesize and release a TrpE polypeptide.

The tryptophan production rate, the first term on the right-hand side of Eq. (10), is assumed to be proportional to the concentration of active (uninhibited) enzymes (E_A) with a proportionality constant K. Anthranilate synthase is inhibited when two tryptophan molecules bind the TrpE subunits in a sequential cooperative reaction with a Hill coefficient of n_H≃1.2. The enzyme feedback inhibition reaction rate constants support a quasi-steady state assumption, from which

All of the model parameters were estimated from the available experimental data. The details are given in Reference [19] supplementary web material. In Table 1 the value of all these parameters is shown.

The model was solved numerically, and the results compared with derepression experiments performed by Yanofsky and Horn [20] on wild-type bacterial cultures as well as trpL29 and trpL75 mutant strains of E. coli. In these experiments, a bacterial culture is first grown in a medium with a high tryptophan concentration during a period of time long enough for the trp operon to reach a steady state. Then the bacteria are washed and shifted to a medium without tryptophan, and the response of the anthranilate synthase activity is measured as a function of time after the nutritional shift.

The trpL29 mutant strain has a mutation A to G at bp 29 in the leader region of the trp operon. This change replaces the leader peptide start codon by GUG, and decreases operon expression in cells growing in the presence or absence of tryptophan. This mutation is simulated by Santillán and Mackey by decreasing the rate constant k_p to 0.004 times its normal value.

The strain trpL75 of E. coli has a mutation of G to A at bp 75 in the leader region of the trp operon. This change decreases the stability of the transcription antiterminator structure, and increases transcription termination at the attenuator. Consequently, it decreases operon expression in cells growing in the presence or absence of tryptophan. This mutation was simulated by increasing the value of parameter b, which determines the probability of transcriptional attenuation at high tryptophan levels. The normal value of this parameter is b≃0.85, and to simulate the mutation it was increased to b≃0.9996.

The experiments of Yanofsky and Horn were simulated numerically. In all of the three cases (wild type and two mutants), the model results show a reasonable qualitative agreement with the experiments, given the simplifying assumptions inherent to the model. The results of these simulations are shown in Fig. 2.

The lac operon is the classical bacterial example of an inducible system which encodes the genes for the transport of external lactose into the cell and its conversion to glucose and galactose. Many experimental studies have examined the control of the genes involved in lactose metabolism in E. coli. One of the most recent mathematical modeling studies of the lac operon was carried out by Wong and coworkers [21]. Although many of the relevant biological details were taken into account in their model, they failed to treat the transcriptional and translational delays. Conversely, Mahaffy and Savev [22] modeled lac operon dynamics and included transcriptional and translational delays, but ignored the conversion of internal lactose to allolactose by β-galactoside and the spontaneous production rate of mRNA together with the individual degradation terms for the proteins. In this section we give a brief summary of the work on a recent model developed by Yildirim et al. [23] on the lac operon including much of the relevant biological detail considered by Wong and coworkers [21] as well as the transcriptional and translational delays considered by Mahaffy and Savev [22]. A more complete mathematical analysis of the model presented here can be found in [23] and [24].

The lac operon consists of a promoter/operator region and three larger structural genes lacZ, lacY, and lacA together with a preceding regulatory operon responsible for producing a repressor (R) protein. To understand the mechanism under this control network easily, it is helpful to refer to Fig. 3. In the absence of glucose available for cellular metabolism, but in the presence of external lactose (L_e), lactose is transported into the cell by a permease (P). Intracellular lactose (L) is then broken down into allolactose (A) first and then glucose and galactose by the enzyme β-galactosidase (B). The allolactose feeds back to bind with the lactose repressor and enable the transcription process to proceed. In the absence of allolactose (A) the repressor R binds to the operator region and prevents the RNA polymerase from transcribing the structural genes. However, if allolactose is present, a complex is formed between allolactose and the repressor that makes binding of the repressor to the operator region impossible. In that case, the RNA polymerase bound to the promoter is able to initiate transcription of the structural genes to produce mRNA. When a mRNA chain, long enough for a ribosome to bind it, has been produced, the process of translation is initiated. The lacZ gene encodes for the mRNA responsible for the production of β-galactosidase (B) and transcription of the lacY gene produces mRNA ultimately responsible for the production of a membrane permease (P). The mRNA produced by transcription of the lacA gene encodes for the production of thiogalactoside transacetylase which is thought to not play a role in the regulation of the lac operon [25] and will not be further considered here.

There is another control mechanism known as catabolite repression. The presence of glucose in the external medium inhibits the production of cAMP in the bacterium. cAMP is necessary for transcription initiation to take place efficiently. In the present work we ignore these regulatory mechanism since we compare with experimental in which there is no glucose in the culture medium and thus, cAMP is always produced at a high rate.

The lac operon regulatory pathway is fairly complex: it involves three repressor binding sites (operators), and two CAP binding sites. Moreover, one repressor can bind two different operators simultaneously [26]. A detailed model of all these regulatory elements is beyond the scope of this work. Nonetheless, according to Yagil and Yagil [27], the fraction of free operators as a function of allolactose (A) can be modeled via a Hill-type function of the form:

The dynamics of the concentration (B) of β-galactosidase are described by Eq. (13) in which τ_B is the time required for β-galactosidase production through mRNA translation. We assume that the rate of production of B is proportional to the concentration of M a time τ_B ago times a dilution factor due to cell growth during the translation process, $\left({\alpha }_{\mathrm{B}}{\mathrm{e}}^{-\mu {\tau }_{\mathrm{B}}}{M}_{{\tau }_{\mathrm{B}}}\right)$. The loss rate of B is given by (γ_B+μ)B.

For the allolactose concentration (A) dynamics, the first term in Eq. (14) gives the β-galactosidase mediated gain in allolactose from the conversion of lactose. The second term accounts for allolactose loss via conversion to glucose and galactose (again mediated by β-galactosidase). Since they are enzymatic reactions, both of these expressions are taken to be in Michaelis–Menten form. The last term takes into account the degradation and dilution of allolactose.

The lactose concentration (L) dynamics are more complicated and given by Eq. (15). The first term in Eq. (15) accounts for the increase of intracellular lactose through the permease facilitated transport of external lactose (L_e). The second term deals with intracellular lactose loss to the extracellular fluid because of the reversible nature of the permease mediated transport [28–31]. The third term accounts for the conversion of lactose to allolactose as well as the hydrolysis of lactose to glucose and galactose via β-galactosidase (B). The fourth term is the decrease in internal lactose concentration due to degradation and dilution.

The permease concentration (P) dynamics are given by Eq. (16). The first term reflects the assumption that permease production is directly proportional to the mRNA concentration a time τ_P in the past, where τ_P is the time it takes to translate a permease polypeptide. The exponential factor e^−μτ_P accounts for dilution of mRNA concentration due to cell growth. The second term accounts for the degradation and dilution of permease.

We have carried out an extensive search of the existing literature for data that would allow us to estimate the model parameters in Eqs. (12)–(16) which are summarized in Table 2, and the details of how we arrived at these parameter values can be found in [24].

The model as formulated in Eqs. (12)–(16) can have one, two, or three steady states depending on the values of the parameters (μ,L_e). The details of how these steady states are determined are contained in [23] and [24]. The results of these considerations are presented in Fig. 4. There we show in $({\mathrm{L}}_{\mathrm{e}},\overline{A})$ space the region where a non-negative steady state can exist. Note in particular that for a range of L_e values there may be coexisting two stable and one unstable steady state values of the intracellular allolactose levels, $\overline{A}$ and, consequently, of $(\overline{M},\overline{B},\overline{A},\overline{L},\overline{P})$.

In other words, analytical and numerical studies of the model predict that for physiologically realistic values of external lactose and the bacterial growth rate, a regime exists where there may be a bistable steady state behavior, and that this corresponds to a cusp bifurcation in the model dynamics. This means that, as the external lactose concentration L_e is increased, the bacterial culture must switch from a low-internal-allolactose-concentration uninduced steady state to a high-internal-allolactose-concentration induced steady state. For very low L_e values, only the uninduced state exists. Conversely, for very high L_e concentrations, only the induced state is available for the system. All these facts are in agreement with the experimentally-observed dynamic behavior of E. coli cultures. On the other hand, the model also predicts a range of L_e values in which both the uninduced and the induced stable steady states coexist. If this bistability behavior is real, an hysteresis phenomenon must be observed in inducing and uninducing the lac operon of a given bacterial culture by changing the external lactose concentration. This prediction is qualitatively confirmed by the observations in [32] and [33].

Although a full analysis of the stability properties of the model is not possible due to its complexity, we have found that the basic properties contained in a reduced version of this model [23] (which considered only the dynamics of M, B, and A) are apparently retained in this much more complicated model as far as we are able to ascertain analytically and numerically.

Given the parameter values determined in Table 2, we numerically solved the model equations to compare the predicted behavior with three distinct experimental data sets. The first data set is from [34]. In there, the activity of β-galactosidase in an E. coli culture (strain ML 30) was measured after the bacteria were switched from a glucose-rich non-lactose to a non-glucose lactose-rich medium. The second data set is from [35]. In this paper, Pestka et al. studied specific inhibition of translation of single mRNA molecules and gave data for the specific activity of β-galactosidase versus time for E. coli 294 in the presence of IPTG. These two data sets and the model simulation determined using Matlab's dde23 [36] routine are shown in Fig. 5.

For this simulation, initial values for the variable were chosen close to the steady state values when L_e=8.0×10⁻² mM. (With this value of L_e there is a single uninduced steady state.) To compare these two sets of experimental data with the model simulation predictions, the data were scaled so the steady-state values of measured β-galactosidase activities and those produced by the simulation were equal. As seen in Fig. 5, there is relatively good agreement between both sets of experimental data and the model predicted temporal approach of β-galactosidase activity to its steady-state value.

In the experiments of Pestka et al. [35], the bacterial culture medium has no glucose all of the time. This means that catabolite repression is not a factor in the regulation of the lac operon and thus, the comparison with our simulations is straightforward. This is not the case with the experiments of Knorre [34] since, in such case, the bacterial were originally growing in a medium rich in glucose. Notice however that the time course of both experimental data is quite similar. This does not necessarily justifies the comparison of such experimental data with our model results, but highlights the necessity of including catabolite repression in further mathematical models of the lac operon in order to test its influence on the system dynamic behavior.

As a third test of the model, a data set from Goodwin [37] was used. In this paper, the dynamic behavior of β-galactosidase was studied in chemostat cultures of E. coli synchronized with respect to cell division by periodic phosphate feeding at a period equal to the bacterial doubling time. Experimentally, oscillations in β-galactosidase concentration were observed with a period equal to the feeding period. Since in these experiments, there is no glucose present in the environment, we can attempt to simulate them with our model.

To mimic the periodic phosphate feeding in our simulation, we assumed that the bacterial growth rate varies as a function of time in manner given by

The maximal value of μ(t) corresponds the times that phosphate was added. The minimal value of this function represents the minimal amount of nutrient left in the vessel. Assuming no nutrient is left ϵ minutes before the addition of phosphate and letting ϵ→0, α can be estimated if the doubling time (T) is known: $\alpha \simeq \overline{\mu }/\mathrm{T}$. In Fig. 6, we compare the data on β-galactosidase activity from the forced culture as a function of time with the model predictions for a feeding period of T=100 min which is the population doubling time. Again, there is a reasonable good agreement between the model predicted curve and this dataset.