6468
Ind. Eng. Chem. Res. 2006, 45, 6468-6476
Wiring Diagrams for Complex Reaction Networks Ilie Fishtik, Caitlin A. Callaghan, and Ravindra Datta* Department of Chemical Engineering, Fuel Cell Center, Worcester Polytechnic Institute, Massachusetts 01609
A graph-theoretic approach to developing reaction route (RR) networks is described that is a powerful new methodology for investigating the integrated architecture and dynamics of complex chemical and biological reaction systems. A distinct feature of our RR graphs is that while the edges represent reaction steps comprising the mechanism, the nodes simply represent their interconnectivity (not necessarily individual species, as has been the convention so far). As a result, all pathways can be traced simply as trails or walks on the RR graphs. Further, these networks are directly converted into wiring diagrams, allowing the use of electrical circuit theory, including Kirchhoff’s current and voltage laws. Thus, they not only provide the reaction network topology but can be used to rigorously analyze network dynamics and the dominant pathways and bottleneck steps. We highlight the approach with three illustrations: (1) an enzymatic reaction system, (2) a heterogeneous catalytic process, and (3) an atmospheric chemistry example. The last example involves 41 reaction steps and provides the most comprehensive chart of atmospheric ozone chemistry to date, while the catalysis illustration with 17 steps provides the most complete view of water-gas shift chemistry. 1. Introduction Chemical and biological systems are replete with examples of highly complex reaction networks,1 many involving hundreds of elementary steps and chemical species. Examples include metabolic2-4 and regulatory5 cellular reaction networks, heterogeneous catalytic processes,6 combustion,7 and atmospheric chemistry.8 In their full glory, these networks, or maps, aim to provide a comprehensive look at the web of interconnections among the various chemical species and their reactions, although this has so far been satisfactorily accomplished in only a handful of cases.6,9 Nonetheless, it is abundantly clear that the reductionist’s approach of considering only one, or a portion of a, pathway at a time is no longer adequate. Rather, the full network must be considered to shed light on the relative significance of the various pathways and their interactions as, for instance, when investigating the effects of a new drug, or one’s genetic predilection.1 Thus, rapid advances are being made in network chemistry and biology,1,2-5,9-16 which has become an important research focus in the data-rich postgenomic era.1,14,6 However, the field is still in its infancy,18 and there is a need to clearly enunciate the fundamental characteristics and organizing principles of such networks, so that they are useful not just as a visualization tool in comprehending the myriad interactions but also in quantitative analysis. A knowledge of both the network structure and its dynamics is the ultimate goal. However, these are independent of each other. This is significant, since the kinetics of elementary steps involved are not yet well-known for most complex systems, with some exceptions,7 although ab initio modeling17 offers hope. The network architecture can, however, be obtained simply from its chemistry19 and is an exceedingly useful tool in its own right. In common vernacular, “reaction networks” represent the graphical depiction of underlying chemistry and reaction pathways as, e.g., in wall charts showing cellular metabolism, ozone chemistry, etc. that grace the hallways of academia. In these, the reaction intermediates, or metabolites, are individually drawn at hubs, or nodes, interconnected via arrows representing their reactions. However, typically, only simplified versions of * To whom correspondence should be addressed. E-mail: rdatta@ wpi.edu.
the full chemistry are represented, since depicting all the various interactions produces jumbled schematics. Also, there are no accepted guidelines for such drawings; these decisions are left largely to the flair of the graphic artist. While a useful teaching tool, these charts are not useful for quantitative network analysis. Graph theory is being increasingly used for translating such renderings into graphs with quantifiable characteristics.5,20-22 Unfortunately, the resulting “reaction graphs”20,21 also adhere to the universal practice of depicting reaction intermediates indiVidually as nodes, while the edges represent chemical reactions. This limits their quantitative application to linear reaction networks.22 Nonetheless, a number of advances have been made, including elucidation of some key graph-theoretic structural aspects of metabolic networks,1,15,16 for instance, their scale-free (rather than single-scale or random) topology, i.e., some molecules such as pyruvate form large hubs, while most species undergo fewer reactions. Metabolic networks are also found to have small-world architecture,11 i.e., most species are connected via a small number of steps and possess a hierarchical network organization. However, a comprehensive approach to depicting and analyzing complex reaction networks is still distant.18 We have recently developed a graph-theoretic approach23-25 that promises to advance the network approach to biology and chemistry. The purpose of this paper is to demonstrate its characteristics and utility in different areas. The resulting “reaction route” graphs, or networks, are distinct from conventional reaction graphs in that the nodes do not necessarily represent a given species but, rather, they simply show how the elementary reaction steps, represented as graph edges, are interconnected to depict the various reaction routes (RRs) or reaction paths, as trails, or walks, on the graphs. Further, the RR graphs can be directly converted into equivalent electrical circuits, or wiring diagrams, wherein each reaction branch represents an appropriate electrical element, e.g., resistance, and reaction rate is akin to branch current and reaction affinity, or Gibbs free energy change, is like branch voltage, thus allowing well-developed methods of electric circuit theory, including Kirchhoff’s laws, to be utilized for a rigorous network analysis.
10.1021/ie050814u CCC: $33.50 © 2006 American Chemical Society Published on Web 08/15/2006
Ind. Eng. Chem. Res., Vol. 45, No. 19, 2006 6469
Figure 1. Conventional (A) and RR graph (B) for enzyme example.
2. Network Topology from Mechanism We describe, for simplicity, the topological construction of an RR network for a reaction system with a single overall n reaction (OR), represented by ∑i)1 νiTi ) 0, where νi is the stoichiometric coefficient of terminal species Ti (i ) 1, 2, ..., n), i.e., the reactants and products of the OR. As an example, we consider the conversion of 7,8-dihydrofolate (A) and NADPH (B) into 5,6,7,8-tetrahydrofolate (C) and NADP+ (D), catalyzed by the enzyme dihydrofolate reductase (E). This system has been well-investigated experimentally and theoretically24,26-28 and is represented by the OR:
-A - B + C + D ) 0 or A + B ) C + D The OR is, of course, a culmination of the various elementary reaction steps sF (F ) 1, 2, ..., p) comprising its mechanism. These steps involve both terminal, Ti, and intermediate species, Ik (k ) 1, 2, ..., l), i.e., those species that do not appear in OR and may be written in the generic form l
s F:
∑
k)1
n
RFkIk +
βFiTi ) 0 ∑ i)1
(F ) 1, 2, ..., p)
(1)
where the stoichiometric coefficient of intermediate species Ik in reaction sF is RFk, while that for the terminal species Ti is βFi. For the enzyme example, the (p ) 13) reaction steps sF are listed in Figure 1. The intermediate species are E, EA, EB, EC, ED, EAB, ECD, EBC, and EAD (l ) 9), while the terminal species are A, B, C, and D (n ) 4). Only q ) 8 of the l intermediates are independent, however, because of the enzyme balance. 2.1. Reaction Routes. The first step is to determine the various possible RRs in the reaction system. In fact, there are two types of RRs of interest: (1) those that result in the OR and are called full routes (FRs), and (2) those that produce no reaction or a zero OR (i.e., all stoichiometric coefficients are zero) and are, hence, called empty routes (ERs).25 For the FRs, since only terminal species are present in the OR, this is akin to finding all of the various linear combinations of steps sF that eliminate all the intermediates. In other words, the gth full reaction route (FRg) is a linear combination of the given elementary reaction steps sF that eliminates all the intermediates and provides the OR
6470
Ind. Eng. Chem. Res., Vol. 45, No. 19, 2006 p
FRg:
∑σgFsF ) OR F)1
(2)
where σgF is called the stoichiometric number29 of reaction step sF in the gth RR. As an illustration, for the enzyme example (Figure 1), FR2: (+1)s3 + (+1)s4 + (+1)s5 + (+1)s6 + (+1)s7 ) OR, the stoichiometric numbers for the other steps being equal to zero, i.e.,
σ2F s3: E + A ) EA
1
s4: EA+ B ) EAB
1
s5: EAB ) ECD
1
s7: ECD ) EC + D 1 s6: EC ) E + C
1
OR: A + B ) C + D This FR may be also presented graphically as a trail or a walk
the given mechanism. However, this results in infinite possibilities if no further constraints are stipulated, since the number of equations is less than the number of unknowns. This nonuniqueness of the RRs is avoided by requiring them to be “direct”,29,31 or minimal, in the sense that if a step is omitted, there is no linear combination of the remaining steps that will produce the OR. Thus, Milner30 defined direct RRs as those that involve a number of elementary steps e rank r + 1 ) q + 1, where r is the matrix of stoichiometric coefficients of the intermediates. All the direct RRs are thus obtained by picking various combinations of q + 1 steps from the given set of p. However, not all the resulting RRs are distinct. Thus, only the set of distinct RRs (FRs and ERs) is retained, forming the RR matrix σ ) [σgF]. These direct RRs are equivalent to the “elementary flux modes” in the metabolic literature,15 of which the so-called “extreme pathways” are an independent subset.14 The set of RRs in σ, while distinct, is not independent, since rank σ ) p - q. Thus, for the enzyme example, p - q ) 13 8 ) 5. As a result, a submatrix of σ comprising a set of linearly independent RRs may be chosen, called the fundamental RR matrix, σf.23 This is, of course, not unique, as different sets of linearly independent RRs in σ may be chosen. For instance, in the enzyme example, we may select the following set: which
In addition, there exist linear combinations of steps that eliminate all of the species, both intermediate and terminal, resulting in ERs p
ERg:
∑σgFsF ) 0 F)1
(3)
An example is ER1 in Figure 1, i.e., (+1)s1 + (+1)s2 + (-1)s4 + (-1)s3 ) 0,
σF s1: E + B ) EB
+1
s2: EB + A ) EAB +1 s4: EA+ B ) EAB
-1
s3: E + A ) EA
-1
Net: 0 ) 0 or
The complete set of RRs (FRs and ERs) for this system is listed in Figure 1 and was enumerated as follows,24 based on the wellestablished stoichiometric theory of reaction routes.29-33 Thus, for a given FRg, the combination of eqs 1 and 2 provides the p RFkσgF ) 0 (k )1, 2, ..., l), which can be set of relations ∑F)1 solved to obtain σgF in terms of RFk, i.e., the stoichiometry of
is, in fact, like selecting a tree for the RR graph with the indicated twigs and links.23,34 For simple systems, a suitable σf may alternately be obtained simply by inspection of the mechanism. However, a systematic development of the complete list of RRs requires use of the stoichiometric algorithm.23,25 With the RRs and an appropriate σf thus in hand, we next construct an RR graph, GR, on which all the RRs can be visualized as trails, i.e., σf, corresponds, in the parlance of conventional graph theory, to the “cycle” matrix of GR.34 2.2. Reaction Route Graphs. This is accomplished by next defining the RR graph, GR, as a connected digraph composed of b branches, or edges, representing the reactions sF as well as the OR, and a set of n nodes, or vertices, representing the interconnectivity of edges or reactions, so that all RRs may be traced as trails, or walks, on the GR. We distinguish between two types of nodes, namely, intermediate and terminal nodes. The intermediate nodes (INs) connect only elementary reactions, while the terminal nodes (TNs) interconnect the OR and elementary reactions. The reactions in GR are assembled so that (a) any cycle is an empty route (ER), (b) any walk between the terminal nodes is a full route (FR), and (c) the RR graph includes, as walks, the complete list of direct FRs and ERs enumerated by the stoichiometric algorithm. The gth walk, or trail, given by wg ) ∑nsfne σgFsF, represents an RR. If the branch sF is oriented along the gth trail, σgF ) +1, and otherwise, σgF ) -1. If a branch is not included in the trail, σgF ) 0. Each reaction branch sF(ni, nj) is incident from node ni and incident to node nj, i.e., it is drawn as a line with an arrow pointing from ni to nj showing its assumed direction, although the reaction may actually proceed in either direction in a given RR. Even though a branch may occur more than once within a GR, they are distinguished by their end nodes, i.e., location. An
Ind. Eng. Chem. Res., Vol. 45, No. 19, 2006 6471
RR graph is minimal, GR,min, if each reaction occurs in it only once, i.e., b ) p + 1 and n ) p - q for GR,min. Our RR graphs depart significantly from the current practice20-22 in that the nodes in RR graphs are not constrained to represent individual species. Rather, nodes are defined based only on reaction connectivity. It turns out, however, that they each can, in fact, be construed to represent either a single or a group of species and their properties, e.g., thermodynamic potential, or their mass balance condition. In fact, we use the quasi-steady-state (QSS) condition of the species to enumerate a finite and unique set of direct INs and TNs, from which one can select an appropriate subset of nodes to construct the RR graph.25 2.3. RR Graph Matrices. Once the graph is drawn, the socalled incidence matrix, M ) [mjF] may be obtained,23 in which the rows correspond to the nodes nj and columns to the branches sF of the GR. The elements of the incidence matrix are mjF ) +1 if the branch sF is incident from node nj, mjF ) -1 if sF is incident to node nj, and mjF ) 0 if the sF is not incident at nj. Thus, the jth node nj is represented by ∑F mjFsF. The number of edges incident on a node nj is its degree d(nj). In the spirit of direct RRs, we define direct nodes as those with the minimum degree.25 The independent submatrix of M is the reduced incidence matrix Mf. For GR,min, rank M ) n - 1. If the elements of the σf and Mf matrices are all +1, -1, or 0, then the mechanism is amenable to representation by a minimal RR graph. This is drawn in the most straightforward manner by starting with a σf which contains one of the shortest FRs, the remaining RRs being among the shortest ERs from the enumerated list. One begins by sketching the FR by drawing all the edges and connecting them via nodes. The sequence of steps in the FR drawn may be chosen so as to be mechanistically meaningful, although the order is inconsequential from a kinetics standpoint. Then, the ERs are added, one at a time, around the existing edges, such that no reaction is repeated. Such a procedure for the enzyme example, starting with σf given by eq 4, readily results in the RR graph shown in Figure 1, where it is also compared to a conventionally drawn reaction graph. The appropriateness of the RR graph is confirmed by checking that all of the enumerated FRs and ERs, as listed in Figure 1, can be traced on this graph as trails and that all the INs and TNs are among the enumerated direct nodes (Figure 1). It may be noted that the nodes in this simple example do, in fact, represent individual intermediates. However, this is not true in general.
Figure 2. Wiring diagram for enzyme example. Table 1. Microkinetic Model for the WGS reaction on a Copper Catalyst s1 s2 s3 s4 s5 s6 s7 s8 s9 s10 s11 s12 s13 s14 s15 s16 s17
CO + S ) COS H2O + S ) H2OS CO2S ) CO2 + S HS + HS ) H2S + S H2S ) H2 + S H2OS + S ) OHS + HS COS + OS ) CO2S + S COS + OHS ) HCOOS + S OHS + S ) OS + HS COS + OHS ) CO2S + HS HCOOS + S ) CO2S + HS HCOOS + OS ) CO2S + OHS H2OS + OS ) 2OHS H2OS + HS ) OHS + H2S OHS + HS ) OS + H2S HCOOS + OHS ) CO2S + H2OS HCOOS + HS ) CO2S + H2S
OR
H2O + CO ) CO2 + H2
the difference between the forward and reverse rates, i.e.,
rF ) b rF - a rF ) B kF
∏i ai-νA
-A kF
∏i aiAν
pi
(5)
where ai is the activity of species i and b VFi and a VFi are the stoichiometric coefficients of species as reactants and products, respectively. The affinity, AF ) ∑i(-V bFiµi) - ∑ia VFiµi, when combined with the chemical potential, µi ) µ0i + RT ln ai, leads to dimensionless affinity AF () AF/RT) for the elementary step in the form
3. Wiring Diagrams and Kirchhoff’s Laws A key advantage of the RR graphs obtained as above is that they can be directly converted into an equivalent electric circuit, or wiring diagram. Thus, each edge in the RR graph may be replaced by its impedance or “resistance” RF (or by a combination of resistance, capacitance, and inductance, for the general case), while the edge representing the OR is replaced by a “voltage” source, AOR, namely, the affinity of the OR. The resulting wiring diagram for the enzyme example is shown in Figure 2. This allows utilization of the well-established methodologies of circuit analysis,34-36 including Kirchhoff’s laws, in the rigorous analysis of such reaction networks. These conservation laws impose important topological constraints on edge rates and affinities, namely, kinetic and thermodynamic22 constraints. 3.1. Reaction Kinetics as Ohm’s Law. Each branch, sF, is characterized by its affinity, AF, akin to branch voltage, and its net rate, rF, which is analogous to edge current. The net rate is
pi
AF ) ln
b rF a rF
(6)
The reaction resistance, RF, i.e., the inverse of reaction conductance, is defined23 as the mean value of the 1/rF between its limiting values, namely, b rF and a rF
RF ≡
1 b rF - a rF
r ) ln(b r /a
∫arbr r1F drF ) br F -F arFF ) F
F
AF rF
(7)
which, hence, provides a linear relation between AF and rF, in the form of Ohm’s law, albeit RF changes with reaction conditions. 3.2. Kirchhoff’s Current Law (KCL). Assuming the nodes to have a zero capacity, the net rates of reactions or edges (akin to edge current) incident at the node nj (INs and TNs) sum to zero
6472
Ind. Eng. Chem. Res., Vol. 45, No. 19, 2006
Figure 3. RR graph for the WGS reaction on a copper catalyst.
∑F mjFrF ) 0
(j ) 1, 2, ..., N)
(8)
Equation 8 is the result of a mass balance (QSS) around nj and is the equivalent of Kirchoff’s current law (KCL). For example, for the central node in Figure 1, -r1 - r3 + r6 + r9 ) 0. 3.3. Kirchhoff’s Voltage Law (KVL). Being a state function, the sum of reaction affinities, akin to edge voltages, around the gth RR (FR or ER) is zero, i.e.,
∑F σgFAF ) 0
(9)
Combining eqs 5, 6, and 9 furthermore provides
(k BF/k AF)σ ∏ ER g
gF
) 1 and
(k BF/k AF ) σ ∏ FR
gF
g
) KOR
(10)
where KOR is the equilibrium constant for the OR. Thus, KVL represents a thermodynamic constraint on kinetics of steps in an RR, with significant consequences. For instance, from Figure 1, for ER2, i.e., s1 + s6 + s10 - s11 ) 0, the edge affinities follow A1 + A6 + A10 - A11 ) 0 or
( )( )( )( )
B k1 B k6 B k 10 A k 11 )1 A k1 A k6 A k 10 B k 11
(11)
The original data provided in the literature was inconsistent with this statement.24 It should be mentioned that the concept of “resistance distance” in a graph is closely related to RR graphs.37-39 3.5. Overall Kinetics. With the reaction network wiring diagram thus available, one can calculate the overall resistance ROR of the circuit in terms of RF, following the standard electrical circuit theory,35 and then obtain the overall rate from rOR )
Ind. Eng. Chem. Res., Vol. 45, No. 19, 2006 6473
AOR/ROR. Further, by using the KCL or, equivalently, the QSS condition and/or the pseudo-equilibrium hypothesis (PEH) for the intermediates, the intermediate concentrations can be obtained in terms of those of the terminal species. Thus, ROR and rOR can be evaluated. It is noteworthy, that the wiring diagram makes use of the entire network, unlike conventional procedures wherein only a single RR is analyzed at a time and ad hoc assumptions, e.g., rate-determining step (RDS) and PEH, are the norm.40 The alternate approach involves numerical solution of the governing differential equations representing the species mass balance.6,8,41 Furthermore, clear conclusions can be deduced from our approach regarding the importance of the various pathways and rational simplifications can be made based on a rigorous comparison of parallel and series resistances as discussed before by us.24
Table 2. Elementary Reactions Comprising the Mechanism of Ozone Destruction/Formation in the Atmosphere
4. Water-Gas Shift Catalysis A 17-step microkinetic model for the water-gas shift (WGS) reaction on Cu(111) considered earlier by us is presented in Table 1,42 where S denotes a catalyst surface site. This model consists of adsorption (s1 and s2) and desorption (s3 and s5) steps, along with several surface reactions. The rate constants for the elementary reactions were predicted a priori.42,43 A more complete mechanism for this reaction is currently under development and will be reported elsewhere.44 Using the developed stoichiometric algorithm,25 the complete list of FRs, ERs, INs, and TNs were enumerated. It was evident from this list that this is a nonminimal mechanism, as there were a number of FRs with σgF) (2. Thus, the nonminimal RR graph for the WGS reaction was constructed as per the following algorithm, and the procedure and the results are shown in Figure 3. This symmetrical RR graph has all of the RRs enumerated as trails and shows how the various mechanistic steps are hooked up. A step for which σgF ) (2 must occur at more than one location, in accordance with the definition of RRs as trails. In fact, the RR graph then turns out to be symmetrical and involves every elementary step as well as the OR exactly twice. This preserves the rate rF and the affinity AF of a step sF irrespective of location. We start, as before, with a σf matrix containing one of the shortest FRs with unit σgF’s, the remaining being the shortest ERs. The FR is drawn with an appropriate sequence of steps, along with another version in which the sequence of steps, but not their direction, is reversed. The two FRs are placed as subgraphs relative to each other so that they have rotational, or point, symmetry of 180°. Thereafter, the ERs are placed, one at a time, once each in both subgraphs such that no step is repeated within a subgraph. The two subgraphs are then interconnected by fusing the appropriate nodes from each subgraph only when the remaining ERs cannot be contained entirely within a subgraph. Two nodes are said to be fused34 when they are replaced by a single node containing all edges incident on the original nodes. Of course, the fused node must also be direct. When all of the remaining ERs in σf have thus been placed and the resulting INs and TNs are among the enumerated nodes, the RR graph is complete and is checked to ensure that the complete list of RRs can be traced on it as trails. It is also evident that the resulting graph is nonplanar (Figure 3). By comparing the resistances in the various pathways shown in Figure 3, we were able to reduce the network to one comprising 11 elementary reaction steps and only 3 parallel pathways.28 Further, a comparison of the sequential resistances in the three parallel pathways revealed the bottleneck steps
s1 s2 s3 s4 s5 s6 s7 s8 s9 s10 s11 s12 s13 s14 s15 s16 s17 s18 s19 s20 s21 s22 s23 s24 s25 s26 s27 s28 s29 s30 s31 s32 s33 s34 s35 s36 s37 s38 s39 s40 s41
OH + O ) H + O2 H + O2 ) HO2 HO2 + O ) OH + O2 H + O3 ) OH + O2 OH + O3 ) HO2 + O2 HO2 + O3 ) OH + 2O2 H + HO2 ) OH + OH O + H2O ) OH + OH OH + HO2 ) H2O + O2 H + HO2 ) H2O + O OH + H2 ) H2O + H H + HO2 ) H2 + O2 O + H2 ) OH + H OH + H2O2 ) H2O + HO2 HO2 + HO2 ) H2O2 + O2 O + H2O2 ) OH + HO2 O + O 3 ) O 2 + O2 O + O2 ) O 3 O3 + hV ) O + O2 O + O ) O2 O2 + hV ) O + O O3 + hV ) O(1D) + O2 O(1D) + M ) O + M O(1D) + O3 ) O2 + O2 Cl + O3 ) ClO + O2 ClO + O ) Cl + O2 Cl + HO2 ) ClO + OH ClO + hV ) Cl + O NO + O3 ) NO2 + O2 NO2 + O ) NO + O2 NO + HO2 ) NO2 + OH NO2 + hV ) NO + O ClO + NO ) NO2 + Cl Cl + HO2 ) HCl + O2 Cl + H2 ) HCl + H Cl + H2O2 ) HCl + HO2 OH + HCl ) H2O + Cl O + HCl ) Cl + OH O + HOCl ) OH + ClO HOCl + hV ) OH + Cl HO2 + ClO ) HOCl + O2
OR
2O3 ) 3O2
(RDS) in each of the three parallel routes.24 Thus, the remaining steps could be treated via PEH. Then, the following explicit rate expression for WGS reaction results, where the step equilibrium constant KF ) B kF/k AF
[
pH21/2
rOR ) B k 6K2pH2Oθ0 (k B8 + B k 10)K1pCO + B k 15 (K4K5)1/2 2
( (
)] / [ )]
B k6 +B k 15 K6
B k 7K2K5K15pCO B k 7K2K5K15pCO + B k 15pH2
pH21/2 B k 7K2K5K15pCO + B k 7K2K5K15pCO + B k 15pH2 (K K )1/2 4 5
(
k 10)K1pCO 1 (k B8 + B
pCO2pH2 KORpH2OpCO
)
(12)
Where, the free catalyst site fraction
θ0 )
1 1 + K1pH2O + K2pCO + (K4K5)-1/2pH21/2
(13)
which includes the three parallel pathways and predicted rate and equilibrium constants (Figure 3) and agrees with predictions of the full network as well as with the experimental rate data very well. In other words, the WGS RR graph (Figure 3) not only affords an unmatched pictorial representation of the
6474
Ind. Eng. Chem. Res., Vol. 45, No. 19, 2006
Figure 4. RR graph for the mechanism of ozone destruction/formation in the atmosphere.
mechanism and the numerous possible reaction pathways but also allows a predictive quantitative analysis and comparison of the resistances of the various pathways, which, in turn, allows identification of the bottleneck steps and a simplified but explicit predictive rate expression. 5. Atmospheric Ozone Chemistry In addition to the biochemical and catalysis examples above, an illustration is provided in the area of atmospheric chemistry, where the concept of catalytic cycles is widely used.8 These catalytic cycles are mainly drawn on an individual basis. No attempt has so far been made to enumerate a global mechanism comprising the complete set of available catalytic cycles describing the ozone chemistry in the atmosphere. Due to a very large number of elementary reactions implicated in the mechanism of ozone chemistry in the atmosphere, this, of course, is not a simple task, since the number of catalytic cycles dramatically increases with an increase in the number of elementary reactions. On the other hand, enumeration of only individual catalytic cycles without a clear notion of their interaction and dominance is, obviously, not satisfactory. Here, we present a RR network describing the ozone production/destruction process involving 41 elementary reactions, as listed in Table 2. The nonminimal RR network was constructed employing the algorithm described above and the final result is presented in Figure 4. In constructing the network, it was assumed that the ozone production/destruction process is described by the OR: 2O3 ) 3O2. Hence, all the other species were considered intermediate species. From the total of l )14 intermediates (H2O, H2, H2O2, OH, O, H, HO2, O(1D), Cl, ClO, HCl, HOCl, NO, and NO2),
only q ) 11 are linearly independent due to three mass balance conditions for the catalytic species H, Cl, and N. As expected, the intermediate nodes in the RR graph do not necessarily represent individual intermediates. In fact, 7 out of 11 INs represent individual intermediates, while 4 INs represent a linear combination of the intermediates (Figure 4). It may further be seen that the RR network is “doubled” and symmetric, i.e., each intermediate node and each elementary reaction are involved twice in the RR network. In addition to nonunit stoichiometric numbers in the enumerated RRs, this is a consequence of the fact that in order to balance the nodes in the RR network it is necessary to add two ORs. Indeed, from the mechanism (Table 2) and the stoichiometry of the OR as well as the RR graph, it follows that the rate of the OR rOR is expressed via the rates of the elementary reactions as
2rOR ) r4 + r5 + r6 + r17 + r18 + r19 + r22 + r24 + r25 + r29 (14) An inspection of the network reveals that besides a few RRs, or catalytic cycles, that are commonly discussed in the literature in order to explain the ozone destruction/production process in the atmosphere, there is in fact a huge number of alternate RRs. The network also explicitly shows the interrelation among the RRs associated with various subsystems, e.g., the interrelation among HOx, ClOx, NOx, etc. subsystems. It may furthermore be observed that the network is nondirect, that is, the network involves both direct and nondirect RRs. The list of direct RRs from the network is complete in the sense that this list may be enumerated by using the soichiometric RR theory. Here, in fact, we have considered only a subset of elementary reactions implicated in the ozone destruction/production mech-
Ind. Eng. Chem. Res., Vol. 45, No. 19, 2006 6475
anism. Other elementary reactions and species, e.g., involving methane, bromine, fluorine, etc., have been proposed in the ozone chemistry and, in principle, may be added to the network. However, this is a complex task, which is in progress. 6. Conclusions This paper provides the basis and algorithm for drawing a new type of reaction network developed by us called reaction route graph, or RR graph, that not only pictorially provides the complete web of interconnectivity of a complex reaction network but also is amenable to quantitative thermodynamic and kinetic analyses. In our approach, each branch or edge represents a reaction step, while the nodes simply represent branch or reaction connectivity, so that all reaction routes or mechanisms can be traced on it simply as walks. The RR graph can be directly converted into an equivalent electrical network, or wiring diagram, with each elementary step represented by an electrical element, e.g., a resistance, with reaction rate being analogous to current and reaction affinity analogous to voltage. Further, the overall network resistance and rate may be calculated using conventional electrical network analysis, i.e., Kirchhoff’s laws. This allows use of the entire network in the calculation of the reaction rate, rather than a single RR at a time, as is the conventional approach. Furthermore, a direct comparison of the resistances allows simplification of the network and elucidation of the dominant routes. The approach is illustrated here with three different examples: a biochemical example involving an enzyme reaction network, a catalysis example for the case of water-gas shift reaction, currently of considerable interest in the context of hydrogen generation, and an atmospheric chemistry example involving ozone chemistry. It is our hope that the approach will be utilized by researchers to unravel and analyze a wide variety of reaction networks. Nomenclature AF ) affinity of the elementary reaction AOR ) affinity of the overall reaction AF ) dimensionless reaction affinity b ) number of branches in a RR graph Ik ) intermediate species kF ) rate constants of the elementary reactions KF ) equilibrium constants of elementary reactions l ) total number of intermediate species M ) incidence matrix Mf ) reduced incidence matrix mjF ) element of incident matrix n ) total number of terminal species nj ) node n ) number of nodes in RR graph p ) total number of elementary reaction steps q ) independent number of intermediate species rF ) rate of the elementary reaction b rF ) rate of the forward reaction a rF ) rate of the reverse reaction rOR ) rate of the overall reaction RF ) resistances of the elementary reactions ROR ) resistance of the overall reaction sF ) elementary reaction S ) catalyst surface site Ti ) terminal species wg ) walk representing RRg
Greek Letters σgF ) stoichiometric number RFk ) stoichiometric coefficient of the intermediate species in an elementary step βFi ) stoichiometric coefficient of the terminal species in an elementary step νi ) stoichiometric coefficient of terminal species in the overall reaction r ) stoichiometric matrix of the intermediate species σ ) reaction route matrix σf ) fundamental reaction route matrix AbbreViations OR ) overall reaction RR ) reaction route FR ) full route ER ) empty route IN ) intermediate node TN ) terminal node QSS ) quasi-steady-state WGS ) water-gas shift PEH ) pseudo-equilibrium hypothesis RDS ) rate-determining step KCL ) Kirchhoff’s current law KVL ) Kirchhoff’s voltage law Literature Cited (1) Baraba´si, A. L. Linked. The New Science of Networks; Perseus Pub.: Cambridge, MA, 2002. (2) Stephanopoulos, G. N.; Aristidou, A. A.; Nielsen, J. Metabolic Engineering. Principles and Methodologies; Academic Press: San Diego, CA, 1998. (3) Jamshidi, N.; Edwards, J. S.; Fahland, T.; Church, G. M.; Palsson, B. O. Dynamic Simulation of the Human Red Blood Cell Metabolic Network. Bioinformatics 2001, 17, 286. (4) Delgado, J.; Liao, J. C. Identifying Constraints on Bioreaction Systems. Chem. Eng. Sci. 1992, 47, 1277. (5) Fell, D. A. Understanding the Control of Metabolism; Portland, London, 1997. (6) Stoltze, P. Microkinetic Simulation of Catalytic Reactions. Prog. Surf. Sci. 2000, 65, 65. (7) Peters, N.; Rogg, B. Reduced Kinetic Mechanisms for Applications in Combustion Systems; Springer: New York, 1993. (8) Brasseur, G. T.; Orlando, J. P.; Tyndall, G. S. Atmospheric Chemistry and Global Change; Oxford University Press: New York, 1999. (9) Cakir, T.; Tacer, C. S.; Uelgen, K. O. Metabolic Pathway Analysis of Enzyme-Deficient Human Red Blood Cells. BioSystems 2004, 78, 49. (10) Watts, D. J.; Strogatz, S. H. Collective Dynamics of ‘Small-World’ Networks. Nature 1998, 393, 440. (11) Wagner, A.; Fell, D. A. The Small World Inside Large Metabolic Networks. Proc. Biol. Sci. 2001, 268, 1803. (12) Stephanopoulos, G.; Alper, H.; Moxley, J. Exploiting Biological Complexity for Strain Improvement Through Systems Biology. Nat. Biotechnol. 2004, 22, 1261. (13) Kell, D. B. Metabolomics and Systems Biology: Making Sense of the Soup. Curr. Opin. Microbiol. 2004, 7, 296. (14) Papin, J. A.; Price, N. D.; Wiback, S. J.; Fell, D. A.; Palsson, B. O. Metabolic Pathwaysin the Post-Genome Era. Trends Biochem. Sci. 2003, 28, 250. (15) Schuster, S.; Fell, D. A.; Dandekar, T. A General Definition of Metabolic Pathways Useful for Systematic Organization and Analysis of Complex Metabolic Networks. Nat. Biotechnol. 2000, 18, 326. (16) Barabasi, A. L.; Oltvai, Z. N. Network Biology: Understanding the Cell’s Functional Organization. Nat. ReV. Genet. 2004, 5, 101. (17) Broadbelt, L. J.; Snurr, R. Q. Applications of Molecular Modeling in Heterogeneous Catalysis Research. Appl. Catal. A, Gen. 2000, 200, 23. (18) Bray, D. Molecular Networks: The Top-Down View. Science 2003, 301, 1864. (19) Schuster, S.; Hilgetag, C.; Woods, J. H.; Fell, D. A. Reaction Routes in Biochemical Reaction Systems: Algebraic Properties, Validated Calculation Procedure and Example from Nucleotide Metabolism. J. Math. Biol. 2002, 45, 153.
6476
Ind. Eng. Chem. Res., Vol. 45, No. 19, 2006
(20) Balaban, A. T. In Reaction Graphs; Bonchev, D., Mekenyan, O., Eds.; Kluwer Academic Publishers: Dordrecht, The Netherlands, 1994; p 137. (21) Temkin, O. N.; Zeigarnik, A. V.; Bonchev, D. G. Chemical Reaction Networks: A Graph-Theoretical Approach; CRC Press: New York, 1996. (22) Beard, D. A.; Babson, E.; Curtis, E.; Qian, H. Thermodynamic Constraints for Biochemical Networks. J. Theor. Biol. 2004, 228, 327. (23) Fishtik, I.; Callaghan, C. A.; Datta, R. Reaction Route Graphs. I. Theory and Algorithm. J. Phys. Chem. B 2004, 108, 5671. (24) Fishtik, I.; Callaghan, C. A.; Datta, R. Reaction Route Graphs. II. Examples of Enzyme- and Surface-Catalyzed Single Overall Reactions. J. Phys. Chem. B 2004, 108, 5683. (25) Fishtik, I.; Callaghan, C. A.; Datta, R. Reaction Route Graphs. III. Non-Minimal Kinetic Mechanisms. J. Phys. Chem. B 2005, 109, 2710. (26) King, E. L.; Altman, C. A Schematic Method of Deriving the Rate Laws for Enzyme-Catalyzed Reactions. J. Phys. Chem. 1956, 60, 1375. (27) Fierke, C. A.; Johnson, K. A.; Benkovic, S. J. Construction and Evaluation of the Kinetic Scheme Associated with Dihydrofolate Reductase from Escherichia Coli. Biochemistry 1987, 26, 4085. (28) Happel, J.; Otarod, M. New Treatment of Enzyme Kinetics Applied to Human Dihydrofolate Reductase. J. Phys. Chem. B 2000, 104, 5209. (29) Horiuti, J. Theory of Reaction Rates as Based on the Stoichiometric Number Concept. Ann. N. Y. Acad. Sci. 1973, 213, 5. (30) Milner, P. C. The Possible Mechanisms of Complex Reactions Involving Consecutive Steps. J. Electrochem. Soc. 1964, 111, 228. (31) Temkin, M. I. The Kinetics of Some Industrial Heterogeneous Catalytic Reactions. AdV. Catal. 1979, 28, 173. (32) Happel, J.; Sellers, P. H. Analysis of the Possible Mechanisms for a Catalytic Reaction System. AdV. Catal. 1983, 32, 273. (33) Fishtik, I.; Datta, R. A Thermodynamic Approach to the Systematic Elucidation of Unique Reaction Routes in Catalytic Reactions. Chem. Eng. Sci. 2000, 55, 4029. (34) Deo, N. Graph Theory with Applications to Engineering and Computer Science; Prentice Hall: Englewood Cliffs, NJ, 1974.
(35) Balabanian, N.; Bickart, T. A. Electrical Network Theory; John Wiley: New York, 1969. (36) Oster, G. F.; Perelson, A. S.; Katchalsky, A. Network Thermodynamics: Dynamic Modelling of Biophysical Systems. Quart. ReV. Biophys. 1973, 6, 1. (37) Klein, D. J.; Randic, M. Resistance Distance. J. Math. Chem. 1993, 12, 81. (38) Balaban, A. T.; Klein, D. J. (Co-Authorship) Erdo¨s Numbers, and Resistance Distances in Graphs. Scientometrics 2002, 55, 59. (39) Bonchev, D.; Balaban, A. Liu, T. X.; Klein, D. J. Molecular Cyclicity and Centricity of Polycyclic Graphs. Part 1. Cyclicity based on Resistance Distances and Reciprocal Distances. Int. J. Quantum Chem. 1994, 50, 1. (40) Boudart, M.; Djega-Mariadassou, D. Kinetics of Heterogeneous Catalytic Reactions; Princeton University Press: Princeton, NJ, 1984. (41) Dumesic, J. A.; Rudd, D. F.; Aparicio, L. M.; Rekoske, J. E.; Trevino, A. A. The Microkinetics of Heterogeneous Catalysis; ACS: Washington, DC, 1993. (42) Callaghan, C.; Fishtik, I.; Datta, R.; Carpenter, M.; Chmielewski, M.; Lugo, A. An Improved Microkinetic Model for the Water Gas Shift Reaction on Copper. Surf. Sci. 2003, 541, 21. (43) Shustorovich, E.; Sellers, H. The UBI-QEP Method: A Practical Theoretical Approach to Understanding Chemistry on Transition Metal Surfaces. Surf. Sci. Rep. 1998, 31, 1-119. (44) Fishtik, I.; Callaghan, C. A.; Datta, R., in preparation.
ReceiVed for reView July 11, 2005 ReVised manuscript receiVed February 22, 2006 Accepted June 26, 2006 IE050814U
