An Asymptotic Analysis of the Malonyl-CoA Route to 3-Hydroxypropionic Acid in Genetically Engineered Microbes.

There has been recent interest in creating an efficient microbial production route for 3-hydroxypropionic acid, an important platform chemical. We develop and solve a mathematical model for the time-dependent metabolite concentrations in the malonyl-CoA pathway for 3-hydroxypropionic acid production in microbes, using a combination of numerical and asymptotic methods. This allows us to identify the most important targets for enzyme regulation therein under conditions of plentiful and sparse pyruvate, and to quantify their relative importance. In our model, we account for sinks of acetyl-CoA and malonyl-CoA to, for example, the citric acid cycle and fatty acid biosynthesis, respectively. Notably, in the plentiful pyruvate case we determine that there is a bifurcation in the asymptotic structure of the system, the crossing of which corresponds to a significant increase in 3-hydroxypropionic acid production. Moreover, we deduce that the most significant increases to 3-hydroxypropionic acid production can be obtained by up-regulating two specific enzymes in tandem, as the inherent nonlinearity of the system means that a solo up-regulation of either does not result in large increases in production. The types of issue arising here are prevalent in synthetic biology applications, and it is hoped that the system considered provides an instructive exemplar for broader applications.

Introduction

3-Hydroxypropionic acid (3HP) is a platform chemical which can be converted into several valuable chemicals, for example acrylic acid and acrylamide (Werpy et al. 2004). There has been significant recent interest in introducing metabolic pathways to microorganisms in order to produce 3HP at industrially viable levels (Borodina et al. 2015; Chen et al. 2014, 2017; Kim et al. 2014; Son et al. 2017; Matsakas et al. 2018). The manipulation of metabolic pathways in microorganisms is an exciting avenue of research arising in synthetic biology. As enzyme production within an organism is genetically controlled, introducing specific genes allows an organism to produce chemicals that it could not previously. One of the most successful examples of a chemical synthesized by this procedure is amorpha-4,11-diene (a precursor to artemisinin, an antimalarial) from Escherichia coli (Martin et al. 2003). It can be expensive and time-consuming to introduce a new chemical pathway to a microorganism and optimize the production of a given metabolite. Moreover, adding and blocking pathways can have unintended consequences for the metabolism of the microorganism, and the experimental parameter space is extensive. Mathematical modelling allows systematic progress to be made in understanding a pathway and can significantly reduce the experimental parameter space that needs to be searched.

In Kumar et al. (2013), three thermodynamically feasible pathways from pyruvate to 3HP are suggested. In previous work, we investigated 3HP production via the -alanine route using mathematical modelling (Dalwadi et al. 2017). Here, we are interested in mathematically modelling 3HP production via the malonyl coenzyme A (malonyl-CoA) route, with the aim of determining the most appropriate enzyme targets for regulation. We consider the reactions where (1f) and (1g) represent the loss of acetyl-CoA and malonyl-CoA to other metabolic pathways, such as the citric acid cycle or fatty acid biosynthesis. We do not track any metabolites once they reach the sink. We show a schematic representation of this pathway in Fig. 1.

Fig. 1

A schematic network diagram for the metabolic pathway we consider in this paper, from pyruvate to 3HP via malonyl-CoA. The arrows denote the direction of the reactions and are labelled by their respective catalytic constants (also referred to as turnover numbers). We only track the metabolites included in this figure and, specifically, not any involved in the acetyl-CoA or malonyl-CoA sinks. The subscript for a given enzyme (for ) corresponds to the subscript of the related maximal reaction rate. The exception here is for the reversible reaction between malonic semialdehyde and 3HP, where both reactions use the enzyme . Hence, corresponds to the pyruvate dehydrogenase complex (EC 1.2.4.1, EC 2.3.1.12, and EC 1.8.1.4), corresponds to acetyl-CoA carboxylase (EC 6.4.1.2), corresponds to malonyl-CoA reductase (EC 1.2.1.75), corresponds to 3-hydroxypropionate dehydrogenase (EC 1.1.1.298), and corresponds to malonate semialdehyde dehydrogenase (acetylating) (EC 1.2.1.18)

Our general aim is to determine how the system behaves as a function of its parameters, with our main goal being to understand how to maximize 3HP production while minimizing the levels of malonic semialdehyde, a toxic intermediate, where possible. To derive and solve our mathematical model, we make several modelling assumptions. We consider a system that is well mixed and thus spatially independent. This means that we are able to formulate a mathematical system in terms of ordinary differential equations in time, rather than partial differential equations in time and space. These differential equations require initial conditions, and we consider the case where pyruvate is instantaneously introduced to a system containing all of the relevant enzymes, but none of the intermediate metabolites. This assumption facilitates a mathematical analysis by reducing the number of unknown parameters in the system. Moreover, we follow the approach of Dalwadi et al. (2017) and Dalwadi et al. (2018) and investigate two simplified cases of pyruvate replenishment, which model the extremes of the actual time-dependent pyruvate replenishment. The first of these is continuous pyruvate replenishment, where the pyruvate is held at a constant concentration, which could represent a continuous culture, and the second is no pyruvate replenishment which could represent a batch culture. Understanding these extreme cases allows us to determine the key targets for enzyme regulation. Finally, we assume that the formation rate of enzyme complex production is much quicker than the rate of substrate consumption, and thus, the reaction rates are governed by Michaelis–Menten-type laws, the specific form of which we obtain from the literature.

To analyse the nonlinear governing equations that we derive, we use a combination of numerical and asymptotic methods. The latter enhances our physical insight into the underlying system and allows us to derive closed-form expressions for how the metabolite concentrations vary as functions of the experimental parameters. Asymptotic techniques (see, for example, Holmes 2012; Kevorkian and Cole 2013) allow us to largely bypass the issue of uncertainty in the parameters, as deriving analytic approximations of the dynamic metabolite concentrations only requires an understanding of the relative order of magnitude of each experimental parameter. Moreover, the nonlinear nature of the system means that a broad understanding of the system behaviour cannot be obtained by simply varying one parameter at a time and collecting system outputs, so asymptotic solutions allow for a quicker and more comprehensive understanding of the system. We note that there is a bifurcation in the asymptotic structure of the system for the case of continuous pyruvate replenishment, and we investigate this in terms of our goal of maximizing 3HP production.

In Sect. 2, we introduce a mathematical model to describe the reaction kinetics. We solve this system numerically and asymptotically, in Sect. 3 for the case with a continuous replenishment of pyruvate, and in Sect. 4 for the case with no replenishment of pyruvate. Finally, in Sect. 5 we discuss our results.

Model Description

The method we use to set up our governing equations is similar to that of Dalwadi et al. (2018), but for a different metabolic network. Thus, to describe the dynamics of the network from pyruvate to 3HP, shown in Fig. 1, we obtain corresponding Michaelis–Menten-type reaction velocities from the literature. The resulting dimensional governing system is where each variable is defined in Table 1. We will consider two extreme cases of system (2), corresponding to continuous and no replenishment of pyruvate, respectively, noting that the reality will lie somewhere between these two extremes. In (2), most of the reaction velocity terms are in standard Michaelis–Menten form, and we obtain the kinetic parameters for these reaction velocities from the corresponding references in Table 2. A lower-case k represents a catalytic constant, and an upper-case K represents a Michaelis (or Michaelis-type) constant. We include a modified Michaelis–Menten term for reaction (1a), which encompasses the uncompetitive inhibitory effect that acetyl-CoA has with pyruvate (Kresze and Ronft 1981a). Therefore, we use the standard form for uncompetitive inhibition for this reaction velocity (Yung-Chi and Prusoff 1973). Reactions (1f) and (1g) correspond to the aggregated loss of acetyl-CoA and malonyl-CoA, respectively, to other metabolic pathways present in the microorganism, such as the citric acid cycle or fatty acid biosynthesis. We model each of these aggregated losses as a sink with first-order reaction kinetics, using the parameters and to represent the strength of each sink for acetyl-CoA and malonyl-CoA, respectively. To get around the issue that these parameters are unknown, we will later consider a distinguished asymptotic limit (in a sense to be made formal later) where both of these sink reactions balance the other reactions. The parameters , where , denote the initial enzyme concentrations for the reactions they control, noting that the reversible reaction with maximal reaction rates and are both controlled by .

Table 1

Dimensional and dimensionless variable definitions

Original variable	Description	Dimensionless variable
[S_1]	Pyruvate	[S_1]=S_0S_1
[S_2]	Acetyl-CoA	[S_2]=S_0S_2
[S_3]	Malonyl-CoA	[S_3]=S_0S_3
[S_4]	Malonic semialdehyde	[S_4]=S_0S_4
[P]	3HP	[P]=S_0P
τ	Time	τ=(S_0/k_1E_1)t

Table 2

Parameter values

Dimensional	Organism	Range	Dimensionless
k_1=15s^-1	Lactococcus lactis (Snoep et al. 1992)	4–30s^-1 (Kresze and Ronft 1981a, b; Snoep et al. 1992)
k_2=12s^-1	Acinetobacter baumannii (Alves et al. 2011)	12–30s^-1 (Alves et al. 2011; Hügler et al. 2003)	k¯_2=k_2E_2/ϵk_1E_1=11.4
k_3=28s^-1	Sulfolobus tokodaii (Alber et al. 2006)	3–50s^-1 (Hügler et al. 2002; Alber et al. 2006)	k¯_3=k_3E_3/k_1E_1=0.19
k_4=100s^-1	Nitrosopumilus maritimus (Otte et al. 2015)	100–120s^-1 (Otte et al. 2015; Kockelkorn and Fuchs 2009)	k¯_4=k_4E_4/k_1E_1=0.67
k_-4=2.1s^-1	Nitrosopumilus maritimus (Otte et al. 2015)	0.2–2.1s^-1 (Otte et al. 2015; Yao et al. 2010)	k¯_-4=k_-4E_4/ϵk_1E_1=2
k_5=20s^-1	Pseudomonas fluorescens (Yamada and Jakoby 1960)	9–20s^-1 (Yamada and Jakoby 1960; Anderson and Magasanik 1971)	k¯_5=k_5E_5/k_1E_1=0.13
K1M=1mM	Lactococcus lactis (Snoep et al. 1992)	0.13–1mM (Kresze and Ronft 1981a; Pronk 1996; Snoep et al. 1992)	K¯1M=K1M/S_0=0.5
K1i=0.014mM	Saccharomyces cerevisiae (Kresze and Ronft 1981a)	0.014–0.018mM (Kresze and Ronft 1981a; Pronk 1996)	ϵ=K1i/S_0=0.007
K2M=0.01mM	Acinetobacter baumannii (Alves et al. 2011)	0.01–0.06mM (Alves et al. 2011; Hügler et al. 2003)	K¯2M=K2M/ϵS_0=0.71
K3M=0.04mM	Sulfolobus tokodaii (Alber et al. 2006)	0.03–0.1mM (Hügler et al. 2002; Alber et al. 2006)	K¯3M=K3M/ϵS_0=2.86
K4M=0.11mM	Nitrosopumilus maritimus (Otte et al. 2015)	0.07–0.11mM (Otte et al. 2015; Kockelkorn and Fuchs 2009)	K¯4M=K4M/ϵS_0=7.86
K-4M=7.9mM	Nitrosopumilus maritimus (Otte et al. 2015)	7.9–17mM (Otte et al. 2015; Yao et al. 2010)	K¯-4M=K-4M/S_0=3.95
K5M=0.03mM	Pseudomonas fluorescens (Yamada and Jakoby 1960)	0.03–0.06mM (Yamada and Jakoby 1960; Anderson and Magasanik 1971)	K¯5M=K5M/ϵS_0=2.14
A[s^-1]			A¯=AS_0/k_1E_1=ϵ
B[s^-1]			B¯=BS_0/k_1E_1=1/ϵ

To give an idea of the range of these parameter values, we provide data from several organisms. As we use an asymptotic analysis, it is the relative magnitude of these values, rather than their exact values, which are important. As discussed in the main text, we use the value , assume that for , and we scale extreme parameter ratios with the small dimensionless parameter , defined as the ratio of and . As there is a distinguished asymptotic limit when and , we use these scalings in our asymptotic analysis, choosing and in our simulations

The initial levels of each metabolite present in the system are relatively unknown. Thus, we make a modelling choice to reduce the number of uncertain parameters in the system to facilitate a simplified analysis of the system. We use initial conditions corresponding to the case in which pyruvate is instantaneously introduced to a system containing all of the relevant enzymes, but none of the intermediate metabolites. That is, we use , where represents the prescribed initial or typical level of pyruvate present in the system, with none of remaining metabolites initially present: .

We now nondimensionalize the system and exploit various small dimensionless parameter ratios to introduce a small parameter into the system, thus allowing us to write the system in terms of constants and one small parameter to quantify the relative size of each term. This is useful because kinetic parameters can vary between different environments, as shown in Table 2. To side-step the issues associated with this parameter uncertainty, we seek to quantitatively understand how the system behaves as these parameters vary by interrogating the system using an asymptotic analysis (see, for example, Holmes 2012; Kevorkian and Cole 2013).

We provide our nondimensional variable scalings in Table 1. Essentially, we scale each dimensional metabolite concentration with , the initial concentration of pyruvate, and we scale time with , the characteristic time of the first reaction, which occurs between pyruvate and acetyl-CoA. To form the dimensionless parameters in the system, we first scale each rate constant with the rate constant of the first reaction, each Michaelis constant with the initial pyruvate concentration, and each enzyme concentration with the concentration of the first enzyme. Then, using the typical dimensional values in Table 2, we note that the parameter ratio is very small. Hence, we scale the other extreme parameter ratios in our system with the small dimensionless parameter , allowing us to write each additional dimensionless parameter as , where c is an parameter (in practice, we take this to be between and ), and j is an integer. This allows us to quantify the “smallness” of each parameter while allowing each parameter c to vary over approximately three orders of magnitude, yielding a system in terms of a small parameter. We take the typical concentration of pyruvate within a microorganism to be , and we assume that for . We make this assumption since we expect the levels of pyruvate dehydrogenase complex (PDC), corresponding to , to be more abundant than the other enzymes in this pathway since PDC is fundamental to cellular respiration in that it links glycolysis to the citric acid cycle. However, we emphasize that different values of can be considered if required by varying the appropriate dimensionless parameter. The resultant dimensionless parameters we use in our system are given in Table 2.

The dimensionless system is then given by where is a small parameter. The dimensionless initial conditions are: , , , , and .

We emphasize that the asymptotic sizes of and are not currently set, and that we have not yet scaled our dependent variables with . As we are not able to accurately estimate the typical sizes of and , we instead consider the distinguished asymptotic limit in which they balance the core reaction velocities over the same timescale: from pyruvate to acetyl-CoA, from acetyl-CoA to malonyl-CoA, from malonyl-CoA to malonic semialdehyde, and from malonic semialdehyde to 3HP. The reaction velocity from pyruvate to acetyl-CoA has asymptotic size . The reaction velocity from acetyl-CoA to malonyl-CoA has asymptotic size . Balancing these two reaction velocities, and noting that we start with , we may deduce that the size of these reaction velocities must be of , and hence we obtain the scaling . For these reaction velocities to balance with , we must consider the distinguished limit . Balancing these terms with the reaction velocity from malonyl-CoA to malonic semialdehyde, with asymptotic size , we obtain the scaling . Then, the balance with yields the distinguished limit . Therefore, we use the scalings and , where . This appears to be the only distinguished limit of interest.

We now solve this system for two different cases of pyruvate replenishment. The first is continuous pyruvate replenishment, where the pyruvate dynamics are not governed by (3a), and instead we impose . The second is no pyruvate replenishment, and (3a) does hold. We will deduce that the two systems are equivalent for , and only diverge when . In both cases, we are interested in determining how to maximize 3HP production. For the continuous-replenishment case, we are also interested in trying to minimize the long-time levels of malonic semialdehyde, a toxic intermediate compound, present in the cell over a long timescale. As the no-replenishment-of-pyruvate case eventually tends to no metabolites in the system, the long-time levels of malonic semialdehyde in that case are trivial. We now consider the continuous-pyruvate-replenishment case.

Continuous Replenishment of Pyruvate

We first consider the case where the pyruvate is continuously replenished, modelling the bacteria being harvested in a continuous culture. A numerical simulation of system (3b–e) with shows that the 3HP concentration is bounded above as (Fig. 2). For our goal of maximizing 3HP production, we would like to understand how to improve the 3HP production by regulating the enzymes in the system. In this section, we solve the nonlinear system (3b–e) with through an asymptotic analysis, exploiting the small parameter . Our analysis reveals that there is a bifurcation in the asymptotic structure for the large-time behaviour of 3HP with important implications for 3HP production. Our results allow us to quantify this bifurcation in terms of the system parameters, providing physical insight and suggesting which reactions should be the key targets for genetic manipulation.

Fig. 2

The dynamics of each metabolite in the system when pyruvate is continually replenished. The solid grey lines are the numerical solutions from system (3b–e) with , the dashed blue lines are the early-time asymptotic results from (6), and the dotted black lines are the late-time results, where we use the asymptotic results (11), (12), and (15) for , , and , respectively. These asymptotic results yield a single reduced ODE for , given in (9), and it is the numerical solution to this that we plot for as a dotted black line. We use the parameter values given in Table 2 (Color figure online)

Asymptotic Structure

We now discuss the asymptotic structure of the continuous-replenishment case, with reference to the metabolic network shown in Fig. 1. In this case, there are two important timescales in the system: early time, where , and late time, where . Over the early timescale, the levels of malonyl-CoA and malonic semialdehyde reach a steady and quasi-steady state, respectively. The remaining metabolite concentrations increase through this early timescale, with a change in the power law governing their growth. Over the late timescale, the levels of acetyl-CoA reach a steady state, but the behaviour of the 3HP depends on the system parameters. There is a bifurcation in the asymptotic structure for the large-time behaviour of 3HP, the levels of which are bounded above unless some critical parameter ratio is reached. We investigate this bifurcation when we consider the late-time behaviour.

Early Time:

We start our analysis over the distinguished timescale , and thus, we make the scaling . By balancing reaction velocities, we find that the appropriate scalings for the dependent variables over this timescale are , which converts system (3b–e) into The leading-order version of (4) iswhere for . We may solve (5) to obtain We see that the early-time results (6) (dashed black lines) agree very well with the numerical results (solid grey lines) when in Fig. 2. However, the metabolite concentrations stray from these asymptotic results (with the exception of ) as t becomes larger than of . This is because new terms in the system become of leading order, and so we must explore a new asymptotic region to understand the system fully. We investigate this in the next section. As such, it is helpful to state the large- limits of solutions (6), as these will allow us to match appropriately into the next timescale. In this limit, we obtain the following expressionsand hence we see that and reach a steady state over , while and continue to grow over the same timescale. In the next subsection, we consider the remaining important timescale .

Late Time:

The second and final distinguished timescale for continuous replenishment is , so we now use the scaling . This is the timescale over which there is a balance between the source from pyruvate and the sinks from acetyl-CoA and malonyl-CoA. Over this timescale, and are both of , and the scalings for the remaining metabolites are still , which follow from (7), the large-time limits of the early-time solutions. Using the scalings mentioned above, we obtain the late-time version of (3b–e): From (8), we obtain the leading-order differential–algebraic system, given by The “initial” conditions of the system are obtained by matching with the large-time results of the early-time system (7), to yieldFrom (9c), we can immediately deduce that is constant over this timescale. That is, over the late timescale, takes the constant valueFrom (9)–(11), we can write in terms of as follows:Substituting (12) into (9), we obtain the following closed ODE for :where the critical parameter is defined asand we shall justify shortly the critical nature of . The differential equation (13) is solved implicitly byThe asymptotic late-time results (11), (12), and (15) (dotted black lines) agree well with the numerical results (solid grey lines) when in Fig. 2. The remaining nonlinear differential equation (9a) decouples from the rest of the system, and solving this equation is not required to understand the dynamics of the rest of the system. However, we note that the remaining equation will be important when there is no replenishment of pyruvate, which we consider in the next section. To understand the dynamics of , (9a) must be solved numerically, and the solution to this agrees well with the full numerical solution of in Fig. 2.

Given that the argument of the logarithm in (15) will vanish for some positive value of if and only if , we may deduce that the large-time behaviour of has a critical dependence on the sign of . Thus, is the bifurcation point for the parameter , in terms of changing the asymptotic structure of the solution. Specifically, when , corresponding to a weak reaction from malonic semialdehyde to 3HP, we obtain the bounded large-time behaviour however when , corresponding to a strong reaction from malonic semialdehyde to 3HP, we obtain the unbounded large-time behaviour We show this bifurcation behaviour in Fig. 3a, including a comparison between the numerical and asymptotic results. We see that our numerical results show the bifurcation behaviour predicted by our asymptotic results. Moreover, our asymptotic results are good predictions of both the location of the bifurcation and the long-time behaviour.

Fig. 3

The large-time limits of a 3HP production and b malonic semialdehyde varying over the critical parameter . The grey lines denote the asymptotic predictions [from (16) and (18)], and the black crosses denote numerical results (solving (3b–e) with ). The subcritical and supercritical regimes are to the left and right, respectively, of the dashed lines in each subfigure. To obtain the large-time numerical limits, we run the simulations until ; we approximate and . Apart from , which we vary in this figure, we use the parameter values given in Table 2, and these correspond to an asymptotic value of the critical parameter

It is of mathematical interest to briefly note that when , the solution to (13) isHence, as passes through the bifurcation point, the long-time behaviour of is to increase with the square root of time.

From (12), we see that will tend to a constant in the far-field; for , we have the subcritical regime and for , we have the supercritical regime We show this bifurcating behaviour in Fig. 3b, including a comparison between the numerical and asymptotic results. As before, we see that our numerical results display the bifurcation behaviour predicted by our asymptotic results, and that the latter are a good predictor of the bifurcation behaviour. Recall that our goal for the continuous-pyruvate-replenishment case is to maximize 3HP production while minimizing the large-time levels of malonic semialdehyde, the latter being a toxic intermediate compound. Therefore, as large-time 3HP production is zero and malonic semialdehyde levels are higher in the subcritical regime, Fig. 3b tells us immediately that the subcritical regime is a bad regime for industrially viable 3HP production, and that the parameter regime is much better for our goal. From Table 2, we note that the critical parameter , but . Therefore, the simulations in Fig. 2 occur in the subcritical regime, as can be deduced by noting that the 3HP concentration is bounded above. Our model predicts that one can achieve large gains in 3HP production if one is able to move into the supercritical regime.

Physical Implications

To frame our discussion around parameters that can be varied experimentally (essentially the levels of each enzyme, regulated by the over- or under-expression of the genes that control their production), we rewrite our results in terms of dimensional quantities. In dimensional terms, the bifurcation occurs when , whereFrom (19), we see that an up-regulation of and a down-regulation of have the most significant effect on lowering , the critical bifurcation parameter, and thus hopefully moving into the supercritical regime for significantly increased 3HP production. Additionally, an up-regulation of will also have an effect in lowering , but with diminishing returns. Notably, our model suggests that regulation of neither nor has a significant effect on this critical bifurcation parameter (see below for a potential important role for , however).

The supercritical regime occurs when , and this corresponds to a nonzero large-time production of 3HP, which is the target if we are to attain the goal of industrially viable 3HP production. In this case, the dimensional large-time production of 3HP isand the dimensional large-time levels of malonic semialdehyde in the system areAs our goal for the continuous-replenishment case is to maximize 3HP production while minimizing the large-time levels of malonic semialdehyde, we consider the ratio of the large-time levels of to , to obtainWhile (22) appears to have a convoluted form, it provides a quantitative value that we wish to maximize. In particular, we are able to see that the most significant gains can be made by up-regulating and in tandem. Our model suggests that an increase in only one of these enzymes will increase ratio (22), but with diminishing returns (Fig. 4). However, increasing both enzymes will result in significant increases in the 3HP production to malonic semialdehyde ratio (Fig. 4). Additionally, we note that 3HP production can also be increased by up-regulating and down-regulating , though this will have diminishing returns. We emphasize that maximizing (22) should not be considered a definitive metric, since the toxic effect of malonic semialdehyde () will have its own nonlinear effects for which to account. Other relevant metrics could involve keeping the maximum value of malonic semialdehyde below a critical value or keeping the cumulative presence of malonic semialdehyde low. Our multiscale approach allows for the efficient analysis of any such metric.

Fig. 4

The effect of regulating enzymes in the pathway shown in Fig. 1 on the key ratio (22), obtained by solving (3b–e) with . We use the parameter values in Table 2, but with increased enzyme concentrations as specified on the x-axis. While up-regulating yields improvement, this has diminishing returns. However, even though up-regulating by itself does not appear to have a significant effect on the ratio, up-regulating and in tandem does have a significant effect (Color figure online)

No Replenishment of Pyruvate

We now consider the case where the pyruvate is never replenished. Due to the connectivity of the metabolic network and the sinks in the system, this problem results in a long-term decay of each metabolite. As the only metabolite initially present in the system is pyruvate, the remaining metabolites will each reach a maximum level at some (generally different) point in time. Our main goal in this subsection is to obtain an asymptotic expression for the maximum level of 3HP present in the system, and the time at which this occurs. In a batch culture system, where nothing is added or removed to the system until the run is stopped for harvesting, this time provides an approximation for when to stop the run and harvest the product. We will summarize the system that governs this optimal time and product yield, but we relegate the details to “Appendix A” for brevity.

We now discuss the asymptotic structure of the no-replenishment case, with reference to the metabolic network shown in Fig. 1. The no-replenishment-of-pyruvate case is equivalent to the continuous-replenishment case at leading order until we reach the timescale (as shown in Fig. 5 for , , , and P), at which point the depletion of pyruvate occurs as a leading-order effect.

Fig. 5

A comparison between the metabolic dynamics in the continuous- and no-replenishment-of-pyruvate cases. The solid grey lines are the numerical solutions from the continuous-replenishment-of-pyruvate system (3b–e) with , and the dashed blue lines are the numerical solutions from the no-replenishment-of-pyruvate system (3). We see that the two cases are essentially equivalent until around , after which the metabolite levels drop in a sharp manner for , , and , and at a slower rate for P. The maximal level of 3HP occurs around the time of the sharp drop. We use the parameter values given in Table 2 for these figures (Color figure online)

Over this longer timescale, the rest of the system remains unchanged until the levels of pyruvate are of , which we define to occur at . The pyruvate is then depleted to a negligible level at a relatively quicker rate, over a timescale of , as can be seen in Fig. 6, while the remaining metabolites are unchanged at leading order. After this rapid depletion of pyruvate, the reaction from pyruvate to acetyl-CoA is negligible and so the levels of acetyl-CoA start to deplete, again over the timescale of . When the levels of acetyl-CoA become of , which we define to occur at , the Michaelis–Menten reaction from acetyl-CoA to malonyl-CoA becomes unsaturated, and the lower levels of acetyl-CoA are felt through the rest of the system. After this point, all the remaining metabolites start to deplete over a timescale of , and so marks the time at which the level of 3HP in the system is maximal. In the next section, we present the reduced system required to be solved to determine .

Fig. 6

The dynamics of pyruvate () when it is never replenished, obtained from solving (3). We see that it becomes exponentially small shortly after . We use the parameter values given in Table 2

Reduced System for

The time at which the level of 3HP in the system is maximal is given by , where . To determine at leading order, it suffices to solve the reduced system and we provide a derivation of this reduced system in “Appendix A”. In (23a), the nonlinear function f is given byAdditionally, since the remaining metabolite concentrations are equivalent to those in the continuous-replenishment case until acetyl-CoA depletion occurs, the concentration of malonic semialdehyde is given by the implicit representation (12), (14)–(15). As shown in “Appendix A”, the asymptotic leading-order time at which pyruvate depletion occurs is , where satisfies the implicit equationThe time at which there is a maximal amount of 3HP in the system is when acetyl-CoA depletion occurs, since this is when the remaining metabolites in the system start to feel the lack of the previous metabolites. Asymptotically, this occurs when or, equivalently, for satisfying the implicit equation which is obtained by integrating (23a) in time and using (23c). It is a simple task to solve this system numerically for any specified parameter values. Then, the maximal amount of 3HP in the system can be obtained by substituting into (15) to yield the implicit equationComparing this reduced system to the full system in Fig. 7, we see that the maximal 3HP produced is very similar between the asymptotic and numerical results, and that there is an difference in the predicted time at which this maximal 3HP is produced between the two, as predicted by the error in our asymptotic results. We show this excellent agreement for a range of relative concentrations of (Fig. 8), the enzyme which recycles malonic semialdehyde back into acetyl-CoA. Moreover, we see that more 3HP is produced when is down-regulated.

Fig. 7

A comparison of numerical solutions for the full system in (3) and the reduced system presented in Sect. 4.2. Our solutions for the reduced system are valid for . Here, is the asymptotic result for the point at which the 3HP in the system is maximal, and in (24) gives the asymptotic result for the maximal level of 3HP. We mark the maximal value of 3HP in the system with a cross in the appropriate colour for the numerical and asymptotic solutions, respectively. We use the parameter values given in Table 2 (Color figure online)

Fig. 8

The maximum amount of 3HP in the system as the concentration of is varied, for the full numerical system (3) and the reduced system presented in Sect. 4.2. We use parameter values from Table 2, with modified as specified on the x-axis. The maximum discrepancy between numerical and asymptotic solutions for the maximum value of 3HP is around 4%, which occurs for smaller relative concentrations of

An alternative implicit equation for is given by (45), which uses the analytic solution for when , stated in (44).

Asymptotic Expression for when

Our model predicts that 3HP production can be improved by under-expressing the enzyme that controls the reaction from malonic semialdehyde to acetyl-CoA, which has maximal reaction rate (Fig. 8). If this reaction is completely removed from the metabolic network, we are able to provide an analytic expression for and . This corresponds to taking the limit , equivalent to . In this limit, there is no further production or depletion of 3HP when , and the total amount of 3HP produced is given by .

As , there are several scalings that reduce the parameter dependence of the system. We scale and , to turn (23a) into where , , the nonlinear function is defined asand has the form where , which will be useful later for notational purposes.

When , the system (25a) is solved implicitly byWe are able to obtain an explicit expression for by re-writing (25c) aswhere we use (25a) for . Re-arranging (27), we obtain We can then obtain an explicit expression for by substituting (27)–(28a) into (26), to yield

When , system (25a) is solved explicitly by and as satisfies , we can state in terms of as follows:where each term in (30) can be determined explicitly from (28).

As we provide the general implicit result for maximal 3HP in (24), we can derive a reduced version for the limit we consider in this subsection by taking the same limit, , noting that also scales with . Using (30), the reduced result for the time at which the levels of 3HP are maximal, the total amount of 3HP produced is given explicitly by To visualize (31), we provide a plot of in Fig. 9, but as a function of and , as these are proxies for the strength of the sink and source in the system, respectively. As F is bounded above by 1 and this limit is attained as , decreasing (and hence ) to zero will have a more significant effect than any finite increase in (and hence ), as long as the organism is able to survive with low values of .

Fig. 9

The asymptotic result for , the total 3HP produced with no replenishment of pyruvate in the limit of , given in (31). We use axes and , which are proxies for the strength of the sink and source in the system, respectively (Color figure online)

Physical Implications

In a similar manner to the analysis in Sect. 3.4, we now discuss the physical implications of our results from this subsection and we reintroduce dimensional variables to our key results.

The first key result is the time at which the 3HP levels are maximal, , whereThe effective production rate of 3HP iswhich can be increased by increasing or , though the latter has diminishing returns. The total amount of 3HP produced is given bywhere F is defined in (31). Definingas proxies for and , we see that up-regulating results in higher levels of 3HP, whereas up-regulating results in lower levels of 3HP, with significant diminishing returns as increases (Fig. 10). The negative effect of up-regulating may appear counter-intuitive, since is required for a connected pathway. However, it can be explained by noting that the effect of increasing is to increase the amount of acetyl-CoA in the system, which has the result of greatly increasing the amount of acetyl-CoA that is taken up by the sink, but only slightly increasing the amount converted into malonyl-CoA. This is because we have assumed that we are in a regime where the sink has first-order kinetics, whereas the reaction from acetyl-CoA to malonyl-CoA follows Michaelis–Menten-type kinetics. This effect may be reduced if we were in a regime where the acetyl-CoA sink was saturating. Hence, our model suggests that a slow conversion of pyruvate into acetyl-CoA is better for 3HP production if it is not possible to over-regulate . We note that this down-regulation is only important for the maximum amount of 3HP produced when the time taken to produce this 3HP is taken into account, as there is no dependence on in the effective production rate (33). We emphasize that is a singular limit in the system, which appears because we are in the no-replenishment-of-pyruvate regime. As , the early-time kinetics of the problem would vary, and we would not have the same late-time problem that we consider here. We do not analyse this limit further, since it is not of particular physical interest.

Fig. 10

A plot showing , the asymptotic result for the total 3HP produced with no replenishment of pyruvate in the limit of , given in (31). In contrast to Fig. 9, we present F in terms of a variation in the parameters and (as shown in (35)), corresponding to regulation of enzymes (Color figure online)

Discussion

We develop and solve a mathematical model for the reaction kinetics of the 3-hydroxypropionic acid (3HP) pathway from pyruvate via malonyl-CoA. We consider two main cases, continuous and no replenishment of pyruvate, to model the metabolite dynamics in the two extreme cases of plentiful and scarce pyruvate. Our main goal is to understand how to maximize 3HP production, and in the case of a continuous pyruvate replenishment to additionally minimize the maximum levels of the toxic intermediary malonic semialdehyde in the system, where possible. We summarize our enzyme regulation recommendations for each case in Table 3.

Table 3

A summary of the enzyme regulations we recommend based on the results of our model

Enzyme [EC number(s)]	Continuous replenishment	No replenishment
Pyruvate dehydrogenase complex [1.2.4.1, 2.3.1.12, 1.8.1.4]	No effect	Down-regulate
Acetyl-CoA carboxylase [6.4.1.2]	Up-regulate	Up-regulate
Malonyl-CoA reductase [1.2.1.75]	Up-regulate	Up-regulate
3-Hydroxypropionate dehydrogenase [1.1.1.298]	Up-regulate	Up-regulate
Malonate semialdehyde dehydrogenase (acetylating) [1.2.1.18]	Down-regulate	Down-regulate

We give results for the two cases of continuous and no replenishment of pyruvate

These are the extremes of the actual time-dependent pyruvate replenishment in the cell and can model continuous and batch culture, respectively

We note that while our regulation recommendations suffer from diminishing returns if implemented on a single enzyme at a time, we can overcome this issue in the continuous-replenishment-of-pyruvate case if we up-regulate acetyl-CoA carboxylase (EC 6.4.1.2) in tandem with 3-hydroxypropionate dehydrogenase (EC 1.1.1.298). Even though a solo up-regulation of either of these enzymes has diminishing returns on improving the ratio of 3HP production rate to maximum malonic semialdehyde, our analysis demonstrates that the up-regulation of both enzymes at the same time has the nonlinear effect of strongly increasing this ratio with no diminishing returns, as shown in Fig. 4.

We also emphasize that the regulation of pyruvate dehydrogenase complex (EC 1.2.4.1, EC 2.3.1.12, and EC 1.8.1.4) is only important when the pyruvate is scarce. In this case, the total 3HP produced can be increased, perhaps counterintuitively, by a slight down-regulation of pyruvate dehydrogenase complex. This is because there are competing destinations for acetyl-CoA—either a sink out of the system or onwards to malonyl-CoA—and the route out of the system becomes favoured as the amount of acetyl-CoA in the system increases. Hence, a slower production of acetyl-CoA from pyruvate is preferable to channel more acetyl-CoA towards malonyl-CoA rather than to the sink, which means a down-regulation of pyruvate dehydrogenase complex. We note that this effect disappears when we account for the time taken to produce the maximal 3HP; the effective production rate of 3HP in this case has no dependence on pyruvate dehydrogenase complex. Moreover, we emphasize that this down-regulation effect is a singular limit of the system. Hence, one cannot obtain maximal 3HP by completely removing pyruvate dehydrogenase complex from the system—in fact the amount of 3HP produced would be zero if this enzyme were completely removed from the system.

To deal with an inherent uncertainty in the parameter values, we maximize our analytic progress by using an asymptotic analysis to solve the system, exploiting the significant difference in typical parameter ratios. This allows us to systematically understand the effect of up- and down-regulating different enzymes without resorting to an expensive parameter sweep. In the continuous-replenishment-of-pyruvate case, our asymptotic analysis reveals a bifurcation in the asymptotic structure of the problem; we quantify the effect of this bifurcation for 3HP production and derive an analytic expression for the critical surface in parameter space at which this occurs. Physically, this bifurcation corresponds to a sudden change in the large-time 3HP production across the critical parameter.

For reactions between the metabolites we tracked, we use modified Michaelis–Menten reaction velocities obtained from the literature, including the effect of inhibition in the reaction from pyruvate to acetyl-CoA from both of these metabolites. We add the effect of acetyl-CoA and malonyl-CoA loss to other metabolic pathways by including an aggregate sink from these metabolites, governed using first-order reaction kinetics. The initial conditions we used corresponded to an instantaneous introduction of pyruvate to a well-mixed solution of enzymes containing no other metabolites. We choose these conditions for mathematical convenience, as the real cells are likely to contain an unknown level of the other metabolites in the system when production is initiated. We do not expect this to be an issue for the case with a continuous replenishment of pyruvate, in that we expect a small initial amount of intermediate metabolites to change the early-time problem only. However, this is not true for the case with no replenishment of pyruvate, as the depletion dynamics will be sensitive to the initial conditions. Obtaining these measurements for the initial conditions would be difficult in vivo, though can be achieved by rapidly quenching the cellular metabolism before using mass spectrometry to examine the metabolites (Bennett et al. 2009). In addition, there may be cases where the assumption of a well-mixed system does not hold, in particular over larger lengthscales, e.g. a colony of bacteria. In such cases, one could determine effective reaction rates by coupling the nonlinear reaction dynamics we consider here to spatial transport processes over longer lengthscales. This procedure could be carried out in a computationally efficient manner if one exploited the extreme ratios of the different lengthscales in the system, as performed in Dalwadi et al. (2018) for diffusive transport and linear reactions.

In this paper, we have accounted for the toxicity of malonic semialdehyde by constructing appropriate metrics to evaluate our results. That is, the toxic effect is not directly included in our model assumptions, only analysed post hoc. Since this toxicity has a negative effect on biomass production, it could be included as a negative feedback in (1) if appropriately quantified.

Finally, we note how this works highlights how mathematical modelling and asymptotic techniques can be used to understand a biological system, and to address the key questions facing experimentalists. In this case, to identify a combination of enzyme regulation with a greater theoretical output than could be obtained by measuring outputs from regulating just one enzyme at a time. This key insight into the system behaviour was possible due to the asymptotic analysis allowing us to overcome uncertainty in parameter values. We hope that these methods can be used to understand other biological systems, and to reduce the time taken to explore their experimental parameter space.