Morphisms of reaction networks that couple structure to function
 Luca Cardelli^{1, 2}Email author
DOI: 10.1186/17520509884
© Cardelli; licensee BioMed Central Ltd. 2014
Received: 25 March 2014
Accepted: 4 July 2014
Published: 15 August 2014
Abstract
Background
The mechanisms underlying complex biological systems are routinely represented as networks. Network kinetics is widely studied, and so is the connection between network structure and behavior. However, similarity of mechanism is better revealed by relationships between network structures.
Results
We define morphisms (mappings) between reaction networks that establish structural connections between them. Some morphisms imply kinetic similarity, and yet their properties can be checked statically on the structure of the networks. In particular we can determine statically that a complex network will emulate a simpler network: it will reproduce its kinetics for all corresponding choices of reaction rates and initial conditions. We use this property to relate the kinetics of many common biological networks of different sizes, also relating them to a fundamental population algorithm.
Conclusions
Structural similarity between reaction networks can be revealed by network morphisms, elucidating mechanistic and functional aspects of complex networks in terms of simpler networks.
Keywords
Morphisms Chemical reaction networks Influence networks Biological networksBackground
Chemical reaction networks
Chemical reaction networks provide a compact language for describing complex dynamical systems of the kind found in inorganic chemistry, biochemistry, and systems biology. They can be presented as certain graphs or as lists of reactions over a set of species. Unlike general formulations of dynamical systems in terms of differential equations, reaction networks explicitly represent mechanism: they present the algorithms that produce certain behaviors by a description of molecular interactions. Implicit in the simple syntax of chemical reactions are (depending on circumstances) stochastic or deterministic kinetic laws that can be used to determine the evolution of systems over time. Unravelling the exact behavior of chemical systems from the kinetic laws can be in general quite demanding; hence, attention has been dedicated to identifying functional properties of reaction networks from their structure or motifs, including questions of multistability and oscillation, and methods for transferring properties of a network to a reduced network (a vast area including[1–11]).
The aforementioned literature is focused on properties of individual reaction networks or their subnetworks. Another way to try to understand the properties of a network is to relate it to another network, perhaps a better known one, either by comparing graph structures[12–15], or more deeply by preserving kinetic features[3, 10]. In this work, we identify kinetic relationships between networks that arise from network mappings, or morphisms. In particular, we explore the notion of network emulation, which allows a complex network to behave like one or more simpler networks. We show how that relationship can be determined from structural properties alone, and how it can be used to transfer system properties. As an application we obtain analytical justification of empirical relationships that have been observed in conjunction with cell cycle switch models[16].
Influence networks
Consider the simple influence network MI in Figure 1(A), where y activates itself and inhibits z, and where conversely z activates itself and inhibits y. It should be fairly intuitive that y and z are competing for dominance, and if y is ever able to fully activate itself and fully inhibit z, then z is forever inhibited, and vice versa (subject to a suitable reaction kinetics). Mutual inhibition networks arise in may areas of biology[18, 20–24]; not all are this simple, and not all are reducible to the particular MI mutual inhibition pattern, but many are routinely summarized in this fashion.
The function of the network QI in Figure 1(B), however, is much harder to interpret: in MI y and z are in mutual inhibition, while in QI we have a kind of quad inhibition. We use QI as an example of a more complicated network of the kind that could occur in biology and whose function might not be obvious (it is in fact a simplified version of[25]). In QI, the y, z species seems to interact in a similar pattern as in MI; for example z is still (indirectly) activating itself and inhibiting y, and conversely. The network structure is thus suggestive of similar functionality, and one could ask whether MI and QI are in fact functionally related.
To know for sure, we need to ask if there is a similarity between the kinetics of the two networks. This question can be typically investigated by studying the kinetic equations of the networks, which can be obtained from their sets of chemical reactions. The approach we take here is instead to study the structure of the networks themselves without, at first, apparent knowledge of their kinetics. In particular we look at morphisms (mappings) between chemical reaction networks, including influence networks such as these. We show that these morphism can characterize functional properties and provide an explanation of kinetic similarity based on structural similarity.
Results and discussion
Mass action interpretation of influence networks
A chemical reaction network is given by a set of irreversible reactions R over a set of species S. Each reaction is written ρ → ^{ k } π, where ρ are the reagents, k > 0 is the rate constant (we assume mass action kinetics), and π are the products. Both ρ and π assign a stoichiometric number to each species in S. For example the reaction 2A + B → ^{ k } B + C, has reactant stoichiometric number 2 for species A, 1 for species B, and 0 for species C, hence ρ _{A} = 2, ρ _{B} = 1, ρ _{C} = 0, and similarly for π on the products side; ρ and π and are called complexes.
Common characterizations of influence networks use sigmoidal (Hill or Reinitz) functions for the transitions between activated and inhibited states[19]. Here we choose an interpretation based on mass action kinetics, which results in a Hill function. Each influence species x (Figure 2, left) is modeled as a triplet of chemical species, denoted x _{0} , x _{1} , x _{2}, connected in a triplet motif of four catalytic chemical reactions (Figure 2, center) with associated rates. The activated and inhibited states of an influence species are represented by separate chemical species: by convention x _{0} is the species whose concentration is high when x is activated and x _{2} is the species whose concentration is high when x is inhibited; x _{1} is an intermediary species that introduces nonlinearity in the transition, and is never otherwise connected to the rest of the network. If, for example, a species i is connected to the inhibition input, then the transition from x _{0} to x _{1} (arrow with hollow circle on top) represents the catalytic reaction${x}_{0}+i{\to}^{{k}_{01}}\phantom{\rule{0.12em}{0ex}}i+{x}_{1}$. A duality is introduced by defining ~x as the species such that ~x _{0} = x _{2}, ~x _{1} = x _{1}, and ~x _{2} = x _{0} (Figure 2, right).
Each influence node x therefore corresponds to a motif of chemical reactions, and it is thus possible to translate any influence network into a chemical reaction network. Solving the mass action equations for the triplet motif at steady state yields a generalized Hill function of coefficient 2. Depending on the four reaction rates, this motif is capable of producing a range of quasilinear, quasihyperbolic, and sigmoidal activation and inhibition responses (see Additional file1).
In summary, we interpret influence networks as an unambiguous, restricted class of mass action chemical reaction networks, in a way that is not too far from common practice because it is based on an explicit mechanism that yields Hill functions. The results that we derive in Methods hold for arbitrary chemical reaction networks, and hence for all influence networks. In the main body of this paper, for exposition purposes, we concentrate on examples of influence networks, while in Additional file2 we provide many examples of small reaction networks. The class of chemical reaction networks that we admit consists of finite sets of irreversible reactions over finite sets of species, where the reaction rates are constants and are interpreted with mass action kinetics. There are no further restrictions: in order to model open chemical systems and systems out of equilibrium, there are no assumptions about conservation of mass or energy, or about detailed balance. Note that even a chemical reaction network of this general kind can be systematically physically realized with DNA nanotechnology[26, 27].
Network emulation
Consider the mapping of species and reactions from QI to MI described in Figure 1 according to equal colors. It will not be clear at this point how this mapping was chosen, and obviously different ones are possible. It satisfies certain properties that will be clarified later, but for now we are just interested in observing some of its consequences. The s and y species of QI are mapped to the y species of MI, and similarly r and z are mapped to z. The mapping of reactions is in this case straightforward: any reaction between species of QI is mapped to a similar reaction between the corresponding species of MI according to the species mapping just described; this is called a homomorphism. As a result, for each reaction of MI we have two reactions of QI that map to it.
We have thus given a mapping between networks based on their structure, and we can now observe a surprising phenomenon about their kinetics. In Figure 1(C) we perform a numerical simulation of the kinetics of MI: there are 6 trajectories, 3 for y (the triplet y _{0}, y _{1}, y _{2}) and 3 for z, as described in the previous section. We have chosen essentially random parameters for reaction rates and initial concentrations, and the 6 trajectories are all distinct. The MI system is inherently bistable, hence (except in degenerate cases) each trajectory converges to a maximum or a minimum. In Figure 1(D) we then simulate the QI network with a particular choice of parameters, and it appears to produce identical kinetics as MI. Actually, there are twice as many traces (and species) in QI: they exactly overlap in pairs, with each pair retracing a corresponding trajectory of MI.
The surprising fact is that, in general, QI can emulate MI (retrace all its trajectories) for any choice of parameters of MI (rates and initial concentrations). The parameters of QI that achieve this matching performance can be systematically extracted from those of MI and from the mapping of species and reactions of Figure 1(AB). In this example, the parameter selection is straightforward, although it can be a bit more complex in general. The initial concentration of each species of QI is taken equal to the initial concentration of the corresponding species of MI under the species mapping, and the rate of each reaction of QI is taken equal to the rate of the corresponding reaction of MI under the reaction mapping. As a cautionary note, this network mapping is not a instance of coarsegraining, at least in the sense that we do not take sums of concentrations.
The fact that QI emulates MI does not preclude QI from having a richer set of behaviors outside of the initial conditions that can be derived from MI: QI does in fact have more degrees of freedom. Still, successful emulation expresses the fact that QI can be regarded as a more complicated version of MI, giving at least a partial insight in its kinetics. And if QI can emulate several unrelated networks, then multiple ‘facets’ of QI can be revealed.
The kinetic emulation property must obviously be a consequence of the kinetic equations of MI and QI, but this is not (directly) what we do here. Instead, we establish the emulation property as a consequence of the existence of a structural morphism between the networks. Network emulation has a number of consequences, which we expand on in the Conclusions. But one should already be evident: a nontrivial but purely structural mapping between networks somehow guarantees that one network can exactly, kinetically, replace another network in all possible circumstances.
Emulation among antagonistic networks
First (although this is stronger than necessary for emulation), all the morphisms in Figure 3 are homomorphic projections: they are obtained by collapsing certain species (as indicated under the arrows) onto species of the target network, and by letting reactions correspond according to the species mapping. In some cases we need to dualize the nodes: for example, in the morphisms leading from MI to AM we collapse ~y onto x, meaning that we map y _{2} onto x _{0}, etc.; see Additional file3 for some detailed network mappings.
Second, and most important, all the morphisms satisfy a stoichiometric relationship (stoichiomorphism), discussed later in this section and in Methods, that can be computed over the stoichiometric matrices and rate constants of the networks irrespectively of initial conditions.
AM, AMr, AMs
The Approximate Majority (AM) network[16] is the quintessential population switch, with asymptotically optimal convergence speed to one of two stable steady states: convergence is achieved in O(log n) with high probability, where n is the number of molecules in a stochastic interpretation of the kinetics[30, 31](A.4)]. Moreover, the steady states are robust to large perturbations, and they are reached quickly even in metastable conditions[30]. AMr and AMs can be generated by separately introducing indirections in the autocatalytic AM reactions, converging towards the characteristic core of a cell cycle switch network, CCR.
The exact interaction pattern of AM can be found in epigenetic switches (Figure 4(A)), where DNA histones can be in one of three states: (M)ethylated, U(nmodified), or (A)cetylated. A contiguous stretch of DNA consists of a population of histones that should be uniformly methylated or acetylated. This is achieved by the M and A states activating two proteins each that catalyze transitions between M,U,A states through the whole population. In Figure 4(A) we have expanded the AM network from Figure 3 into triplets, so that three chemical species are visible and labeled M,U,A. The resulting autocatalytic network reproduces Figure one from[28] in our notation. The known kinetics of AM implies robust uniform settling of the whole histone population into either M or A states, which is the conclusion reached in[28].
MI
The MI network is the basic mutual inhibition network discussed in Background along with QI; its pattern can be found in many biological networks, at least in simplified form[18, 20, 23, 24]. In genetic toggle switches, for example, the selfactivation loops are usually replaced by inducers or constitutive transcription. The morphisms from QI to MI and from MI to AM are detailed in Additional file3.
SI
The SI network is another mutual inhibition network between two species, but with a different algorithm than MI. Two antagonists, instead of promoting themselves, are doubly active in opposing their antagonist. The SI network has exactly the same steady states as AM, while MI has an additional class of unstable steady states (see Additional file4). QI morphs into SI, but by a less obvious mapping than into MI.
In Figure 4(B) we expand the SI network from Figure 3 into triplets. The resulting network largely matches Figure one A from[21], which is a septation initiation network: the ellipses represent the old and new spindle pole bodies that separate, and the other species are in the cytosol. Differences from[21] include the grey links, which are missing in a minimal model, or replaced by other mechanisms in more detailed models.
CCR, GW, NCC
These networks are related to the G2M cell cycle switch[32]. Their common structure consists of the two righthandside feedback loops around the species x,s,r (or z,s,r), where x (or z) is the Cdk1 protein: a cyclindependent kinase that has an essential role in the progress of the cell cycle. Those two loops by themselves give a close but imperfect match to AM kinetics[16], and do not technically emulate AM.
Since the basic feedback loops do not achieve optimal performance, it was suggested in[16] to consider additional known feedbacks involving the Greatwall kinase (Gwl): this gives the GW network, which was shown to perform better in simulations. We can now show analytically that GW emulates AM, and hence that it has its same switching performance. Moreover, it was independently shown[29] that the Greatwall kinase is in fact necessary for the proper biological function of the cell cycle switch. Figure 4(C) reproduces the influence network and protein assignment from[29], Figure seven: apart from an extra feedback around PPA2 that is necessary to reset the switch, that is exactly the GW network. The CCR network is a simplified version of GW where some feedbacks are shortcircuited: it too emulates AM.
Figure 4(D) shows the influence network from[25], Figure three, again in our notation (a species S that is missing here just represents downstream targets of Cdk1). This is the NCC network, which has been proposed as a more complete model of the cell cycle switch, refining for example the interactions of GW. Even this rather complex network can exactly emulate AM. Moreover, the influence interactions are modeled in[25] by phosphorylation/dephosphorylation dynamics, therefore compatibly with our triplet interpretation.
NCC is a highly symmetric network, and it can emulate the equally symmetrical QI, and through it also CCR and AM. Note however that symmetry is not necessary to achieve emulation: GW is not nearly as symmetrical, and neither are AMr and AMs. It is also possible to go from NCC to QI in two steps, resulting in two less symmetrical intermediate networks before symmetry is restored (only the required collapsing of species is indicated).
NCC and GW disagree as models of the cell cycle switch (they do not emulate each other), but they agree on basic functionality. They can both emulate MI, which embodies the essence of mutual inhibition, and indirectly also AM, which embodies the essence of fast switching. Therefore, even through biological uncertainties that may be reflected in conflicting models, and even through the simplifying assumptions of modeling, we can mathematically justify what is believed to be the functional kernel of these networks. This kind of insight was used to determine that the basic cell cycle model in[16, 32] may be missing some feedbacks, because it does not emulate AM and thus does not have optimal performance.
QI
The QI network that we discussed in Background has at least two ‘facets’: it can emulate both CCR and MI, while those networks cannot emulate each other. Through CCR, QI can also emulate SI, but again MI and SI cannot emulate each other. We have no direct biological analog to offer for QI. But it can be seen as a more symmetric variation on GW where the antagonism from z to y is carried out indirectly through s and r, or else as a version of MI where selfactivation is replaced by mutual activation: these are all options that are biochemically available. Both QI and GW can emulate MI, but not each other.
DN
The DN network is a schematic DeltaNotch configuration between two neighboring cells (top and bottom halves). The tight coupling of each two nodes is due to, e.g., lowNotch (s) inducing highDelta (z) in the same cell (top half) and highNotch conversely inducing lowDelta because of degradation[33]. Although the basic functionality of DeltaNotch is well represented, this is an example where there is no close match with models from the literature, which all have differences in detail and miss some of the interactions. In these instances one can follow the empirical path from[16], to see whether the more realistic models still approximate the behavior of DN and AM.
SCR, CCR’, SCR’
These are variations on other small networks, indicating that a rather large set of possibilities exists in network connectivity even after fixing the kinetics of the species. The SCR network is a version of AM where the two direct feedbacks of x onto itself are replaced by indirect feedbacks trough s and r. The point here is that indirections, which are common in biological networks, can (sometimes) be introduced while maintaining the emulation property. The CCR’ and SCR’ networks are just variations in the connectivity of CCR and SCR respectively, with the same number of species and the same ability to emulate AM.
Relationships
Not all networks, even closely related ones, are connected by emulations. For example in reference to Figure 3, we can go from NCC to QI and from QI to SI in two separate steps each, but we cannot go from QI to MI in two steps. (There is no suitable morphism from CCR to MI. E.g., if we map x to z there is no target for the reactions from x _{2}. If we merge s with y in QI and we take the homomorphic projection, we obtain a network unrelated to any in the figure.) Still, if we start from QI and we merge s with y and r with z at once, then the homomorphic projection is an emulation that leads directly to MI. This shows that there may be gradual ways to connect two networks via emulations, or there may be ‘jumps’ in complexity that, for example, would not be discovered by algorithms that attempt to identify a pair of nodes at a time. And there may be no way to relate two networks except through a common simpler network. The question of how to reach a network from another network through a sequence of emulations, is essentially the questions of how complex networks can arise from simple networks without upsetting functionality.
Discussion: how common, useful, and brittle are network emulations?
Although we have shown many meaningful emulations morphisms, this is not the same as saying that such morphisms are common: how can we find them? The networks in Figure 3 all involve two or more antagonistic species, so that certain symmetries become apparent, and quick simulations can be used to falsify suspected symmetries. Emulations that work for some parameters can be extended to other parameters (see Methods), hence a simple strategy is to first test potential emulations with unitrate parameters and heterogeneous initial conditions. This heuristic is imperfect, but is sufficient for all the networks in Figure 3, which are also all homomorphisms (the simplest kind of morphisms). It would seem feasible to use a tool to check all possible unitrate homomorphisms between two networks, since we can perform emulation checks by simple matrix calculations rather than simulations (see Methods, and examples of such calculations in Additional file3). More subtle heuristics could be devised for more subtle morphisms, and more principled algorithms for network matching could also be studied, similar to the techniques for the analysis of transition systems (see[34] and the related work referenced in[10]). It is not clear at this point how such techniques would fare on largescale biological networks, and some notion of approximate matching would likely be required[13–15]. The examples in Additional file2, at the level of chemical reaction networks, give some further hints about what networks may or may not be related by emulations.
Once found, the simple existence of a network emulation can reveal biologicallyrelevant properties of complex networks, by deriving them from known properties of simpler networks. For example, a few of the networks in Figure 3 are related to cell cycle switches, and it has been shown that the speed of cell cycle switching is important to avoid genetic instability during replication[24]. By the existence of those emulations, we immediately obtain that cell cycles switches are capable of achieving asymptotically optimal switching performance, because that is a known property of AM trajectories. This is a nontrivial consequence of emulation that speaks about the kinetics (speed of stabilization), robustness (stability of steady states), and reliability (likelihood of metastable states) of the networks, rather than just their steady state landscape. These are consequences that would be arduous to derive analytically for the larger networks, but through emulation morphisms we can take advantage of the fact that they have been derived analytically for AM[30].We have defined the notion of emulation as exact matching of trajectories between two networks, and in our examples we have chosen network patterns that produce exact kinetic emulations. This way we can be mathematically definite about the nature of the relationship, and we can derive a theory and computational methods for such morphisms. However, exact matching of trajectories is not likely to happen in practice. Even if the reaction rates were to cooperate, realistic networks always have uncertainties and minor details in network connectivity that would cause some divergence: we have discussed examples in conjunction with Figure 4 of more or less exact correspondences with biological networks. Even the notion of activation and inhibition, although widely used, is itself an approximation of underlying mechanisms that vary from case to case. And even if the true networks enjoyed a (deterministic) emulation property, stochastic mechanisms could still differentiate them.
The reality is that exact network emulation is uncommon and susceptible to perturbations: to work perfectly, it requires some adequate amount of mathematical abstraction and deviation from exhaustive biological detail. Nonetheless, the idea can be applied to imperfect situations: if the true networks do not deviate too much from the ideal networks, it is reasonable to expect that they retain their fundamental features, including the emulation relationships. This assumption should of course be tested, either with a theory of approximate network emulation that can take perturbations into account, or in absence of it, by validating approximate emulations in each particular case, for example by simulations. The latter approach was taken in[16], where it is shown that an imperfect but close emulation exists between the classical cell cycle switch[32] and AM, which is already sufficient to establish the nearoptimal performance of the cell cycle switch. Moreover, a perfect emulation exists between the GW cell cycle switch[29] and AM, suggesting that a GWlike network would perform better. In general, ideal emulation relationships may suggest similar, possibly less than perfect, connections between networks. Thus, network morphisms and emulations provide a new perspective on network structure and similarity that may give helpful insights even when not perfectly realized.
Compositionality and modularity of network emulation
The solid arrows in Figure 3 are transitive in the sense that all the relevant morphism properties formally compose (see Additional file5). For example, QI has in general 12 distinct trajectories (4 nodes with 3 species each), and they can be made to emulate any 9 trajectories of CCR, and any 6 trajectories of SI and MI. The 6 trajectories of MI, also, can emulate any 3 trajectories of AM. It then follows by compositionality that the 12 trajectories of QI can emulate any 3 trajectories of AM.
Summary of formal results
We now give an overview of key definitions and theorems that analytically characterize the notions of network morphism and emulation: the detailed presentation is found in Methods.
 1.A reactant morphism is such that the reactants of each reaction are mapped by m _{ R } according to the mapping m _{ S } on species, but the reaction rates and products are unconstrained. Since network structure is related to stoichiometry, that condition turns out to be equivalent to a condition over reactant matrices:${{\mathit{m}}_{\mathit{S}}}^{\mathbf{T}}\xb7\mathit{\rho}=\widehat{\mathit{\rho}}\xb7{{\mathit{m}}_{\mathit{R}}}^{\mathbf{T}}\phantom{\rule{.3em}{0ex}}\text{(Definition: Reactant morphism)}$
 2.A stoichiomorphism relates the stoichiometry of the networks in a deeper way. The instantaneous stoichiometry of a species s in a reaction ρ → ^{ k } π is defined as φ(s, ρ → ^{ k } π) = k · (π _{ s }  ρ _{ s }), that is the net stoichiometry in the reaction multiplied by the reaction rate. A stoichiomorphism is then defined to satisfy:$\mathit{\phi}\xb7{\mathit{m}}_{\mathit{R}}={\mathit{m}}_{\mathit{S}}\xb7\widehat{\mathit{\phi}}\phantom{\rule{.3em}{0ex}}\text{(Definition: Stoichiomorphism)}$
 3.An emulation is a morphism that relates the differential systems of two networks. The differential system $F\in {\mathit{\mathbb{R}}}_{+}^{S}\to {\mathit{\mathbb{R}}}^{S}$ of (S, R) gives for each state $\mathit{v}\in {\mathit{\mathbb{R}}}_{+}^{S}$ (where ν assigns concentrations to species), and for each species s ∈ S, the derivative F(v)(s) ∈ ℝ of the concentration of s in v. F is defined according to the law of mass action, and similarly for $\widehat{F}$ over $\left(\widehat{S},\phantom{\rule{0.12em}{0ex}}\widehat{R}\right)$. A morphism is an emulation if:$\forall \widehat{\mathit{v}}\in {\mathit{\mathbb{R}}}_{+}^{\widehat{S}}\phantom{\rule{0.12em}{0ex}}F\left(\widehat{\mathit{v}}\circ {m}_{S}\right)=\widehat{F}\left(\widehat{\mathit{v}}\right)\circ {m}_{S}\phantom{\rule{.3em}{0ex}}\text{(Definition: Emulation)}$
This says that the derivatives of the two systems are related by m _{ S }. In particular, the derivatives coincide when the initial states coincide under m _{ S }, guaranteeing by determinism that whole trajectories coincide. This definition characterizes the coincidence of trajectories observed earlier in simulations.
Given any network morphism, the reactant morphism and stoichiomorphism conditions can be checked purely on the connectivity, stoichiometry, and rate constants of the networks. Our main theorem states that those structural conditions guarantee that the network morphism is a kinetic emulation:
Theorem (Emulation): If a morphisms (m _{ S }, m _{ R }) is a reactant morphism and a stoichiomorphism, then it is an emulation.
That is, for any choice of initial conditions$\widehat{\mathit{v}}$ of$\left(\widehat{S},\phantom{\rule{0.12em}{0ex}}\widehat{R}\right)$ we can pick initial conditions$\widehat{\mathit{v}}\circ {m}_{S}$ for (S, R) such that the trajectories of the two systems coincide. The rates of the two networks are coupled by the stoichiomorphism condition, but a second theorem then guarantees that we can achieve emulation also after changing rates. We show that if there is a stoichiomorphism from (S, R) to$\left(\widehat{S},\phantom{\rule{0.12em}{0ex}}\widehat{R}\right)$ and we change the rates of$\left(\widehat{S},\phantom{\rule{0.12em}{0ex}}\widehat{R}\right)$, then we can correspondingly change the rates of (S, R) so that we have again a stoichiomorphism, and hence the previous theorem applies.
Discussion: network ‘structure’ and rate independence
The ‘structure’ of a reaction network can be understood as its connectivity, without considerations of reaction rates: in this view the network structure is a graph with no kinetic information. Interesting notions of morphisms can be developed based on this level of information[10, 12, 13] as well useful tools[14, 15]. It is even possible to draw interesting conclusions about network kinetics just from the graph structure, but usually under some minimal assumptions about the underlying kinetics[1]. Reaction rates are usually present in the development of mathematical results, and are necessary to be able to state a rateindependence property.
More broadly, the structure of a network can be understood as the whole presentation of the network: the stateindependent information that does not change over time and that is, in particular, independent of the initial state. Rate constants, when provided, are structural in that sense: they are part of the ‘syntax’ or raw presentation of a network, just like stoichiometric constants, before any behavioral questions are considered. They are also technically part of the mathematical structure of chemical reaction networks. It may later be possible to show that rate information does not matter for certain network properties, but note that, similarly, some connectivity information can be shown not to matter (e.g., a reaction s → s).
As already mentioned, certain behavioral properties of individual networks, such as the number of steady states, may be characterized independently of reaction rates[1]. Here, however, we consider morphisms between two networks and their kinetic implications. It does not seem possible to assume that we can assign arbitrary rates to both networks and draw kinetic conclusions. But we can assign arbitrary rates to one network, and the rates of the other will follow in a systematic way (see the Change of Rates Theorem in Methods). This may be as rateindependent as we can hope for, when relating two networks.
Conclusions
The main question that we study is: what are the conditions under which a network can emulate another one? We have shown relatively complex examples of emulation, and outlined technical answers that rely only on network structure. We also explained how our main question is of interest in the study of chemical and biological networks. The formal results tell us that network emulation relates the mass action kinetics of two networks, and that it can be characterized by morphisms defined only over network connectivity, stoichiometry, and rate constants. There are several ways in which we can interpret these results more informally.
Network morphisms that are emulations provide an explanation of network structure, in that they reveal structural connections between networks that entail kinetic connections. For example, we may suspect that the main purpose of a networks is to stabilize a system in one of two states. An emulation from that network to the AM network can confirm that suspicion, as a dynamicalsystem analysis could also reveal. Moreover, the mapping of reactions that entails emulation explains how stabilization is achieved mechanistically, and because of known results about the speed of AM convergence to steady state, how fast it can happen.
Network emulation can be understood in terms of redundant implementation of a particular network kinetics. Redundancy, in our examples, is not just simple replication: even when QI is exactly reproducing the kinetics of MI, it is doing so through an intricately interconnected set of reactions. Redundant implementation may seem to be wasteful, but there are situations in which it arises naturally. Biological networks, for example, are not known for their minimality. Redundancy there may be due to material constraints that prevent the minimal realization of a network but allow a more complex one. Redundancy also implies robustness: for identical kinetics, the more complex realization will be less perturbed by link or node deletion than the minimal realization (a precise characterization requires a theory of approximation). In these circumstances, criteria that allow us to check whether a network is capable of implementing another one may lead to a better understanding of the essential and inessential aspects of the more complex networks. For example, we can now understand which aspects of the cell cycle switches GW and CCR lead them to emulate AM exactly or approximately[16].
If we have reason to believe that a complex network is implementing a simpler one, we may also use this knowledge for model reduction. For example, we may know that concentrations of certain species follow similar trajectories, and therefore they may be identified. This resembles the notion of abstraction or coarsegraining, which has been widely studied for discrete systems as well as continuous ones[35, 36]. Except that in that approach the abstract (simpler) network and the concrete (more complex) network should have, as much as possible, the same behavior for any parameter range. In our work, on the contrary, we seek an emulation or finegraining of a simple network yielding a more complex network that retains the original behavior in appropriate conditions[37], but that may well diverge from it in general. Emulation is less constraining than abstraction in that it has to work only in specific contexts, namely in the connectivity context of the simpler network.
Another aspect of emulation is its kinetic neutrality. Neutral drift in RNA landscapes[38] allows RNA systems to explore alternative organizations, including more complex ones, without at first affecting functionality. Similarly, we can look at a sequence of emulations, such as the ones in Figure 3, as tracing (backwards) a neutral path in network space. If an evolutionary event accidentally produces a new network that, perhaps approximately, emulates the previous one, then the new network will not be immediately selected against. We show that the conditions for emulation need not be demanding, operating with the same parameters as the old network and providing for rate compensation. Emulation steps compose, so they can be repeated many times, resulting in very distant and possibly much more complex networks that can later evolve new functionality due to their extra degrees of freedom. Another possibility is a lateral jump: we can go first neutrally from AM to AMr, and then introduce another species to obtain CCR (which does not emulate AMr) while still emulating AM. Hence, CCR is reached in two gradual steps (one of which is not an emulation) instead of a single complex jump. The diagram in Figure 3, which is certainly not exhaustive, thus emphasizes the richness of network space: kinetic neutrality must be relatively rare, yet so many meaningful connections can be found.
In conclusion, we have shown that kinetic connections between networks can be established via network morphism defined only on stateindependent attributes. In contrast, when studying directly the kinetic equations of a network, the structural properties are reduced to functional properties and disappear from sight, or at least become much more implicit. Therefore, network morphisms aid in understanding kinetic relationships between networks that are structural and not accidental: morphisms provide a structural reason for kinetic similarity.
Methods
We define morphism between chemical reaction networks (CRNs) based on static network structure: connectivity, stoichiometry, and rate constants. These are therefore morphisms of the syntax of the networks; we then study implications about the kinetics. We prove two main theorems: an Emulation Theorem stating structural conditions under which two networks can have identical traces for any choice of initial conditions, and a Change of Rates Theorem generalizing that result to any choice or reaction rates as well. In Additional file5 we discuss conditions under which the Emulation Theorem has an inverse.
Let ℕ be the natural numbers, ℤ be the integers, ℝ be the reals, ℝ _{+} be the nonnegative reals, and ℙ be the strictly positive reals. A set A has cardinality A. We write A → B and B ^{ A } for the functions from A to B. When f ∈ A → B and a ∈ A we use f(a) and f _{ a } for function application. A function f ∈ A → B has images f(X ⊆ A) = {b ∈ B ∃ a ∈ X f(a) = b} and fibers f ^{ 1}(b ∈ B) = {a ∈ Af(a) = b} (inverse images of singleton sets). A function f ∈ A × B → ℝ is a matrix of dimensions A × B, with matrix multiplication f · g and matrix transpose f ^{T}(i, j) = f(j, i).
Chemical reaction networks (CRNs)
where s _{ i } are the species,${\rho}_{{s}_{i}}$ are the multiplicities (stoichiometric numbers) of the reactant species,${\pi}_{{s}_{i}}$ are the multiplicities of the product species, and k is the reaction rate constant. We use 1^{ st }(ρ, π, k) = ρ.
(Two reactions ρ → ^{ k } π, ρ → ^{ k '} π are kinetically equivalent in mass action to a single reaction ρ → ^{ k + k '} π.)
As common[1, 39, 40], this definition of CRN considers only irreversible reactions, with reversible reactions modeled as pairs of opposite reaction. Moreover, it does not enforce conservation or detail balance laws so that it can cover the description of open chemical systems that, e.g., take energy or materials from the environment.
Species maps and complex maps
With these preliminaries, we can now define various morphisms between CRNs.
Morphisms between CRNs
A CRN morphism from (S, R) to$\left(\widehat{S},\phantom{\rule{0.12em}{0ex}}\widehat{R}\right)$ is a pair of maps${m}_{S}\in S\to \widehat{S}$ and${m}_{R}\in R\to \widehat{R}$. We write$m\in \left(S,\phantom{\rule{0.12em}{0ex}}R\right)\to \left(\widehat{S},\phantom{\rule{0.12em}{0ex}}\widehat{R}\right)$ for$m=\left({m}_{S},\phantom{\rule{0.12em}{0ex}}{m}_{R}\right)\in \left(S\to \widehat{S}\right)\times \left(R\to \widehat{R}\right)$, given that (S, R) and$\left(\widehat{S},\phantom{\rule{0.12em}{0ex}}\widehat{R}\right)$ are CRNs. See Figure 6 for some simple examples of morphisms between CRNs.
Linear algebra notation is sometimes useful for compactness. We write${\mathit{m}}_{\mathit{S}}\in S\times \widehat{S}\to \left\{0,\phantom{\rule{0.12em}{0ex}}1\right\}$ for the characteristic matrix of m _{ S }, such that$\forall s\in S\forall \widehat{s}\in \widehat{S}\phantom{\rule{0.5em}{0ex}}{\mathit{m}}_{\mathit{S}}\left(s,\phantom{\rule{0.12em}{0ex}}\widehat{s}\right)=\left(\mathit{if}\phantom{\rule{0.12em}{0ex}}{m}_{\mathsf{S}}\left(s\right)=\widehat{s}\phantom{\rule{0.24em}{0ex}}\mathit{then}\phantom{\rule{0.24em}{0ex}}1\phantom{\rule{0.24em}{0ex}}\mathit{else}\phantom{\rule{0.24em}{0ex}}0\right)$, and similarly for${\mathit{m}}_{R}\in R\times \widehat{R}\to \left\{0,\phantom{\rule{0.12em}{0ex}}1\right\}$. We write φ _{(S,R)} ∈ S × R → ℝ for the (instantaneous) stoichiometric matrix such that ∀ s ∈ S ∀ r ∈ R φ _{(S,R)}(s, r) = φ(s, r). We write ρ _{(S, R)} ∈ S × R → ℕ for the reactant matrix such that ∀ s ∈ S ∀ ρ → ^{ k } π ∈ R ρ _{(S, R)}(s, ρ → ^{ k } π) = ρ _{ s }.
(This property does not imply homomorphism, e.g., m(s → s) = 2s → 2s, or m(s _{0} → s _{0} + s _{1}) = 2s _{0} → 2s _{0} + s _{1}, determine maps between singlereaction CRNs that are not homomorphisms.)
Although homomorphisms produce the most natural examples of network morphisms, it is useful for our results to consider weaker morphisms that preserve just the reactant component of a reaction (which is helpful because the mass action of a reaction is defined on the reactants).
In the sequel we often omit the subscripts on m _{ S } and m _{ R }.
We can say that a stoichiomorphism preserves stoichiometry over reaction fibers, meaning that the sum of the stoichiometry of any s ∈ S in the reactions that map to any$\widehat{r}\in \widehat{R}$, must equal the stoichiometry of m(s) in$\widehat{r}$.
The combination of homomorphism and stoichiomorphism forms a natural notion of network correspondence, strongly preserving both the graph structure and the stoichiometric structure of a network, while not requiring in general injectivity or surjectivity of the mapping. Many typical examples fall into this combined class. However, weaker notions of network morphisms are also useful, so we investigate stoichiomorphisms separately from homomorphisms, and in particular we combine them with reactant morphisms.
Remark – homomorphic projection
Given a CRN (S, R) and a species map${m}_{S}\in S\to \mathcal{S}$, the homomorphic projection of (S, R) via m _{ S } is the CRN (m _{ S }(S), m _{ R }(R)) where m _{ R }(ρ → ^{ k } π) = m _{ S }(ρ) → ^{ k } m _{ S }(π). By construction (m _{ S }(S), m _{ R }(R)) is a CRN, because R ⊆ ℛ_{ S } and ρ → ^{ k } π ∈ R imply ρ ∈ ℕ ^{ S } and${m}_{S}\left(\rho \right)\in {\mathbb{N}}^{{m}_{S}\left(S\right)}$, and similarly for π, hence${m}_{\mathit{\mathcal{R}}}\left(\rho {\to}^{k}\pi \right)\in {{\mathit{\mathcal{R}}}_{{m}_{S}}}_{\left(S\right)}$ and${m}_{R}\left(R\right)\subseteq {\mathit{\mathcal{R}}}_{{m}_{S}\left(S\right)}$. Moreover, the homomorphic projection is a CRN epimorphism: m _{ S } and m _{ R } are surjective. So if$\widehat{\rho}{\to}^{\widehat{k}}\widehat{\pi},\phantom{\rule{0.12em}{0ex}}\widehat{\rho}{\to}^{\widehat{k}\text{'}}\widehat{\pi}\in {m}_{R}\left(R\right)$ then there is some ρ → ^{ k } π, ρ → ^{ k '} π ∈ R with$\widehat{\rho}={m}_{S}\left(\rho \right)$,$\widehat{\pi}={m}_{S}\left(\pi \right)$,$\widehat{k}=k$,$\widehat{k}\text{'}=k\text{'}$. Since (S, R) is a CRN we have k = k’ and hence$\widehat{k}=\widehat{k}\text{'}$.
Remark – mathematical and graphical representations of CRNs
Following[40], a CRN is defined above as a finite set of species S ⊆ S and a reaction relation R ⊆ ℛ_{ S } = ℕ ^{ S } × ℕ ^{ S } × ℙ over complexes ℕ ^{ S }. In the CRNT approach[1], a CRN is instead defined as finite set of complexes C ⊆ ℕ ^{ S } and a reaction relation R ⊆ ℛ_{ C } = C × C × ℙ. We can of course translate between these two representations: our reaction triples ρ → ^{ k } π = (ρ, π, k) are already in the CRNT format. We even obtain the same graphical representation by drawing the directed graphs of the relations ℛ_{ S } and ℛ_{ C } with complexes as nodes and reactions as edges (labeling the edges with the rate ℙ), which is the most common depiction of CRNs.
The morphisms that arise from those graphs, however, are different from ours. In the CRNT case we would naturally map complexes to complexes, but it is hard to see what kinetic properties could be preserved by such mappings, unless we required further relationships between the complexes. Our morphisms instead map species to species, which constrains the relationships between the complexes being mapped. The notion of morphism that we study is therefore predicated on the (S, R) representation.
To better visualize these CRNs and their morphisms, we draw graphs where nodes represent species, not complexes. A reaction then becomes a ‘directed multiedge’ with sets of species as source and target. A good way to visualize such edges is as Petri nets: a CRN is drawn as a directed bipartite graph between (round) species nodes and (square) reaction nodes (many bipartite variants exist[9], but we conform to Petri nets[34]). A species that occurs as a reactant in a reaction is connected to it by a directed edge, with as many such edges as the stoichiometric number of the species. A reaction that produces a species as a product is connected to it by a directed edge, with as many such edges as the stoichiometric number of the species. If we omit an arrowhead on an edge, the edge is by convention directed from the species node to the reaction node. Reaction rates are affixed to the reaction nodes; if omitted they are equal to 1. A CRN morphism is then represented as a mapping between species nodes and a mapping between reaction nodes, between two Petri nets (Figure 6).
CRN kinetics
We now give a formulation of standard mass action kinetics, for the purpose of next presenting results about the connection between CRN morphisms and kinetics.
A species s has dimension mol (amount of substance); its concentration has dimension molarity M = mol · l^{‒ 1}. A reaction r = ρ → ^{ k } π has order r = ∑ _{ s ∈ S } ρ _{ s } and its rate k has dimension s^{‒ 1} · M^{1 ‒ r}.
For each initial state v the differential system F has a unique maximal (for some t ≤ + ∞) differentiable solution $f\in \left[0,\phantom{\rule{0.12em}{0ex}}\left(\right)close=")">t,\to ,{\mathit{\mathbb{R}}}_{+}^{S}\right.$ such that f(0) = v and$\dot{f}=F\circ f$, where$\dot{f}$ is the (time) derivative of f [CauchyLipschitz]. Since F arises from a CRN, the solution is everywhere nonnegative[1].
Note that it is traditional to group the three factors of the differential system as (π _{ s }  ρ _{ s }) · (k · v ^{ ρ }), where k · v ^{ ρ } are the rate functions, which are interpreted as the kinetics of the system: kinetics other than mass action can be investigated by varying the rate functions. In mass action, however, the grouping (k · (π _{ s }  ρ _{ s })) · v ^{ ρ } is also natural because k · (π _{ s }  ρ _{ s }) = φ(s, r) is the stateindependent factor and v ^{ ρ } = [r]_{ v } is the statedependent factor, so that the syntactic structure of the network is gathered in φ(s, r). By representing reactions as triples (ρ, π, k) we have already committed to a rate function that depends on a single rate parameter k. Other common kinetics can be approximated in mass action, such as MichaelisMenten enzyme kinetics, and Hill kinetics as in the triplet motif of Figure 2.
CRN emulation
Starting from the bottom left towards the top right, if we apply the emulating F to a state$\widehat{\mathit{v}}\circ m$ of S which is a copy under m ^{1} of a state$\widehat{\mathit{v}}$ of$\widehat{S}$, we should obtain derivatives that are copies under m ^{1} of the derivatives of the emulated$\widehat{F}$ in state$\widehat{\mathit{v}}$.
Hence both$f,\phantom{\rule{0.24em}{0ex}}\widehat{f}$ and$\dot{f},\phantom{\rule{0.12em}{0ex}}\dot{\widehat{f}}$ coincide at 0, and therefore every strajectory of F will coincide with the m(s) trajectory of$\widehat{F}$ until the maximal time of$\widehat{f}$. (If m is not surjective, then$\widehat{f}$ as a whole could stop sooner than f, e.g., if$\widehat{S}$ has an extra species s′ and$\widehat{R}$ has an extra reaction 2s′ → 3s′).
In summary, an emulation$m\in \left(S,\phantom{\rule{0.12em}{0ex}}R\right)\to \left(\widehat{S},\phantom{\rule{0.12em}{0ex}}\widehat{R}\right)$ is such that for any initial condition$\widehat{\mathit{v}}$ for$\left(\widehat{S},\phantom{\rule{0.12em}{0ex}}\widehat{R}\right)$ there is some initial condition for (S, R), namely$\widehat{\mathit{v}}\circ m$, such that (S, R) can exactly reproduce the kinetics of$\left(\widehat{S},\phantom{\rule{0.12em}{0ex}}\widehat{R}\right)$.
Emulation theorem
We can now relate CRN stoichiomorphisms to CRN emulations.
Lemma  mass action
 1)
Let $m\in \left(S,\phantom{\rule{0.12em}{0ex}}R\right)\to \left(\widehat{S},\phantom{\rule{0.12em}{0ex}}\widehat{R}\right)$ be any CRN morphism. For any state $\widehat{\mathit{v}}\in {\mathit{\mathcal{R}}}_{+}^{\widehat{S}}$ and complex ρ ∈ ℕ ^{ S }, we have ${\left(\widehat{\mathit{v}}\circ m\right)}^{\rho}={\widehat{\mathit{v}}}^{m\left(\rho \right)}$.
 2)
Let $m\in \left(S,\phantom{\rule{0.12em}{0ex}}R\right)\to \left(\widehat{S},\phantom{\rule{0.12em}{0ex}}\widehat{R}\right)$ be a CRN reactant morphism. For any reaction r ∈ R and state $\widehat{\mathit{v}}\in {\mathit{\mathbb{R}}}_{+}^{\widehat{S}}$, we have ${\left[r\right]}_{\widehat{\mathit{v}}\circ m}={\left[m\left(r\right)\right]}_{\widehat{\mathit{v}}}$.
 1)By the Fiber Lemma (see Additional file 5), for any function $m\in S\to \widehat{S}$ with $g\in S\times \widehat{S}\to {\mathit{\mathbb{R}}}_{\mathit{+}}$ we have $\prod}_{\widehat{s}\in \widehat{S}}{\displaystyle {\prod}_{s\in {m}^{1}\left(\widehat{s}\right)}g\left(s,\phantom{\rule{0.12em}{0ex}}\widehat{s}\right)}={\displaystyle {\prod}_{s\in S}g\left(s,m\left(s\right)\phantom{\rule{0.12em}{0ex}}\right)$. Then:$\begin{array}{ll}{\widehat{\mathit{v}}}^{m\left(\rho \right)}={\prod}_{\widehat{s}\in \widehat{S}}{\widehat{\mathit{v}}}_{\widehat{s}}^{m{\left(\rho \right)}_{\widehat{s}}}& \mathrm{notational}\phantom{\rule{0.25em}{0ex}}\mathrm{definition}\\ ={\prod}_{\widehat{s}\in \widehat{S}}{\widehat{\mathit{v}}}_{\widehat{S}}^{{\sum}_{s\in {m}^{1}\left(\widehat{s}\right)}{\rho}_{s}}& \mathrm{by}\phantom{\rule{0.25em}{0ex}}\mathrm{definition}\phantom{\rule{0.25em}{0ex}}\mathrm{of}\phantom{\rule{0.5em}{0ex}}\mathit{m}\left(\rho \right)\\ ={\prod}_{\widehat{s}\in \widehat{S}}{\prod}_{s\in {m}^{1}\left(\widehat{s}\right)}{\widehat{\mathit{v}}}_{\widehat{s}}^{\phantom{\rule{0.05em}{0ex}}{\rho}_{s}}& \mathrm{by}\phantom{\rule{0.25em}{0ex}}\mathrm{distribution}\phantom{\rule{0.25em}{0ex}}\mathrm{of}\phantom{\rule{0.25em}{0ex}}\mathrm{exponents}\\ ={\prod}_{s\in S}{{\widehat{\mathit{v}}}_{m\left(s\right)}}^{{\rho}_{s}}& \mathrm{by}\phantom{\rule{0.5em}{0ex}}\mathrm{the}\phantom{\rule{0.5em}{0ex}}\mathrm{Fiber}\phantom{\rule{0.5em}{0ex}}\mathrm{Lemma}\phantom{\rule{0.5em}{0ex}}\mathrm{with}\phantom{\rule{0.5em}{0ex}}g\left(x,\phantom{\rule{0.5em}{0ex}}y\right)\phantom{\rule{0.1em}{0ex}}=\phantom{\rule{0.2em}{0ex}}{{\widehat{\mathit{v}}}_{y}}^{{\rho}_{x}}\\ ={{\prod}_{s\in S}\left(\widehat{\mathit{v}}\circ m\right)}_{s}^{{\rho}_{s}}& \mathrm{by}\phantom{\rule{0.25em}{0ex}}\mathrm{definition}\phantom{\rule{0.25em}{0ex}}\mathrm{of}\phantom{\rule{0.5em}{0ex}}\circ \\ ={\left(\widehat{\mathit{v}}\circ m\right)}^{\rho}& \mathrm{notational}\phantom{\rule{0.25em}{0ex}}\mathrm{definition}\end{array}$
 2)Let r = ρ → ^{ k } π. Then:$\begin{array}{ll}{\left[m\left(r\right)\right]}_{\widehat{\mathit{v}}}={\left[m\left(\rho \right){\to}^{\widehat{k}}\widehat{\pi}\right]}_{\widehat{\mathit{v}}}& \phantom{\rule{1.4em}{0ex}}\mathrm{by}\phantom{\rule{0.25em}{0ex}}\mathrm{definition}\phantom{\rule{0.25em}{0ex}}\mathrm{of}\phantom{\rule{0.25em}{0ex}}\mathrm{reactant}\phantom{\rule{0.25em}{0ex}}\mathrm{morphism}\phantom{\rule{0.25em}{0ex}}\mathrm{for}\phantom{\rule{0.25em}{0ex}}\mathrm{some}\phantom{\rule{0.5em}{0ex}}\widehat{\pi},\phantom{\rule{0.12em}{0ex}}\widehat{k}\\ ={\widehat{\mathit{v}}}^{m\left(\rho \right)}& \phantom{\rule{1.4em}{0ex}}\mathrm{by}\phantom{\rule{0.25em}{0ex}}\mathrm{definition}\phantom{\rule{0.25em}{0ex}}\mathrm{of}\phantom{\rule{0.25em}{0ex}}\mathrm{mass}\phantom{\rule{0.25em}{0ex}}\mathrm{action}\phantom{\rule{0.25em}{0ex}}\mathrm{of}\phantom{\rule{0.25em}{0ex}}\mathrm{a}\phantom{\rule{0.25em}{0ex}}\mathrm{reaction}\\ ={\left(\widehat{\mathit{v}}\circ m\right)}^{\rho}& \phantom{\rule{1.4em}{0ex}}\mathrm{by}\phantom{\rule{0.25em}{0ex}}\left(1\right)\\ ={\left[\rho {\to}^{k}\pi \right]}_{\widehat{\mathit{v}}\circ m}={\left[r\right]}_{\widehat{\mathit{v}}\circ m}& \phantom{\rule{1.4em}{0ex}}\mathrm{by}\phantom{\rule{0.25em}{0ex}}\mathrm{definition}\phantom{\rule{0.25em}{0ex}}\mathrm{of}\phantom{\rule{0.25em}{0ex}}\mathrm{mass}\phantom{\rule{0.25em}{0ex}}\mathrm{action}\phantom{\rule{0.25em}{0ex}}\mathrm{of}\phantom{\rule{0.25em}{0ex}}\mathrm{a}\phantom{\rule{0.25em}{0ex}}\mathrm{reaction}\end{array}$
■
Theorem  emulation
Proof
■
Note that a stoichiomorphism need not be a homomorphism for the theorem to hold: both rates and products may be allowed to vary as long as the stoichiomorphism property is satisfied. However the stoichiomorphism needs to be a reactant morphism because that is the basis of the Mass Action lemma. Because of the form of the differential system F(v)(s) = ∑ _{ r ∈ R } φ(s, r) · [r]_{ v }, the stoichiomorphism condition supports the φ(s, r) factor, while the reactant morphism condition supports the [r]_{ v } factor. As matrix equations, the assumptions for this theorem are${{\mathit{m}}_{\mathit{S}}}^{\mathbf{T}}\xb7\mathit{\rho}\mathbf{=}\widehat{\mathit{\rho}}\xb7{{\mathit{m}}_{\mathit{R}}}^{\mathbf{T}}$ (reactant morphism) and$\mathbf{\phi}\xb7{\mathit{m}}_{\mathit{R}}\mathbf{=}{\mathit{m}}_{\mathit{S}}\xb7\widehat{\mathbf{\phi}}$ (stoichiomorphism). (See the supporting examples 10, 11, 12, in Additional file2.)
To conclude this section, the following proposition shows that any emulation, and hence any morphism that is both a reactant morphism and a stoichiomorphism, preserves steady states from the target CRN to the source CRN. Therefore the existence of a stoichiomorphism/reactant morphism can be used to determine, syntactically, that a CRN has at least certain steady states inherited from another CRN whose steady states are known.
Proposition – steady states
Proof
We have that$\forall \widehat{\mathit{v}}\in {\mathit{\mathbb{R}}}_{+}^{\widehat{S}}\forall s\in S\phantom{\rule{0.24em}{0ex}}\widehat{F}\left(\widehat{\mathit{v}}\right)\left(m\left(s\right)\right)=F\left(\widehat{\mathit{v}}\circ m\right)\left(s\right)$.
Hence, if$\widehat{\mathit{v}}\in {\mathit{\mathbb{R}}}_{+}^{\widehat{S}}$ is a (partial) steady state of$\left(\widehat{S},\phantom{\rule{0.12em}{0ex}}\widehat{R}\right)$ via m, that is if$\forall s\in S\phantom{\rule{0.24em}{0ex}}\widehat{F}\left(\widehat{\mathit{v}}\right)\left(m\left(s\right)\right)=0$, then$\forall s\in S\phantom{\rule{0.24em}{0ex}}F\left(\widehat{\mathit{v}}\circ m\right)\left(s\right)=0$, so$\widehat{\mathit{v}}\circ m\in {\mathit{\mathbb{R}}}_{+}^{S}$ is a (full) steady state of (S, R). The same holds if$\widehat{\mathit{v}}$ is a full steady state.
Conversely, if m is a bijection on species and$\mathit{v}\in {\mathit{\mathbb{R}}}_{+}^{S}$ is a steady state of F, then take$\widehat{\mathit{v}}\mathit{=}\mathit{v}\circ {m}^{\u20121}$. We have that for all s ∈ S,$\widehat{F}\left(\widehat{\mathit{v}}\right)\left(m\left(s\right)\right)=F\left(\widehat{\mathit{v}}\circ m\right)\left(s\right)=F\left(\mathit{v}\right)\left(s\right)=0$; hence$\widehat{\mathit{v}}$ is a steady state of$\widehat{F}$.
■
The converse direction fails if m is not injective: if v is any steady state of (S, R), then there might not be any$\widehat{\mathit{v}}$ with$\widehat{\mathit{v}}\circ m=\mathit{v}$ (if ‘all v _{ S } differ’) and hence no corresponding steady state in$\left(\widehat{S},\phantom{\rule{0.12em}{0ex}}\widehat{R}\right)$.
Change of rates theorem
We now generalize the Emulation Theorem, which allows us to choose arbitrary initial conditions, to also allow choosing arbitrary rates for the target CRN. (See the supporting Example 10 in Additional file2.)
A change of rates is a CRN morphism ι ∈ (S, R) → (S, R′) such that ι _{ S } ∈ S → S is the identity and ι _{ R } ∈ R → R′ is a bijection such that:${\iota}_{R}\left(\rho {\to}^{k}\pi \right)={\rho}^{\prime}{\to}^{{k}^{\prime}}{\pi}^{\prime}\Rightarrow \rho ={\rho}^{\prime}\&\pi ={\pi}^{\prime}$. Then ι ^{ 1} ∈ (S, R′) → (S, R) is also a change of rates. As usual, we omit the subscripts on ι _{ S }, ι _{ R }.
Theorem  change of rates
Let$m\in \left(S,R\right)\to \left(\widehat{S},\widehat{R}\right)$ be a stoichiomorphism, and$\widehat{\iota}\in \left(\widehat{S},\widehat{R}\right)\to \left(\widehat{S},{\widehat{R}}^{\prime}\right)$ be any change of rates. Then there is a change of rates$\iota \in \left(S,R\right)\to \left(S,{R}^{\prime}\right)$ such that$\widehat{\iota}\circ m\circ {\iota}^{1}\in \left(S,{R}^{\prime}\right)\to \left(\widehat{S},{\widehat{R}}^{\prime}\right)$ is a stoichiomorphism.
Proof
Note that without the condition on CRNs that$\rho {\to}^{{k}_{1}}\pi ,\phantom{\rule{0.12em}{0ex}}\rho {\to}^{{k}_{2}}\pi \in R\Rightarrow \phantom{\rule{0.48em}{0ex}}{k}_{1}={k}_{2}$, we might have that ι is not injective, if we had two such reactions in R with k _{1} ≠ k _{2}, and for some values${k}_{1}\xb7\frac{{\widehat{k}}_{1}^{\prime}}{{\widehat{k}}_{1}}={k}_{2}\xb7\frac{{\widehat{k}}_{2}^{\prime}}{{\widehat{k}}_{2}}$. This is the one place where we take advantage of that condition, so that ι is injective and hence is a change of rates.
Take any s ∈ S and any${\widehat{r}}^{\prime}={\widehat{\rho}}^{\prime}{\to}^{{\widehat{k}}^{\prime}}{\widehat{\pi}}^{\prime}\in {\widehat{R}}^{\prime}$.
Let$\widehat{r}=\widehat{\rho}{\to}^{\widehat{k}}\widehat{\pi}={\widehat{\iota}}^{1}\left({\widehat{r}}^{\prime}\right)$, implying$\widehat{\rho}={\widehat{\rho}}^{\prime}\phantom{\rule{0.36em}{0ex}}\&\phantom{\rule{0.36em}{0ex}}\widehat{\pi}={\widehat{\pi}}^{\prime}$ because$\widehat{\iota}$ is a change of rates.
Which means that$\widehat{\iota}\circ m\circ {\iota}^{1}$ is a stoichiomorphism.
■
Corollary  emulation under reactant morphism
If$m\in \left(S,R\right)\to \left(\widehat{S},\widehat{R}\right)$ is a CRN reactant morphism and a CRN stoichiomorphism, then for any change of rates$\widehat{\iota}\in \left(\widehat{S},\widehat{R}\right)\to \left(\widehat{S},{\widehat{R}}^{\prime}\right)$, there is some change of rates ι ∈ (S, R) → (S, R′) such that$\widehat{\iota}\circ m\circ {\iota}^{1}\in \left(S,{R}^{\prime}\right)\to \left(\widehat{S},\widehat{{R}^{\prime}}\right)$ is a CRN reactant morphism and a CRN emulation.
Proof
By the Change of Rates Theorem for any$\widehat{\iota}$ there is a ι such that$\widehat{\iota}\circ m\circ {\iota}^{1}$ is a stoichiomorphism. Moreover, if m is a reactant morphism then$\widehat{\iota}\circ m\circ {\iota}^{1}$ is too, because$\widehat{\iota}$, ι are changes of rates and so:$\left(\widehat{\iota}\circ m\circ {\iota}^{1}\right)\left(\rho {\to}^{{k}^{\prime}}\pi \right)=\left(\widehat{\iota}\circ m\right)\left(\rho {\to}^{k}\pi \right)=\widehat{\iota}\left(m\left(\rho \right){\to}^{\widehat{k}}\widehat{\pi}\right)=m\left(\rho \right){\to}^{\widehat{{k}^{\prime}}}\widehat{\pi}$, where$m\left(\rho \right)=\left(i{d}_{\widehat{S}}\circ m\circ i{d}_{S}\right)\left(\rho \right)=\left(i\circ m\circ {\iota}^{1}\right)\left(\rho \right)$ because$\widehat{\iota}$, ι are identities on species. Therefore, by the Emulation Theorem,$\widehat{\iota}\circ m\circ {\iota}^{1}$ is an emulation.
■
Corollary  emulation under homomorphism
If$m\in \left(S,R\right)\to \left(\widehat{S},\widehat{R}\right)$ is a CRN homomorphism and a CRN stoichiomorphism, then for any change of rates$\widehat{\iota}\in \left(\widehat{S},\widehat{R}\right)\to \left(\widehat{S},{\widehat{R}}^{\prime}\right)$, there is some change of rates ι ∈ (S, R) → (S, R′) such that$\widehat{\iota}\circ m\circ {\iota}^{1}\in \left(S,{R}^{\prime}\right)\to \left(\widehat{S},{\widehat{R}}^{\prime}\right)$ is a CRN homomorphism and a CRN emulation.
Proof
By the Change of Rates Theorem for any$\widehat{\iota}$ there is a ι such that$\widehat{\iota}\circ m\circ {\iota}^{1}$ is a stoichiomorphism (and such that${k}^{\prime}=k\xb7\frac{\widehat{{k}^{\prime}}}{\widehat{k}}$). Moreover, if m is a homomorphism and$\widehat{\iota}$, ι are changes of rates, we have:$\left(\widehat{\iota}\circ m\circ {\iota}^{1}\right)\left(\rho {\to}^{{k}^{\prime}}\pi \right)=\left(\widehat{\iota}\circ m\right)\left(\rho {\to}^{k}\pi \right)=\widehat{\iota}\left(m\left(\rho \right){\to}^{\widehat{k}}m\left(\pi \right)\right)=m\left(\rho \right){\to}^{\widehat{{k}^{\prime}}}m\left(\pi \right)$, where$m\left(\rho \right)=\left(i{d}_{\widehat{S}}\circ m\circ i{d}_{S}\right)\left(\rho \right)=\left(i\circ m\circ {\iota}^{1}\right)\left(\rho \right)$ and m(π) = (i ∘ m ∘ ι ^{ 1})(π) because$\widehat{\iota}$, ι are identities on species. Moreover${k}^{\prime}=k\xb7\frac{\widehat{{k}^{\prime}}}{\widehat{k}}$ where$k=\widehat{k}$ by homomorphism and hence${k}^{\prime}=\widehat{{k}^{\prime}}$. Therefore$\widehat{\iota}\circ m\circ {\iota}^{1}$ is a homomorphism, and hence a reactant morphism. Therefore, by the Emulation Theorem,$\widehat{\iota}\circ m\circ {\iota}^{1}$ is an emulation.
■
In fact, as exemplified in Additional file3, all the emulation morphism in Figure 3 (which are all homomorphism) can be checked with this strategy, by checking the net stoichiomorphism condition between unit rate networks.
In summary, a stoichiomorphism$m\in \left(S,\phantom{\rule{0.12em}{0ex}}R\right)\to \left(\widehat{S},\phantom{\rule{0.12em}{0ex}}\widehat{R}\right)$ that is also a reactant morphism, determines an emulation for any choice of rates of$\left(\widehat{S},\phantom{\rule{0.12em}{0ex}}\widehat{R}\right)$. Those emulations can match any initial conditions of any choice of rates of$\left(\widehat{S},\phantom{\rule{0.12em}{0ex}}\widehat{R}\right)$ with some initial conditions of some choice of rates of (S, R).
Abbreviations
 AM:

Approximate majority (algorithm)
 CRN:

Chemical reaction network.
Declarations
Acknowledgements
David Soloveichik made key observations that directed this work onto its analytical path. Attila CsikászNagy provided deep insights on cell cycle switches and other biological networks.
Authors’ Affiliations
References
 Craciun G, Tang Y, Feinberg M: Understanding bistability in complex enzymedriven reaction networks. Proc Natl Acad Sci. 2006, 103 (23): 86978702. 10.1073/pnas.0602767103.PubMed CentralView ArticlePubMedGoogle Scholar
 Joshi B, Shiu A: Atoms of multistationarity in chemical reaction network. J Math Chem. 2013, 51 (1): 153178.Google Scholar
 Feliu E, Wiuf C: Simplifying biochemical models with intermediate species. J R Soc Interface. 2013, 10 (87): 2013048410.1098/rsif.2013.0484.PubMed CentralView ArticlePubMedGoogle Scholar
 Okino MS, Mavrovouniotis ML: Simplification of mathematical models of chemical reaction systems. Chem Rev. 1998, 98 (2): 391408. 10.1021/cr950223l.View ArticlePubMedGoogle Scholar
 Orth JD, Thiele I, Palsson BØ: What is flux balance analysis?. Nat Biotechnol. 2010, 28: 245248. 10.1038/nbt.1614.PubMed CentralView ArticlePubMedGoogle Scholar
 Alon U: Network motifs: theory and experimental approaches. Nat Rev Genet. 2007, 8: 450461. 10.1038/nrg2102.View ArticlePubMedGoogle Scholar
 Tyson JJ: Classification of instabilities in chemical reaction systems. J Chem Phys. 1975, 62: 10101015. 10.1063/1.430567.View ArticleGoogle Scholar
 Soulé C: Graphic requirements for multistationarity. Com Plex Us. 2003, 1: 123133.Google Scholar
 Mincheva M, Roussel MR: Graphtheoretic methods for the analysis of chemical and biochemical networks. I. Multistability and oscillations in ordinary differential equation models. J Math Biol. 2007, 55: 6186. 10.1007/s0028500700991.View ArticlePubMedGoogle Scholar
 Naldi A, Remy E, Thieffry D, Chaouiya C: Dynamically consistent reduction of logical regulatory graphs. Theor Comput Sci. 2011, 412 (21): 22072218. 10.1016/j.tcs.2010.10.021.View ArticleGoogle Scholar
 Soliman S: A stronger necessary condition for the multistationarity of chemical reaction networks. Bull Math Biol. 2013, 75 (11): 22892303. 10.1007/s1153801398937.View ArticlePubMedGoogle Scholar
 Gay S, Soliman S, Fages F: A graphical method for reducing and relating models in systems biology. Bioinformatics. 2010, 26 (18): i575i581. 10.1093/bioinformatics/btq388.PubMed CentralView ArticlePubMedGoogle Scholar
 Gay S, Fages F, Martinez T, Soliman S, Solnon C: On the subgraph epimorphism problem. Discret Appl Math. 2014, 162: 214228.View ArticleGoogle Scholar
 Tian Y, McEachin RC, Santos C, States DJ, Patel JM: SAGA: a subgraph matching tool for biological graphs. Bioinformatics. 2007, 23 (2): 232239. 10.1093/bioinformatics/btl571.View ArticlePubMedGoogle Scholar
 Ferro A, Giugno R, Pigola G, Pulvirenti A, Skripin D, Bader GD, Shasha D: NetMatch: a Cytoscape plugin for searching biological networks. Bioinformatics. 2007, 23 (7): 910912. 10.1093/bioinformatics/btm032.View ArticlePubMedGoogle Scholar
 Cardelli L, CsikászNagy A: The cell cycle switch computes approximate majority. Sci Rep. 2012, 2: 656PubMed CentralView ArticlePubMedGoogle Scholar
 Fages F, Soliman S: From reaction models to influence graphs and back: a theorem. Formal Methods Syst Biol. 2008, LNCS 5054: 90102.View ArticleGoogle Scholar
 Verdugo A, Vinod PK, Tyson JJ, Novak B: Molecular mechanisms creating bistable switches at cell cycle transitions. Open Biol. 2013, 3 (3): 12017910.1098/rsob.120179.PubMed CentralView ArticlePubMedGoogle Scholar
 Tyson JJ, Novák B: Functional motifs in biochemical reaction networks. Annu Rev Phys Chem. 2010, 61: 219240. 10.1146/annurev.physchem.012809.103457.PubMed CentralView ArticlePubMedGoogle Scholar
 Motegi F, Seydoux G: The PAR network: redundancy and robustness in a symmetrybreaking system. Philos Trans R Soc Lond B Biol Sci. 2013, 368 (1629): 2013001010.1098/rstb.2013.0010.PubMed CentralView ArticlePubMedGoogle Scholar
 Bajpai A, Feoktistova A, Chen JS, McCollum D, Sato M, CarazoSalas RE, Gould KL, CsikászNagy A: Dynamics of SIN asymmetry establishment. PLoS Comput Biol. 2013, 9 (7): e100314710.1371/journal.pcbi.1003147.PubMed CentralView ArticlePubMedGoogle Scholar
 Oates AC, Morelli LG, Ares S: Patterning embryos with oscillations: structure, function and dynamics of the vertebrate segmentation clock. Development. 2012, 139: 625639. 10.1242/dev.063735.View ArticlePubMedGoogle Scholar
 Gardner TS, Cantor CR, Collins JJ: Construction of a genetic toggle switch in Escherichia coli. Nature. 2000, 403 (6767): 339342. 10.1038/35002131.View ArticlePubMedGoogle Scholar
 Yang X, Lau KY, Sevim V, Tang C: Design principles of the yeast G1/S switch. PLoS Biol. 2013, 11 (10): e100167310.1371/journal.pbio.1001673.PubMed CentralView ArticlePubMedGoogle Scholar
 Fisher D, Krasinska L, Coudreuse D, Novak B: Phosphorylation network dynamics in the control of cell cycle transitions. J Cell Sci. 2012, 125 (Pt 20): 47034711.View ArticlePubMedGoogle Scholar
 Soloveichik D, Seelig G, Winfree E: DNA as a universal substrate for chemical kinetics. PNAS. 2010, 107 (12): 53935398. 10.1073/pnas.0909380107.PubMed CentralView ArticlePubMedGoogle Scholar
 Chen YJ, Dalchau N, Srinivas N, Phillips A, Cardelli L, Soloveichik D, Seelig G: Programmable chemical controllers made from DNA. Nat Nanotechnol. 2013, 8: 755762. 10.1038/nnano.2013.189.PubMed CentralView ArticlePubMedGoogle Scholar
 Dodd IB, Micheelsen MA, Sneppen K, Thon G: Theoretical analysis of epigenetic cell memory by nucleosome modification. Cell. 2007, 129 (4): 813822. 10.1016/j.cell.2007.02.053.View ArticlePubMedGoogle Scholar
 Hara M, Abe Y, Tanaka T, Yamamoto T, Okumura E, Kishimoto T: Greatwall kinase and cyclin BCdk1 are both critical constituents of Mphasepromoting factor. Nat Commun. 2012, 3: 1059PubMed CentralView ArticlePubMedGoogle Scholar
 Angluin D, Aspnes J, Eisenstat D: A simple population protocol for fast robust approximate majority. Distrib Comput. 2008, 21 (2): 87102. 10.1007/s004460080059z.View ArticleGoogle Scholar
 Soloveichik D: Robust stochastic chemical reaction networks and bounded Tauleaping. J Comput Biol. 2009, 16 (3): 501522. 10.1089/cmb.2008.0063.PubMed CentralView ArticlePubMedGoogle Scholar
 Novak B, Tyson JJ: Numerical analysis of a comprehensive model of Mphase control in Xenopus oocyte extracts and intact embryos. J Cell Sci. 1993, 106 (Pt 4): 11531168.PubMedGoogle Scholar
 Ghosh R, Tomlin CJ: Lateral inhibition through DeltaNotch signaling: A Piecewise affine hybrid model. HSCC. 2001, LNCS 2034: 232246.Google Scholar
 Jörg D, Gabriel J: What is a Petri Net? Informal answers for the informed reader. Unifying Petri Nets. 2001, LNCS 2128: 125.Google Scholar
 Iancu B, Czeizler E, Czeizler E, Petre I: Quantitative refinement of reaction models. Int J Unconv Comput. 2012, 8 (5–6): 529550.Google Scholar
 Danos V, Feret J, Fontana W, Harmer R, Krivine J: Abstracting the Differential Semantics of RuleBased Models: Exact and Automated Model Reduction. Proceedings of the 25th Annual IEEE Symposium on Logic in Computer Science. 2010, Los Alamitas California, Washington, Tokyo: IEEE Computer Society, Conference Publishing Services, 362381.Google Scholar
 Tschaikowski M, Tribastone M: Exact fluid lumpability in Markovian process algebra. Theor Comput Sci. in press,http://dx.doi.org/10.1016/j.tcs.2013.07.029,
 Fontana W, Schuster P: Continuity in evolution: on the nature of transitions. Science. 1998, 280: 14511455. 10.1126/science.280.5368.1451.View ArticlePubMedGoogle Scholar
 Hárs V, Tóth J: On the Inverse Problem of Reaction Kinetics. Colloquia Math. Soc. Jànos Bolyai. 30: Qualitative Theory of Differential Equations. 1979, Szeged: Bolyai János Matematikai Társulat, Hungary, 363379.Google Scholar
 Soloveichik D, Cook M, Winfree E, Bruck S: Computation with finite stochastic chemical reaction networks. Nat Comput. 2008, 7 (4): 615633. 10.1007/s110470089067y.View ArticleGoogle Scholar
Copyright
This article is published under license to BioMed Central Ltd. This is an Open Access article distributed under the terms of the Creative Commons Attribution License (http://creativecommons.org/licenses/by/4.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly credited. The Creative Commons Public Domain Dedication waiver (http://creativecommons.org/publicdomain/zero/1.0/) applies to the data made available in this article, unless otherwise stated.