Memory functions reveal structural properties of gene regulatory networks

The first term on the r.h.s. of (11) represents binding of protein p to DNA conformation (j, n − ep) to make DNA conformation (j, n). Here ep is a unit vector whose i-th entry is δip, i.e. a 1 at position p and zeros elsewhere. Similarly the third term on the right of (11) describes unbinding of protein p. The other two terms capture how DNA in conformation (j, n) can be lost by these two processes.

We next establish the relationship between the rate constants k in the mass action representation (11) and the affinities w of the thermodynamic form. The rate constants need to be chosen to ensure the QSS (2) of the DNA concentrations. The simplest assumption is detailed balance, which requires that in this steady state the net change in DNA concentration from any individual pair of binding and unbinding reactions between two conformations vanishes. Detailed balance means, in particular, that where there are multiple reaction paths from one DNA conformation to another, which might be distinguished by the order in which the various proteins bind, each path carries zero reaction flux rather than having a preferred forward or backward reaction direction. Balancing in this way the net change from binding of protein p to DNA in conformation (j, n) with the reverse unbinding transition, for the QSS concentrations (2), gives

The protein concentrations cancel from this as they should and we are left with a constraint on the rate constants:

To finish our construction, we need to ensure that the DNA concentrations reach the QSS values that we have imposed above with the detailed balance constraint (13). This requires that we make all DNA binding and unbinding reactions fast. Formally we introduce a fast rate factor γ and write all rate constants in the form k(j,n)p±=γk˜(j,n)p±, with the understanding that we will take the limit γ → ∞ at fixed k˜. This construction is analogous to the one previously used for Michaelis-Menten reactions [21]. As there, we also need to ensure that the changes to protein concentrations from binding to DNA are negligible, as such effects are not described by the original thermodynamic eq (1). (We do, however, include them in the finite-γ numerics of Fig 2B.) This is achieved by scaling all DNA concentrations as x(j,n)=x˜(j,n)/γ′. The intuition for this construction is that the amount of protein bound to DNA is, at most, of the order of the total concentration of DNA, which vanishes for γ′ → ∞. Finally, to compensate for the DNA concentration scaling, protein production rates have to be scaled as α(j,n)=γ′α˜(j,n).

With the above rescalings, the limit γ → ∞ ensures that the “fast” species, i.e. the DNA concentrations are always in QSS with respect to the “slow” protein concentrations. The limit γ′ → ∞, on the other hand, means that the effects of binding to, and unbinding from, DNA can be neglected in the equations for the protein concentrations. These statements can be shown mathematically along the lines of the arguments in [21]. After dropping all tildes on rates and concentrations again—this is the notational convention we adopt for the rest of the paper—the time evolution of the concentrations in our expanded network is given by (1) for the protein equations, while for the DNA species one has

Note that in the above argument we do not require a specific relation between γ and γ′. In reality the rate for producing a single protein from the appropriate DNA, which is of order γ′, is lower than the DNA binding/unbinding rates, which are O(γ). However, quantifying this would be non-trivial as the model described here simplifies many biological steps into each of the two elementary processes of protein production and binding/unbinding. In any event, for our argument to hold we only need both γ and γ′ to be large, and for simplicity we take then γ = γ′ as was done in [21].

We now proceed to use the expanded network to derive linearised projected equations. After the expansion, this requires us to keep track of deviations from steady state of all protein and DNA concentrations. We label these observables by “s” and “b” for subnetwork and bulk proteins as before, and “a” for DNA. Below we will refer to the DNA species in the system as “fast”, as expressed mathematically by the factor γ on the r.h.s. of (14), while protein species are “slow”. Note that the rate constants for protein production are also fast, of order γ′, but that the resulting changes to protein concentrations are slow because of the low concentration of (protein-producing) DNA.

We construct our L matrix by placing the DNA in the bulk, such that S = {s} and B = {b, a}, giving the block form

The expression (7) for the resulting memory function now explicitly involves the fast rate factor γ and the rate constants k˜(j,n)p± used in the construction of our expanded network. These parameters appear in the third column in the second expression for L in (15), which encodes the time evolution equations for the DNA species. Our remaining task is to take the limit γ → ∞, with the aim of obtaining an expression for the memory function that only involves parameters of the original thermodynamic eqs (1 and 2). To take this limit, one notes that all blocks in the third column of (15) are proportional to γ. As the memory function contains an exponential of L, it will then have contributions decaying for time differences of order 1/γ, arising from the DNA dynamics, as well as slow contributions decaying on timescales of order unity.

For γ → ∞ the fast memory function contributions become effectively instantaneous and add to the rate matrix, and only the slow contributions remain in the memory function. This separation into fast and slow pieces of the memory can be performed using the method in [21] because L has exactly the same division into fast and slow blocks as there. It leads to the simple result that the final rate matrix and (slow) memory function for γ → ∞ can be calculated from an effective L-matrix no longer involving the fast degrees of freedom. To state the expression for the effective L, we first isolate the slow piece of our original L:

(The zero blocks here arise from the fact that subnetwork and bulk proteins do not interact directly in our GRN setting, but only via the DNA species.) The effective L-matrix that gives us the rate matrix and memory function for γ → ∞ is then of the form Leff = L∖a + ΔL∖a, which written out reads:

Here the second term including the minus is the additional contribution ΔL∖a from the fast degrees of freedom. The matrix Leff can be inserted directly into (6 and 7) to obtain the projected equations for linearised dynamics, with the fast rate limit already taken.

As in [21], the result (17) has a simple interpretation: Leff can be obtained by eliminating the fast degrees of freedom from the linearised time evolution equations using a QSS condition. Using the block form (5), these equations can be written as

where δxsT, δxbT and δxaT are vectors collecting the concentration deviations of subnetwork proteins, bulk proteins and DNA species, respectively. All blocks of L appearing in (18c) are proportional to γ. This justifies elimination of the fast variables using a QSS condition for these fast (DNA) degrees of freedom, effectively setting the r.h.s. of (18c) to zero. The fast species concentrations are then expressed as:

If we now substitute this back into the time evolution eqs (18a and 18b) we obtain:

Writing these time evolution equations in matrix form gives exactly the effective L-matrix Leff defined in (17), as claimed. Two terms have been highlighted here by brackets for later comparison.

So far we have shown that the rate matrix and memory functions for the linearised dynamics can be found by applying standard projection results to a set of time evolution equations for only the slow degrees of freedom. This reduced set of equations is obtained by first expanding the thermodynamic equations into a set of mass-action equations, linearising, and finally eliminating the fast variables using a QSS assumption. We next show that this procedure is equivalent to directly linearising the original thermodynamic equations, for any case where there is a timescale separation. This generalises previous work [21] that only considered separate (non-interacting) fast species. To verify the equivalence, it is useful to write the expanded mass-action equations in the generic form

where xs, xb and xa are generic components of xs, xb and xa, respectively. This expanded description is constructed so that for γ → ∞, when the xa can be replaced by their QSS values xa*, one recovers the original thermodynamic equations. These can therefore be written as

Expanding now the thermodynamic equations to linear order around a steady state one has

Here all derivatives are evaluated at the steady state, and the terms in the second line arise from the variation of the QSS values of the fast variables. The required coefficients of the type ∂xa*∂xs can be found by differentiating (24) with respect to the relevant slow variable, here xs:

We can now show that (25) is identical to (20) derived above. To see this, we note that by linearising the expanded mass action eq (21) one can identify the entries of L: for example, Lb,s has entries ∂Rs/∂xb, Ls,a has entries γ ∂Ra/∂xs etc. The QSS coefficients arising from (26) can then be written in matrix form as ∂xa*/∂xs=-(Ls,a(La,a)-1)sa, with a similar expression for ∂xa*/∂xb. Inserting these expressions into (25) does indeed lead directly to (20), with matching terms indicated by the brackets. An exactly analogous argument demonstrates the equivalence for ∂t δxb.

With the above arguments we have justified the heuristic method of Sec. Projected equations for linearised dynamics for obtaining linearised projected equations for the dynamics of a subnetwork within a larger GRN written in thermodynamic form: expand the GRN equations to linear order in protein concentrations around a steady state. Then construct from this the matrix Leff and partition this into “s” and “b” blocks according to the chosen subnetwork–bulk split. Finally determine the rate matrix and memory function matrix from (6 and 7). We invoked an expanded network of binary reactions involving DNA conformations to derive that this is the correct method, but note that its application only requires the original thermodynamic equations as input. Therefore the resulting projected equations are, as they should be, independent of the details of the rates k˜(j,n)p± used in the construction of the expanded network.

To capture more of the nonlinearity inherent in GRNs one can extend the approach described in the previous sections, to include terms that are quadratic in the deviations from a steady state. The memory function then contains fast and slow contributions, with the former being transformed into additional rate matrix terms when γ → ∞. The resulting rate matrix and slow memory functions can be obtained from an effective L-matrix involving only slow (linear and quadratic) observables. As before, this Leff can be viewed as arising from a QSS elimination of the fast observables. Finally, one can show that this explicit elimination of the fast observables is equivalent to directly expanding the original thermodynamic equations to quadratic order in concentration deviations from the steady state: this justifies the heuristic method given in the Sec. Nonlinear projected equations.

Except for the last step, the above chain of argument is analogous to the linearised dynamics case, being based entirely on the identical structure of the partition—compare (5) and (S1 Text, Eq. S.1)—of L into fast and slow blocks. The final reduction to a quadratic expansion of the thermodynamic equations is more subtle because the projection approach treats quadratic observables as distinct and not directly tied to products of linear observables. We defer the details to S1 Text and note only that our derivation there both simplifies and generalises that in [21], from enzyme species that do not interact with each other to arbitrary interactions of the fast degrees of freedom.

We applied the Zwanzig-Mori projection method to a cross repressive motif to illustrate our approach. Fig 2A shows the expansion procedure for this very simple regulation network: each of two genes is repressed through one binding site for the other factor. This system can be expanded to represent all possible DNA conformations, bound and unbound, as shown graphically by the four outer nodes in the network in the centre of Fig 2A. As the arrows indicate, in the expanded network unbound DNA can bind protein to form bound DNA; the reverse process is also possible. In this way one creates a set of reaction equations that has been extended and brought into mass action form by including the concentration of DNA conformations. The equations used for Fig 2 are of the form of (1) & (2):

and conversely for g2Prot. In numerical simulations of the expanded network we first checked the fast rate limit, which we control by setting the ratio, γ, of typical DNA reaction rates to protein production rates. Fig 2B shows that already with a moderate value of γ = 10 the predicted time courses of protein concentrations are almost identical to those predicted by the thermodynamic equations. As γ → ∞ the trajectories from the mass action equations and the thermodynamic equations become identical.

We chose to separate the system so that one protein species, g1Prot, is in the subnetwork and the other, g2Prot, is in the bulk. After expanding into protein and DNA species we proceeded to place all DNA species along with g2Prot in the bulk. Because we have performed the expansion around a steady state, the method is more precise when close to the steady state as expected (Fig 2D). Note that the cross-repressive motif with one binding site is not bistable. It could be made bistable by adding binding sites; with the current approach the behaviour around each steady state would then have to be analysed separately, as we also do in the neural tube application in the Sec. Application to neural tube network. More general projection methods can overcome this limitation as detailed in the Discussion, although it is clear that there are limits of principle here. For example, where the dynamical choice of steady state depends crucially on the initial concentration of a bulk species, this effect cannot be captured by a subnetwork description that by definition does not track this species.

Next we analyse the memory functions in our cross-repressive example, which can be thought of as describing the strength of a signal returning from the bulk a specific time after the original signal entered the bulk from the subnetwork (Fig 2E). The amplitude of the leading linear memory function—where the signal is δx1(t′)—is positive. This makes sense as, in the cross repressive motif, a protein is inhibiting its inhibitor, thus promoting its own production. The second order memory function governs the effects of the signal δx12(t′). It is, in this case, negative and acts as a correction that captures nonlinearities beyond the linear memory. We observe that the decay in time of the memory functions is determined by the decay rate of the bulk protein. This is consistent with the memory function describing a protein repressing its own repressor, hence the effect only lasts as long as the repressor protein is present. Overall, the memory functions provide compact and readable descriptions of how the bulk modifies and influences the activity of the subnetwork.

Reducing part of a network into memory functions offers a way to simplify and analyse the effect of the factors in this bulk part on the remaining subnetwork. To illustrate this approach we chose to analyse a four gene network involved in the embryonic patterning of the vertebrate neural tube [29, 30]. The developing neural tube is a well characterised example of developmental pattern formation. In this tissue, the secreted molecule, Sonic-Hedgehog (Shh) forms a ventral to dorsal gradient [31]. This in turn generates discrete domains of gene expression that define the progenitors of the distinct neuronal subtypes that comprise spinal cord circuitry (Fig 3A and 3B). A transcriptional network has been identified that is controlled by graded Shh signaling and is responsible for specifying ventral progenitor domains [30, 32]. Genetic and molecular experiments have determined the activity of the components of this network (Fig 3C) and documented the temporal dynamics of pattern formation in vitro and in vivo [24, 30, 33]. For the three most ventral domains of progenitors, four transcription factors are essential, Nkx2.2, Olig2, Pax6 and Irx3. The most ventral domain (p3 domain) is characterised by high Nkx2.2 and low levels of the other three proteins; the domain adjacent to this (pMN domain) has high levels of Olig2, medium levels of Pax6 and low levels of the other two proteins. Dorsal to this, the p2 domain expresses high levels of Pax6 and Irx3 and low levels of the other two proteins.

Mathematical models have been formulated based on the experimental data and these are able to replicate key aspects of cell patterning [24, 33] (Fig 3D). We take advantage of these to develop a Zwanzig-Mori projection of the system (using equations as described in [24] & S3 Text). To this end, we chose Nkx2.2 and Olig2 to be the subnetwork species, given that they are receiving direct input from Shh (Fig 3C), and replaced Irx3 and Pax6 with memory functions.