Perfect adaptation and optimal equilibrium productivity in a simple

Subscriber access provided by University of Rochester | River Campus & Miner Libraries

Article

Perfect adaptation and optimal equilibrium productivity in a simple microbial biofuel metabolic pathway using dynamic integral control Corentin Briat, and Mustafa H. Khammash ACS Synth. Biol., Just Accepted Manuscript • DOI: 10.1021/acssynbio.7b00188 • Publication Date (Web): 17 Jan 2018 Downloaded from http://pubs.acs.org on January 18, 2018

Just Accepted “Just Accepted” manuscripts have been peer-reviewed and accepted for publication. They are posted online prior to technical editing, formatting for publication and author proofing. The American Chemical Society provides “Just Accepted” as a free service to the research community to expedite the dissemination of scientific material as soon as possible after acceptance. “Just Accepted” manuscripts appear in full in PDF format accompanied by an HTML abstract. “Just Accepted” manuscripts have been fully peer reviewed, but should not be considered the official version of record. They are accessible to all readers and citable by the Digital Object Identifier (DOI®). “Just Accepted” is an optional service offered to authors. Therefore, the “Just Accepted” Web site may not include all articles that will be published in the journal. After a manuscript is technically edited and formatted, it will be removed from the “Just Accepted” Web site and published as an ASAP article. Note that technical editing may introduce minor changes to the manuscript text and/or graphics which could affect content, and all legal disclaimers and ethical guidelines that apply to the journal pertain. ACS cannot be held responsible for errors or consequences arising from the use of information contained in these “Just Accepted” manuscripts.

ACS Synthetic Biology is published by the American Chemical Society. 1155 Sixteenth Street N.W., Washington, DC 20036 Published by American Chemical Society. Copyright © American Chemical Society. However, no copyright claim is made to original U.S. Government works, or works produced by employees of any Commonwealth realm Crown government in the course of their duties.

Page 1 of 47 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

ACS Synthetic Biology

Perfect adaptation and optimal equilibrium productivity in a simple microbial biofuel metabolic pathway using dynamic integral control Corentin Briat and Mustafa Khammash∗ Department of Biosystems Science and Engineering, ETH Z¨ urich, Basel E-mail: [email protected]

Abstract The production of complex biomolecules by genetically engineered organisms is one of the most promising applications of metabolic engineering and synthetic biology. To obtain processes with high productivity, it is therefore crucial to design and implement efficient dynamic in-vivo regulation strategies. We consider here the microbial biofuel production model of Dunlop et al. (2010) for which we demonstrate that an antithetic dynamic integral control strategy can achieve robust perfect adaptation for the intracellular biofuel concentration in presence of poorly known network parameters and implementation errors in certain rate parameters of the controller. We also show that the maximum equilibrium extracellular biofuel productivity is fully defined by some of the network parameters and, in this respect, it can only be achieved when all the corresponding parameters are perfectly known. Since this optimum is a network property, it cannot be improved by the use of any controller that measures the intracellular biofuel

1 ACS Paragon Plus Environment

ACS Synthetic Biology 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

concentration and acts on the production of pump proteins. Additional intrinsic fundamental properties for the process are also unveiled, the most important ones being the existence of a conservation relation between the productivity and the toxicity, a low sensitivity of the optimal productivity with respect to a poor implementation of the set-point for the intracellular biofuel, and a strong intrinsic robustness property of the optimal productivity with respect to poorly known parameters. Taken together, these results demonstrate that a high and robust equilibrium rate of production for the extracellular biofuel can be achieved when the parameters of the model are poorly known and those of the controllers are poorly implemented. Finally, several advantages of the proposed dynamic strategy over a static one are also emphasized.

Keywords Antithetic integral control, perfect adaptation, metabolic engineering

Introduction Living microorganisms have the ability to turn cheap carbon sources, such as sugars, into highly complex biomolecules at a low energy cost (1 ). The tools developed in metabolic engineering and the field of synthetic biology (2 –4 ) allow us to take advantage of this ability for producing certain biomolecules of interest in low-cost industrial-scale bioreactors. Recent examples include drugs or drug precursors (5 –7 ), biofuels (8 –10 ) and antibodies (11 ). In order to optimize the efficiency of chemical reactors, control methods, such as model predictive control (12 ), have been traditionally used. In this case, control systems measure key process variables (such as temperature, pressure, and the concentration of various chemical species) using sensors and act back on the process through actuators after 2 ACS Paragon Plus Environment

Page 2 of 47

Page 3 of 47 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


a processing phase performed by the controller. This type of control algorithm is usually implemented in a computer, and is therefore referred to as in silico control. In the case of bioreactors (13 ), this type of control is also possible, with the difference being that reactions do not take place homogeneously, but are instead confined inside cells, thereby adding a layer of complexity to the problem (14 ). More recently, approaches relying on the use of genetically engineered cells and optogenetics (15 –18 ) or microfluidics (19 , 20 ) have been proposed as new ways for performing in silico control. There is, however, another control opportunity that synthetic biology allows for. Indeed, recent advances in this field suggest the possibility of implementing controllers in vivo through circuits that are engineered into living cells (21 –26 ) . While the implementation of such controllers is still an elusive task, the nascent general mathematical theory pertaining to the design and the analysis of such controllers at the molecular level is developing much faster (21 , 22 ). Whether in silico or in vivo control is used, the analysis and design of controllers that act at the genetic scale has given rise to a rich area of study – cybergenetics (22 ). In addition to synthetic in vivo and in silico regulation motifs, several endogenous motifs have also been discovered and analyzed both theoretically and experimentally; see e.g. (27 – 32 ). A recurring theme in many of these controllers is that the controller network is not clearly separated from the network it actually controls, as the boundary between process and controller is unclear. The point of view of control theory is to view control systems as an interconnection of a controller that regulates a separate system, and to precisely characterize the different controller networks that can be used and the control tasks they can achieve – a pure control theoretic perspective. This point of view is notably developed in (21 , 22 ) in the context of set-point regulation and perfect-adaptation in both the stochastic and the deterministic settings. In in vivo metabolic control, two general control strategies are often opposed: static



regulation and dynamic regulation (33 –36 ). The former does not use any information on the controlled metabolic networks and can be viewed as an open-loop control strategy; e.g. constitutive production. Hence, a simple strategy would be to simply overexpress the molecule of interest. However, it is known that overexpression leads to high toxicity and growth inhibition arising, for instance, from the overexpression of membrane proteins (37 ) or from the intrinsic toxicity of some metabolites (e.g. biofuel). Dynamic regulation allows one to circumvent this difficulty: by measuring some metabolites using biosensors, the dynamic controller implementing, for instance, a negative feedback loop, can act back on the metabolic pathway in order to optimize its overall productivity while, at the same time, reducing harm to the cell. Such control strategy is a closed-loop one and has been shown to be necessary if one wants to be robust to network uncertainties, cell-to-cell variability and environmental fluctuations; see e.g. (36 ). Here we propose to theoretically study a variation of the model of a biofuel production process considered in (8 ) where four different types of nonlinear static/dynamic regulation strategies were proposed and numerically analyzed via simulations. A major limitation of these controllers lies in their inability to ensure robust set-point tracking and perfect adaptation for the intracellular biofuel concentration. The perfect adaptation property is here of great importance as it will allow for the automatic compensation of unknown parametric uncertainties and of disturbances, which are systematically present due to model uncertainties, parameter fluctuations and unmodeled interactions. This motivates the consideration of a homeostatic controller implementing an integral action, a crucial element of proportionalintegral-derivative (PID) controllers that have been used to solve regulation problems in wide variety of contexts (38 ). The integral action is known for its (structural) tracking and perfect adaptation properties which are the main objectives of those kind of problems. Biological versions of this structure have also been discovered to be present in calcium regulation (29 ), glucose regulation (39 , 40 ) and in bacterial chemotaxis (28 ). It has also been posited in 4 ACS Paragon Plus Environment

Page 4 of 47

Page 5 of 47 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


(21 ) that the σ/anti-σ factors system in bacteria may endogenously implement such an integral control structure. On a more theoretical level, certain negative feedback and incoherent feedforward motifs have also been shown to exhibit (perfect) adaptation (31 ). The approach we consider here is a cybergenetic one, where we connect a controller network to another network (a metabolic network) with the goal of controlling it. The considered in vivo integral controller, referred to as antithetic integral controller (21 ), is different from the structure commonly used in control engineering, as it is nonlinear, takes only nonnegative state values, and is implementable in terms of biomolecular reactions with mass-action kinetics. It is assumed that the controller can directly sense the concentration of the intracellular biofuel and that it can act on the production rate of the pump proteins. Using this controller, several important properties for the closed-loop network are established. It is proved that the equilibrium point for the controlled dynamics that achieves tracking of the intracellular biofuel concentration to the desired set-point can be made locally exponentially stable provided that the controller parameters belong to a certain region which we partially characterize. The closed-loop network is also demonstrated to exhibit perfect adaptation for the intracellular biofuel concentration when the network is subject to changes in its parameters and to constant disturbances. We prove next that there exists a value for the intracellular biofuel set-point for which the equilibrium rate of production of the extracellular biofuel is maximum, which is interesting from an optimization standpoint. However, as this value depends on some of the network parameters, their perfect knowledge is required for this maximum to be achievable. Note also that this maximum value is a property of the system which cannot be improved through the use of any other types of controllers that measure the intracellular biofuel and act on the pump proteins expression. A robustness analysis is then carried out in order to consider more realistic scenarios where system parameters are poorly known and when the controller parameters are inaccurately implemented. In the former case, a robust optimization approach (41 ) allows us to compute the value for the set-point for which the 5 ACS Paragon Plus Environment


worst-case equilibrium productivity is maximized. In practice, this means that the equilibrium productivity will always be higher than this value, whence the term “worst-case”. In the latter case, a sensitivity analysis can be carried out in order to evaluate the performance deterioration. We numerically show that this deterioration is rather limited for the considered model and the considered parameters. In the process of proving all these results, a conservation relation between the equilibrium toxicity and the equilibrium productivity is unveiled and shown to be independent of the controller; i.e. it is an invariant of the system. A consequence of this relationship is twofold. Firstly, it implies that the value of the intracellular biofuel set-point that maximizes the equilibrium productivity minimizes the overall equilibrium toxicity. Secondly, it indicates that the productivity cannot be improved by considering a different controller structure since this productivity/toxicity relationship does not depend on the controller. Finally, we theoretically prove the existence of an intrinsic robustness property for the equilibrium productivity when considering a set-point that is proportional to the ratio of the normalized growth rate to the biofuel toxicity, a quantity that is easy to determine experimentally. This result notably states that the relative performance deterioration is at most 43.75% but also that this deterioration is quite low for a wide range of values for the parameters. At last, we discuss about the benefits and drawbacks of the proposed dynamic approach when compared with a static one in terms of sensitivity, disturbance rejection and robustness properties. It is notably shown that, unlike the static strategy, the proposed dynamic control approach has interesting robustness properties with respect to unknown leakiness, unknown pump expression models, unknown input channel structure and to constant disturbances within the pump and channel dynamics.


Page 6 of 47

Page 7 of 47 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


Results Biofuel production model The biofuel production model proposed in (8 , 9 ) and schematically depicted in Fig. 1 is given by n˙

= αn n(1 − n) − δn bi n −

αn np p + γp

b˙ i = αb n − δb pbi p˙

(1)

= u − βp p

b˙ e = V δb pbi n where n is the normalized cell density, bi is the concentration of intracellular biofuel, p is the concentration of efflux pump proteins, u is the rate of production of pump proteins, which will also be our control input, and be is the extracellular concentration of biofuel. The parameter αn is the specific growth rate of the cells, δn is the biofuel toxicity coefficient, αb is the biofuel production coefficient, δb is the biofuel export rate per pump, βp is the pump protein degradation rate, V is ratio of intra to extracellular volume and γp is the constant that sets the pump toxicity threshold. It is worth commenting on this model. First of all, it is a simplified model that does not include a full metabolic pathway but rather a simplified production pathway where the rate of production of intracellular biofuel is simply proportional to the cell-density n. Secondly, it does not capture the dynamics of the substrates, their intake (e.g. through diffusion and/or permeases) and the effect their potential shortage has on the pathway. Therefore, this model cannot be used to carry out a yield analysis which would require the computation of the ratio of the quantity of produced biofuel to the quantity of consumed substrate (34 ). As a consequence, additional tradeoffs that may arise, for instance, from the competition of permeases and pumps for space on the membrane, are not taken into account; see e.g. 7 ACS Paragon Plus Environment


(10 ) for a discussion on the simultaneous use of multiple types of efflux pumps. It is also important to note that the rate of production of the extracellular biofuel differs from the model in (8 ) and has been taken from (9 ). The motivation behind this modification is that bi needs to represent the concentration of biofuel inside a given cell (and not the overall quantity of intracellular biofuel across the entire population) since it will have to be locally monitored by the controller in order to regulate its concentration at the single-cell level. In this regard, the overall rate of production will be proportional to the export rate of a cell multiplied by the normalized cell-density (assuming that all cells are identical) and the ratio of intra to extracellular volume. Finally, the term αb n modeling the rate of production of the intracellular biofuel is here to, as stated in (9 ), capture ”the impact of cell viability on biofuel production” and should be understood on a phenomenological level rather than a mechanistic one. In spite of these limitations, this model remains interesting as it is simple and allows one to carry out an in-depth theoretical analysis of some of the intrinsic properties and limitations of biofuel production, many of which hold in more complete and accurate models, although they may be encapsulated in a different mathematical form. For example, an important feature of this model is that it captures the fundamental tradeoff between the cell-growth and pump/biofuel toxicity through the second and third terms in the first equation in (1). Modifying this model to incorporate more complex/accurate reactions including, for instance, saturations through Hill kinetics will simply impose constraints on the levels of intracellular biofuel and pump proteins that can be maintained within the cell. In a similar way, the competition between various membrane proteins will introduce additional tradeoffs that will impose, in turn, restrictions on the set of equilibrium points that can be reached by the system. In this regard, this model, although simple, can still be used to obtain interesting theoretical results and draw relevant conclusions which can then be used to gain intuition when facing more complex models. Finally, it is important to anticipate over the upcoming 8 ACS Paragon Plus Environment

Page 8 of 47

Page 9 of 47 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


results and mention now that the considered controller will be robust and will work even when more complex systems, unmodelled dynamics, etc. are present.

The control problem We define here the control problem that will be considered in the paper. We formulate this problem in control theoretical terms as control theory has the precise vocabulary for defining such problems. See Box 1 in (21 ) for more details about control theory in a biological context. See also (22 ) for a similar setup. The overall idea is to find a controller consisting of a set of additional species and reactions which, by interacting with the network described by the model (1), achieves certain desired properties for the closed-loop network (which consists of the interconnection of the network described by the model (1) and the chosen controller). This control problem is defined below: Problem. Find an in-vivo controller (i.e. that can be implemented in terms of biomolecular reactions) that measures the intracellular biofuel concentration and acts on the pump proteins production rate in a way that • [Asymptotic stability] makes the closed-loop network locally asymptotically stable (i.e. the trajectories of the closed-loop network starting sufficiently near the desired equilibrium point converge to it); • [Set-point tracking and perfect adaptation] ensures that the concentration of intracellular biofuel tracks a desired set-point µ > 0 and exhibits perfect adaptation properties (i.e. the equilibrium value of the intracellular biofuel concentration must be equal to the set-point µ even in presence of unknown certain constant inputs and changes in the parameters of the network);



• [Maximum productivity] maximizes the equilibrium rate of production of extracellular biofuel (i.e. the set-point µ should be chosen in a way that the equilibrium productivity – which will defined later – is maximum); • [Low metabolic load] limits the metabolic load on the host (i.e. the controller parameters should be chosen in a way that makes the energy cost of the controller reactions – which will also be defined later – low). • [Robustness properties] guarantees crucial robustness properties with respect to a poor knowledge of the parameters of the network and/or a poor implementation of the controller parameters (i.e. the design should ensure that despite poorly known network parameters and poorly implemented controller parameters, the closed-loop network still behaves in the desired way and that the performance deterioration is limited);

We now explain the rationale for defining such a control problem. The requirement on the asymptotic stability of the equilibrium stems from the fact that we would like the trajectories of the closed-loop network to be well-behaved. To this stability objective, we add the constraint that the equilibrium point of the closed-loop network be such that the equilibrium intracellular biofuel concentration equals a desired set-point value µ that we can choose. The closed-loop network should also exhibit the perfect adaptation property for the intracellular biofuel concentration. Note that we do not require here that the other molecular species or the cell-density adapt. Finally, the optimality objective is here to ensure that the controlled metabolic process is efficient while robustness is considered to account for process uncertainties and a poor implementation of the controller parameters.


Page 10 of 47

Page 11 of 47 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


Proposed integral controller Various controllers have been proposed in (8 ), namely, control through the constitutive production of pumps, biofuel-responsive (i.e. negative Hill-type nonlinear feedback) control, control through a repression cascade and control through a feedforward loop. The first one is static while the others are dynamic control strategies. The main problems with all these control strategies is that they do not ensure perfect adaptation and are difficult to tune to obtain optimal biofuel production. The latter issue stems from the fact that the controller networks involve a large number of network parameters whose effect on the dynamics and the equilibrium points is quite intricate. We claim here (and prove later) that the proposed network controller below alleviates these limitations: θ

Bi −−−→ Z2 + Bi ,

ν

∅ −−−→ Z1 ,

η

Z1 + Z2 −−−→ ∅,

k

Z2 −−−→ Z2 + P

(2)

where Z1 and Z2 are the controller species, Bi is the intracellular biofuel species and P is the pump protein species. The first reaction is the measurement reaction that catalytically creates one molecule of species Z2 from Bi with rate θ > 0. The second one is the reference reaction that implements part of the set-point value µ = ν/θ > 0 for the concentration of the intracellular biofuel Bi . The third reaction is the comparison reaction that correlates the species Z1 and Z2 and closes the overall control loop. Finally, the fourth and last reaction is the actuation reaction as it is the reaction that acts on the system in order to control it. The topology of the closed-loop network is depicted in Fig. 1. Note that the species Z2 plays the role of both the sensing and actuating species unlike in (21 ) where Z1 is the actuating species. In accordance with the number of reactions, this controller involves four parameters that can be tuned in order to meet the specified control objectives. The mathematical model (mass action kinetics) corresponding to the controller network



Page 12 of 47

(2) is given by z˙1 = ν − ηz1 z2 , z˙2 = θbi − ηz1 z2 , u

(3)

= kz2

which is the deterministic counterpart of the antithetic integral controller introduced in (21 ) with the only difference being that the control input u is now equal to kz2 . The reason for this stems from the fact that the (local) gain of the process around any of the possible equilibrium points is negative in the sense that when the input is increased, the output is decreased. Indeed, increasing u will increase the concentration of pump proteins which will, at the same time, a) decrease the cell-density by toxicity and hence decrease the production rate of the intracellular biofuel as it is proportional to the cell-density, and b) directly increase the extraction rate of the intracellular biofuel. In order to connect the above model to the usual integral controller, it is enough to observe that the difference system given by z˙1 − z˙2 = ν − θbi exactly behaves like the desired integrator; see (21 , 42 ). A potential implementation of the controller (2) relies on a network based upon σ-factors and anti σ-factors where Z2 would denote the σ-factor and Z1 the anti σ-factor (see Fig. 2). In this design, the anti σ-factor molecules would be constitutively produced at a rate equal to ν, a rate which could be easily modified using chemical inducers or optogenetics. Conveniently, the comparison reaction would be automatically implemented by exploiting the strong complementarity of the σ-factor and the anti σ-factor molecules known to have very high affinity (dissociation constant of 0.01 nM (43 , 44 )). Finally, the actuation reaction is performed by the RNA polymerase holoenzyme complex, which consists of a RNA polymerase molecule bound to a σ-factor molecule, that can specifically bind to the pump genes to initiate transcription. Even though the model does not include all the details of the overall control mechanism,


Page 13 of 47 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


it captures all its key elements. For instance, the dynamics of the RNA polymerase or the unbinding reactions between pairs of σ- and anti σ-factor molecules is ignored but it can be shown that they do not alter the integral action of the controller (see Section S10 of the SI). Generally speaking, the integral action is present as long as the species Z1 and Z2 do not degrade/dilute and tracking will be achieved as long as the set-point can be reached by the system (and that the closed-loop network is stable). The assumption of nondegrading/diluting species may not be realistic, especially for bacteria which have a short doubling time. However, a sensitivity approach will show that despite destroying the integral action of the controller, dilution has, in general, a mild impact on the performance of the controller. Finally, dilution can also be further compensated by considering a high-gain and a high-annihilation constant η; see Section S10.3 of the SI and (45 ).

Stability, tracking and perfect adaptation Stability of the equilibrium point. It is shown in Section S2.2 of the SI that the closedloop network (1)-(3) depicted in Fig. 1 admits a single positive equilibrium point provided that the set-point µ := ν/θ satisfies the growth-toxicity condition given by

µ < ρ :=

αn . δn

(4)

When this is not the case, no positive equilibrium point exists (e.g. negative or zero pump concentration) for the closed-loop dynamics, which results in a regime where the model has no real physical meaning. Interestingly, this inequality formulates an intrinsic limit on the maximum equilibrium intracellular biofuel concentration that can be reached, independently of the controller. The parameter ρ is defined here as the tolerance constant as it measures the relative tolerance of the organism for the biofuel. Note that this parameter can be easily



characterized using independent experiments where the intracellular biofuel concentration is fixed (8 ). Moreover, since this condition is independent of the controller, it can not be overcome by any other controller structure. When this condition is met, the positive-density equilibrium point can be made locally asymptotically stable provided that the controller parameters k, η > 0 are chosen to satisfy certain constraints which depend on the choice of the set-point parameters ν, θ > 0. In particular, if the parameter k, which is related to the actuation reaction speed, is chosen sufficiently small, then the equilibrium point of the closed-loop network will be locally asymptotically stable for any values of the parameter η > 0. If, on the other hand, the value of the parameter k lies beyond a certain threshold, then unstable/oscillatory trajectories for the closed-loop network will arise if the value of the parameter η also lies beyond a certain threshold, as illustrated in Fig. 3A (see Section S3 in the SI for more details about these conditions). This result is particularly important when the controller aims to be implemented using sigma and anti-sigma factors, as suggested in (21 ), since, in this case, the parameter η will be large due to the strong affinity between sigma and anti-sigma factor molecules. Note that the issue of having a large η may not arise when using other implementations of the antithetic controller, such as those involving small RNAs; see e.g. (45 –47 ). To conclude, the parameter k should be chosen such that it is not too small (which would result in a slow response due to a slow actuation reaction) and not too large (in order to preserve the stability of the controlled network). Even if not fully explored here, a potentially suitable value for the gain k would be the one that minimizes the spectral abscissa of the Jacobian matrix of the system; i.e. that minimizes the real-part of the right-most eigenvalues of the Jacobian matrix (see Proposition S3.1 in the SI for an explicit expression of the Jacobian matrix).

Tracking and perfect adaptation. When the positive-density equilibrium is locally asymptotically stable, the integral controller ensures that the equilibrium value of the in14 ACS Paragon Plus Environment

Page 14 of 47

Page 15 of 47 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


tracellular biofuel will be equal to the set-point value µ (see Fig. 4), regardless of the value of the network parameters. In this respect, no prior information about them is needed to achieve this property, which is a great advantage of integral controllers that many other controllers do not have; e.g. proportional ones (see Section S1.3.5 of the SI). When the network parameters change values, the integral controller will automatically react and bring back the intracellular biofuel concentration to the set-point µ, as shown in Fig. 4. This will also happen when the control input u is affected by constant disturbances of sufficiently small amplitude. Similar to tracking, perfect adaptation cannot be achieved using proportional controllers or Hill-type controllers (8 , 21 , 38 ) as they are very sensitive to parameter changes. For more details, see Section S4 in the SI.

Maximizing the biofuel production while limiting the metabolic load Productivity and toxicity measures. We propose here to consider a slightly different performance measure than is traditionally used in metabolic engineering and control. The yield is defined in (34 ) as the “mass of product formed per mass of substrate consumed (gproduct /gsubstrate ).” This performance measure cannot be considered here as the model does not take into account substrates. On the other hand, in the same reference, the productivity is defined as the “overall rate of production for the entire batch, the concentration of product per unit of time (g/(L·h))” (here, the quantity be (T )/T ) whereas the titer is defined as the “concentration of product at the end of a batch (gproduct /L)” (here, the quantity be (T )). These measures are related to each other in the sense that the former considers the rate of synthesis of the product (or metabolic flux) whereas the latter considers the final concentration of the product. The titer, equal here to be (T ) for some final time T , was considered in (8 ) as a 15 ACS Paragon Plus Environment


Page 16 of 47

comparative performance measure between several controllers using numerical simulations. Unfortunately, it is difficult to theoretically analyze it as its value is sensitive to initial conditions and to systems parameters. This is the reason why we consider here slightly different productivity concepts: the instantaneous productivity

Pi (t) := V δb bi (t)p(t)n(t) = b˙ e (t)

(5)

and the equilibrium productivity

Pe∗ (µ) := V δb b∗i (µ)p∗ (µ)n∗ (µ) =

V 2 ∗ 2 ∗ 2 δ b (µ) p (µ) αb b i

(6)

where b∗i (µ) = µ and p∗ (µ) are the equilibrium concentrations for the intracellular biofuel and pump proteins when the set-point is µ. When the positive-density equilibrium point is locally asymptotically stable and the system trajectories converge to it, then we get that Pi (t) → Pe∗ (µ) as t → ∞. Additionally, we also have the following relationship (assuming be (0) = 0) 1 be (T ) = T T

Z

T

T →∞

Pi (s)ds −−−→ Pe∗ (µ)

(7)

0

which reconciles the concepts of titer (i.e. be (T ) ≈ T Pe∗ (µ) for any sufficiently large time T ), productivity be (T )/T , instantaneous productivity Pi (t) and equilibrium productivity Pe∗ (µ) all together. An immediate advantage of the equilibrium productivity is that it only depends on the set-point µ (which makes it generic since there is no dependence on the transient behavior nor initial conditions of the model) and easy to analyze since it has a closed-form solution (see equation (S5.4) in the SI). In the light of the above discussion, it is clear that the equilibrium productivity is the most adapted concept for carrying out a theoretical analysis for the productivity of the process.


Page 17 of 47 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


Similarly, we define here the instantaneous toxicity (see also (8 )) as

Ti (t) :=

δn bi (t) + | {z } Biofuel toxicity

αn p(t) . γp + p(t) | {z } Pump toxicity

(8)

and the equilibrium toxicity as

Te∗ (µ) := δn b∗i +

α n p∗ . γp + p ∗

(9)

Under the assumption that the positive-density equilibrium point is locally asymptotically stable and the system trajectories converge to it, then we have that Ti (t) → Te∗ (µ) as t → ∞ where Te∗ (µ) only depends on µ and the network parameters.

A productivity-toxicity conservation law. An interesting underlying characteristic of the process, which does not depend on the controller, is the following productivity-toxicity conservation law: Te∗ (µ) + αn

Pe∗ (µ) V αb

1/2 = 1, µ ∈ (0, ρ)

(10)

which describes a nonlinear relationship between the equilibrium productivity and the equilibrium toxicity (see Proposition S5.6 in the SI). The interpretation of this formula is that low toxicity is necessarily associated with high productivity, and vice-versa. Though not surprising, such a relationship does not seem to have been captured analytically in the literature. Note, however, that the dependence on µ is nonlinear and that the evolution of the point (Pe∗ (µ), Te∗ (µ)) on the curve defined by (10) may not be monotonic. In fact, it will be indicated later that the function Pe∗ (µ) is non-monotonic. Finally, it is important to stress that this relation does not hold for the instantaneous productivity/toxicity as shown in Fig. 3B. where the evolution of the instantaneous productivity Pi (t) is plotted against the



Page 18 of 47

instantaneous toxicity Ti (t) for a piecewise constant trajectory for µ.

Optimal biofuel production. We are now interested in establishing whether it is possible to maximize the equilibrium productivity to optimize the process. It is important to stress right away that this value only depends on the process parameters and is independent of the controller. In this regard, it cannot be improved through the use of any other type of controllers measuring the intracellular biofuel concentration and acting on pump production. The only degree of freedom we have for the equilibrium productivity is through the parameter µ, hence we consider from now on the equilibrium productivity Pe∗ (µ) as a function of µ. It is shown in Section S5 of the SI that this function exhibits the global maximum

∗ Popt :=

V αb

αb δb γp ρ 1/2 (αb

!2

+ (δb γp ρ)1/2 )2

(11)

which is attained with the set-point value

µopt =

(αb δb γp ρ3 ) 1/2

(αb

1/2

+ (δb γp ρ)1/2 )2

∈ (0, ρ/4],

(12)

∗ hence Popt = Pe∗ (µopt ). A consequence of this result is that, by virtue of the productivity-

toxicity conservation law, this value also globally minimizes the equilibrium toxicity of the process, as shown in Fig. 3D. Hence, when the system trajectories converge to the positivedensity equilibrium point, then the instantaneous productivity converges to its optimal equilibrium value. An interesting property illustrated in Fig. 3D is that µopt ≤ ρ/4, which indicates some asymmetry in the equilibrium productivity Pe∗ (µ) as its maximum is always attained in the interval (0, ρ/4]. Finally, it is interesting to mention the following remarkably simple relationship between the maximum equilibrium productivity and the corresponding


Page 19 of 47 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


set-point ∗ Popt αb δb γp = . 2 µopt ρ

(13)

This expression will be turn out to be useful when evaluating the metabolic cost of the controller in the next section. Metabolic cost. We saw that by adjusting the set-point µ, it was possible to maximize the equilibrium productivity of the process. However, this should not be at the expense of a high metabolic load on the host which can be seen as another, yet unmodeled, toxicity source. It is hence necessary to evaluate the metabolic burden of the controller onto the host to show how to minimize it. Following the ideas (21 , 22 ), the metabolic cost can be evaluated through the concept of equilibrium power consumption defined for the current circuit as We∗ (ν, θ) = κr ν + κm θb∗i + κc ηz1∗ z2∗ + κa kz2∗ =

> `∗ (µ) := κa βp p∗ . κa βp p∗ (κr + κm + κc )ν + {z } | | {z } cost of adaptation constitutive limit

(14)

where κr , κm , κc and κa (unit is [J/M] (joule per molar)) are the elementary metabolic costs associated with the reference reaction, the measurement reaction, the comparison reaction and the actuation reaction, respectively (see Section S6 in the SI for more details). A low equilibrium power consumption is equivalent to saying that the consumed energy at equilibrium will increase at a low rate. By choosing a sufficiently small ν while preserving the ratio ν/θ = µ constant, we can approach the constitutive limit `∗ (µ) arbitrarily closely without changing the productivity as this property only depends on the ratio ν/θ. Interestingly, the constitutive limit is equal to the cost of the actuation reaction and coincides with the cost of a static constitutive control law that would lead to the same set-point for the equilibrium intracellular biofuel. The rest of the cost consists of the cost of adaptation and comes from the dynamic part of the controller. The function `∗ (µ) is a decreasing function of µ (see



Page 20 of 47

Section S6 of the SI) and a lower equilibrium power consumption is obtained by increasing µ, which is intuitive as more intracellular biofuel means less pumps. When µ = µopt , hence at optimality, we have that r

∗

` (µopt ) = κa βp

αb γp ρδb

(15)

which is the cost of a static controller that would achieve optimal equilibrium productivity (but with no adaptation property). Finally, it seems important to stress that even the above analysis is not exact when the reactions are implemented as multiple more elementary reactions, then the computed lower bound still holds as it only depends on the actuation reaction.

Robustness results Robustness and maximum worst-case equilibrium productivity. When the parameters of the process are poorly known and belong to some interval of values, the biofuel production can still be optimized using a robust optimization approach (41 ). If we indeed assume that the parameters (ρ, αb , δb , γp ) are poorly known and belong the box

P := [ρ, ρ] × [αb , αb ] × [δb , δb ] × [γp , γp ],

(16)

then the formula (12) can intriguingly still be used by simply substituting the parameter values by their lower bounds in order to get a “robustified” version of µopt , denoted by µrob , which is proved in Theorem S8.1 in the SI to be given by µrob :=

3

αb δb γp ρ

1/2

αb 1/2 + (δb γp ρ)1/2

2 ∈ (0, ρ/4].


(17)

Page 21 of 47 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


The above result is a consequence of the fact that the equilibrium productivity is an increasing function of the network parameters. In this regard, µrob is such that the worst-case productivity is maximum (see for instance Fig. 5A). From a practical viewpoint, this means that the equilibrium productivity will always be higher than this worst-case value. For more details, see Section S8 of the SI.

Sensitivity analysis with respect to set-point implement errors. Optimum equilibrium productivity is only possible when the set-point is exactly µopt , which is not practically feasible. It is therefore interesting to quantify the performance deterioration, in terms of loss of equilibrium productivity, that is associated with a set-point implementation error. It is interesting to mention that the same sensitivity analysis can be used to quantify the deterioration effect of dilution terms on the controller species which would destroy the integral action and induce some steady-state error for the equilibrium biofuel concentration. Let ε be the set-point relative implementation error (i.e. ε = (µ − µopt )/µopt ) or the relative steady-state error (i.e. ε = (b∗i − µopt )/µopt ) in the case of a leaky integrator having degrading/diluting species. The numerical results for the sensitivity are depicted in Fig. 5B where we can see that, for the considered parameters, in the process is quite robust with respect to set-point implementation errors since an implementation error of 80% roughly preserves 50% of the maximum performance.

A simple robust suboptimal set-point. The main issue with the implementation of µopt in (12) or µrob in (17) is that they require some knowledge about four parameters. The former requires the exact knowledge of the parameters ρ, αb , δb , γp whereas the latter requires the knowledge of their respective lower-bound. However, some of these parameters or even their range may be very difficult to infer from experiments. In this regard, it would be



Page 22 of 47

convenient to have a value for the set-point µ that has less dependency on the parameters while preserving at the same time most of the equilibrium productivity. We propose the following reference value ρ µρ = , 4

(18)

which obviously only depends on the parameter ρ = αn /δn , which is known to be easily identifiable (8 ). The rationale for considering this value is that, given ρ, it is the maximum value for µopt . To evaluate the performance deterioration in terms of the equilibrium productivity, we define the relative error as

ρ :=

∗ − P ∗ (µρ ) Popt ≥ 0. ∗ Popt

(19)

It is proved in Section S9 of the SI that this error satisfies the following inequality

0 ≤ ρ
0. In other words, by implementing µρ in place of µopt , we can only lose at most 43.75% of the performance. This fact is illustrated in Fig. 6. When the parameters are uncertain, we simply need to implement the value

µρ,rob =

ρ 4

(21)

and we will only lose at most 43.75% of the maximum worst-case performance P ∗ (µrob ). In particular, we can see in Fig. 5 (see also Section S9 of the SI) that for the considered values for the parameters, the performance degradation is less than 1%. This demonstrates that high equilibrium productivity can even be achieved with little knowledge about the system.


Page 23 of 47 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


Comparison with static (open-loop) control It is interesting to compare the proposed dynamic control strategy with a static one where the control input u is simply set to a constant value. In this case, the expression for the open-loop equilibrium productivity is given by Pe (u) = V αb

δb αn βp γp u (u + βp γp )(αn δb u + βp αb δn )

2 (22)

and is maximized with the control input value r uopt = βp

γp α b δ n = βp αn δb

r

γp α b . ρδb

(23)

For comparison purposes, we select the parameters θ = 1, ν = µopt = 0.1762, k = 0.1, η = 0.1 and u = uopt = 0.1297 and we obtain the simulation results depicted in Fig. 7. In Fig. 7A. we can see that the extracellular biofuel production is higher when considering our dynamic control strategy, mostly because of the transient phase of the dynamic control strategy. In this regard, we cannot really assess what strategy is the best with this performance measure. However, we can observe in Fig. 7B. that the sensitivity of our dynamic control strategy is lower for negative error values than in the static case, but the opposite effect occurs for positive errors. However, it is worth mentioning that, in practice, the set-point µ or the control input u could be more or less precisely adjusted using, for instance, optogenetics and, perhaps to a lesser extent, chemical inducers. In this case, large errors may be ruled out and, if we assume that the maximum error is a 50% error, then we can immediately see that the performance degradation of both control strategies remain quite close to each other. On the other hand, the proposed dynamic control strategy is robust with respect to the value of the pump degradation rate βp as shown in Fig. 7C. where it can be seen that the dynamic control strategy can automatically restore the maximum performance after a change in βp . 23 ACS Paragon Plus Environment


Page 24 of 47

This can not be achieved by the static control law as it explicitly depends on the value βp . Interestingly, this property generalized to more complex input functions such as saturated Hill functions or more complex functions arising from the presence of intermediary reactions between the output of the controller and the effective pump production rate; e.g. binding with RNA polymerase, production of mRNA molecules, pump protein maturation to an active form, etc. This can be incorporated in the model using the steady-state relationship p = h(u) where h(0) = 0 and h is a monotonically increasing function, which may be bounded or not. While the integral controller will be able to be robust with respect to all the parameters involved in the expression of the function h, the open-loop control will not since the equilibrium productivity Pe (u) = V αb

δb αn γp h(u) (h(u) + γp )(αn δb h(u) + αb δn )

2 (24)

is only maximized with the optimal input

−1

r

uopt = h

γp α b ρδb

(25)

which, obviously, depends on the parameters involved in the expression of the function h. Along the same lines, we can see in Fig. 7D. that the dynamic control strategy will be able to compensate for the presence of leaky pump promoter by suitably adapting the value of the control input. The static controller will not be able to achieve this property. Besides these differences in terms of performance between the two controller structures, it is interesting to mention that the two approaches lead to very different ways to view the considered network. For instance, working directly with the control input would have made much more complicated the discovery of the approximate optimal set-point µρ which only depends on ρ. The worst-case analysis leading to a value for µrob would have also be made


Page 25 of 47 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


more complicated.

Discussion Integral control has been a dominant control strategy in many engineering disciplines as it can control a process in a way that ensures output tracking and constant disturbance rejection properties. In spite of the slight differences with the usual integral control, we proved that the proposed one possesses the same properties, thereby demonstrating its potential usefulness in synthetic biology and metabolic control. As demonstrated in this paper, this type of controllers has the potential of leading to optimized single-cell behavior in the deterministic setting, which parallels previously obtained results in the case of stochastic cell populations (21 ). We have notably shown that the intracellular biofuel concentration setpoint can be optimally chosen so that the equilibrium productivity is maximum. Notably, the value of the set-point that yields this maximum only depends on the parameters of the network, demonstrating that the type of considered controller is irrelevant and that those parameters need to be perfectly known to achieve optimal productivity. This indicates why the numerically optimized controllers considered in (8 ) are also leading to similar values for the equilibrium productivity. Note, however, that this fact was not reported in (8 ) as the considered performance criterion was the extracellular biofuel concentration at some stopping time T . This emphasizes the fact that such controllers can be considered in order to yield high equilibrium productivity. Unfortunately, they cannot ensure perfect adaptation properties and are difficult to tune due to a large number of intertwined parameters. The proposed controller has been shown to ensure robust perfect adaptation for the intracellular biofuel concentration and to be tunable in a systematic and simple way. The implementation can be also made robust in the sense that implementation errors and poorly known system parameters – some issues that are very likely to occur in biology – may be tackled 25 ACS Paragon Plus Environment


in a rigourous way so as to reduce the performance deterioration. In particular, the proposed dynamic controller can compensate for the leakiness of promoters and the presence of unknown reactions between the output of the controller and the effective production rate of pump proteins. These properties are not achieved by any other controller structure different from those involving an integral action. Even if precise tuning could be done a single-cell level when scaling the problem up to the population level, cell-to-cell variability will cause serious issues as no single common set-point that maximizes the equilibrium productivity for all the cells at the same time. In this case, a complementary in-silico control strategy using population-level measurements should be considered in order to find the right common set-point for the overall cell-population. We have also shown that many properties of the model can be stated in terms of the tolerance constant ρ which can be easily determined experimentally (8 ). The maximum admissible set-point is indeed equal to ρ and the maximum value for the set-point µopt that maximizes the equilibrium productivity is equal to ρ/4. Even more surprisingly, by implementing this upper-bound in place of µopt , we preserve at least 56.25% of the optimal equilibrium productivity. It is unclear at this point whether these robustness properties are generic to metabolic network models or are specific to the considered model. The worst case scenario is rather extreme and in the present case, the deterioration of performance is less than 1%. It is also unclear what is the source for such robustness properties or, in other words, what mechanisms in the model yield such robustness properties. All the proven results suggest the existence of multiple design procedures for the setpoint and the parameters of the controllers. These procedures are represented in Fig. 8. All of them start with a parameter identification phase which can be either total and exact (all the parameters are perfectly known), total and approximate (lower bounds for all the parameters are known) or even partial and exact (the exact value of the tolerance constant


Page 26 of 47

Page 27 of 47 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


ρ) is known. Note that this tolerance constant can be directly obtained from some simple experiments as in (8 ). Once this is done, we can select the associated set-point with the considered scenario. If energy consumption is an issue then ν should be chosen small and θ such that θ = ν/µ. Note, however, that when ν is small, then system may be slow and there will be likely a tradeoff between energy consumption and speed of the system. Finally, the parameters k, η should be chosen so that stability for the closed-loop network is achieved. It is possible that η is not tunable because it may be depend on the affinity of two molecules such as a σ factor and its associated anti-σ factor. In this case, k can always be chosen small enough so that stability is achieved. In the process of proving the main results of the paper, an interesting fundamental property of the model, the productivity-toxicity conservation law, has been unveiled. This relationship indicates that a low toxicity is always associated with a high productivity through a nonlinear relationship. The existence of such a relationship describing a tradeoff between the productivity and the toxicity was expected but was never theoretically reported so far. It is therefore tempting to posit that an analogous relationship exists for other metabolic networks, perhaps in a more complicated form. For instance, assume that instead of having a linear relationship between the equilibrium cell-density n∗ and the equilibrium productivity Pe∗ (µ) we have the nonlinear relationship Pe∗ (µ) = αb f (n∗ )

(26)

where f (n) is any nonnegative monotonically increasing function. Such a relationship may be obtained, from instance, by solving for the equilibrium point of a cascade of reactions with Hill kinetics. Then, in such a case, the conservation law becomes

f

−1

Pe∗ (µ)

1/2 !

V αb

Te∗ (µ) + =1 αn


(27)


where f −1 is the inverse function of f , which is also monotonically increasing. Note that when f (n) = n, then we retrieve the nonlinear relationship obtained in this paper. Hence, in this more general case also, maximum equilibrium productivity will be associated with minimum equilibrium toxicity. It is important to stress here that this result is essentially of theoretical nature and that it would be interesting to have an experimental confirmation of the existence of such a relationship. As far as theoretical models are concerned, a novel, simple and generic class of models based on this concept of conservation law which may be useful for the study and the design of metabolic networks could be defined. This class of models would be the metabolic network analogue of the class of gene regulatory network models which have been extensively studied over the past recent years. Those generic models could be referred to as productivity-toxicity models in which toxicity terms limit growth whereas productivity depends on cell-density, either directly or indirectly through intermediate metabolites. It is possible that this class of models would lead to a better mathematical understanding of the performance tradeoffs and the robustness properties of metabolic networks. Perhaps more optimistically, the understanding of the limitations on the productivity imposed by productivity-toxicity conservation laws could lead to a rational design of synthetic heterologous metabolic networks in order to maximize the achievable maximum productivity (48 ). This can be done by choosing host organisms with high tolerance (i.e. high ρ in the current model) but by also designing networks that do not saturate too quickly. It is also possible that this can be combined with automated computational tools for the design of heterologous metabolic networks; see e.g. (49 ). It is also important to mention here that the proposed analysis is preliminary as many phenomena and elements have been neglected in the model. For instance, butanol may diffuse in and out the cell directly through the membrane and, hence, the extracellular butanol will


Page 28 of 47

Page 29 of 47 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


also be toxic to the cell, a phenomenon that is neglected here. We also assume here that the controller can directly sense the butanol and acts on the production of the pump proteins. However, this does not seem to be practically feasible since it is unclear whether a biosensor can be an integrator and a transcription factor at the same time. Adding a biosensor is clearly necessary here and the mathematical model should be updated to account for its presence in order to capture its dynamics, its sensitivity and its range. The dilution of the controller species has been also shown to induce a steady-state error and, as a consequence, a productivity deterioration. However, a sensitivity analysis of the productivity demonstrated that this deterioration is, in general, quite mild, even for large steady-state errors. A steadystate error of 50% only yields a 10% decrease of the maximum productivity. A well-known limitation of integral control lies on its destabilizing effect that may lead to oscillatory trajectories and high overshoot. A possible solution to this problem is to adjoin a proportional action (one of the motifs proposed in (8 ) for instance) in order to improve the behavior of the trajectories. This control strategy would be analogous to a proportional-integral (PI) controller and would combine the advantages of both controllers: a fast and stable response (proportional action) and the properties of set-point tracking and perfect adaptation (integral action). Interestingly, proportional action has also been shown to reduce variability in stochastic models (50 , 51 ) and may also be helpful when controlling cell-populations. Note, finally, that since all the properties of the considered metabolic network do not depend on the type of integral controller that is considered, it may be substituted by any other integral controller.

Associated content Supporting Information. Supplemental methods, discussion, and figures.



Author information All the authors are with the Department of Biosystems Science and Engineering (ETHZ¨ urich), Mattenstrasse 26, 4058 Basel, Switzerland. Corresponding author. Email: [email protected] Author contribution. C.B. and M.K devised and performed the research, and wrote the paper. Notes. The authors declare no competing financial interest.

Acknowledgement The authors acknowledge funding support from the Swiss National Science Foundation grant 200021-157129.

References 1. Alberts, B., Johnson, A., Lewis, J., Morgan, D., Raff, M., Roberts, K., and Walter, P. Molecular Biology of the Cell, 6th ed.; Garland Science, 2014. 2. Stephanopoulos, G. (2012) Synthetic Biology and Metabolic Engineering. ACS Synth. Biol. 1(11), 514–525. 3. Zhao, H., Ed. Synthetic Biology. Tools and Applications; Elsevier, 2013. 4. Jullessona, D., Davida, F., Pflegerb, B., and Nielsena, J. (2015) Impact of synthetic biology and metabolic engineering on industrial production of fine chemicals. IEEE Trans. Autom. Control 33(7), 1395–1402. 30 ACS Paragon Plus Environment

Page 30 of 47

Page 31 of 47 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


5. Martin, V. J. J., Pietra, D. J., Withers, S. T., Newman, J. D., and Keasling, J. D. (2003) Engineering a mevalonate pathway in Escherichia coli for production of terpenoids. Nat. Biotechnol. 21(7), 796–802. 6. Ro, D. K., Paradise, E. M., Ouellet, M., Fisher, K. J., Newman, K. L., Ndungu, J. M., Ho, K. A., Eachus, R. A., Ham, T. S., Kirby, J., Chang, M. C. Y., Withers, S. T., Sarpong, Y. S. R., and Keasling, J. D. (2006) Production of the antimalarial drug precursor artemisinic acid in engineered yeast. Nature 440(7086), 940–943. 7. Galanie, S., Thodey, K., Trenchard, I. J., Filsinger Interrante, M., and Smolke, C. D. (2015) Complete biosynthesis of opioids in yeast. Science 349(6252), 1095–1100. 8. Dunlop, M. J., Keasling, J. D., and Mukhopadhyay, A. (2010) A model for improving microbial biofuel production using a synthetic feedback loop. Syst. Biol. Synth. Biol. 4(2), 95–104. 9. Harrison, M. E., and Dunlop, M. J. (2012) Synthetic feedback loop model for increasing microbial biofuel production using a biosensor. Front Microbiol 3, 360. 10. Turner, W. J., and Dunlop, M. J. (2015) Trade-offs in improving biofuel tolerance using combinations of efflux pumps. ACS Synth. Biol. 4(10), 1056–1063. 11. Lienert, F., Lohmueller, J. J., Garg, A., and Silver, P. A. (2014) Synthetic biology in mammalian cells: next generation research tools and therapeutics. Nat. Rev. Mol. Cell Biol. 15, 95–107. 12. Nikolaou, M. Model predictive controllers: A critical synthesis of theory and industrial needs; Advances in Chemical Engineering; Academic Press, 2001; Vol. 26; pp 131–204. 13. Bastin, G., and Dochain, D., Eds. On-line Estimation and Adaptive Control of Bioreactors; Elsevier, Amsterdam, 1990. 31 ACS Paragon Plus Environment


14. Elowitz, M. B., Levine, A. J., Siggia, E. D., and Swain, P. S. (2002) Stochastic gene expression in a single cell. Science 1183, 1183–1186. 15. Milias-Argeitis, A., Summers, S., Stewart-Ornstein, J., Zuleta, I., Pincus, D., ElSamad, H., Khammash, M., and Lygeros, J. (2011) In silico feedback for in vivo regulation of a gene expression circuit. Nat. Biotechnol. 29, 1114–1116. 16. Briat, C., and Khammash, M. Computer control of gene expression: Robust setpoint tracking of protein mean and variance using integral feedback. 51st IEEE Conference on Decision and Control. Maui, Hawaii, USA, 2012; pp 3582–3588. 17. Briat, C., and Khammash, M. Integral population control of a quadratic dimerization process. 52nd IEEE Conference on Decision and Control. Florence, Italy, 2013; pp 3367– 3372. 18. Milias-Argeitis, A., Rullan, M., Aoki, S. K., Buchmann, P., and Khammash, M. (2016) Automated optogenetic feedback control for precise and robust regulation of gene expression and cell growth. Nat. Commun. 7(12546), 1–11. 19. Uhlendorf, J., Miermont, A., Delaveau, T., Charvin, G., Fages, F., Bottani, S., Batt, G., Bottani, S., Batt, G., and Hersen, P. (2012) Long-term model predictive control of gene expression at the population and single-cell levels. Proc. Natl. Acad. Sci. USA. 109(35), 14271–14276. 20. Menolascina, F., Fiore, G., Orabona, E., De Stefano, L., Ferry, M., Hasty, J., di Bernardo, M., and di Bernardo, D. (2014) In-Vivo Real-Time Control of Protein Expression from Endogenous and Synthetic Gene Networks. PLoS Comput. Biol. 10(5), e1003625. 21. Briat, C., Gupta, A., and Khammash, M. (2016) Antithetic Integral Feedback ensures robust perfect adaptation in noisy biomolecular networks. Cell Systems 2, 17–28. 32 ACS Paragon Plus Environment

Page 32 of 47

Page 33 of 47 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


22. Briat, C., Zechner, C., and Khammash, M. (2016) Design of a synthetic integral feedback circuit: dynamic analysis and DNA implementation. ACS Synthetic Biology 5(10), 1108– 1116. 23. Del Vecchio, D., and Murray, R. M., Eds. Biomolecular Feedback Systems; Princeton University Press, 2015. 24. Franco, E., Giordano, G., Forsberg, P.-O., and Murray, R. M. (2014) Negative Autoregulation Matches Production and Demand in Synthetic Transcriptional Networks. ACS Synth. Biol. 3(8), 589–599. 25. Hsiao, V., de los Santos, E. L. C., Whitaker, W., Dueber, J. E., and Murray, R. M. (2015) Design and Implementation of a Biomolecular Concentration Tracker. ACS Synth. Biol. 4(2), 150–161. 26. Del Vecchio, D., Abdallah, H., Qian, Y., and Collins, J. J. (2017) A Blueprint for a Synthetic Genetic Feedback Controller to Reprogram Cell Fate. Cell Systems 4(1), 109– 120. 27. Alon, U., Surette, M. G., Barkai, N., and Leibler, S. (1999) Robustness in bacterial chemotaxis. Nature 397, 168–171. 28. Yi, T.-M., Huang, Y., Simon, M. I., and Doyle, J. (2000) Robust perfect adaptation in bacterial chemotaxis through integral feedback control. Proc. Natl. Acad. Sci. USA. 97(9), 4649–4653. 29. El-Samad, H., Goff, J. P., and Khammash, M. (2002) Calcium homeostasis and parturient hypocalcemia: An integral feedback perspective. J. Theor. Biol. 214, 17–29. 30. El-Samad, H., Kurata, H., Doyle, J. C., Gross, C. A., and Khammash, M. (2005) Sur-



viving heat shock: Control strategies for robustness and performance. Proc. Natl. Acad. Sci. USA. 102(8), 2736–2741. 31. Ma, W., Trusina, A., El-Samad, H., Lim, W. A., and Tang, C. (2009) Defining network topologies that can achieve biochemical adaptation. Cell 138, 760–773. 32. Shinar, G., and Feinberg, M. (2010) Structural sources of robustness in biochemical reaction networks. Science 327, 1389–1391. 33. Berkhout, J., Teusink, B., and Bruggeman, F. J. (2013) Gene expression requirements for regulation of metabolic gene expression to a desired state. Sci. Rep. 3, 1417. 34. Venayak, N., Anesiadis, N., Cluett, W. R., and Mahadevan, R. (2015) Engineering metabolism through dynamic control. Curr. Opin. Biotechnol. 34, 142–152. 35. Cress, B. F., Trantas, E. A., Ververidis, F., Linhardt, R. J., and Koffas, M. A. G. (2015) Sensitive cells: enabling tools for static and dynamic control of microbial metabolic pathways. Curr. Opin. Biotechnol. 36, 205–214. 36. Chukubov, V., Mukhopadhyay, A., Petzold, C. J., Keasling, J. D., and Mart´ın, H. G. (2016) Synthetic and systems biology for microbial production of commodity chemicals. Systems Biology and Applications 2, 16009. 37. Wagner, S., Baars, L., Ytterberg, A. J., Klussmeiert, A., Wagner, C. S., Nord, O., Nygren, P.-˚ A., van Wijk, K. J., and de Gier, J.-W. (2007) Consequences of Membrane Protein Overexpression in Escherichia coli. Mol. Cell. Proteomics 6.9, 1527–1550. 38. ˚ Astr¨om, K. J., and H¨agglund, T. PID Controllers: Theory, Design, and Tuning; Instrument Society of America: Research Triangle Park, North Carolina, USA, 1995.


Page 34 of 47

Page 35 of 47 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


39. Topp, B., Promislow, K., DeVries, G., Miura, R. M., and Finegood, D. T. (2000) A Model of β-Cell Mass, Insulin, and Glucose Kinetics: Pathways to Diabetes. J. Theor. Biol. 206, 605–619. 40. Karin, O., Swisa, A., Glaser, B., Dor, Y., and Alon, U. (2016) Dynamical compensation in physiological circuits. Mol. Syst. Biol. 12, 1–7. 41. Bental, A., El Ghaoui, L., and Nemirovski, A. Robust Optimization; Princeton University Press, 2009. 42. Oishi, K., and Klavins, E. (2010) Biomolecular implementation of linear I/O systems. IET Systems Biology 5(4), 252–260. 43. Maeda, H., Fujita, N., and Ishihama, A. (2000) Competition among seven Escherichia coli σ subunits: relative binding affinities to the core RNA polymerase. Nucleic Acids Res 28(18), 3497–3503. 44. Bakshi, S., Siryaporn, A., Goulian, M., and Weisshaar, J. C. (2012) Superresolution imaging of ribosomes and RNA polymerase in live Escherichia coli cells. Mol. Microbiol. 85(1), 21–38. 45. Qian, Y., and Del Vecchio, D. Realizing “integral control” in living cells: How to overcome leaky integration due to dilution? bioRxiv 46. Levine, E., Zhang, Z., Kuhlman, T., and Hwa, T. (2013) Quantitative characteristics of gene regulation by small RNA. PLoS Biol 5(9), e229. 47. Yoo, S. M., Na, D., and Lee, S. Y. (2013) Design and use of synthetic regulatory small RNAs to control gene expression in Escherichia coli. Nat. Protocols 8(9), 1694–1707.



48. Yadav, V. G., De Mey, M., Lim, C. G., Ajikumar, P. K., and Stephanopoulos, G. (2012) The future of metabolic engineering and synthetic biology: Towards a systematic practice. Metab. Eng. 14(3), 233–241. 49. Carbonell, P., Fichera, D., Pandit, S. B., and Faulon, J.-L. (2012) Enumerating metabolic pathways for the production of heterologous target chemicals in chassis organisms. BMC Syst. Biol. 6(10), 1–18. 50. Milias-Argeitis, A., Engblom, S., Bauer, P., and Khammash, M. (2015) Stochastic focusing coupled with negative feedback enables robust regulation in biochemical reaction networks. J. R. Soc., Interface 12, 20150831. 51. Oyarzun, S. A., Lugagne, J.-B., and Stan, G.-B. V. (2015) Noise Propagation in Synthetic Gene Circuits for Metabolic Control. ACS Synth. Biol. 4(2), 116–125.

Figure legends • [Legend Fig. 1]The controlled metabolic pathway. The metabolic pathway (in the cloud) consumes sugars in order to produce some butanol (not included in the model). Since butanol is toxic to the cell, we need to export it using pumps that are located on the membrane of the cells. The production of pumps is controlled by an antithetic integral controller that measures the current concentration of biofuel and acts back on the pump expression accordingly. • [Legend Fig. 2] Practical implementation of the antithetic integral controller using σfactors and anti σ-factors (denoted by σ). The reference reaction is implemented as the constitutive production of anti σ-factor molecules. This reaction can be tuned using chemical inducers or optogenetics. The production of σ-factor molecules is activated 36 ACS Paragon Plus Environment

Page 36 of 47

Page 37 of 47 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


by the intracellular biofuel, possibly through the use of an intermediary biosensor. The comparison reaction is automatically realized through the complementarity property of the σ-factor and the anti σ-factor molecules. The free σ-factor molecules can then bind to the RNA polymerase to form a holoenzyme complex that specifically binds to the promoter region of the pump genes to initiate transcription. • [Legend Fig. 3] A. Evolution of the intracellular biofuel concentration bi (t) with initial condition x(0) = (0.01, 0, 0, 0, 0), k = 0.3289, θ = 1 and ν = µopt = 0.1762 (in black dots) for different values of η. When η is smaller than η ≈ 2.99, then the trajectories converge to the equilibrium value while when η is larger than η, then a stable limit cycle arises. B. Evolution of the instantaneous productivity versus the instantaneous toxicity for some piecewise constant trajectory for ν which successively takes the values 0.0145 → 0.1813 (section 1 in blue), 0.1813 → 0.3626 (section 2 in red), 0.3626 → 0.5440 (section 3 in orange) and 0.5440 → 0.6527 (section 4 in purple). We can see that after some transient phase, the trajectories always converge back to the curve described by the conservation law (10). The fact that the trajectories move in both directions illustrate the non-monotonic behavior of the productivity and toxicity as functions of µ. C. Bifurcation curve in the (k, η)-plane for θ = 1, ν = µopt = 0.1762. When the pair of parameters (k, η) is chosen to lie beneath the curve, the equilibrium point is locally asymptotically stable but comes oscillatory when the pair is located above the bifurcation curve. D. Evolution of the normalized equilibrium productivity ∗ ∗ Pe∗ (µ)/Popt and the normalized equilibrium toxicity T ∗ (µ)/Topt as a function of µ. We

can clearly see that when µ = µopt = 0.1762 the productivity is maximized whereas the toxicity is minimized, which is an immediate consequence of the productivity-toxicity formula (10). • [Legend Fig. 4] Response of bi (t) to a variation in the parameter αb (A), δb (B), γp



(C) and βp (D) at t = 200h with initial condition x(0) = (0.01, 0, 0, 0, 0), θ = 0.3, ν = 0.0529 (i.e. µ = 0.1762), k = 0.12 and η = 0.1. The red curve corresponds to a variation of −50% of the parameter while the blue one corresponds to a variation of +100% of the parameter. In all these scenarios, the process shows perfect adaptation, thereby demonstrating the importance of integral control. • [Legend Fig. 5]A. Evolution of the equilibrium productivity Pe∗ (µ) as a function of µ and for different values of ρ (the other parameters are fixed). We can see that, in accordance with the robustness result (Theorem S8.1 in the SI) and the formula (17), the value µrob maximizes the worst-case productivity. B. Evolution of the positive and − negative maximum relative errors ε+ σ and εσ in the implementation of µ = µopt as a

function of the performance deterioration σ. We can see that an error of 50% in the implementation of µ = µopt = 0.1762 results in a performance deterioration of less than 10%. An error of 80% still retains 50% of the performance. Note that same conclusion holds in the case of a nonzero steady-state error in the case of a leaky integrator. • [Legend Fig. 6] Relative error ρ . We observe that, as predicted by the results of Section S9 of the SI that the upper-bound is at least 0 and at most 43.75. We can also observe that in the central conic region, the performance degradation is very small. The white circle indicates the parameter values obtained in (8 ) which is clearly in the small error region (less than 1% error). This demonstrates the potential of implementing µ = µρ = 0.1813 instead of µ = µopt = 0.1762. • [Legend Fig. 7] Comparison of the proposed dynamic strategy with a static one for the parameters θ = 1, ν = µopt = 0.1762, k = 0.1, η = 0.1 and u = uopt = 0.1297. A. Evolution of the extracellular biofuel for both the dynamic and static control strategy. B. Evolution of the normalized equilibrium productivity as a function of the implementation error in ν = µopt and u = uopt . C. Evolution of the normalized 38 ACS Paragon Plus Environment

Page 38 of 47

Page 39 of 47 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


instantaneous productivity under the influence of a one-fold change in the value of βp . D. Evolution of the normalized instantaneous productivity under the appearance of leakiness term equal to uopt /2 in the production of pump proteins. • [Legend Fig. 8] Diagram describing the procedure for designing the parameters of the controller. The first step is the parameter identification step where either all the parameter are exactly identified, or their respective lower bounds are identified or simply the tolerance ρ is directly measured. Then, the set-point µ is chosen from which ν and θ can be obtained by, for instance, choosing a small ν to limit the power consumption of the controller. The last step is the stability check step where k, η need to be chosen such that the system is stable.

Tables No tables.



Figures

Figure 1: The controlled metabolic pathway. The metabolic pathway (in the cloud) consumes sugars in order to produce some butanol (not included in the model). Since butanol is toxic to the cell, we need to export it using pumps that are located on the membrane of the cells. The production of pumps is controlled by an antithetic integral controller that measures the current concentration of biofuel and acts back on the pump expression accordingly.

Figure 2: Practical implementation of the antithetic integral controller using σ-factors and anti σ-factors (denoted by σ). The reference reaction is implemented as the constitutive production of anti σ-factor molecules. This reaction can be tuned using chemical inducers or optogenetics. The production of σ-factor molecules is activated by the intracellular biofuel, possibly through the use of an intermediary biosensor. The comparison reaction is automatically realized through the complementarity property of the σ-factor and the anti σ-factor molecules. The free σ-factor molecules can then bind to the RNA polymerase to form a holoenzyme complex that specifically binds to the promoter region of the pump genes to initiate transcription.


Page 40 of 47

Page 41 of 47 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

A


B

#10-3

5 4 3 2

2

1 0 0.5

C

D

0.52

3 1 0.54

4 0.56

0.58

0.6

0.62

Figure 3: A. Evolution of the intracellular biofuel concentration bi (t) with initial condition x(0) = (0.01, 0, 0, 0, 0), k = 0.3289, θ = 1 and ν = µopt = 0.1762 (in black dots) for different values of η. When η is smaller than η ≈ 2.99, then the trajectories converge to the equilibrium value while when η is larger than η, then a stable limit cycle arises. B. Evolution of the instantaneous productivity versus the instantaneous toxicity for some piecewise constant trajectory for ν which successively takes the values 0.0145 → 0.1813 (section 1 in blue), 0.1813 → 0.3626 (section 2 in red), 0.3626 → 0.5440 (section 3 in orange) and 0.5440 → 0.6527 (section 4 in purple). We can see that after some transient phase, the trajectories always converge back to the curve described by the conservation law (10). The fact that the trajectories move in both directions illustrate the non-monotonic behavior of the productivity and toxicity as functions of µ. C. Bifurcation curve in the (k, η)-plane for θ = 1, ν = µopt = 0.1762. When the pair of parameters (k, η) is chosen to lie beneath the curve, the equilibrium point is locally asymptotically stable but comes oscillatory when the pair is located above the bifurcation curve. D. Evolution of the normalized equilibrium ∗ ∗ as a function of productivity Pe∗ (µ)/Popt and the normalized equilibrium toxicity T ∗ (µ)/Topt µ. We can clearly see that when µ = µopt = 0.1762 the productivity is maximized whereas the toxicity is minimized, which is an immediate consequence of the productivity-toxicity formula (10).


0.64


A

B

C

D

Figure 4: Response of bi (t) to a variation in the parameter αb (A), δb (B), γp (C) and βp (D) at t = 200h with initial condition x(0) = (0.01, 0, 0, 0, 0), θ = 1, ν = µopt = 0.1762, k = 0.12 and η = 0.1. The red curve corresponds to a variation of −50% of the parameter while the blue one corresponds to a variation of +100% of the parameter. In all these scenarios, the process shows perfect adaptation, thereby demonstrating the importance of integral control.


Page 42 of 47

Page 43 of 47 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


A

B

Figure 5: A. Evolution of the equilibrium productivity Pe∗ (µ) as a function of µ and for different values of ρ (the other parameters are fixed). We can see that, in accordance with the robustness result (Theorem S8.1 in the SI) and the formula (17), the value µrob maximizes the worst-case productivity. B. Evolution of the positive and negative maximum relative errors − ε+ σ and εσ in the implementation of µ = µopt as a function of the performance deterioration σ. We can see that an error of 50% in the implementation of µ = µopt = 0.1762 results in a performance deterioration of less than 10%. An error of 80% still retains 50% of the performance. Note that same conclusion holds in the case of a nonzero steady-state error in the case of a leaky integrator.



Figure 6: Relative error ρ . We observe that, as predicted by the results of Section S9 of the SI that the upper-bound is at least 0 and at most 43.75. We can also observe that in the central conic region, the performance degradation is very small. The white circle indicates the values of the parameters obtained in (8 ) which is clearly in the small error region (less than 1% error). This demonstrates the potential of implementing µ = µρ = 0.1813 instead of µ = µopt = 0.1762.


Page 44 of 47

Page 45 of 47 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


A

B

C

D

Figure 7: Comparison of the proposed dynamic strategy with a static one for the parameters θ = 1, ν = µopt = 0.1762, k = 0.1, η = 0.1 and u = uopt = 0.1297. A. Evolution of the extracellular biofuel for both the dynamic and static control strategy. B. Evolution of the normalized equilibrium productivity as a function of the implementation error in ν = µopt and u = uopt . C. Evolution of the normalized instantaneous productivity under the influence of a one-fold change in the value of βp . D. Evolution of the normalized instantaneous productivity under the appearance of leakiness term equal to uopt /2 in the production of pump proteins.



Figure 8: Diagram describing the procedure for designing the parameters of the controller. The first step is the parameter identification step where either all the parameter are exactly identified, or their respective lower bounds are identified or simply the tolerance ρ is directly measured. Then, the set-point µ is chosen from which ν and θ can be obtained by, for instance, choosing a small ν to limit the power consumption of the controller. The last step is the stability check step where k, η need to be chosen such that the system is stable.


Page 46 of 47

Page 47 of ACS 47Synthetic Biology 1 2 3 4 5

ACS Paragon Plus Environment

Perfect adaptation and optimal equilibrium productivity in a simple

Recommend Documents