Understanding system dynamics of an adaptive enzyme network from globally profiled kinetic parameters

Background A major challenge in mathematical modeling of biological systems is to determine how model parameters contribute to systems dynamics. As biological processes are often complex in nature, it is desirable to address this issue using a systematic approach. Here, we propose a simple methodology that first performs an enrichment test to find patterns in the values of globally profiled kinetic parameters with which a model can produce the required system dynamics; this is then followed by a statistical test to elucidate the association between individual parameters and different parts of the system’s dynamics. Results We demonstrate our methodology on a prototype biological system of perfect adaptation dynamics, namely the chemotaxis model for Escherichia coli. Our results agreed well with those derived from experimental data and theoretical studies in the literature. Using this model system, we showed that there are motifs in kinetic parameters and that these motifs are governed by constraints of the specified system dynamics. Conclusions A systematic approach based on enrichment statistical tests has been developed to elucidate the relationships between model parameters and the roles they play in affecting system dynamics of a prototype biological network. The proposed approach is generally applicable and therefore can find wide use in systems biology modeling research.

}, , , where s is the concentration vector of A, B and C, f 1 , f 2 , and f 3 are the rate equation of A, B, and C, respectively, and p is the parameter of interest. At steady state, f 1 = f 2 = f 3 = 0, we can solve the sensitivity indices (change of concentration with respect to change of parameter) at steady state analytically as follows, p f s f dp  2  2   2  2   2   3  3  3  3  3  3  3  3  3  3   2  2  2  2  2  2  2  2  2  2   1  1  1  1  1  1

II. Additional analysis for the NFBLB model
In this section, we computed the sensitivity indices of the NFBLB model using analytical derivation (equations S2-S4) and numerical estimations of AMIGO ( Figure   S1); compared the bootstrap results of AMIGO with our kinetic motifs ( Figure S2); presented the results of additional simulations ( Figure S3), of using different PR and SN thresholds ( Figure S4) and of using additional value classes ( Figure S5); determined the number of sampling/simulations required for converging kinetic motifs ( Figure S6).   A is identical to Figure 2 in the main text. B, C, D, E, F are results of five replicates each, like A, was carried out for a total of N=10 5 sampling/simulations. G: the results of 10 6 sampling/simulations. M is the number of kinetic solutions; all other notations are the same as those described in Figure 2 in the main text. As the thresholds became more stringent, motifs may reduce their value range or even disappear (e.g. k BA , K AC in D, E; k BA , K AC and K BA in F), but the main features of the motifs and the resulting functionality networks remained unchanged from those described in the main text for SN > 1 and PN > 10 (i.e, Figure 3). All the notations are the same as those described in Figure 3 in the main text.  Thus, when the number of solutions for a specific parameter (a specific color line) is greater than that of the black line, this parameter exhibits a motif with p-value ≤ 10 -3 .
(B) The resulting kinetic motifs as a function of sampling size (number of simulations). The motifs exhibited no more changes in the last three sets of simulations (boxed by dash lines); in other words, the minimal number of sampling/simulations to obtain stabilized kinetic motifs was 8x10 4 for the NFBLP model. Note that these motifs are identical to those presented in Figure 3 in the main text.

III. An analysis of computational cost
The three steps of our method are: 1) sampling and simulation, 2) identification of kinetic motifs, and 3) functional association test. Let N be the sample size, k the number of parameters, v the number of value classes, m the number of kinetic motifs, and f the number of functional elements to be associated with parameters (e.g. 2, sensitivity (SN) and precision (PR) of adaptation, in this study). Our computation cost can be estimated as follows.
1) For the sampling/simulation step, the computation cost is N model evaluations, i.e. N simulations for a total of N simulated samples (N sets of parameters).
2) For the identification of kinetic motifs step, the computation cost is kv evaluations of the enrichment test (equation (10) in the main text), because each kinetic parameter is tested for each of the v value classes.
3) Upon identifying m kinetic motifs, each motif is tested by the Mann-Whitney U-test (see Method in the main text) for association with each of the f functional elements in both the motif and non-motif group. The computation cost for this step is estimated to be 3mf evaluations (for each motif, 2 calculations of group means for each functional element plus 1 calculation of the U-test).
The total computation cost of the three sequential steps is therefore N+kv+3mf evaluations. Since m ≤ k, by replacing m with k, the maximum computation cost is N+k(v+3f) evaluations. Because N >> k, v and f, the computational complexity of our method is essentially equivalent to the order of, and therefore dominated by, N.

IV. Simulation for a model of two linked NFBLPs (NFBLP 2 )
To investigate the computational cost of our method for a larger network, the following mock model was simulated where two NFBLB modules were linked such that the output of module 1 was the input to module 2. This resulted in a network with twice the number of kinetic parameters to be sampled than the single NFBLB model studied in the main text. For simplicity, the same 'perfect adaptation' and all its thresholds of parameters as defined in the main text were followed to find kinetic solutions and, in the subsequent statistical analysis, kinetic motifs. Note that these motifs showed considerable similarity but were not identical to those of the NFBLP model because the input to the second NFBLP module of the NFBLP 2 network was different from that to the first (or single) NFBLP module.

V. Results for GA-augmented simulations of the NFBLP model.
In this section, we present simulation results from a LHS-GA hybrid approach, in which 1000 parameter sets sampled by the LHS method [35] were subsequently augmented by 100 generations of GA optimization. A transformed SN+PR score used in the computer program of Ma et al. [33] for adaptation dynamics was used as the objective function for the GA optimization. A subroutine from GAlib (a genetic algorithm library; http://lancet.mit.edu/ga/) was implemented for this hybrid approach. For each GA generation, the best 1000 individuals (parameter sets) were selected to create a new generation by operations of genetic crossover (at 0.9 rate) and mutation (at 0.01 rate). This GA cycle was repeated until 100 generations were reached, at which point the maximal score of the objective function had reached a plateau, indicative of convergence (see below). The results (of the kinetic motifs) obtained from three independent runs of this hybrid approach are presented in Table S1. In square brackets are the intervals of parameter values for the indicated class. b.
Out of the kinetic solutions, x is the number of solutions with the value of the indicated parameter belonging to the indicated value class. c.
Out of a total of N (=10 5 ) parameter sets sampled (100 generations upon each of 1,000 initial parameter sets), y is the number of sets with the value of the indicated parameter belonging to the indicated value class. d.
An enrichment test is considered statistically significant if its p-value<10 -3 and is highlighted in boldface. Based on p-values (using 10 -3 as threshold), an enrichment state, i.e. motif, was assigned. Shaded p-values are those determined as a motif in Table 1 but not here in the GA-augmented hybrid method; note that these are the motifs exhibiting marginal significance statistically in Table 1 in the main text. The boxed p-value is the only new one found to be statistically significant (for Run 2 only) in the GA-augmented simulations.

VI. Results for the IFFLP model
The mathematical equations for the IFFLP model ( Figure S8) are: where all the notations are the same as in equation (7) in the main text. For the IFFLP model, 10 4 sets of parameter values were sampled and 6,073 of them were found to exhibit adaptation dynamics, of which 131 satisfied the criteria for perfect adaptation (see Methods). The same procedures of analysis and statistical tests as described in the main text were followed, which produced a plot of kinetic parameter distributions ( Figure S9), a kinetic-functionality network ( Figure S10), a plot of correlation between kinetic parameters (Figure S11), and results from the enrichment (Table S2) and association (Table S3) statistical tests. These results indicated that catalytic rate constant k AB was biased towards the 3 rd values class, while k AC and k BA were both biased towards the 4 th values class. The Michaelis-Menten constants K AC and K BA were biased towards the first two values classes, and K AB and K CB were biased towards one or both of the last two values classes. The other parameters, k BC , k CB and K BC , did not show any biases. In addition, the four kinetic parameters, k AC , k BA , K AC and K BA , were shown to be highly significant in improving the sensitivity scores, while k AB , K AB and K BC were significant in improving the precision scores. Again, like in the NFBLB model, a kinetic-functionality network highlighting the relationship between kinetic parameters and functionalities of dynamics can be derived ( Figure   S8).
Below, we provide possible explanations for the observed associations between kinetic motifs and perfect adaptation for the IFFLP model. Here, the introduction of an external signal decreases the concentration of activated node C by the negative regulating effect from node A ( Figure S8). For this to occur, the deactivation rate of node C must be high (the second term (v C2 ) of the rate equation dC/dt from equation (S5)), and this can be achieved with a large value for k AC , consistent with its 4 th value class motif (Table S2). According to Goldbeter and Koshland [45], in order to ensure ultra-sensitivity, the system must be in the condition of saturation, which is in accord with our finding that both parameters K AC and K BA must take on small values (i.e. < 10 0 , or value class 1 st , 2 nd , or 3 rd ).
From the topology of the network ( Figure S8), one can also observe that node A decreases directly the concentration of node C, and after some delay it also increases the concentration of node C via node B. The latter route thus serves as means for returning node C to its pre-stimulated level, and to do so the inhibition mechanism of node B via node C must be kept at a slow rate, which justifies our numerical finding that K CB takes on large values (4 th and 5 th value class) in order for the system to attain precision (as K CB appears in the denominator of the inhibition rate of node B). Note that, like in the NFBLB model, we took values within the first three value classes as small and within the last two as large, and the relative magnitudes of these parameters as deduced from analytical analysis are in accord with their corresponding value class motifs obtained from the enrichment tests (Table S2).
However, like in the NFBLB model, our numerical results do not completely agree with all our intuitive explanations. For example, k BA was found to be biased towards large values (4 th value class), which can result in a high deactivation rate for node A and should damp the effect of input signal on node C. From this, we speculate that a large k BA should cause high precision and should have no effect on the sensitivity part of the systems dynamics, but our numerical analysis showed the contrary (Table S3) values to produce perfect adaptation. This may imply a less stringent constraint on K AB than on some of the other parameters, and the observed motif of K AB represents a sufficient condition but not a necessary condition for producing adaptation dynamics.
Here, we should note that the intuitive deductions do not provide necessary proofs of the numerical simulation results. More sophisticated mathematical treatments or experimental investigations are required to resolve any discrepancies.
Finally, like in the NFBLB model, we checked whether kinetic parameters work in a cooperative way with each other when contributing to the system's adaptation dynamics. Figure S11 shows the results of correlation tests performed on all pairs of the seven kinetic parameters exhibiting value class biases. We identified two significant positive correlations (p-value < 0.05) in (k AB , K AB ) and (K AC , K CB ) pairs, and three significant negative correlations in (k AC , K AB ), (k AB , K AC ) and (k AB , K CB ) pairs.
-25 -Again, like in the NFBLB model, with the exceptions of k AC and K AC , all the correlated parameters contributed to system's precision ( Figure S10), suggesting that they work in a cooperative manner in the precision mechanism of adaptation dynamics. In contrast, kinetic parameters contributing to the system's sensitivity (i.e. k AC , k BA , K AC and K BA , see Figure S10) were not correlated with each other, implying that they function independently in the sensitivity mechanism.    b. "Pr m " ("Sn m ") is the mean logarithm value of precision (sensitivity) scores for the motif group, and "Pr~m" ("Sn~m") is the same but for the non-motif group.