Mathematical description of biochemical reactions is deeply rooted in chemical kinetics. At its very core, a system of chemical reactions may be written in the form of the *Law of Mass Action*

where

*α *_{i} β_{i}, are the kinetic rates, and

*P*_{i m},

*Q*_{i m} are the

*matrices of stoichiometric coefficients* in the direct and inverse reactions, respectively [

57,

58]. Depending on the nature and complexity of the system under investigation, the quantities{

*x*_{i}} may represent concentration of various biochemical constituents participating in the process, including individual molecules or their aggregates. There is no unique way of representing the biochemical machinery in mathematical form: depending on the level of

*structural granularity* and temporal resolution, the same process may be seen either as an individual chemical reaction or as a complicated system of reactions. Formally,

equations (1.1) are appropriate for a system in which each constituent is generated by one direct and one reverse reaction. The reality of large biochemical systems is, of course, far more complex. In particular, there may be several competitive reactions producing and degrading the same constituents but following different intermediate pathways. For these cases, a more appropriate form of the equations would be

known as the law of

*Generalized Mass Action.* Here

*L*_{i},

*M*_{i} are the numbers of concurrent reactions of production and degradation,

*α*_{ni} β_{ni}, are the

*matrices of rates*, and

*P*_{n i m},

*Q*_{n i} m are the

*tensors of stoichiometric coefficients*. In principle, however, this more complex system is reducible to the form (

1.1) by appropriate redefinition of chemical constituents [

59].

The

*Power Law Formalism* expresesed by the

equation (1.1) has been introduced by Savageau in 1969 as a method for studying complex biochemical phenomena [

57,

60,

61]. In ensuing years, this method has been applied to numerous problems in biology [

62–

69]. Later on, the Power Law Formalism has been applied in a more general context of

*organizationally complex systems*, beyond purely biological meaning of this word [

70,

71]. The term

*S-Systems* appeared later for denoting the

*synergy-saturation systems*. In addition to its subject matter meaning, this term was useful for distinguishing the S-Systems methodology from other methodologies containing the key words

*power law* in their names. Extensive practical work with S-Systems revealed many useful analytical properties which put them into a position of a much more potent instrument than it was originally anticipated. In particular, it turned out that similarly to

*neural networks* the S-Systems may serve as

*universal approximators*, which means that

*in principle* any system of multidimensional nonlinear differential equations may be recast into the form of S-System [

72]. Furthermore, it was discovered that S-Systems can be transformed into the

*Lotka-Volterra system* known in

*population dynamics*; these two nonlinear dynamical models have been developed within two completely independent contexts and were originally thought to be unrelated to each other [

73]. In addition, the S-Systems has shown their usefulness in dealing with purely computational problems, such as those in probability theory [

74], in numerical integration of differential equations[

75], in finding the roots of algebraic equations [

76]. A comprehensive summary of these developments, as well as some newer results, are given in [

58].

In the S-System approach, the state of a system is characterized in terms of time-dependent *constituents*, {*x*_{i} (*t*)}. These constituents, as well as kinetic rates and stoichiometric matrices, should be appropriately defined at each level of biological organization. To make this premise clearer, let us imagine that on a certain level of abstraction it could be quite sufficient to describe gene expression by the following pair of biochemical reactions

where TF and RNAP stand for transcription factor and RNA Polymerase, respectively. The structure of stoichiometric matrices is quite obvious from this notation. As to the kinetic rates, they should be defined and measured on this particular level of abstraction. A more close look at the processes symbolically depicted as a chemical reaction (

1.3) would reveal a much more complicated picture involving hundreds of biomolecules and thousands of elemental steps, each representing a separate chemical reaction [

31,

77]. If one would like to create a model accounting for this level of details then the definitions of

*constituents* may drastically change, and the model would become much more complicated and rich in properties. However, in the S-Systems methodology, the

*functional form of the model* (

1.1) will remain the same. That’s why the S-System is not just another model, one of many possible. It is a

*unified framework of thinking* about big and complex hierarchical systems. The originators of the S-System touched upon this fundamental property in many of their works. They refer to the process of convolution of a detailed model of lower level of abstraction into a more concise but more coarse-grained model at a higher level of abstraction as

*aggregation*. They write in [

78]: “There are alternative hierarchical levels at which biochemical systems can be described by this formalism. The lowest level might correspond to the elemental chemical kinetic description for each of the steps in each of the enzyme-catalyzed reactions. At this level, the power-law description is exactly the same as that of conventional chemical kinetics... At a higher level, one can aggregate system components, and the net rates of increase and decrease of these aggregate variables are again mathematical expressions that can be represented by power-law functions. At a still higher level, one can describe the growth of organisms and their interactions at a biochemical level by means of this same power-law formalism.”

Obviously, any hierarchical model of reality, however complex it is, is always open at the top towards further generalization, and at the bottom towards further elaboration and enrichment. In the S-System, however, the

*functional form of dynamics* remains the same while the focus of attention moves from the bottom to the top in this hierarchy. Voit and Savageau called this property of S-Systems a

*telescopic property* [

70,

78].

A major limitation of ordinary differential equations (ODE), such as (

1.1) and (

1.2), in application to intracellular biochemistry is that they are valid only for the

*well stirred* systems. In such systems, any chemical constituent produced anywhere in the system becomes immediately available for all other chemical transformations. As discussed in Section 2, in the tight nuclear environment this assumption may be questionable because of slow delivery of transcription factors to regulatory sites. A way of overcoming this drawback is in accounting for spatial distributions of biochemical constituents and compartmentalization of the model. These measures automatically lead to replacement of the system of ODEs by a system of partial differential equations (PDE) with diffusion terms. A detailed review of these approaches may be found in [

79].

As with any discipline seeking support in mathematics, a fundamental limitation comes from the fact that any mathematical model leaves out of its scope a vast universe of

*unmodeled realities* thus introducing uncontrollable distortions into any theoretical description. During centuries-long history of this question, a number of powerful techniques have been developed to address the problem. In mathematical biology, existence and influence of the unmodeled realities is usually taken into account through introducing the concept of

*intrinsic noise*. Formally, this measure leads to the replacement of the ODE of chemical kinetics by the system of stochastic differential equations (SDE) thus enabling reformulation of dynamics in probabilistic terms [

13,

20,

80]. It should be noted, however, that the very usage of the word

*noise* has a connotation that

*in principle* the system allows for deterministic description and it is only the nuisance factors beyond our control that make such a description impossible. Such a premise is not always justifiable even in low-dimensional nonlinear systems (with Lorenz attractor being a celebrated example [

81]) in which the dynamics may assume the form of

*deterministic chaos*. This means that the motion in the system, being in principle deterministic and fully predictable, is nevertheless so complex and tangled that may formally satisfy the criteria of stochasticity; hence, for all practical purposes they may be regarded as random [

82]. This situation is especially likely in nonlinear systems of very high dimension with strong interactions between components. Equations of chemical kinetics represent a perfect example of such systems. Many aspects of stochastic behavior in strongly nonlinear high dimensional systems has been discussed earlier in the works [

23,

24,

83,

84] by SR.

Serious limitations of the basic model of chemical kinetics in applications to intracellular biomolecular dynamics forced researchers to look for alternative approaches. One such approach has emerged comparatively recently and represents a powerful combination of two, largely independent, mathematical disciplines: graph theory and nonlinear dynamics. It has been termed

*Chemical Reaction Networks* (CRN) theory, and in recent decades has been slowly but steadily percolating into intra-cellular biochemistry [

85,

86]. A major motivation for introducing the CRN and other graph-theoretical techniques is a regrettable, yet inescapable, fact that

*quantitative parameterization* of any large-scale biochemical model may be prohibitively time-, labor- and money-consuming, if feasible at all. Hence, in the vast majority of models such a parameterization may be simply unavailable. In addition, there is no guarantee that

*in vitro* measurements of kinetic rates can be truly representative of

*in vivo* dynamics. At the same time,

*qualitative* relations between the major biochemical players expressed in the form of graphical networks, such as the charts of metabolic pathways [

87], are often known much better. Therefore, prior to labor-consuming effort of dynamical modeling, it seems natural to derive as much as possible from the careful examination of the network’s topological structure. CRN theory provides a theoretical framework for such an approach and often allows for acquiring surprisingly rich information. In particular, topological analysis allows for making judgments regarding stability, vulnerability to damage, mechanisms of self-repair, scenarios of self-organization and other properties. Moreover, some versions of the network dynamics are capable of elucidating the mechanisms of prebiotic evolution through spontaneous birth of small

*autocatalytic molecular sets* followed by rapid self-organization into complex biomolecular structures [

88]. In general, topological analysis is capable of elucidating

*systemic properties*, that is, the ones that characterize the network’s

*global behavior*. These global structural properties are hardly possible to derive directly from the equations of chemical kinetics (see {Rosenfeld, 2008 721/id} for more detail.) The power of graph-theoretical methods is seen from the simple fact that neither kinetic rates nor actual concentrations of biochemical constituents are used in topological analyses. This means that the conclusions derived from the topological analysis are applicable to a wide class of systems possessing similar topologies, but otherwise entirely different [

87].