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. 2014Chen et al. , 2017Son 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 Acetyl-CoA Malonyl-CoA Malonic semialdehyde Malonic semialdehyde 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. 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.

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 E i (for i = 1, . . . , 5) 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 E 4 . Hence, E 1 corresponds to the pyruvate dehydrogenase complex (EC 1.2.4.1, EC 2.3.1.12, and EC 1.8.1.4), E 2 corresponds to acetyl-CoA carboxylase (EC 6.4.1.2), E 3 corresponds to malonyl-CoA reductase (EC 1.2.1.75), E 4 corresponds to 3-hydroxypropionate dehydrogenase (EC 1.1.1.298), and E 5 corresponds to malonate semialdehyde dehydrogenase (acetylating) (EC 1.2.1.18) 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  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 , 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 A and B 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 E i , where i = 1, . . . , 5, denote the initial enzyme concentrations for the reactions they control, noting that the reversible reaction with maximal reaction rates k 4 and k −4 are both controlled by E 4 . 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 [S 1 ](0) = S 0 , where S 0 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 O(1) 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 S 0 , the initial concentration of pyruvate, and we scale time with S 0 /(k 1 E 1 ), 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  (Snoep et al. 1992) 0.13-1 mM (Kresze and Ronft 1981a;Pronk 1996;Snoep et al. 1992)K 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 S 0 = 2 mM, assume that E  Table 2, we note that the parameter ratio := K i 1 /S 0 = 7 × 10 −3 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 c j , where c is an O(1) parameter (in practice, we take this to be between 1/2 and −1/2 ), 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 S 0 = 2 mM, and we assume that E i /E 1 = 0.1 for i ∈ {2, . . . 5}. We make this assumption since we expect the levels of pyruvate dehydrogenase complex (PDC), corresponding to E 1 , 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 E i 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: We emphasize that the asymptotic sizes ofĀ andB 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Ā andB, 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 O(min(1, S 1 , /S 2 )). The reaction velocity from acetyl-CoA to malonyl-CoA has asymptotic size O(min( , S 2 )). Balancing these two reaction velocities, and noting that we start with S 1 = O(1), we may deduce that the size of these reaction velocities must be of O( ), and hence we obtain the scaling S 2 = O(1). For these reaction velocities to balance withĀS 2 , we must consider the distinguished limitĀ = O( ). Balancing these terms with the reaction velocity from malonyl-CoA to malonic semialdehyde, with asymptotic size O(min(1, S 3 / )), we obtain the scaling S 3 = O( 2 ). Then, the balance withB S 3 yields the distinguished limitB = O(1/ ). Therefore, we use the scalingsĀ = Â andB =B/ , whereÂ,B = O(1). 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 S 1 (t) ≡ 1. The second is no pyruvate replenishment, and (3a) does hold. We will deduce that the two systems are equivalent for t = O(1), and only diverge when t = O(1/ ). 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 noreplenishment-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 S 1 (t) ≡ 1 shows that the 3HP concentration is bounded above as t → ∞ (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 S 1 (t) ≡ 1 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.

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 t = O( ), and late time, where t = O(1/ ). 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. 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 S 1 ≡ 1, 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 S 3 , S 4 , and P, respectively. These asymptotic results yield a single reduced ODE for S 2 , given in (9), and it is the numerical solution to this that we plot for S 2 as a dotted black line. We use the parameter values given in Table 2 (Color figure online)

Early Time: t = O( )
We start our analysis over the distinguished timescale t = O( ), and thus, we make the scaling t = t. By balancing reaction velocities, we find that the appropriate scalings for the dependent variables over this timescale are The leading-order version of (4) is where 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 t = O( ) = O(10 −2 ) in Fig. 2. However, the metabolite concentrations stray from these asymptotic results (with the exception of S 3 ) as t becomes larger than of O( ). 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-t limits of solutions (6), as these will allow us to match appropriately into the next timescale. In this limit, we obtain the following expressionŝ and hence we see thatŜ 3 andŜ 4 reach a steady state over t = O( ), whileŜ 2 and P continue to grow over the same timescale. In the next subsection, we consider the remaining important timescale t = O(1/ ).

Late Time: t = O(1/ )
The second and final distinguished timescale for continuous replenishment is t = O(1/ ), so we now use the scaling t = T / . 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, S 2 and P are both of O(1), and the scalings for the remaining metabolites are still (S 3 , S 4 ) = ( 2Ŝ 3 , 2Ŝ 4 ), 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 yield , P(0) = 0. (10) From (9c), we can immediately deduce thatS 3 is constant over this timescale. That is, over the late timescale,S 3 takes the constant valuẽ From (9)-(11), we can writeS 4 in terms of P as follows: Substituting (12) into (9), we obtain the following closed ODE for P(T ): where the critical parameterk * is defined as and we shall justify shortly the critical nature ofk * . The differential equation (13) is solved implicitly bȳ The asymptotic late-time results (11), (12), and (15) (dotted black lines) agree well with the numerical results (solid grey lines) when t = O(1/ ) = O(10 2 ) 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 S 2 , (9a) must be solved numerically, and the solution to this agrees well with the full numerical solution of S 2 in Fig. 2.
Given that the argument of the logarithm in (15) will vanish for some positive value of P if and only ifk 4 <k * , we may deduce that the large-time behaviour of P has a critical dependence on the sign ofk 4 −k * . Thus,k * is the bifurcation point for the parameterk 4 , in terms of changing the asymptotic structure of the solution. Specifically, whenk 4 <k * , corresponding to a weak reaction from malonic semialdehyde to 3HP, we obtain the bounded large-time behaviour . 3 The large-time limits of a 3HP production and b malonic semialdehyde varying over the critical parameterk 4 =k * . The grey lines denote the asymptotic predictions [from (16) and (18)], and the black crosses denote numerical results (solving (3b-e) with S 1 ≡ 1). 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 t = t end := 10 12 ; we approximate lim t→∞ d P/dt ≈ P(t end )/t end and lim t→∞ S 4 ≈ S 4 (t end ). Apart fromk 4 , 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 parameterk * ≈ 1.3951 however whenk 4 >k * , 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. It is of mathematical interest to briefly note that whenk 4 =k * , the solution to (13) is Hence, ask 4 passes through the bifurcation point, the long-time behaviour of P is to increase with the square root of time. From (12), we see thatS 4 will tend to a constant in the far-field; fork 4 <k * , we have the subcritical regimẽ and fork 4 >k * , we have the supercritical regimẽ 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 continuouspyruvate-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 regimek 4 >k * is much better for our goal. From Table 2, we note that the critical parameterk * ≈ 1.3951, butk 4 = 0.67. 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 k 4 = k * , where From (19), we see that an up-regulation of E 2 and a down-regulation of E 5 have the most significant effect on lowering k * , the critical bifurcation parameter, and thus hopefully moving into the supercritical regime for significantly increased 3HP production. Additionally, an up-regulation of E 3 will also have an effect in lowering k * , but with diminishing returns. Notably, our model suggests that regulation of neither E 1 nor E 4 has a significant effect on this critical bifurcation parameter (see below for a potential important role for E 4 , however). The supercritical regime occurs when k 4 > k * , 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 is and the dimensional large-time levels of malonic semialdehyde in the system are As 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 d[P]/dτ to [S 4 ], to obtain While (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 E 2 and E 4 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 E 3 and down-regulating E 5 , 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 ([S 4 ]) 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 S 1 ≡ 1. We use the parameter values in Table 2, but with increased enzyme concentrations as specified on the x-axis. While up-regulating E 2 yields improvement, this has diminishing returns. However, even though up-regulating E 4 by itself does not appear to have a significant effect on the ratio, up-regulating E 2 and E 4 in tandem does have a significant effect (Color figure online)

Asymptotic Structure
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 t = O(1/ ) (as shown in Fig. 5 for S 2 , S 3 , S 4 , and P), at which point the depletion of pyruvate occurs as a leading-order effect. Over this longer timescale, the rest of the system remains unchanged until the levels of pyruvate are of O( ), which we define to occur at t = T d / = O(1/ ). The pyruvate is then depleted to a negligible level at a relatively quicker rate, over a timescale of t − T d / = O(1), 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 t = O(1/ ). When the levels of acetyl-CoA become of O( ), which we define to occur at t = T * / = O(1/ ), 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 t = O(1/ ), and so t = T * / 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 T * .

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 S 1 ≡ 1, and the dashed blue lines are the numerical solutions from the no-replenishmentof-pyruvate system (3). We see that the two cases are essentially equivalent until around t = 5, after which the metabolite levels drop in a sharp manner for S 2 , S 3 , and S 4 , 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)

Reduced System for T *
The time at which the level of 3HP in the system is maximal is given by To determine T * 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 by Fig. 6 The dynamics of pyruvate (S 1 ) when it is never replenished, obtained from solving (3). We see that it becomes exponentially small shortly after t = 12. We use the parameter values given in Table 2 f Additionally, since the remaining metabolite concentrations are equivalent to those in the continuous-replenishment case until acetyl-CoA depletion occurs, the concentration of malonic semialdehydeŜ 4 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 T = T d , where T d satisfies the implicit equation The time at which there is a maximal amount of 3HP in the system T = T * 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 S 2 (T * ) = 0 or, equivalently, for T * 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 T * into (15) to yield the implicit equation 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 0 < t < T * / . Here, t = T * / is the asymptotic result for the point at which the 3HP in the system is maximal, and P(T * ) 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) Comparing 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 O(1) 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 E 5 (Fig. 8), the enzyme which recycles malonic semialdehyde back into acetyl-CoA. Moreover, we see that more 3HP is produced when E 5 is down-regulated. An alternative implicit equation for T * is given by (45), which uses the analytic solution for S 2 when T d < T < T * , stated in (44).

Asymptotic Expression for T * whenv 5 → 0
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 ratek 5 (Fig. 8). If this reaction is completely removed from the metabolic network, we are able to provide an analytic expression for T * and P(T * ). This corresponds to taking the limitk 5 → 0, equivalent tov 5 → 0. In this limit, there is no further production or depletion of 3HP when t > T * / , and the total amount of 3HP produced is given by P(T * ).
Asv 5 → 0, there are several scalings that reduce the parameter dependence of the system. We scale T =T /Â and S 2 (T ) =k 2 Y (T )/Â, to turn (23a) into The maximum amount of 3HP in the system as the concentration of E 5 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 E 5 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 E 5 where α =Â/k 2 2 , T * =T * /Â, the nonlinear functionf is defined as and T d =T d /Â has the form where β =k 2 /2, which will be useful later for notational purposes. When 0 <T <T d , the system (25a) is solved implicitly by We are able to obtain an explicit expression for Y (T d ) by re-writing (25c) as where we use (25a) forẎ ≡ dY /dT . Re-arranging (27), we obtain We can then obtain an explicit expression forT d by substituting (27)-(28a) into (26), to yieldT WhenT d <T <T * , system (25a) is solved explicitly by and asT * satisfies Y (T * ) = 0, we can stateT * in terms ofT d 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, v 5 → 0, noting thatk * also scales withv 5 . 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 F(α, β) in Fig. 9, but as a function of αβ 2 =Â/4 and β =k 2 /2, 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 α → 0, decreasing α (and henceÂ) to zero will have a more significant effect than any finite increase in β (and hencek 2 ), as long as the organism is able to survive with low values ofÂ.

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. Fig. 9 The asymptotic result for F(α, β), the total 3HP produced with no replenishment of pyruvate in the limit ofv 5 → 0, given in (31). We use axes 4αβ 2 and 2β, which are proxies for the strength of the sink and source in the system, respectively (Color figure online) The first key result is the time at which the 3HP levels are maximal, τ = τ * , where The effective production rate of 3HP is which can be increased by increasing E 2 or E 3 , though the latter has diminishing returns. The total amount of 3HP produced is given by where F is defined in (31). Defininĝ  (α,β), the asymptotic result for the total 3HP produced with no replenishment of pyruvate in the limit ofv 5 → 0, given in (31). In contrast to Fig. 9, we present F in terms of a variation in the parameters E 1 and E 2 (as shown in (35)), corresponding to regulation of enzymes (Color figure online) as proxies for E 1 and E 2 , we see that up-regulating E 2 results in higher levels of 3HP, whereas up-regulating E 1 results in lower levels of 3HP, with significant diminishing returns as E 1 increases (Fig. 10). The negative effect of up-regulating E 1 may appear counter-intuitive, since E 1 is required for a connected pathway. However, it can be explained by noting that the effect of increasing E 1 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 E 2 . 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 E 1 in the effective production rate (33). We emphasize that E 1 → 0 is a singular limit in the system, which appears because we are in the no-replenishment-of-pyruvate regime. As E 1 → 0, the earlytime 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.

Discussion
We develop and solve a mathematical model for the reaction kinetics of the 3hydroxypropionic 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. 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  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. material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons.org/licenses/by/4.0/.

A Depletion Dynamics
In the case where the pyruvate is not replenished, the system we consider is (3), and we no longer impose S 1 ≡ 1. However, in the no-replenishment case, the leading-order system is equivalent to the leading-order system for the continuous-replenishment case until t = O(1/ ). Thus, it is helpful to work in the long timescale T = t = O(1), over which the full system is given by (8b)-(8d), along with (36b)

A.1 Slow Depletion of Pyruvate: 0 < T < T d
Over the timescale T = O(1) and until the pyruvate becomes of O( ), the leadingorder system for S 2 ,S 4 ,S 3 , and P is still given by (9), the leading-order system for the continuous-replenishment case. The leading-order equation for S 1 is obtained from (36a) to yield which we can integrate to obtain the leading-order solution in terms of S 2 . The solution (38) is valid until S 1 = O( ), at which point a new asymptotic balance occurs. This occurs when T = T d (equivalently t = T d / ), where with an O( ) correction (equivalently O(1) in terms of t). Hence, T d can be obtained by solving dS 2 dT = 1 S 2 −k 2 +v 5Ŝ4 (T ) −ÂS 2 , with S 2 ∼ (2T ) 1/2 as T → 0 + , forward in time, using (12), (14)-(15) forŜ 4 , and terminating the calculation when S 1 = 0. Thus, this analysis is valid for T < T d .

A.2 Rapid Depletion of Pyruvate: t − T d / = O(1)
The next timescale is an interior layer (Holmes 2012), which occurs when t − T d / = (T −T d )/ =t = O(1). Over this timescale, S 1 varies from being of O( ) to exponentially small, while the remaining dependent variables are unchanged at leading order. Making the scaling S 1 = Y , the leading-order equation for the pyruvate concentration in this interior layer is given by with appropriate matching conditions ast → −∞. This depletion results in the exponential decay of Y , and hence S 1 , and is a slightly modified version of the coupled pyruvate and acetyl-CoA depletion dynamics considered in Appendix B in . As these depletion dynamics are not required to calculate details of 3HP production, we do not consider the dynamics of S 1 any further here for brevity; the analysis will be similar to that in  and require first-order correction terms to match to an appropriate accuracy. Since S 1 is exponentially small as we move forward in time out of this interior layer, we henceforth treat it as zero in the remaining timescales.

A.3 Slow Depletion of Acetyl-CoA: T d < T < T *
After the depletion of pyruvate, we remain in the timescale T = t/ = O(1), but now with T > T d . The system is thus given by (8b)-(8d), replacing (8a) with As each non-pyruvate dependent variable remains unchanged over the previous interior layer, the correct matching conditions are obtained by imposing continuity of each dependent variable across T = T d . The leading-order system is given by (9), replacing the first equation with dS 2 dT = −k 2 +v 5Ŝ4 −ÂS 2 .
This system is solved by (11)-(12) and (15) and these solutions are valid until S 2 = O( ), at which point a new asymptotic balance will occur. This occurs when T = T * (equivalently t = T * / ), where with an O( ) correction (equivalently O(1) in terms of t). After this point, the remaining metabolites will all start to deplete, so P(T * ) represents the maximal level of 3HP in the system with no replenishment of pyruvate. Physically, t = T * / gives the asymptotic "optimal" time at which to terminate a batch run and harvest the 3HP that has been produced.

A.4 The Remaining Depletion Dynamics
The remaining depletion dynamics are unimportant to the task of determining the maximal 3HP in the system, so we only state them briefly here for completion. The exact details of the transition of S 2 from O(1) to O( ) occur in another interior layer over the timescale t − T * / = (T − T * )/ = O(1), whereŜ 3 andŜ 4 also vary but do not change in asymptotic order. Over this timescale, the levels of P do not vary. In the forward-in-time far-field of this interior layer, the levels of each metabolite tend to a constant. The final depletion dynamics occur over the timescale T = t/ = O(1), but with T > T * . Over this timescale, the metabolites decay exponentially, governed by an ODE for P with the remaining metabolites coupled to the levels of P in a quasi-steady manner.