Studying system-level properties of bio-molecules is essential to systems biology [1
]. But most studies are based on either network topology that is not working very well at the protein level, or flux analysis that lacks in system level perspective [13
]. To overcome such drawbacks, we propose a method of criticality characterization on the basis of kinetic modeling. In a kinetic system, every interaction is expressed by a kinetic rate equation. How a component influences the system is determined by both its position and the kinetic parameters. Position is equivalent to topological property, while kinetic parameters encode specific biochemical/biological functions. Both kinds of information are integrated in modeling and revealed by dynamical simulation [15
]. According to the typical formulism of biochemical systems, the kinetic rate equations constitute the deterministic part of the complex system dynamics and they can be viewed as the "driving force" of the system [23
]. Thus theoretically, the criticality characterization proposed in our method is the study of structural factors built into the "driving force" of a system.
Differing from topology-based methods, our method characterizes system-level properties on the quantitative basis. But unlike the conventional sensitivity analysis, we employ the network structure information by calculating the distances from the deleted spot to the affected entities besides computing the fluctuations. Moreover, unlike conventional flux-analysis approaches, we explore the system stability and retrieve system dynamics structure. Incorporating the network/dynamics structure information allows us to reveal the simultaneous/collateral influences and the overall impact on the system. Another major difference from the sensitivity analysis is that we use in silico
deletions instead of mild perturbations (e.g. 5% or 10%, as most flux-analysis approaches do). Because a well-casted biological network usually has parametric properties favouring the robustness in dynamics, critical components may well tolerate mild perturbations (i.e. parameters exhibiting the Lyapunov stability). Therefore, individual sensitivity analysis often fails to identify such critical spots, and its inability to reveal simultaneous influences worsens the situation. That is why we develop the "criticality characterization". In silico
deletion is equivalent to investigating how the system would be if the component is forcefully assumed to be absent, eliminating the parametric properties stated earlier. Furthermore, our method's capacity of revealing simultaneous/overall impacts at the system level enables it to distinguish real critical spots from uncritical ones more effectively. In addition, utilizing kinetic model as the analytical basis is a superiority over the stoichiometric flux-balance modeling in traditional flux-analysis methods, enabling us to appropriately explore system behaviours in the real-time scale [15
]. For example, both traditional topology-based and flux-analysis approaches regard TIS as peripheral as it is not highly connected and it is not on any uni-directional or rate limiting reactions. However, there were experimental studies showing that knocking out tpiA
(i.e. the gene encoding TIS) attenuated the cell growth by about 70%. And our method appropriately revealed that TIS could exert large impacts on the system if deleted, because of the designs we made (mentioned above). Hence methodologically, our method creates a different angle from topology-based methods and can be viewed as an improvement of conventional flux-analysis approaches.
After in silico
deleting a protein, the residual system is actually a virtual structure. We assume that this structure encodes important information about whether the mutant can maintain its functionality and how it would dynamically behave/evolve provided that it stills operates on its own. The residual system fails to maintain functionality when its kinetics goes beyond regular ranges (e.g. occurring negative values, or soaring to extreme values exceeding regular intracellular molecular concentrations), or its dynamics is trapped in a mode where the stable equilibrium is sabotaged, as stable equilibrium gives rise to robustness and is an essential prerequisite for valid mathematical formulations of living cell dynamics [17
]. Either case indicates that deleting the protein makes the system so ill-suited that it cannot run on its own.
By applying our method to E. coli
central carbon metabolism, we find that deleting enzymes such as PGK, GAPDH, etc. causes the system to become a very ill-suited structure as some state values soaring to levels beyond the normal range and the trajectories are highly divergent throughout the state space (Figure and ). Likewise, deleting enzymes such as TKb, ALDO, etc. also causes relatively large impacts on both system kinetics and qualitative dynamics (Figure and ). On the contrary, knocking out enzymes such as PGI, PK, etc. exerts very small influences (Figure and ). We also find enzymes can mediate large influences on distant metabolites or enzymes. For instance, TKb, PGK, PGI, etc. can all exert the largest impacts on entities of distances other than 1 (Figure , and Additional file 2
). This is because bio-systems have complex structures consisting of branches, alternative pathways and loops, as well as various kinetic parameters differing in orders of magnitudes [6
]. Such structure acts as a special leverage, determining special ways of interactions and influence propagations. Only kinetic modeling can reveal such knowledge, and such analyses can give us more clues on selecting potential regulatory targets for use in drug development, metabolic engineering, etc.
Figure 6 Remote impact and local flux compensation. A system component can affect distant entities even greater than its closest neighbours, which is illustrated by TKb. Moreover, asymmetry in the criticalities of components can result in local flux compensation, (more ...)
By utilizing the power of kinetic model for approaching real-time events, we simulated fluxes after enzyme deletions and compared the results with a previously validated study of metabolic essentiality [6
]. The comparison shows that our characterization of criticality is well supported by functional essentiality. Interestingly, we discovered that the asymmetry in criticalities of building blocks might give rise to local flux compensation. For instance, multiple metabolites (e.g. ribulose-5-phosphate, sedoheptulose-7- phosphate, etc.) in the pentose phosphate pathway have increased flux-sums after PGI knockout (Figure and ). The cutoff of PGI induces the two alternative pathways for generating the essential metabolite fructose-6-phosphate, TKb and TA, to operate at a greater volume. Thus fluxes through relative reactions are compensated, resulting in local amplified fluxes. This is a likely result in accordance to the MOMA mechanism [17
]. Although MOMA can compensate system fluxes/states to some degrees, our results show that the effects caused by deletions of critical components such as TKb, TA, PGK, etc. cannot be smoothed by such compensations (Figure , Additional file 5
). This is because such compensations are mainly mediated by alternative pathways [6
]. When a critical component is deleted, leaving inferior components as backup to rely on, the system cannot work efficiently. On the contrary, deleting PGI leaves its two alternative pathways that are of superior properties at the "ON" state and the system still works, thus fluxes/states can be efficiently compensated. This gives a hint on how criticality characterization can help in bio-system modifications such as in metabolic engineering. We can delete some system components with inferior properties, leaving alternative pathways with superior properties to work. And phenotypes in the local areas relating to such alternative pathways might be compensated due to the leverage of system structure and the MOMA mechanism. Therefore, comprehensive methods of exploring system-level properties can help us make use of bio-complexity in engineering, as well as in knowledge discovery.
It is noteworthy that functionally important components are not necessarily critical, as studies suggest that the more important a reaction is in function, the more likely that it has a backup pathway [6
]. For example, PK connects very fundamental chemical compounds but it is regarded as uncritical at the system level, because there are alternative paths (e.g. the phosphotransferase system - PTS, in bacteria glycolysis; Additional file 1
) that can prevent large impacts on system kinetics/dynamics. This exemplifies that bio-system components have dichotomy. They have "importance" as biochemical molecules, and they also have "criticality" to the system as constitutive building blocks. Actually, our method does not aim to find the "functionally important" molecules, but those "critical" to the system, i.e. components that cannot be absent, or the system will be severely aberrant. Since the criticality of an enzyme depends on many factors (e.g. kinetic parameters, substrates inhibiting/activating other reactions, degree of the effects, etc.), the assignments of system boundaries in modeling might affect prediction results. As the enzymes located on the boundary might have incomplete interplay structure, the above factors may not occur properly in the kinetic equations. Therefore, accurate criticality characterization is facilitated by appropriate system inclusiveness in modeling. For example, glucose-1-phosphate adenyltransferase (G1PAT) only connects the external polysaccharide synthesis pathway, with few interactions with large-capacity reactions both in the system and outside pathways. Thus as the boundary is assigned up to it, the validity of the results are enhanced (Table ). Furthermore, fundamental, common and conserved pathways must be chosen for comparison with genome-scale gene essentiality studies that regard to global cellular functionality. For instance, the bacteria central carbon metabolism here is an appropriate example [5
], thus various predictions of protein criticality are well consistent with global gene essentiality characterizations [13
Although we used a metabolic system as the working model, the application of our method is not confined to metabolic systems. For instance, we can model gene transcription dynamics by deriving gene transcription rate with the power-law formulism, the Hill equation, or equations of chemical kinetic actions [27
]. Or we can describe ligand-receptor and protein-protein binding actions with the mass action law and build models for signaling networks [22
]. We even do not have to obtain exact parameters fitting the modeled solutions to assay measurements when analyzing the generic behavioral potential of the system, e.g. in what parameter ranges the system exhibits certain dynamics and how they change with parameters. Such qualitative predictions are also useful in revealing general principles governing complex bio-systems. Naturally, complicated bifurcation dynamics will be harder to analyze; but the idea of our method can be well applied once the coexisting dynamical characteristics in bifurcation are associated with biological implications [28
]. By integrating knowledge and using theoretical generic forms of models [15
], kinetic modeling will be eventually feasible for more organisms. Hence, instead of the traditional approaches, we propose that complex systems be studied by casting the network into kinetic equations and computing the system-level properties with respect to system kinetics/dynamics (criticality). Overall, our method may provide a new viewpoint in revealing constitutive/functional properties of building blocks in a biological system.