Figures and data
(A). A graphical summary of the three models used in this study. The solid circles represent the metabolites and the lines represent the metabolic reactions. The dashed circles and lines are the subsystems not implemented in the model. (B-D). Example time courses for each model. The top and bottom panels are examples of the strong and weak response, respectively. (E-G). The response coefficient distribution P (χ) of each model (bottom). The average response coefficient distribution of the random catalytic reaction network (RCRN) model is overlaid (green line). The response coefficient distribution of each instance of the RCRN model is shown on the top panel. The model schemes, time courses, and distributions are aligned in each column, i.e. the leftmost network of (A), (B), and (E) are from the model developed by Khodayari et al. Khodayari et al. (2014).
The impact of setting each metabolite’s concentration to its constant, steady-state value is depicted on the metabolic network. This impact is quantified by the relative change in the average response coefficient, as detailed in (1) and (2). The filled point and filled square symbols represent the decrease 
(A). A schematic figure of the network expansion simulation. Reactions are randomly added to the original metabolic network, and then, the corresponding ODE model is constructed. In the illustration, the additional reaction is highlighted as the dashed red lines. The perturbation response simulation is performed on the expanded models to obtain the response coefficient distribution. The average response coefficient is calculated as a scalar indicator of the responsiveness for each response coefficient distribution. (B). The distribution of the average response coefficient is plotted against the number of additional reactions for the three models. The dashed line indicates the average response coefficient of the original model with no additional reaction. For each Nadd, we generated 256 extended network, and 128 trajectories were computed for each extended network.
The distribution of the average response coefficient is plotted against the number of reactions, R, for the model with N = 64 metabolites plus 3 cofactors (the panels labeled as “Cofactor”) and N = 67 metabolites without cofactors (the panel labeled as “Catalytic”). For each combination of R and f, we generated 512 random network, and computed 128 trajectories for each of these networks. The inverse temperature β is set at 10 for the cofactor version and 20 for the catalytic version. This is because the maximum chemical potential difference in the cofactor version is 2.0, but it is 1.0 in the catalytic version. The total concentrations of the cofactors are set to be unity.
The list of the substrate-level regulations in the Khodayari model
The list of the substrate-level regulations in the Boecker model
The list of the substrate-level regulations in the Chassagnole model
The list of the additional reactions used for generating Fig. S7.
Scatter plot of the response coefficient. Each point corresponds to a single initial condition, while the horizontal and vertical axes represent the response coefficient χ computed by the original and the linearized model, respectively, around the steady state. The linearized model is given by dx/dt = J x, where J is the Jacobi matrix of the model at steady state. If any metabolite concentrations become negative, the response coefficient is incomputable and “0” is assigned to such trajectories (note that the minimum value of the response coefficient is unity). Since the response coefficient is based on the logarithm of the concentrations, as the metabolite concentrations approach zero, the response coefficient becomes larger. The high response coefficient in the Boecker and Chassagnole model would be explained by this artifact. The linearized Khodayari model shows either χ ≈ 1 or χ = 0 (one or more metabolite concentrations become negative).This could be due to the number of variables in the model. For the response coefficient to have a larger value, the perturbation should be along the eigenvector that leads to oscillatory dynamics with long relaxation time (i.e., the corresponding eigenvalue has a small real part in terms of absolute value and a non-zero imaginary part). However, since the Khodayari model has about 800 variables, if perturbations are along such directions, there is a high probability that one or more metabolite concentrations will become negative.
The response coefficient distribution of each model obtained by the perturbation-response simulation with different perturbation strengths (20%, 40% and 60% relative perturbation). Each distribution is obtained by simulating 4096 trajectories. The qualitative features (double-peaks for the Khodayari model and the plateau for the Boecker and Chassagnole models) are unvaried, while the scales of the response coefficient of the Khodayari model differ among the perturbation strengths.
The time course of each model with the strongly responding metabolites being highlighted. A metabolite is highlighted if its concentration changes more than X-fold from its steady concentration. X is set to be 5 for the Boecker and Chassagnole models, while X = 25 is chosen for the Khodayari model because too many metabolites are highlighted with X = 5.
The average response coefficients of the models with the concentration of indicated metabolite fixed. The dashed line in each panel is the average response coefficient of the original model.
The Sobol’ total sensitivity index 
The response coefficient distributions of models with weakened substrate-level regulation of enzyme activity. The relaxations of the inhibition strength of AKG on CSICD and OAA on MDH increase the responsiveness of the model, while those of FUM on MDH and OAA on FRD have no effect. To relax the inhibition, the equilibrium constant for the corresponding competitive inhibition krxn_i_met is set to 106 × krxn_i_met (the inhibition term is implemented with the form xmet /krxn_i_met in the denominator of the reaction rate function), where rxn, i, and met represent the reaction name, “inhibition”, and the metabolite name, respectively. Each distribution is computed from 4096 trajectories.
The violin plot of the average response coefficient of the network expansion of the Khodayari model. The horizontal axis represents the number of additional reactions. The additional reactions are chosen from the reactions registered in the metabolic model database (see Table. S4).
The response coefficient distribution of 215 realizations of the PYK model by the MASSpy framework Haiman et al. (2021) are overlaied (each model has the same network structure, kinetic rate law, but the parameters are different). For each model realization, we have computed 256 trajectories. Trajectories had negative value for any metabolites are excluded for the computation of the distribution.
