{'fraction_non_alphanumeric': 0.03800098419815189, 'fraction_numerical': 0.006479304499972661, 'mean_word_length': 4.541433116194516, 'pattern_counts': {'":': 0, '<': 1, '<?xml version=': 0, '>': 3, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 1, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Boolean constraint satisfaction problems for reaction networks

We define and study a class of (random) Boolean constraint satisfaction problems representing minimal feasibility constraints for networks of chemical reactions. The constraints we consider encode, respectively, for hard mass-balance conditions (where the consumption and production fluxes of each chemical species are matched) and for soft mass-balance conditions (where a net production of compounds is in principle allowed). We solve these constraint satisfaction problems under the Bethe approximation and derive the corresponding Belief Propagation equations, that involve 8 different messages. The statistical properties of ensembles of random problems are studied via the population dynamics methods. By varying a chemical potential attached to the activity of reactions, we find first order transitions and strong hysteresis, suggesting a non-trivial structure in the space of feasible solutions.


I. INTRODUCTION
Biological networks map out the complex set of interactions that may occur among different units (genes, proteins, signalling molecules, enzymes, etc.) in cells [1,2]. Their structure is thought to reflect, at least in part, the specific physiologic function(s) they are meant to control and, with the topology-to-function mapping mostly still unclear, many important questions can be formulated about the optimality, robustness and evolvability of these systems [3][4][5][6][7][8][9][10]. On the other hand, functional control in cells is achieved through the physical dynamics that takes place on the networks, which is usually much more complicated than the network structure by itself would suggest. To make an example, consider transcriptional regulatory networks. While their structure only encodes for the possible protein-DNA interactions by which the transcription of RNA can be turned on or off, regulation results from the reciprocal adjustment of transcriptional activity and protein levels. This process however involves a variety of regulated steps, like DNA-binding and unbinding events by multiple proteins (possibly preceded by the formation of protein complexes), RNA polymerization (by specifically recruited molecular machinery) and transport, post-transcriptional modification events and, finally, translation. Each node in this network therefore lumps together a number of molecular species and elementary processes, all of which are subject to noise. In such a complex interacting environment, the overall patterns of activity may be hard to uncover even if one is only interested in steady states.
It is tempting, then, to implement a coarse-grained approach and explore the possibility of characterizing the operation of biological networks through simpler, perhaps Boolean, dynamical rules or, at an even more basic level, through elementary feasibility constraints [11][12][13][14][15][16]. This type of scheme is especially suited (a) to identify robust attractors of the dynamics and/or groups of nodes that are likely to behave in a highly correlated way (viz. the emergence of network motifs discussed in [13]), and (b) to evaluate 'degrees of activity' for the different nodes, by which one may, for example, guide more refined techniques that simulate the full dynamics of the system towards physiologically relevant states. Besides, from a purely theoretical viewpoint, very often the problems thus defined present phase structures and algorithmic challenges that suffice by themselves to attract a considerable statistical mechanical effort (see e.g. [17,18]).
In this paper we define and study a Boolean constraint-satisfaction problem (CSP) designed to represent minimal operational and stability requirements for the non-equilibrium steady states (NESS) of biochemical reaction networks, like the metabolic networks that relate enzymes to the substrates and products of the reactions they catalyze in any given cell type [19,20]. In essence, we shall enforce feasibility constraints that link enzyme activity to substrate and product availability, and vice-versa, similarly to the approach defined in [21][22][23]. From a physical viewpoint, the model describes, in different limits, different types of NESS, and therefore different physiological scenarios. The corresponding CSPs, on the other hand, turn out to be of a novel type, requiring ad hoc message-passing methods to be analyzed in detail.
This article focuses on the properties of CSPs defined on ensembles of artificial (random) reaction networks -introduced in Section II -which will be studied by the statistical physics tools sketched, together with the corresponding results, in Section III (and fully exposed in the Appendix). A subsequent work will be concerned with the analysis of solutions for single random networks and real metabolic networks.

A. Random Reaction Networks
We define a random reaction network (RRN) to be a bipartite random graph with two types of nodes, representing respectively chemical species (or metabolites) and enzymes (or reactions). We shall denote by N and M , respectively, the number of reactions and that of metabolites. Both N and M will be taken to be large, i.e. N, M 1. For sakes of simplicity, we shall assume here that each reaction has a well defined operational direction, so that the bipartite graph is directed. Its topology will be encoded in an adjacency matrix ξ, with entries ξ m i = 1 if reaction i produces metabolite m, ξ m i = −1 if reaction i consumes metabolite m, and ξ m i = 0 otherwise. We furthermore define ∂m in (resp. ∂m out ) as the set of reactions producing (resp. consuming) m; likewise, for each reaction i, ∂i in (resp. ∂i out ) will denote the set of its substrates (resp. products).
The topology of the RRN is specified by the probability distributions of the degrees of the two node types. For metabolite nodes we shall assume that the in-and out-degrees in ≡ |∂m in | and out ≡ |∂m out | are independent random variables, both distributed according to a Poissonian with parameter λ, i.e.
Metabolites having ( in , out ) = (0, 0) are disconnected from the network and will be ignored in what follows. We shall generically assume that λ ≥ λ p = 1 (the percolation threshold), ensuring the existence of a 'giant' connected subgraph. Metabolites with ( in , out ) = (0, ≥ 1) represent the substrates that the reaction network derives from the environment (the 'nutrients'), whereas metabolites with ( in , out ) = ( ≥ 1, 0) will be considered to be the final products or sinks (e.g. excreted compounds or molecules that are employed in processes other than chemical reactions) of the network. The fraction of such 'leaves' (nutrients or sinks) is given by e −λ (1 − e −λ ) e −λ for large enough λ. Likewise, for reaction nodes, the quantities |∂i in | ≡ d in and |∂i out | ≡ d out will be assumed to be independent random variables, both distributed according to where 0 ≤ q ≤ 1 is a parameter. In other words, reactions can be of four different types according to their in-and out-degrees ((d in , d out ) ∈ {(1, 1), (1, 2), (2, 1), (2, 2)}) and q weights the relative number of bi-component reactions (as inputs, outputs or both). The only structural control parameters that we shall use in the following are the mean degrees of metabolites (λ) and of reactions (q). A sketch of the network is given in Fig. 1. j k n m p o i Figure 1: Sketch of a random reaction network of the type discussed in the text. m, n, o and p denote metabolites (squares), i, j and k are instead reactions (circles). Red (resp. blue) links carry substrate-like (resp. product-like) couplings with ξ m i = −1 (resp. ξ p i = 1).

B. Constraints: Hard versus Soft Mass Balance
In order to define a CSP embodying realistic operational constraints, we focus on the characterization of the NESS induced by non-zero in-and out-fluxes of nutrients and sinks, respectively [24,25], following two different (but related) schemes. In Flux-Balance-Analysis (FBA) it is assumed that fluxes in NESS ensure mass balance at each metabolite node in the network [26]. (We do not consider here the optimization schemes that are typically coupled to such constraints in biological implementations of FBA [19,27].) If we denote by J i the flux of reaction i (with J i ≥ 0 for an irreversible reaction), this amounts to solving the system where σ m i is the stoichiometric coefficient of metabolite m in reaction i (such that sgn(σ m i ) = ξ m i ). One easily understands that the above conditions are equivalent to Kirchhoff's node laws for the flow of matter through metabolite nodes and describe NESS with constant (time-independent) levels for each metabolite. A soft version of this model [28,29] assumes instead that intracellular concentrations may be allowed to increase linearly over time at constant rate, e.g. because some metabolites have to be available for processes outside of metabolism strictly defined. This simply leads to replacing (3) More formally, the above conditions can be seen to derive from Von Neumann's optimal growth scenario [30] and provide a useful means of characterizing a reaction network's production capabilities [31][32][33]. The general problem posed by (3) and (4) consists, given the matrix σ = {σ m i }, in retrieving the flux vectors J = {J i } satisfying the M linear conditions. An interesting feature that is observed in the solutions of the above models is that a sizeable fraction of reactions carries a null flux in each solution [32,34]. This suggests that, to a first approximation, if one is interested in capturing certain aspects of NESS within a coarse-grained description it might suffice to just distinguish, for each reaction, the inactive state from the active one. We shall then introduce, for each reaction, a variable ν i ∈ {0, 1} (inactive/active). Similarly, we shall link to every metabolite a variable µ m ∈ {0, 1} that characterizes whether that particular chemical species is available (µ m = 1) or not (µ m = 0) to enzymes that process it. Our next task is to devise Boolean CSPs that embed the basic features underlied by (3) and (4), respectively.
Starting from (3), it is simple to understand that a minimal necessary requirement that is encoded in the massbalance conditions is that, for each metabolite which is produced by an active reaction, there must be at least one active reaction consuming it, and vice-versa. This means that, for each compound m, all assignments of ν i 's are acceptable except those for which all active reactions either produce or consume it. We can therefore define the number of active reactions producing and consuming chemical species m as and, in turn, introduce an indicator function Γ m ≡ Γ m (µ m , {ν i }) for every m as Given a configuration {ν i }, metabolite m will be said to be SAT when Γ m = 1, i.e. when either no reaction in which it is involved is active (x m = 0 and y m = 0) and the metabolite is unavailable (µ m = 0), or when the metabolite is available (µ m = 1) and at least one reaction produces it (x m > 0) and at least one reaction consumes it (y m > 0). Similarly, we define a reaction to be SAT when the indicator function ∆ i ≡ ∆ i (ν i , {µ m }), given by with ∂i = ∂i in ∪ ∂i out , equals 1. That is, i can be active only if all its neighbouring metabolites (including both substrates and products) are available. We note that (7) can actually be re-cast as according to which i can be active only if all of its inputs are available: it is indeed clear that if a reaction is active but one (say) of its products is unavailable, then the constraint imposed on the metabolite will either be violated or force that metabolite to become available. Notice that ∆ i = 1 does not imply that i is active when all of its substrates are available.
The CSP corresponding to (3) can then be formulated as follows: find a non-trivial assignment of ν i 's (ν i 's not all zero) such that all reactions and all metabolites are SAT, i.e. Γ m = 1 ∀m and ∆ i = 1 ∀i, with Γ m and ∆ i given by (6) and (8), respectively. We shall call this CSP Hard Mass Balance, or Hard-MB for brevity.
In order to get a Boolean representation of (4), we note that the main difference between this case and that of (3) is that, because of the soft constraint, it is no longer necessary that production fluxes are balanced by consumption fluxes. Therefore, while constraint (8) remains valid, (6) has to be replaced by In other terms, metabolite m can be available as soon as at least one reaction producing it is active (x m > 0). It is convenient to re-write Γ m for this case as so that the constraint at each metabolite node only includes incoming degrees of freedom, making the directionality inherent in the corresponding CSP explicit. It is straightforward to see that (9) or (10), together with (8), which retains validity, return the same configurations.
The CSP corresponding to (4) is then the following: find a non-trivial assignment of ν i 's such that all reactions and all metabolites are SAT, i.e. Γ m = 1 ∀m and ∆ i = 1 ∀i, with Γ m and ∆ i given by (10) and (8), respectively. We shall call this CSP Soft Mass Balance, or Soft-MB for brevity.
Note that the constraints behind the two problems can be written compactly as where α = 1 for Hard-MB and α = 0 for Soft-MB. We shall be interested in solutions obtained upon fixing the probability that a nutrient is available, which we denote below as ρ in . For the moment, no specific assumption will be made on sinks.

III. RESULTS
In a nutshell, the above setup aims at retrieving Boolean patterns of activity of reactions (or of metabolite availabilities) induced, on network architectures defined by q and λ, by the fact that a certain set of metabolites (nutrients) is available from the outset. Ideally, one would like to devise a method to sample configurations (ν = {ν i }, µ = {µ m }) with a controlled probability given by which forbids states that don't satisfy all constraints. The 'chemical potential' θ appearing above can be tuned externally in order to concentrate the measure around configurations with a different average fraction N −1 i ν i of active reactions, where angular brackets represent the average with respect to the measure (13).

A. Soft Mass-Balance
In order to find the configurations of reaction and metabolite variables that solve the above CSPs one may resort to statistical mechanics techniques. In particular, we have employed a cavity theory to devise a belief propagation/population dynamics algorithm to sample the probability distribution (13). Details about the theory and the algorithms are reported in the Appendix. We shall concentrate here on the scenario that emerges for different q and λ upon varying two parameters, namely the chemical potential θ and the probability ρ in that nutrients are available. In specific, we have computed the average reaction activity and the average metabolite availability following two protocols: first, by gradually reducing θ starting from a large, positive value, and, second, by doing the reverse. Averages obtained in these ways will be denoted, respectively, by · · · + and · · · − . These averages (that we call magnetizations, using a statistical physics jargon) need not coincide, in which case the two quantities will display hysteresis when plotted against the chemical potential. Generally, the presence of hysteresis is a main characteristic of a discontinuous (first order) phase transition, while for continuous (second order) ones no hysteresis is observed, as also happens in cases where no phase transition takes place.
The average fractions of available compounds (metabolites) and active reactions obtained upon varying θ at fixed ρ in = 0.5 for Soft-MB is displayed in Fig. 2. One sees that, expectedly, larger values of θ lead, on average, to larger fractions of available metabolites and of active reactions. For large enough values of λ and q, however, as the · · · + and . . . − averages become steeper functions of θ, the curves obtained by increasing and decreasing θ no longer coincide. Notice that, while for lower λ and q solutions can be found over a broad range of values of the magnetizations, when λ and q increase the average metabolite availability seems to concentrate in small ranges close to the extremes 0 and 1, distinguishing solutions with few available metabolites from solutions with a large fraction of available compounds. This type of picture is however not observed for reactions (we shall return to this point later on). A simple way to quantify the onset of hysteresis is by measuring the quantity (we focus for simplicity on metabolites) which vanishes when µ + = µ − and generically differs from 0 in presence of hysteresis. A map of the values of ∆µ in the parameter space (λ, q) is presented in Figure 3 for the limiting choices ρ in = 1 and ρ in = 0. While hysteretic behaviour can be found practically all throughout the (λ, q) plane, it becomes stronger at high enough λ and q, where an abrupt jump in the magnetizations takes place. The presence of such a large hysteresis, and the coexistence of low and high magnetization solutions, signal a non trivial structure in the space of solutions to the CSP. Such a non-trivial structure appears also in many other well-known CSP, as the random k-XORSAT [35] and random k-SAT [36,37], and is the origin of the onset of long range correlations, that have important consequences on the behavior of searching algorithms [38]. Away from the hysteretic portion, µ + and µ − vary smoothly with θ, allowing one to sample easily solutions with any magnetization not in the jump.
On the other hand, the overall structure of the solutions (in terms of ∆µ) appears to vary weakly with ρ in . This strongly suggests that main observed effects (e.g. the jump and the hysteresis) are essentially due to topology of the network, rather than to the boundary conditions. It is interesting to observe that the hysteretic region shrinks as ρ in increases, suggesting that, within the constraints imposed by Soft-MB, a larger repertoire of available nutrients stabilizes the output by allowing to achieve higher values of the magnetization for smaller values of θ.

B. Hard Mass-Balance
The ∆µ-map for the Hard-MB case is displayed in Figure 4. In contrast with the Soft-MB case, Hard-MB solutions display strong hysteresis for all choices of λ, q and ρ in . Furthermore, comparing the results at ρ in = 0 and ρ in = 1, it is clearly seen that, again, changing ρ in (i.e. increasing the number of available nutrients) has little influence on the overall structure of the phase space. Rather, its main effect is that of reducing the magnitude of hysteresis cycles. It is interesting to note that the maximum value of ∆µ in Hard-MB is more than double than the one in Soft-MB.
The presence of strong hysteresis markedly distinguishes the solution space of the two CSP problems. A comparison between the behaviour of the magnetization obtained in the Soft-and Hard-MB cases for selected parameter values is displayed in Figs 5 and 6. In first place, one sees that the limiting value of the average magnetization for θ → ±∞ in the Hard-MB problem is identical to that of the Soft-MB problem, suggesting that in specific cases the Hard-MB CSP may actually acquire a strong directional nature (like the Soft-MB case), despite the fact that in Hard-MB substrates and products are highly correlated between each other. Secondly, the increasing-θ protocol appears to be unable to identify active solutions in the Hard-MB case, suggesting that the Hard-MB constraints bias solutions towards activating a large fraction of metabolite nodes. In the Hard-MB case, it seems that it is possible to start from the all-on configuration and gradually switch off the network, but it is very difficult to switch on part of the network starting from the all-off configuration: for this reason the all-off solution is very stable in the Hard-MB case.    A more quantitative view of this is given in Fig. 7, where we display the distribution of values of the magnetization for metabolites and reactions obtained for a value of θ at the transition, where the difference in behaviour between Soft-and Hard-MB is more striking.
From the distribution of metabolite availabilities one clearly sees that, generically, fluctuations are larger in Soft-MB than in Hard-MB, implying that, while Soft-MB sustains non-trivial solutions over a wide range of values of the magnetizations, Hard-MB only admits solutions with a large and tightly constrained value of the average metabolite availability. Interestingly, the overall structure of the distributions changes when one considers reactions, for which both Soft-and Hard-MB can lead a large variability (much larger, in turn, than what occurs for metabolites). This is consistent with our constraints, which do not impose to activate a reaction even when all of its neighbouring metabolites are available. Note that both for reactions and metabolites Soft-MB allows for solutions with very low magnetization that are generically absent in Hard-MB.
Finally, we notice that not all of the solutions to Hard-MB would be able to carry non-vanishing fluxes in the linear problem defined by (3), which is only possible if the number of available metabolites does not exceed that of active reactions. To see this, one can compare the quantities M µ and N ν , see Figure 8 (left panel), which are respectively the number of equations and the number of unknowns in the FBA problem.
It is clear that only for sufficiently large values of λ will Boolean configurations correspond to realizable flux states in FBA. This confirms the intuition that redundant network structures (larger λ's) confer flexibility (i.e. the possibility of operating in different states) to a reaction network. What looks counterintuitive in Figure 8 (left panel) is that small q values are also to be preferred. An explanation to this fact can be obtained by plotting µ versus ν (right panel in Figure 8), and noticing that data with different q values fall on the same curve [42] Since the data in the left panel of Figure 8 are obtained by multiplying the x values in the right panel by N/M = λ/(1 + q), large q data are more keen to cross the line at the boundary of the feasible solutions region.
In the right panel of Figure 8 we also notice that µ spans a rather limited range (roughly 0.8 < µ ≤ 1) which is mostly independent on the topology (i.e., on λ), while the range of valid ν values becomes very broad for redundant networks (i.e., for large values of λ). In other words, solutions to the Boolean constrained problem on a RRN do exist only if a very large fraction of metabolites are present, while the fraction of active reactions can be made small only if the topology is redundant enough.

IV. OUTLOOK
In this work we have defined and studied a class of random Boolean CSPs representing minimal feasibility and operational constraints for reaction networks. We focused specifically on two sets of conditions: (a) mass balance for metabolites, corresponding to NESS with time-independent concentrations, and (b) mass unbalance allowing for a global net production of chemical species, corresponding to NESS with distinct metabolic output profiles linearly dependent on time. In both cases, we have been interested in computing statistical properties of the solutions induced by exchanges of the network with the environment, modelled by allowing for a fraction of chemical species (with precise topological characteristics) to be available from the outset, either as nutrients (basic inputs) or as sinks (basic outputs) of metabolism. We have displayed results obtained via a cavity-derived algorithm, whose theory is fully described in the Appendix.
A non-trivial dependence of the 'magnetizations' (i.e. of the average reaction activity and the average metabolite availability) on the chemical potential has been uncovered, characterized by hysteresis and, hence, by a first-order transition (in θ). While this phenomenology is present in both models, it is more marked in the Hard-MB case. Remarkably, the overall structure of the solution space (in terms of hysteresis) appears to be only weakly dependent on the value of ρ in , characterizing the size of the repertoire of available nutrients, suggesting a considerable robustness with respect to environmental changes, at least in the case of RRNs.
We argue here that many properties observed in this paper are typical of random networks, which have very peculiar topological properties: above the percolation threshold a random network is made of a unique connected component, having no modularity. In other words, the fact that neighbours in a RRN are chosen at random does not allow the network to have correlated structures at short scales (the so-called modules), and in this sense a random network is very different from real ones. On a modular network we expect the first-order-like transition to occur in each module separately, thus producing a more complicates pattern when the chemical potential θ is varied. Nevertheless in this work we are interested in the influence of the constraints on the exploration of the fixed point.
The method developed here will be employed in the analysis of the feasible patterns of activity of real metabolic networks. A (minor) adjustment required in order to port the tools described here to a realistic setting concerns the fact that, in physiological conditions, reactions may be reversible. It is possible to account for reversibility in the RRN model described here by simply extending the single reaction state space from {0, 1} to {0, ±1}. The theory carried out in the Appendix is indeed very easily generalized to the latter case. While such a generalization does not appear to bring out significant qualitative differences for RRNs, it is crucial to deal with genome-scale reconstructions of cellular metabolic networks. On the other hand, real instances of biological networks require algorithms that are capable of extracting individual solutions besides statistical properties. The methods developed here will serve as the necessary basis for accomplishing such a task. [42] The reasons why curves in the right panel in Figure 8 depend on λ but not on q rely on a particular property of a random graph, whose detailed derivation is outside the scope of the present paper. We leave this fact as an observation in the numerical data.
Appendix: Cavity approach for Soft-and Hard-MB problems

Cavity equations
In general a CSP, as Soft-MB or Hard-MB, can be solved efficiently on random networks by the belief propagation algorithm [39] or equivalently by the replica symmetric cavity method [40]. In this method, the marginal of a variable is computed by creating a "cavity" inside the system, removing a subpart of the network. Thus it is possible to obtain a "cavity marginal" and then reintroduce the variables removed. Finally the complete marginal of the variables follows directly from the cavity marginals.
In this kind of approach the system is represented by a graph made of "variable" and "function" nodes. Variable nodes are both metabolites and reactions, while a function node exists for each constraint. In the following we will use letters a, b, .. for the metabolite constraints and e, f, .. for the reaction constraints. Furthermore we introduce the condensed notations: if a is the constraint of metabolite m, then ∂a R = ∂a\m is the set of reactions involving metabolite m; if e is the constraint of reaction i, then ∂e M = ∂e\i is the set of metabolites involved in reaction i. Moreover by dividing in two groups reactions producing and consuming a given metabolite, we call ∂a R i the set of reactions in the same group as i, excluding i, and ∂a R ¬i the opposite group. Using this notation we can rewrite the constraints (8) and (10) for the Soft-MB case in a simpler form that will be useful in the computation of the equations of the system: where we have substituted δ xm,0 = j∈∂a R in (1 − ν j ), with a R in being the set of reactions producing the metabolite whose constraint is a.
In the same way, for the Hard-MB case, the constraints in the new notation can be written as: The representation of the graph that we obtain with the variables and the functions nodes is given in Figure 9 for Soft-MB and Hard-MB. As we can see the main difference between the two approaches is that the constraint a of the metabolite m changes its form.
From an algorithmic point of view, it is possible to search for the solutions to the cavity equations by an iterative procedure, in which "messages" are exchanged between variable and function nodes. For the two CSP defined in the present work, eight type of messages are required: ψ a→m µm , ψ m→a µm , η e→m µm , η m→e µm , ψ a→i νi , ψ i→a νi , η e→i νi , η i→e νi . Each message represents the belief that a variable (function) has about its neighbouring function (variable) state. The messages can be divided in two classes: from function nodes to variable nodes and from variable nodes to function nodes. The first class of messages, e.g. ψ a→m µm , is the probability that metabolite m is in state µ m when there is only the function metabolite a. While the second class, e.g. ψ m→a µm , is the probability that metabolite m is in state µ m when the edge (am) is not present. a a Figure 9: Summary of the cavity method messages for Soft-MB (left) and Hard-MB (right) constraints.
Introducing a parameter α to interpolate between Soft-MB (α = 0) and Hard-MB constraints (α = 1), the equations to be satisfied by the messages in the two CSPs under study can be written as follows: The equations we have just written hold in Soft-MB or Hard-MB case, with the difference that in Soft-MB case reaction nodes are connected only to output metabolite functions, while in Hard-MB all metabolite functions are connected to reaction nodes. Another caution we have to take is that the reaction function node is connected only to the input metabolites [see equations (15) and (16)] regardless of the constraints used. Writing explicitly the constraints we can compute the cavity equations obtaining, for the metabolite constraints, and, for the reaction constraints, Using these equations, we can iterate until convergence the algorithms presented in Appendix 3 and 4, finding solutions that satisfies the constraints and obtaining the cavity marginals. We can then compute the real marginals of the variable nodes as: where: (18) The main assumption behind the cavity method is that the messages coming from two neighbouring nodes are independent. This clearly depends on the length of loops in the network: if the length of typical loops grows with the system size (as in RRN), then the above assumption can be valid, at least in the thermodynamic limit. As we can see clearly from Figure 9 short loops are not present in the Soft-MB, but they arise in the Hard-MB problem. If the assumption breaks down, then message passing algorithms may fail to converge, although it has been observed that also in networks with many short loops (as e.g. regular lattices) it is possible to find solutions to the cavity equations by iterative message passing algorithms [41]. In the two CSPs under study we have observed a better convergence of the algorithm for Soft-MB constraints (α = 0) than for Hard-MB constraints (α = 1), and to ensure convergence in any case we have reached the Hard-MB limit by varying smoothly the α parameter.

Nutrients and Outputs
A B a Figure 10: Representation of the external metabolites in our network. A is the product while B is the nutrient.
How to deal with the nutrients (metabolites with in-degree 0 and out-degree larger than 0) and the outputs (metabolites with in-degree larger than 0 and out degree 0) is probably the trickiest part of the network analysis. Indeed looking at the cavity equations derived in 1, we immediately see that nutrients and outputs are automatically switched off because in these cases ψ a→m µm = δ µm,0 , while in real systems these variables are usually active, as they represent the interaction with the environment. To overcome this limitation we will consider in the following that nutrients are external variables fixed by the environment and thus have a probability of being present ρ in . Furthermore these variables send a message to the neighbouring reaction-constraint of the type:  On the other hand the products are internal variables with no reaction constraint node associated, a probability of being present p(µ m ) = ψ a→m µm (taken from 17) and send a message: to the neighbouring metabolite-constraint node. Furthermore the metabolite-constraint has to be a Soft-MB constraint otherwise the outputs would be always off in the Hard-MB case. A schematic view of the form of the network for the external metabolites is presented in Figure 10.
In principle for the outputs it is possible to define another parameter ρ out as: to force the network to switch on a fraction of outputs. Nevertheless it is then required to check at convergence of the algorithm that this value is consistent with the value of p(µ m ). A simple way to check this is by measuring µ out = i∈outputs µ i /N out and checking if this value is consistent with the value of ρ out given as a parameter. In Figure 11 this check is done for a particular case. In this Figure we see that there is only one value of ρ out consistent with µ out (this result holds similarly for Soft-MB and Hard-MB, and for any q and λ). We then verified that this solution is exactly the same as the one obtained by using equation (20), hence showing that the parameter ρ out is not necessary to explore all the possible solutions of the outputs. As a consequence of this setting on the inputs and outputs, the presence of the nutrients is determined by the parameter ρ in while the presence of the outputs is determined at convergence depending on the state of the network.

Belief Propagation Algorithm
In order to solve the equations we can run an iterative algorithm until convergence of the value of the messages. Thus the outcome of this algorithm is a set of probabilities of the states in which the variables are. This kind of algorithm is called Belief Propagation (BP).
The BP algorithm that we used to sample the solutions of the system is the same for Soft-MB and for Hard-MB. In this algorithm, first we generate a graph with a given q and λ. Then we initialize the messages with a random value and we iterate the equations until convergence. The simplest way to sample the solutions is by fixing one of the two free variables remained: θ or ρ in . By changing ρ in we can see how the configuration of the solutions changes when the nutrients have a probability ρ in of functioning. Whereas by changing θ we can observe what happens if we constrain the system to switch on (or off) the reactions. Each behaviour is interesting to understand how the system is organized. In each case we will compute the mean over the metabolites, µ , and the reactions, ν , where x is the average over the measure P (µ, ν), (13).

Population dynamics
BP is an algorithm for inferring the marginal probabilities on a specific graph. However, when one is interested in the behavior of typical samples of the RRNs with given parameters q and λ, then the equations presented in Appendix 1 can be solved using population dynamics [40]. The idea behind this approach is that, instead of computing the messages on a given graph, one considers the probabilities, P (ψ) and Q(η), of having a message ψ or η in the system. Self-consistency equations for these probabilities can be written as follows: Q(η) = E λ,q dη Q(η ) dψ P (ψ)δ(η − G(η , ψ)) , where the product is over the neighbours and the functions F (η, ψ) and G(η, ψ) are given by the equations in Appendix 1. These population dynamics equations can be solved iteratively and once the fixed point has been reached, averages over the RRN ensemble can be directly computed.
In the population dynamics algorithm, we start by initializing the system with a random population of messages and by fixing the parameters of the RRN, q and λ. Then we iterate using the equations of Appendix 1 where the neighbours are extracted at random, using the distributions (1) and (2). This is done until convergence of the mean of the messages in the system. At convergence we can compute the mean value of the metabolites, µ , and reactions, ν , with respect to the ensemble of RRNs and over the measure (13).