
				 understanding metabolic activity from the underlying genotype is one of the most addressed problems in computational biology .
of particular interest is the issue of the identification of the reactions that are indispensable for an organism to survive , grow or perform a specific function in a given growth medium or , conversely , of the potentially most harmful knock - outs or enzymopathies .
several experimental protocols are able to assess the essentiality of gene products ( and hence of the corresponding metabolic reactions ) , ranging from individual knock - outs to transposon mutagenesis and rna interference [ 15 ]
. computational approaches on the other hand might provide important clues on the system - level organization by investigating genome - scale network reconstructions .
the functional modularity of metabolic networks suggests that topological aspects may provide a key to identify a class of essential pathways [ 6 , 7 ] .
however the metabolic genotype only constitutes the frame on the top of which the dynamic phenotype is built .
the essentiality of a metabolic pathway will in general depend on both structural considerations based on the network reconstruction from genomic information , and on the  model of metabolism  defined on it , for example , on the corresponding steady state fluxes . in   e.coli
, phenotypical essentiality of metabolic genes has been associated with a reduced allowed variability of the corresponding fluxes , suggesting that dynamically stiff reactions may constitute an evolutionarily robust backbone of metabolism conserved over different species . here
we attempt a more thorough integration of topological and dynamical views to obtain a more comprehensive insight into a metabolic network 's organization , efficiency , and ability to respond to perturbations .
we will first associate the essentiality of a reaction with a measure of its topological replaceability by enumerating the alternative paths from substrate to product along the network edges , with the rationale that from a purely structural viewpoint more replaceable reactions are less likely to be crucial nodes of the network .
then we will validate and compare the essentiality map thus obtained with the metabolic phenotype resulting from the definition of a general constraint - based model for metabolic flux prediction .
we shall see that dynamical and structural measures of essentiality may offer complementary views of a reaction network 's robustness .
we carry out our analysis on the metabolic network of the human red blood cell ( hrbc ) , one of the most studied complexes in systems biology , from the earliest mathematical models of single biochemical pathways [ 9 , 10 ] to the currently available genome - scale reconstructions .
the reason for this choice lies essentially in its limited size . on the one hand , it allows to compute reaction replaceabilities   exactly by a suitable modification of johnson 's algorithm for counting loops in a directed graph . on the other
, it allows for the efficient application of various sampling methods to the space of viable flux states [ 8 , 13 ] .
the latter is a vital ingredient to address many important properties of erythrocytes . indeed for some organisms under certain conditions
it is reasonable to assume that the metabolic activity is aimed at maximizing a subset of the metabolic reactions ( or a function of them ) associated with a certain biological function . in such cases
for example , e.   coli 's metabolism has been shown to maximize biomass production under evolutionary pressure , but after a genetic knockout it responds with a minimum rearrangement of fluxes . while the production of the cofactors ensuring the maintenance of osmotic balance and the release of oxygen may be argued to be their metabolic goal
information - rich directions in flux space must be retrieved by coupling the underlying constraints on fluxes with other types of analyses .
much understanding has indeed been obtained from the uniform sampling of feasible states [ 13 , 16 , 17 ] and by functional studies , like the computation of extreme pathways , of metabolic regulatory structures [ 19 , 20 ] and of metabolic pools .
these aspects combined make hrbcs a key benchmark for both theories of metabolism and computational tools .
it is worth noting that the detailed structural information we derive ( i.e. , the full map of alternative paths for each substrate / product pair ) can not be retrieved by other methods .
unluckily , computation times still prevent scaling the approach we employ up to networks larger than a few hundred nodes .
given a reaction network , we want to compute , for any pair of metabolites a and b that are , respectively , substrate and product in a reaction i ( this situation will be indicated by aib ) , the number ab( ) of alternative pathways , excluding reaction i , of length  allowing for the conversion a  b , see figure 1 .
the rationale is that a reaction performing a metabolite conversion aib for which ab( ) ( or , more properly , ab( ) ) is large will be more easily substituted , in case of an enzymopathy or a knockout , than one for which the above quantity is small .
finding paths connecting two points of a directed network is a long - studied problem in computer science .
the focus is usually on locating the shortest paths or the fastest way to find any path .
enumerating   all the distinct paths between two vertices is however a less confronted issue . in our case
it is crucial to avoid overcounting , for example , due to self - intersecting paths .
we will identify the substitutive paths using the following trick : for each pair ( a , b ) of metabolites such that aib , revert i fictitiously .
this results in a new graph where an auxiliary edge bia replaced the edge aib , see again figure 1
. counting the number of alternative reaction chains producing b from a then comes down to computing the number of directed cycles , that is , non self - intersecting directed closed paths along the edges of the new graph , passing through the fictitious edge bia . thanks to the limited size of the hrbc network it is possible to solve this enumeration problem
,  will denote here the number of reactions in the alternative pathway (  = 4 in figure 1 ) .
the space of viable fluxes will be defined through a constraint - based approach which relies on more general assumptions than flux - balance analysis ( fba , ) .
fba is the standard method to model steady - state reaction networks where mass balance constraints are imposed to every metabolite . for a reaction network with n reactions and m metabolites ,
let us denote by a and b , respectively , the m  n matrices of output and input stoichiometric coefficients .
the stoichiometric matrix is given by s = a  b. letting  = ( i)i=1 denote a vector of fluxes ( with properly chosen bounds i  i  i ) , the concentrations c = ( c)a=1 of metabolites vary in time according to c=s-u , where u = ( u)a=1 stands for the net cellular uptake of metabolite a ( u > 0 if a is a global output of metabolism , u < 0 if a is consumed by the organism , u = 0 if a is mass - balanced ) . assuming a steady state , the concentrations are constant in time ( i.e. , c=0 ) and vectors  satisfying s = u , or 
 ( 1)(ab)=u , 
							 represent flux configurations ensuring that each metabolite meets its production or consumption constraints at fixed concentrations .
as n is typically larger than m , the system is underdetermined and feasible flux states form a convex set of dimension n - rank(s ) embedded in the n - dimensional space of fluxes . in absence of a selection criterion that allows to pick one solution out of this set ( as e.g. , a maximum biomass principle ) , a uniform sampling of the solution space should be carried out .
when n is sufficiently small ( as for hrbcs ) , this can be achieved effectively , albeit at a considerable computational cost , by monte carlo methods [ 13 , 16 ] or by message - passing procedures . here
we will consider a different but related flux scheme based on von neumann 's ( vn ) model of reaction networks . in the vn framework ,
one fixes the environment through a small set of intakes on nutrients and defines a self - consistent flux problem where the network chooses , given a target growth rate , how much of the nutrients to use and which metabolites are globally produced .
the equations describing the vn model have been studied by statistical mechanics methods in [ 23 , 24 ] . for an intuitive derivation ,
note that the quantities a and b represent , respectively , the total output and the total input of each metabolite for a given flux vector . then a flux vector such that a  b , with some constant  > 0 , describes a network state where metabolites are being produced at a rate at least equal to  , since for each of them the total output is at least  times the total input . it is simple to see that as  increases the volume of such flux vectors shrinks continuously ( for  = 0 every flux vector is a solution ) . in particular , there exists a value  of  , representing the maximum metabolic production rate compatible with the stoichiometric constraints , above which no suitable flux vectors exist .
the presence of conserved metabolic pools   implies  = 1 , so that in metabolic networks optimal steady state fluxes correspond to the solutions of the solutions of ( 2 ) do not coincide with those of ( 1 ) even for u = 0 .
interestingly , a finite volume of ( optimal ) flux states turns out to satisfy the above constraints .
this trait is at odds with both the behavior of the solutions of ( 2 ) for a random reaction network ( where a single solution survives at  ) and with the optimization that is usually coupled to fba ( where typically a single flux state maximizes the objective function ) , and points to the robustness of metabolic phenotypes . for   e.coli , in particular , the solutions of ( 2 )
have been shown to reproduce both the large - scale organization of fluxes and the individual measured rates .
in addition , fluxes with the smallest solution - to - solution fluctuations , representing the most susceptible parts of the network , turn out to be strongly correlated with   e.coli's phenomenologically essential genes .
the main technical advantage in using the vn scheme lies in the fact that its solution space can be sampled uniformly at very modest computational costs even for genome - scale models .
the algorithm allowing for this , which has been recently applied to sample   e.coli's solution space , is detailed in the following section .
we consider the hrbc metabolic network studied in , a map of which is shown in figure 2 ; table 1 lists reactions and the corresponding abbreviations .
the network comprises three main pathways , namely , glycolysis ( reactions 113 ) , the pentose phosphate ( pp ) pathway ( 1421 ) and the adenosine metabolism , with a total of m = 39 metabolites linked by n = 59 reactions : 49 internal reactions ( 34 of which come from the splitting of 17 reversible processes ) , 3 auxiliary fluxes to maintain the osmotic equilibrium and the redox state of the cell ( atpase , nadhase , nadphase ) and 7 uptake reactions to guarantee the intake of the necessary nutrients ( glu , ade , ado , ino ) , and of the cytosol elements ( h2o , h , pi ) .
the forward and backward parts of reversible reactions are treated separately throughout this study , both in the structural and in the flux analysis .
structural vulnerabilities are identified by analyzing the loop structure of a modified metabolic reaction network , created from the original one by inverting  in turn  the direction of the single reaction for which we want to compute the replaceability , as explained in figure 1 .
the fastest known exact algorithm ( for the worst case scenario ) of this cycle enumeration problem for a directed graph was introduced by johnson .
we shall now shortly describe its key ideas , referring to   for a pseudocode .
given a directed graph with n vertices and e edges , the algorithm is designed to build non self - intersecting paths from a root vertex r to itself , loading them onto stacks .
the main ingredients allowing for an optimal exploration of the graph are ( a ) a smart choice of the root vertex , and ( b ) an efficient method to avoid duplicating cycles and repeating searches on the same portions of the graph . to achieve this ,
vertices are initially ordered in a lexicographic sequence , and the algorithm only selects as roots those nodes that are the 
vertex ( in the initial ordering ) of at least one cycle . the algorithm described in   guarantees to find such vertices in o(n + e ) operations .
moreover , to avoid self - intersections , each time a node is loaded onto a stack it is also given a  blocked 
it was proven by johnson that if a vertex v stays blocked as long as every path from v to the root vertex r intersects the current path at a vertex other than r , the algorithm outputs all cycles exactly once . by sufficiently delaying the unblocking of each of these vertices and by keeping track of the portions of the graph that have been searched holding the current stack , the maximum time that can elapse between two consecutive cycle outputs can be reduced to o(n + e ) .
the same holds for the time window before the first cycle is delivered and for the one after the output of the last cycle .
hence , the total time needed to list the , say , c cycles of the graph is o((n + e)(c + 1 ) ) . in our case
, each fictitious reaction reversal generates a new graph , so that computing the complete substitutability map for a network of n reactions requires a time of the order o(n(n + e)(c + 1 ) ) . for practical reasons
, we perform this analysis on the bipartite metabolic network ( as in figure 1 ) rather than the reduced network of figure 2 .
this implies that in our case n = n + m. one can in principle consider different measures of replaceability of a metabolic conversion aib .
the quantity ab = ab( ) , counting the total number of paths alternative to i from a to b of any length , is perhaps the most obvious option .
taking into account the fact that , typically , longer detours can be less convenient than shorter ones from an energetic viewpoint one could instead consider -weighted functions like ab = exp  ( )ab( ) , with the caveat that shorter pathways might require more atp than longer ones
we limit ourselves to identifying three key reaction groups that are   independent of the replaceability measure used :    the group of reactions such that each substrate - product pair involved in them can be substituted ( this is putatively the part of the network that is most robust to enzymopathies ) ;   the group of reactions that can not be substituted , corresponding to the most harmful enzymopathies ;    the group of reversible reactions that are only replaceable in one direction , corresponding to the situation in which a conversion ab can only be substituted in one direction in case of a knockout .  
					
optimal flux vectors , that is , solutions of ( 2 ) , are computed by the algorithm introduced in   based on .
the idea is to modify fluxes iteratively until all inequalities in ( 2 ) are satisfied .
 ( with  = 1 in our case ) define ( ) = a  b and let ( ) denote the rows of ( ) , for a  { 1 ,  , m}. let also , for each iteration step t ,
(t ) be the flux vector at step t and 
 ( 3)m(t)=arg  min  a  a()(t ) , 
							 at each t , the algorithm runs as follows .
if ( )  (t ) < 0 , update fluxes according to 
 ( 4)i(t+1)=max  { 0,i(t)+im(t)( ) } , 
							 and iterate in t. else , if ( )  (t )  0 stop , that is , (t ) is a solution .
convergence to a solution is rigorously ensured for all 0   <  , and  can be approximated with the desired resolution by iterating the above process for increasing values of  . to guarantee that solutions are well defined
one can either resort to setting fixed upper bounds on i 's or , as we do , impose a linear constraint of the form ii(t ) = n on the solutions ( this is equivalent to singling out one flux as the reference unit for the other fluxes ) .
different initial points generate trajectories to different solutions at  and the sampling of the solution space thus obtained turns out to be uniform .
contrary to fba , the solution space of vn 's model is generically not a polytope .
indeed much useful information can be retrieved from its shape . as a means to characterize it we employ the average overlap between different optimal flux vectors , defined as follows .
let  and  denote two distinct solution vectors of ( 2 ) and , for each flux i , let 
 ( 5)q(i)=2iii2+i2 . 
							
 and  , equals 1 if flux i takes on the same value in solutions  and  and decreases as the values differ more and more . averaging q(i ) over different pairs of solutions provides a measure of the allowed variability of flux i ( smaller variability corresponds to larger average overlap ) , complementary to the standard deviation of the resulting flux distribution .
the complexity of the solution space can then be roughly understood by distinguishing narrower directions with larger overlap or less variable fluxes from broader ones .
it is reasonable to think that a cell will be more sensitive to perturbations ( e.g. , knockouts ) of fluxes with larger overlap . analyzing the susceptibility of the solution space to perturbations along the directions identified by different fluxes then allows to extract a list of the potentially more deleterious perturbations , in analogy with previous work on   e.coli . 
					 in order to test the hrbc network against enzymopathies ,
one can first employ a structural criterion : the knockout of a metabolic conversion a  b that is less easily  substituted  is more likely to be deleterious for the cell than the knockout of a highly replaceable conversion .
as said above , we concentrate here on a coarse - grained view of replaceability based on classifying reactions into the groups ( a ) , ( b ) , and ( c ) defined above , with groups ( b ) and ( c ) containing potentially essential reactions .
the second criterion is based on fluxes : fluxes with smaller allowed variability ( i.e. , larger overlap ) in the healthy cell are more likely to be essential links of the network than fluxes whose value can be changed over a larger range without losing optimality .
as is to be expected , the essentiality maps produced in these ways have a large degree of similarity , and reactions in the group ( b ) discussed above coincide with the physiologically most critical parts of hrbc 's metabolism . the simplest way to simulate an enzymopathy on flux
deficiencies can be partial , that is , of a smaller degree , the closer i is to the upper limit of the allowed range in the healthy cell , or total if i=0 .
such constraints cause in principle a modification of the solution space along the direction i which in turn cascades on the entire volume , modifying the optimal states of the metabolic network .
( for the sake of simplicity we exclude the highly replaceable currency exchange fluxes from this discussion . )
the most replaceable core of the network lies in the pp pathway ( reactions 1721 ) , which constitutes the main source of nadph , the key metabolite that in erythrocytes limits the accumulation of peroxides protecting the cell from hemolysis .
the high reliability coming with replaceability partly explains the reason why this group of reactions plays a central role not just as an auxiliary pathway for glycolysis , see the following analysis of fluxes .
unreplaceable reactions are instead lined up along glycolysis ( numbers 1,6,813 ) , in the bridge between glycolysis and the pp - pathway ( 14 and 16 ) or in auxiliary modules ( 22 , 27 , 29 ; the adeamp conversion in 25 is also not replaceable being directly linked to the ade uptake ) .
the physiologically most deleterious knockouts ( hk , pk , and g6pdh ) all belong to this group .
for instance , deficiency in the level of g6pdh is the basis of different types of hemolytic anemias , including favism , and is also linked to malaria resistance .
finally , there is a group of reversible reactions ( numbers 2 , 4 , 7 , 15 , 23 , 26 ) that can be replaced only in one direction .
note however that the last three of these could still be replaced in case of an enzymopathy if a proper medium is selected .
for instance , if reaction 15 is removed , it could be substituted by an alternative chain of reactions provided 6pgc is externally supplied .
this is instead not possible for reaction 4 and possibly 23 ( depending on the directionality of reaction 26 ) , as a knockout in these cases would necessarily result in a net production of fdp and r1p .
the flux distribution corresponding to optimal states in the healthy and enzyme deficient hrbc are displayed reaction by reaction in figure 4 , obtained by sampling 10000 solutions of ( 2 ) , while a pictorial representation of the optimal flux states is given in figure 5 .   for the healthy cell ( black line in figure 4 and top left panel in figure 5 ) the large flux backbone is formed by the second part of glycolysis ( crucial for atp , nadh , and 23dpg production ) and the pp pathway ( nadph production ) .
the latter gives a substantial contribution to the former , not just as salvage way .
in addition to glu , which is the fundamental substrate for hrbcs , the ino uptake plays an important role as an alternative way to the pp pathway .
it is worth stressing that these solutions imply a net production of 23dpg , the crucial regulator for oxygen release , which is obtained without any imposed constraint .
this picture is strongly reminiscent of the first eigenpathway obtained by extreme pathways analysis in , though the thermodynamic constraints and production requirements used in , including one on 23dpg , are more strict than the self - consistent analysis presented here . comparing the distributions with fba studies on the same system , one notices instead a general rearrangement of fluxes in the network apart from glycolysis .
a close inspection reveals that such a rearrangement is mostly quantitative , as preferred reaction directions are generically preserved , the noteworthy exception being the rpi flux , that in the vn solution strengthens the pp pathway with respect to the fba solution .
this scenario is not surprising in view of the basic difference between fba and the vn approach . it should be kept in mind however that the a priori constraints on flux variability are quite more strict in fba than they are in the vn model , and the flux distribution appear to be particularly sensitive to the assumed upper and lower bounds for the fluxes . in figure 6
one sees that the large overlap backbone ( signaling dynamically stiff fluxes ) coincides to a large degree with the structurally most vulnerable parts of the network .
note that the overlap of reactions 2 , 4 , 15 and 26 is larger in the direction that can not be replaced , further pointing to a higher susceptibility , and that currency reactions ( 3135 ) belong to the most constrained part of the network .
revealingly , however , topological and dynamical characterizations prove to be complementary in some cases .
this is seen , for example , from reaction 3 , which is flux - constrained but also highly replaceable , so that the damage due to removal is limited even in presence of a small allowed dynamical range .
( a similar picture holds for reaction 23 . ) to conclude , we remind that in our framework uptake fluxes are optimized variables not fixed by boundary conditions . in the optimal state five of the uptakes have a limited allowed variability , implying rather severe constraints on the cell 's environment .
we have simulated the most studied enzymopathies by constraining the flux of the corresponding reaction .
indeed , we have observed appreciable differences in relevant cellular functions compared to the nondeficient case only under full enzyme deficiencies , as also observed in   within a standard fba optimization approach . even under the most serious enzyme deficiencies the network appears to be able to maintain the production of atp , nadh and nadph almost constant , see also .
as shown in figure 4 , the alterations in the flux distributions are not particularly striking and indeed we do not observe global flux rearrangements on the network 's scale .
the g6pdh enzymopathy appears to only cause local changes , confirming the structural predictions , the overlap calculations and also in agreement with clinical observations .
the synoptic analysis of figure 5 shows that in general the response to the perturbation consisted in a drop of the glu uptake , and in a reduction of the glycolytic flux , while the rapoport - leubering shunt ( reactions 8 - 9 ) for the production of 23dpg remains particularly stable , as does the adenosine metabolism . for the glycolytic deficiencies pk and
hk we further observe an increase of the ino uptake to sustain the pp pathway and allow for the second part of glycolysis , and with it the production of atp and nadh , to take place .
detailed flux configurations corresponding to the next most severe enzymopathies ( hk , en , pgk and pgm ) are available from http://chimera.roma1.infn.it/sysbio .
in this work we compare two robustness measures for biochemical networks , one based on structural properties ( the reaction replaceability ) , the other based on dynamical stiffness ( the overlaps ) .
the former can be exactly assessed by enumerating the alternative paths joining substrates and products of a given reaction in a network .
the latter depends on both the network topology and the model defined on it . within vn 's frame
, we found that unreplaceable reactions mostly correspond to processes with a smaller allowed flux variability .
in such directions , reaction removals as well as constraints on the fluxes are expected to be generically harmful .
reactions with limited ( but non zero ) replaceability tend to have instead smaller overlap , so that while the reaction is difficult to substitute still its flux can be largely adjusted . in an evolutionary perspective
, the former pathways appear as  frozen  , and perturbations at these nodes will require large - scale flux rearrangement , while a mutation affecting the latter group may be neutral and could be preserved across generations .
these , albeit structurally robust , are dynamically constrained and should be considered as essential pathways of the metabolism as well . integrating dynamical and structural characterizations may thus provide a rather complete picture of the emerging network robustness .
the fact that topological and dynamical essentiality may not coincide could also prove to be important in view of the present challenges to understand the dynamical basis of topological modularity .
extended flux state sampling was achieved here by an algorithm that is easily scalable to larger networks , see .
exhaustive structural analysis instead was made possible by the small size of hrbc 's metabolic network , which has served here as a model system to test basic concepts and algorithms .
the use of the same procedure on a larger network , such as e.   coli 's , is likely to be prevented by cpu time growth .
however , message - passing algorithms , designed specifically to solve combinatorial optimization or counting problems  albeit approximately  on graphical models , may be a suitable replacement .
 characterizing the capabilities , key dependencies , and response to perturbations of genome - scale metabolic networks is a basic problem with important applications . 
 a key question concerns the identification of the potentially 
 most harmful reaction knockouts . 
 the integration of combinatorial methods with sampling techniques to explore the space of viable flux states may provide crucial insights on this issue . 
 we assess the replaceability of every metabolic conversion in the human red blood cell by enumerating the alternative paths from substrate to product , obtaining a complete map of 
  he potential damage of single enzymopathies . sampling the space of optimal steady state fluxes in 
 the healthy and in the mutated cell reveals both correlations and complementarity between topologic and dynamical aspects . 
1. Introduction
2. Approach
3. Methods
4. Results
5. Final Remarks