The effective deficiency of biochemical networks

The deficiency of a (bio)chemical reaction network can be conceptually interpreted as a measure of its ability to support exotic dynamical behavior and/or multistationarity. The classical definition of deficiency relates to the capacity of a network to permit variations of the complex formation rate vector at steady state, irrespective of the network kinetics. However, the deficiency is by definition completely insensitive to the fine details of the directionality of reactions as well as bounds on reaction fluxes. While the classical definition of deficiency can be readily applied in the analysis of unconstrained, weakly reversible networks, it only provides an upper bound in the cases where relevant constraints on reaction fluxes are imposed. Here we propose the concept of effective deficiency, which provides a more accurate assessment of the network’s capacity to permit steady state variations at the complex level for constrained networks of any reversibility patterns. The effective deficiency relies on the concept of nonstoichiometric balanced complexes, which we have already shown to be present in real-world biochemical networks operating under flux constraints. Our results demonstrate that the effective deficiency of real-world biochemical networks is smaller than the classical deficiency, indicating the effects of reaction directionality and flux bounds on the variation of the complex formation rate vector at steady state.


The effective deficiency of biochemical networks
Damoun Langary 1,2 , Anika Küken 1 & Zoran Nikoloski 1,2* The deficiency of a (bio)chemical reaction network can be conceptually interpreted as a measure of its ability to support exotic dynamical behavior and/or multistationarity.The classical definition of deficiency relates to the capacity of a network to permit variations of the complex formation rate vector at steady state, irrespective of the network kinetics.However, the deficiency is by definition completely insensitive to the fine details of the directionality of reactions as well as bounds on reaction fluxes.While the classical definition of deficiency can be readily applied in the analysis of unconstrained, weakly reversible networks, it only provides an upper bound in the cases where relevant constraints on reaction fluxes are imposed.Here we propose the concept of effective deficiency, which provides a more accurate assessment of the network's capacity to permit steady state variations at the complex level for constrained networks of any reversibility patterns.The effective deficiency relies on the concept of nonstoichiometric balanced complexes, which we have already shown to be present in real-world biochemical networks operating under flux constraints.Our results demonstrate that the effective deficiency of real-world biochemical networks is smaller than the classical deficiency, indicating the effects of reaction directionality and flux bounds on the variation of the complex formation rate vector at steady state.
Chemical reaction network theory (CRNT) has made seminal contributions to establishing necessary and/or sufficient conditions that a (bio)chemical network exhibits particular properties, such as: presence/absence of multiple steady state, stability of steady states, and robustness of steady-state concentrations 1 .The deficiency of a network is of central importance in this theory, and many steady-state concentration properties are guaranteed or precluded for networks of particular deficiency [2][3][4] .These results often hold for (bio)chemical networks whose reactions are endowed with (generalized) mass action kinetics.However, a notable feature of this theory is that it does not impose any physico-chemical constraints that real-world networks obey.As a result, the notion of deficiency is not concerned with the (ir)reversibility of the considered biochemical reactions that may arise in practice, and the extent to which this may affect the properties of the corresponding steady states the network supports.
In contrast to CRNT, which almost exclusively deals with concentration-based properties, the constraintbased modeling framework 5,6 , imposes physico-chemical constraints as lower and/or upper bounds on reaction rates (i.e.fluxes) in making predictions about macro-level phenotypes, such as growth or yield of a chemical product of interest.The constraint-based modeling framework has found numerous applications, largely due to its capacity to model large-scale networks by building on approaches for convex optimization.
A natural question then arises of how to integrate any imposed constraints on reaction fluxes into the notion of deficiency, and if the consideration of such constraints leads to smaller values such that one may expand the usability of the classical results from CRNT 2,3 .Here, we rely on the recently introduced classification of balanced complexes 7 , in particular the class of nonstoichiometric balanced complexes, to define the notion of effective deficiency of a network.Nonstoichiometric balanced complexes have been shown to arise as a combined result of several factors, including the algebraic and graphical structure of the network as well as operational bounds on reaction kinetics 7 .A remarkable feature of the effective deficiency is that, like that of the classical, so-called structural, it is defined and can be calculated for networks of arbitrary kinetics.
Our theoretical results show that decrease in the structural deficiency of a network, as quantified by the effective deficiency, is due to the presence of reactions whose fluxes are fixated at their lower bounds-i.e. they are robust irrespective of the environment.In other words, given the irreversibility and flux bound constraints, if the network contains nonstoichiometric balanced complexes, some reactions are guaranteed to exhibit absolute flux robustness.It is then not surprising that the effective deficiency captures the reduced capacity for variation of the complex formation rate at steady state.We conjecture that the effective deficiency, like the structural deficiency,

Background on CRNT
, each of which is a multiset of species C j ∈ N S , and a set of r reac- tions,R = R q r q=1 ⊂ C × C , which symbolize the potential conversion of complexes into each other in the network 2,3,8,9 .
We denote the standard basis in R n by e j n j=1 , where is an n− vector with a unit value at the jth entry and zero entries elsewhere.Let us next assume some arbitrary ordering on the sets of species Given this order- ing, any complex C j ∈ C can be associated with a vector e j ∈ R n indexing its position in the ordered set, and at the same time with a unique vector y j ∈ R m , which represents its species content.This defines the following mapping Consequently, any given CRN is associated with a stoichiometric map, defined by the matrix Similarly, the ordering associates each reaction R q : C p′ ⇀ C p , converting complex C p′ to C p , with a vector a q = e p − e p′ ∈ R n that represents this conversion in the complex space, and a vec- tor n q = y p − y p′ ∈ R m , which represents this conversion in the species space.Thus, the CRN is at the same time associated with the complex-reaction incidence matrix , known as stoichiometry matrix.It follows immediately from the definition of the stoichiometric map that N = YA.The column span of N is called the stoichiometric subspace, the dimension of which,s , is termed rank of the CRN, that is,s = rank(N).
Assuming the reaction network is endowed with some kinetics, at any state of the system, denoted by the vector of concentrations x ∈ R m ≥0 , the dynamics is governed by the following ODE where v(x) ∈ R r returns the reaction rates determined by system kinetics as a generally nonlinear function of the state vector x.
Blocked reactions and the steady state flux space.The vector v ∈ R r , referred to as a flux distribution for the network G = (S, C, R) , is said to be at steady state, if In addition, the flux vector v is often assumed to be bounded by box constraints as follows where ≤ operates element-wise, and v l , v u ∈ R r specify the individual lower-and upper bounds on flux through reactions, respectively.By convention,v u is strictly positive for all reactions , v u > 0 .Going forward, we assume a network denoted G also encodes information about upper and lower bounds on reaction fluxes.4).The set of all feasible flux distributions is an intercept of some hyperplanes and halfspaces, which defines a convex polyhedron, referred to as the steady state flux space, denoted The set of reversible reactions R rev ⊆ R is defined by negative entries of v l , that is where v l,R denotes the entry of v l corresponding to reaction R .The set of irreversible reactions R irr ⊆ R is then defined as the complement of R rev , that is,R irr = R \ R rev .Hence, irreversible reactions are associated with v l,R ≥ 0 .By contrast to irreversible reactions, a reversible reaction R q : C p′ ⇋ C p portrays the mutual conversion of complexes C p′ and C p .Similarly, let R zl and R nl denote the sets of reactions with zero and nonzero lower bounds, respectively.It follows from the definition that,R zl ∩ R nl = ∅; R zl ∪ R nl = R; and R zl ⊆ R irr .Given the arbitrariness of index- ing, one can always order the reactions in R such that elements of R zl and R irr precede those of R nl and R rev , respectively.As a result, the incidence matrix A can be block-partitioned as either The network is said to be operating under a canonical flux regime, if R zl = R irr .It is said to be operating under an unbounded flux regime, if v u = ∞ and v l,R = −∞, ∀R ∈ R rev .Under a canonical and unbounded flux regime, the set of feasible flux distributions forms a polyhedral convex cone, referred to as the steady state flux cone.The flux regime is called bounded, if it is not unbounded.
For any reaction R ∈ R , we say R is a blocked reaction in G , if for all flux distributions v ∈ F (G) , the flux through R is zero, namely , v R ≡ 0 .This is a recurring phenomenon in metabolic networks, in particular when restrictive flux bounds and/or optimization of particular objectives are imposed.As long as we only concern ourselves with steady state analysis, all blocked reactions can be safely removed from the network.
More generally, we say a reaction R ∈ R is fixated (at some flux value f ), if for all flux distributions v ∈ F (G) , the flux going through R is unchanged, namely, v R = f , ∀v ∈ F (G) .Obviously, the set of blocked reactions is contained in the set of fixated reactions.
Linkage structure and deficiency.Two complexes C, C′ ∈ C are said to be directly linked, , each of which is directly linked to the immediate preceding and succeeding ele- ments.The equivalence relation ∼ partitions C into a family of ℓ equivalence classes {L l } ℓ l=1 , called the linkage classes of the network.
For any two complexes C, C′ ∈ C , we say C directly converts to C′ , denoted C → C′ , if (C, C′) ∈ R or (C′, C) ∈ R rev .We say C converts to C′ , denoted C ⇒ C′ , if there exist a sequence of complexes A terminal strong linkage class is a strong linkage class , no complex of which converts to any complex in C\ .We denote the number of terminal strong linkage classes by ⊔ .It is trivial to show that, in general,ℓ s ≥ ⊔ ≥ ℓ .A reaction network is said to be weakly reversible, If C ∼ C′ implies C ≈ C′ , that is, any two linked complexes convert to each other.In graph-theoretic terms, it means each component of the CRN graph-defined by complexes as vertices and reactions as directed edges -is strongly connected.For weakly reversible networks, every linkage class is a terminal strong linkage class, hence ℓ s = ⊔ = ℓ.
The deficiency of a reaction network is the nonnegative integer, denoted δ , defined by 2 Alternatively, it can be defined as 10,11 As Eq. ( 6) demonstrates, the deficiency of a network has to do with the image of stoichiometric flux modes under linear mapping A .In this light, for any network G , we define the deficiency space D(G) as the set of all feasible complex formation rate vectors, given as follows Any nonzero vector a ∈ D(G) points at a nonzero complex conversion, i.e.net production or consumption of some complexes at steady state, which is unobserved on the species level, because all species concentrations are constant.Thus, one may interpret the dimension of the deficiency space as a measure of the capacity of the network to permit unobserved complex conversions at steady state.This is why the notion of deficiency stands at the center of several seminal results in CRNT, which determine whether special classes of networks may exhibit multistationarity and/or exotic dynamical behavior 3 .
Given the definition of A , it is not surprising that the linkage structure of the network is closely interlinked with the left nullspace of A .Let G be a closed network comprising ℓ linkage classes {L l } ℓ l=1 .An obvious choice of a basis for the left nullspace of A , namely kerA T , is the set u (l) ℓ l=1 , where u (l) = j:C j ∈L l e j , 1 ≤ l ≤ ℓ .Therefore, the column span of the complex -linkage class incidence matrix Vol:.( 1234567890 www.nature.com/scientificreports/coincides with the left nullspace of A .In a similar fashion, we can define the complex -strong linkage class inci- dence matrix as follows Since every linkage class is a disjoint union of some strong linkage classes, it follows that for any arbitrary network, im(U) ⊆ im(U s ) , with equality holding only for weakly reversible networks.

Mathematics of balanced complexes
Let us next give a formal definition of a balanced complex.For a given matrix X , let X :i and X j: denote the ith column and jth row of X , respectively.A complex C j ∈ C is referred to as a balanced complex (BC) for network G , if A j: v ≡ 0 for all v ∈ F (G) .Given the properties of the incidence matrix A , this definition complies with the notion that the algebraic sum of fluxes entering complex C j must be equal to the sum of fluxes leaving it for all feasible distributions v ∈ F (G) , which means this complex has a zero complex formation rate at all steady states.
The contents of this section build heavily upon earlier works in 12 and 7 .We refer the readers interested in a more detailed discussion to these studies for more information.

Factorizations of balanced complexes.
As was shown in 7 , a complex C j ∈ C is a BC, if and only if where variables The set of all balanced complexes in G , denoted B = B(G) , can then be defined as A large subset of these balanced complexes can be characterized as follows In a canonical flux regime, the above two sets coincide 7 , i.e.B = B 1 .Note that in that case,A nl = A rev and Furthermore, in a canonical flux regime, if the network is void of any blocked reactions, the set of balanced complexes reduces to a particular subset of which is given as follows The equalities in Eqs.(10), (11), ( 12) and ( 13), referred to as factorizations, explain the formation of a balanced complex C j ∈ C as a combined effect of a number of underlying factors, namely, the stoichiometry ( Y ), linkage structure ( U ), irreversibility patters ( R zl ⊆ R irr ) and flux bounds ( v l , v u ).

Nonstoichiometric balanced complexes. It is worth noting that
A complex C j ∈ C is called a strictly stoichiometric BC, if C j ∈ B 3 ; the equality in ( 13) is referred to as a strictly stoichiometric factorization for complex in ( 12) is a stoichiometric factorization for complex C j .By contrast, a complex The notion stoichiometric arises from the fact that due to the factorizations given in ( 12) and ( 13), the balancing property for such complexes relies heavily on the stoichiometric structure of the network, and not on other imposed constraints, such as: irreversibility and flux bounds.The notions are borrowed from 7 , to which we refer the interested reader for more details.
While all nonstoichiometric BCs have an implicit factorization of the form (10), we make a distinction between those which also have an explicit factorization of the form (11), and those which do not.A complex that is, it has a factorization of the form (11), but no stoichiometric factorization.A complex C j ∈ C is called a type-II nonstoichiometric BC, if C j ∈ B\B 1 , that is, it has a factorization of the form (10), but none of the form (11).A complex C j ∈ C is unbalanced, if it has no factorization of the form (10), that is , C j ∈ C\B.
The emergence of nonstoichiometric BCs is an intriguing phenomenon in (bio)chemical reaction networks.As will become clear in the next section, it reflects on some key structural properties of the reaction network.Furthermore, it also reflects on the reduced capacity of a network to permit variations of the complex formation rate vector at steady state.

Properties of nonstoichiometric BCs
We hereby aim to seek conditions under which nonstoichiometric BCs may emerge in a network, and to see what implications their existence may have on structural properties of a reaction network.We begin by analyzing type-I nonstoichiometric BCs, which is the only viable type in canonical flux regimes.However, it shall be noted that they may also emerge in non-canonical flux regimes.The following statement was proven in 7 .
Proposition 4.1 Given a network G , if the set B 1 \ B 2 is nonempty, then G contains at least two irreversible reac- tions which are blocked at steady state.
The proof follows from the fact that nonzero (positive) entries in vectors A zl T θ 1 and A zl T θ 2 must correspond to blocked reactions in R zl .Existence of C j ∈ B 1 \ B 2 requires that the vectors A zl T θ 1 and A zl T θ 2 are both nonzero and also not collinear.Proposition 4.1 shows how the presence of a nonstoichiometric BC reflects qualitatively on the steady state properties.In fact, the following statements show that it also contains information about the graphical structure of the CRN.Proposition 4.2 Let the complex C j ∈ B 1 have an explicit factorization of the form (11) For weakly reversible networks, U s = U, and hence , im(U s ) = im(U) .Therefore, the parameter θ t can be removed by merging its value with the term Uξ t , which yields a stoichiometric factorization for each C j ∈ B 1 .As a result, type-I nonstoichiometric BCs do not emerge in weakly reversible reaction networks.
In particular, it can be shown that blocked irreversible reactions indexed by strictly positive entries of A zl T θ t , t = 1, 2 are "bridge reactions" that connect distinct strong linkage classes of the network.In fact, removal of all such blocked reactions increases the number of linkage classes by at least two.Proposition 4.3 shows an interesting contrast to a well-established result in chemical reaction network theory, which states that (full) complex balancing can be obtained at a positive steady state, only if the network is weakly reversible [Proposition 16.5.7 in 1 ].
Interestingly, once the blocked reactions predicted by Proposition 4.1 are removed, a number of strong linkage classes will be detached from the rest of the network and form separate linkage classes.
The following corollary follows immediately from Proposition 4.1.

Corollary 4.5
Let G be a network with all its blocked reactions have been removed; then we have B 1 = B 2 and B \ B 1 = B \ B 2 ; that is, all nonstoichiometric BCs are of type II.
By contrast to type-I nonstoichiometric BCs, the type-II BCs do not necessarily require the network to contain blocked reactions at steady state; hence, they do not automatically follow from modified topological features of the network at steady state.However, a rather similar trend could be observed: For a type-II nonstoichiometric BC to exist, a number of reactions must be fixated at nonzero lower-or upper bounds.Proposition 4. 6 Let G be a network, all blocked reactions of which have been removed.Suppose B \ B 1 is nonempty.Then G contains at least three reactions fixated at a corresponding nonzero lower-or upper bound, for all v ∈ F (G).

Nonstoichiometric BCs and the integration of phantom species
Balanced complexes have been shown to play a key role in the reduction of large-scale metabolic networks [see 12 and references therein].They shall be viewed as intrinsic properties of the network that shed light on its steady state characteristics regardless of the underlying kinetics governing the conversion of species.The nonstoichiometric BCs have an additional interesting facet to themselves: They carry extra information about the network structure, which is not contained in the steady state Eq. ( 3) or in the linkage closedness condition U T Av = 0 , but emerges as a joint product of the stoichiometry, graphical structure, and flux constraints.
In what follows, we introduce transformations on a given reaction network G , which encodes these extra pieces of information into the steady state equation of a modified network G * , which has the same graphical but a different algebraic structure.This will enable us to exploit G * to make observations about steady-state characteristics of the original CRN G .Even though only nonstoichiometric BCs carry extra pieces of informa- tion to encode, we shall next begin by introducing these transformations in the more general case, i.e. for any arbitrary BC C b ∈ C.

Insertion of phantom species.
Let G be a CRN and C b ∈ B(G) be any balanced complex in G , i.e.

e b
T Av = 0, ∀v ∈ F (G) .Suppose we 'amend'G by inserting one molecule of some phantom species σ into complex C b .This yields the modified network G * = (S * , C * , R * ) , where S * = S ∪ {σ }, Accordingly, the stoichiometric map and incidence matrix for G * are given by We refer to this transformation as injection of vector e b into network G .As this procedure is only applied to some balanced complex C b and it involves insertion of some new species (not already in the network), the fol- lowing result is trivially obtained.Lemma 5.1 Let G be a network,C b ∈ B(G) and G * the modified network obtained by injection of e b into G .Then G and G * have identical steady state flux distributions, that is, F (G) = F (G * ) .In particular, any balanced complex of one is a balanced complex of the other, B = B * .Moreover, we have the inclusion While G and G * have identical balanced complexes, the stoichiometric nature of BCs may differ across the two networks.In particular, for the modified complex, we have C b ∈ B * 3 , regardless of its original categorization in G .For any chosen C b ∈ B , this transformation preserves the steady state flux space, and the information carried by any such BC is hereby encoded into the steady state equation.Clearly, the number of nonstoichiometric BCs strictly decreases in this process, in case of Characteristics of the modified network.The next question to address is how the stoichiometric subspace develops under the introduced transformation.On the face of it, this seems to contradict the statement of Lemma 5.1.However, the ambiguity is easily sorted out, because not every vector in ker(N) lies in F (G) , simply due to irreversibility patterns and flux bounds.The appar- ent difference is only due to the fact that the extra piece of information carried by nonstoichiometric BC C b is now encoded in the steady state equation The following result follows readily from the definition of deficiency: The deficiency is often characterized as a measure of how tightly the complex formation rate vector, i.e. the image of the steady state flux space under mapping A , is constrained at steady state 1 , that is, the algebraic dimension of the deficiency space.Let us next focus on cases where C b is a nonstoichiometric BC.Given that G and G * have the exact same flux distributions and the exact same incidence matrix, the contrast in their deficiency values is intriguing.
It is worth noting that the deficiency, as defined in (5) or (6), completely neglects of the fine details of the directionality of reactions, or the active bounds on flux through reactions.The reaction arrows exert their influence only to the extent that they serve to partition the complexes into linkage classes 1 .For the rather abstract notion of unconstrained networks, it is trivial to show that the deficiency value coincides with the dimension of the deficiency space, that is However, in more practical settings, e.g. the constrained framework presented in flux balance analysis of metabolic networks, the deficiency -as defined-should be understood only as an upper bound on the dimension of the deficiency space, that is, the following more general relation holds In the case of G and G * , the two networks form identical complex formation rate vectors at steady state, that is, have identical deficiency spaces.However, the nominal deficiency of G,δ, is larger than that of G * , i.e.δ * = δ − 1 .This sug- gests that, due to abovementioned factors, the network G is in fact more constrained at steady state than is reflected by δ = ker(Y) ∩ im(A) .Those factors cause the deficiency of G to be effectively less than δ , and no larger than δ * .

Effective deficiency
Let G be a network and be the set of nonstoichiometric BCs for G .Let G 1 be the network obtained by injection of vector e b 1 into G .Given the discussion in Section "Insertion of phantom spe- cies", δ 1 = δ − 1 .Moreover, the equality is a nonstoichiometric BC for G 1 .Without loss of generality, let it be labelled C b 2 .Let G 2 be the network obtained by injection of vector e b 2 into G 1 .It follows that δ 2 = δ 1 − 1 .
We then proceed further by searching the set C b j |B\B 2 | j=3 for any complex that may lie in B \ B 2 (G 2 ) .If it exists, we construct G 3 by injection of the corresponding vector, and so on.
As long as the most recent network in the sequence G → G 1 → G 2 → • • • has any nonstoichiometric BC, one can basically repeat this process of inserting distinct phantom species and obtain the next modified network in the sequence.After repeating this process consecutively d times, for some d , we will obtain B 2 (G d ) = B ; that is, there remains no nonstoichiometric BC left to carry on with the process.
Our first observation is that, by virtue of Lemma 5.1, all networks G, G 1 , • • • , G d have identical steady state flux distributions.Since they all have the same incidence matrix, they also share the same deficiency space, i.e. the steady state variations of the complex formation rate vector is also identical across all such networks.However, by virtue of Proposition 5.4, It follows that, due to constraints on the flux space imposed by irreversibility patterns and flux bounds, the deficiency of network G is effectively no larger than δ − d.
In order to be able to define the notion of effective deficiency unambiguously, the question of uniqueness needs to be addressed: Had one chosen a different order for the sequence of nonstoichiometric BCs to inject and construct G → G 1 → G 2 → • • • , would one have still obtained the same value for d ?The next statement addresses this question.

Theorem 6.1 The maximum length of the sequence
that is,d is independent of the choice and order of nonstoichiometric BCs used to construct it.Theorem 6.1 paves the way for the definition of effective deficiency for networks under irreversibility constraints and flux bounds.Definition Let G be a network and B \ B 2 the set of nonstoichiometric BCs for G .Let G 1 , • • • , G d be any arbi- trary sequence of modified networks constructed iteratively by means of injecting nonstoichiometric BCs, and suppose B 2 (G d ) = B .The effective deficiency of G , denoted δ eff , is defined as follows.
The following result facilitates the computa- tion of effective deficiency, without having to explicitly construct any sequence of modified networks.
is the set of nonstoichiometric BCs for G .Let us construct the matrix One can seamlessly replace the set B \ B 2 in Proposition 6.2 with the set of all BCs B (and replace E B\B 2 by E B , which contains all balancing vectors as columns).Nevertheless, the exact same equality will still hold: However, the same value of δ eff would be obtained at the expense of a higher computational effort.This is due to the fact that injection of stoichiometric BCs has neither an impact on the deficiency of the network nor on its flux distributions, simply because it incorporates no additional information into the stoichiometric structure.
We remark that, by the virtue of Eqs.(18) or (19), once the balanced complexes have been identified, obtaining the effective deficiency reduces to rank computation, which can be conducted by performing Gaussian elimination to obtain the row echelon form.Hence, the computation of effective deficiency has a time complexity of O n 3 , at most.In more technical terms, it can be shown that the effective deficiency has to do with the codimension of im Y T + im(U) , as stated in the following proposition.
where V /W denotes the quotient space of V by W.
Considering that the effective deficiency reflects the true dimension of the deficiency space in a network under functional and operational constraints, it is conceivable that it is in fact the effective deficiency, and not the nominal value of (structural) deficiency, that determines the steady state characteristics of a metabolic network, e.g.allowing unique equilibrium and excluding exotic dynamic behavior.We present this intuition as the following conjecture.Conjecture 6.4 Let G be a network of arbitrary deficiency, potentially larger than one.Suppose the effective defi- ciency of G satisfies either the conditions of Deficiency Zero Theorem or those of Deficiency One Theorem.Then, the results of the corresponding theorem, such as uniqueness of the equilibrium also apply to G.While we do not prove this statement here, we would like to emphasize that all the results presented here point at this property.In particular, we have shown that there is a one-to-one correspondence between the flux space of G and that of the modified networks obtained by insertion of the phantom species, which reflects itself in equivalent deficiency spaces.It seems only logical that if a modified network does not exhibit exotic steady state behavior, the same property should hold for G.

Some examples
In Section "Effective deficiency", we introduced the notion of effective deficiency, and established its relation to the existence of nonstoichiometric BCs in the network.In the following, we illustrate this phenomenon via a number of examples.Here, we present toy networks, the design of which does not necessarily reflect real-world chemical or metabolic pathways, but is actually tailored for the purpose of illustration.On that account, the unrealistic nature of these toy networks shall not be viewed critically, as the goal is to simply show a couple of nonstoichiometric BCs and their consequences within a small and representable setting.
As far as real-world applications are concerned, the presence of nonstoichiometric BCs across a wide range of metabolic models has already been confirmed in a recent study 7 , which implies the exact same consequences for the effective deficiency of those metabolic networks.
Toy network I: a type-I nonstoichiometric BC.Let us consider the toy network in Fig. 1.
The basic conversion diagram presented in Fig. 1 (including the reactions in brown) portrays a network operating in a canonical flux regime with n = 9 complexes,ℓ = 2 linkage classes and of rank s = 5 .As a result, the network is of deficiency δ = n − ℓ − s = 2 .It contains no stoichiometric BCs.However, it can be shown that the complex ( 2A ) has a nonstoichiometric factorization of the form Eq. ( 11), hence, it is a type-I nonstoichiometric BC.The detailed parameter values for the factorization are presented in the Supplementary.
As predicted by Proposition 4.1, the presence of a type-I nonstoichiometric BC is accompanied by having two irreversible reactions blocked at steady state, both of which are highlighted in brown.Moreover, (2A) is the only balanced complex, hence, d is bound from below and above to be exactly one.Therefore, while the network's nominal deficiency is two, it follows from the above results that the network is of effective deficiency Let us next consider the reduced network obtained by removing all blocked reactions, but including all reactions shown in black.The reduced network will contain no nonstoichiometric BCs, since the balanced complex (2A) has turned into a stoichiometric BC in this network.For the reversible reactions, the larger arrow size depict the direction of the flux associated with a positive sign.In line with the prediction of Propositions 4.6 and 4.7, the network contains a number of reactions fixated at corresponding lower-and upper bounds, shown in blue and red, respectively.
The toy example in Fig. 2 portrays a network operating in a non-canonical and bounded flux regime.The three complexes (A + B) , (A + C) and (D + E) are strictly stoichiometric BCs, due to the fact that species B,C , and E do not appear elsewhere in the toy network.The network has n = 6 complexes,ℓ = 1 linkage class and rank s = 4 .Hence, the nominal deficiency of this network is δ The other three complexes, that is,(2A), (A + D) and (2D) are type-II nonstoichiometric BCs.One can basi- cally associate each with a nonstoichiometric factorization of the form (10) (details are available in the Supplementary).The existence of a nonstoichiometric BC implies that the effective deficiency δ eff of this network must be strictly smaller than its nominal deficiency.Coupled with the knowledge that the deficiency only takes nonnegative values, it follows that δ eff = 0 .In other words, this is effectively a deficiency zero network, and thereby it is complex-balanced.
In line with the statements of Propositions 4.6 and 4.7, the network contains a number of reactions fixated at their lower-and upper bounds, shown in blue and red, respectively.It is worth noting that the formation of nonstoichiometric BCs in this network rely on the irreversibility of R 1 and R 2 , the non-canonical flux regime imposed by strictly positive lower bounds on fluxes of of R 1 and R 2 , as well as the unbounded flux regime imposed by finite upper bounds on fluxes of of R 3 and R 4 .These imposed constraints force the network to perform like a deficiency-zero network and hence be complex-balanced, despite the fact that its underlying structure had the capacity to perform like a network of higher deficiency.
Effective deficiency of large-scale metabolic networks.We use twelve genome-scale metabolic networks [13][14][15][16][17][18][19][20][21][22][23][24] from all kingdoms of life obtained from Küken et al. 12 to investigate deficiency in networks with two sets of constraints (i) imposing reaction irreversibility constraints as specified in the original model reconstruction, and (ii) imposed optimality of specific growth rate in addition to reaction reversibility constraints.Note that blocked reactions are removed from the networks before BC detection and, therefore, the networks do not contain type-I nonstoichiometric BCs.In other words, all type-I nonstoichiometric BCs are transformed into stoichiometric BCs as a result of these removals (see Proposition 4.4).In this setting, we compare the structural deficiency and effective deficiencies obtained under the two different sets of constraints.The structural deficiency ranges from 57 for T. maritima to 1097 for S. cerevisiae (Fig. 3).We find that the effective deficiency obtained from scenario (i) is the same with the structural deficiency for the networks of A. thaliana, M. musculus, N. pharaonis, and P. putida.For the remaining eight networks, we find a reduction of the effective deficiency in comparison to the structural deficiency, ranging from 1.8% in E. coli and 35.3% in M. barkeri (Fig. 3).Considering the additional constraint on the specific growth rate, fixed at its optimum, in scenario (ii) we observe the effective deficiency to be smaller than the structural deficiency in all networks, with the smallest decrease in networks of A. thaliana (0.8%), E. coli (1.8%), M. musculus (3.6%) and A. niger (4.2%) and the largest decrease in the networks of M. barkeri (35.3%),M. acetivorans (11%), C. reinhardtii (9.9%) and T. maritima (9.6%).

Conclusion
Here we introduced the notion of effective deficiency that takes into account not only the structure of the network but also its operational constraints in the calculation of deficiency.The effective deficiency relies on the presence of nonstoichiometric balanced complexes, which we have shown to be present in large-scale metabolic networks across kingdoms of life.In addition, our results point at a subtle relation between effective deficiency and robustness of some reaction fluxes.Future work will aim at employing these findings in characterizing classes of networks that exhibit particular dynamical properties that can be ensured or precluded based on the notion of effective deficiency.
C and C′ are strongly linked, denoted C ≈ C′ , if C ⇒ C′ and C′ ⇒ C .The equivalence relation ≈ partitions C into a family of ℓ s equivalence classes {� l } ℓ s l=1 , called the strong linkage classes of the network.

Proposition 5 . 2 Corollary 5 . 3
Let G be a network,C b ∈ B(G) and G * the modified network obtained by injection of e b into G .Let us denote the stoichiometry matrices of G and G * by N and N * , respectively.Then, if and only if C b ∈ B\B 2 .Moreover, rank(N * ) = rank(N), if and only if C b ∈ B 2 .rank N * = rank(N) + 1, It is only the injection of nonstoichiometric BCs that increases the rank of the stoichiometric subspace.This yields the following corollary.Let G, G * be as defined in Proposition 5.2, and C b ∈ B\B 2 (G) .Then, dimker(N * ) = dimker(N) − 1.

Proposition 5 . 4
Let G be a network,C b ∈ B(G) and G * the modified network obtained by injection of e b into G .Let us denote the deficiency of G and G * by δ and δ * , respectively.Then.if and only if

Proposition 6 . 3
Suppose B = C b j |B| j=1 be the set of BCs for G , and E B defined as E B = e b 1 e b 2 • • • e b |B| .Then, the drop in deficiency,d , is the codimension of It has n′ = 9 complexes,ℓ′ = 4 linkage classes and rank s′ = 4 .It follows that the reduced network is of deficiency δ′ = n′ − ℓ′ − s′ = 1 .This is in accord with the value calculated for the effective deficiency of the original network.Toy network II: a type-II nonstoichiometric BC.Next, let us consider the conversion diagram shown in Fig. 2.

Figure 1 .
Figure 1.A type-I nonstoichiometric BC in a toy network.The basic conversion diagram depicts a network with m = 6 species, n = 9 complexes, and r = 8 reactions.The network does not contain any stoichiometric BC.However, it contains exactly one type-I nonstoichiometric BC (complex 2A ).In line with the prediction of Proposition 4.1, the network also contains two blocked irreversible reactions (highlighted in brown).

Figure 2 .
Figure 2. Type-II nonstoichiometric BCs in a toy network.The conversion diagram portrays a network with m = 5 species, n = 6 complexes, and r = 7 reactions.The network operates in a non-canonical and bounded flux regime, where the fluxes of irreversible reactions R 1 and R 2 have strictly positive lower bounds ( v l,1 = v l,2 = 100 ), while the fluxes of reactions R 3 and R 4 have finite upper bounds ( v u,3 = v u,4 = 200 ).The network contains three type-II nonstoichiometric BCs, highlighted in yellow.The rest of complexes are stoichiometric BCs.

Figure 3 .
Figure 3. Structural and effective deficiency of real-world metabolic networks.Balanced complexes are identified in networks of twelve species from all kingdoms of life under two different scenarios: (i) imposing reaction irreversibility constraints, and (ii) considering in addition specific growth rate optimality.The effective deficiency is smaller than the structural deficiency when the additional constraint on optimality of specific growth rate is imposed for all networks.
One can actually come up with a slightly more informative statement.For a network G , suppose B \ B 1 is nonempty.Then G contains a reaction R ∈ R irr \ R zl fixated at a positive lower bound.Moreover, there exists another reaction R′ ∈ R fixated at a positive upper bound, or a reaction R′ ∈ R rev fixated at a negative lower bound.Note that none of these results imposes any specific restrictions on the graphical structure of the network, unlike the type-I case.For type-II nonstoichiometric BCs, it is not the topological feature but the imposed flux constraints that set the stage for the emergence of additional balanced complexes.The next result follows immediately from Proposition 4.7.Given a network G , for the set B \ B 1 to be nonempty,G must operate in a bounded and non-canonical flux regime.Conversely, suppose G is a network operating under an unbounded and/or canonical flux regime, and it contains no blocked reactions.Then, B = B 1 = B 2 .
) Scientific Reports | (2023) 13:14589 | https://doi.org/10.1038/s41598-023-41767-1www.nature.com/scientificreports/ ,If the conjecture holds, it can have interesting implications even for large metabolic networks.Consider a metabolic network whose linkage classes are (i) independent subnetworks, and (ii) each of effective deficiency zero or one.Such a network then, by this conjecture, cannot admit more than a single steady state within any positive stoichiometric compatibility class, regardless of the values of the rate constants 1,2 .