Large, naturally evolved biomolecular networks typically fulfil multiple functions. When modelling or redesigning such systems, functional subsystems are often analysed independently first, before subsequent integration into larger-scale computational models. In the design and analysis process, it is therefore important to quantitatively analyse and predict the dynamics of the interactions between integrated subsystems; in particular, how the incremental effect of integrating a subsystem into a network depends on the existing dynamics of that network. In this paper we present a framework for simulating the contribution of any given functional subsystem when integrated together with one or more other subsystems. This is achieved through a cascaded layering of a network into functional subsystems, where each layer is defined by an appropriate subset of the reactions. We exploit symmetries in our formulation to exhaustively quantify each subsystem’s incremental effects with minimal computational effort. When combining subsystems, their isolated behaviour may be amplified, attenuated, or be subject to more complicated effects. We propose the concept of mutual dynamics to quantify such nonlinear phenomena, thereby defining the incompatibility and cooperativity between all pairs of subsystems when integrated into any larger network. We exemplify our theoretical framework by analysing diverse behaviours in three dynamic models of signalling and metabolic pathways: the effect of crosstalk mechanisms on the dynamics of parallel signal transduction pathways; reciprocal side-effects between several integral feedback mechanisms and the subsystems they stabilise; and consequences of nonlinear interactions between elementary flux modes in glycolysis for metabolic engineering strategies. Our analysis shows that it is not sufficient to just specify subsystems and analyse their pairwise interactions; the environment in which the interaction takes place must also be explicitly defined. Our framework provides a natural representation of nonlinear interaction phenomena, and will therefore be an important tool for modelling large-scale evolved or synthetic biomolecular networks.
Complex biochemical reaction networks serve a broad variety of different tasks within the cell. Systems Biology researchers apply a range of systems analysis techniques to these networks to identify and model functional subsystems and their interaction structure. In the context of biomolecular networks, the subsystems that can be identified often have biological interpretations: for example, the heat shock response and the chemotaxis pathways represent two functional subsystems within a model describing the complete biomolecular reaction network of Escherichia coli; the synthesis pathways of individual products represent distinct functional subsystems within a metabolic network; and so on. Other functional subsystems may also have system-theoretic interpretations: for example, interacting, distributed feedback control mechanisms; or subsystems which can sense, compute, or actuate on the cell and its environment. This tangle of different objectives within the same network leads to trade-off situations: evolutionary or synthetic changes to one functional subsystem can lead to declining performance or unexpected side effects with respect to another. A fundamental challenge of Systems Biology is to not only establish the behaviour of each functional subsystem in isolation, but also to understand how they dynamically influence one another. This problem is particularly acute when applying the modelling and analysis tools of Systems Biology to adapt and redesign modular biomolecular networks in Synthetic Biology [1–3]. The dynamics of many functional subsystems, whether evolved biochemical networks or synthetic devices, often do not proceed as modelled when integrated into a cell due to their interactions with one another and the environment of their cellular host. Possible sources of nonlinear interactions between pairs of functional subsystems and the cellular environment include retroactivity in genetic [4–6] and signalling [7] networks, crosstalk between parallel signalling pathways [8,9], and the coupling of multiple transcription or translation rates through competition for shared resources [10–12]. In each of these settings, the change in input–output behaviour of a given subsystem upon integration with its context is examined. There are two complementary goals of this paper. First, we investigate the behaviour of each subsystem between the two extremes of ‘isolated’ or ‘integrated’, when it is integrated with any subset of the other subsystems. The second goal, which is achieved as a consequence of the first, is to then systematically quantify each of the pairwise interactions between the network’s subsystems. The approach we will take in this work is to define a functionality as a group of reactions which corresponds to an identified functional subsystem of a biomolecular network. The reactions that determine each functionality can be selected either through biological insight, or by applying existing computational approaches such as elementary flux mode (EFM) analysis [13–15] (see also Results). We characterise the behaviour, or effect, of a functionality as the solution of an ordinary differential equation (ODE) model determined by the particular group of reactions. This approach exploits the recently-introduced decomposition technique known as layering [16,17]. As depicted in Fig 1, such an approach is distinct from established modular approaches to network decomposition, which are characterised by identifying sets of species with a high connectivity inside the module, and significantly lower connectivity to species in other modules [18–26]. While often many species and reactions in a given network are implicated in multiple network functions, these modular approaches generally do not allow for such a high degree of overlap between modules. For example, if a network of two pathways responds to two external signals with a single output species, a modular decomposition of this network requires the common output species to be assigned to a module representing exactly one of the pathways, or potentially to an additional separate module. Either way, the input–output behaviour of both pathways cannot be easily defined. However, in the layered framework, the common output is associated with both layers, and hence the output of each layer can be defined in terms of its biological function. Thus, in some cases, layers are preferable to modules for defining the functional subsystems of the network, since the layered framework explicitly allows for overlap in species and reaction subsets [16,17], as will be illustrated further in ‘Mapping Functionalities to Layers’ below. Most importantly, we make it explicit that the behaviour of any functionality also depends on the other functionalities with which it is integrated, to which we refer as the context of the functionality. This contextual dependence is formalised by developing a notational framework that will unambiguously define a functionality’s behaviour in a particular context. The resulting concept of conditional dynamics will be key to our understanding of each functionality as being defined only in the context of others, allowing us to systematically investigate the interdependence of an entire network’s functionalities. The subsequent aim of this framework is to characterise all of the interactions between each pair of functionalities, each of which is also context-dependent. Our approach is a formalisation and extension of previous investigations into additive (i. e. independent), synergistic, or antagonistic subsystem interactions. Examples of these phenomena include the cytokine secretion by macrophages in response to stimulation with different sets of ligands [27], the response of bacteria to different combinations of drugs [28], or calcium signalling responses to different stimuli [29]. Importantly, we demonstrate how the strength and the type of interactions between functionalities depends on mediated indirect interactions with the other functionalities comprising their context. The relationship of our approach to the concepts in [27–29] is further discussed in the section “Calculating With Layers”. In addition to the previous literature on context-dependent dynamics, our theoretical framework is also related to steady-state methods. For instance, the third of our examples will exploit EFMs, a technique designed to analyse the steady-state flux distribution in metabolic networks [13–15]. Furthermore, modular and hierarchical control/response analysis is concerned with the different behaviour of subsystems both in isolation and integrated in larger systems, with particular reference to the steady-state responses of biochemical networks to parameter perturbations [6,30–34]. The key distinction between our method and these is that we analyse the dynamics of kinetic models [35,36] represented by sets of ODEs, rather than steady states. Moreover, our method is not based on linearisation, allowing us to adequately capture nonlinear interactions between functionalities. Quantifying the dynamic interactions of each pair of functionalities in all possible contexts requires multiple ODE simulations; for its practical applicability, it is important to minimise the computational effort involved. This paper is structured as follows: in the Methods section we show how to use a layered decomposition to identify the incremental effect of a functionality, making its context-dependence explicit. We continue by defining the interdependence, or mutual dynamics, between any two functionalities. We summarise this interdependence by the incompatibility and the cooperativity between functionalities. In the final part of the Methods section we describe how to analyse all functionalities and their dependencies with minimal computation. We demonstrate our method on three familiar biomolecular networks in the Results section. The first example is of two signalling pathways with two crosstalk mechanisms, in which we use our approach to quantify the nonlinear interactions between crosstalk mechanisms. In the second example we analyse an unstable pathway stabilised by two integral feedback loops, finding the interactions between each controller and the pathway, and also between the controllers. Finally, we consider the glycolytic pathway in Saccharomyces cerevisiae, with functionalities defined by an EFM analysis. We apply our approach to compare how different knock-out strategies in metabolic engineering influence the yield of ethanol, industrially relevant in biofuel production.
Consider a biochemical reaction network with NX species Xi of time-varying concentrations xi (t), taking part in NR reactions R1, …, RNR, each of which proceeds at the concentration-dependent rate vj (x1, …, xNX) for j = 1, …, NR. Forming vectors v = (v1, …, vNR) T and x = (x1, …, xNX) T gives an ODE model of the system x ˙ (t) = S v (x (t) ), x (0) = x 0, (1) where the stoichiometric matrix S maps reaction rates to the rate of change of concentrations. The layered decomposition strategy [16,17] defines NL new stoichiometric matrices S 1, …, S N L such that S = S 1 + ⋯ + S N L, and defines NL associated state variables xl taking values x l (t) ∈ ℝ N X (see Fig 1B). Each layer’s state xl has dynamics x ˙ l (t) = S l v (x 0 + x l (t) + ∑ k ≠ l x k (t) ), (2) from initial conditions xl (0) = 0, for l = 1, …, NL. The original state’s dynamics are recovered by summing the layers’ states x (t) = x 0 + ∑ l = 1 N L x l (t). Denote by r = rank (S) the dimension of the original system, which defines the dimension of the manifold in ℝ N X in which x (t) evolves. It follows that rl = rank (Sl) defines the dimension of the state space of each layer. Hence, even though the state space of each layer is also embedded in ℝ N X, each layer is a lower-dimensional system than the original system if rl < r. In our previous work, we have applied this decomposition strategy by choosing the matrices Sl to reflect timescale separation [17], and to reflect the propagation of steady-state responses to parametric perturbations [16]. A feature of both of these approaches was that, in the form (2), each layer’s dynamics depend on all other layers’ states (as in Fig 1B, for example). Consequently, all layers had to be numerically integrated together, and the effect of one specific layer on all others could not be easily determined. Also, the approach was constrained to define layers by strict partitions of the reaction set, somewhat limiting its flexibility to capture the widest possible range of functional subsystems. In this article, we significantly extend the layering framework in two ways. First, we introduce the concept of functionalities, which are possibly overlapping sets of reactions working together for a common purpose. To enforce a cascade structure between the functionalities, we adapt the layered dynamics corresponding to each functionality depending on its position in the cascade. The following section will use this cascade structure to define the incremental dynamic effect of each functionality. Let a functionality Fi of a network be defined for i = 1, …, NL to be a subset of N R i reactions Fi ⊆ {R1, …, RNR} necessary to fulfil a given task of the network, where superscript integers index functionalities and their properties. It is assumed for the remainder of this section that these subsets are given, and that all reactions take part in at least one functionality. The question of how to choose each subset Fi ⊆ {R1, …, RNR} remains out of the scope of this work. Nevertheless, there are numerous non-modular decomposition strategies taken in recent related research that we can use to justify this definition of a functionality. For example, Oishi and Klavins [37] identify control blocks as specific groups of reactions, connected by shared species. Kurata et al. [38] identified reaction groups forming ‘flux modules’ in the Escherichia coli heat shock response system. Similarly, the decomposition of signalling networks into component pathways exhibiting crosstalk [8] also identifies functionalities as groups of reactions. Finally, we can also consider elementary flux modes (EFMs) of metabolic networks [13] as being sets of reactions with the specific ‘task’ of converting one or more substrates into given products. Several of these examples are explored further in the Results section of this paper. In this section, we will assume that the functionalities are ordered by their index F 1, …, F N L. We first identify the dynamics of the isolated functionality F1 as the dynamics of a biomolecular network consisting of only the reactions associated with F1. We then identify the conditional dynamics of the next functionality in the cascade as the effect of extending the pre-existing network with the reactions in the new functionality. First consider, without loss of generality, the network defined by only the subset of reactions making up functionality F1 ⊆ {R1, …, RNR}, taken in isolation from the other reactions. For given initial conditions x0, we now identify the isolated dynamics of this functionality as the solution to the layer x ˙ 1 (t) = S 1 v (x 0 + x 1 (t) ), x 1 (0) = 0. (3) Here, the stoichiometric matrix S1 is defined S j k 1 = { S j k R k ∈ F 1,0 otherwise, by copying the columns of the original stoichiometric matrix S in (1) corresponding to the reactions in F1 and setting the other columns to zero. We will denote this trajectory x1 = L (F1), where the notation L represents a map from the functionality F1 to the solution x1 of the dynamics (3) from initial conditions x1 (0) = 0. Note that L (F1) depends on the specific initial condition x0 of the network (3), which is in general distinct from the initial condition x1 (0) = 0 of the layer’s state. To make this dependence explicit, it is sometimes helpful (see Examples 2 and 3) to define a ‘zero layer’ F0 with constant trajectory L (F0) = x0. We can then make clear that L (F1) is dependent on the initial conditions by writing it as L (F1∣F0). The layered framework also implies that the absolute concentrations in this network are modelled by the translated trajectory x0 + L (F1∣F0). We next consider extending the functionality F1 by combining it with the reactions in F2. The extended network can be simulated through a similar process to the original network above, as follows. Define S1,2 as S j k 1,2 = { S j k R k ∈ F 1 ∪ F 2,0 otherwise, considering only the reactions in at least one of F1 or F2. Using S1,2 we can then simulate the layer corresponding to the extended network x ˙ 1,2 (t) = S 1,2 v (x 0 + x 1,2 (t) ), x 1,2 (0) = 0, (4) the solution of which can be written L (F1, F2∣F0) = x1,2. This denotes the trajectory of the combined functionalities F1 and F2. The fact that each of (3) and (4) are layers with states in ℝ N X implies that we can calculate the difference between each of the trajectories. This difference is clearly interpreted as the incremental effect of extending a network made up of the initial conditions F0 and the isolated functionality F1, by also including F2. We thus define L (F 2 | F 1, F 0) = L (F 1, F 2 | F 0) - L (F 1 | F 0) (5) as the conditional dynamics of F2, given the specified context of F1 and the initial condition layer F0. However, rather than simulating the layer (4) representing the combined functionalities, we may further exploit the layered framework described above to directly simulate L (F2∣F1, F0). Suppose we already have x1 = L (F1∣F0), found as the solution to the dynamics (3). We now define the layer x ˙ 2 (t) = S 1,2 v (x 0 + x 1 (t) + x 2 (t) ) - S 1 v (x 0 + x 1 (t) ), x 2 (0) = 0, (6) with S1 and S1,2 given above. Note that this layer is downstream of (3), since it depends on the state x1. It is clear from summing the vector fields in (3) and (6) that the sum (x1 + x2) of the layers’ states follows exactly the same dynamics as the combined network’s state x1,2 in (4). Thus, since x2 = x1,2 − x1 it follows that the dynamics (6) directly simulate L (F2∣F1, F0), with the input L (F1∣F0) simulated by (3). We can rewrite the dynamics (6) corresponding to the simulation of L (F2∣F1, F0) as x ˙ 2 (t) = S 2 v (x 0 + x 1 (t) + x 2 (t) ) + S 1 v a l t (x 0 + x 1 (t), x 2 (t) ) (7a) where S2 = S1,2 − S1 corresponds to the reactions in F2\F1, and v a l t (x 0 + x 1, x 2) = v (x 0 + x 1 + x 2) - v (x 0 + x 1) (7b) are the rates of the ‘altered reactions’: the rates of reactions in F1 which are modified by the presence of F2 (shown as broken green arrows in Fig 1C). This description allows us to see the degree to which F2 is ‘downstream’ of F1. For example, if valt = 0, then we can say that the reactions in F1 are independent of those in F2 and that F2 is strictly downstream of F1. Note that, especially for larger networks, many of the altered reaction rates in valt are zero and can be omitted (see Example 3), simplifying simulation. Given that the trajectory of x1 is already determined from simulating (3), we can simulate either (6) or (7a), using x1 (t) as a time-dependent input to obtain the conditional dynamics L (F2∣F1). The latter approach is taken in our examples (see Results section). The definitions above easily extend to larger combinations of functionalities. In full generality, we can consider the network defined by the combination of n1 functionalities F 1, …, F n 1, and its extension through the additional n2 functionalities F n 1 + 1, …, F n 1 + n 2. By writing F ‾ 1 = F 1 ∪ ⋯ ∪ F n 1 and F ‾ 2 = F n 1 + 1 ∪ ⋯ ∪ F n 1 + n 2, the definitions above can apply in the simulation of L (F 1, …, F n 1 | F 0) = L (F ¯ 1 | F 0), L (F n 1 + 1, …, F n 1 + n 2 | F 0, F 1, …, F n) = L (F ¯ 2 | F 0, F ¯ 1). Here we have defined a notation for the trajectory L (F 1, …, F n 1 ∣ F 0) of the biochemical network made up of the reactions which comprise an arbitrary combination of functionalities F 1, …, F n 1. We have also defined the change in trajectory L (F n 1 + 1, …, F n 1 + n 2 ∣ F 0, F 1, …, F n 1) incurred by extending that network with the additional reactions in F n 1 + 1, …, F n 1 + n 2. Finally, we have shown how to identify the dynamical systems that can simulate these trajectories. Consider the two layers L (F2∣F1) and L (F2) that both describe the effect of the functionality F2. This effect is different depending on the presence or absence of F1. The difference between these two trajectories defines how the presence of F1 changes the behaviour of F2; that is, the dependence of F2 on F1. We will now demonstrate how our layered analysis allows us to define the interdependence between two functionalities, thereby capturing the nonlinear effects arising from modelling a biomolecular network as being constructed from a combination of functional subsystems. In order to quantify the interactions between functionalities, we can exploit the layered formulation above. For simplicity, from this point on we suppress the F0 notation, with the acknowledgement that all of the trajectories depend on the system’s initial conditions L (F0) = x0. The definition of conditional dynamics in (5) implies that L (F1, F2) = L (F2∣F1) + L (F1). This represents a layered cascade, where the dynamics of an integrated network are the linear combination of the conditional dynamics of its functionalities. There are two natural questions associated with this approach. First, how is the contribution of functionality F2, considered in isolation, different from the conditional dynamics of F2 when integrated with F1? Secondly, how is the behaviour of the integrated F1, F2 network different from the linear combination of the isolated functionalities? That is, how different is L (F1, F2) from L (F1) + L (F2)? The answers to these two questions are the same. We denote the error incurred by approximating the integrated system as the linear combination of the isolated dynamics by the quantity M (F1; F2), defined as M (F 1; F 2) = L (F 1) + L (F 2) - L (F 1, F 2), which we call mutual dynamics. This can be interpreted as the nonlinearity that arises from integrating the two functionalities together. Note that, since M is defined symmetrically, we can use (5) to rewrite M as M (F 1; F 2) = L (F 2) - L (F 2 | F 1) = L (F 1) - L (F 1 | F 2). (8) Therefore M measures how the function of F2 (or F1) is changed when considered in the context of F1 (or F2). Thus M (F1; F2) is a symmetric measure of the interdependence between the two functionalities. We use (8) to calculate the the mutual dynamics between two functionalities, which requires us to first obtain either both trajectories L (F1) and L (F1∣F2), or alternatively both trajectories L (F2) and L (F2∣F1). These trajectories can be simulated, as described in the previous section, or calculated by the methods described in ‘Reducing Computational Burden’ below. We have been careful to make explicit through our Bayesian-style notation that the dynamics of all functionalities are context-dependent. It is also the case that the interdependence between any two functionalities is context-dependent. Therefore, we need to extend the definition of mutual dynamics to consider how the interdependence between F1 and F2 is dependent on the wider context of the network, which we denote by another functionality, F3. Similarly to the definition above, we can define the conditional mutual dynamics between F1 and F2, given F3, with the formula M (F 1; F 2 | F 3) = L (F 1 | F 3) + L (F 2 | F 3) - L (F 1, F 2 | F 3), to quantify the difference between the dynamics of the integrated and isolated functionalities F1 and F2, in the context of F3. As before, this can also be expressed in terms of layered dynamics as M (F 1; F 2 | F 3) = L (F 2 | F 3) - L (F 2 | F 1, F 3) = L (F 1 | F 3) - L (F 1 | F 2, F 3), to quantify how the effect of F2 on its context changes with the presence of F1, and vice versa. The geometric intuition underlying the conditional mutual dynamics can be seen in Fig 2. The key interpretation of M is that it captures the nonlinearities that arise from combining F1 and F2 into a single network, conditioned on F3 if necessary. The conditional mutual dynamics M (F1; F2∣F3) is a time-varying, vector trajectory. We can base on M the following time-varying scalar, which we call the incompatibility and denote I (F1; F2∣F3) with formula I (F 1; F 2 | F 3) = ∥ M (F 1; F 2 | F 3) ∥ ∥ L (F 1 | F 3) + L (F 2 | F 3) ∥. (9) Here, ‖. ‖ represents the Euclidean norm. In this paper we use the unweighted Euclidean norm, but in certain cases it might be appropriate to introduce a weight, for example if the concentrations of the species in a network are at different orders of magnitude. One might also decide to set the weight of certain intermediate species of limited interest to zero (see below). This incompatibility measures the relative size of the error made by approximating the integration of two functionalities as the sum of their individual behaviour. To gain some intuition about this number, we can consider a number of special cases. If I (F1; F2∣F3) = 0, then this indicates that the trajectory of the integrated functionalities is simply the sum of the isolated functionalities’ trajectories: L (F1, F2∣F3) = L (F1∣F3) + L (F2∣F3). If the incompatibility is nonzero but small then L (F1, F2∣F3) ≈ L (F1∣F3) + L (F2∣F3) is a reasonable approximation, since the incurred error is relatively small. However, if I is of significant size, then the dynamics of the integrated functionalities can be expected to significantly differ from their individual behaviours. Besides I, which measures the relative size of the mutual dynamics M, the direction of M is also important, since this determines if the integration of two functionalities together enhances or attenuates their individual dynamics, or causes other effects. We define the cooperativity C as the cosine of the angle between −M and the sum of the isolated layers: C (F 1; F 2 | F 3) = - M (F 1; F 2 | F 3) · (L (F 1 | F 3) + L (F 2 | F 3) ) ∥ M (F 1; F 2 | F 3) ∥ ∥ L (F 1 | F 3) + L (F 2 | F 3) ∥, (10) with ⋅ denoting the scalar product. See Fig 2 for a geometric representation of C. Note that when using a weighted norm, the scalar product should be weighted accordingly. Again, we consider a number of special cases. Suppose that C (F1; F2∣F3) equals or is close to minus one, so that M is approximately parallel to and pointing in the same direction as L (F1∣F3) + L (F2∣F3). In this case the integrated dynamics can be approximated as an attenuation of the isolated behaviour L (F 1, F 2 | F 3) = L (F 1 | F 3) + L (F 2 | F 3) - M (F 1; F 2 | F 3) ≈ (1 - I (F 1; F 2 | F 3) ) (L (F 1 | F 3) + L (F 2 | F 3) ). Conversely, if the cooperativity equals or is close to one, the two functionalities enhance each other, in the sense that we can approximate the integrated behaviour as an amplification of the isolated behaviours L (F 1, F 2 | F 3) ≈ (1 + I (F 1; F 2 | F 3) ) (L (F 1 | F 3) + L (F 2 | F 3) ). However, once the cooperativity C equals or is close to zero, the mutual dynamics M are orthogonal to L (F1∣F3) + L (F2∣F3). This means that when the functionalities are integrated, both isolated functionalities are maintained, but there are also additional interactions (with an effect of strength I) in directions orthogonal to the summed isolated dynamics. For example, suppose that the functionalities F1 and F2 correspond to sets of reactions responsible for mediating the cellular responses to two different input signals. By setting the influence of all but the common output (measured) species to zero, an incompatibility I (F1; F2) close to zero corresponds to an additive interaction of the input signals. If I (F1; F2) is larger, it may correspond to either a synergetic (for C (F1; F2) = 1) or antagonistic (for C (F1; F2) = −1) interaction (see e. g. [28]). Note that for a scalar output, orthogonal dynamics are not possible. For systems composed of many different functionalities, one might also take the time averages ⟨I⟩ and ⟨C⟩ of the incompatibility, respectively the cooperativity, over the simulation time ΔT to obtain single measures quantifying the interactions between layers: ⟨ I ⟩ (F 1; F 2 | F 3) = 1 Δ T ∫ 0 Δ T I (F 1; F 2 | F 3) d t ⟨ C ⟩ (F 1; F 2 | F 3) = 1 Δ T ∫ 0 Δ T C (F 1; F 2 | F 3) d t Although they are useful for obtaining a first impression of how functionalities interact, time averages should be carefully applied. They can hide transient interactions between functionalities, including potential sign changes of the state-dependent cooperativities (see Example 1). We have now identified how to measure the interdependence between two functionalities, which we have defined as the change in the dynamics of one functionality when the other is present. We have made explicit how this interdependence is itself dependent on the context of the rest of the network. In the remainder of this section we will describe how to minimise the computational burden incurred when calculating all possible interactions between functionalities. The notation L (F1) and L (F2∣F1) describing the map from a functionality (or set of functionalities) to the resulting trajectory simplifies the calculations we may wish to carry out to understand how the functionalities combine. Using the key definition of ‘conditional dynamics’ given by (5), we can prove a number of rules for combining layers which appear analogous to those known from Information Theory [39]. For example, for two random variables X and Y, the well-known quantities of joint entropy H (X, Y) = H (X) + H (Y∣X) and mutual information I (X; Y) = H (X) − H (X∣Y) are each definitions of the same form as those given above of conditional dynamics (5) and mutual dynamics (8) respectively. However, it is important to note that this similarity is only superficial, and any intuition gained by seeking analogies between our work and information theoretic concepts should be applied carefully. This caveat applies in particular to the two results below. Two lemmas allowing the quick combination of layer dynamics can be easily proved directly from the definitions of L (F1) and L (F2∣F1), their extensions to larger combinations of functionalities, and Eq (5). The first is an analogue of Bayes’ Rule, given by L (F 1 | F 2, F 3) = L (F 2 | F 1, F 3) + L (F 1 | F 3) - L (F 2 | F 3). (11) We will demonstrate how this rule can be used for quickly deducing the incremental effects of layers when combined in a different order. This is fundamental, since a natural ordering of the layers is generally not given. A second rule, which is analogous to Bayes’ Factor, is given by L (F 3 | F 1) - L (F 2 | F 1) ︸ Posterior Dynamics = L (F 1 | F 3) - L (F 1 | F 2) ︸ Bayes Factor B 32 | 1 + L (F 3) - L (F 2) ︸ Prior dynamics. (12) This rule applies when we have a choice between integrating two functionalities to the F1-only network. The difference in their effects is decomposed into the difference L (F3) − L (F2) between their isolated behaviours, summed with the difference L (F1∣F3) − L (F1∣F2) in the incremental effect of F1 on each. We can use these rules to compute all possible functionality combinations with a minimal amount of simulation. In order to answer particular biological questions, we may be interested in the incremental effect of a given functionality on a specific ‘base’ network, such as those described in Examples 1 and 2 in the Results section. In other cases we may be interested in all possible interactions between the functionalities, such as the situation in Example 3. In a biochemical network whose reactions are decomposed into NL functionalities, the latter case suggests that we must simulate all NL layers for each of the NL! different orderings of the functionalities, resulting in (NL + 1)! − NL! layers to be numerically solved. This burden can be significantly reduced using the calculation rules deduced above. The cascaded layers representing all possible orderings of functionalities can be arranged in an acyclic directed graph, shown in Fig 3. Each node represents the trajectory arising from the incremental addition of a new functionality, given those already present. The graph is organised into levels, corresponding to the position of the new layer in the sequence. The root of the layering graph (referred to as Level 0) represents the given initial conditions x0. Each of the subsequent levels l = 1, …, NL consists of (NL−l+1) (NLl−1) nodes. Each node represents the dynamics of a functionality Fi conditioned on a subset of size l − 1 of the remaining functionalities Fj, j ≠ i. A directed edge from a node in Level l to a node in Level l + 1 exists if the node in Level l + 1 is conditioned on all functionalities taking part in the node in Level l. Each directed path from Level 0 to Level NL (i. e. from the root to a leaf) represents one of the NL! possible orderings of functionalities. By adding up the layer dynamics corresponding to the nodes in each path, the trajectory of the complete network is obtained. Each node in Level 1 represents the dynamics of an isolated functionality. The dynamics represented by the nodes at the leaves of the layering graph are also of specific interest. Multiplying the layer dynamics in a given node in Level NL (corresponding to a particular functionality Fi) by −1 gives the change in the dynamics that results from removing Fi from the system while keeping all other functionalities intact (see Example 3). To obtain the dynamics of all layers for all orderings of functionalities, we can use rules (11) and (12) to exploit certain symmetries in the layering graph and reduce the number of numerical integrations of layer ODE systems. If all of the dynamics in Level l − 1 are already known, it is only necessary to numerically solve (N L l) layer dynamics in Level l. The rest of the level can then be deduced using Bayes’ rule (11). For example, suppose we have the trajectories of all nodes in Layer 1, and have simulated L (F1∣F0, F2) in Level 2. Using (11) we can deduce L (F2∣F0, F1) without having to simulate again. In fact, when considering a network with NL functionalities, only ∑ l = 1 N L (N L l) = 2 N L - 1 layers have to be numerically integrated (see Example 1), compared to (NL + 1)! − NL! integrations necessary when analyzing all orderings separately. Although the number of integrations still exponentially grows with NL, it nevertheless becomes possible to analyse systems with up to ten or more distinct functionalities in reasonable time. For NL = 10, integrating every layer would require more than 3. 5 ⋅ 104 times as much computational time as is required by exploiting (11). Furthermore, the symmetry of the layering graph allows a certain degree of freedom in choosing which layers to simulate, and which to deduce from (11). In particular, one should always choose to simulate the simplest layers, with respect to the number of states or reactions in the respective ODE system. Fig 3 shows the situation where NL = 3: we need to simulate 23 − 1 = 7 layers. We have decided to, wherever possible, simulate F3, ahead of F2, ahead of F1. After simulating the isolated behaviour L (F1∣F0) of the functionality F1, we can calculate its behaviour in any other context without having to simulate it again. For example, L (F1∣F0, F2) can be calculated, using (11), as L (F2∣F0, F1) + L (F1∣F0) − L (F2∣F0). Similarly, substituting the resulting trajectory into (11) once more means that we can calculate L (F1∣F0, F2, F3) as the linear combination L (F3∣F0, F1, F2) + L (F1∣F0, F2) − L (F3∣F0, F2). As can be seen by the indicated nodes in Fig 3, only one layer corresponding to F1 had to be numerically simulated, whereas we simulated layers corresponding to F2 twice, and layers corresponding to F3 four times. Thus we can reduce computation time even further by avoiding the repeated simulation of high-dimensional layers.
The high osmolarity glycerol and the pheromone response mitogen-activated protein (MAP) kinase pathways in S. cerevisiae share the common species Ste11 [40]. Such a common species can constitute a mutually excitatory crosstalk mechanism by cross-activating one pathway upon activation of the other. However, McClean et al. [8] observed that only one pathway responds, and deduced that a second, mutually inhibitory crosstalk mechanism exists. They constructed a model, available as Model 115 in the BioModels Database [41], that includes both crosstalk mechanisms. As an initial example of our approach we implemented a layered decomposition of this model, shown in Fig 4. The mathematical description of the complete model and of all layers can be found in the Supporting Information. We applied our framework to systematically investigate the function of each of the crosstalk mechanisms. We first considered the effects of each crosstalk mechanism on the signalling pathways. This gives some indication of the function of each crosstalk mechanism, but our framework takes this analysis further. The effect of excitatory crosstalk on the network is altered by the presence of inhibitory crosstalk, and vice versa. The conditional mutual dynamics between the two crosstalk mechanisms, integrated with the crosstalk-free network, can be used to quantify their interdependence. In the model, we set the strength of the mutually excitatory crosstalk to ka = 0. 1 and the strength of the mutually inhibitory crosstalk to kd = 1, corresponding to a monostable network (see Fig. 1 in [8]). We assume zero initial conditions. The three reactions R1, R2, and R3 represent the first signal transduction pathway, which affects species concentrations X1–X3, so that we define the functionality F1 = {R1, R2, R3}. Similarly, we choose F2 = {R4, R5, R6} to represent the second signal transduction pathway, affecting species Y1–Y3. The mutually inhibiting crosstalk between the two pathways comprises F3 = {R7, R8}. Finally, the mutually excitatory crosstalk is given by F4 = {R9, R10}. We first calculated all of the possible incremental effects of each functionality Fi integrated with each possible subset {Fj ∣ j ≠ i} of the others. Since we have NL = 4 functionalities we only need to simulate 2 N L − 1 = 15 layers to be able to calculate all NL ⋅ NL! = 96 layer dynamics for all orderings of the functionalities. Fig 5 shows the trajectories corresponding to several of these nodes. The plots were generated for time-varying inputs S1 and S2 given by [ S 1 (t), S 2 (t) ] = { [ 0,0 ] 0 ≤ t < 20, [ 5,0 ] 20 ≤ t < 40, [ 5,5 ] 40 ≤ t < 80, [ 0,5 ] 80 ≤ t < 100, [ 5,5 ] 100 ≤ t. 10. 1371/journal. pcbi. 1004235. g005 Fig 5 Layer dynamics of the crosstalk example. Dynamics of the layers for switching between different input combinations: (-/-) S1 = S2 = 0; (+/-) S1 = 5, S2 = 0; (+/+) S1 = S2 = 5; and (-/+) S1 = 0, S2 = 5. The top left figure corresponds to the dynamics of the complete network, the central and right top figure to the dynamics of the two isolated signalling pathways without any crosstalk. If the mutually inhibitory crosstalk L (F3∣F1, F2) is added first, the network becomes bistable, then monostable again by then including the mutually excitatory crosstalk L (F4∣F1, F2, F3). Adding first the excitatory crosstalk L (F4∣F1, F2) leads to strong cross-activation of the pathways, which is significantly weakened by the mutual inhibitory crosstalk L (F3∣F1, F2, F4). The switching times between the input combinations were chosen such that, at the end of each period, the species of the network have converged to their corresponding steady-state concentrations. The top-left plot in Fig 5 shows the response L (F1, F2, F3, F4) of the entire network to this input pattern. Each of the other plots depicts the incremental effects of integrating a certain functionality with a group of others. Fig 6 depicts the trajectories of the mutual dynamics between several pairs of functionalities, also calculated from the 15 simulations of layered ODE systems. Below these trajectories we plot their corresponding incompatibilities I and cooperativities C. We consider the basic network comprised only of F1 and F2, where neither crosstalk mechanism is active. From L (F1) and L (F2) in Fig 5 we see that after t = 20 (respectively, t = 40) the concentrations of species X1–X3 (respectively, Y1–Y3) quickly saturate. It is easily shown that the mutual dynamics between the two isolated pathways M (F1; F2) = 0. This means that there is no interdependency between the pathways, and their integrated dynamics equal the sum of their isolated dynamics. Of course, this conclusion is intuitively clear because we are considering a network with neither crosstalk mechanism active. We are now in a position to investigate what each crosstalk mechanism does to this basic network, by identifying the effects of each of F3 and F4 in turn. We first identify the effect of the mutual inhibitory crosstalk F3 on the basic network. This is depicted in Fig 5 by L (F3∣F1, F2). The input signals S1 = S2 = 5 are both ‘on’ during two time intervals t ∈ [40,80] and t ≥ 100. In each of these time intervals, we can observe in L (F3∣F1, F2) two different steady states of X1–X3 and Y1–Y3. Thus one effect of mutual inhibition is that the resulting network is bistable. We can also observe that the X1–X3 values of L (F3∣F1, F2), which show the effect of F3 on the Xi concentrations, become sufficiently negative during t ∈ [80,100] that they cancel out the positive values of X1–X3 in L (F1, F2). Hence, we can conclude that, by integrating mutual inhibition with the basic network, the removal of S1 at t = 80 now causes the corresponding pathway to deactivate. We can also analyse the interaction between mutual inhibition and the basic network in terms of the mutual dynamics M (F3; F1, F2), as depicted in Fig 6. For t ≥ 40, the cooperativity C (F3; F1, F2) is always negative, and the incompatibility I (F3; F1, F2) is high. Thus, the effect of integrating F3 with (F1, F2) is to strongly attenuate the levels of all species from their saturated state. We next consider the effect on the basic network of mutually excitatory crosstalk F4, depicted by L (F4∣F1, F2) in Fig 5. We will focus in particular on the time interval t ∈ [20,40], where S1 is first activated. For this time interval, the concentrations of Y2 and Y3 increase, while the concentrations of X1–X3 are transiently reduced. Thus the effect of including mutual excitation is that a non-zero input S1 is sufficient to activate the second pathway, as well as the first. As in the previous case, we can also analyse the interaction between mutual excitation and the basic network in terms of the mutual dynamics M (F4; F1, F2), as depicted in Fig 6. In this time interval, the cooperativity C (F4; F1, F2) is briefly negative before returning to zero, while the incompatibility I (F4; F1, F2) increases monotonically and converges to around 0. 81. Negative cooperativity implies that mutual excitation attenuates the isolated dynamics, although only marginally since I is small. This attenuation is also only transient; C approaches zero and I approaches approximately 0. 81 as t→40. Thus, by the end of this time interval, integrating mutual excitation creates a significant additional effect in a direction orthogonal to the isolated dynamics. This is intuitively clear, since the incremental effect of mutual excitation is the additive excitation of Y2 and Y3, orthogonal to the saturation of X1–X3 during t ∈ [20,40]. We have identified the effect of each crosstalk mechanism on the basic network. However, we can see from the plots of L (F3∣F1, F2, F4) and L (F4∣F1, F2, F3) that the function of each crosstalk mechanism (when defined as its incremental effect on a network) is very different when the other crosstalk mechanism is present. We can see from the first of these trajectories that one effect of integrating inhibition into the cross-activated network is to effectively insulate the Y1–Y3 pathway, since the excitation shown by L (F4∣F1, F2) during t ∈ [20,40] is cancelled out by the values of L (F3∣F1, F2, F4). Furthermore, on comparing L (F3∣F1, F2, F4) with L (F3∣F1, F2) we can see that mutual inhibition has remarkably different incremental effects depending on whether or not mutual excitation is present. This can be quantified by observing the trajectory of M (F 3; F 4 | F 1, F 2) = L (F 3 | F 1, F 2) - L (F 3 | F 1, F 2, F 4) and the associated C and I values in the right-hand plots of Fig 6. For t ∈ [20,40], the value of I (F3; F4∣F1, F2) is large, and C (F3; F4∣F1, F2) is close to −1. This confirms the assertion above that, in this time interval, the interdependence between crosstalk mechanisms causes them to approximately cancel each other out, so that the behaviour of the entire system is much like that of the isolated pathways. We can observe that, in this time interval, the effect of mutual inhibition on the basic network L (F3∣F1, F2) is zero. However, non-zero mutual dynamics M (F3; F4∣F1, F2) in this time interval means that the presence of mutual excitation causes the inhibition crosstalk to have a non-zero function. To summarise, we have demonstrated how to use the concepts of conditional dynamics L (Fi∣F1, F2) and mutual dynamics M (Fi; F1, F2) to analyse the effect of each crosstalk mechanism on a pair of isolated pathways. We extended this analysis by using M (F3; F4∣F1, F2) to quantify how, by integrating the crosstalk mechanisms together, they influence one another’s isolated functions. This example illustrates the application of our layered framework to a toy network, comprised of a biomolecular cascade containing two integral feedback motifs (Fig 7). This toy network was adapted from a motif identified in [42] as sufficient to exhibit perfect adaptation to changing external input concentrations, and underlying the ability of the chemotaxis pathway in E. coli to show near-perfect adaptation to extracellular chemoattractant concentrations [43–45]. The network consists of two intermediates, Y1 and Y2. The species Y1 is produced with a time-varying rate k (t) depending on the concentrations of one or more external species. Intermediate Y1 is irreversibly converted to Y2 at rate V1 (y1) and Y2 is consumed with rate V2 (y2), each following Michaelis-Menten kinetics (see Fig 7). For certain production rates of Y1, V1 is close to saturation so that the concentration of Y1 eventually becomes unstable for high rates k (t). To stabilise the network and achieve perfect adaptation of the concentrations of Y1 and Y2 to slowly changing production rates k (t), the network needs to be controlled. Such controllers can be implemented through additional Michaelis-Menten reactions Va (y1) and Vb (y2) forming species A and B from Y1 and Y2, respectively. A and B are consumed with Michaelis-Menten rates V−a (a) and V−b (b) respectively, both of which are assumed to be close to saturation. The two control loops are closed by non-competitive inhibition of the production of Y1 and of the conversion of Y1 to Y2 modelled by adapting the production rates with the factors ϕa (a) and ϕb (b) respectively, as shown in Fig 7. The dynamics of the system can be modelled by the following ODEs: [ y ˙ 1 y ˙ 2 a ˙ b ˙ ] ︸ =: x ˙ = [ 1 - 1 0 - 1 0 0 0 0 0 0 1 - 1 0 0 0 - 1 0 0 0 0 0 0 1 - 1 0 0 0 0 0 0 0 0 0 0 1 - 1 ] ︸ =: S [ k (t) ϕ a (a) V 1 (y 1) ϕ b (b) V 2 (y 2) V a (y 1) V a (y 1) V - a (a) V b (y 2) V b (y 2) V - b (b) ] ︸ v (x), with V i (x) = k i x m i + x, ϕ i (x) = 1 p i + x. Reactions Va and Vb have two distinct effects on the pathway: first, they decrease the concentration of Yi independently of any interaction of the controller; and second, they increase the concentration of A and B and are thus the means by which the controller observes the state of the network. To distinguish these two effects mathematically, in the dynamics above we have distributed each of Va and Vb into two reactions with the same rates, one only affecting the controller, and the other only the pathway. To analyse the system, we define three functionalities, corresponding to the uncontrolled pathway F1 (columns 1–4 and 7 in S) and the two controllers F2 (columns 5 and 6) and F3 (columns 8 and 9). The initial condition x0 was set to represent the steady state concentrations of the species for k = 3. Note that, in order to make explicit the dependency of the layered dynamics on this initial condition, we created a ‘zero layer’ F0 with constant trajectory L (F0) = x0 as described in the Methods section. The ODEs describing the layers’ dynamics for each ordering of the functionalities, the parameter values, as well as the precise initial conditions for this example can be found in the Supporting Information. We then determined the dynamics of the layers when applying a step change of the Y1 production rate from k = 3 to k = 10 at t = 50. Fig 8 displays the dynamics of the pathway, the conditional dynamics of the controllers given the pathway for both orderings of the controllers, and the complete dynamics of the system. Similarly to the previous example, we will consider a basic network consisting of only the uncontrolled pathway, F1. We are then able to identify the effect of each controller as the conditional dynamics of F2 and F3 given the basic network. However, the effect on the basic network of both controllers will be different to the sum of each controller’s individual effect. We will thus compute the mutual dynamics between the two controllers in order to analyse how the controllers interact with one another when both are integrated in the context of the basic network. We first consider the isolated dynamics L (F1∣F0) of the uncontrolled pathway which makes up the basic network, depicted in Fig 8. After the step increase at t = 50, the two reactions decreasing Y1 quickly saturate. The isolated system is therefore unstable, and the Y1 component of L (F1∣F0) increases to infinity at a positive constant slope. We will now use this basic network to define the function of each of the controllers. Integrating the controller F2 with the pathway stabilises the network. This is indicated by the trajectory of L (F2∣F0, F1) in Fig 8. In particular, the Y1 component of this trajectory decreases to minus infinity with approximately −1 times the rate of increase of Y1 in the isolated pathway. Another way to observe this stabilising effect is to consider the mutual dynamics M (F1; F2) between the pathway and the controller, shown in Fig 9. The inconsistency I (F1; F2∣F0) between pathway and controller approaches one, while the cooperativity C (F1; F2∣F0) minus one. We interpret these values as the two functionalities cancelling one another out; in other words, F2 stabilising F1. We next consider the effect of the second controller F3, depicted in Fig 8 by L (F3∣F0, F1). The Y1 component of L (F3∣F0, F1) also diverges to positive infinity, albeit at a much slower rate than that of the basic network. Intuitively, the second controller cannot stabilise the basic network, since it has no means to observe the unstable Y1 component. Instead, its interventions further destabise the pathway. We can quantify this through the mutual dynamics between the controller and pathway. Indeed, since C (F1; F3∣F0) ≈ 1, the second controller amplifies the unstable dynamics. However, since I (F1; F3∣F0) ≪ 1, this effect is only minimal. We now consider the interaction between the two controllers when combined with the uncontrolled pathway. Interestingly, the effect of the second controller when combined with the first is only transient, as shown by L (F3∣F0, F1, F2) in Fig 8. However, the mutual dynamics M (F2; F3∣F0, F1) between the two controllers have a Y1 component increasing slowly to infinity, and the cooperativity between the controllers C (F2; F3∣F0, F1) approaches 1. This positive cooperativity indicates that the first controller increases its interventions to stabilise the network when the second controller is present. That the two controllers given the pathway act cooperatively might initially appear surprising, since the first controller stabilises the network while the second destabilises the network. This indicates the utility of our notation in allowing a rigorous quantification of how functionalities interact, given their environment. Consider the following three cases for how the two controllers F2 and F3 interact. First consider the interaction M (F2; F3∣F0, F1) between the two controllers when integrated with the pathway. Since the presence of F3 increases the control action by F2 necessary to stabilise the pathway, the cooperativity given the pathway C (F2; F3∣F0, F1) approaches 1, as described above. Next, consider the interaction M (F2; F3∣F0) between the two controllers without the pathway. They cannot influence each other since the respective subnetworks are not connected. Consequently they do not interact and hence M (F2; F3∣F0) = C (F2; F3∣F0) = I (F2; F3∣F0) = 0. Finally, consider the interaction M (F2; (F1, F3) ∣F0) between the first controller F2 and the unstable network comprised of the pathway F1 extended with F3. The first controller stabilises the extended pathway, which is represented by the cooperativity C (F2; (F1, F3) ∣F0) approaching minus one, as F2 acts to attenuate the instability of F1, F3. Thus, we have demonstrated that the rigorous definition of interdependence using our cascaded layering framework allows us to identify which components of a network amplify or attenuate each other, by defining the components which interact and the context in which they do so. Glycolysis is a central ten-step process in most organisms responsible for the production of energy in form of ATP and NADH by catabolism of glucose and other sugars (see [46], p. 88 ff.). Besides being one of the best studied metabolic pathways, it is also an important starting point for biotechnological processes [47]. For non-growing S. cerevisiae cells, a kinetic model of the glycolytic pathway was created by Hynne et al. [48], including fermentation, glycerol production, lactonitrile and glycogen formation, and cellular import and export processes for glucose and other metabolites. We will use this model (available at the BioModels Database [41], model 61) to exemplify how to apply our layering approach in the context of metabolic engineering. In this section we use Elementary Flux Modes (EFMs) [13], minimal functional pathways which can carry non-zero fluxes at steady state and which fulfil positivity constraints for irreversible reactions. From the unique set of EFMs of a network, all steady-state flux distributions can be obtained by non-negative linear combinations of EFMs. Due to their simplicity, EFMs can often be associated with certain elementary ‘tasks’ of a network, like the production of one or more final products from various available extracellular substrates. Thus EFMs are natural candidates to represent functionalities. The model in [48] includes the dynamics of the cofactors NAD+, NADH, AMP, ADP, and ATP, with the consumption of ATP by the rest of the cell modelled by first order mass action kinetics. The cofactor concentrations can be interpreted as control inputs to the junctions of the glycolytic pathway dynamically channelling the metabolite flux into the different branches depending on cellular requirements, e. g. during hyperosmotic conditions [49]. The feedback loop is closed by an integral controller ‘observing’ the consumption and production of the cofactors by the glycolytic pathway and other cellular processes. We established a mixed layer structure for the glycolysis model (Fig 10). This mixed layer structure consists of the EFMs (not taking into account mass balance of the cofactors) as functionalities in a cascaded layer structure. The dynamics of cofactors, together with ATP consumption (ATP→ADP, Reaction 23 in [48]) and the AK reaction (ATP + AMP ↔ 2ADP, Reaction 24) form a control layer that communicates with all layers in the cascade without being part of the cascade. Using the software efmtool [50], we identified eight EFMs for the modified network where the cofactor dynamics were removed by setting the appropriate rows in the stoichiometric matrix to zero. Three of the identified EFMs were non-negative linear combinations of other EFMs, a consequence of treating each reversible reaction as two separate irreversible reactions. Based on the concept of simplicity (condition C3 in [13]), we removed the linearly-dependent EFMs with the lowest number of zero entries. The remaining five EFMs can be interpreted as representing elementary ‘tasks’ of the network: Glycogen buildup: production and storage of glucose-6-phospate G6P Production and excretion of glycerol Glyc Fermentation: production and excretion of ethanol EtOH Production and excretion of acetaldehyde ACA Lactonitrile lacto formation. Since the cofactors participate in many of the reactions, our approach to not take them into account when calculating the EFMs is conceptually similar to classifying metabolites taking part in more than a threshold number of reactions as ‘external’, as proposed in [51]. Both approaches result in the same significantly simpler and biologically interpretable set of EFMs. However, unlike the subnetworks identified in [48], the EFMs are not necessarily redox neutral. We will briefly discuss the consequences of this at the end of this section. In Fig 10 we represent the common reactions of the five EFMs by a graph. Surprisingly, the graph is a binary tree, indicating that, while there are junctions in the reaction pathway of the metabolites, there are no joins. Furthermore, there exists one main branch, consisting of Reactions 1 − 11 and Reaction 18, and each EFM can be represented by a unique junction from this main branch: EFM F1 is the junction after Reaction 3, F2 after Reactions 6 and 7, F3 after Reaction 11, and F4 and F5 after Reaction 18. This property provides a natural order for the EFMs, with F1 at the top of the cascade and F5 at the bottom. The ODEs describing the layers’ dynamics for this ordering of the functionalities, the dynamics of the control layer, the parameter values, and the initial conditions can be found in the Supporting Information. The layers and their interconnection can be represented as a species-reaction graph (SR graph, see [18]) depicting the reactions, the non-zero altered reactions, the species with non-zero dynamics in each layer, and the inter-layer communication (Fig 10). Recall from (7b) that the vector of altered reactions for layer i is defined by v a l t (x 1 + ⋯ + x i - 1, x i) ≔ v (x 1 + ⋯ + x i) - v (x 1 + ⋯ + x i - 1). In this example we identified reactions j where vj (x1 + … + xi) = vj (x1 + … + xi−1), meaning the rate of the reaction is not altered by the interconnection of functionality Fi. A sufficient condition [18] to conclude that altered reaction j does not change the dynamics of Layer i and thus can be omitted is if the jth element of I (δ δ x v) exp (I (∑ k = 1 i - 1 S k) I (δ δ x v) ) · I (S i) · 1 (13) is zero. Here, exp (.) represents the matrix exponent, Sk the stoichiometric matrix of layer k, I the element-wise indicator function being one if the respective element is non-zero and zero otherwise, and 1 the NR × 1 vector with all elements being one. Similarly, the integration of a species can be omitted in a layer if it is not affected by any of that layer’s reactions. Fig 10 shows that each layer is comparatively small and includes few altered reactions. Furthermore, since reactions R13 and R14 are linear (in the metabolites), information about the metabolite concentrations involved in these reactions is not required to be transmitted between the layers (for example, EtOH3 and EtOHX3 are not transmitted to Layer 4 and 5). This demonstrates that, for many biological networks, layers and the interfaces between them are comparatively simple. We equilibrated the model for a mixed flow glucose concentration of 10 mM and initialised the dynamics of the initial condition layer L (F0) = x0 to these steady-state values. This corresponds to low glucose conditions in the range of the Km values of high-affinity glucose transporters [52], and slightly below the mixed flow glucose concentration for which the parameters of the model were identified (18. 5 mM) to prevent glycogenic oscillations (see [48]). All other parameters and medium conditions were kept unchanged as compared to the version of the glycolysis model [48] available at the BioModels Database. We then calculated the layer dynamics for each ordering of the EFMs into all possible cascade structures by simulating the layering graph for 1000 min as described in Methods. This requires the simulation of 25 − 1 = 31 layers to populate the 5 ⋅ 5! = 600 nodes of the graph. Based on the layering graph we then analysed all pairwise interactions between Fi and Fj, given all possible combinations of other functionalities. In Fig 11, we display the steady-state incompatibilities and cooperativities between the EFMs after convergence at the end of the simulation. Note that although the incompatibilities and cooperativities converge, this is not necessarily true for the corresponding layer and mutual dynamics, which diverge for some EFMs (compare Fig 9). Several pairs have an incompatibility close to one and a cooperativity close to minus one. This combination is typical if one of the functionalities (given its environment) is unstable but is stabilised by the other. For example, in the isolated layer L (F3∣F0), pyruvate increases to infinity since the PDC reaction (Reaction 11 in [48]) is saturated and becomes rate limiting. Integrating either F1 or F2 with F3 reduces the production rate of pyruvate and thus stabilises the network, since both L (F1, F3∣F0) and L (F2, F3∣F0) are stable. On the other hand, in the isolated layer L (F5∣F0) both pyruvate and acetaldehyde (ACA) concentrations are unstable, the latter due to limited cyanide (CNX) availability. Integrating F1 or F2 with F5 only stabilises the pyruvate concentration, but not acetaldehyde. The prediction of infinite growth of species concentrations is unlikely to be observed experimentally. Instead, the instability might indicate that the conversion of pyruvate to acetaldehyde catalysed by the pyruvate decarboxylase (pdc) becomes a rate limiting step when synthetically increasing the flow through the respective EFM. Indeed, it was experimentally validated [53] that under certain conditions over-expression of pdc after reducing glycerol synthesis can lead to increased growth rates and ethanol yield. The incompatibilities between two EFMs, given the three other EFMs, are of specific interest (rightmost ten interactions in Fig 11). They indicate that F1 and F3 are significantly more incompatible than F2 and F3. Recall that F3 corresponds to the fermentation capability of the network. To increase biofuel production one would intuitively expect that the highest yields could be obtained by knocking out or down the EFMs which are most incompatible to F3, thereby removing EFMs which act to attenuate its function. Recall that the effect of knocking-out a certain functionality Fi on the overall dynamics of the network is equal to −L (Fi∣Fj, j ≠ i). Thus, the effect of knocking out GPD (Reaction 15) and, thus, glycerol production F2 is −L (F2∣F0, F1, F3, F4, F5). The effect of a double knock-out of glycerol production together with glycogen build-up F1 corresponds to −L (F1, F2∣F0, F3, F4, F5). The theoretical effect on fermentation efficiency of all possible combinations of knock-outs (Fig 12) can thus be directly derived, given the layering graph. Indeed, this analysis confirms that fermentation F3 is significantly more incompatible with glycogen build-up F1 than with glycerol production F2 at the given experimental setup. The effect of knocking out more than one functionality is non-linear, and a double knock-out of the two functionalities yielding the highest gains when knocked out separately is not necessarily optimal, as shown in Fig 12. Thus, the layering graph provides an effective tool to analyse such potential interactions. It is out of the scope of this article to take side effects of the proposed genetic modifications into account, such as possible viability issues when knocking out some of the pathways. Our intention was to provide a proof of principle of how to apply our layering framework in the context of metabolic engineering. Notably, our analysis is based on a model of a non-industrial strain of S. cerevisia cells grown under glucose starvation, and which was created mainly to explain glycogenic oscillations at low glucose conditions [48] rather than being optimised for metabolic engineering purposes. Furthermore, we essentially avoided the challenge of balancing the redox state (NAD+/NADH ratio) by our mixed layering approach (Fig 10) in which the control layer was always receiving the cofactor consumption and production rates of all cascaded layers while simulating the layering graph. Nevertheless, our analysis indicates that ethanol yield can be increased by approximately 6% by knocking out glycerol formation (F2). This prediction is in good agreement with a reported [54] experimental increase of about 7% after inhibition of glycerol formation from dihydroxyacetone phosphate (DHAP) by knock-down of glycerol-3-phosphate dehydrogenase (GPD). In [54] the redox state was balanced by synthetically engineering an alternative route in the endogeneous Embden–Meyerhof–Parnas pathway based on NADP+ and NADPH instead of NAD+ and NADH.
We have described a cascaded layering approach which can be used to systematically identify the interactions of functionalities (i. e. functional subsystems) in a decomposed network. Functionalities are defined as sets of reactions which together accomplish a given purpose. We have identified the dynamics of an isolated functionality as the dynamics of the subnetwork made up of only those reactions. However, the dynamics of a functionality are different when it is in the context of others. We have formulated the conditional dynamics of a functionality as the incremental change to the dynamics caused by integrating it with its context. We have also demonstrated that the computational burden associated with this layered framework can be minimised by exploiting symmetries within the definitions of interdependence. We have used the conditional dynamics to define the mutual dynamics, which describe the context-dependent interaction between any pair of functionalities. We can thus identify if the interaction between them is strong or weak (incompatibility), and to what extent they amplify or attenuate each other (cooperativity). This interdependence is also context-dependent, so that two functionalities may interact in vastly different ways depending on the wider system in which they are integrated. Our framework allows the unambiguous quantification of nonlinear interactions between functionalities. Finally, we illustrated our layering framework with three examples. The first considered signalling pathways interconnected by two distinct crosstalk mechanisms, a mutually inhibiting and a mutually excitatory crosstalk. We not only quantified how each mechanism separately influences the dynamics, but also how the crosstalk mechanisms interact with each other. Our second example concerned a pathway stabilised by two integral feedback mechanisms. For certain input signals, the pathway alone becomes unstable. It gets stabilised by the first integral feedback, but the second feedback further destabilises the system, as displayed by the negative and positive cooperativity respectively. Thus each feedback has an opposite effect on the pathway to the other when each is considered in isolation. However, when they are integrated together they interact cooperatively, due to the first controller increasing its strength to stabilise the network when the second controller is present. Third, we showed how to apply our layering framework in the context of metabolic engineering to analyse the dynamic interactions between elementary flux modes in glycolysis. Interestingly, the effect of abolishing an EFM by knocking out its enzymes can be easily described as the negative of the conditional dynamics of that EFM, given the rest of the network. The incompatibility between two EFMs, on the other hand, gives information about the expected increase in the flux through one EFM if the other one is knocked out. These interpretations allowed us to efficiently calculate the expected increase of ethanol yield in biofuel production when knocking out arbitrary combinations of EFMs. Interestingly, this yield is not additive, and a double knock-out of the two EFMs giving the highest increases in ethanol yield alone is not necessarily optimal. Instead, our method to calculate the integrated effect of all combinations of knock-outs allowed us to assess the best metabolic engineering strategy with a minimum amount of simulations required. As discussed in the Introduction, the layering framework [16,17] differs from modularization approaches [18–26]. Each of these frameworks allows different network architectures to be considered, reflecting a distinction between vertical and horizontal network decomposition [55]. Layered networks represent the overlay of multiple (possibly competing) functional subnetworks, while the modular framework interprets biological networks in terms of engineered interconnections of input–output systems. Although it is a modular phenomenon, previous work on quantifying retroactivity [4,5] can be related to our layered formulation. A goal of both techniques is to measure the difference in a given subsystem’s isolated and integrated behaviours. In particular, an upper bound on the Euclidean norm of the difference in a module’s isolated and integrated trajectories has been derived in terms of the system parameters in [5]. In the Supporting Information to this paper, we show in detail how this norm relates to our formulation of the incremental dynamics L (F2∣F1), and the mutual dynamics M (F1, F2). Although our layering framework and the concept of retroactivity were developed to quantify the interactions between different kinds of subsystems, this comparison suggests that our additional measures of incompatibility and cooperativity as defining an interaction direction may also be applicable to further understand retroactivity in multi-dimensional modular systems. The computational method employed to quantify the interaction strength between a given functionality pair is clearly dependent on the kinetic constants and other parameter values in the model, including the initial conditions. A repetition of this computation across a range of such values may give a more informative picture of how uncertain biological functionalities interact, but will be difficult to visualise and expensive to evaluate. Possible directions for future work may include the derivation of analytical estimates in terms of the system’s parameter values for the mutual dynamics, incompatibility, or cooperativity, of the type produced for layered steady state perturbations in [16] or retroactive perturbations in [5]. It would also be interesting to consider if any estimates of our measures can be computationally derived to hold for all parameters, similarly to the structural results on the direction of steady state responses to parameter perturbations given in [34]. Alternatively, semi-definite programming may be used to calculate worst-case estimates of the difference between the trajectories L (F1) and L (F1∣F2) without resorting to simulation, similarly to previous work on model reduction error estimation [56]. An important assumption made in this paper was that the network functionalities Fi are given. The examples in this paper provide two different strategies for how the functionalities necessary for our layering framework can be defined. In the first two examples this was done by intuition and prior knowledge, while in the third example we applied an already available method to group functionally related reactions in metabolic pathways, namely EFMs. Clearly, this is not a comprehensive list of how functionalities can be defined, and (depending on the precise application) other strategies might be more promising. One such example may arise in the context of Synthetic Biology, from the specification of synthetic biomolecular devices such as toggle switches, oscillators, and so on, the dynamics of which are significantly affected when combined into large-scale systems [2]. Such an application of our framework to Synthetic Biology would require these biomolecular devices–typically designed and modelled as modules–to be mathematically expressed as layers, respectively functionalities. If this is possible without requiring modifications in the biomolecular implementation remains a question for future research. Our layering approach provides a general and concise framework for the quantification of the nonlinear interactions between functionalities of signalling and metabolic pathways. It offers a great flexibility since it only requires that each reaction must be part of a functionality, but allows reactions to be part of more or even all functionalities. Besides the mathematical definition of the model, the only input data needed is the problem-specific definition of the sets of reactions which make up a functionality. For many signalling pathways these definitions are typically given by biological insight, and for metabolic networks implementations of efficient algorithms to determine the sets of EFMs are available (e. g. [50]). To allow a quick assessment of our approach, we have implemented a MATLAB (The MathWorks, Natick, MA) toolbox, available under the GNU General Public License from sysos. eng. ox. ac. uk/control/sysos/index. php/User: Prescott/Code. This toolbox provides algorithms to automatically construct the models of the layers from SBML files, to derive the layering graph with minimal numerical simulation, and to easily assess the mutual dynamics, the incompatibility and the cooperativity between functionalities. We hope that our implementation can serve as a starting point to integrate our layering framework into standard software solutions for the analysis of biomolecular networks.
To better understand the dynamic behaviour of cells and their interaction with the environment, mathematical models describing the interplay between proteins, metabolites or signalling molecules are used extensively in Systems Biology. Typically, such models focus on single functional subsystems and neglect the rest of the biochemical reaction network. However, the behaviour of multiple functional subsystems when integrated together can differ significantly from each subsystem’s isolated behaviour. In this article we describe a methodology for assessing the nonlinear effects of combining multiple functional subsystems of a biological system. This is key for answering questions related to Systems and Synthetic Biology as well as Metabolic Engineering. For example, if we can identify the isolated behaviours of two subsystems, we can determine if they persist when the subsystems interact. Similarly, we can show how modifications to single functional subsystems (such as increasing particular metabolic yields) have different effects in the context of the integrated system.
Quantification of Interactions between Dynamic Cellular Network Functionalities by Cascaded Layering