A Molecular Dynamic Model of Tryptophan Overproduction in Escherichia coli

Abstract

Several deterministic models simulate the main molecular biology interactions among the numerous mechanisms controlling the dynamics of the tryptophan operon in native strains. However, no models exist to investigate bacterial tryptophan production from a biotechnological point of view. Here, we modified tryptophan models for native production to propose a biotechnological working model that incorporates the activity of tryptophan secretion systems and genetic modifications made in two reported E. coli strains. The resultant deterministic model could emulate the production of tryptophan in the same order of magnitude as those quantified experimentally by the genetically engineered E. coli strains GPT1001 and GPT1002 in shake flasks. We hope this work may contribute to the rational development of biological models that define and include the main parameters and molecular components for designing and engineering efficient biotechnological chassis to produce valuable chemicals.

1. Introduction

L-tryptophan (tryptophan) is an essential amino acid for humans, being one of the three aromatic amino acids (AAAs) that are part of proteins (proteinogenic amino acids) [1]. It is a neutrally charged molecule whose final biosynthesis in bacteria and plants involves condensing an indole molecule and an L-serine. The three aromatic amino acids, L-phenylalanine, L-tyrosine, and L-tryptophan, share a common biosynthetic pathway known as the AAAs general pathway or the pre-chorismate pathway [2]. This common precursor pathway begins with D-erythrose 4-phosphate and phosphoenolpyruvate. It ends with synthesizing chorismate, which gives place to any of the three AAAs [3]. The specific pathway of tryptophan biosynthesis is highly regulated, so the higher the intracellular concentration of tryptophan, the stronger the inhibition of tryptophan synthesis. Three mechanisms are sensitive to the intracellular concentration of tryptophan and exert specific inhibitions: the transcription factor tryptophan repressor (TrpR), transcription attenuation (att), and enzyme feedback inhibition [4] (Figure 1). This over-regulation is probably related to the chemical complexity of the tryptophan molecule, which, requires considerable resources to be biosynthesized [5].

The tryptophan operon is one of the most extensively studied biological control systems, maybe only behind the lac operon. Several models describe the dynamics of the tryptophan operon [6,7,8,9,10,11,12]. Most of them include one or several of the regulatory mechanisms of the tryptophan operon, while few include the dynamics of tryptophan transport from and to the extracellular milieu [13,14,15,16]. Outstandingly, the models proposed by Santillán and Mackey [17] and Santillán and Zeron [18] include all the known mechanisms of regulation of the tryptophan operon. Both models successfully reproduced the enzymatic activity of anthranilate synthase during derepression experiments. The anthranilate synthase activity, TrpE, represents the six reactions for tryptophan synthesis, starting from chorismate. This assumption is reasonable since TrpE is the only enzyme subject to feedback inhibition and thus is the bottleneck of the biosynthetic pathway. Few models include elements for importing tryptophan from the media [7,11] or their secretion from the cell to the milieu [13,16]. The authors usually utilize a previously proposed expression to represent the transport dynamics [6]. Until very recently, we did not know the magnitude of the contribution of passive and active transport to the overall flux of tryptophan. Based on the thermodynamic and chemical characteristics of tryptophan, authors speculated that tryptophan should be able to passively cross the cell membrane and that such flux could be considerable [19,20,21,22,23]. However, more recently, new information suggests that tryptophan transport is a secondary active process [24,25,26].

Thus, in this work, we adopted one of the existing models to account for the dynamics of tryptophan expulsion for overproducing strains. We further adapted the model to include genetic modifications of a couple of overproducing E. coli strains. In such a way, we reproduced the same order of magnitude in the final tryptophan production that was experimentally observed. Below, we briefly address the natural mechanisms of tryptophan regulation. Then, we explain the central elements included in the model proposed by Santillán and Zeron [18] (our base model). Afterward, we describe the approaches and modifications applied to this base model. Finally, we show the results compared to experimentally reported data and make some concluding remarks.

2. Materials and Methods

2.1. The Tryptophan Operon

In E. coli, the tryptophan operon is subject to negative control by three different mechanisms whose grade of repression is proportional to the intracellular tryptophan concentration, either directly or indirectly. The TrpR transcription factor exerts trp operon repression by blocking the beginning of its transcription. TrpR is a homodimer of 25 kDa each that is unbinding to DNA unless they bind to a pair of tryptophan molecules, which act as corepressors [27]. At least five different transcriptional units make up the TrpR regulon, namely trpLEDCBA, mtr, aroL, aroH, and trpR [28]. The homodimer concentration within the cell varies from 120 to 375 dimers per cell at any moment, capable of repressing the expression of the trp operon by about 80-fold when tryptophan is >10 µM [29,30].

The second element regulating tryptophan production is transcriptional attenuation. Kasai introduced this term to refer to a new regulatory mechanism found in the leader region of his operon [31]. By that time, Yanofski also had trouble explaining how it was possible that mutants lacking TrpR still exhibited severe repression over the trp operon. He arrived at the same conclusion as Kasai: there had to be a “regulated transcription termination site…” in the leader region of these operons [32]. The trp operon has a leader sequence of 162 base pairs preceding the structural genes. Within the leader region, another “attenuator“ subregion is responsible for forming the terminator and anti-terminator structures. These secondary structures are mutually exclusive when the mRNA grows during the transcription of the leader region. These loop mechanisms are sensitive to the concentration of charged tRNA^Trp within the cell so that transcription halts in the order of eight folds under circumstances of unlimited tRNA^Trp. Together, transcription and attenuation repress the trp operon expression by about 600 folds [33].

The third element regulating tryptophan biosynthesis is anthranilate synthase feedback inhibition [34]. The terminal pathway of tryptophan biosynthesis begins with this enzyme. Although this is the first reaction of the path, it occurs in conjunction with the second reaction, which consists of the conversion of anthranilate to N-(5-phosphoribosyl)-anthranilate. The same heterotetramer comprises four polypeptides, two units of each trpE and trpD [35]. The tryptophan molecule induces structural changes affecting the aggregate, thus inhibiting the anthranilate synthase at a concentration of >10 µM tryptophan.

2.2. The Tryptophan Biosynthesis Mathematical Models

As previously mentioned, all the models proposed to date that describe the dynamics of the tryptophan operon do not consider or simplify tryptophan transport. Thus, we started working with a model previously proposed by Santillán and Zeron [18]. The overall model comprises three differential Equations (1)–(3), with Equation (2) containing a slight delay (M_τE) (see below).

$$\frac{dM}{dt}=k_M\left[O_{tot}\right]\frac{\frac{P}{K_P}}{1+\frac{P}{K_P}+\frac{R_{tot}}{K_R}{\left(\frac{T}{T+K_T}\right)}^2}\frac{1+2\alpha \frac{T}{K_G+T}}{{\left(1+\alpha \frac{T}{K_G+T}\right)}^2}-\left({\gamma }_M+\mu \right)M$$

(1)

$$\frac{d{E}_{tot}}{dt}=\frac{1}{2}{k}_E{M}_{\tau E}-\left({\gamma }_E+\mu \right)E_{tot}$$

(2)

$$\frac{d{T}_{tot}}{dt}=k_{T}E-\rho \frac{T}{K_{\rho }+T}-\mu T_{tot}$$

(3)

$$E=E_{tot}{\left(\frac{K_I}{T+K_I}\right)}^2$$

(4)

Equation (1) models the change in concentration of mRNA from the tryptophan operon. In Component A in Equation (1), “M” is the concentration of mRNA from trpE. k_M is the rate at which RNA polymerase (RNAP) molecules join the promoter to start transcription. O_tot is the total concentration of operator-promoter regions. P represents the concentration of RNAP. K_P is the dissociation constant of the binding reaction between RNAP and the tryptophan operon promoter. R_tot is the total concentration of active repressor molecules. T is the intracellular tryptophan concentration. K_T is the dissociation constant of the binding reaction between TrpR and its corepressor tryptophan. O_tot, P, and R_tot are assumed to be constant. Thus, Component B in Equation (1) accounts for TrpR repression of the tryptophan operon and the occupation of the tryptophan operon promoter by RNAP molecules. Component C in Equation (1) was formulated by Santillán and Zeron [18] from theoretical considerations and accounts for transcription attenuation. α is the proportion of the concentrations of total tRNA^trp and the dissociation constant of the binding reaction between a charged tRNA^trp and one of the tryptophan codons in the trpL chain. K_G is the dissociation constant of the binding reaction between tRNA and tryptophan. Lastly, in Component D in Equation (1), γ_MM is the rate of mRNA degradation and μM is the dilution of mRNA due to cell volume variation.

Equation (2) models the change in intracellular TrpE concentration. E represents the concentration of TrpE. The “1/2” comes from the fact that two polypeptides catalyze the first reaction of the pathway. Hence, modeling the dynamics of one of those polypeptides (TrpE) accounts for half of the enzymes required to catalyze the first reaction. K_E is the translation initiation rate of trpE. Since the half-life of trpE mRNA is 1 min, M_τE is the concentration of mRNA with a 1 min delay. This delay represents the time in which each mRNA molecule can be translated before degrading. γ_EE_Tot is the rate of TrpE degradation and μE_tot stands for the dilution of the enzyme. In Equation (3), T_tot is the concentration of intracellular tryptophan. k_T is the rate of tryptophan production per anthranilate synthase. This model assumes that tryptophan production is proportional to the availability of free TrpE. r is the maximum rate of metabolic tryptophan consumption, as Bliss initially determined [36]. K_r is analogous to the Michaelis constant and is equal to or less than one-tenth of the intracellular concentration of tryptophan at a steady-state [13]. μT_tot is the dilution of tryptophan due to variations in the intracellular volume. Finally, Equation (4) describes the concentration of free, active, anthranilate synthase (E), i.e., the proportion of enzymes that have not been inhibited by interacting with tryptophan molecules.

2.3. Model Implementation

We started by implementing the model proposed by Santillán and Zeron [18] in MATLAB R2020b (education license) and named this our base model. To solve the system of equations, we used the fourth order Runge–Kutta method with a step size of 0.01. The complete implemented code in this work is available in File S1.

2.4. Inclusion of TnaAB Enzymes and Transmembranal Transporters in the Model

We modified the base model code to fit it with the model proposed by Orozco-Gómez et al. [37], mainly to introduce the function of the tnaAB operon. TnaA is the tryptophanase that degrades tryptophan and TnaB is an importer system that introduces tryptophan from the milieu. Additionally, we incorporated the collective activity of the various endogenous tryptophan export systems. We considered that they work with the mixed-activation model consistent with those proposed by Cornish-Bowden [38,39]. In addition, we included the effects of increased carbon flux in the aromatic amino acid pathway through the (C/(C + k_mc)) term (see below).

2.5. Adjusting Parameters of Trp Transmembranal Transport

We determined that the effects predicted by the model refer to the degradation and export of tryptophan by the transmembrane transport fluxes and the activity of the TrpE, TrpR, and TnaA enzymes. An additional differential equation was added to the model to quantify changes in the concentration of external tryptophan (see below).

2.6. Model Validation with Experimental Reports

To simulate tryptophan synthesis, we used data reported for the GPT1001 and GPT1002 strains [40] by adjusting the respective equations. The GPT1001 E. coli strain was subjected to various genetic modifications to increase tryptophan production. To estimate the cell concentration in the cultures of GPT1001 and GPT1002, we made use of the cell calculator (OD600 to Bacterial Cell Number) provided by Agilent^® (Santa Clara, CA, USA) (https://www.agilent.com/store/bioCalcs.jsp#ncc (Last accessed on 19 March 2021)), supplemented with data from Sezonov et al. [41]. From these data, we calculated the specific growth rate for the culture of each mutant strain according to the methodology used by Maier [42]. The code of the model can be reached at: https://gitlab.com/amalib/trp-biot-dynmod (Last accessed on 19 March 2021).

3. Results

3.1. Modifications to the Based Model

The first change introduced to the base model was to replace Component B in Equation (1), which accounts for tryptophan attenuation, with an expression that represents a model of attenuation already used by Santillán and Mackey [17] (Equation (5)). Both models return similar outputs, but using the term in Equation (5) returns results slightly closer to those experimentally observed. Equation (5) represents the probability that the transcript is attenuated as a function of the concentration of intracellular tryptophan; hence, the likelihood that the tryptophan operon will be fully transcribed is equal to 1-A(T). The constants “b” and “c” have approximate values of 0.85 and 0.04 μM, respectively. The following modification introduces the contributions of the tnaAB operon to the model, i.e., the effects of tryptophanase and the tryptophan-uptake permease TnaB. So, we propose the equations proposed by Orozco-Gómez et al. [37], describing the tnaAB operon dynamic. That original model consists of six equations. The tnaAB operon is subject to catabolite repression via extracellular glucose, while intracellular tryptophan induces their expression. Furthermore, tryptophanase activity is subject to forming aggregates (foci) inside the cell [43]. Thus, Equation (6) describes the inhibition of the tnaAB operon by extracellular glucose. Equation (7) is analogous to Equation (6); however, it describes tnaAB induction by intracellular tryptophan. Equations (8) and (9) illustrate the behavior of tryptophanase, including the fraction of active tryptophanase, that is, tryptophanase out of the foci. Equations (10) and (11) represent the rate of change in the concentration of tryptophanase and TnaB, respectively. All the values for the constants used in this study are shown in

Table 1. Value parameters that are used in the model.

Estimated by	Kinetic Constants	Rate Constants	Dimensionless Constants
Santillán and Zeron [18]	Rtot ≈ 0.8 μM	ρ ≈ 2.4 × 102 μM min−1
	KP ≈ 4.5 × 10−2 μM	KM ≈ 5.1 min−1
	KI≈ 4.1 μM	γM ≈ 0.69 min−1
	P ≈ 3 μM	kE ≈ 30 min−1
	Kρ ≈ 10 μM	γE ≈ 0 min−1
	Otot ≈ 4 × 10−3 μM	kT ≈ 7.4 × 104 min−1
	KR ≈ 2 × 10−4 μM
	τE ≈ 1 min
	KT ≈ 40 μM
Santillán and Mackey [17]	c ≈ 0.04 μM		b ≈ 0.85
Santillán et al. [45]	KG2 = 2.6 μM	# kmlac = 0.18 min−1	Pp = 0.127
			kpc = 30
			nh = 1.3
Orozco-Gómez et al. [37]	KG =11 μM		nG = 4
			nW = 4
			K1 = 14
			λ = 7
			K2 = 9
This work	C ≈ 8.5 μM	ktnaA ≈ 8.7066 min−1
	K’x ≈ 168 μM	ktnaB ≈ 1750 min−1
	K2′x = 762.6 mM	Vmax ≈ 161.96 mM min−1
	K3′x ≈ 135.3 mM	Vmax2 ≈ 522.2 M min−1
	Kx ≈ 700 μM	Vmax3 ≈ 21,609 mM min−1
	K2x ≈ 41 mM
	K3x ≈ 16.4 mM
	km1′ ≈ 0.2 μM
	km2′ ≈ 10 mM
	km3′ ≈ 12.3 mM
	OP = 0.0125 μM
	OAB ≈ 4 × 10−3 μM
	KW = 60 μM

. The γ_A, γ_B, K_A, and K_B values are unknown; hence, we estimated them. γ_A and γ_B are the degradation rates of tnaAB. Since it is a single transcriptional unit, both values must be equal; therefore, γ_A = γ_B. According to Orozco-Gómez et al. [37], (γ_A + μ) ≈ μ; thus, γ_A = γ_B ≈ 0. With similar reasoning, K_A = K_B. Both constants are unknown and comparable to the transcription product and the translation rates of tnaAB. The wild-type promoter of the tryptophan operon is relatively strong. Thus, we suppose that the values of K_A and K_B should not be higher than that of the tryptophan operon, which turns into an upper limit. Bhartiya et al. [11] adjusted an equivalent value for the tryptophan operon to 65 min⁻¹; consequently, K_A = K_B ≤ 65 min⁻¹. Orozco-Gómez et al. [37] set the value of K_w to 60,000 μM; however, to fit the simulation values with the flux values through TnaB reported by Zhao et al. [44], the value of K_w was set to 60 μM.

$$A\left(T\right)=b\left(1-e^{\frac{-T}{c}}\right)$$

(5)

$$P_G{G}_e=\frac{K_G^{n_G}}{K_G^{n_G}+G_E^{n_G}}$$

(6)

$$P_{w}w=\frac{T^{n_w}}{K_w^{n_w}+T^{n_w}}$$

(7)

$$x=\lambda G_e+W_e$$

(8)

$$P_A\left(G_e,W_e\right)=1-{(\frac{x}{K_1})}^3 e^{\frac{-x}{K_2}}$$

(9)

$$\frac{dA}{dt}=K_A{O}_{AB}{P}_G{G}_e{P}_{w}w-\left({\gamma }_A+\mu \right)A$$

(10)

$$\frac{dB}{dt}=K_B{O}_{AB}{P}_G{G}_e{P}_{w}w-\left({\gamma }_B+\mu \right)B$$

(11)

As previously mentioned, the original model assumes that tryptophan production is only conditional on the concentration of active TrpE. However, Michaelis–Menten kinetics indicate that enzymatic activity also depends on substrate concentration: chorismate, in this case. Thus, the first element in Equation (3) (k_TE) should be further multiplied by (C/C + k_mc), where C is the chorismate concentration, and k_mc is the Michaelis–Menten constant of the reaction carried out by TrpE (hereafter represented as (C/C + k_mc) = ω (Equation (12)). The value of k_mc is 5.5 µM, as previously estimated [46]. No model has been proposed for the dynamics of chorismate concentration. Accordingly, for simplicity, we decided to keep “C” as a constant, indicating the steady-state concentration of chorismate. To the best of our knowledge, that value is also unknown; however, a value that fits the results with experimental data is C = 8.5 µM. Substituting “C” and k_mc for their numerical values, ω = 0.607.

Furthermore, it is necessary to add transport elements to the original model. We added the elements of transport used by Santillán and Mackey [17] and Sun et al. [16] to the model (Equation (12), components A and B). Water transport flux through TnaB was modeled using pure Michaelis–Menten kinetics (Equation (12), Component C). TnaB’s tryptophan transport (β) rate has not been experimentally defined; however, a value of 1750 min⁻¹ fits the experimental data. A similar gap happens with the tryptophanase rate for indole production. However, judging by the specific activity reported for commercial tryptophanase (t.ly/g37Q), this value (k_tnaA) should be less than the maximum rate of metabolic tryptophan consumption; the upper bound should be 240 µM min⁻¹. Consequently, k_tnaA ≤ 240 µM min⁻¹. The original model does not include any term or equation accounting for extracellular tryptophan dynamics; therefore, we added Equation (13) to the model. This equation represents a mass balance of tryptophan in the extracellular media; it includes all transport elements in Equation (12). Signs are inverted because of a loss of tryptophan in bacteria (Equation (12)) and represent a gain of tryptophan for the extracellular media (Equation (13)). Even with all the elements previously mentioned, the model cannot reproduce an extracellular tryptophan production of 0.58 µM, a value previously reported for wild-type E. coli [40]. Thus, we needed to consider additional transport elements to be included in the model. Modeling transport systems that follow pure Michaelis–Menten kinetics could not reproduce experimental data, and systems modeled uncompetitively inhibited transport systems, as did components A and B in Equation (12). We introduced three terms accounting for transport systems displaying mixed-activation kinetics [38,39] in Equations (12) and (13). Components D and E are activated by extracellular tryptophan, while Component F is activated by intracellular tryptophan. Only in this way was it possible to reproduce the experimental results. The values of the constants in Components D, E, and F were determined by trial and error, looking for the result that best fits the experimental observations (Table S1).

$$\begin{array}{cc}\frac{d{T}_{tot}}{dt}=k_{T}E\mathsf{\omega }-\hfill & \rho \frac{T}{K_{\rho }+T}+d\frac{W_e}{e+W_e\left(1+\frac{T}{f}\right)}-d\frac{T}{e+T\left(1+\frac{W_e}{f}\right)}-k_{TnaA}{P}_A\left(G_e,W_e\right)AT+k_{TnaB}\beta \frac{W_e}{W_e+K_{mb}}\hfill \\ & -V_{max}^{\prime }\frac{T}{T\left(1+\frac{K_x^{\prime }}{We}\right)+k{m}_1{}^{\prime }\left(1+\frac{K_x}{We}\right)}-V_{max2}^{\prime }\frac{T}{T\left(1+\frac{K{2}_x^{\prime }}{We}\right)+k{m}_2{}^{\prime }\left(1+\frac{K{2}_x}{We}\right)}\hfill \\ & -V_{max3}^{\prime }\frac{T}{T\left(1+\frac{K{3}_x^{\prime }}{T}\right)+k{m}_3{}^{\prime }\left(1+\frac{K{3}_x}{T}\right)}-\mu T_{tot}\hfill \end{array}$$

(12)

$$\begin{array}{cc}\frac{d{W}_e}{dt}=-d\hfill & \frac{W_e}{e+W_e\left(1+\frac{T}{f}\right)}+d\frac{T}{e+T\left(1+\frac{W_e}{f}\right)}-k_{tnaB}B\frac{W_e}{W_e+K_{mb}}+V_{max}\frac{T}{T\left(1+\frac{K_{We}^{\prime }}{We}\right)+k{m}_1{}^{\prime }\left(1+\frac{K_{We}}{We}\right)}\hfill \\ & +V_{max1}\frac{T}{T\left(1+\frac{K{1}_{We}^{\prime }}{We}\right)+k{m}_2{}^{\prime }\left(1+\frac{K{2}_{We}}{We}\right)}+V_{max2}\frac{T}{T\left(1+\frac{K{1}_T^{\prime }}{T}\right)+k{m}_3{}^{\prime }\left(1+\frac{K{2}_T}{T}\right)}\hfill \end{array}$$

(13)

In Equation (12), A refers to $d\frac{W_e}{e+W_e\left(1+\frac{T}{f}\right)}$; B refers to $d\frac{T}{e+T\left(1+\frac{W_e}{f}\right)}$; Tryptophanase refers to $k_{TnaA}{P}_A\left(G_e,W_e\right)AT$; C refers to $k_{TnaB}\beta \frac{W_e}{W_e+K_{mb}}$; D refers to $V_{max}^{\prime }\frac{T}{T\left(1+\frac{K_x^{\prime }}{We}\right)+k{m}_1{}^{\prime }\left(1+\frac{K_x}{We}\right)}$; E refers to $V_{max2}^{\prime }\frac{T}{T\left(1+\frac{K{2}_x^{\prime }}{We}\right)+k{m}_2{}^{\prime }\left(1+\frac{K{2}_x}{We}\right)}$; F refers to $V_{max3}^{\prime }\frac{T}{T\left(1+\frac{K{3}_x^{\prime }}{T}\right)+k{m}_3{}^{\prime }\left(1+\frac{K{3}_x}{T}\right)}$.

Then, we chose two overproducer E. coli strains as a reference to adapt the equations developed so far. These strains were named GPT1001 and GPT1002; their tryptophan production, rate of operon expression, specific growth rate, and other characteristics were reported by the authors of [40]. The strain GPT1001 carries a deletion of the trpR gene (ΔtrpR). Furthermore, it was transformed with the plasmid pTAT, containing extra copies of the genes trpE, aroG, and tktA. In vitro mutagenesis was also completed to modify the TrpE and AroG proteins to eliminate enzyme feedback inhibition. All of these genes were controlled by the lactose operon promoter (lac). Researchers made further modifications to this strain. The ptsG and tnaA genes were deleted. In addition, the expression of genes that catalyze pathway precursors, such as PRPP (phosphoribosyl pyrophosphate) and L-serine, were increased. Finally, the product of the ptsG gene, the IIBC subunit of the PTS system, consumes around 50% of the PEP produced by the cell [47]. Therefore, in this strain, ptsG elimination is accompanied by an increase in the expression of tktA, which increases the flux through the general aromatic amino acid pathway. The strain GPT1002, in addition to all the characteristics described for GPT1001, carries a deletion of the sequence encoding the tryptophan attenuator. Moreover, the native promoter of the genomic tryptophan operon of GPT1002 was replaced by a tandem arrangement of five Ptac synthetic strong promoter sequences.

These modifications were introduced to the dynamic model as follows. We simulated the increased flux through the pathway by increasing the chorismate concentration (Component C in Equation (12)) from 8.5 to 297.5 µM. That resulted from assuming that PEP that was not being spent by the PTS system would be redirected to the pre-chorismate pathway. The plasmid pTAT was built using pCL1920 as the backbone; according to their creators, this plasmid contains the lac operon promoter [48]. Consequently, the lac promoter controls the genes cloned into the pTAT plasmid. Thus, we introduce the effects of the pTAT plasmid to the original model through a set of equations that were previously proposed [45], which model the dynamics of the lac operon (Equations (14) and (15)). Equation (14) models the carbon catabolite repression, exerted by extracellular glucose, affecting the expression of the lac promoter. Such repression is related to the catabolite activator protein (CAP) effects. Equation (15) is a Hill function describing a reaction inhibited by glucose. Equations (14) and (15) were introduced in Equation (1), leading to Equation (16). K_mlac is the rate of lac promoter expression, which is 0.18 min⁻¹ [49]. O_P stands for the concentration of pTAT. According to Lerner and Inouye [48], E. coli cells transformed with pCL1920 contain five copies per cell. Thus, O_P ≈ 0.0125 μM. PR (A) in Equation (16) is a term initially introduced to account for the probability of the lac operon not being repressed by the lac repressor. However, this term was originally modeled by Santillán [45] to describe the behavior of a wild-type lac promoter, that is, a promoter including an operator sequence together with two operator sequences. However, pTAT lacks operator sequences. Therefore, this term was simplified considering the following: (1) A wild-type lac promoter is repressed by a factor of 1300 due to the effects of the lac repressor. (2) Tryptophan production by GPT1001 and GPT1002 was carried out under the presence of IPTG; according to previous results, lac promoter activity is increased by a factor of 6 under IPTG concentrations similar to those used by Gu et al. [40]. (3) Repression exerted over a lac operon lacking the two pseudo-operator sequences is 1.4% of that of the wild-type operon. Under the culture conditions described by the author, the lac promoter in pTAT is repressed by a factor of 3 (1300 × 0.014 × 1/6). Thus, the probability that the lac operon will not be repressed is reversed to the repression factor; here, PR (A) has a constant value of 1/3.

$$PD\left(G_e\right)=\frac{{{\rm P}}_p\left(1+P_c\left(G_e\right)\right)\left(k_{pc}-1\right)}{1+{{\rm P}}_p{P}_c\left(G_e\right)\left(k_{pc}-1\right)}$$

(14)

$$P_c\left(G_e\right)=\frac{K_{G2}^{nh}}{K_{G2}^{nh}+G_e^{nh}}$$

(15)

$$\frac{dM}{dt}=k_m\left[O_{tot}\right]\frac{\frac{P}{K_P}}{1+\frac{P}{K_P}+\frac{R_{tot}}{K_R}{\frac{T}{T+K_T}}^2}\left(1-A\left(T\right)\right)+k_{mlac}\left[O_P\right]P_D\left(G_e\right)PR\left(A\right)-\left({\gamma }_M+\mu \right)$$

(16)

Equation (16) represents the rate of change in trpE mRNA concentration for the GPT1001 mutant. In the case of mutant GPT1002, Equation (16) was further modified, eliminating the term for transcription attenuation (1-A(T)). Then, the effect of the 5Ptac promoter tandem was introduced, considering that the author reported an overall two-fold increase, compared to GPT1001, in the expression of the entire transcriptional unit of the tryptophan operon after such modification. TrpE production in GPT1001 is equal to the product of k_m and (1 − A(T)). A(T) is 0.85 for almost any concentration of tryptophan, except for tryptophan concentrations near starvation. Thus, the actual value of expression of the tryptophan operon in GPT1001 is 0.765 min⁻¹ (15% times 5.1 min⁻¹). Therefore, we replaced the original value of k_m for GPT1002 by two times 0.765 (Equation (17)).

$$\frac{dM}{dt}=1.53 \left[O_{tot}\right]\frac{\frac{P}{K_P}}{1+\frac{P}{K_P}+\frac{R_{tot}}{K_R}{\frac{T}{T+K_T}}^2}+k_{mlac}\left[O_P\right]P_D\left(G_e\right)PR\left(A\right)-\left({\gamma }_M+\mu \right)$$

(17)

3.2. Comparing the Model to Experimental Data

In a complex medium, Gu et al. [40] reported a final tryptophan titer of 0.58 μM by wild-type E. coli. However, they did not register the culture-specific growth rate and cell concentration. Thus, the optical density reported elsewhere [50] for this same E. coli growing on the same complex media and with the same conditions (Table S1) was used to calculate the culture-specific growth rate and cell concentration using the Agilent^® “OD₆₀₀ to bacterial cell number” calculator (https://www.agilent.com/store/bioCalcs.jsp#ncc (Last accessed on 19 March 2021)). We calculated the specific growth rate from the cell-density data sample points [41]. These data in Supplementary Tables S1–S3 were introduced into the model to obtain growth and glucose consumption curves. Parts 8–10 of the code in Supplementary File S1 state how we include the specific growth rates to simulate tryptophan production (also shown in the final parts of the program codes for GPT1001 and GPT 1002 at (https://gitlab.com/amalib/trp-biot-dynmod (Last accessed on 19 March 2021)). The system representing a wild-type strain growing on 50 mL media consists of Equations (2), (3) and (10)–(13). We simulated the culture for 8 h (Figure 2), which is the time in which culture in complex media reaches a steady state (Table S1). As seen in Figure 2, the final tryptophan production indicated by the model is close to 0.58 μM. It was remarked that we could only reproduce exogenous tryptophan dynamics by introducing transport systems that follow mixed-activation dynamics. When simulations were carried out using transport systems following pure Michaelis–Menten dynamics or uncompetitively inhibited dynamics, the system displayed high-frequency oscillations even at steady state. At the same time, external tryptophan concentration never reached more than 0.1 μM without showing high-frequency oscillations (not shown).

After introducing a mixed-activation transport, we simulated the models describing tryptophan overproducing in mutants GPT1001 and GPT1002. The GPT1001 model includes Equations (2), (3), (10), (11), (13) and (16). In this case, we used the optical density data provided by the authors to calculate cell density and specific growth rate using the same strategy described before. A dynamic for extracellular glucose concentration was also introduced to the model, according to Table S2. The model was adjusted to simulate a 300 mL culture. The resulting curve describes extracellular tryptophan concentration close to 1.3 g/L obtained by Gu et al. [40] after a 36 h culture (Figure 3). For GPT1002, the model includes the same equations, except for Equation (16), which is replaced by Equation (17). We calculated the specific growth rate and cell concentration from the reported data (Table S3) and we inferred a dynamic for glucose concentration (not shown). The resulting curve fits the first stages of the tryptophan dynamic (Figure 4); however, the last two points differ significantly from that experimentally defined. In the discussion section, we speculate about the possible causes of this discrepancy.

4. Discussion

In this work, we adapted a working dynamic model that includes the primary molecular components involved in bacteria’s last part of the tryptophan biosynthetic pathway. Our approach was from a biotechnological point of view; thus, we reflect in the model mRNA, attenuator, transport, regulators, and enzymatic genes whose activity is known to control tryptophan production. The model introduces the genetic and regulatory changes experimentally described in developing two E. coli strains (GPT1001 and GPT1002) [40]. Although the authors performed additional genetic modifications in those strains, such as sugar transport and pathway precursors, these were not directly included in our model since their effect could be more general in bacterial physiology. Another critical declaration is that we did not intend to reproduce what happens in a controlled bioreactor. Although it increases metabolite production, the way the production rises is possibly at the physiological level or due to the extensive residency time of bacteria in bioreactors, giving place to adaptive modifications in the population. For instance, in the referred study [40], when culturing in a 5 L bioreactor, a 20 h reaction occurred before bacteria grew and produced tryptophan. Hence, we focus on reproducing that which was obtained in that study in shake flasks. An interesting approach that may provide insight into physiology’s contribution is coupling the effects of oscillations in the central metabolism and tryptophan production [51,52].

Even with our effort, we had a discrepancy in the last two points of the culture, particularly on the GPT1002 strain (Figure 4), where the latest points fell and did not recuperate as in GPT1001 (Figure 3). We speculate that the difference between the model and the experimental observation seen in Figure 4 could be due to several reasons. Since this is a simplified model, many things may be happening in the overall E. coli metabolism that may cause a decrease in tryptophan production. A possible cause is that, according to previous results [53], intracellular RNAP and ribosomal concentrations decrease at low specific growth rates, as could happen in entering the stationary phase [54]. The last two points of the simulation in Figure 4 are also the slowest growing points of the bacterial population. Hence, it is possible to think that such diminution might be related to lower transcription and translation rates of the genes in the tryptophan operon. That condition is especially true for a strain with an increased metabolic burden. Thus, decreasing the value of the average number of ribosomes translating each mRNA (k_E), from 30 min⁻¹ to 15 min⁻¹ starting in the minute 1440, results in extracellular tryptophan dynamics that fit closer to those experimentally observed, i.e., 1.7 g L⁻¹ (Figure 5A).

Furthermore, as mentioned earlier, previous tryptophan operon models assume that tryptophan production is proportional to the concentration of free TrpE since it is considered the operon bottleneck. However, it was necessary to take a different perspective in the case of mutant GPT1002. As mentioned, Gu et al. [40] reported an overall two-fold increase in the expression of the tryptophan operon transcriptional unit. However, such an increase was not equal for all the genes in the transcriptional unit since trpE displayed a nine-fold increase, while the expression of trpD showed a two-fold increase. The differentiated expression of genes in operons has been previously reported [55]. Thus, if we had stuck to the considerations of the original model, we should have added a nine-fold increase in the expression of TrpE for mutant GPT1002. However, we reasoned that since an enzymatic aggregate carries out the first reaction of the tryptophan pathway (TrpED) and TrpE cannot catalyze the path’s first reaction in the absence of TrpD [35], the effect of the increased expression of trpE would be minimal if it is not accompanied by an increase of the same magnitude of trpD. Since such an increase did not occur, the expression of trpD would be the new bottleneck. That is why a two-fold increase in operon expression was applied in Equation (17). It appears that it was a correct deduction since a nine-fold increase in the operon expression used in Equation (17) results in tryptophan dynamics that fall far from the experimental results (Figure 5B).

Other phenomena also may not be considered here, such as the diminution of mRNA concentration during the different growth phases. However, summing the possible causes that may produce discrepancies may not be the best approach when working with a deterministic model, neglecting the stochastic nature of these interactions. So, further analysis will be needed to determine the precise causes of this discrepancy. In cultures near or on the stationary phase, the apparition of alternative phenotypes, such as dormant or persistent cells, may contribute to a less productive population [56,57].

5. Conclusions

This work shows the first approach to achieving a dynamic working model for biotechnological tryptophan production in E. coli. This model includes the contributions of the main genetic modifications and speculates the performance of the bulk of transport reactions to fit the quantity of this metabolite as measured in the culture media. There are many opportunity areas to improve the model; we think one main grey area to experimentally discern is the specific contributions of the different transport systems to obtain tryptophan in the extracellular concentrations as experimentally measured. Other unspecific or unconsidered transport systems may participate in the secretion of tryptophan when produced at high cellular rates. There is also consideration of how to include physiological changes and apparitions of additional genetic modifications in long-time cultures. Controlled conditions and feedback bioreactors can maintain optimal conditions that contribute to optimal physiology and sustain highly productive populations. We hope this work can contribute to more precise working biotechnological models to produce valuable chemicals.