Benito Hernández-Bermejo

Renormalization group approach to power-law modeling of complex metabolic networks

Keywords: GMA models, S-systems, self-similarity, scaling, reaction rates

In the modeling of complex biological systems, and specially in the framework of the description of metabolic pathways, the use of power-law models (such as S-systems and GMA systems) often provides a remarkable accuracy over several orders of magnitude in concentrations, an unusually broad range not fully understood at present. In order to provide additional insight in this sense, this article is devoted to the renormalization group analysis of reactions in fractal or self-similar media. In particular, the renormalization group methodology is applied to the investigation of how rate-laws describing such reactions are transformed when the geometric scale is changed. The precise purpose of such analysis is to investigate whether or not power-law rate-laws present some remarkable features accounting for the successes of power-law modeling. As we shall see, according to the renormalization group point of view the answer is positive, as far as power-laws are the critical solutions of the renormalization group transformation, namely power-law rate-laws are the renormalization group invariant solutions. Moreover, it is shown that these results also imply invariance under the group of concentration scalings, thus accounting for the reported power-law model accuracy over several orders of magnitude in metabolite concentrations.

Power-law modeling

Power-law models constitute a highly structured nonlinear representation of many complex biological systems. Such models, originally developed on the basis of a first order

Taylor series in log-log space, have provided a very powerful modeling framework, specially (but not only) in the domain of biochemical pathways. For instance, see [START_REF] Savageau | Biochemical systems analysis: a study of function and design in molecular biology[END_REF]; [START_REF] Voit | Computational analysis of biochemical systems: a practical guide for biochemists and molecular biologists[END_REF]; [START_REF] Voit | Introduction to S-systems[END_REF] for a general and detailed description of the formalism and its applications. As indicated, one important advantage of such formalism is that the mathematical description of complex systems is highly structured, thus allowing a remarkable mathematical tractability, let it be in the form of S-systems

n i x b x a dt dx n j h j i n j g j i i ij ij , , 1 , 1 1 … = - = ∏ ∏ = =
or of generalized mass-action (or GMA, in what follows) systems:

n i x b x a dt dx q j n k h k ij p j n k g k ij i ijk ijk , , 1 , 1 1 1 1 … = - = ∑ ∏ ∑ ∏ = = = =
In both cases, and in the biochemical context, the system variables x i describe metabolite concentrations. Therefore, both S-systems and GMA systems are purely reactional models, the power-law functions describing the kinetic rate-laws (also termed velocity functions) of the reactions being modeled. An additional and very significant feature of modeling based on either S-systems or GMA systems is that, typically, such models are valid over several orders of magnitude in concentrations. This behavior has been widely reported in the literature [START_REF] Savageau | Biochemical systems analysis: a study of function and design in molecular biology[END_REF]Sorribas and Savageau 1989a, 1989b, 1989c;[START_REF] Voit | Accuracy of alternative representations for integrated biochemical systems[END_REF][START_REF] Voit | Introduction to S-systems[END_REF]. In fact, it is not uncommon that such models provide accurate predictions over variations of two or three orders of magnitude in concentrations, namely it is possible to find sizes as large as 100-or 1000-fold for the range of validity of A c c e p t e d m a n u s c r i p t 3 the models. In spite that power-law models were in principle conceived as precise nonlinear approximations, such success seems to suggest that more fundamental reasons might exist in such a way that power-law models could, in fact, be the natural description for many processes, specially in a biochemical framework. In the last context, an equivalent question is: to what extent power-law rate-laws do implement some fundamental properties of reactions, that could explain the accuracy of power-law models? In the literature, different approaches and explanations have been proposed in order to account for such features. In this sense, we can mainly mention the mathematical nature of the formalism, the existence of systemic regulations that maintain the concentrations within narrow limits, or the advantages of the agregation of interactions in the case of S-systems [START_REF] Voit | Introduction to S-systems[END_REF].

Other arguments are derived from considerations based on approximation theory (Savageau, 1979a(Savageau, , 1979b) ) as well as on the role of fractal kinetics [START_REF] Aon | The fractal architecture of cytoplasmatic organization: Scaling, kinetics and emergence in metabolic networks[END_REF][START_REF] Kopelman | Rate processes on fractals: theory, simulations and experiments[END_REF][START_REF] Kopelman | Reaction kinetics in restricted spaces[END_REF][START_REF] Savageau | Influence of fractal kinetics on molecular recognition[END_REF][START_REF] Savageau | Michaelis-Menten mechanism reconsidered: implications of fractal kynetics[END_REF][START_REF] Savageau | Development of fractal kinetic theory for enzyme-catalysed reactions and implications for the design of biochemical pathways[END_REF]. In particular, it is interesting to recall the important role of fractals in the context of in vivo chemistry, as far as many reactions are confined to two-dimensional membranes, one-dimensional channels, or fractal domains of non-integer dimension, and an alternative to the classical framework is necessary for the analysis of such phenomena. All the contributions just mentioned provide very valuable insights that help understanding the naturalness of the power-law function as a tool for modeling, in spite that the ultimate reasons accounting for this success seem to be not yet understood. In this context, the methods of group theory in general (and of the powerful tools of renormalization group in particular) have not been exploited in order to analyze the fundamental reasons that confer on power-laws such a special role. Actually, the different studies about power-law systems based on group-theoretic properties [START_REF] Díaz-Sierra | Graph-theoretic description of the interplay between non-linearity and connectivity in biological systems[END_REF][START_REF] Hernández-Bermejo | Lotka-Volterra representation of general nonlinear systems[END_REF][START_REF] Voit | Symmetries of S-systems[END_REF] are not devoted to the reported accuracy of power-law models. In fact, such papers deal mainly with mathematical properties of power-law systems, such as the existence of first integrals or reduction procedures. Precisely, the purpose of the present work is to make use of renormalization group methods in order to shed some light on the fundamental properties accounting for the success of power-law functions in the modeling through several orders of magnitude in concentrations. In first place, let us recall some general features regarding the renormalization group in order to provide a self-contained presentation. (1992); [START_REF] Fisher | Renormalization group theory: its basis and formulation in statistical physics[END_REF]; [START_REF] Schroeder | Fractals, chaos, power-laws[END_REF]; [START_REF] Takayasu | Fractals in the physical sciences[END_REF] for some clear introductions to renormalization group methods in statistical mechanics, including applications to fractals and self-similarity.

The purpose of the renormalization group is to treat quantitatively the change of a physical magnitude when the geometric scale is changed. Here the term "physical" is to be understood in a broad sense. For instance, let q be a certain physical quantity measured at a certain scale of coarse-graining, and let q * be the same quantity measured with the scale of coarse-graining changed by a factor λ. The rescaled value q * can be related to the original value of q by an appropriate transformation q * =f λ (q). This transformation f λ is the renormalization group transformation for magnitude q. Now the basis of the renormalization group analysis proceeds to the search of a fixed point of the renormalization group transformation, namely the search of a critical value q c such that q c =f λ (q c ). The value q c thus corresponds to a renormalization group invariant solution, as far as it remains unaltered after the application of the renormalization group transformation. In statistical mechanics, such solutions usually correspond to critical values of different processes such as phase transitions, percolation, etc. The reason is clear: after successive application of the renormalization group transformation, we find q c =f λ (f λ (q c )). This explains also the close relationship of the renormalization group invariant solutions and fractality. In fact, historically the notions of fractal and renormalization group appeared usually f λ is not the exact description of how q changes after the scaling. Instead, it is in general a simplified approximation in order to (i) have the possibility of explicitly writing a renormalization group transformation f λ , and (ii) find a transformation f λ that is amenable to analysis. This implies also that the solutions found by means of renormalization group techniques often provide a hint about the system properties, usually qualitative (and only approximate from a quantitative point of view) rather that leading to a precise analysis. Of course, the previous considerations imply also that a given problem can be analyzed by means of the renormalization group methods with progressively refined degrees of approximation.

In a biological context, the origin of power-law or allometric relationships has deserved a significant interest from many decades (e.g. see [START_REF] Huxley | Problems of relative growth[END_REF]. In this sense, very diverse contributions can be found in the literature, for instance those based on statistical approaches [START_REF] Kaitaniemi | Testing the allometric scaling laws[END_REF][START_REF] Packard | On the use of logarithmic transformations in allometric analysis[END_REF][START_REF] Wu | A statistical model for the genetic origin of allometric scaling laws in biology[END_REF], network theory [START_REF] Furusawa | Evolutionary origin of power-laws in a biochemical reaction network: embedding the distribution of abundance into topology[END_REF], and on physical arguments, mainly based on fractality and scaling principles (see [START_REF] Aon | The fractal architecture of cytoplasmatic organization: Scaling, kinetics and emergence in metabolic networks[END_REF][START_REF] Auffray | Scale relativity theory and integrative systems biology: 1 Founding principles and scale laws[END_REF][START_REF] Demetrius | The origin of allometric scaling laws in biology[END_REF][START_REF] West | The origin of universal scaling laws in biology[END_REF][START_REF] West | Allometric scaling of metabolic rate from molecules and mitochondria to cells and mammals[END_REF][START_REF] West | The origin of allometric scaling laws in biology from genomes to ecosystems: towards a quantitative unifying theory of biological structure and organization[END_REF]. In addition, some authors have specifically addressed the relevance of group theory and the renormalization group approach [START_REF] Derome | Biological similarity and group theory[END_REF][START_REF] West | Comments on the renormalization group, scaling and measures of complexity[END_REF] in such framework. However, as indicated in the previous subsection, the renormalization group methodology has not yet been applied in the context of powerlaw models (such as S-systems and GMA systems) and power-law rate-laws in order to account for the accuracy of the resulting models. As indicated, this is just the purpose of the present work. More precisely, in this article the renormalization group analysis of reactions in fractal or self-similar media is considered. In particular, the renormalization group methodology is applied to the investigation of how rate-laws describing such reactions are transformed when the geometric scale is changed. This will be done in two successive steps corresponding to different degrees of refinement in the renormalization group approximation. Such steps shall be termed Zero-Order renormalization group approximation (Section 2) and First-Order renormalization group approximation (Section 3). The precise purpose of such analysis is to investigate whether or not power-law rate-

One-variable case

Let us first consider the simplest situation in which a single metabolite A reacts in a given region of volume V in which the reacting species is restricted to a region of fractal dimension d. In the following, we shall denote by M the initial number of molecules in that volume, and therefore the concentration of the metabolite is x=M/V. In addition, let N be the increment in the number of molecules of A in that region per time unit due to the reaction. Thus, N>0 implies a net production of A, while N<0 means that A is depleted.

Accordingly, the rate of the reaction is given by v=N/V. Assume now that due to the underlying fractality, the system is self-similar after a scaling, and let λ be the scaling factor, namely we postulate the invariance of the system in a scale λ times larger, and let us denote with a star superscript the system parameters in the new scale. Obviously, now we have

V * = λ 3 V but notice that it is M * = λ d M.
Consequently, for the concentration we find

x * =M * /V * =λ d-3 x. Now, it is postulated that the invariance of the system in the new scale can be quantitatively expressed as:

* * * V M N MV N = (1)
In other words, the system invariance is assumed to mean that the fraction of molecules reacting per unit of time and per unit of volume remains invariant after the scaling. Thus, after (1) we readily find that

N * =NM * V * M -1 V -1 = λ d+3 N, as well as v * =N * /V * = λ d v.
Since the reactional rate-law v(x) is actually a function of the concentration, the scaling condition v * = λ d v must be explicitly written as an identity between functions, namely v * (ξ)= λ d v(ξ)

for any argument ξ. This identity for v is specially important for what is to follow, as far as it will be necessary for the application of the renormalization group condition. In fact, our main goal is to characterize the renormalization group invariant rate-laws. As indicated in the Introduction, the renormalization group condition assumes the scale invariance of the system, in our case in the context of the invariance of the rate-law. Consequently, in our To begin with, let us consider the case of a power-law functional form: v(x)=v 0 x g .

We thus have for the renormalization group condition:

v * (x * ) = λ d v(λ d-3 x) = λ d+g(d-3) v 0 x g = v 0 x g = v(x).
Then, after equating both functions, it is found that the renormalization group condition amounts to: λ d+g (d-3) =1. This equation admits a solution that actually provides invariance for all scales (i.e. for all values of λ). Such solution relates the renormalization group value of the kinetic order g with the fractal dimension d of the region. Thus the renormalization group value is:

3 0 , 3 < < - = d d d g (2)
This result presents several interesting features. First of all, notice that g is positive, but can take arbitrarily large values. It is worth emphasizing that the renormalization group method is an approximation. Let us ignore for the moment the formal divergence at d=3, which is a consequence of the approximate character of the renormalization group approach in the Zero-Order case (note that for d=3 the only solution corresponds to the transformation for which λ=1). Such divergence will dissapear in the First-Order renormalization group approach to be presented later. The important feature at this stage is that the kinetic order increases when the fractal dimension grows, and such kinetic order can take a wide range of possible values. This is a reasonable property since the dimensional restriction of the medium should limit the ability of the metabolite to react. In fact, we see that g tends to zero when d also tends to zero, which is to be expected. The most important conclusion here is that, in this approximation, the power-law rate-law is an exact solution of the renormalization group condition for all values of λ.

Let us now compare this result with other possible functional forms. For instance, let us first turn to kinetic functions of the kind v(x)=v 0 x g /(K+ x g ), with g>0. This form Once the renormalization group condition v * (x * )= v(x) is imposed, we readily find:

g d g d g d g x K x K ) 3 ( ) 3 ( 1 - - + + = + λ λ
Clearly, this outcome does not allow the renormalization group invariance in general.

Anyway, such invariance can be approximately found in the limit case K>> x g . In this situation, i.e. if x g can be considered negligible when compared to K (for instance, in the case of sufficiently low concentration) then the kinetic function can be approximated by a power-law: v(x)≈(v 0 /K)x g . Of course, this functional dependence was already considered before. Accordingly, kinetic functions of the form v(x)=v 0 x g /(K+ x g ) do not satisfy the renormalization group invariance condition.

As an additional possibility, we shall consider the case of quasi-polynomial kinetic functions, namely:

∑ = = q i g i i x a x v 1 ) (
, where the g i are in general real numbers. We thus impose the renormalization group invariance condition v * (x * )= v(x) and the outcome is that the following identity must be verified:

⎟ ⎟ ⎠ ⎞ ⎜ ⎜ ⎝ ⎛ = ∑ ∑ = - = q i g d g i d q i g i i i i x a x a 1 ) 3 ( 1 λ λ
Consequently, renormalization group invariance exists if and only if d+g i (d-3)=0 for every i=1,...,q. In turn, this implies that g 1 = g 2 =...= g q =d/(3-d), namely all exponents are equal and then the quasi-polynomial is actually a power-law. Otherwise, renormalization group invariance is not present for any quasi-polynomial not being a power-law.

We thus conclude that the power-law seems to be the only evident renormalization group invariant kinetic function in the single-variable case. Clearly, the same property

Generalization to the n-variable case

Let us now consider the generalized situation in which several metabolites A 1 ,...,A n react in a region of volume V and fractal dimension d. We denote by M 1 ,...,M n the initial number of molecules of each kind in that volume, the concentrations of the metabolites thus being x i =M i /V. As before, N i will describe the variation in the number of molecules of A i that are present in that region per unit of time due to the reactions taking place, and then the reaction velocities are v i =N i /V. Let us assume again that the system is self-similar after a scaling of factor λ, and denote with a star superscript the system parameters in the new scale. Now we have

V * = λ 3 V but M i * = λ d M i .
Then, for the concentrations we find

x i * =M i * /V * =λ d-3 x i .
As we did in the previous subsection, it is postulated that the dynamical invariance of the system in the new scale can be quantitatively expressed as:

* * * V M N V M N i i i i = (3)
for all i=1,...,n. In other words, now the system invariance is assumed to mean that the fraction of molecules of each species reacting per unit of time and per unit of volume remains invariant after the scaling. Thus, after (3) we again find that N i * = λ d+3 N i , as well as

v i * =N i * /V * = λ d v i .
Recalling that the rate-law v i (x 1 ,...,x n ) is actually described as a function of the concentrations, the scaling condition v i * = λ d v i must be understood as an identity between functions, namely

v i * (ξ 1 ,...,ξ n )= λ d v i (ξ 1 ,...,ξ n ) for any n-variable argument (ξ 1 ,...,ξ n ).
Again, our main goal is to characterize the renormalization group invariant rate-laws. In first place, we consider the case of power-law functions:

n i x v x x v ij g j n j i n i , , 1 , ) , , ( 1 0 1 … … = ∏ = =
After some algebra, the joint application of the renormalization group invariance criterion and the scaling condition leads to the following compatibility equations:

n i d d d g i n j ij , , 1 , 3 0 , 3 1 … = < < - = Γ = ∑ = (4)
For future convenience, we shall define

n Γ = = Γ ≡ Γ … 1 .
Obviously, equation ( 4) is a generalization of equation ( 2). As it was the case in the single-metabolite situation, equation ( 4) constitutes a solution that actually provides an exact invariance for all scales (i.e. for all values of λ). This result presents several interesting features additional to those found for a single variable. First of all, note that condition (4) admits an infinity of different compatible solutions for the kinetic orders g ij , including the possibility of simultaneous negative and positive values. Again, the average kinetic order increases when the fractal dimension grows, showing that the dimensional restriction of the medium should limit the ability of the metabolites to react. Once again, the most relevant conclusion is that, in this approximation, the power-law rate-law is an exact solution of the renormalization group condition for all values of λ (perfect self-similarity).

Due to the fact that equation ( 4) admits an infinity of solutions for the kinetic orders g ij , now the renormalization group renormalization group invariance extends also to the case of quasi-polynomial velocity functions:

n i x v x x v i ijk q j n k g k ij n i , , 1 , ) , , ( 1 1 0 1 … … = = ∑ ∏ = =
* n i n i x x v x x v … … =
, after some calculations we find the following compatibility conditions:

j i d d d g n k ijk , , 3 0 , 3 1 ∀ < < - = ∑ = (5)
This is obviously a generalization of (4). Thus, from (5) we see that every power-law making up every quasipolynomial velocity function v i must verify independently condition (4). Consequently, in particular ( 5) implies that all quasi-polynomial functions v i must be homogeneous of the same degree d/(3-d). Therefore, when modeling the reactions of n metabolites, the renormalization group invariance criterion is compatible with the use of the actual power-law models (GMA systems in general, including S-system models in particular).

For the sake of brevity, calculations involving alternative functional forms (such as n-variable rational functions generalizing Michaelis-Menten or Hill rate-laws) are not presented here. Such calculations are direct generalizations of those given for a single reacting metabolite. Moreover, the conclusions are entirely similar, as far as those functions are not compatible with renormalization group invariance. We thus conclude that powerlaw and quasipolynomial rate-laws appear as the natural renormalization group invariant solutions in the Zero-Order renormalization group approximation.

Alternative interpretation of the Zero-Order renormalization group exact solutions

Let us consider a general power-law based model of a metabolic pathway, namely a description in terms of a system of differential equations of GMA form, which for simplicity we can write as:

A c c e p t e d m a n u s c r i p t

13 n i x v x x v x x v dt dx n k g k q j ij q j n ij n i i ijk i i , , 1 , ) , , ( ) , , ( 1 1 0 1 1 1 … … … = = = = ∏ ∑ ∑ = = = (6)
Therefore in (6) we have that ) , , ( 1

n ij x x v
… is a power-law for all i, j. In addition, we assume that the renormalization group condition ( 4) is verified by every power-law function v ij or equivalently, that identity (5) holds for all functions v i conforming the righthand side of equations ( 6). Consider now that a scaling transformation is applied over the GMA system (6). Such transformation is defined in terms of a one-parameter Lie group transforming the system variables as: y i =μx i for all i=1,...,n. Since the x i are the system concentrations, actually such scaling is not a geometric one (as it was before, in the application of the renormalization group criterion) but a scaling in the metabolite concentrations. This explains that the group parameter now is termed μ instead of λ. In the new, transformed variables, the GMA system becomes:

∑ ∑ ∑ = Γ - = - - = = = = = i i i q j n ij q j n ij q j n ij i i y y v y y v x x v dt dx dt dy 1 1 1 1 1 1 1 1 1 ) , , ( ) , , ( ) , , ( … …
… μ μ μ μ μ μ
Precisely, the last step in the previous equation is the direct consequence of identities (4) and ( 5): note that if (5) holds for a power-law function v ij , then after some direct calculations we find that ) , , ( ) ,,( 1

1 1 1 n ij n ij y y v y y v … … Γ - - - = μ μ μ
. Therefore, in the new variables y i the GMA system remains the same with the only exception of a constant multiplicative factor Γ -1 μ . This means that if the GMA system velocities v i in (6) verify the renormalization group condition (5), then the phase-space trajectories of the rescaled flow are identical to the original ones, namely such trajectories are exactly the same. In fact, this can be shown in a simple way if we perform a time reparametrization of the form dτ = Γ -1 μ dt, where t is the initial time variable and τ is the new time. We thus see that:

n i y y v dt dy dt dy d dt d dy i q j n ij i i i , , 1 , ) , , ( 1 1 1 … … = = = ⋅ = ∑ = - Γ μ τ τ (7)
A c c e p t e d m a n u s c r i p t 14 As indicated, a time reparametrization does not modify the system trajectories, its only effect being a rescaling of the "speed" at which the phase state vector moves over such trajectories. Therefore, we see in ( 7) that the phase-space of the GMA system remains invariant after the scaling in the concentrations, and this happens if and only if the Zero-Order renormalization group invariance condition is verified. Mathematically, such condition in fact means that the GMA system is homogeneous, as anticipated, and this property accounts for the invariance found. Note that this property is independent of the actual value of parameter Γ, actually. Then, we have shown that the renormalization group invariance implies also phase-space invariance under scaling of concentrations and, in particular, this means that the GMA model remains valid over an infinite range of concentrations. In fact, this result accounts for the property, already mentioned in the Introduction, involving the reported validity of GMA and S-system models over several orders of magnitude in concentrations. In the Zero-Order approximation we obtain an idealized result in which the theoretical renormalization group range for concentrations is infinite, something that will be polished in the next section. We thus see how the renormalization group methods provide a direct link between self-similarity (geometric scaling) and the validity of the model over broad variations in the system variables (scaling in metabolite concentrations). Therefore the renormalization group approach leads to a unified perspective of both aspects of power-law modeling. self-similarity at all scales. This is clearly an idealization. It seems thus a good idea to proceed to a higher degree of approximation, in order to refine and validate different aspects of the results previously found. As it was done in the last section, now two cases (a single metabolite and more than one metabolite) are distinguished for the sake of clarity.

One-variable case

As in Subsection 2.1, we consider a single metabolite A reacting in a volume V and fractal dimension d. With the same notation, M is the initial number of molecules in that volume, and N is the increment per unit of time in the number of molecules of A in that region. As before, the velocity is v=N/V. After a geometric scaling of factor λ, we again have V * = λ 3 V as well as M * = λ d M. These are general features that remain unaltered. Now the refinement to be presented arises from the evaluation of the concentration x. In the Zero-Order renormalization group approach, it was x=M/V. However, M is the initial number of molecules in the region, namely such number ignores the variation due to the reaction taking place. In order to refine the evaluation of x, this time we shall consider the evolution due to the ongoing reaction during an infinitesimal time interval δ. At the initial time t, it is true that the concentration in the region is M/V. However at the time t+δ, such concentration will be M/V + δv. We can thus make use of an averaged value, and write in what follows that:

( ) ε δ δ + ≡ ⎟ ⎠ ⎞ ⎜ ⎝ ⎛ + = + = 1 2 1 2 V M M N V M V N V M x (8)

A c c e p t e d m a n u s c r i p t

16 where ε<<1 is a small dimensionless parameter. The smallness of ε arises from the fact that typically N<<M and, in addition, the time interval δ is infinitesimal and for practical purposes can be made arbitrarily small. It is worth recalling that the kind of approximation based on taking an average value at the central point t+δ/2 of the time interval is very common in numerical analysis, for instance in the domain of finite element methods [START_REF] Young | A survey of numerical mathematics[END_REF]. Notice also from (8) that in the limit case 0 → ε , we remain in the Zero-Order renormalization group approximation. As before, we now assume that the geometric and dynamic invariance of the system in the new scale can be quantitatively expressed as:

* * * V M N MV N =
Taking into account this invariance condition, we again obtain N * = λ d+3 N and v * =N * /V * =λ d v. Finally, we can use (8) in order to compute x * :

( ) ε λ λ δ 3 3 * * * * * 1 2 1 + = ⎟ ⎟ ⎠ ⎞ ⎜ ⎜ ⎝ ⎛ + = - V M M N V M x d
Let us now impose the one-variable renormalization group condition for velocities, namely

v * (x * )= v(x).
In what follows, we shall focus on the case of a power-law rate-law, v(x)=v 0 x g . Then, after some simplifications we arrive to the condition: ( ) (

)

g d g d g ε λ λ ε 3 ) 3 ( 1 1 + = + - + .
Note that this equation is identically satisfied when λ=1 for all values of g, d and ε. In addition, in the limit 0 → ε , we consistently retrieve the Zero-Order renormalization group identity (2), as expected. From the previous First-Order renormalization group relationship, by taking logarithms we obtain that:

⎟ ⎠ ⎞ ⎜ ⎝ ⎛ + + + - = ε λ ε λ 3 1 1 log log 1 ) 3 ( d d g (9)
A c c e p t e d m a n u s c r i p t 17 As before, the limit 0 → ε amounts to the Zero-Order renormalization group identity (2).

However, now some important differences arise. First of all, notice that now the Zero-Order divergence observed for d=3 is avoided. In fact, g can take a wide range of values for 3 ≤ d without a formal divergence. As an illustration of this result, Figure 1 9), namely that such kinetic order is dependent on λ. As explained before, an exact renormalization group invariant solution being independent of λ, as it was the case in the Zero-Order approach, reflects a perfect and exact invariance at all scales. However now this is not the case, as far as ( 9) is λ-dependent. This means that, in fact, a given power-law ratelaw now can be an exact solution of the First-Order renormalization group equations for a single value λ 0 , but not for all values of λ. In other words, in this approximation a powerlaw with a fixed kinetic order g does not comply to a perfect self-similarity at all scales. However, this is not the end of the story. Actually, power-law functions do present an

A c c e p t e d m a n u s c r i p t

18 approximate First-Order renormalization group invariance for a finite range of scales: the reason is that condition (9) depends on the small parameter ε. As indicated, if ε=0 we have an exact invariance for all scales (Zero-Order approximation). Thus, in the First-Order approach the power-law is not exactly renormalization group invariant for every λ, but it is invariant to a high degree of approximation over an interval of values of λ, which depends on the small parameter ε. This is the case because function 9) is close to zero in absolute value and varies slowly as λ changes. For instance, Figure 3 presents the relative deviation between the First-Order kinetic order (9) and the Zero-Order one (2). Let us denote both kinetic orders by g F and g Z , respectively. Then the relative deviation is defined as R=(g Z -g F )/g Z . Accordingly, Figure 3 depicts the variation of R for a fixed value of the fractal dimension d and for broad variations in the scaling parameter λ and in ε. As expected, the relative variation remains small.

⎟ ⎠ ⎞ ⎜ ⎝ ⎛ + + ε λ ε 3 1 1 log in (

Figure 3.

Therefore, it is worth recalling that the degree of approximate First-Order renormalization group invariance of a power-law is dependent on the precise value of λ, not being equal at all scales. As we have seen, now g is a function of (d, λ, ε). This means that the renormalization group invariant exponent g now is not one and the same for all scales.

Accordingly, a given power-law can only provide exact invariance for a given scaling, namely for a given value λ 0 of parameter λ, as well as a very approximate invariance for values of λ close to λ 0 (for constant values of d and ε). Emphasizing once again the approximate character of the renormalization group techniques, the previous conclusions seem to provide a justification about the prominent role of the power-law functions in the description of metabolic pathways in self-similar media: briefly speaking, it can be concluded that such functions approximate very precisely the property of self-similarity over a broad (but necessarily finite) range of scales.

A c c e p t e d m a n u s c r i p t 19

In what follows, let us analyze the generalization to several reacting metabolites.

Generalization to the n-variable case

As it was the case in Subsection 2.2, we again consider the more general situation in which several metabolites A 1 ,...,A n react in a region of volume V and fractal dimension d.

As before, M i is the initial number of molecules of A i in that volume for every i=1,...,n, respectively. Also, N i is the increment in the number of molecules of A i in that region per unit of time, and then v i =N i /V. Accordingly, if the system remains self-similar after a scaling of factor λ, we again have

V * = λ 3 V and M i * = λ d M i .
As in the previous arguments, it is postulated that the dynamical invariance of the system in the new scale can be quantitatively expressed as:

* * * V M N V M N i i i i =
for all i=1,...,n. After such relationship we again find N i * = λ d+3 N i , as well as

v i * = λ d v i .
Recall that this condition must be written in functional terms as

v i * (ξ 1 ,...,ξ n )= λ d v i (ξ 1 ,...,ξ n )
for any n-variable argument (ξ 1 ,...,ξ n ). To complete the preliminaries, we need an expression for the concentrations x i . Following a similar reasoning to that in Subsection 3.1, we introduce an infinitesimal time interval δ and write an averaged value for the concentrations:

( ) i i i i i i i i V M M N V M V N V M x ε δ δ + ≡ ⎟ ⎟ ⎠ ⎞ ⎜ ⎜ ⎝ ⎛ + = + = 1 2 1 2
Consequently, after the scaling the concentrations are transformed as:

i i d i i i i V M M N V M x ε λ λ δ 3 3 * * * * * 1 2 1 + = ⎟ ⎟ ⎠ ⎞ ⎜ ⎜ ⎝ ⎛ + = -
Once again, we wish to investigate the renormalization group invariance. As we know, the renormalization group criterion amounts to ) , , ( ) , , (

1 * * 1 * n i n i x x v x x v … …
= for all i=1,...,n.

Specifically, we shall analyze the case of a power-law functional dependence, i.e. we investigate the invariance of the following set of rate-laws:

n i x v x x v ij g j n j i n i , , 1 , ) , , ( 1 0 1 … … = ∏ = =
After some calculations based on the same procedures used in previous sections, eventually we arrive to the following set of compatibility conditions:

n i d g d ij n j j j , , 1 , 1 1 log log 1 3 1 3 … = = ⎥ ⎥ ⎦ ⎤ ⎢ ⎢ ⎣ ⎡ ⎟ ⎟ ⎠ ⎞ ⎜ ⎜ ⎝ ⎛ + + + - ∑ = ε λ ε λ (10)
Notice that these equations admit an infinity of possible solutions for the kinetic orders g ij , including in particular the possibility of solutions having negative components. In addition, in the limit 0 → i ε the Zero-Order renormalization group result (4) is retrieved. The nature of relationships ( 10) is more evident if we consider the particular case in which ε 1 = ...=ε n ≡ε. If this is the situation, equation (10) becomes simplified as:

n i d d g n j ij , , 1 , 1 1 log log 1 3 3 1 … = ⎟ ⎠ ⎞ ⎜ ⎝ ⎛ + + + - = Γ ≡ ∑ = ε λ ε λ (11)
It is evident in (11) that it is a direct generalization of (9). It is thus natural that the conclusions obtained in the one-variable First-Order renormalization group approach remain valid.

A c c e p t e d m a n u s c r i p t 21

As before, the limit 0 → ε amounts to the Zero-Order renormalization group identity (4).

However, now the Zero-Order divergence observed for d=3 is again avoided, and the kinetic orders can take a broad range of values also in the case d=3. The second relevant property of the First-Order renormalization group result ( 11) is that such equations are again dependent on λ, thus showing that power-law rate-laws can be exact solutions of the First-Order renormalization group equations for a given value λ 0 of the scaling parameter λ, but not for all values of λ. In addition, we now have an approximate renormalization group invariance for values of λ close to λ 0 , and such approximate invariance is due to the dependence on the small parameter ε. Then, as it was the case in the one variable First-Order renormalization group analysis, the power-law functional form is not exactly renormalization group invariant over a finite interval of λ values, but it is invariant to a high degree of approximation over such interval which depends on the small parameter ε. In addition, the dependence on λ implies that the First-Order renormalization group exponents g ij now are not the same for all scales. Accordingly, a given set of power-laws can only provide an exact invariance for a certain scaling parameter λ 0 as well as an approximate invariance for values of λ close to λ 0 (provided constant values of d and ε). Therefore, we again conclude that power-law rate-laws implement very precisely the property of selfsimilarity over a limited (not infinite) range of scales. As shown, in the limit 0 → ε the exact (infinite) invariance is retrieved. Therefore, the property of self-similarity over a limited range of scales is implemented with an accuracy which can be as precise as desired if ε takes a value sufficiently close to zero in absolute value. Obviously, following a reasoning similar to the one displayed in Subsection 2.2, analogous conclusions are valid for quasi-polynomial rate-laws. However, the calculations are omitted for the sake of conciseness.

Alternative interpretation of the First-Order renormalization group exact solutions

As it was shown in the Zero-Order renormalization group analysis, when applied to a general GMA model (6), conditions (4) and ( 5) amount to the fact that the GMA differential equations are homogeneous. This allowed an additional interpretation of the

A c c e p t e d m a n u s c r i p t

22 Zero-Order renormalization group condition as the GMA system invariance with respect to the scaling group y i =μx i , i=1,...,n. As indicated, since the x i are the system variables (which are not lengths), actually such scaling is not a geometric one (as it was before, in the application of the renormalization group criterion) but a scaling in the metabolite concentrations. In the First-Order renormalization group approximation, now conditions (4) are replaced by their generalization (10). As we have seen, such generalization depends on the small parameters ε i in such a way that (4) is retrieved in the limit 0 → i ε . Due to the smallness of parameters ε i it is clear that conditions (10) imply that the associated GMA models are almost homogeneous, the exact homogeneity arising only in the limit 0

→ i ε .
Accordingly, the First-Order renormalization group criterion leads to an approximate GMA system invariance with respect to the scaling group y i =μx i , i=1,...,n. In other words, we have shown that GMA models verifying (10), namely almost homogeneous GMA systems, remain valid over a broad (but this time finite) range of concentrations. According to the original renormalization group hypothesis this result provides a link for GMA systems between (i) an approximate geometric self-similarity of the system over a finite range of geometric scaling factors, and (ii) an approximate phase-space invariance under scaling of concentrations over a finite range of variation of such concentrations. Consequently, we retrieve the property, mentioned in the Introduction, involving the validity of GMA and Ssystem models over several orders of magnitude in concentrations. In the Zero-Order approximation we obtained an idealized result in which the theoretical renormalization group range for concentrations is infinite, something that now has been polished. Therefore the First-Order renormalization group approximation allows establishing a direct link between finite-range self-similarity (geometric scaling) and the validity of the model over broad variations of the system variables (scaling in metabolite concentrations).

Example: linear chain with feedback inhibition

In what follows a detailed example is presented for the sake of illustration. Specifically, the metabolic pathway considered is a linear chain with feedback control by

Γ - + Γ - = - = 2 2 1 2 2 1 1 2 0 1 1 x b x a dt dx x b x x a dt dx γ γ (12)
Recall that this Zero-Order model is necessarily homogeneous, as discussed in Section 2.

We see that the resulting GMA model ( 12) is also an S-system. Notice that, as usual in this case, the concentration x 0 of the initial precursor X 0 is considered to have a fixed value that can be independently modified by experimental means. The resulting system is thus twodimensional, a very convenient feature for illustrative purposes, as we shall see. In what follows, the values Γ = 2.5, d = 15/7 ≈ 2.14 and γ = 0.5 will be the ones employed. In addition, for the sake of clarity it will be very convenient to have the steady-state of the Zero-Order model in the coordinates (1,1). Accordingly, we set . 1

2 2 1 0 1 = = = = + Γ b a b x a γ
This is just a convenient choice of parametric values suitable for presentation purposes, but any other set of values would be equally good. The dynamics of the resulting system correspond to a stable steady-state, as depicted in Figure 5. Let us now define the First-Order system associated to this pathway. For this, we must take into account equation ( 10). For convenience, let us introduce the following definition:

2 , 1 , 0 , 1 1 log log 1 3 ) , , ( 3 = ⎟ ⎟ ⎠ ⎞ ⎜ ⎜ ⎝ ⎛ + + + - = j d d w j j j j ε λ ε λ ε λ
We then arrive to the First-Order system, which is given by:

2 1 1 2 0 / 2 / 1 2 / 1 / ) ) ( ( 2 0 1 1 w d w d w d w w d x x dt dx x x x a dt dx - = - = + Γ - + Γ γ γ (13)
As it can be seen, in the previous equations the kinetic order of X 0 is maintained as Γ+γ for clarity of illustration as in the Zero-Order case (recall that by definition it is 1

0 1 = + Γ γ x a
).

Note that the expression γ + Γ 0 1 x a takes the value 1 for convenience, in order to place the steady-state at point (1,1), as mentioned before. Such factor is explicitly written in Eq. ( 13) for clarity, because from condition (10) the value of the exponent of x 0 determines the value of the exponent of x 2 in the same monomial. Alternatively, the factor γ + Γ 0 1 x a can be replaced by its constant value 1. Note that in general, the rescaling is applied to both independent and dependent variables: however, independent variables are actually constant parameters and, in fact, they remain unaltered when the rescaling is applied to them. This can be verified easily, and x 0 in the present model is an example of it. In addition we shall set ε 0 = 0 and consistently w 0 = 3-d. Also, the values ε 1 = 10 -5 and ε 2 = -10 -7 are defined for the rest of the example.

A c c e p t e d m a n u s c r i p t 25

In first place, we shall consider the effect of geometric scalings over the First-Order system (13). According to (10), mathematically this amounts to the variation of parameter λ. Of course, this variation modifies the values of the kinetic orders in model ( 13). The effect of progressively higher scalings is displayed Figure 6. In particular, the phase plots are displayed for λ=2 (Figure 6a), λ= 50 (Figure 6b) and λ=100 (Figure 6c). In addition, the trajectories displayed in each plot are equivalent, namely for the sake of comparison the same initial conditions have been chosen for the trajectories plotted in Figure 5 and in Figures 6a, 6b and 6c. The presentation of this kind of graphical comparison is possible as far as the steady-state remains by construction located at (1,1), as anticipated. It is evident that even with a 100-fold geometric scaling (namely, over several orders of magnitude) the system trajectories remain essentially unaltered, in full agreement with the approximate invariance established for the First-Order renormalization group power-law models.

Figure 6 (Figures 6a, 6b, 6c).

In order to complete the example, we shall now turn to examine the effect of scalings in the metabolite concentrations in the First-Order system (13). Let us recall that such scalings are given by the one-parameter Lie group y i = μx i for all i=1,2,3. When applied to system (13), and taking into account that 1

0 1 = + Γ γ x a
, the outcome is: More precisely, the phase plots in Figure 7 are displayed for the values μ=2 (Figure 7a), μ= 50 (Figure 7b) and μ=100 (Figure 7c). In all cases, the geometric scaling parameter has the value λ=2. Note that, consistently, for every value of parameter μ now the steady-state is approximately centered at the coordinates (μ,μ). Again, even with a 100-fold concentration scaling it is clear that near-invariance is present as far as the stable steady-state topology of the solutions remains unaltered. The "straight-line" aspect of trajectories in Figures 7b and7c is due to the magnification efect produced on the phase-plot by this kind of transformation. Consequently, we verify also in the case of concentration scalings the approximate invariance established for the First-Order renormalization group power-law models over several orders of magnitude.

2 2 1 1 1 1 2 0 2 0 / 2 / 1 / 1 / 1 2 / 1 / 1 / ) ) ( ( 2 / ) ) ( ( 1
A c c e p t e d m a n u s c r i p t 27

Discussion

As shown throughout the previous section, power-law functions can be regarded as the most basic renormalization group invariant rate-laws. This has several implications that deserve some additional comments.

In first place, it is worth recalling that geometric invariance is actually present only over a limited (not infinite) range of scales. Precisely, the significance of power-laws is that they approximate such invariance over a limited range of scales, which is the actual property to be described. As we have seen, such feature is not present in other common rate-laws such as Michaelis-Menten, Hill, or their generalizations. However, the previous results are also useful in order to account for the ability of many functions to provide accurate predictions when employed in the construction of differential models: a possible explanation is that such functions are often relatively close to power-laws [START_REF] Sorribas | Cooperativity and saturation in biochemical networks: a saturable formalism using Taylor series approximations[END_REF]); accordingly, those functions could, in some situations, reproduce precisely the actual invariance properties of the system which is being analyzed.

It is also convenient to provide some additional considerations regarding the assumptions being made in the present work. This is necessary in order to put in perspective the present results with respect to previous approaches [START_REF] Kopelman | Rate processes on fractals: theory, simulations and experiments[END_REF][START_REF] Kopelman | Reaction kinetics in restricted spaces[END_REF][START_REF] Savageau | Influence of fractal kinetics on molecular recognition[END_REF][START_REF] Savageau | Michaelis-Menten mechanism reconsidered: implications of fractal kynetics[END_REF][START_REF] Savageau | Development of fractal kinetic theory for enzyme-catalysed reactions and implications for the design of biochemical pathways[END_REF]. In first place, it is worth noting that the existence of kinetic equilibrium is not assumed, as far as an ongoing reaction is considered in the previous calculations. As a consequence, the concentrations are in fact not taken as constant in the first-order analysis. The fact that the average kinetic order is a decreasing function of the fractal dimension d, as well as the precise functional form of this dependence, is one the most direct and novel consequences obtained from this point of view. Another difference with previous approaches is that usually the effects of a dimensionally limited diffusion are the ones analyzed, while here we have focused on purely reactional effects (namely excluding diffusion). In other words, our approach deals with kinetics taking place in a dimensionally limited environment, that is to say a region of fractal dimension d < 3, on which only reaction kinetics is considered. In the present work, the specific analysis of presented, which is to the author's knowledge new in the literature. In addition, it seems that focusing on reactional effects is natural in this context, as far as power-law models such as S-systems and GMA systems are purely reactional.

Another aspect of interest involves the practical use of GMA models and the evaluation of kinetic orders. As we have seen in the First-Order renormalization group approach, according to (10) the invariant GMA models should be almost homogeneous. In practice, it is worth recalling the difficulties in the construction of models (e.g. see [START_REF] Voit | Computational analysis of biochemical systems: a practical guide for biochemists and molecular biologists[END_REF]. In fact, even the fitting of a single power-law from data is today a controversial issue [START_REF] Kaitaniemi | Testing the allometric scaling laws[END_REF][START_REF] Packard | On the use of logarithmic transformations in allometric analysis[END_REF]. Due to the relevance of different uncertainties and sources of errors (experimental and numerical) in parameter estimation for model determination, specially when modeling in vivo processes, probably it is not to be expected that the homogeneity or almost homogeneity properties should be apparent in practice.

As recalled in the previous analysis, the renormalization group methodology is an approximation. In spite that its quantitative predictions are often not precise, renormalization group equations usually provide qualitative hints of great value in a workable way. This explains the successes of renormalization group theory in many different domains. In our case, the most important results obtained are those providing a link among three aspects of modeling: (i) the use of power-laws for the description of ratelaws; (ii) the geometric (and according to relationships (3) also physical) self-similarity under spatial scalings; and (iii) the self-similarity under concentration scalings (accounting for the reported GMA model validity over several orders of magnitude in metabolite concentrations). In particular, the last item is a consequence of the system invariance characterized which provides a mathematical background for the observed model validity First-order renormalization group kinetic orders. A typical behavior of the first-order renormalization group kinetic order (9) for some selected values. In this graph, the value λ=5 is the one selected for the scaling parameter. It can be seen that the kinetic order g is a decreasing function of the fractal dimension d.

Figure 2:

First-order renormalization group kinetic orders. Alternative representation of a typical behavior of the first-order renormalization group kinetic order (9) for some selected values. In this graph, the value d=2.3 of the fractal dimension is fixed, while λ and ε are allowed to vary.

Figure 3:

Relative deviation between the renormalization group kinetic orders. Behavior of the relative deviation R for the kinetic order g between the first and Zero-Order renormalization group approximations. R is displayed for a fixed value of the fractal dimension d=2.3 and for broad variations in the scaling parameter λ and in ε. As expected, the relative deviation remains small, in this case within limits of ±5%. Trajectories of the First-Order system (13) for different values of the scaling parameter λ. The phase plots are displayed for λ=2 (Figure 6a), λ= 50 (Figure 6b) and λ=100 (Figure 6c). The trajectories displayed in each plot are equivalent, namely for the sake of comparison the same initial conditions have been chosen for the trajectories plotted in Trajectories of the First-Order system ( 14) for different values of the concentration scaling parameter μ. The phase plots are displayed for μ=2 (Figure 7a), μ= 50 (Figure 7b) and μ=100 (Figure 7c). In all cases, the geometric scaling parameter has the value λ=2.

  independently, but both were intended to analyze what is invariant under the change of scale of observation: fractal for geometrical objects, and renormalization group for physical quantities. It is also important to stress that, except for very simple problems, the is generally an approximation. In other words,

  remarkable features accounting for the successes of power-law modeling. As we shall see, according to the renormalization group point of view the answer is positive, as far as power-laws are the critical solutions of the renormalization group transformation, namely power-law rate-laws are the renormalization group invariant solutions. Moreover, it is shown that these results also imply the system invariance under the group of concentration scalings, thus accounting for the reported power-law model accuracy over several orders of magnitude in metabolite concentrations. The work concludes in Section 4 with a discussion of the results found.

  group condition amounts to v * (x * )= v(x). This identity can be combined with the scaling condition v * (ξ)= λ d v(ξ) in order to explore the implications of renormalization group invariance in several cases.

  cases the well-known Michaelis-Menten (g=1) and Hill rate-laws.

  valid in the generalization to several reacting metabolites. The investigation of this issue is the aim of the next subsection.

Following

  i=1,...,n. Let us now analyze the implications of this condition for several functional dependences.

  if we now apply the renormalization group invariance condition,

  renormalization group approximation: exact solutionAs indicated in the Introduction, the renormalization group is essentially an approximation, as are the corresponding results. The Zero-Order renormalization group calculations just developed provide a good starting point in order to understand the significant role of power-law models in the implementation and description of selfsimilarity. Actually, the results are valid for all values of λ, therefore describing a perfect

  displays a typical behavior of the First-Order renormalization group kinetic order (9) for a fixed value of λ and varying values of d and ε.

Figure 1 .

 1 Figure 1.

  Figure 2.

Figure 2

 2 Figure 2 also illustrates another relevant property of the First-Order renormalization group result (9), namely that such kinetic order is dependent on λ. As explained before, an exact

  Figure 5.

  means of system (14) we can check the additional effect of scaling in concentrations, which is displayed in Figure7.

  Figure 7 (Figures 7a, 7b, 7c).

  allowed the use of the renormalization group methods in order to develop the combined perspective of geometric scalings and concentration scalings just

  FIGURE CAPTIONS

Figure 4 :

 4 Figure 4: Metabolic pathway for the example in Section 3.4. The example considers a linear chain with feedback control by inhibition. The end product X 2 acts acts as an allosteric effector causing inhibition of the first reaction in the sequence.

  Figure 6:

Figure 5

 5 Figure 5 and in Figures 6a, 6b and 6c.

  Figure 7:

  Fig 1

  Renormalization group techniques constitute a significant tool of physics. Actually, renormalization group methods are very diverse, ranging from quantum field theory to statistical mechanics. The last perspective, naturally associated to real space transformations, will be the one of interest in what follows. For instance, seeCreswick et al.

A c c e p t e d m a n u s c r i p t 4 1.2. Brief outline of the renormalization group methods

inhibition, a well-known mechanism[START_REF] Savageau | Biochemical systems analysis: a study of function and design in molecular biology[END_REF]. Such pathway is schematically displayed in Figure4.Figure 4.In terms of power-law modeling, and in agreement with equations (4) and (5), the Zero-Order model for the system is described by the following set of equations, in which for every index i the variable x i accounts for the concentration of the metabolite X i . Γ Γ

Acknowledgements

This work has been supported by the Spanish MEC, Project Ref. MTM2006-10053.

In addition, the author acknowledges two anonymous referees for their advice and fruitful comments.