(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 and the increase of the response coefficient, respectively. The size of the symbols is scaled using the maximum and minimum values of computed for each model. The cross symbol indicates that the original steady-state attractor becomes unstable when the concentration of the corresponding metabolite is fixed. If a dot, square, or cross is absent, the concentration of the corresponding metabolite is not a variable in the model. Only selected, representative metabolites are displayed in the figure. Comprehensive results are provided in Fig. S3.

(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 broken 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 broken 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 additional reactions used for generating Fig. S4.

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 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. S1).