Engineering metabolic regulation for random environments: best compromises between enzyme cost and homeostasis

Engineering metabolic regulation for random environments: best compromises between enzyme cost and homeostasis Arthur Lequertier1, Wolfram Liebermeister1, and Alberto Tonda2, 3 1Universit´e Paris-Saclay, INRAE, MaIAGE, Jouy-en-Josas, France 2UMR 518 MIA-PS, INRAE, Universit’e Paris-Saclay, Palaiseau, France 3UAR 3106 Institut des Syst`emes Complexes, Paris, France Abstract The metabolic fluxes in cells are regulated by metabolites that bind to enzymes and modulate their activities. Metabolites can stabilize their own concentrations by exerting a negative feedback on their production pathways. Although metabolic networks have been studied extensively, many regulation arrows remain unknown. To model the costs and benefits of direct enzyme regulation, we apply multi-objective optimization with one loss function describing how regulation stabilizes the metabolic state against random perturbations and another loss function scoring the extra enzyme amounts required by regulation arrows. Using the number of arrows as a third objective for ease of interpretation, we study how activating and inhibiting arrows should be arranged in the system. Using an evolutionary multi-objective approach, simulating an evolution under biological trade-offs, we explore optimal arrow configurations. Our framework can also be used with other biological objectives, including optimal adaptation or information transmission in metabolic networks, to study their trade-offs with objectives such as energy or enzyme investment in cells. Keywords: Metabolic model, Enzyme regulation, Robustness, Enzyme cost, Multi-objective optimization, NSGA-II 1 Introduction The metabolic fluxes in cells are constantly adapted to environmental changes and internal demands. An important regulation mechanism behind this is the activation or inhibition of enzymes by small effector molecules. Such regulation arrows in metabolic networks can stabilize metabolic states in changing environments. For example, some metabolites stabilize their concentrations by negatively regulating enzymes that catalyze their own producing reactions. Although the metabolic network structures of many microbial species have been uncovered, the regulation arrows in these networks are often unknown. Different species tend to share similar metabolic network structures, but appear to implement regulation in different ways. However, there are typical arrangements of regulation arrows, reminiscent of network motifs in gene regulation [1]. In contrast to these motifs, the ”motifs” of metabolic activation and inhibition arrows are additions to existing metabolic networks, for example, to existing synthesis pathways. To explain these regulation motifs, we may claim that they are evolutionarily selected for specific benefits and costs, for example, providing homeostasis, but also requiring costly investments in enzymes. Assuming such compromises, we may ask whether evolution favors ”wiring economy”, that is, whether regulation should be implemented by few, strong arrows or rather by larger numbers of weaker arrows [2]. The stabilizing and destabilizing effects of enzyme regulation have been explored in [3, 4] using computational models. Structural Kinetic Modeling [3, 5] is a method to study the dynamic effects of network structures, for example of the presence and signs of regulation arrows, regardless of quantitative parameters. This is done by considering model ensembles with a special model parameterization and randomly distributed model parameters. Here we employ this method to optimize the arrangement of regulation arrows in cells that live in stochastic environments. By adding an arrow, the cell can stabilize its metabolic states under the influence of external perturbations. Assuming that deviations from a given reference state are penalized by a quadratic fitness function, a given arrangement can be scored by the resulting fitness loss due to non-robustness, where stabilizing errors will lower this loss. However, there is also a second effect, which makes arrows costly: any regulation arrow will make its target enzyme less efficient on average, 1 and the cell needs to compensate for this by spending more enzyme, which we count as a second loss term. We postulate that the existing regulation arrows in cells can be seen as the best compromise between the two objectives. Both effects play a role for the cells’ Darwinian fitness, but they are hard to compare quantitatively, and their relative importance may vary between cell types and cell environments. Therefore, we study them here by multi-objective optimization, trying to explore different possible compromises, and to get possible explanations for the different implementations of regulation in different microbial species. In our model, we consider a cell with a given metabolic reference state. A stochastic environment perturbing the system leads to a random ensemble of metabolic states. We consider possible arrangements of activating or inhibiting regulatory arrows, each with a specific strength (a number between 0 and 1, from no effect to maximal effect), and apply multiobjective optimization to explore arrangements providing optimal compromises between favorable control properties and low enzyme cost. To study whether our two combined objectives favor parsimonious arrangements of arrows, as an example of ”wiring economy”, we also consider the number of arrows as a third minimization objective. For simplicity, we consider a simple branched metabolic pathway with six reactions and seven metabolites, where the effects of individual regulation arrows can be intuitively understood. Optimization is performed using the evolutionary multiobjective optimization algorithm NSGA-II [6]. Our optimization framework could be used more generally to study compromises between cellular objectives describing adaptation, robustness, and information transmission in metabolic networks and the burden of enzyme production and maintenance. 2 Background This section summarizes our research question as well as the methods employed: metabolic models, SKMs, and multi-objective evolutionary optimization. 2.1 Design principles for enzyme kinetics and metabolic regulation Metabolic networks contain recurrent patterns of regulation. For example, in many cases the end product of a metabolic pathway inhibits the first enzyme in the pathway. Sometimes, the same control pattern is even realized twice: once via direct enzyme regulation by the metabolite and another time, indirectly, via transcriptional control of the enzyme level. In evolution, each regulation arrow is ”reinvented” (or implemented) separately, so if typical regulation patterns exist, this hints at a selection for functionally favorable patterns. A simple design principle We hypothesize that the arrangement of the regulation arrows can be explained as the result of an optimal design, achieved by random mutations and selection in evolution. More specifically, if regulation arrows increase the performance of the metabolic system, we expect that they will be selected and preserved. We further propose that the optimality principle behind this is not based on a single objective but on a compromise between robust metabolic behavior and costs for cellular enzymes. We study whether compromises between these two objectives lead to typical known configurations such as product inhibition or substrate activation, and whether they lead to ”wiring economy”, that is, a parsimonious usage of arrows. Modeling arrangements of regulation arrows as an optimal compromise In our model, we assume that the metabolic system is exposed to a changing environment, represented by fluctuating external metabolites; their values are represented by statistical distributions, which describe a random ensemble of steady states. System performance is scored by a fitness function that trades a metabolic benefit (e.g. stable fluxes at low concentrations of intermediates) against an enzymatic cost (e.g. average concentration of enzymes). The fluctuations and the associated fitness values are computed using a local approximation around a metabolic reference state; this enables us to use analytical formulas based on metabolic control theory. We then start from a model without regulation and evaluate the fitness gained by adding regulations. To keep the different models comparable, we require that the regulated system still shows the same reference state (compensated regulation, as described below). As a consequence, each regulation arrow will increase the necessary amount of enzyme; the additional enzyme investment (to support the regulation arrows) will be treated as a second cost function. 2.2 Metabolic network models Our framework for metabolic models is based on Structural Kinetic Modeling and follows the description in [5]. A metabolic model describes the biochemical compounds and reactions within a cell. Network nodes irepresent chemical species called metabolites, while hyperedges lrepresent chemical reactions. The metabolites consumed and produced in each reaction are described by a stoichiometric matrix Nwhose rows represent the metabolites and whose columns correspond to the reac2 tions. A matrix element nji represents the stoichiometric coefficient between a metabolite and a reaction, with negative elements for reaction substrates and positive elements for reaction products. Considering the mass balance for each metabolite, we can relate the temporal variation of metabolite concentrations (in the vector c) to the reaction rate v dc dt=N·v(c),(1) with reaction rates described by unknown rate laws vi=eiκi(c). Assuming a given reference state of the system, we describe how small changes in metabolite concentrations would influence the reaction rates. The scaled elasticities of a reaction, defined as Eli = ∂ln vl/∂ ln ciand parameterized by enzyme saturation coefficients ϕli, quantify how sensitive the reaction rate vlis to small changes in the concentrations ciof metabolites. Reaction elasticities The unscaled elasticity matrix can be written as E= diag(v∗)Ediag(c∗)−1(2) where Eis the scaled version of the elasticity matrix. To study the effects of regulation arrows in a network, we write the elasticity matrix Eas the sum of two sparse matrices: E=Ekin +Ereg.(3) The matrix of kinetic elasticities Ekin, describing the effects of the enzymes’ own reactants, is assumed to be known. The matrix Ereg represents additional regulatory effects that metabolites can have on enzyme catalytic activities. If an enzyme catalyzes a reaction, the reaction rate vldepends on the activity of its enzyme as described in Section 2.2, and metabolites regulating the enzyme will lead to non-zero elasticities Ereg li . Therefore, positive (or negative) regulation coefficient Ereg li indicate activation (or inhibition) arrows in the network. The matrix Ereg will later serve as the only modified input of our model for the optimization process. Structural kinetic modeling In the Structural Kinetic Modeling (SKM) approach, we construct a linearized metabolic model in which reaction elasticities are treated as model parameters [3]. The SKM approach uses matrices Nand E, as well as fixed vectors vand cin the reference state, to represent the metabolic dynamic. The stoichiometric matrix Nrepresents the well-known topology of the metabolic network, and plausible reference states can be guessed, the elasticity matrix is typically unknown. In SKM, one may bypass this problem by sampling this vector at random [3, 7, 5]. Below, instead, we will take some of the elasticities as given and determine the others as choice variables in a multi-objective optimization. Metabolic steady-state responses to external perturbations To model variability in metabolic states, we describe all concentrations and reaction rates with random variables [8]. From given random distributions of the external variables (external metabolite concentrations and enzyme levels), we obtain the distributions of all state variables (internal metabolite concentrations and reaction rates). All model variables are described on a logarithmic scale. Assuming small variations, we replace the system dynamics by a linear approximation using concepts from Metabolic Control Analysis (MCA) [9, 10]: the linearized steady-state responses δz=Rδxof logarithmic state variable vectors to logarithmic external perturbation vectors are described by the scaled response matrix R[10], which can be computed from the elasticity matrix E. To compute the variances and covariances, we assume that all variables (on logarithmic scale) follow normal distributions centered around the reference state. For logarithmic external variables (in a random vector x), we predefine a given diagonal covariance matrix Cov(x). In our linear approximation, the logarithmic state variables (in a random vector z) will then follow a multivariate normal distribution with covariance matrix [8] Cov(z) = RCov(x)R>.(4) If all variables are considered on logarithmic scale (and collected in a vector zcomprising metabolite concentrations and fluxes), : Rz x=Rc x Rv x(5) where Rc x=−dg(c)−1L J−1NRExdg(x) Rv x= dg(v)−1(EcRc x+Ex) dg(x).(6) Here the vectors cand vrefer to state variables in the unperturbed state, the vector xcontains the corresponding external variables, ”dg” converts a vector into a diagonal matrix, the matrix J=NRE L is the Jacobian matrix of the system, the matrix Exis an elasticity matrix with respect to the external variables in x, and the matrices L and NR, with N=L NRand NRhaving full row rank, have been introduced to handle the fact that the stoichiometric matrix Nmay be rank-deficient [9]. 3 Figure 1: Metabolic network model considered in this study: a branched metabolic pathway. As an example, the figure shows two regulation arrows (red: inhibition; blue: activation). Through product inhibition and substrate activation, the arrows tend to stabilize the concentration of metabolite C. Correlated variation of metabolic variables To model random environments, we treat external perturbations, and therefore all steady-state variables, as random variables [8]. To compute the correlated variation of metabolic variables, we assume that all variables, described on a logarithmic scale, follow normal distributions centered around the reference state. For logarithmic external variables (in a random vector x), we assume a given diagonal covariance matrix Cov(x). In our linear approximation, the logarithmic state variables (in a random vector z) then follow a multivariate normal distribution with covariance matrix [8] Cov(z) = Rz xCov(x)Rz x >.(7) Eq. (7) shows how an uncertain (or variable) cell environment leads to uncertainty (or variability) in the steady state of the metabolic system. The covariance matrix Cov(z) depends on the covariance matrix Cov(x) of external variables and on the scaled response matrix R, which in turn depends on the saturation matrix, and therefore on the choice and strengths of the regulation arrows. 2.3 Multi-objective evolutionary optimization In traditional optimization methods, conflicting objectives can be addressed by constructing a composite objective function, typically as a weighted sum of individual objectives. In contrast, multiobjective optimization [11] does not require predefined weights, and instead explores the possible trade-offs among conflicting objectives. The result is not just single solution, but a set of solutions that are not Pareto-dominated by others. Evolutionary algorithms (EAs) are well suited for multiobjective optimization and are widely seen as the state of the art in this domain [12]. For this study, we chose the well-established Non-Dominated Sorting Genetic Algorithm II (NSGA-II) [6]. Despite not being the most recent method, NSGA-II remains highly competitive for optimization problems involving up to three objectives and is available across various programming languages, including C++ and Python. 3 Proposed approach To model metabolic dynamics, we start from a metabolic network with a given reference state, add some regulation arrows, and score the system by objectives describing the total enzyme demand, the resulting (non)robustness of the metabolic steady state, and the number of regulation arrows. Using multi-objective optimization, we determine arrangements of regulation arrows that represent optimal compromises between the objectives. An overview is given in Figure 2. As shown in the figure, a model instance is defined by a stoichiometric matrix N, a known reference state (v∗ and c∗), known kinetic reaction elasticities Eli, and unknown regulation elasticities Eli, parameterized by a binary vector band a vector xof real values. While bspecifies only whether the control arrow is activated or not, xrepresents the strength of the regulation and is associated with a vector of maximum and minimum elasticity values. Typically if xi= 0, it is a maximum inhibition, if xi= 0.5, it is a zero regulation, and xi= 1 is a maximum activation. A model instance yields an input-output relationship described by a response coefficient matrix R. Based on an assumed random distribution of external variables (”environment”), the model generates a covariance matrix of all the model variables. The aim is to determine the placement and strengths of arrows that provide an optimal compromise between low enzyme cost and high robustness of the metabolic steady state against perturbations. In our estimation procedure for saturation values, we, therefore, consider two biological objectives. Our loss function f1represents additional enzyme investments due to the existence of regulation arrows. Our loss function f2scores the nonrobustness of the metabolic steady state, which can be improved by the presence of regulation arrows. For convenience, the number of regulation arrows is considered as a third loss function f3in the Pareto optimization. 3.1 Regulation arrows, arrow strengths, effector elasticities, and response coefficients Introducing a new regulation arrow from metabolite ito enzyme lmeans that the rate law of enzyme 4 Figure 2: Multi-objective problem for arrangements of regulation arrows in metabolic models (description see text). lobtains a prefactor freg li (ci) = ci/KA li 1 + ci/KA li for an activation freg li (ci) = 1 1 + ci/KA li for an inhibition (8) where ciis the concentration of the effector metabolite (in mM), and KAand KIis the respective activation or inhibition constant (in mM). To account for these terms in the elasticity matrix, we add to the elasticity matrix a second matrix of scaled ”effector elasticities” Ereg li . This matrix contains a positive element for each activation arrow, a negative element for each inhibition arrow, and zero elements otherwise. The non-zero elements have absolute values in ]0,1[ (called ”arrow strengths”), which physically depend on the choices of KA li and KI li. In our model, a pattern of regulation arrows can be encoded by a sparse regulation matrix qli. In an SKM approach, where the model is parameterized by elasticities directly (and not by the underlying kinetic parameters), we can therefore sample, for each existing arrow (l,i), an arrow strength ϕli ∈ ]0,1[ (independently) and set Ereg li =qli ϕli. Adding this to the scaled ”kinetic” elasticity matrix yields the total elasticity matrix, from which the response matrices Rc xand Rv xcan be computed. 3.2 Candidate solutions In our optimality problem, each candidate solution represents an arrangement of arrows defining a regulatory elasticity matrix Ereg. For convenience, a candidate solution is described by a bit vector b (denoting the presence or absence of each possible arrow) and a real-valued vector aencoding the arrow strengths and signs (activating or inhibiting arrow) b∈ {0,1}d a∈[0,1]d,(9) where dis the number of possible arcs (counting activating and inhibiting arrows separately). The vector of binary elements bdenotes the existence of arcs between different metabolites and reactions in the pathway; the real-valued elements of aspecify the arrow signs as q= sign(a−0.5) (a < 0.5 =⇒ inhibition and a > 0.5 =⇒activation)and the arrow strengths (saturation values) as ϕ=|2a−1|), ranging between 0 (no effect) and 1 (maximal effect). The latter are directly related to kinetic constants KAor KI, as described above. Therefore, each individual can have a positive or negative arrow from a metabolite to a reaction, but not both. This choice is motivated by the fact that the opposite could lead to the generation of too many individuals with combinations of arrows canceling each other out, but with a high enzyme cost. To focus the simulations on the relevant part of the Pareto front, we capped enzyme cost function at a limit 50 (arbitrary choice). This choice prevent the algorithm from adding too many regulation arrows to improve robustness of the system, we introduce a limit to enzyme cost function of 50 (arbitrary choice) to create pressure on individuals with many arrows and thus make them mutate more intensively to reach configurations with fewer arrows. Furthermore, for calculation reasons, we 5 limit the range of possible absolute strength of a regulation to 0.99 instead of 1. Finally, by studying the sign of the higher eigenvalues of a model configuration defined by an individual, it can be tested whether the entire system is stable or not. In the case of an individual that represents an unstable model configuration, we set the 3 fitness values to their maximum to be sure that this individual will be eliminated at the next generation. 3.3 Objective functions If a new regulation arrow is added to a cell’s metabolic network, this will have two effects. First, it decreases the average activity of the regulated enzyme, and to restore the activity (and keep the reference state unchanged) the enzyme level must be increased by a factor that depends on the strength of the arrow. Hence, each arrangement of arrows and arrow strengths leads to an increase in the cell’s enzyme demand. Second, a new regulation arrow also changes the way state variables respond to external perturbations. In a random environment, each additional arrow will reshape the probability distribution of all state variables. Assuming that variability of state variables may incur fitness losses, we can characterize each arrangement of arrows by an average fitness loss. The three loss functions are computed as follows. Loss function f1: Enzyme cost The ”enzyme cost” loss of an arrow arrangement (described by vectors band a) describes the extra amount of enzyme that a cell needs to invest to maintain all steady-state fluxes and metabolite concentrations at their reference values. It is given by f1=X l""Y i 1 1−ϕli #−1#el(10) where eldenotes the reference enzyme levels and ϕli is the regulation strength of enzyme lwith respect to effector metabolite i, determined by the respective element in a. Without regulation arrows, the enzyme cost would be 0. Using arguments based on whole cell models, an increase in enzyme amounts could be translated into a decrease in cell growth rate, a proxy for Darwinian fitness [13, 14]. The formula (10) itself is easy to show. According to Eq. (8), a regulation arrow decreases the efficiency of enzyme lby a factor fli, given by 1−ϕli = 1 − |Ereg li |(this holds for activation and inhibition arrows, and also for ”nonarrows”, if we set ϕli = 0 in this case; the proof takes a few lines). For an enzyme l(with all possible regulation arrows targeting this enzyme), the necessary compensating enzyme level is, therefore, ecomp l=Y i 1 1− |Ereg li |el.(11) By subtracting the original enzyme level eland summing the costs for all enzymes, we obtain formula (10). Loss function f2: Non-robustness The ”nonrobustness” loss function penalizes variability of metabolic state variables, described by the variances from Eq. (7). At the same time, the cell is subject to an uncertain or varying environment: the external variables (e.g. external metabolite concentrations and enzyme levels) follow a random distribution (with a given mean vector and covariance matrix), leading to a random distribution of all cell variables. Therefore, also the cell fitness is not a fixed number, but a random variable. Changes in the system dynamics (e.g. a stabilization of certain cell variables by regulatory feedback loops) will change this distribution, and in particular the expectation value, that is, the average fitness. To see this, let us consider the mean vector of external variables, and the resulting mean vector of all cell variables, and let us assume that this means vector is exactly the fitness optimum. In this case, any deviation from the mean will lead to a loss of fitness and a random distribution of external variables will lead to an average loss of fitness. This fitness loss, evaluated for a set of regulation arrows that shape the system dynamics, is our second minimization objective. For the formula, we do not have to assume that the mean vector is the fitness optimum. Instead, we just need to assume our usual multivariate Gaussian distribution (with mean vector hziand covariance matrix Cov(z) and expand the fitness function quadratically around the mean: g(hzi+δz)≈f(hzi)+0+1 2Tr(Gzz Cov(z)) (12) where ”Tr” denotes the trace and zis the vector of logarithmic state variables. The first term depends only on the mean values and is not affected by the addition of regulation arrows. The second term (due to the local gradient of the fitness function) vanishes because of symmetry. Only the third term matters, and since we try to maximize fitness, we need to minimize Tr(GzzCov(z)). To see how it depends on our regulation arrows, we expand the covariance matrix (of cell variables) to Cov(z) = Rz xCov(x) (Rz x)T. We further assume that the fitness curvature matrix has strictly negative eigenvalues (which would be the case in the vicinity of a stable local fitness optimum, assuming the fitness landscape is sufficiently smooth), and set 6 M=−Gzz (now with strictly positive eigenvalues). We obtain our loss term f2= Tr[MCov(z)] (13) where the covariance matrix Cov(z) = Rz xCov(x) (Rz x)>follows from the response coefficients matrix Rz xand therefore from the arrangement of arrows and arrow strengths. In general, a model without arrows will have a non-zero loss, and arrows are expected to either improve it or make it worse. For simplicity, we may assume that biological fitness ghas a sum form g(z) = Pkgk(zk) or even depends only on a single state variable zk,g(z) = gk(zk). In these cases, the loss function simplifies, respectively, to f2=X k mkVar(zk) (14) with positive weights mk. In a network without arrows, there will be a positive loss. As arrows are added, this loss may increase or decrease. In our experiment below, we consider only this simple case, with only a single element state to be stabilized: f2= Var(zk).(15) Loss function f3: Number of regulation arrows To study wiring economy – that is, a possible parsimonious usage of arrows – we treat the number of regulation arrows as a third loss function f3. 3.4 Evolutionary operators Given the structure of a candidate solution, we designed structure-aware operators that generate valid offspring and preserve semantic information. The operators are applied with probability pmfor structure-aware mutations and pcfor structureaware crossovers. Structure-aware mutation With uniform probability, the mutation operator can change an individual by (i) flipping a bit in the vector bor (ii) performing a Gaussian mutation on a single random element of a. For a Gaussian mutation, we add a Gaussian random number to one of the values aj(choosing an element jfor which bj= 1) and bound the result within the allowed domain [0,1] ⊂R. This mutation has two hyperparameters: a mean µmand a standard deviation σm. Structure-aware crossover The crossover operator used is a one-point crossover, choosing randomly with uniform probability a cut point for both band a. This preserves the semantic information linking the existence of regulation arrows (encoded by b) with the regulation signs and strengths (encoded by a). 3.5 A test case: metabolic branchpoint model For our experiments we used the simple pathway model shown in Figure 1. Among the 7 metabolites, Aext, Eext, and Gext are external metabolites with concentrations treated as sources of variability. The concentrations of the internal metabolites B, C, D, and F and the reaction rates vl(for reaction indices l= 1..6) are state variables. As a simple objective to describe homeostasis, we use a loss function f2= Var(cC) that penalizes variability in the concentration of metabolite C. In cells, stabilizing the concentration of intermediates can be important in reducing toxicity, changes in osmotic pressure, or losses caused by membrane leakage. To focus our simulations on the most relevant part of the Pareto front, we capped the enzyme cost function f1at an upper value of 50 (arbitrary choice). For the following experiment, we assume that concentrations of metabolites involved in a reaction half saturate the reaction enzymes. Such hypotheses allow one to express the kinetic elasticity by the stoichiometric matrix, such that Ekin =−1 2NT. To introduce an asymmetry between the two branches or consuming reactions, we assume a flux of 2 in the upper consuming branch and only 1 in the lower branch, implying a flux of 3 in the producing branch. For clarity of the example, all model parameters and variables are given in arbitrary units and with values of 1 (if not stated otherwise). Rerunning our analysis with biologically realistic networks and plausible numbers would be straightforward. 4 Experimental evaluation 4.1 Set-up for multi-objective optimization As our multi-objective evolutionary algorithm for experimental evaluation, we chose the Non-Sorting Genetic Algorithm II (NSGA-II) [15], an established choice for problems with two or three objectives. For all experiments, the hyperparameters of the algorithm were set as follows: population size µ= 100, offspring size λ= 100, maximum number of generations G= 1000, pc= 0.8, pm= 0.5, µm= 0.0, σm= 0.1. In our case study, each individual (a possible arrow arrangement) was encoded by vectors band aof length d= 6 ·7 = 42. With this configuration, a single run took around 4 hours on an end-user laptop. 7 Figure 3: Pareto front for the model in Figure 1, obtained after 230 generations of evolutionary optimization with 3 objectives: enzyme cost f1(x-axis, capped at 50), non-robustness loss f2(y-axis), and arrow count number f3(shown in color, where ”6+” represents numbers ≥6). The plot represents a projection of the 3d Pareto front on the f1/f2plane. The scripts for the experiments, coded in Python, use the inspyred [16] library1for the implementation of NSGA-II. All the code and data necessary to reproduce the experiments are available in a public GitHub repository2. 4.2 Results Our optimization results are shown in Figures 3 and 5 and 4. Figure 3 shows the Pareto front of our 3-objective problem projected onto the plane of f1 (enzyme cost) and f2(non-robustness of the central metabolite C). The third objective f3, representing the number of regulation arrows, is shown in color. By separating the regulation arrows into two types, activation (blue) and inhibition (red), Figure 5 shows the distribution between individuals of the presence of this regulation. Figure 4 shows, for each regulation arrow, the Pareto front with the individual colored blue if this arrow is present among this individual as an activation, red as an inhibition, and white if the arrows are not present. The overall shape of the front reveals features that are in line with expectations, such as the location of the model instance without any regulation arrows (”empty arrangement of arrows”) on the upper left (with no extra enzyme cost at all, f1= 0 on the x-axis, and a positive non-robustness cost f2) that we used to seed the initial population. In 1Inspyred, https://github.com/aarongarrett/ inspyred 2https://github.com/albertotonda/ evolutionary-optimization-cell-models principle, adding regulation arrows could make this second cost better or worse. However, since these arrows would also increase the enzyme cost f1, all points on the front must have values of f2better than this empty arrangement; therefore, the empty arrangement forms the upper left end of the front. In the 2d projection in Figure 3, the 3d Pareto front consists of a set of separate fronts corresponding to different total numbers of arrows: for each count number f3we observe a simple curved front with a trade-off between f1and f2. If we see f1 and f2as fitness relevant (and f3as a mere feature of the solutions), the optimization does not seem to favor parsimonious solutions. At a given enzyme investment (x-axis), better homeostasis (that is, lower non-robustness) requires larger numbers of arrows. However, for count numbers larger than 6, the resulting extra benefits are negligible. Limiting the f1function to a value of 50 creates enormous pressure on individuals with a large number of arrows and a poor f1function. We can see that the decrease in enzyme costs continues well beyond this limit, because a discontinuity exists at this limit with points not aligned with the Pareto front. The predicted favorable arrows are as expected: the most common arrows to stabilize the metabolite C are feedback inhibitions from C to its producing reactions 1 and 2, or forward activations of the consumption pathways, like the arrows displayed in Figure 1. Surprisingly, there were only a few regulation arrows between the metabolite Cand the reaction 3. However, we can see in Figure 4 that the absence of this regulation was compensated by regulating reaction 3 via downstream metabolites, especially Eand F. So, if Cincreases, the concentrations of these downstream metabolites also increase, activating the inhibition of reaction flux 3, in the same way that Cwould have done if it had its own arrow for self-inhibition. The question that then arises is: Although preferential configurations for stabilizing the metabolite C exist with few regulatory arrows, why does a large number of solutions continue to have a large number of arrows? An answer may lie in the fact that these solutions are in a kind of equilibrium: as explained earlier, the introduction of an arrow can have both a beneficial and a detrimental effect on f2, and in this case, several arrows could have a compensatory effect on each other. So, when such an individual is randomly mutated, the removal of one of these arrows could have such a negative effect on f2that it would not be selected in the next generation. How can we interpret these results? Our aim was to stabilize the central metabolite in our model by minimizing our loss function f2, and to see what arrangements of arrows would achieve this at a 8 Figure 4: Occurrences of each possible regulation arrow along the Pareto front. Each subplot in the grid corresponds to one arrow (similar to Figure 5, also compare Figure 1) and contains a picture of the Pareto front. Blue (red) points denote solutions in which the arrow in question appears as an activating (inhibiting) arrow. low enzyme cost (as implemented by our loss function f1. Among the 84 possible arrows (between 7 metabolites and 6 reactions, each possible arrow being either activating or inhibiting), the intuitively expected arrows start from metabolite C itself (as the element to be stabilized) and inhibit its producing reaction or activate its consuming reactions. This is exactly what we found. Of course, this simple result is due to the simple assumptions made in our model, where variables (reference fluxes, metabolite concentrations, enzyme levels, and kinetic reaction elasticities) were set to simple numerical values and all external variables (external metabolite concentrations and enzyme levels) were treated as sources of noise. Having gained confidence in our approach, it will be interesting and easy to explore how all these factors affect the arrangements of arrows. A more general question concerns the idea of parsimonious regulation, or ”wiring economy”. Here, we found that, instead of employing only one regulation arrow, the cell should employ several arrows to reach the same degree of robustness at a lower enzyme cost (or a higher degree of robustness at the same enzyme cost). This contradicts the idea of a parsimonious usage of regulation. Of course, this result depends on the assumptions made in the model (in particular: no fixed enzyme costs for a regulation arrow but increasing costs depending on the strength of an arrow, with no cost for zero strength). Understanding the specific reasons for this will require more experiments with variants of the model, for example, variants with different sources of noise (e.g. noise only in the external substrate or only in one of the enzymes). Models with other assumptions could make different predictions, but here we can already see that a simple model with reasonable assumption does not predict wiring economy ”for free”. 5 Conclusions We developed a framework for assessing the benefit and cost of enzyme regulation in cells. Regulatory interactions were scored by loss functions that reflect biological fitness objectives: regulation incurs an enzyme cost, as higher enzyme levels are needed to compensate for reduced enzyme activities, but it can also improve the robustness of metabolic variables. The evolutionary algorithm used here should not just be seen as a computational tool, but resembles a model of biological evolution, potentially explaining diversity in metabolic regulation across microbial species. It can simulate, for example, how cell populations might evolve new regulation arrows under changing environmental selection pressures. Although some biological regulation systems show wiring economy [2] – a preference for sparse connections – our model predicts a preference for multiple, weaker arrows, at least up to a certain number. This result depends on our cost function, which penalizes large regulation strengths very strongly and makes distributed regulation more favorable – although above a certain number of arrows the extra benefits become very small. Even if our cost function is biologically realistic, this prediction is not fully conclusive: considering other costly effects in our model, giving rise to cost functions similar to a L1norm term, may have led to a different result. Moreover, our results reflect many other model assumptions made: for example, enzyme levels were treated as simple sources of noise, while in reality they are controlled by transcriptional feedback regulation, which comes at its own benefits and costs. In real cells, only a limited number of enzyme regulations are known, although many more may exist. If wiring economy is not a general principle, then weak, numerous, undetected regulatory interactions might play a larger role in metabolism than previously assumed. Parsimonious solutions might be expected if each arrow had a fixed cost, independent of the arrow strength. Such a cost function would make sense to describe allosterically regulated enzymes, which need to be larger to accommodate effector binding sites. However, here we deliberately studied a separate cost effect that would concern any kind of direct regulation (including competitive inhibition, which does not require larger enzymes): a regulation makes enzymes less efficient, requires higher enzyme levels for compensation, and thereby in9