Michaelis-Menten dynamics in protein subnetworks

# Michaelis-Menten dynamics in protein subnetworks

## Abstract

To understand the behaviour of complex systems it is often necessary to use models that describe the dynamics of subnetworks. It has previously been established using projection methods that such subnetwork dynamics generically involves memory of the past, and that the memory functions can be calculated explicitly for biochemical reaction networks made up of unary and binary reactions. However, many established network models involve also Michaelis-Menten kinetics, to describe e.g. enzymatic reactions. We show that the projection approach to subnetwork dynamics can be extended to such networks, thus significantly broadening its range of applicability. To derive the extension we construct a larger network that represents enzymes and enzyme complexes explicitly, obtain the projected equations, and finally take the limit of fast enzyme reactions that gives back Michaelis-Menten kinetics. The crucial point is that this limit can be taken in closed form. The outcome is a simple procedure that allows one to obtain a description of subnetwork dynamics, including memory functions, starting directly from any given network of unary, binary and Michaelis-Menten reactions. Numerical tests show that this closed form enzyme elimination gives a much more accurate description of the subnetwork dynamics than the simpler method that represents enzymes explicitly, and is also more efficient computationally.

## I Introduction

Biological networks are often complex and it is frequently necessary to focus on subnetworks of a larger system, e.g. to enable a better understanding of the network properties Bhalla (2003); Ackermann et al. (2012); Conradi et al. (2007). Such subnetworks may be of interest because they carry out important biological functions or because they capture parts of the system where there is less uncertainty in the network structure or dynamical parameters. The choice of subnetwork may also be dictated by experimental data only being available for a limited number of molecular species.

There are many different methods of model reduction that have been used to simplify a large model down to one for a subnetwork Okino and Mavrovouniotis (1998); Radulescu et al. (2012). We have previously used projection techniques to find a systematic description of the dynamics of a subnetwork embedded in a “bulk” network Rubin et al. (2014), identifying the occurrence of memory terms as one of the key features. The resulting methods for calculating memory functions were, however, restricted to networks involving only unary and binary reactions with concentration-independent reaction rates. While this class of networks is large, it excludes networks with reactions following Michaelis-Menten kinetics, which is commonly used to represent e.g. enzyme reactions. Our aim in this paper is to remove this restriction, by providing an explicit method for constructing the description of the dynamics of any subnetwork within a network consisting of unary, binary and Michaelis-Menten reactions.

An approximate way of achieving the above aim is to extend a network with Michaelis-Menten reactions to one that explicitly involves enzymes and enzyme complexes. The Michaelis-Menten terms in the kinetics are thus “unfolded” into unary and binary reactions so that the original projection approach Rubin et al. (2014) can be applied. However, this is inconvenient both conceptually – we need to include extra species not present in the original system of reaction equations – and numerically, because the fast rates of enzyme reactions typically create a “stiff” system of differential equations that has to be integrated using small timesteps. The approach nevertheless provides the inspiration for the method we adopt here: we follow the above route analytically, taking the limit of fast enzyme reaction rates in such a way as to retrieve the original Michaelis-Menten kinetics exactly. The main achievement of our analysis is to show that the limit can be taken in closed form. This leads to a procedure for constructing the subnetwork dynamics directly from the original reaction network, without reference to any extra species.

In Section II.1 we give a summary of Michaelis-Menten dynamics and the conditions under which it is retrieved as the limit of a fast formation and dissociation of an enzyme complex; in essence these conditions are that the enzyme reactions must be fast, and the enzyme concentration low. Section III details our approach to obtaining the projected equations that describe subnetwork dynamics, for systems of reaction equations that include Michaelis-Menten reactions: we temporarily add enzymes to represent these reactions, derive the projected equations, and then take the limit of fast enzyme rate and low enzyme concentration in closed form. We describe the approach separately for linearised (Sec. IV) and nonlinear (Sec. V) dynamics, as the nonlinear case is more complicated technically but follows the same conceptual route. Remarkably, even though Michaelis-Menten terms are generally nonlinear, we find that in the memory terms that are characteristic of the projected equations no additional nonlinearities appear, i.e. the memory terms involve linear concentration fluctuations for linearised dynamics, and linear and quadratic concentration fluctuations for nonlinear dynamics. Finally in Section VI we compare predictions from the original reaction equations with Michaelis-Menten dynamics to the projected equations with either enzymes added explicitly or eliminated in closed form using the method derived in this paper. We find that the closed form elimination is both faster to evaluate computationally, and gives a more accurate approximation to the original reaction equations.

## Ii Michaelis-Menten dynamics

Many biochemical reactions are catalysed by enzymes. Generally each enzyme will enable a particular reaction without being consumed. A simple model of an enzyme reaction Henri (1902); Michaelis and Menten (1913) is written

 u+ek+ue,c⇌k−c,ueck−c,pe−−→e+p (1)

where is a substrate (we use “u” not “s” here as we will need the letter “s” to denote “subnetwork” later), is an enzyme, is an enzyme-substrate complex and is a product. The reaction rates are denoted . These reactions describe the binding of a free enzyme with a substrate to form a substrate-enzyme complex. This complex can then dissociate into enzyme and a product. In the traditional model substrate binding is reversible but product formation is not; we will consider also the more general case below, where both reactions are reversible.

### ii.1 Derivation of Michaelis-Menten equations

Let be the concentration of species . Then the set of mass action equations for system (1) is

 ∂∂txu=k−c,uexc−k+ue,cxuxe∂∂txe=k−c,uexc−k+ue,cxuxe+k−c,pexc∂∂txc=−k−c,uexc+k+ue,cxuxe−k−c,pexc∂∂txp=k−c,pexc (2)

From these equations it follows that there is a conservation law between the enzyme and complex such that

 ∂∂txe+∂∂txc=0⟹xe+xc=xtote (3)

The Michaelis-Menten description of the dynamics is obtained by exploiting the fact that enzyme reactions are typically fast. This allows one to reduce the system (1) to a simpler description where the enzyme and enzyme complex no longer appear explicitly.

To achieve this simplification one uses the quasi-steady state assumption Briggs and Haldane (1925). We assume all enzyme reactions are fast so that the enzyme complex is in quasi-steady state at any time. We use the qualifier “quasi” because this steady state depends on the concentrations of both substrate and product, which themselves generally vary in time. For the reaction flux one then finds

 v=k−c,pextotexu(k−c,ue+k−c,pe)/k+ue,c+xu=VmaxxuKm+xu (4)

where

 Vmax=k−c,pextote (5)

is the maximum reaction flux that can be achieved and

 Km=k−c,ue+k−c,pek+ue,c (6)

is the Michaelis constant. The simplified description of the original reaction system can now be written in terms of the reaction flux in (4), as

 −∂∂txu=VmaxxuKm+xu=∂∂txp (7)

In the limit where the complex dissociates into enzyme and substrate at a much higher rate than for enzyme and product, one can achieve a further simplification known as the “rapid equilibrium assumption” Michaelis and Menten (1913). However, this is just a limiting case of the quasi-steady state assumption where is neglected against in determining the Michaelis constant .

### ii.2 Reversible Michaelis-Menten dynamics

We discuss briefly how the above analysis is modified when there is a back reaction from the enzyme and product to the complex Haldane (1930); Sauro (2013). In such cases one has the more general reaction scheme

 u+ek+ue,c⇌k−c,ueck−c,pe⇌k+pe,ce+p (8)

with the additional rate constant . The equation for the enzyme complex concentration now has an extra contribution , while the equations for enzyme and product contain the same additional term with a negative sign. Using the quasi-steady state assumption one then obtains a reaction flux of the form

 v=Vuxu/Ku−Vpxp/Kp1+xu/Ku+xp/Kp (9)

where and are the maximum reaction rates for the forward and reverse reactions respectively and are given by

 Vu=k−c,pextote,Vp=k−c,uextote (10)

Similarly and are the Michaelis constants for the forward and reverse reactions and are written in terms of the mass action rate constants as

 Ku=k−c,ue+k−c,pek+ue,c,Kp=k−c,ue+k−c,pek+pe,c (11)

### ii.3 Quantitative accuracy of Michaelis-Menten approximation

We outlined above how the Michaelis-Menten dynamics (7) – or its generalisation (9) – emerges from a quasi-steady state approximation. We now consider under what conditions on the enzyme reaction parameters these approximations become exact, so that we can later take the appropriate limit for the enzyme parameters in the construction of the projected equations. Intuitively the enzyme reactions need to be fast; from the definition of the maximal reaction rate (5), which in a reaction specified from the outset in Michaelis-Menten form is given, the total enzyme concentration then needs to be small. This makes sense as the Michaelis-Menten kinetics always gives equal and opposite reaction fluxes for substrate and product, neglecting any substrate or product captured in the enzyme complex. A more formal analysis by Segel and Slemrod (1989), based on a singular perturbation approach to a dimensionless form of the reaction equations, gives as the condition for validity of the quasi-steady state approximation. As can also be regarded as fixed by the specification of a Michaelis-Menten reaction, this is consistent with the intuitive requirement of small enzyme concentration.

The discussion so far suggests that we should rewrite the enzyme reaction rate constants and total enzyme concentrations as

 k+ue,c=¯k+ue,cγ, k−c,ue=¯k−c,ueγ, k−c,pe=¯k−c,peγ, xtote=¯xtote/γ (12)

where is a dimensionless fast rate scale parameter. The definitions of (5) and (6) then allow us to write the remaining parameters as

 ¯xtote=Vmax/¯k−c,pe, ¯k+ue,c=¯k−c,ue+¯k−c,peKm (13)

We now have three parameters, , and , that parameterise the possible enzyme kinetics underlying a given Michaelis-Menten reaction. From the criterion of Segel and Slemrod (1989) the Michaelis-Menten description will then become exact for , irrespective of the values of and Van Slyke and Cullen (1914). The same exactness statement applies to the more general case where there is also a back reaction from product and enzyme to form an enzyme complex Li and Li (2013); Kollar and Siskova (2015).

## Iii Enzyme reactions in the projected equations

The projection method as applied to protein interaction networks Rubin et al. (2014) works with mass action kinetics for unary and binary reactions. Therefore if we are given an interaction network that includes Michaelis-Menten reactions then a priori we need to represent these explicitly in such a form, to allow us to compute the projected equations. As explained in the introduction, this is a disadvantage both conceptually and computationally. Our aim here is to implement this approach analytically instead: we add enzyme species to the network to get mass action kinetics and take the large enzyme rate and low enzyme concentration limit in which the mass action dynamics becomes exactly identical to Michaelis-Menten dynamics. The challenge is to understand what happens in this limit to the projected equations.

The aim of the projection method generally is to provide a description of the dynamics of a protein interaction subnetwork embedded in a bulk network. The Zwanzig-Mori projection method Mori (1965) can be used to obtain such a description, specifically equations for the time evolution of the protein concentrations in a chosen subnetwork. Full details of the projection method applied to protein interaction networks are given in Rubin et al. (2014). We summarise below the features necessary for the analysis of enzyme dynamics.

The natural specification of the state of a given network is the vector of concentrations for all molecular species in the network, which we call . In the steady state will fluctuate around some mean and it will be useful to switch variables to , which has zero mean in steady state. The protein concentrations (concentration deviations from the mean, more precisely) evolve in time according to a Fokker-Planck equation, where the stochasticity arises from copy number fluctuations. The variance of these fluctuations scales as the inverse reaction (e.g. cell) volume . While a nonzero is needed initially to apply the Zwanzig-Mori formalism, we consider the limit of small throughout.

The time evolution of any observable is given by , where is the adjoint Fokker-Planck operator whose drift term encodes the mass action kinetics. We showed in Rubin et al. (2014) that if we focus on a set of observables containing and all products like and , then the operator can be written in matrix form such that

 ∂∂tzα=∑βzβLβα+δx3+O(ϵ) (14)

where indicates terms of third or higher order in while the terms vanish in the low noise limit.

If we now let be a set of observables from the subnetwork, the projected equations have the form

 ∂∂taα(t)=∑βaβ(t)Ωβα+∫t0dt′∑βaβ(t′)Mβα(t−t′)+rα(t) (15)

We choose specifically as subnetwork observables all the subnetwork concentrations and their products, i.e. the entries of only involving the subnetwork. We denote these collectively by S, and the remaining species – all bulk concentrations (b) and concentration products involving the bulk (sb and bb) – by the letter B.

With the observables ordered appropriately, we can then decompose the matrix formed from the entries into subnetwork and bulk blocks:

 L=(LS,SLS,BLB,SLB,B) (16)

The four blocks here can be broken down further according to the specific type of subnetwork and bulk observable, for example

 LB,S=⎛⎜ ⎜⎝Lb,s0Lsb,sLsb,ssLbb,sLbb,ss⎞⎟ ⎟⎠ (17)

Here “s” denotes linear subnetwork observables, “ss” product observables from the subnetwork, “sb” cross-product observables between subnetwork and bulk, and so on. In this way contains the linear coefficients of bulk concentrations in the equations of motion for subnetwork concentrations, while the coefficients of subnetwork-bulk products in these equations are in the block .

In terms of the blocks of the -matrix defined as above, the “rate matrix” in (15) is simply Rubin et al. (2014)

 Ω=LS,S=(Ls,s0Lss,sLss,ss) (18)

It contains terms from the subnetwork dynamics that are local in time. The memory function (matrix) can be written as Rubin et al. (2014)

 M(Δt)=LS,BeLB,BΔtLB,S (19)

where . The entries of this matrix, which appear in (15), are the memory functions : they determine how strongly the past values of the observable affect the present rate of change of . For later use we note here the Laplace transform version of the memory function, which is

 ^M(z)=LS,B(z−LB,B)−1LB,S (20)

For brevity we will use the term “memory function” also for when it is clear from the context that the Laplace transform is meant.

One important property of the memory functions is their boundary structure: if we define a boundary species as a subnetwork species that is involved in a reaction with a bulk species, then among the projected equations for subnetwork concentrations only those for boundary species contain memory terms. We also note that linear memory functions are only nonzero for boundary species influencing other boundary species. Similarly, only concentration products involving at least one boundary species appear in nonlinear memory terms Rubin et al. (2014).

The final term in (15), , is what is known as the random force. It accounts for the fact that because of the interaction between subnetwork and bulk, the time evolution of the subnetwork observables cannot be closed. When using the projected equations to make predictions in numerical examples, we will drop the random force as there is no simple way of calculating it. The otherwise exact projected equations thus become an approximation, but one that in previous work Rubin et al. (2014) we showed to be remarkably accurate.

Our analysis starts from a given reaction network involving unary and binary protein reactions, and enzyme reactions described by Michaelis-Menten terms. We assume for simplicity that all enzyme reactions are reversible; the irreversible scenario can be obtained from this as the limiting case where the rate of dissociation into enzyme and product is much larger than the rate for formation of enzyme complex in the reverse direction. The mass-action kinetics for each enzyme reaction is

 ∂∂tδxu=−fue,c+…∂∂tδxe=−fue,c−fpe,c∂∂tδxp=−fpe,c+… (21)

where the dots indicate fluxes from other reactions. The reaction fluxes from substrate and enzyme to complex, and from product and enzyme to complex, read respectively (compare (2))

 fue,c=k−c,ue(−xtote+ye+δxe)+k+ue,c(yu+δxu)(ye+δxe)fpe,c=k−c,pe(−xtote+ye+δxe)+k+pe,c(yp+δxp)(ye+δxe) (22)

Here we have used the enzyme conservation law (3) to eliminate the complex concentration via . The enzyme steady state concentration can be found by requiring that in the steady state, where , the two fluxes must sum to zero to ensure . This gives

 ye=xtotek−c,ue+k−c,pek−c,ue+k−c,pe+k+ue,cyu+k+pe,cyp (23)

We now write the enzymatic reaction rate constants in terms of a fast rate scale as in (12), and similarly the steady state enzyme concentration, which from (23) must scale as the inverse of , i.e. , like the total enzyme concentration . This gives

 fue,c=¯k−c,ueγ(−¯xtote/γ+¯ye/γ+δxe)+¯k+ue,cγ(yu+δxu)(¯ye/γ+δxe)fpe,c=¯k−c,peγ(−¯xtote/γ+¯ye/γ+δxe)+¯k+pe,cγ(yp+δxp)(¯ye/γ+δxe) (24)

The scaled rate constants and steady state enzyme concentration are related to the Michaelis-Menten parameters as explained after (8), i.e.

 Vu=¯k−c,pe¯ye(¯k−c,ue+¯k−c,pe+¯k+ue,cyu+¯k+pe,cyp)¯k−c,ue+¯k−c,pe,Ku=¯k−c,pe+¯k−c,ue¯k+ue,cVp=¯k−c,ue¯ye(¯k−c,ue+¯k−c,pe+¯k+ue,cyu+¯k+pe,cyp)¯k−c,ue+¯k−c,pe,Kp=¯k−c,pe+¯k−c,ue¯k+pe,c (25)

In the above representation it is not obvious which terms have to be regarded as “fast” in the remainder of the analysis, and which as slow. We therefore switch to dimensionless concentration variables . In terms of these we have

 ∂∂tδ~xu=−fue,c/yu+…∂∂tδ~xe=−γ(fue,c+fpe,c)/¯ye∂∂tδ~xp=−fpe,c/yp+… (26)

with enzymatic fluxes

 fue,c=¯k−c,ue¯ye(−¯xtote/¯ye+1+δ~xe)+¯k+ue,cyu¯ye(1+δ~xu)(1+δ~xe)fpe,c=¯k−c,pe¯ye(−¯xtote/¯ye+1+δ~xe)+¯k+pe,cyp¯ye(1+δ~xp)(1+δ~xe) (27)

Here one sees clearly that the enzyme evolution equation contains only fast terms that scale with , while the equations for substrate and product only contain slow terms. We will therefore use dimensionless concentrations throughout, and to lighten the notation we will in the following drop the tildes, as well as the bars indicating rescaling with . Note also that as at steady state () the total flux into or out of any molecular species must be zero, we can drop the constant pieces from and : they have to cancel against each other in the equation for , and against the other steady state fluxes in the equations for and .

Our task now is to take a mass-action (unary and binary) reaction system where every enzyme reaction is represented as in (27), and to find closed form expressions for the rate matrix and memory functions in the limit where we know that this mass-action description becomes identical to Michaelis-Menten dynamics. The effect of the enzyme reactions depends on where they are located relative to subnetwork and bulk, with three potentially distinct cases as shown in Fig. 1. If both substrate and product are in the subnetwork, also the enzyme will be located there and, as is clear from Fig. 1a, away from the boundary of the subnetwork. We will therefore find “fast” equations of motion for such enzymes without any memory terms. These enzymes can therefore be kept explicitly in a first stage of our analysis, and eliminated in a second stage following the standard logic that leads to the Michaelis-Menten description.

If both substrate and product are within the bulk then the entire enzyme reaction takes place there (Fig. 1b), contributing fast reaction rate constants to . Accordingly we will find that such reactions only give contributions to the memory functions, not the rate matrix.

The third case is the one where the substrate is on the boundary and the product is in the bulk (or vice versa, but for reversible enzyme reactions the two cases are mathematically equivalent). We then choose to assign also the enzyme and the enzyme complex to the bulk, and refer to this situation as a (bulk) enzyme on the boundary. Reactions involving such enzymes, whose rate constants appear in and , will contribute fast terms in the memory functions that decay on a timescale of order . In the limit these terms become local in time and so turn into contributions to the rate matrix.

## Iv Linearised Dynamics

In linearised dynamics we only consider terms in the mass-action kinetics up to linear order in . The dimensionless scaled reaction equations for a reversible Michaelis-Menten reaction are then, from (26) and (27),

 ∂∂tδxu=−k−c,ue(ye/yu)δxe−k+ue,cye(δxu+δxe)+…∂∂tδxe=−γ[(k−c,ue+k−c,ep)δxe+k+ue,cyu(δxu+δxe)+k+ep,cyp(δxp+δxe)]∂∂tδxp=−k−c,ep(ye/yp)δxe−k+ep,cye(δxp+δxe)+… (28)

Let us partition the matrix form of the adjoint Fokker-Planck matrix operator (16) so that the bulk species are split into fast and slow blocks. If e and e represent the collection of subnetwork and bulk enzymes respectively and s and b represent the other molecular species in the subnetwork and bulk, we partition as

 L=(LS,SLS,BLB,SLB,B)=⎛⎜ ⎜ ⎜ ⎜ ⎜ ⎜⎝Ls,sLs,e′Ls,bLs,eLe′,sLe′,e′Le′,bLe′,eLb,sLb,e′Lb,bLb,eLe,sLe,e′Le,bLe,e⎞⎟ ⎟ ⎟ ⎟ ⎟ ⎟⎠=⎛⎜ ⎜⎝c|ccmw1f1w2w3f2w4w5f3⎞⎟ ⎟⎠ (29)

where are slow terms and are fast terms; the top left block denoted contains a mixture of fast and slow terms. In writing the last equality above we have grouped s and e together; the resulting specific block structure of slow and fast terms is one that we will find again in the case of the full nonlinear dynamics. Note that because subnetwork enzymes only have interactions with subnetwork species (s and e), and are zero. Similarly, because bulk proteins or enzymes do not interact with subnetwork enzymes, and vanish. This means that subnetwork enzymes do not feature at all in the calculation of the memory function (20), which makes intuitive sense. The vanishing of and is important also as these blocks contain rate constants for the time evolution of (subnetwork) enzymes, which by our construction scale with : if these blocks were nonzero, it would change the character of and from slow to fast.

To analyse the memory function (19) that results from (29), we note that has both slow and fast sub-blocks. As a result the memory function should in general have both slow contributions that decay on timescales, and fast contributions that decay for time differences of . As the memory function appears as a weight in an integral over the past (15), the fast contributions only matter for if their amplitude is proportional to so that the integral over all time differences remains finite. Accounting also for subleading terms in the amplitude dependence then suggests the following decomposition of the memory function:

 M(Δt)=(γM0f(Δ¯t)+M1f(Δ¯t)+…)+(M0w(Δt)+1γM1w(Δt)+…) (30)

where . In principle an arbitrary constant can be added to e.g.  and subtracted from ; we fix this constant by requiring that all fast contributions decay to zero for large . Fig. 2 shows an example of the above decomposition.

The leading fast and slow contributions can now be extracted relatively simply from the Laplace transform of the memory function (20), which has the decomposition

 ^M(z)=(^M0f(¯z)+1γ^M1f(¯z)+…)+(^M0w(z)+1γ^M1w(z)+…) (31)

where . The leading fast term can be obtained by taking the limit

 limγ→∞^M(z)∣∣¯z=const=^M0f(¯z) (32)

because the subleading fast terms are down by powers of , and in the slow terms so that the Laplace transforms etc. vanish. The leading slow part can then be found from a limit at constant , namely

 limγ→∞[^M(z)−^M0f(0)]∣∣z=const=^M0w(z) (33)

We could have equivalently written inside the square brackets, as when at fixed .

Using the method above, we can now find the fast and slow pieces of the memory function derived from the adjoint Fokker-Planck (matrix) operator in equation (29). From (20) this memory function reads

 ^M(z)=(w1f1)(z−w3−f2−w5z−f3)−1(w2w4) (34)

and after working out the inverse using standard block inversion identities and simplifying we obtain

 ^M(z)=(w1−f1(z−f3)−1(−w5))(z−w3+f2(z−f3)−1w5)−1(w2+f2(z−f3)−1w4)+f1(z−f3)−1w4 (35)

If we now write all fast blocks as , we can see that application of (32) identifies the fast part of the memory function as

 ^M0f(¯z)=¯f1(¯z−¯f3)−1w4 (36)

while (33) gives for the slow part

 ^M0w(z)=(w1−¯f1(¯f3)−1w5)(z−w3−¯f2(¯f3)−1w5)−1(w2−¯f2(¯f3)−1w4) (37)

Here we have used the fact that in the combination , the first term can be neglected when at constant .

The fast part of the memory decays on an ever shorter timescale as increases. In the limit, and when it is used inside a memory function integral, it becomes equivalent to a delta function multiplied by the area under the fast piece, which is just . The rate matrix that we obtain from (29) for is therefore

 Ω=m−¯f1(¯f3)−1w4 (38)

while the memory function in the same limit is given by (37). We can now compare to the rate matrix and memory function that would result from an adjoint Fokker-Planck operator like (29) without the fast bulk variables e, i.e.

 L∖e=(mw1w2w3) (39)

This would give and . Looking at (37) and (38), we conclude that the large- limit gives results that can equivalently be obtained by using an adjoint Fokker-Planck matrix without the fast bulk variables that is modified from to , where

 ΔL∖e=−(¯f1(¯f3)−1w4¯f1(¯f3)−1w5¯f2(¯f3)−1w4¯f2(¯f3)−1w5) (40)

This is the key result of our first stage of elimination, which has removed the fast bulk variables.

### iv.1 Bulk enzyme elimination as quasi-steady state method

Before moving on to the second stage of also eliminating the subnetwork enzymes, we pause briefly to give a simpler form of our last result. While (40) gives a closed form for the effective -matrix we obtain after eliminating the bulk enzymes, this form is not very intuitive. We show next that there is a much simpler statement of the result, namely that can be obtained by treating all bulk enzymes as in quasi-steady state with the other molecular species.

To see this, we reinstate in the generic block structure of (29) the specific notation for the linearised dynamics considered here, i.e.

 L=⎛⎜ ⎜⎝mw1f1w2w3f2w4w5f3⎞⎟ ⎟⎠=⎛⎜ ⎜⎝LS,SLS,bLS,eLb,SLb,bLb,eLe,SLe,bLe,e⎞⎟ ⎟⎠ (41)

where, as before, S collects all subnetwork variables, i.e. subnetwork proteins s and subnetwork enzymes e. The dynamics of the system can then be written as

 ∂∂tδxST=δxSTLS,S+δxbTLb,S+δxeTLe,S∂∂tδxbT=δxSTLS,b+δxbTLb,b+δxe% TLe,b∂∂tδxeT=δxSTLS,e+δxbTLb,e+δxe% TLe,e (42)

If we now impose a quasi-steady state condition for the bulk enzymes , this gives

 δxeT=−(δxSTLS,e+δxbTLb,e)(Le,e)−1 (43)

Substituting this back into the equations of motion for the subnetwork species and the bulk proteins gives

 ∂∂tδxST=δxST(LS,S−LS,e(Le,e)−1Le,S)+δxbT(Lb,S−Lb,e(Le,e)−1Le,S)∂∂tδxbT=δxST(LS,b−LS,e(Le,e)−1Le,b)+δxbT(Lb,b−Lb,e(Le,e)−1Le,b) (44)

This is exactly the dynamics that is defined by the effective -matrix derived above (see (39) and (40)), hence proving our claim that this matrix can be constructed by imposing a quasi-steady state condition for the bulk enzymes.

### iv.2 Michaelis-Menten terms as effective unary reactions

The above bulk enzyme elimination can be carried out in closed form. Each bulk enzyme can be eliminated by setting the time derivative of its concentration to zero. Here and throughout we assume that each enzyme only catalyses one reaction. (In the matrix formulation, this implies that the block is diagonal.) We then find for the effective dynamical equations of any substrate and product the following form:

 ∂∂tδxu=−λupδxu+λpu(yp/yu)δxp+…∂∂tδxp=−λpuδxp+λup(yu/yp)δxu+… (45)

where

 λup=k+ue,cye(k−c,pe+k+pe,cyp)k−c,ue+k−c,pe+k+ue,cyu+k+pe,cypλpu=k+pe,cye(k−c,ue+k+ue,cyu)k−c,ue+k−c,pe+k+ue,cyu+k+pe,cyp (46)

are the rates for effective unary reactions converting substrate to product and back, respectively. The factors of in (45) arise because we are using dimensionless concentration variables.

The bulk enzyme elimination thus has the simple effect of replacing all bulk Michaelis-Menten reactions by unary conversion reactions with constant rates. From (25) one sees that these effective rates can be expressed directly in terms of the Michaelis-Menten parameters, as

 λup=Vu/Ku+(Vu+Vp)/Ku(yp/Kp)(1+yu/Ku+yp/Kp)2λpu=Vp/Kp+(Vu+Vp)/Kp(yu/Ku)(1+yu/Ku+yp/Kp)2 (47)

Comparing with (9) shows that the rates are obtained by linearising the Michaelis-Menten reaction flux around the steady state concentrations of substrate and product. This is the closed-form procedure for bulk enzyme elimination we were after: it only requires as input the Michaelis-Menten parameters of the original network, and its steady state.

Note that the above discussion includes enzyme reactions both entirely in the bulk, or on the boundary of subnetwork and bulk (cf. Fig. 1b and c). The only difference between these two cases is that for the latter group of enzymes, the effective unary reactions we have derived are between a subnetwork boundary species and a bulk species and so will contribute to the rate matrix, while for enzyme reactions entirely in the bulk the effective reactions only affect the memory function.

### iv.3 Elimination of subnetwork enzymes

So far we have described how bulk enzymes can be eliminated, replacing them by effective unary conversion reactions. This allows the rate matrix and memory functions to be calculated from an effective -matrix . These quantities determine the projected equations of motion for the subnetwork proteins s and the subnetwork enzymes e.

The second, final stage of the elimination procedure is now to eliminate the subnetwork enzymes. Elimination of the bulk enzymes does not affect the equations of motion for the subnetwork enzymes, i.e. these do not acquire memory terms nor are the local-in-time terms changed. Looking at the general formula (19) for the memory function, read in terms of the effective -matrix, then confirms that the equations for the e species do not contain memory terms. This makes sense: the subnetwork enzymes are not boundary species to start with, and this is not changed by the effective unary reactions from bulk enzymes.

The projected equation for each subnetwork enzyme thus looks exactly as the original equation in (26). Because of the fast rate scale in this, in the limit each subnetwork enzyme will be in quasi-steady state with its substrate and product. Substituting the quasi-steady state enzyme concentration into the equations for substrate and product then gives again effective unary conversion reactions, with rates as given in (47) for the case of bulk enzymes.

### iv.4 Summary of enzyme elimination procedure for linearised dynamics

The final procedure we have arrived at for constructing projected equations for reaction systems with Michaelis-Menten terms, within linearised dynamics, is remarkably simple: replace each Michaelis-Menten term, whether in the subnetwork, the bulk or on the boundary, by its linearisation around the steady state. This gives effective rates for unary conversion reactions among each substrate-product pair (see (46)).

## V Nonlinear Dynamics

We now want to extend the above approach of eliminating enzymes from the projected equations to the full nonlinear dynamics. Directly transplanting the results from the linearised dynamics is not possible, however: if we use the quasi-steady state assumption for the enzymes as in Section IV, then we get back the full Michaelis-Menten nonlinearities. As these go beyond second order in , they cannot be used directly in our construction of the projected equations, which starts from reaction equations with only linear and quadratic terms as appropriate for a mass action description of unary and binary reactions.

We take as our starting point the nonlinear -matrix, as shown in equation (16), but subdivide this into smaller blocks below in order to single out contributions from enzymes. Focussing though for now just on the distinction between linear and quadratic observables, we have two new kinds of entries. Firstly, mixed linear-quadratic elements as contained in e.g. : these are coefficients of quadratic terms in equations of motion for the concentrations (linear observables), so can be read off directly from the mass action equations. The quadratic-quadratic elements as in are coefficients from equations of motion of concentration products. For a generic product they are of the form

 ∂∂tδxiδxj=δxj∂∂tδxi+δxi∂∂tδxj (48)

Because we are only considering terms up to quadratic order on the r.h.s., we need to insert only the linearised equations of motion for and . All quadratic-quadratic elements of the -matrix are therefore “inherited” from the linearised dynamics. In particular, the above structure of the equations of motion for quadratic observables means that the time evolution of any product containing at least one enzyme factor will contain fast terms.

For the purpose of eliminating the fast degrees of freedom, the nonlinear -matrix can be split into four blocks mirroring the structure of the linear -matrix in equation (29), viz.

 L=(LS,SLS,BLB,SL