Understanding system dynamics of an adaptive enzyme network from globally profiled kinetic parameters
- Austin WT Chiang^{1, 2, 3},
- Wei-Chung Liu^{4},
- Pep Charusanti^{5} and
- Ming-Jing Hwang^{1, 2, 3}Email author
https://doi.org/10.1186/1752-0509-8-4
© Chiang et al.; licensee BioMed Central Ltd. 2014
Received: 14 August 2013
Accepted: 6 January 2014
Published: 15 January 2014
Abstract
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.
Keywords
Background
Systems biology aims to unravel the design principles of living organisms from a network perspective [1, 2]. Advances in experimental studies have generated a large amount of data on several key biological processes [3–8], and networks of interactions between molecular species have been hypothesized [9–14]. Despite these advances, one unresolved challenge in systems biology is to understand how the hypothesized molecular interactions can lead to the observed biological phenomenon for complex biological systems. One way of pursuing this is via mathematical modeling of biological processes, which can also generate testable hypothesis for future experiments.
Biological processes are often complicated, and the complexity of their mathematical models usually increases with the amount of parameters involved. This generally gives rise to two fundamental problems in mathematical modeling. First, it is possible to have multiple sets of parameter values that are equally likely to produce the observed data and finding the “best” parameter set might be insufficient to fully characterize a biological system if such a parameter set is not the only set with biological relevance [15]. Second, understanding the role of individual parameters on different aspects of observed systems dynamics can be difficult for parameter-rich models as it might be too time consuming to explore this systematically and exhaustively. Although a number of approaches, for instance the genetic algorithm-based method [16] and others [15, 17], have been developed to search for parameter solutions in high dimensional spaces, they have not been extended to make inferences on the contribution of individual parameters to specific components of the system dynamics. As such, in this paper we propose a general framework for addressing the aforementioned problems in mathematical modeling of biological systems. Our methodology can be summarized as follows. For a given mathematical model of a biological process, one first defines the system’s dynamics that the model is required to reproduce; one then searches through the parameter space by a sampling method, and keeps those sets of parameters with which the model can produce the required dynamics. These then form a matrix Z of functional parameter values, with each element z _{ ij } denoting the value of the j- th parameter in the i-th parameter set. This matrix then undergoes statistical analysis to test whether a particular parameter is biased towards a certain value (or certain range of values) for the model to produce the target dynamics. After this is done for all parameters, the results can be compiled to identify recurrent parameter values and any patterns they might form. Finally, if the system’s dynamics can be decomposed into different components or parts, then further analysis can be performed to associate a parameter with a particular aspect of the dynamics.
Methods
Model of the adaptive enzyme network
The original E. coli chemotaxis model was proposed as a two-state model [29] and was later expanded to include the phosphorylation cascade [34]. In this study, the model we used is essentially the same as that used by [29] and [34], in which we consider the CheA-bound receptor complex as a single entity (node A in Figure 1C) and assume that this receptor complex exists in either a CheA active (A ^{ m }) or a CheA inactive (A) state. The superscript m denotes the methylated form of the receptor complex.
where I denotes the concentration of chemo-attractant, K _{ IA } is the ligand dissociation constant and k _{ IA } is the ligand catalytic constant.
where B ^{ P } denotes the concentration of phosphorylated demethylase CheB, K _{ BA } is the dissociation constant of phosphorylated CheB and k _{ BA } is the catalytic constant of phosphorylated CheB.
where K _{ AB } and k _{ AB } are the dissociation constant and the catalytic constant of active CheA, respectively.
where k _{ FBB } is the dephosphorylation constant.
where K _{ AC } and k _{ AC } are the dissociation constant and the catalytic constant of active CheA, respectively.
where K _{ ZC } and k _{ ZC } are the dissociation constant and the catalytic constant of CheZ, respectively.
where I denotes the input signal (i.e. the concentration of chemo-attractant); A, B and C denote the concentrations in active states; (1-A), (1-B) and (1-C) denote the concentrations in inactive forms; F _{ B } and F _{ C } are the concentration of deactivating enzymes (assumed to be a constant value of 0.5 as in [33]) that transform the active states of B and C into inactive states. Kinetic parameters k _{ ij } and K _{ ij } are respectively the catalytic rate constant and Michaelis-Menten constant for the catalytic reaction between substrate i and its regulator (activating or deactivating) enzyme j, where i and j = A, B, C, F _{ B }, or F _{ C }. Note that for each node A, B or C, the first term (v _{ 1 }) and the second term (v _{ 2 }) of the equation represent its activation and deactivation rates respectively.
Measuring adaptation
where O_{1} and O_{2} represent two steady-state values, respectively corresponding to the two input values I_{1} and I_{2} (I_{1} = 0.5 and I_{2} =0.6 following [33]), and O_{p} is the peak value of a transient pulse in response to the input change (see Figure 1A).
Sampling parameter values and numerical simulations
As in [33], Latin hypercube sampling [35] was used to sample uniformly at random the values of kinetic parameters on a logarithmic scale, with the catalytic rate constant k being in the range of [10^{-1}, 10^{1}] and the Michaelis-Menten constant K in the range of [10^{-3}, 10^{2}]. These two parameter ranges were chosen by Ma et al. [33] to, presumably, minimize computational cost while covering previously reported values used to model the E coli chemotaxis system [36–38]. In order for our results to be comparable with those from previous studies, we opted to use the same parameter ranges in this paper. For each parameter set, the model (equation (7) for the NFBLB model and equation (S1) in the supplemental data for the IFFLP model) was numerically simulated and those producing the desired adaptation dynamics were identified. The original 10^{4} parameter sets sampled by Ma et al. for the NFBLB model produced only eight kinetic solutions [33], which are insufficient for discovering parameter motifs with any statistical significance. To remedy this, a 10-fold greater sampling size was used in this study; for the IFFLP model, the original 10^{4} sampling produced 131 solutions, enough for the subsequent enrichment tests to be carried out. Results from additional sets of simulation and also from a run increasing the sample size to 10 fold, which increased the number of kinetic solutions to roughly 10 fold, indicated that the parameter motifs reported here had been reliably deduced (see Discussion).
Following Ma et al. [33], we discarded those parameter sets that render the model to produce extremely small steady-state values, persistent oscillations, weakly damped oscillations, and exceedingly long transient dynamics. For each of the remaining parameter sets (46,715 for NFBLB and 6,073 for IFFLP), the sensitivity and precision scores (i.e. equation (8) and equation (9) respectively) were calculated. We said a particular parameter set was a kinetic solution to the model if its sensitivity score was greater than 1 and its precision score greater than 10, as these criteria have been used to define perfect adaptation [33]. It can be shown from equation (8) and (9) that these thresholds were chosen to ensure a stimulus of at least 20% of the initial steady-state value, and the system can return back to this value within an error of 2%, in consistence with those used in experimental measurements (Khan et al. [39] and Alon et al. [40]). More stringent thresholds reduced the value range of the parameter motifs, but the resultant minor changes did not significantly affect the main findings (see Discussion). We used computer programs from Ma et al. [33] and Matlab software (version 7.6.0.324, release R2008a) [41] available at http://www.mathworks.com, to implement a computational pipeline to simulate both the NFBLB model and the IFFLP model numerically. We validated our simulation pipeline by reproducing the numerical results of [33] and also the steady-state solution of a four-node transcriptional regulation cascade of [42].
Enrichment test
In order to see if there exists an underlying pattern in the values of kinetic parameters for the NFBLB model to yield perfect adaptation dynamics, we obtained the kinetic solutions and plotted the distribution of parameter values for each kinetic parameter. To test whether the resulting distributions of parameter values were enriched in the sets of kinetic solutions among those 10^{5} (NFBLB) parameter sets sampled, we adopted an enrichment test that has found many uses in genomics sciences [43, 44]. For a given kinetic parameter, let N be the total number of parameter values sampled (here, N = 10^{5} for NFBLB model), and let y _{ i } be the number of those in N belonging to the i ^{th} value class in the logarithm scale (i = 1, 2, …, 5); furthermore, let M be the number of kinetic solutions (here, M = 74 for NFBLB model), and x _{ i } be the number of those in M belonging to the same i ^{th} value class. Parameter values belonging to the 1^{st} value class are within the interval [10^{-3}, 10^{-2}], those belonging to the 2^{nd} value class are within [10^{-2}, 10^{-1}], and so on and so forth, with the 5^{th} value class containing parameter values within [10^{1}, 10^{2}]. The five value classes of kinetic parameters may correspond to varying strengths of enzymatic reactions that can be measured and classified experimentally. Doubling the number of value classes resulted in a similar, albeit finer, map of parameter motifs, but did not alter the conclusions reached (see Discussion).
As in the enrichment test employed in genomics sciences [43, 44], we can then compute the p-value to measure the statistical significance of the likelihood for observing x _{ i } when the null distribution (equation (10)) is assumed to be the true count distribution. In this study, we used a p-value threshold of 10^{-3} to decide the statistical significance of the enrichment test.
Thus, for each Michaelis-Menten constant K, we carried out five independent enrichment tests, each for each value class, and for each catalytic rate constant k, we carried out two independent enrichment tests, one for the 3^{rd} value class (i.e. [10^{-1}, 10^{0}]) and the other for the 4^{th} value class (i.e. [10^{0}, 10^{1}]), due to the smaller range of parameter values for k (see above). If an enrichment test was statistically significant (p-value < 10^{-3}), a motif in the form of value class was assigned to the kinetic parameter tested. To sum up, the outcome of this analysis produced motifs of kinetic parameters, which tells us whether a particular kinetic parameter is biased towards any specific value class on the logarithm scale, or none at all.
Functional association test and parameter inter-dependence
For each kinetic parameter, values from all the kinetic parameter sets (46,715 for NFBLB and 6,073 for IFFLP) were partitioned into two groups. The first group, called the motif group, consists of parameter values belonging to the biased value class(es) as identified by the above mentioned enrichment test on kinetic solutions; and the second group, called the non-motif group, comprises those parameter values that are not in the motif group.
We then tested whether the motif group exhibited higher sensitivity or precision scores than the non-motif group by comparing the score medians of the two groups. Since the scores were not normally distributed, we used the Mann–Whitney U-Test [45] for the analysis. The test produced a z-score, and we said a particular kinetic parameter is positively associated with the sensitivity or precision parts of the adaptation dynamics if the corresponding z-score is greater than 3.29 (i.e. the upper boundary of the critical value of the 99.9% confidence intervals). In this study our focus was on finding parameters that can significantly improve the function, although some of the parameters may exhibit a large negative z-score indicative of a negative role in the function. A bipartite kinetic-functionality network can then be constructed to display the associations between the kinetic motifs and the functionalities identified.
Finally, we investigated the cooperation between every pair of kinetic parameters. Specifically, we took the kinetic parameters from those 74 kinetic solutions (131 for IFFLP) and performed Pearson correlation test between the values of a pair of kinetic parameters that appeared in the kinetic motifs identified. We said two parameters are correlated if the p-value of the correlation test is less than 0.05.
Results
Detecting kinetic motifs
Results of parameter enrichment test for the NFBLB model
Parameter | 1^{st}class | 2^{nd}class | 3^{rd}class | 4^{th}class | 5^{th}class | ||||||||||
---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|
[10^{-3}, 10^{-2}]^{ a } | [10^{-2}, 10^{-1}] | [10^{-1}, 10^{0}] | [10^{0}, 10^{1}] | [10^{1},10^{2}] | |||||||||||
x ^{ b } | y ^{ c } | p ^{ d } | x | y | p ^{ d } | x | y | p ^{ d } | x | y | p ^{ d } | x | y | p ^{ d } | |
k _{ AB } | 32 | 49,999 | 8.5E-01 | 42 | 50,001 | 1.0E-01 | |||||||||
k _{ AC } | 32 | 50,000 | 8.5E-01 | 42 | 50,000 | 1.0E-01 | |||||||||
k _{ BA } | 24 | 50,000 | 1.0E + 00 | 50 | 50,000 | 7.6E-04 | |||||||||
k _{ FBB } | 73 | 50,000 | <2.2E-16 | 1 | 50,000 | 1.0E + 00 | |||||||||
k _{ FCC } | 35 | 50,000 | 6.4E-01 | 39 | 50,000 | 2.8E-01 | |||||||||
K _{ AB } | 44 | 19,998 | 9.4E-11 | 25 | 19,999 | 1.7E-03 | 4 | 20,000 | 1.0E + 00 | 1 | 20,000 | 1.0E + 00 | 0 | 20,000 | 1.0E + 00 |
K _{ AC } | 26 | 20,000 | 7.4E-04 | 17 | 19,999 | 2.1E-01 | 17 | 20,000 | 2.1E-01 | 6 | 20,000 | 1.0E + 00 | 8 | 20,000 | 9.7E-01 |
K _{ BA } | 33 | 19,994 | 4.3E-07 | 28 | 20,000 | 1.1E-04 | 11 | 20,000 | 8.3E-01 | 2 | 20,000 | 1.0E + 00 | 0 | 20,000 | 1.0E + 00 |
K _{ FBB } | 57 | 19,997 | <2.2E-16 | 17 | 19,999 | 2.1E-01 | 0 | 19,999 | 1.0E + 00 | 0 | 20,000 | 1.0E + 00 | 0 | 20,000 | 1.0E + 00 |
K _{ FCC } | 9 | 19,999 | 9.4E-01 | 23 | 19,999 | 8.0E-03 | 32 | 20,000 | 1.5E-06 | 9 | 20,000 | 9.4E-01 | 1 | 20,000 | 1.0E + 00 |
Functional roles of kinetic motifs
For each of the seven kinetic parameters showing bias to at least one value class with statistical significance as determined in Table 1, we carried out Mann–Whitney U-test [45] to find out whether a motif (i.e. a preferred value class) of that parameter is involved in the sensitivity function of the adaptation dynamics, or, correspondingly, whether the median of the sensitivity scores for the motif group is significantly larger than that for the non-motif group (see Methods). The same analysis was then repeated for precision.
Functional association test results for the NFBLB model
Parameter | Motif [m]^{a} | Non-motif [~m] | Precision test^{b} | Sensitivity test^{b} | Function^{c} | ||||
---|---|---|---|---|---|---|---|---|---|
Pr_{m} | Pr_{~m} | z^{d} | Sn_{m} | Sn_{~m} | z^{d} | ||||
k _{ BA } | [10^{0},10^{1}] | [10^{-1},10^{0}] | 1.81 | 1.35 | 30.88 | −1.18 | −0.76 | −8.75 | PR |
k _{ FBB } | [10^{-1},10^{0}] | [10^{0},10^{1}] | 1.68 | 1.45 | 16.71 | −1.16 | −0.80 | −6.67 | PR |
K _{ AB } | [10^{-3}, 10^{-2}] | [10^{-2},10^{2}] | 1.59 | 1.42 | 6.60 | −1.40 | −0.90 | −6.65 | PR |
K _{ AC } | [10^{-3}, 10^{-2}] | [10^{-2},10^{2}] | 1.57 | 1.52 | 9.23 | 1.23 | −1.38 | 26.57 | PR,SN |
K _{ BA } | [10^{-3}, 10^{-1}] | [10^{-1},10^{2}] | 1.76 | 1.29 | 35.34 | −1.41 | −0.68 | −12.72 | PR |
K _{ FBB } | [10^{-3}, 10^{-2}] | [10^{-2},10^{2}] | 1.47 | 2.05 | −34.32 | −0.61 | −1.06 | 3.95 | SN |
K _{ FCC } | [10^{-1}, 10^{0}] | [10^{-3}, 10^{-1}] and [10^{0}, 10^{2}] | 1.63 | 1.34 | 16.98 | 3.30 | −2.31 | 28.37 | PR,SN |
Cooperation of kinetic parameters
Experimental data for kinetic parameters of E. coli chemotaxis
Experimental data for the kinetic parameters of E. coli chemotaxis
Description | Reaction^{a} | Eq7^{b} | k(s^{-1}) | K | Notes^{c}[Reference] |
---|---|---|---|---|---|
Receptor Complex Demethylation | ${B}^{P}+{A}^{m}\stackrel{{K}_{\mathit{BA}}}{\leftrightarrow}{B}^{P}{A}^{m}\stackrel{{k}_{\mathit{BA}}}{\to}{B}^{P}+A$ | v _{ A2 } | 1.2 | 0.08 | k ^{ exp } = 1.2 s^{-1}[46] and K ^{ exp } = 0.39 μ M [46]; k _{ BA } = k ^{ exp } and K _{ BA } = K ^{ exp }/[A_{t}]. |
CheB | ${A}^{m}+B\stackrel{{K}_{\mathit{AB}}}{\leftrightarrow}{A}^{m}B\stackrel{{k}_{\mathit{AB}}}{\to}{A}^{m}+{B}^{P}$ | v _{ B1 } | 3.2 | 0.281 | k ^{ exp } = 3.2 s^{-1}[37] and K ^{ exp } = 1.405 μ M [37]; k _{ AB } = k ^{ exp } and K _{ AB } = K ^{ exp }/[B_{t}]. |
Phosphotransfer | |||||
CheB | ${B}^{P}\stackrel{{k}_{{F}_{B}B}}{\to}B$ | v _{ B2 } | 0.35 | - | k ^{ exp } = 0.7 s^{-1} at 35°C, or 0.35 s^{-1} at 25°C [47]; k _{ FBB } = k ^{ exp } . |
Dephosphotransfer | |||||
CheY | ${A}^{m}+C\stackrel{{K}_{\mathit{AC}}}{\leftrightarrow}{A}^{m}C\stackrel{{k}_{\mathit{AC}}}{\to}{A}^{m}+{C}^{P}$ | v _{ C1 } | 650 | 0.36 | k ^{ exp } = 650 s^{-1}[48] and K ^{ exp } = 6.5 μ M [48]; k _{ AC } = k ^{ exp } and K _{ AC } = K ^{ exp }/[C_{t}]. |
Phosphotransfer | |||||
CheY | ${F}_{C}+{C}^{P}\stackrel{{K}_{{F}_{C}C}}{\leftrightarrow}{F}_{C}{C}^{P}\stackrel{{k}_{{F}_{C}C}}{\to}{F}_{C}+C$ | v _{ C2 } | 30 | 0.006 | k ^{ exp } = 650 s^{-1}[34] and K ^{ exp } = 0.1 μ M [34]; k _{ FCC } = k ^{ exp } and K _{ FCC } = K ^{ exp }/[C_{t}]. |
Dephosphotransfer |
and the concentration of node C (i.e. CheY) disappears from the rate equation completely.
The kinetic parameter K _{ FBB } is involved in the deactivation of node B in the NFBLB model and K _{ FBB } is also found to be biased to small values (1^{st} value class motif in Table 1) for the model to exhibit sensitivity dynamics. Although node B of the NBFLP model corresponds to CheB in the chemotaxis of E. coli, to the best of our knowledge there are no natural enzymes reported to deactivate CheB^{p} in the manner suggested by the NFBLB model. Perhaps an unrevealed enzyme exists to dephosphorylate CheB^{p}, the phosphorylate form of CheB, or there might be other mechanism of CheB^{p} dephosphorylation that is more complicated than the current knowledge can offer.
As for precision, the empirical estimates for K _{ AB } (normalized with the concentration of CheB), k _{ FBB }, K _{ BA } (normalized with the concentration of receptor complex) and k _{ BA } are 0.281, 0.35, 0.08 and 1.2, respectively (Table 3). These are also in agreement with our findings (Table 1) that K _{ AB }, k _{ FBB } and K _{ BA } are biased towards small values, and k _{ BA } is constrained to large values. K _{ AB } and k _{ FBB } are involved in the phosphorylation and dephosphorylation of CheB, respectively. According to Ma et al. [33], K _{ AB } is constrained to small values by the topological features of the NFBLB model such that the rate equation of node B (i.e. CheB) is independent of the input level I; this then implies the system is in a stable state independent of the initial perturbation and thus is able to maintain high precision level. From equation (7), a small value for k _{ FBB } can reduce the dephosphorylation rate of CheB, this in turn increases the deactivation rate of the receptor complex, thereby ensuing high precision. Finally, K _{ BA } and k _{ BA } play a part in the demethylation of the receptor complex, and intuitively the rate of its demethylation must be great (e.g. with a large k _{ BA } and a small K _{ BA }) such that the perturbation initiated by the input signal can be mitigated in order to maintain high precision.
These observations of key parameters of the NFBLB model are furthermore supported by a number of experimental findings on the chemotaxis system of E. coli: 1) K _{ FCC } having a large effect on the sensitivity part of the system dynamics is in agreement with CheZ playing an important role in adjusting the concentration of CheY [43]; 2) K _{ AC } being an important kinetic parameter in affecting sensitivity agrees well with the receptor complex being a major contributor to sensitivity [44]; 3) CheB mutants being far less sensitive than the wild type due to the functional abnormality of the receptor complex in those mutants [45] implies that K _{ FBB } is an important parameter affecting sensitivity; 4) that phosphorylated CheB can increase the rate of receptor demethylation and thus speeds up adaptation [51] supports the finding that k _{ BA }, k _{ FBB }, K _{ AB } and K _{ BA } were all important parameters in contributing to the system’s precision mechanism.
Discussion
We have developed here a methodology using parameter profiling and enrichment statistical tests to uncover not only the sets of kinetic parameters with which a model produces user-specified system dynamics, but also motifs (i.e. enriched value classes) of these parameters and their associations with specific functional aspects of the system’s dynamics. For these tasks, conventional methods usually focus on identifying the “best” parameter set in fitting the empirical data and on using complicated analytical explorations of models (e.g., sensitivity analysis) [52–54] or by a laborious local approach examining how changing the parameters one at a time would affect systems dynamics [55–57]. Note that the sensitivity analysis of conventional methods (not to confuse with the sensitivity of adaptation (Equation 8) studied in this work) is usually carried out for one specific output dynamics, whereas our profiling approach investigated the sensitivity of functional elements (sensitivity and precision of adaptation) for a collection of outpout dynamics (those that qualified as a perfect adaptation; there were 74 for the NFBLB model). Interestingly, the kinetic parameters showing value class biases (Table 2) are those exhibiting non-negligible sensitivity indices obtained from analytical derivation (Additional file 1: Equation S2) or from the numerical method implemented in AMIGO [58] (Additional file 1: Figure S1), the latter two being computed using the output dynamics and the parameter set of a randomly chosen one of the 74 kinetic solutions. Note that bootstrap-derived distributions of the kinetic parameters and their confidence intervals for this particular kinetic solution were generally in accord with the kinetic motifs deduced from the 74 kinetic solutions (Additional file 1: Figure S2). In summary, by profiling the parameters as a whole our method takes a global view to find not just one but clusters of viable parameter sets, thus moving a step further to account for the complexity of biological systems. Although a number of global-view approaches have recently been developed to sample from large and high dimensional parameter spaces, including a combined global and local exploration [15] and an approach with model checking on partitioned regions of parameters [17], these studies do not make inferences on whether motifs of parameters exist and how they might contribute to specific elements of the system dynamics. Our approach therefore offers a simple framework to systematically characterize and elucidate the functional contribution of kinetic parameters in a biological network.
The kinetic motifs obtained are quite robust in that only minor differences in their resolution were resulted from independent sampling runs (Additional file 1: Figure S3), different thresholds of the adaptation objectives (Additional file 1: Figure S4), and added number of value classes (Additional file 1: Figure S5). Further analysis showed that, for NFBLB and to satisfy the statistical significance of p-value < 10^{-3}, 8 × 10^{4} (rather than 10^{5}) sampling runs were required to converge and stabilize the kinetic motifs (Additional file 1: Figure S6). The computational cost of our method is dominated by the sampling and simulation step (see an analysis Additional file 1: Section III). To investigate the difficulties that will inevitably arise from larger networks for our method, we artificially linked two modules of NFBLB together (see Additional file 1: section IV) while requiring the system to produce the same adaptation dynamics as before. The results showed that 5 times of sampling/simulations (4 × 10^{5} vs. 8 × 10^{4}) were needed to produce stabilized kinetic motifs (Additional file 1: Figure S7) for the twice-sized network (20 kinetic parameters vs. 10). However, it should be noted that the required number of sampling/simulations depends on many factors, including the specified output dynamics, the level of statistical significance desired, and the network topology (e.g., the IFFLP model needed 10 times less number of sampling/simulations than the NFBLB model to exhibit stabilized kinetic motifs, despite they both having 10 kinetic parameters). Note that our approach is quite general and can be integrated with other approaches. For instance, the Latin hypercube sampling (LHS) used here (and as in [33]) can be replaced by other methods, such as those mentioned above, to identify the kinetic solutions needed in the subsequent enrichment tests. Results from an experiment of combining LHS and genetic algorithm (GA) showed that 10^{3} of LHS sampling followed by 100 generations of GA optimization yielded a similar set of kinetic motifs for NFBLB (Additional file 1: Table S1), suggesting that optimization can help to find more solutions from a smaller initial parameter set, but in this case the hybrid approach did not reduce the computational cost (since it also needed a total of 10^{5} model evaluations), and could miss some of the marginally significant motifs (Additional file 1: Table S1). While further research is required to fully address the dimensionality problem of scaling up the system, complex biological networks are known to be composed of simple recurrent structural components [59–61]; a deeper understanding of these network components is an important first step toward a better understanding of the assemblage and functioning of a much larger system.
One interesting observation from the case studies of the NFBLB and the IFFLP (presented in Additional file 1: section VI) adaptation model is that the majority of parameters seem to contribute to only one single functionality (i.e. either sensitivity or precision, see Figure 3, Additional file 1: Figure S10). If the mathematical models are indeed mechanistic, this could have important biological implications as follows. The behavior or dynamics of a biological system is likely to be the result of intricate interactions of many biological processes. If a biological process has a major influence on all aspects of the systems behavior, then any changes to such a process may have a drastic impact on the system. Modularity is ubiquitous in biological systems as it provides an effective mechanism to corral damaging perturbations to local consequences [62]. Thus, to ensure system robustness, evolution might have favored a biological system with a fine division of labors (i.e. modularity) among different biological components or processes [63–65]. Here, in the kinetic-functionality association, we may have uncovered yet another example of nature’s modularity design manifestated in the organization of kinetic parameters.
Previous studies have shown that certain types and arrangements of network structures are required to produce certain types of system dynamics [33, 66–68]; here, we have shown that, for a given network structure, certain types of values, or motifs, also exist for kinetic parameters in order to achieve specific system dynamics. Our results suggest that, as has been noted by others [69, 70], system dynamics can put constraints on the values of kinetic parameters. The discovery of these motifs underscores the intricate inter-relationships between structure (i.e. topology) of the biological network, kinetic parameters of the reactions involved, and the function of the biological system. Delineation of these relationships by methods such as the one developed here, which is general and can be applied to other types of well defined dynamics, will greatly advance our understanding of the design principles of prototype biological systems.
Conclusions
An increasing number of studies have revealed that complicated biological systems often share simple and universal design principles that are more understandable to biologists. The identification of motifs in biological networks is a prime example relating recurrent network structures to biological functions. Many studies have also argued for the importance of kinetic parameters in determining the dynamics of biological networks, but dissecting the association between system dynamics and kinetic parameters has been difficult. In this study, we have developed a methodology, akin to the enrichment analysis of gene expression profiles, to determine whether a preference of kinetic parameters adopting certain parameter values exists. Such preferences, or kinetic motifs, encapsulate the possible roles and functional constraints of kinetic parameters. Our analysis on models for the adaptation dynamics of the chemotaxis of Escherichia coli showed that design principles also exist from the perspective of kinetic parameters. Our methodology, owning to its generality and simplicity, provides a computational framework for understanding the kinetic mechanics of prototype biological networks.
Declarations
Acknowledgements
We thank Dr. Chao Tang of the Department of Bioengineering and Therapeutic Sciences, UCSF, and Dr. Wenzhe Ma of the Department of Systems Biology, Harvard Medical School, for making their computer programs available to us. AWTC is grateful to Prof. Bernhard Ø. Palsson of the Department of Bioengineering, UCSD, and members of his lab, Mr. Aarash Bordbar in particular, for helpful discussions and assistance in this study during his visit.
Authors’ Affiliations
References
- Ideker T, Galitski T, Hood L: A new approach to decoding life: systems biology. Annu Rev Genomics Hum Genet. 2001, 2: 343-372. 10.1146/annurev.genom.2.1.343.View ArticlePubMedGoogle Scholar
- Kitano H: Systems biology: a brief overview. Science. 2002, 295 (5560): 1662-1664. 10.1126/science.1069492.View ArticlePubMedGoogle Scholar
- Weinberger LS, Burnett JC, Toettcher JE, Arkin AP, Schaffer DV: Stochastic gene expression in a lentiviral positive-feedback loop: HIV-1 Tat fluctuations drive phenotypic diversity. Cell. 2005, 122 (2): 169-182. 10.1016/j.cell.2005.06.006.View ArticlePubMedGoogle Scholar
- Sotiriou C, Piccart MJ: Opinion - Taking gene-expression profiling to the clinic: when will molecular signatures become relevant to patient care?. Nat Rev Cancer. 2007, 7 (7): 545-553. 10.1038/nrc2173.View ArticlePubMedGoogle Scholar
- Acar M, Mettetal JT, van Oudenaarden A: Stochastic switching as a survival strategy in fluctuating environments. Nat Genet. 2008, 40 (4): 471-475. 10.1038/ng.110.View ArticlePubMedGoogle Scholar
- Luo J, Emanuele MJ, Li D, Creighton CJ, Schlabach MR, Westbrook TF, Wong KK, Elledge SJ: A genome-wide RNAi screen identifies multiple synthetic lethal interactions with the Ras oncogene. Cell. 2009, 137 (5): 835-848. 10.1016/j.cell.2009.05.006.PubMed CentralView ArticlePubMedGoogle Scholar
- Costanzo M, Baryshnikova A, Bellay J, Kim Y, Spear ED, Sevier CS, Ding H, Koh JL, Toufighi K, Mostafavi S: The genetic landscape of a cell. Science. 2010, 327 (5964): 425-431. 10.1126/science.1180823.View ArticlePubMedGoogle Scholar
- Spencer SL, Sorger PK: Measuring and modeling apoptosis in single cells. Cell. 2011, 144 (6): 926-939. 10.1016/j.cell.2011.03.002.PubMed CentralView ArticlePubMedGoogle Scholar
- Chen KC, Calzone L, Csikasz-Nagy A, Cross FR, Novak B, Tyson JJ: Integrative analysis of cell cycle control in budding yeast. Mol Biol Cell. 2004, 15 (8): 3841-3862. 10.1091/mbc.E03-11-0794.PubMed CentralView ArticlePubMedGoogle Scholar
- Palsson B: Systems biology: properties of reconstructed networks. 2006, Cambridge; New York: Cambridge University PressView ArticleGoogle Scholar
- Alon U: An introduction to systems biology: design principles of biological circuits. 2007, Chapman & Hall/CRC: Boca Raton, FLGoogle Scholar
- Chuang HY, Hofree M, Ideker T: A decade of systems biology. Annu Rev Cell Dev Biol. 2010, 26: 721-744. 10.1146/annurev-cellbio-100109-104122.PubMed CentralView ArticlePubMedGoogle Scholar
- Arkin AP, Schaffer DV: Network news: innovations in 21st century systems biology. Cell. 2011, 144 (6): 844-849. 10.1016/j.cell.2011.03.008.View ArticlePubMedGoogle Scholar
- Palsson B: Systems biology: simulation of dynamic network states. 2011, Cambridge, UK; New York: Cambridge University PressView ArticleGoogle Scholar
- Zamora-Sillero E, Hafner M, Ibig A, Stelling J, Wagner A: Efficient characterization of high-dimensional parameter spaces for systems biology. BMC Syst Biol. 2011, 5: 142-10.1186/1752-0509-5-142.PubMed CentralView ArticlePubMedGoogle Scholar
- Chen BS, Chen PW: GA-based Design Algorithms for the Robust Synthetic Genetic Oscillators with Prescribed Amplitude, Period and Phase. Gene Regul Syst Biol. 2010, 4: 35-52.View ArticleGoogle Scholar
- Batt G, Yordanov B, Weiss R, Belta C: Robustness analysis and tuning of synthetic gene networks. Bioinformatics. 2007, 23 (18): 2415-2422. 10.1093/bioinformatics/btm362.View ArticlePubMedGoogle Scholar
- Berg HC, Brown DA: Chemotaxis in Escherichia coli analysed by three-dimensional tracking. Nature. 1972, 239 (5374): 500-504. 10.1038/239500a0.View ArticlePubMedGoogle Scholar
- Macnab RM, Koshland DE: The gradient-sensing mechanism in bacterial chemotaxis. Proc Natl Acad Sci USA. 1972, 69 (9): 2509-2512. 10.1073/pnas.69.9.2509.PubMed CentralView ArticlePubMedGoogle Scholar
- Berg HC, Tedesco PM: Transient response to chemotactic stimuli in Escherichia coli. Proc Natl Acad Sci USA. 1975, 72 (8): 3235-3239. 10.1073/pnas.72.8.3235.PubMed CentralView ArticlePubMedGoogle Scholar
- Mello BA, Tu Y: Perfect and near-perfect adaptation in a model of bacterial chemotaxis. Biophys J. 2003, 84 (5): 2943-2956. 10.1016/S0006-3495(03)70021-6.PubMed CentralView ArticlePubMedGoogle Scholar
- Zigmond SH, Sullivan SJ: Sensory Adaptation of Leukocytes to Chemotactic Peptides. J Cell Biol. 1979, 82 (2): 517-527. 10.1083/jcb.82.2.517.View ArticlePubMedGoogle Scholar
- Shi W, Zusman DR: Sensory adaptation during negative chemotaxis in Myxococcus xanthus. J Bacteriol. 1994, 176 (5): 1517-1520.PubMed CentralPubMedGoogle Scholar
- Marwan W, Bibikov SI, Montrone M, Oesterhelt D: Mechanism of photosensory adaptation in Halobacterium salinarium. J Mol Biol. 1995, 246 (4): 493-499. 10.1006/jmbi.1994.0101.View ArticlePubMedGoogle Scholar
- Hilliard MA, Apicella AJ, Kerr R, Suzuki H, Bazzicalupo P, Schafer WR: In vivo imaging of C. elegans ASH neurons: cellular response and adaptation to chemical repellents. EMBO J. 2005, 24 (1): 63-72. 10.1038/sj.emboj.7600493.PubMed CentralView ArticlePubMedGoogle Scholar
- Jaasma MJ, Jackson WM, Tang RY, Keaveny TM: Adaptation of cellular mechanical behavior to mechanical loading for osteoblastic cells. J Biomech. 2007, 40 (9): 1938-1945. 10.1016/j.jbiomech.2006.09.010.View ArticlePubMedGoogle Scholar
- Muzzey D, Gomez-Uribe CA, Mettetal JT, van Oudenaarden A: A systems-level analysis of perfect adaptation in yeast osmoregulation. Cell. 2009, 138 (1): 160-171. 10.1016/j.cell.2009.04.047.PubMed CentralView ArticlePubMedGoogle Scholar
- Spehr J, Hagendorf S, Weiss J, Spehr M, Leinders-Zufall T, Zufall F: Ca2+ − calmodulin feedback mediates sensory adaptation and inhibits pheromone-sensitive ion channels in the vomeronasal organ. J Neurosci. 2009, 29 (7): 2125-2135. 10.1523/JNEUROSCI.5416-08.2009.View ArticlePubMedGoogle Scholar
- Barkai N, Leibler S: Robustness in simple biochemical networks. Nature. 1997, 387 (6636): 913-917. 10.1038/43199.View ArticlePubMedGoogle Scholar
- Yi TM, Huang Y, Simon MI, Doyle J: Robust perfect adaptation in bacterial chemotaxis through integral feedback control. Proc Natl Acad Sci USA. 2000, 97 (9): 4649-4653. 10.1073/pnas.97.9.4649.PubMed CentralView ArticlePubMedGoogle Scholar
- Wadhams GH, Armitage JP: Making sense of it all: Bacterial chemotaxis. Nat Rev Mol Cell Biol. 2004, 5 (12): 1024-1037. 10.1038/nrm1524.View ArticlePubMedGoogle Scholar
- Porter SL, Wadhams GH, Armitage JP: Signal processing in complex chemotaxis pathways. Nat Rev Microbiol. 2011, 9 (3): 153-165. 10.1038/nrmicro2505.View ArticlePubMedGoogle Scholar
- Ma WZ, Trusina A, El-Samad H, Lim WA, Tang C: Defining Network Topologies that Can Achieve Biochemical Adaptation. Cell. 2009, 138 (4): 760-773. 10.1016/j.cell.2009.06.013.PubMed CentralView ArticlePubMedGoogle Scholar
- Sourjik V, Berg HC: Binding of the Escherichia coli response regulator CheY to its target measured in vivo by fluorescence resonance energy transfer. Proc Natl Acad Sci USA. 2002, 99 (20): 12669-12674. 10.1073/pnas.192463199.PubMed CentralView ArticlePubMedGoogle Scholar
- Iman RL, Helton JC, Campbell JE: An Approach to Sensitivity Analysis of Computer-Models.1. Introduction, Input Variable Selection and Preliminary Variable Assessment. J Qual Technol. 1981, 13 (3): 174-183.Google Scholar
- Kollmann M, Lovdok L, Bartholome K, Timmer J, Sourjik V: Design principles of a bacterial signalling network. Nature. 2005, 438 (7067): 504-507. 10.1038/nature04228.View ArticlePubMedGoogle Scholar
- Morton-Firth CJ, Shimizu TS, Bray D: A free-energy-based stochastic simulation of the Tar receptor complex. J Mol Biol. 1999, 286 (4): 1059-1074. 10.1006/jmbi.1999.2535.View ArticlePubMedGoogle Scholar
- Rao CV, Kirby JR, Arkin AP: Design and diversity in bacterial chemotaxis: a comparative study in Escherichia coli and Bacillus subtilis. PLoS Biol. 2004, 2 (2): E49-10.1371/journal.pbio.0020049.PubMed CentralView ArticlePubMedGoogle Scholar
- Khan S, Castellano F, Spudich JL, McCray JA, Goody RS, Reid GP, Trentham DR: Excitatory signaling in bacterial probed by caged chemoeffectors. Biophys J. 1993, 65 (6): 2368-2382. 10.1016/S0006-3495(93)81317-1.PubMed CentralView ArticlePubMedGoogle Scholar
- Alon U, Surette MG, Barkai N, Leibler S: Robustness in bacterial chemotaxis. Nature. 1999, 397 (6715): 168-171. 10.1038/16483.View ArticlePubMedGoogle Scholar
- Henson R, Cetto L: The MATLAB bioinformatics toolbox. Encyclopedia of Genetics, Genomics, Proteomics and Bioinformatics. 2005, Natick, MA, USA: The MathWorks, IncGoogle Scholar
- Chiang AWT, Hwang MJ: A computational pipeline for identifying kinetic motifs to aid in the design and improvement of synthetic gene circuits. BMC Bioinformatics. 2013, 14 (Suppl 16): S5-10.1186/1471-2105-14-S16-S5.PubMed CentralView ArticlePubMedGoogle Scholar
- Mootha VK, Lindgren CM, Eriksson KF, Subramanian A, Sihag S, Lehar J, Puigserver P, Carlsson E, Ridderstrale M, Laurila E: PGC-1alpha-responsive genes involved in oxidative phosphorylation are coordinately downregulated in human diabetes. Nat Genet. 2003, 34 (3): 267-273. 10.1038/ng1180.View ArticlePubMedGoogle Scholar
- Subramanian A, Tamayo P, Mootha VK, Mukherjee S, Ebert BL, Gillette MA, Paulovich A, Pomeroy SL, Golub TR, Lander ES: Gene set enrichment analysis: a knowledge-based approach for interpreting genome-wide expression profiles. Proc Natl Acad Sci USA. 2005, 102 (43): 15545-15550. 10.1073/pnas.0506580102.PubMed CentralView ArticlePubMedGoogle Scholar
- Wilcoxon F: Individual Comparisons by Ranking Methods. Biometrics Bull. 1945, 1 (6): 80-83. 10.2307/3001968.View ArticleGoogle Scholar
- Vuppula RR, Tirumkudulu MS, Venkatesh KV: Mathematical modeling and experimental validation of chemotaxis under controlled gradients of methyl-aspartate in Escherichia coli. Mol Biosyst. 2010, 6 (6): 1082-1092. 10.1039/b924368b.View ArticlePubMedGoogle Scholar
- Stewart RC: Activating and inhibitory mutations in the regulatory domain of CheB, the methylesterase in bacterial chemotaxis. J Biol Chem. 1993, 268 (3): 1921-1930.PubMedGoogle Scholar
- Stewart RC, Jahreis K, Parkinson JS: Rapid phosphotransfer to CheY from a CheA protein lacking the CheY-binding domain. Biochemistry. 2000, 39 (43): 13157-13165. 10.1021/bi001100k.View ArticlePubMedGoogle Scholar
- Rogers A, Gibon Y: Chapter 4. Enzyme kinetics: theory and practice. Plant metabolic networks. Edited by: Schwender J. 2009, New York: Springer, 71-103.View ArticleGoogle Scholar
- Goldbeter A, Koshland DE: An amplified sensitivity arising from covalent modification in biological systems. Proc Natl Acad Sci USA. 1981, 78 (11): 6840-6844. 10.1073/pnas.78.11.6840.PubMed CentralView ArticlePubMedGoogle Scholar
- Lupas A, Stock J: Phosphorylation of an N-terminal regulatory domain activates the CheB methylesterase in bacterial chemotaxis. J Biol Chem. 1989, 264 (29): 17337-17342.PubMedGoogle Scholar
- Li ZF, Osborne MR, Prvan T: Parameter estimation of ordinary differential equations. Ima J Numer Anal. 2005, 25 (2): 264-285. 10.1093/imanum/drh016.View ArticleGoogle Scholar
- Chou IC, Voit EO: Recent developments in parameter estimation and structure identification of biochemical and genomic systems. Math Biosci. 2009, 219 (2): 57-83. 10.1016/j.mbs.2009.03.002.PubMed CentralView ArticlePubMedGoogle Scholar
- Schellenberger J, Zielinski DC, Choi W, Madireddi S, Portnoy V, Scott DA, Reed JL, Osterman AL, Palsson B: Predicting outcomes of steady-state 13C isotope tracing experiments using Monte Carlo sampling. BMC Syst Biol. 2012, 6: 9-10.1186/1752-0509-6-9.PubMed CentralView ArticlePubMedGoogle Scholar
- Gonze D, Halloy J, Goldbeter A: Robustness of circadian rhythms with respect to molecular noise. Proc Natl Acad Sci USA. 2002, 99 (2): 673-678. 10.1073/pnas.022628299.PubMed CentralView ArticlePubMedGoogle Scholar
- Zak DE, Stelling J, Doyle FJ: Sensitivity analysis of oscillatory (bio)chemical systems. Comput Chem Eng. 2005, 29 (3): 663-673. 10.1016/j.compchemeng.2004.08.021.View ArticleGoogle Scholar
- Szallasi Z, Stelling J, Periwal V: System modeling in cell biology: from concepts to nuts and bolts. 2006, Cambridge, Mass: MIT PressView ArticleGoogle Scholar
- Balsa-Canto E, Banga JR: AMIGO, a toolbox for advanced model identification in systems biology using global optimization. Bioinformatics. 2011, 27 (16): 2311-2313. 10.1093/bioinformatics/btr370.PubMed CentralView ArticlePubMedGoogle Scholar
- Wolf DM, Arkin AP: Motifs, modules and games in bacteria. Curr Opin Microbiol. 2003, 6 (2): 125-134. 10.1016/S1369-5274(03)00033-X.View ArticlePubMedGoogle Scholar
- Alon U: Network motifs: theory and experimental approaches. Nat Rev Genet. 2007, 8 (6): 450-461. 10.1038/nrg2102.View ArticlePubMedGoogle Scholar
- Tyson JJ, Novak B: Functional motifs in biochemical reaction networks. Annu Rev Phys Chem. 2010, 61: 219-240. 10.1146/annurev.physchem.012809.103457.PubMed CentralView ArticlePubMedGoogle Scholar
- Kitano H: Biological robustness. Nat Rev Genet. 2004, 5 (11): 826-837. 10.1038/nrg1471.View ArticlePubMedGoogle Scholar
- Hartwell LH, Hopfield JJ, Leibler S, Murray AW: From molecular to modular cell biology. Nature. 1999, 402 (6761 Suppl): C47-C52.View ArticlePubMedGoogle Scholar
- Ihmels J, Friedlander G, Bergmann S, Sarig O, Ziv Y, Barkai N: Revealing modular organization in the yeast transcriptional network. Nat Genet. 2002, 31 (4): 370-377.PubMedGoogle Scholar
- Schlosser G, Wagner GP: Modularity in development and evolution. 2004, Chicago: University of Chicago PressGoogle Scholar
- Shen-Orr SS, Milo R, Mangan S, Alon U: Network motifs in the transcriptional regulation network of Escherichia coli. Nat Genet. 2002, 31 (1): 64-68. 10.1038/ng881.View ArticlePubMedGoogle Scholar
- Brandman O, Meyer T: Feedback loops shape cellular signals in space and time. Science. 2008, 322 (5900): 390-395. 10.1126/science.1160617.PubMed CentralView ArticlePubMedGoogle Scholar
- Tsai TY, Choi YS, Ma W, Pomerening JR, Tang C, Ferrell JE: Robust, tunable biological oscillations from interlinked positive and negative feedback loops. Science. 2008, 321 (5885): 126-129. 10.1126/science.1156951.PubMed CentralView ArticlePubMedGoogle Scholar
- Famili I, Forster J, Nielsen J, Palsson BO: Saccharomyces cerevisiae phenotypes can be predicted by using constraint-based analysis of a genome-scale reconstructed metabolic network. Proc Natl Acad Sci USA. 2003, 100 (23): 13134-13139. 10.1073/pnas.2235812100.PubMed CentralView ArticlePubMedGoogle Scholar
- Liu J: Kinetic constraints for formation of steady states in biochemical networks. Biophys J. 2005, 88 (5): 3212-3223. 10.1529/biophysj.104.056085.PubMed CentralView ArticlePubMedGoogle Scholar
Copyright
This article is published under license to BioMed Central Ltd. This is an open access article distributed under the terms of the Creative Commons Attribution License (http://creativecommons.org/licenses/by/2.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.