Mechanistic Modeling of Biochemical Systems without A Priori Parameter Values Using the Design Space Toolbox v.3.0

Summary Mechanistic models of biochemical systems provide a rigorous description of biological phenomena. They are indispensable for making predictions and elucidating biological design principles. To date, mathematical analysis and characterization of these models encounter a bottleneck consisting of large numbers of unknown parameter values. Here, we introduce the Design Space Toolbox v.3.0 (DST3), a software implementation of the Design Space formalism enabling mechanistic modeling without requiring previous knowledge of parameter values. This is achieved by using a phenotype-centric modeling approach, in which the system is first decomposed into a series of biochemical phenotypes. Parameter values realizing phenotypes of interest are subsequently predicted. DST3 represents the most generally applicable implementation of the Design Space formalism and offers unique advantages over earlier versions. By expanding the Design Space formalism and streamlining its distribution, DST3 represents a valuable tool for elucidating biological design principles and designing novel synthetic circuits.


INTRODUCTION
Mechanistic modeling is the current gold standard for quantitative characterization and understanding of complex biological phenomena. However, mathematical analysis of rigorous mechanistic models is not yet the method of choice to address biological questions owing to a persistent gap in knowledge of the involved parameter values. Within this context, modeling without a priori values for the parameters is a radical claim. To see how radical, it helps to contrast the fundamentally new phenotype-centric approach implemented here with the conventional simulation-centric approach based on sampling and simulation (Lomnitz and Savageau, 2016a;Valderrama-Gó mez et al., 2018). Both approaches start with a typical qualitative hypothesis or model consisting of (1) the number of molecular elements, (2) their interactions with one another (connections), (3) the signs of their interactions (+/À), and (4) the stoichiometry of the interactions (number of molecules involved). Such a model would have numerous mechanistic parameters (rate constants, binding constants, equilibrium constants) whose values are unknown and difficult to determine.
In the conventional strategy, one samples parameter values, which can be done by randomly drawing values from a uniform distribution within a pre-defined range (Lee et al., 2014), simulates the non-linear differential equations, and compares the results with the corresponding experimentally observed quantitative phenotypes of the system. This is repeated numerous times to select the values that tend to improve the fit between the simulated results and the experimental data (some examples of this methodology include the work by Forger and Peskin (2003), which implemented a coordinate search algorithm to find optimal parameter values; the work by Zhou et al. (2005), which manually tuned parameter values; and the work by Mirsky et al. (2009), which used an evolutionary strategy to find a set of parameters that most closely matched experimental data). The result might be a good fit, one that is marginally acceptable, or no fit at all.
In contrast to this conventional approach, the phenotype-centric approach provides a fundamental advance that is entirely different. It starts with the same qualitative model as the conventional approach.
However, well-stablished powerful linear analysis methods (Strang, 2005) are used (Savageau, 1969(Savageau, , 1971a(Savageau, , 1971b to identify and enumerate the entire repertoire of biochemical phenotypes for the model . In this first step there is already the ability to discriminate among models/hypothesis; if the experimentally observed phenotype is not present in the repertoire, then the model can be eliminated (Lomnitz and Savageau, 2016b). On the other hand, if the phenotype of interest is present in the repertoire, then the same, powerful, well-stablished linear analytic methods can be used to predict a full set of parameter values that will realize that phenotype (Fasani and Savageau, 2010;Lomnitz and Savageau, 2016a;Valderrama-Gó mez et al., 2018).
Here, we introduce the Design Space Toolbox v.3.0 (DST3) to expand the capabilities of the Design Space formalism . This new version builds on previous iterations of the software (DST1: Fasani and Savageau, 2010, implemented the general data structures; DST2: Lomnitz and Savageau, 2016a, added parallelization routines and a graphical user interface) to allow the automatic identification and mathematical characterization of additional biochemical phenotypes arising from critically important under-determined cases.
These special phenotypes emerge from cycles, metabolic imbalances, and conservation constraints present in many biochemical systems. Examples of such systems include but are not limited to metabolic networks containing multiple reversible reactions (Schuster et al., 2000), a common motif leading to cycles; signaling cascades of chemical species with different regulatory states (e.g., achieved by different conformational or phosphorylation states) that are linked by conservation relationships (Haraldsdó ttir and Fleming, 2016;Klipp and Liebermeister, 2006); and saturable processes, for which a pathway pool is imbalanced with respect its influx and efflux, thus generating a steady accumulation or depletion of material within the system (Levy, 1999;Dahl et al., 2013).
Besides the expanded computational engine of DST3, its C-library that is now able to handle these biologically important singularities automatically when they appear individually and simultaneously, the new version of the Design Space Toolbox (DST3) has further unique features that distinguish it from previous versions. First, the portability and installation of DST3 ( Figure S1) is greatly improved and simplified by using Docker (Boettiger, 2015). Second, DST3 offers an improved and more stable IPython-based user interface for users with limited programming experience (refer to Figures S2 and S5-S11). Third, calculating the product of the global tolerances  for all parameters in log-coordinates, a proxy for a phenotype's volume in parameter space and for its associated global robustness (Valderrama-Gó mez et al., 2018) is now automated within DST3. Fourth, solvers for systems of both ordinary differential (ODE) and differential algebraic (DAE) equations (Petzold, 1983; were incorporated, thus allowing a fully integrated dynamic characterization of the Full System. The analysis of the Full System by these completely different methodologies provides an independent means of confirming the results obtained by the Design Space methodology. The Results section is divided into three sub-sections. The first sub-section reviews briefly key concepts of the Design Space formalism (see Savageau et al., 2009 andSavageau, 2010 for a more detailed theoretical treatment) needed to understand the advances described in this work. In the second sub-section, we build on these concepts to develop general mathematical strategies aimed at resolving matrix singularities arising from system topologies containing cycles, moiety conservations, and metabolic imbalances. In the third sub-section, we illustrate the capabilities of DST3 by analyzing a case study, the protocatechuate metabolic system of Acinetobacter (Dal et al., 2005;Trautwein and Gerischer, 2001), a biochemical system exhibiting multiple, nested singularities. The Methods section provides details on the software architecture of DST3, the different ways to access its computational capabilities, and the installation instructions via Docker.

Review of Key Concepts
Biochemical systems described by the power-law functions of chemical kinetics and the rational functions of biochemical kinetics can be represented by generalized mass action (GMA) kinetics (Savageau and Voit, 1987)  where a ik and b ik represent rate constants and g ijk h ijk are kinetic orders. P i and Q i are the number of positive and negative terms in the i-th equation, respectively. X i represent the variables of interest (typically the concentration of a biochemical species) in a system containing a total of n dependent and m independent variables. In general, dependent variables can be split into two groups: chemical variables, for which a differential equation exists, and auxiliary variables, for which algebraic constraints are defined. Each group contains n c and n À n c members, respectively. The set n c of chemical variables has a direct biological meaning and represents all the chemical/biological entities (i.e., enzymes, metabolites, chemical species, mRNA molecules) of a given system. On the other hand, the set n À n c of auxiliary variables has meaning in the context of recasting the system of ordinary differential equations into its GMA form. Independent variables for which a differential equation or algebraic constraint are not defined are treated as parameters.
For any system, one of the P i positive terms and one of the Q i negative terms in Equation 1 will momentarily dominate over the others in each one of the n equations in the system. This gives rise to a so-called dominant S-System (Savageau, 1969;Savageau et al., 2009), which can be generically described by Equation 2: Here, k represents indices of corresponding non-dominant terms. Steady-state concentrations of the dependent variables can be obtained in three steps (Savageau, 1969;Fasani and Savageau, 2010). By rearranging Equation 3 and taking logarithms, one obtains Equation 6: which can be written in matrix form as: where y j = ln X j , a ij = g ijpi À h ijqi , and b i = lnðb in =a iq Þ. In a second step, dependent (y D ) and independent (y I ) variables are split to obtain: The vector of dependent concentration variables y D can be obtained in a third step by matrix operations:

Strategies for Treating Three Types of Singularities
Earlier versions of the Design Space Toolbox (Fasani and Savageau, 2010;Lomnitz and Savageau, 2016b) have exclusively dealt with cases for which the inverse of the matrix A D exists ( Figure S4). However, a 11 X 1 X 2 X 3 31 12 23 31 33 . Systems Containing Cycles, Conservation Constraints, and Metabolic Imbalances (A) A cycle comprising three metabolic species is shown. Mass action kinetics is used to mathematically describe the flux through each reaction. Rate constants for each one of these reactions are shown. (B) If the fluxes represented by the blue arrows dominate over the others, a singularity arises that prevents the system from having a unique steady-state solution. This singularity is resolved by introducing a mass balance equation around a control volume (blue rectangle) containing the cycle. (C) Three-component system containing one conservation constraint. Fluxes to and from metabolite pools X 1 , X 2 , and X 3 are mathematically described using mass action kinetics. (D) No fluxes enter or leave the control volume (blue rectangle) around metabolites involved in the conservation. (E) Metabolic pool with one input and one output flux. Depending on the numerical values of fluxes v 0 and v 1 , the concentration of X 1 can steadily increase (v 0 >v 1 ), decrease (v 0 <v 1 ), or remain unchanged over time (v 0 = v 1 ).

OPEN ACCESS
4 iScience 23, 101200, June 26, 2020 iScience Article number of important system topologies-cycles, conservations, metabolic imbalances ( Figure 1)-cause the matrix A D to become singular and effectively limit the utility of earlier versions of the Design Space Toolbox. The ubiquity of such topologies in biochemical networks renders their automatic handling highly relevant in the context of the Design Space formalism.
The mathematical framework presented here to deal with these singularities can be seen as an extension of previous theory with ad hoc applications in a more limited context by our group (refer to Savageau, 1969 for metabolic imbalances; Alves and Savageau, 2003 for conservations;Tolla et al., 2015 for cycles). The extension presented in this work generalizes the theory and automates the computational applications thereby expanding the Design Space formalism

Cycles Are Resolved by Considering Global Dominance Equations
Cycles of reactions are a common feature of biochemical systems. They typically have a number of input and output fluxes. A simple example is the fumarate nitrate reduction (FNR) regulator of Escherichia coli that exists in a cycle with three forms having one influx and two effluxes (Tolla et al., 2015). The FNR global regulator is responsible for sensing the environment, more specifically the availability of O 2 . It modifies the global gene expression of the cell to adapt its machinery for the transition from aerobic to anaerobic growth. Let us consider the simple system shown in Figure 1A, in which species X 1 , X 2 , and X 3 interact to form a substrate cycle driven far from thermodynamic equilibrium. This example is deliberately selected to focus on the mathematical details of the singularity contained in the kinetic equations and on the strategy that resolves it. We start by setting up equations to describe the change in the concentration of each chemical species over time. Mass action kinetics are used to generate the rate laws describing the flux through each reaction. The resulting expressions are then combined by means of Kirchhoff's node law to generate balance equations for each metabolite in the system.
The Design Space formalism  can be applied to decompose this set of equations into different cases, each having a unique set of dominant terms and being valid within a specific region in parameter space. One such case for this system is case number 27, with case signature [22 11 21]. This signature contains three pairs of indices, one for each equation, indicating the identity of the positive and negative term dominating in each equation. Visual inspection of the matrix A D for this case reveals the presence of a linear dependency among its rows, and, thus, there is no unique steady-state solution. Nevertheless, the system of algebraic equations is consistent and a solution (or set of solutions) can be found by analyzing global dominance conditions on the influxes and effluxes that describe a mass balance around the cycle present in this system. Once the extended system has been automatically set up, the Design Space formalism can be applied to identify valid sub-cases that resolve the cyclical case. Table S1 shows S-system equations for each one of the six valid sub-cases generated from the extended system. Note the special form of these equations. Ssystems originating from the dominance analysis of the global dominance equation are used to replace the differential equation for the pool with the dominant efflux. For instance, the expression a 11 À b 11 X 1 (obtained when the first positive and first negative term in the global dominance equation are dominant) is used to replace the differential equation for X 1 (refer to sub-case 1 in Table S1). Additionally, this expression is scaled to match the coefficient of the negative term in the Full System (i.e., the original set of equations). Consider for instance the expression a 11 À 2b 33 X 3 , which is obtained when the first positive and third negative term in the global dominance equation are dominant. Since the coefficient of the negative term in the original equation is 1, the scaled expression 1 2 a 11 À b 33 X 3 is used to construct sub-case 3 of Table S1.
Owing to the special way in which terms stemming from the global dominance equation are used to construct equations for sub-cases, a three-digit case signature is introduced. This allows for tracking the origin of terms that make up the S-systems of these sub-cases. In addition to the indices of dominant positive and negative terms contained in the traditional two-digit signature, the extended three-digit signature contains the index of the equation from which its positive dominant term originated. Consider, for instance, the case signature for sub-case 27_5 in Table S1 and its associated S-system. The signature [22 312 21] dictates that differential equations for pools X 1 and X 3 are constructed by picking dominant terms in the traditional way, whereas the differential equation for pool X 2 is made from the first positive term of the third equation and the second negative term of the second equation.

Conserved Moieties Are Handled by Considering the Total Size of Conserved Pools
Pools of metabolites with constant total concentration, on some timescale, are a common feature of complex metabolic and signaling systems. They involve conserved moieties, which are groups of atoms that remain intact in all reactions of a system. AMP, NAD + , and NADP + are prominent examples of conserved moieties in energy metabolism (Haraldsdó ttir and Fleming, 2016). Consider the simple system in Figure 1C with three components linked by a conservation relationship. As in the case of the cycles in section ''Cycles Are Resolved by Considering Global Dominance Equations,'' the presence of conservations causes the matrix A D to be singular, since the three concentrations are not independent. Nevertheless, steady-state solutions can be obtained by discarding one of the differential equations and adding the algebraic constraint that the sum of the three concentrations must equal their conserved amount.
The analysis of the differential-algebraic system using the Design Space formalism involves the usual generation of cases by picking dominant terms for each of the equations of the system. The three cases that result are shown in Table S2. Note that each case is defined by only two differential equations and one algebraic constraint. In order to capture the special way in which the equations are constructed for each case, the indices of the differential equation being deleted are set to zero in the case signature. Case 3, for instance, in which the differential equation for pool X 3 is missing, has a case signature of [11 11 00 13] to reflect this fact.

Metabolic Imbalances Are Treated by Considering Knife-Edge Conditions
Flux imbalances are frequently encountered in the metabolism of engineered microbial strains (Dahl et al., 2013;George et al., 2014;Alonso-Gutierrez et al., 2017) and in inborn metabolic diseases such as phenylketonuria (Levy, 1999) and maple syrup urine disease (Haymond et al., 1973), often by the excretion of some metabolite. Let us consider the simplest example of a metabolic pool as shown in Figure 1E with one input (v 0 ) and one output (v 1 ) flux, where v 0 is a constant and the output flux is described by a Michaelis-Menten . K M and V M represent the Michaelis constant and maximal reaction rate, respectively. The change in concentration of metabolite X 1 over time can be described by the generalized mass action system where D represents an auxiliary variable introduced in the recasting process to describe the denominator of the Michaelis-Menten rate law. Let us now consider the equations for the case with signature [11 21]: which together with its associated dominance condition This system does not have a steady-state solution. Indeed, Equation 16 only provides a consistency condition for the concentration of X 1 to remain unchanged over time: 0 = v 0 À V M . We will refer to this kind of constraint as a knife-edge condition. In general, we are interested in the behavior of the system when knife-edge conditions are not satisfied, i.e., v 0 sV M . Violating the knife-edge condition in a specific direction implies an extreme value for X 1 :X 1 /N or X 1 /0. The validity of either situation is assessed by checking the validity of the associated dominance conditions, as shown in Table S3. Taken together, these results indicate that, for the system shown in Figure 1E, the concentration of the pool X 1 will steadily increase over time, i.e., it will blow up if the system's parameters fulfill the conditions K À1 M X 1 >1 and v 0 V À1 M > 1.
Note that the case with signature [11 11], and associated dominance condition, has a conventional steadystate solution given by iScience Article

Analysis of a Biochemical System Exhibiting Multiple Singularities
For a case study, we have selected the protocatechuate degradation pathway of Acinetobacter sp. strain ADP1. Protocatechuate is a chemical compound derived from phenolic species including p-cresol, 4-hydroxybenzoate, and numerous lignin monomers and is converted to b-ketoadipate by the pathway under study (Harwood and Parales, 1996). The protocatechuate degradation pathway is one of the two branches of the b-ketoadipate pathway and is widely distributed among taxonomically diverse eubacteria and fungi (Harwood and Parales, 1996). The ecological importance of the b-ketoadipate pathway resides in its key role in recycling vast amounts of aromatic material in the natural carbon cycle. From an industrial point of view, the protocatechuate degradation pathway plays a central role in the valorization of lignin for the production of fuels, chemicals, and materials (Linger et al., 2014;Salvachua et al., 2015) as well as in the detoxification of environmental pollutants (Harwood and Parales, 1996). There are many microbes capable of these industrially important processes in the context of toxic environmental hydrocarbons. Although they may have similar if not identical pathways for these functions, their genomic architectures exhibit major differences that are not well understood (Jimé nez et al., 2002;Harwood and Parales, 1996). One of the simplest of these architectures for the transport and catabolism of protocatechuate is found in Acinetobacter sp. strain ADP1. Its genetic system encodes a single polycistronic mRNA for transporter and catabolic enzymes of the protocatechuate specific pathway, as well as several shared enzymes (Trautwein and Gerischer, 2001;Dal et al., 2005). For our purposes, we shall focus only on the protocatechuate specific pathway illustrated in Figure 2A (refer to Table S5 in the Supplemental Information for a summary of symbols used in this section). The system is composed of a signaling cascade, a gene circuit, and a metabolic module. The transcription factor (PcaU) functions as both a repressor (U 1 ) and an activator (U 2 ) and has a conserved total concentration of U T = U 1 + U 2 . It controls the synthesis of a polycistronic mRNA (M 3 ) that encodes the transporter, PcaK (T 4 ), and several enzymes, three of which constitute this pathway: protocatechuate 3,4-dioxygenase, PcaGH (GH 5 ), b-carboxy-cis,cis-muconate cycloisomerase, PcaB (B 6 ), and 4carboxymuconolactone decarboxylase, PcaC (C 7 ). The environmentally supplied substrate protocatechuate (P 0 ) is transported into the cell where it becomes intracellular protocatechuate (P 8 ), which is both a metabolic intermediate and the natural inducer for the transcription factor. The following metabolites in the pathway are b-carboxy-cis,cis-muconate (CM 9 ) and g-carboxymuconolactone (CL 10 ) (Dal et al., 2005).
This system exhibits all three types of singularities discussed in section ''Strategies for Treating Three Types of Singularities,'' thus making it an excellent example for demonstrating the relevance of our work. It contains a conservation relationship from the transcription factor, cycles from the reversible enzymatic reactions within the metabolic pathway, and blow-ups from imbalances within the same pathway. The dynamics of the system can be described by a set of differential algebraic equations involving 30 parameter values, which can be considered unknown for our purposes here. Refer to Equations S34-S45 in the Supplemental Information for details on the mathematical model and to Tables S5 and S6 for an overview of the symbols used and their biological context.

Filtering the Phenotypic Repertoire for Phenotypes of Interest
Enumerating the phenotypic repertoire of a system is typically the first step in the phenotype-centric modeling strategy (Valderrama-Gó mez et al., 2018). Even systems of moderate size can exhibit a surprisingly large number of biochemical phenotypes. Therefore, the second important step is to filter the repertoire for the phenotypes of interest. For example, filtering for cases with 2 eigenvalues with positive real part can be used to identify oscillatory phenotypes (Lomnitz and Savageau, 2014), filtering for cases with 1 eigenvalue with positive real part can be used to identify multi-stability and hysteresis (Fasani and Savageau, 2013), and filtering for a logical function consisting of a pattern of dependent variables that increase, decrease, or remain unchanged in response to a change in an independent variable can be used for model discrimination (Lomnitz and Savageau, 2016a). Since all of the phenotype characteristics can be exported from DST3 to an Excel spread sheet, many types of user-defined filters can be customized to meet the user's needs (see part 3 of the tutorial contained within the DST3 Docker image for an example).
Here, we show how one can progressively filter the repertoire of the protocatechuate system to narrow the focus on phenotypes of interest. If we allow for all possibilities, DST3 shows that the system represented in Figure 2A is capable of exhibiting 3,450 valid phenotypes, of a total of 10,368 potential phenotypes (Figure S7). However, we can progressively filter this list automatically to include only those phenotypes of interest. In the context of this case study, we are interested in the steady states of this system that maximizes The signaling module responds to the inducer signal P 8 , which stimulates the conversion of the transcription factor from the repressor U 1 into the activator U 2 form. The regulator controls the synthesis of a polycistronic mRNA molecule M 3 , from which four proteins T 4 , GH 5 , B 6 , and C 7 are translated. Transporter T 4 catalyzes the import of metabolite P 8 into the cell from an external pool P 0 . Enzymes GH 5 and B 6 then catalyze the reversible conversion of P 8 into CM 9 and CM 9 into CL 10 . The last enzyme C 7 catalyzes the conversion of CL 10 into the end-product b-ketoadipate enol-lactone (not shown).
(C) A Trajectories Plot is used to characterize the operating point shown in (B) The solid gray line represents the evolution of the system starting from an initial state in which all concentrations are set to 0.001. The gray star represents the steady state reached by the system after numerical integration for 500 time units. Refer to Table S7 for further details on the numerical solvers used by DST3. The black dot next to the gray star represents the steady state predicted by DST3 for phenotype 7633 using linear algebra.
(D) The operating point of the system has been modified by decreasing k cat5 from 1.0 to 0.01 so that it is now contained within the region of the pathological phenotype 7718.1. Following this initial screening, the DST3 can be used to characterize automatically the inducible responders by comparing them on the basis of three functional criteria: global robustness to a change in phenotype, energy index (maximum flux with minimum production of protein machinery), and toxicity index (maximum flux with minimum accumulation of toxic intermediates). Global robustness is determined by the product of the global tolerances ) for all of the parameters of the system, which is a proxy for the volume of the phenotype's polytope in the system Design Space-note that the ability to automatically compute phenotypic volumes is exclusive to DST3. We define the energy index as the cost/ benefit determined by the ratio of the logarithmic gain in mRNA, which is a proxy for the increased expenditure of energy for protein production, to the logarithmic gain in the pathway flux produced, energy index = LðM 3 ;P 0 Þ=LðF;P 0 Þ. The toxicity index is the cost/benefit determined by the ratio of the logarithmic gain in the toxic intermediate, protocatechuate (P 8 ), to the logarithmic gain in the pathway flux produced, toxicity index = LðP 8 ; P 0 Þ=LðF; P 0 Þ.
The results summarized in Table 1 show that 32 of the 64 phenotypes have the best global robustness, best energy index, and best toxicity index. The next 28 phenotypes have intermediate values for these three criteria and the remaining 4 phenotypes have the worst global robustness and worst energy and toxicity indices. There is a clear trade-off revealed by this analysis. The pathway flux can be increased by moving from the phenotypes in the first group to those in the third group, but only by sacrificing global robustness, energy efficiency, and toxicity. Note that generation of the results in Table 1 did not require any kind of simulation or optimization. Rather, it involved the application of linear algebra to analyze phenotype-specific properties, such as logarithmic gains and phenotypic volumes, that were used to compare a given set of biochemical phenotypes using three different criteria: global robustness as well as energy and toxicity indices. After values for these criteria were calculated, the set of phenotypes were distributed into three different groups. Then, these groups were sorted by its total normalized volume (global robustness).

Phenotype Analysis Reveals Dynamic Properties of the System
Let us consider the phenotype numbered 7633 as a representative member of the best class, which could conceivably have been selected in nature. We start our analysis by using DST3 to predict specific values for each one of the 30 parameters required to fully define an operating point of the system within phenotype 7633 along with its phenotypic boundaries. For that, we use the Analyze Case tab ( Figure S9A) within the Main Menu of the DST3 user interface. Note that an analogous analysis can be performed for any other phenotype.
Since we have the ability to characterize the Full System exhibiting this phenotype, we can ask if there might be a strategy for further improving its performance (section ''Logarithmic Gains Can Guide the Design of Engineering Strategies'') or for avoiding dysfunction through rational engineering. An obvious dysfunction occurs when there is a violation of one of the most basic design principles, namely, the maximal velocity of a downstream enzyme should be greater than that of the upstream enzymes in the pathway (Savageau et al., The temporal behavior of metabolite pools P 8 and CM 9 is shown. CM 9 reaches a steady state of approximately 0.01, whereas the concentration of P 8 continuously increases over time and does not reach a steady state. The gray star represents the state of the system after numerical integration for 50,000 time units. When integrated for a longer period of time, the gray star continues moving further to the right at a constant CM 9 concentration, till it eventually reaches and passes the location of the black dot. Initial conditions are the same as in (C).

OPEN ACCESS
iScience 23, 101200, June 26, 2020 9 iScience Article 2009). The relation of this pathology to the phenotype 7633 is made evident in a Design Space Plot (Figure 2B) with the turnover number (k cat5 ) for the enzyme converting intracellular protocatechuate (P 8 ) to b-carboxy-cis,cis-muconate (CM 9 ) on the y axis. The relative location of physiologically relevant phenotypes coded by color is shown in Figure 2B. The operating point of the system, which was automatically calculated by DST3, is represented by the white dot within the valid region for phenotype 7633 (red polytope). In this specific view of the Design Space, the white blank region in the lower portion of the plot indicates the location of pathological phenotypes exhibiting metabolic imbalances.
The dynamical behavior of the Full System, when its steady-state operating point is located within the valid region for phenotype 7633, can be studied by simulation using a Trajectories Plot ( Figure 2C), one of three types of plots available in the Full System tab of the Main Menu. This plot shows the evolution of the system starting from a set of initial conditions to reach its steady state, which is marked by a gray star in Figure 2C. The black dot in this figure represents the steady-state prediction made by DST3 for phenotype 7633 using linear algebra. The relative position of the black dot and the gray star in the trajectories plot demonstrates the accuracy of DST3 when approximating steady states. The Full System tab also allows dynamic phenotypes (case number and signature) to be reported at each point in the solution when there is a change in dominance conditions; this is currently implemented as an exportable.xlsx file (data not shown). In the specific example of the trajectory shown in Figure 2C, the initial conditions of the system locate its operating point in phenotype 1515. As the concentrations of the pools in the system evolve toward a steady state, the operating point of the system transits through phenotypes 649 and 721 to finally reach its steady state located within phenotype 7633.
DST3 makes it possible to identify the nature of pathological phenotypes. By clicking the ''Check for Blowups'' option in the construction of the Design Space ( Figure S5), we obtain Figure 2D, which has five additional phenotypes displayed (7718.1, 7754.1, 7758.1, 9030_10.1, and 9066_10.1). Note that the ability to identify pathological phenotypes exhibiting metabolic imbalances is unique to DST3. The dynamical nature of the pathological phenotypes contained in the lower portion of the Design Space Plot can be trivially predicted. Decreasing the numerical value of k cat5 while keeping all other parameter values constant at the operating point of the system would lead to a steady increase of the toxic intermediate P 8 .  Figure S6B) and filtered for desirable phenotypes as described in the main text.

ll
OPEN ACCESS iScience 23, 101200, June 26, 2020 iScience Article A Trajectories Plot ( Figure 2E), with k cat5 decreased from its nominal value of 1 to 0.01, is used to characterize the dynamics of the Full System when the operating point of the system is located within the region of pathological phenotype 7718.1 ( Figure 2D). In Figure 2E, the state of the system after numerical integration using the IDA solver (refer to Table S7) is marked by a gray star, whereas the black dot represents the prediction made by DST3 for phenotype 7718.1 using linear algebra (note that the value of 10 12 is an arbitrarily large value used to signify an approach to infinity

Logarithmic Gains Can Guide the Design of Engineering Strategies
Once a stable and globally robust operating point for a given system has been identified, one might be interested in finding strategies to increase the flux through a specific metabolic pathway or to increase the steadystate concentration of certain intermediate metabolites.
Here, we exemplify how an analysis of logarithmic gains can be used to identify such strategies. Logarithmic gains are amplification factors relating changes in input signals (independent variables) to the resulting changes in output signals (dependent variables). The term parameter sensitivity is used instead of logarithmic gain when the effect of varying a parameter on a dependent variable is analyzed. These parameter (in)sensitivities represent the local robustness of a system, in contrast to the global robustness provided by the volume of the valid region of a phenotype in Design Space. Both logarithmic gains and parameter sensitivities are properties that depend exclusively on the kinetic orders of the system and can be calculated for concentrations or fluxes (Savageau, 1971a). DST3 allows the calculation of logarithmic gains and parameter sensitivities using the tab Analyze Case of the Main Menu in the user interface ( Figure S9A). For simplicity, we will use the term logarithmic gain for both logarithmic gains and parameter sensitivities. Table  2 lists logarithmic gains for the phenotype 7633, the representative phenotype of the first group of phenotypes with desired properties (Table 1).
A number of engineering strategies are contained in Table 2. For instance, interventions increasing the flux through the metabolic pathway without altering the steady-state concentration of the potentially toxic metabolic intermediate P 8 are identified in the third column of Table 2 and are graphically represented in Figure 3A, where each individual arrow represents a different strategy. Note that all these strategies ultimately lead to an increase in the availability of the mRNA molecule M 3 and can be categorized into two groups. The first group contains strategies that directly increase the synthesis (by either increasing a 3max , decreasing the binding constant K 2 or increasing K 1 ) and reduce the degradation (by decreasing the rate constant b 3 ) of M 3 .
The second group encompasses indirect strategies that point at increasing the steady-state concentration of the activator form of the transcription factor U 2 by modifying rate constants (P 8 or M 3 ). An analogous analysis can be done to identify strategies increasing the steady-state concentration of metabolic intermediates without increasing the pathway flux. We use the Full System tab within the Main Menu of the user interface of DST3 ( Figure S11) to demonstrate the validity of these predictions by means of two Titration Plots. In each case, the maximal synthesis rate a 3max is increased 100-fold and decreased 10-fold from its nominal operating value of 1 and the effect on the pathway flux F ( Figure 3B) and on the steady-state concentration of the toxic metabolic intermediate P 8 (Figure 3C) is computed for the Full System. As predicted by a logarithmic gain of LðF; a 3max Þ = 1 for phenotype 7633, increasing a 3max leads to an increase in the steady-state flux through the metabolic pathway in the Full System. On the other hand, and as indicated by a logarithmic gain of LðP 8 ; a 3max Þ = 0, increasing or decreasing a 3max has no effect on the steady-state concentration of P 8 in the Full System. Thus, increasing a 3max , which can be experimentally achieved by engineering the promoter region of the polycistronic mRNA or increasing the copy number of the pca operon, from its nominal operating value can be used as a strategy to increase the flux through the metabolic pathway without increasing the steady-state concentration of P 8 .

DST3 Efficiently Locates Experimentally Observed Induction Patterns of the System
Trautwein and Gerischer (2001)  iScience Article studied in this work. Based on expression experiments under different conditions, Trautwein and Gerischer suggested a bifunctional nature of the transcription factor PcaU, concluded that protocatechuate-dependent regulation of the pca genes only occurs in the presence of a functional PcaU, and quantitatively characterized the degree of induction of the pathway under different growth conditions. The authors measured a 94-fold induction of the pathway, when the pca gene expression was compared between growth on succinate and growth on p-hydroxybenzoate-note that Acinetobacter sp. strain ADP1 degrades p-hydroxybenzoate via protocatechuate, which induces pca gene expression by activating the regulator PcaU. On the other hand, succinate or its degradation products do not induce pca gene expression. This value represented the highest degree of induction observed by the authors among different growth conditions.
Trautwein and Gerischer (2001) also studied the effect that PcaU had on pca gene expression in the absence of the inducer protocatechuate. This analysis involved a wild-type and a mutant strain missing iScience Article two-thirds of the pcaU gene. After growth on succinate, pyruvate, and acetate, the DpcaU derivative strain expressed between 2-and 3.5-fold higher enzyme levels than the wild-type under the same conditions, thus suggesting that PcaU can act as a repressor in the absence of high concentrations of the inducer protocatechuate.
In section ''Filtering the Phenotypic Repertoire for Phenotypes of Interest,'' we filtered the phenotypic repertoire exhibited by the network represented in Figure 2A to identify 64 phenotypes of interest that were inducible and responded to changes in the concentration of the inducer protocatechuate. Here, we ask the question whether these phenotypes are able to exhibit induction patterns that match the experimental observations made by Trautwein and Gerischer (2001). Once again, we focus our analysis around phenotype 7633, which is a representative phenotype for the group 1 in Table 1. Figure 4A shows a Design Space plot centered on this phenotype. The x axis corresponds to the logarithm of the concentration of the environmentally supplied substrate protocatechuate P 0 , whereas the y axis corresponds to the logarithm of A B C Figure 3. Engineering Strategies for Phenotype 7633 (A) Six different strategies to increase the flux through the metabolic pathway without increasing the steady-state concentration of P 8 are shown. Each colored arrow represents an individual strategy. An arrow placed directly over a synthesis or a degradation flux targets its associated rate constant: a for synthesis and b for degradation. Arrows located within gray boxes usually target kinetic properties of an enzyme or a process. Blue upward arrows symbolize increase, whereas red downward arrows represent decrease. All of the strategies represented in the figure are biologically feasible. However, modifying kinetic properties of a given enzyme or process requires, in most of the cases, a greater experimental effort than modifying its synthesis rate. The effects of perturbing the maximal synthesis rate of the mRNA molecule M 3 ða 3max Þ on the flux through the metabolic pathway F and the steady-state concentration of P 8 are shown in (B) and (C), respectively. The operating point of the system is the one depicted in Figure 2B by the white dot. Vertical dashed lines represent the nominal value of a 3max , from which the system is perturbed.  iScience Article the total amount of the regulator PcaU present in the system ðU T Þ. Parameter values are the same as in Figure 2B.
To obtain a global view of the expression pattern exhibited by the system around phenotype 7633, we generated a Steady State Concentration Plot ( Figure 4B) of the mRNA molecule M 3 . This plot has the same axes as Figure 4A but shows the logarithm of the concentration of the mRNA molecule M 3 as a heatmap. Solid black lines in Figure 4B represent boundaries between phenotypes, whose identity can be obtained from Figure 4A. Detailed induction patterns were extracted from Figure 4B by a series of two Titration Plots at different constant levels of the regulator PcaU ( Figures 4C and 4D) and two Titration Plots at constant levels of P 0 (Figures 4E and 4F). The parametric trajectories followed by each one of these titration plots are denoted by the white dashed lines in Figure 4B.
Taken together, these four titration plots correctly capture key experimental observations made by Trautwein and Gerischer (2001). Figure 4C along with Figure 4D clearly demonstrates that the protocatechuatedependent regulation of the pca genes only occurs in the presence of a sufficiently high concentration of the regulator PcaU ðU T Þ. If the concentration of PcaU is too low, or if the regulator is not present in a functional form, increasing the concentration of protocatechuate ðP 0 Þ will not lead to a change from the basal state of the system, as shown in Figure 4D by a constant steady-state concentration of the mRNA molecule M 3 . However, if U T is increased to a sufficiently high value, a 100-fold pca gene induction by P 0 is observed ( Figure 4C). Interestingly, this extent of induction is in agreement with the maximal induction factor experimentally observed by Trautwein and Gerischer (2001).
Evidence for the repressor activity of PcaU at a low protocatechuate level is shown in Figure 4E by an increase in the steady-state concentration of the mRNA molecule M 3 as the total amount of PcaU ðU T Þ is decreased. Interestingly, a calculated 10-fold repression is in the same order of magnitude as the experimentally observed 2-to 3.5-fold repression. Note that, at a high concentration of the inducer, decreasing the total amount of PcaU ðU T Þ leads to a decrease in the steady-state concentration of the mRNA molecule M 3 , thus reinforcing the observation that PcaU acts as a repressor only at low concentrations of the inducer protocatechuate.
The ability of the underlying mathematical model to correctly describe experimentally observed induction patterns is not completely surprising, partly because the potential for these features was incorporated during the construction of the model. Rather, the information extracted from Figures

DISCUSSION
The Design Space Toolbox v.3.0 offers a variety of advantages over its predecessor version. By distributing the software via a Docker image, the installation process of DST3 is reduced to installing Docker itself. All necessary software dependencies and configurations are already contained within the Docker image, so that users can focus on the actual application of the tools to the analysis of biochemical systems. To improve software usability, DST3 comes with a more stable user interface. By integrating ODE and DAE solvers into the user interface (see Table S7 for further details), it is now possible to directly test independently the accuracy of predictions made by the Design Space formalism from within DST3. Additionally, the capabilities of the computational engine of the Design Space Toolbox were extended. DST3 is now able to analyze biochemical systems containing multiple, nested singularities, something that was out of the reach of previous versions of the toolbox. We demonstrated the utility of DST3 by analyzing a case study of an integrated biochemical system consisting of a signaling cascade, a gene circuit, and a metabolic pathway. The system's topology encoded a cycle, a conservation relationship, and the potential to exhibit blow-up behavior. We applied a recently developed phenotype-centric modeling strategy (Lomnitz and Savageau, 2015;Valderrama-Gó mez et al., 2018) to identify a stable and globally robust operating point of the system. This process involved listing the phenotypic repertoire and filtering it for phenotypes of interest. From a total of 3,450 valid phenotypes, 384 were found to be non-pathological, 192 of these physiological phenotypes were responsive to changes in the concentration of environmental substrate protocatechuate (P 0 ), and 64 of the responder phenotypes were found to be inducible. The latter fell into three groups based on the steepness of the induction characteristic: LðM 3 ; P 0 Þ = 2 with 32 phenotypes, LðM 3 ; P 0 Þ = 4 with 28, and LðM 3 ; P 0 Þ = 6 with 4. When compared on the bases of three criteria, global robustness, energy efficiency, and toxicity, the first group was best and the third group was worst. An analysis of the volume of the 64 phenotypes with desired properties revealed that their combined volume only accounted for 5.61 3 10 À5 % of the total volume of all non-pathological phenotypes identified by DST3 (see Table 1). When pathological phenotypes were considered (phenotypes exhibiting a blow-up behavior), this value decreased to 2.93 3 10 À22 %. This suggests that desirable phenotypes will have to be actively selected for by nature, since the vast majority of parameter values chosen at random (increased entropy) would produce few desirable phenotypes.
These figures highlight the power of DST3 and the phenotype-centric modeling strategy it enables. DST3 predictions regarding steady states, stability, and blow-up behavior were accurate, as demonstrated by time course, titration, and trajectory plots generated for the Full System. By finding strategies to increase the flux through the metabolic pathway of the system without increasing the steady-state concentration of an intermediate metabolite, we aimed at showing a glimpse of the potential that the Design Space formalism has to offer to the field of rational Metabolic Engineering (Bailey, 1991). Further potential applications relate to the ability of DST3 to correctly identify and characterize blow-up phenotypes, which are commonly found in metabolic systems. Often, in the process of strain development, intermediate strains are generated, in which a given intermediate metabolite excessively accumulates or is totally consumed, thus generating a metabolic imbalance within the cell. This decreases strain fitness and can ultimately lead to cellular death (Dahl et al., 2013;George et al., 2014;Alonso-Gutierrez et al., 2017). DST3 is able to identify regions in the parameter space leading to metabolic imbalances and to provide clues to rectify these phenotypes. For instance, consider the operating point of the system shown in Figure 2D, which is located within the blow-up phenotype 7718.1. Inspection of the Design Space plot around this phenotype indicates that increasing the value of k cat5 to values larger than 0.1 would place the operating point of the system within phenotype 7705, 7633, 4,177 or 5,473_6, all of which exhibit a stable, non-pathological steady state. The specific location of the operating point within any of these phenotypes, and the resulting behavior, will depend on the extent of the increase. Alternative, and more practical strategies to rectify the pathological behavior of an operating point located within phenotype 7718.1 include increasing the amount of the PcaGH enzyme by cloning its gene sequence on a controllable plasmid or engineering its ribosomal binding site.
The application of mechanistic models for the identification of metabolic engineering strategies has been rather limited. This has been mainly caused by a lack of knowledge of associated parameter values. As a consequence, constraint-based modeling has been the method of choice applied to rationally guide metabolic engineering strategies (Valderrama-Gó mez et al., 2017). By enabling a mechanistic, phenotype-ll OPEN ACCESS iScience 23, 101200, June 26, 2020 iScience Article centric modeling strategy not dependent on parameter values, the Design Space formalism and associated toolbox offers enormous potential for the field of metabolic engineering.
Elucidating biological design principles is another important area for application of the Design Space formalism that was not explored in this work owing to space limitations. In the Design Space, boundaries delimiting biochemical phenotypes are linear functions of the system's parameters in logarithmic coordinates. Thus, design principles can be readily identified in the form of mathematical inequalities involving the parameters of the system. These ideas were applied by Fasani and Savageau (2013) to study properties of toxin-antitoxin systems, which have been linked with the medically relevant persister phenotype exhibited by certain bacterial strains. The study revealed factors affecting the frequency of persisters in the population, such as the overall number of toxin-antitoxin modules and the size and position of the bistable region, a property emerging from the system's architecture.
There are many examples of systems that appear to perform the same function, and yet they exhibit radically different genomic architectures, the reasons for which are poorly understood. An example is provided by the protocatechuate degradation pathway studied in this work. It is one of the two branches of the b-ketoadipate pathway, a chromosomally encoded convergent pathway for aromatic compound degradation that is widely distributed in soil bacteria and fungi. Enzyme studies suggest that the pathway is highly conserved in diverse bacteria; however, its regulation and gene organization differ greatly (Harwood and Parales, 1996). For instance, the pathway genes from Pseudomonas aeruginosa and P. syringae are arranged in three and four different clusters, respectively. By contrast, all genes are arranged in a single cluster in Acinetobacter sp. ADP1 (studied in this work) and in P. fluorescens (Jimé nez et al., 2002). It has been suggested that evolutionary processes have shaped moldable aspects of the b-ketoadipate pathway to optimally serve diverse lifestyles of bacteria (Harwood and Parales, 1996). DST3 could be used to compare and contrast inherent aspects of each system, such as its dynamic properties, induction characteristics, and trade-offs regarding energy and toxicity, thus potentially allowing the elucidation of underlying design principles used by nature to create the alternative genomic architectures observed in organisms with different environments and lifestyles.

Limitations of the Study
Even moderate-sized systems are capable of exhibiting a large number of biochemical phenotypes. The ability to filter the repertoire for specific phenotypes is limited by currently available filtering criteria: hysteretic switches, limit cycle oscillations, global robustness, and logical patterns of system-wide increasing and decreasing responses to a given input stimulus. Although we are developing additional criteria, users may need to develop criteria to characterize phenotypes specifically of interest to them. Although the rigorous definition of boundaries between phenotypes is of primary focus in the novel phenotype-centric modeling strategy, three limitations have their origin in the mathematical transformations used to decompose the Full System into a finite set of dominant biochemical phenotypes (represented by S-Systems) and to mathematically analyze them in the logarithmic space.
(1) Numerical accuracy, which becomes undefined at boundaries between phenotypes in parameter space where there is no clear pattern of dominance. We denote the extreme case when the operating point of the system is directly located on a phenotypic boundary by a phenotype number with an asterisk. This limitation can be easily overcome by placing the operating point of the system far away from a phenotypic boundary. The Design Space Toolbox provides necessary functionalities to modify the operating point of the system accordingly (see the tutorial contained in the Docker image under /Tutorials/Tutorial_DST3).
(2) The logarithmic transformation employed to linearize S-System equations causes trivial solutions to be ignored. In such a solution, at least one of the dependent variables has a steady-state value of zero. We are currently implementing computational routines that address this limitation.
(3) The description of some processes, such as the membrane potential, require mathematical functions (e.g., exponential functions) that cannot be transformed trivially into the GMA form, which in the Design Space formalism is used to represent the Full System. An approximation by rational functions is necessary, and in most cases generating the right approximation requires some expertise. We will be addressing this topic in a future study.

Lead Contact
Further information and requests for resources should be directed to and will be fulfilled by the Lead Contact, Michael A. Savageau (masavageau@ucdavis.edu).
Tolla, D.A., Kiley, P.J., Lomnitz, J.G., and Savageau, M.A. (2015).  Figure S3. Phenotypic Deconstruction of a Biochemical System Using the DST3 Python Module. Related to Figure 1. Three steps are involved in the computational Design Space analysis of any biochemical system. First, an Equations Python object is generated by means of the class dspace.Equations. This process involves recasting ordinary differential equations describing the systems' dynamics into the GMA form, followed by a further transformation into a list of strings according to syntax rules described in the main text. Auxiliary variables stemming from the recasting process or introduced by conservation constraints need to be declared explicitly using the key argument aux_var. In a second step, the Equations object is passed to the class dspace.DesignSpace along with necessary key arguments to inform the computational engine about the presence of conservations, cycles or metabolic imbalances. The output of this second step is a DesignSpace object. Methods associated with this object allow, among other things, the generation of a list of strings containing identifiers of valid cases. In a third step, the DesignSpace object can be used to generate Case objects using valid case identifiers as input. Each Case object contains a respective SSystem object. Methods associated with these two objects allow a comprehensive characterization of each valid case, which includes, but is by no means limited to, calculation of interior parameter sets, logarithmic gains, determination of dynamical stability, etc. Refer to Lomnitz and Savageau (2016)  .remove_algebraic_constraints() .log_gain('X 1 ', 'X 2 ') .eigenvalues(pvals) .positive_roots(pvals) Figure S4. Overview of the computational workflow implemented by DST3. Related to Figure 2. Three main checkpoints define the set of computational algorithms applied to analyze a potential phenotype of a given system. Routines to handle under-determined cases (conserved relationships, metabolic imbalances and cycles) are exclusive to DST3. Refer to Section 2 of the Supplementary Information for details on the mathematical treatment.  Figure S5. Initial Input Window of the DST3 User Interface. Related to Figures 2, 3 and 4. Relevant information defining the biochemical system under study such as its mame, equations, auxiliary variables, parametric constraints and kinetic orders can be passed to the computational engine of DST3 through this window. The conservations field (highlighted by the number 1 in the blue circle), as well as the checkboxes for cycles (highlighted by the number 2 in the blue circle) and metabolic imbalances (e.g., blowing phenotypes, highlighted by the number 3 in the blue circle) can be used to control computational routines used by DST3 to analyze biochemical phenotypes.  Figure S7. System Tab of the DST3 User Interface. Related to Figure 2 and Table 1. General information of the system under analysis such as the maximal number of biochemical phenotypes it can exhibit, its signature and its equations are provided within the System tab.  Table 2. The Analyze Case tab (A) allows the generation of a comprehensive analysis of a given phenotype specified by its case number. A central functionality provided by this action tab is the automatic calculation of parameter values for the realization of the specific biochemical phenotype. This set of parameters allows a more directed exploration and analysis of the parameter space and the behavior of the system using the tabs Create Plot (Fig. S10) and Full System (Fig. S11). The Case Intersections (B) and Co-localizations (C) tabs are used to study the relative position of a given set of phenotypes in the parameter space.   Table S7 for an overview of the solvers available in DST3. Three different types of plots can be generated in this tab: Time Courses, Titrations and Trajectories plots.
Supplemental Tables   Table S1. Cyclic Cases, their Signatures and Case Numbers. Related to Figure 1. Resolving the singularity contained in case 27 (see Fig. 1B) involves a Design Space analysis of the sub-system described by Eqs. S13-S16. This analysis generates six sub-cases, each one of which is generated as dictated by a threedigit signature. Note that the reference system of equations from which dominant terms are picked according to the case signature is the Full System (Eqs. S1-S3) and not the sub-system used to resolve the singularity (Eqs. S13-S16). The case numbers for the valid sub-cases have the parent case number with an underscore followed by a number associated with the sub-case; e.g., 27_3.  Table S2. Conserved Cases and their Signatures. Related to Figure 1. The synthetic network shown in Fig.  1C can be decomposed into three different cases by applying the Design Space formalism. Dominance analysis on the conservation constraint, i.e., Eq. S25, gives rise to each of these cases. Note that each case lacks differential equations, when compared with the Full System (Eqs. S22-S25). The missing differential equation in each case is identified by a pair of zeros in the case signature.  Figure 1. Each sub-case is defined by the violation of its knife-edge condition in one of two possible directions. Subcase 2.1 is valid because fulfills the associated dominance condition . The same is not true for sub-case 2.2 because does not satisfy this dominance condition. Computationally determining the validity of each sub-case involves the solution of a linear program, as described by Fasani and Savageau (2010); the only difference being the incorporation of additional inequalities to account for the violation of the knife-edge conditions and associated extreme values for the chemical pools. Since the linear program is formulated in logarithmic coordinates, a value of or is used when or , respectively. The case numbers for the valid sub-cases have the parent case number with a period followed by a number associated with the sub-case. The total number of sub-cases corresponds to , where refers to the number of knife-edge conditions in the system.    Figure 2. The biological identity and context of symbols used to represent pools modeled in Eqs. S34-S44 is provided.

U1
Repressor form of the transcription factor PcaU PcaU belongs to the IclR family of transcriptional regulators (PobR subfamily) and carries a helix-turn-helix motif at the N terminus. It acts as an activator in the presence of the inducer protocatechuate and as a repressor in its absence (Trautwein and Gerischer, 2001)

U2
Activator form of the transcription factor PcaU

M3
Polycistronic mRNA molecule It encodes pca genes for the degradation of protocatechuate (pcaBKCHG)

4-Hydroxybenzoate transporter
This transporter protein, encoded by the gene pcaK, is also able to transport protocatechuate into the cell.

4-carboxymuconolactone decarboxylase
This protein is encodeded by the gene pcaC. It catalyzes the conversion of γcarboxymuconolactone to β-ketoadipate enol-lactone

Environmentally supplied protocatechuate
Protocatechuate is a (toxic) chemical compound derived from phenolic species including p-cresol, 4-hydroxybenzoate and numerous lignin monomers and is converted to β-ketoadipate by the pathway under study. Intracellular protocatechuate induces the transcription of pca genes via PcaU.

CM9
β-carboxy-cis,cis-muconate Intermediate metabolites of the protocatechuate degradation pathway CL10 γ-carboxymuconolactone Michaelis Menten constant of enzyme B6 for the reverse reaction Turnover number of the enzyme C7 Michaelis Menten constant of enzyme C7 Table S7. Numerical Solvers Available in DST3. Related to Figures 2, 3 and 4. The Full System tab of the DST3 user interface (Fig. S11) offers access to three different numerical solvers: a solver of ordinary differential equations (ODE) and two solvers of differential algebraic equations (DAE). When the Full System contains a conservation relationship in the form of an algebraic constraint, one of the two DAE solvers is used to numerically integrate the Full System. In this work, all numerical integrations reported were done using the IDA solver of the Sundials solver suite. The step size was 10 steps per time unit for all analyses. Default values for other relevant parameters are listed below.

Available in Python Package
Assimulo  Assimulo  Scipy

Integration method
Variable-order, variablecoefficient backward differentiation formulas in fixed-leading-coefficient form. The method order varies between 1 and 5 Runge-Kutta method of order 5 Adams methods for non-stiff problems (orders 1-12). Backward differentiation formulas (orders 1-5 with modified Newton iteration) for stiff problems

Transparent Methods
The theoretical foundation of the Design Space formalism was laid back in the 70's (Savageau, 1969;Savageau, 1971a;Savageau, 1971b;Savageau, 1979). This early work introduced the concept of S-systems and their mathematical characterization regarding dynamic stability of steady states, logarithmic gains for signal amplification, and parameter sensitivities for local robustness. Recently, these concepts were integrated to describe a generic approach to the construction of the Design Space, a structured parameter space in which qualitatively distinct biochemical phenotypes can be identified, counted and located . Linking regions of the parameter space with biochemical phenotypes has allowed the elucidation of design principles and the introduction of a radically new phenotype-centric modeling strategy (Lomnitz and Savageau, 2016;Valderrama-Gómez and Savageau, 2018 (Fasani and Savageau, 2010). Later, Jason Lomnitz introduced the Design Space Toolbox V2 (DST2) (Lomnitz and Savageau, 2016). DST2 consisted of a collection of tools comprised of a stand-alone library, written in the C language, that implements its own symbolic algebra engine and leverages open-source compiled libraries for linear algebra and linear optimization (via the GLPK library). By using multi-threaded concurrent algorithms to speed up calculations, DST2 took advantage of the parallelizable nature of the Design Space approach by analyzing each biochemical phenotype of the system independently.
Here, we show how various computational tools are integrated to create DST3. We go on to explain how Docker images and containers can be used to access DST3 on virtually any operating system. Then, we briefly describe two components of DST3 that can be used to access the computational capabilities of DST3: its user interface and its Python module. Refer to Ipythonnotebooks contained in the Docker image under /Tutorials/Tutorial_DST3 for a detailed description of the user interface.

Design Space Toolbox v.3.0
Innovations contained in DST3 aim at improving three key aspects of the software: utility, usability and portability. By further developing the C library of DST3 to allow for the automatic identification and mathematical characterization of various types of singularities, we increased the scope of systems that can be analyzed by DST3, thus improving its utility. By enhancing stability and functionality of the IPython-based user interface of DST3, we increased software usability for users with limited programming knowledge. For advanced users, we generated a python module that runs on python 3.7.3 and can be integrated into customized programs. Due to its various external software dependencies, DST2 suffered from a limited portability. We address this issue by distributing DST3 via Docker. This effectively renders Docker the only software dependency necessary to run DST3 and guarantees portability across major operating systems.

Technology overview
Three main components make up DST3: A C library, a Python package and a user interface (Fig. S1). All three components are interconnected, with the C library being the computational engine that performs most of the numerical analyses. The Python package was designed to provide high-level access to the C library, making further software development simpler and faster. The user interface was built using IPython widgets and accesses the C library through the Python module.
Four steps are required to install and access all components of DST3: That command will create a container without read and write privileges in the host computer. Files created within the container will be lost after the container is stopped. In order to grant read and write privileges, the following command should be used instead:   (Fig. S8A) and Tables menus (Fig. S8B), respectively. The System tab ( Fig. S7) contains general information about the specific set of differential equations subject to analysis. It includes the name of the system, the total number of potential cases and its system signature.
Action tabs contained in Main Menu include the following: Phenotypic Repertoire: This window allows the user to list and filter biochemical phenotypes -cases -according to their validity, case signature, log-gain values, number of eigenvalues with positive real part and volume (Fig. S6B). The resulting phenotypic repertoire can be exported to a .xlsx file or saved into the Tables menu if desired.
Analyze Case: This tab creates the full analysis of a case referenced by its case number or case signature (Fig. S9A). This analysis includes the S-system equations that mathematically define the case, conditions that need to be fulfilled for its validity, . This characterization is useful to assess the accuracy of predictions made using the S-system approximations (Fig. S11). Analysis of systems of differential algebraic equations (DAE) is performed through the python package Assimulo , while systems of ordinary differential equations (ODE) are numerically solved by the ODEINT solver of the Python package Scipy. Refer to Table S7

Conservations
As discussed before, the system described by Eqs. S22-S25. contains a conservation constraint among its constituent pools, as defined by Eq. S25. According to the DST3 syntax rules, this conservation relationship can be explicitly defined as an algebraic constraint and should be placed in the last position of the string representation of the system (see line 2 of the snippet below).
Since the Python module is being used to analyze this system, an associated auxiliary variable needs to be explicitly defined, as shown in line 3.  ['1', '2', '3'] The computational engine of DST3 is informed about the conservation constraint by using two key arguments: number_conservations=1 and resolve_conservations=True. Valid cases can be printed by means of the method valid_cases() of the ds object.

Metabolic Imbalances
The last example consists of a system containing one differential equation and one algebraic Valid cases can be customarily printed using the method valid_cases() of the ds object. Note that the case identifier 2.1 refers to the sub-case 2.1 of Table S3.

Data and Code Availability
The Docker images used by DST3 are freely available at https://hub.docker.com/r/savageau/dst3

Additional Resources
A tutorial is available as various IPython notebooks within the DST3 docker image under /Tutorials/Tutorial_DST3 X 1

Cycles are resolved by considering global dominance equations.
We start by setting up Eqs. S1 to S3 to describe the change in time for the concentration of each chemical species shown in Fig. 1A .
Necessary conditions for these terms to be dominant are described by Eqs. S7 -S11: . (S11) It should be noted that the notation used in the computational algorithms differs from the conventional notation in the main text, where parameters (b-vector) and independent variables ( vector) are separate, since both are treated here as independent variables ( vector). The bvector, which conventionally represented ratios of parameters, is now used to represent only stoichiometry values. This facilitates the computational matrix operations.
By setting the left-hand side of Eqs. S4-S6 to zero, taking logarithms and rearranging, one obtains Eq. S12, which exhibits the form of Eq. 8 in the main text: ,  (Strang, 2005), and the system of algebraic equations has the potential to contain a single or multiple cycles (Fig. S4). Its solution(s) can be found by applying the Design Space 3. Generate extended sub-system by introducing the global dominance equation(s). Eqs. S13-S16 along with the conditions defined by Eqs. S7-S11 define the extended subsystem that should be used to resolve the cyclical case generated by Eqs. S4-S6.
(S16) DST3 implements steps 1 to 3 in a recursive fashion to resolve multiple and nested cycles. Once the extended sub-system has been set up, the Design Space formalism can be applied to identify valid sub-cases resolving the cyclical case. Table S1 shows S-system equations for each one of the six valid sub-cases generated from the extended system. Note the special form of these equations. S-systems originating from a dominance analysis on the global dominance equation (Eq. S16) are used to replace the differential equation for the pool with the dominant efflux. For instance, the expression , which is obtained when the first positive and first negative term in Eq. S16 are dominant, is used to replace the differential equation for (refer to subcase 27_1 in Table S1). Additionally, this expression is scaled to match the stoichiometric coefficient of the negative term in the Full System (i.e., the original set of equations). Consider for instance the expression , which is obtained when the first positive and third negative term of Eq. S16 are dominant. Since the original stoichiometric coefficient of the negative term is 1 (see Eq. S3), the scaled expression is used to construct sub-case 27_3 of Table   S1.

Conserved moieties are handled by considering the total size of conserved pools.
Consider the simple system in Fig. 1C in the main text with three components linked by a conservation relationship. Mass balance equations can be set up for each metabolite using the rate laws shown in Fig. 1C: (S17) (S18) . (S19) Balance Eqs. S17-S19 can be compactly expressed in matrix form to yield Eq. S20: , (S20) 0 = α 11 + 2α 31 − β 11 X 1 − 2β 22 X 2 − 2β 33 X 3 α 11 − β 11 X 1 X 1 α 11 − 2β 33 X 3 1 2 α 11 − β 33 X 3 dX 1 dt = r 1 − r 2 dX 2 dt = r 2 − r 3 dX 3 dt = r 3 − r 1 dX dt = Sr with being the stoichiometric matrix and a vector of rate laws describing the flux through each reaction. The number of conservations within the system is given by: , which equals 1 for the system under consideration.
Conservation relationships can be mathematically described as linear dependencies among metabolite pools. They can be expressed in matrix form as: .
is a vector of independent variables and expresses the total pool size of each conservation. represents a vector of concentration pools. Applying Eq. S21 to the system defined by Eqs. S17-S19 yields: . DST3 is able to handle biochemical systems with single or multiple conservation relationships. For the system depicted in Fig. 1C, Eqs. S22-S25 represents an appropriate set of equations that can be analyzed by DST3 given the conservation relationship provided: The analysis of this system using the Design Space formalism involves the usual generation of cases by picking dominant terms for each of the equations. To handle the singularity generated by the conservation encoded in the system, DST3 eliminates the differential equation (s) corresponding to the dominant negative term of each conservation relationship. This generates cases with S-systems that are deficient in differential equations. Table S2 contains three cases that result from applying the Design Space formalism to the system defined by Eqs S22-S25. Note that each case is defined by only two differential equations and one algebraic constraint. In order to capture this special way of constructing the S-system equations, the indices of the differential equation being deleted are set to zero in the case signature. Case 3 for instance, S r n cr n cr = col(Null(S T )) dX 1 dt = α 11 X 3 − β 11 X 1 dX 2 dt = β 11 X 1 − β 21 X 2 dX 3 dt = β 21 X 2 − α 11 X 3 0 = CR 1 − X 1 − X 2 − X 3 n cr in which the differential equation for pool is missing, has a case signature of [11 11 00 13] to reflect this fact.
Note the similarity between the system topology of the cyclical case 27 (Fig. 1B) and the conserved system of this section (Fig. 1C). In both instances, metabolites , and interact in a cyclical fashion to introduce a linear dependency among their pools that renders their matrix singular. For cyclical cases (e.g. case 27), this dependency is eliminated by means of a global dominance equation, which represents a mass balance around the cycle. This strategy cannot be applied for conserved moieties, because they are not synthesized, degraded or exchanged with the environment (Haraldsdóttir and Fleming, 2016). From a mass balance perspective, this implies that fluxes entering or leaving the conservation do not exist, as exemplified in Fig. 1D. The singularity is thus eliminated by replacing the differential equation for one of the metabolites involved in the conservation by an algebraic constraint that contains an additional independent parameter ( ).

Metabolic Imbalances are treated by considering knife-edge conditions
Given the system shown in Fig. 1E and described by Eqs. 11-12 in the main text, consider the equations for case 2 with signature [11 21]: and its associated dominance condition: .
Substituting the algebraic constraint (Eq. S27) into the differential equation for (Eq. S26) yields the following dynamical system: This is because in a biological context, these fluxes will never be dominant fluxes defining the dynamics of the pool .
Differential equations: Algebraic constraints: dC 7 dt = α 7 M 3 − β 7 C 7 dP 8 dt = T 4 k cat 4 K M 0 −1 P 0 D 4 −1 + GH 5 k cat5 K M 5 f −1 K eq5 −1 CM 9 D 5 −1 − GH 5 k cat5 K M 5 f −1 P 8 D 5 −1 − β 8 P 8 dCM 9 dt = GH 5 k cat5 K M 5 f −1 P 8 D 5 −1 + B 6 k cat 6 K M 6 f −1 K eq5 −1 CL 10 D 6 The cooperativity (Hill number) for transcription factor binding is represented by the kinetic orders n and p. The combination of transcriptional repression and activation requires the constraints: . The and parameters are rate constants for synthesis and degradation processes, respectively. The and parameters used in the metabolic pathway represent turnover numbers and Michaelis constants, respectively. Auxiliary variables to replace the denominators in the biochemical kinetic expressions during the process of recasting the original ODE system into the equivalent GMA system. Refer to Tables S5 and S6 for the biological context of symbols used in Eqs. S46 to S62. The 17 equations (Eqs. S46-S62) contain a total of 30 parameters (Table S6). As shown in the main text, no previous knowledge (other than the constraints noted above) for any of these parameter values is required by DST3 to characterize the dynamical behavior of the system and identify a robust operating point for the system.