Memory functions reveal structural properties of gene regulatory networks

Gene regulatory networks (GRNs) control cellular function and decision making during tissue development and homeostasis. Mathematical tools based on dynamical systems theory are often used to model these networks, but the size and complexity of these models mean that their behaviour is not always intuitive and the underlying mechanisms can be difficult to decipher. For this reason, methods that simplify and aid exploration of complex networks are necessary. To this end we develop a broadly applicable form of the Zwanzig-Mori projection. By first converting a thermodynamic state ensemble model of gene regulation into mass action reactions we derive a general method that produces a set of time evolution equations for a subset of components of a network. The influence of the rest of the network, the bulk, is captured by memory functions that describe how the subnetwork reacts to its own past state via components in the bulk. These memory functions provide probes of near-steady state dynamics, revealing information not easily accessible otherwise. We illustrate the method on a simple cross-repressive transcriptional motif to show that memory functions not only simplify the analysis of the subnetwork but also have a natural interpretation. We then apply the approach to a GRN from the vertebrate neural tube, a well characterised developmental transcriptional network composed of four interacting transcription factors. The memory functions reveal the function of specific links within the neural tube network and identify features of the regulatory structure that specifically increase the robustness of the network to initial conditions. Taken together, the study provides evidence that Zwanzig-Mori projections offer powerful and effective tools for simplifying and exploring the behaviour of GRNs.

Biological systems are complex, comprising multiple interacting components. In many cases this complexity makes it difficult to identify underlying mechanisms and to understand the function of a system. Gene regulatory networks (GRNs) are an example of this problem [1]. A GRN comprises the set of interacting genes responsible for the development, differentiation or homoeostasis of a tissue and provides a formal system-level, causative explanation for gene regulation. In physical terms, a GRN consists of modular DNA sequences—cis regulatory elements—that bind to specific sets of transcriptional activators and repressors, which control the expression of associated genes. Some of the regulated genes are themselves transcriptional regulators. Thus, at the core of a GRN is a recursive set of regulatory links that forms a transcriptional network, the dynamics of which is responsible for the spatial and temporal patterns of gene expression.

Attempts have been made to map large transcriptional networks, yet even for relatively small networks, the number of links and the feedback within the system make intuitive understanding difficult to obtain. Various computational models have been developed to address this. Logical models, which describe regulatory interactions qualitatively, provide a flexible and simplifying formalism to explore and understand the behaviour of a network [2]. However, these approaches are unable to capture subtler features of a network that depend on specific aspects of the timing or concentration of components of the network. For this, continuous models based on, for example, ordinary differential equations (ODEs) are often employed [3, 4]. These models describe gene regulation in much greater detail, and there are well-developed mathematical theories that provide powerful tools to distill the dynamical details of such systems. These have been used successfully to gain insight into the operation of GRNs and suggest explanations for otherwise difficult to understand behaviours, including emergent phenomena such as self-organisation, oscillations, spatial patterning, and scaling of pattern size [5–8]. A disadvantage of this approach, however, is that ODE models often require a large number of parameters, including binding affinities, degradation rates, production rates, etc. While it is often possible to numerically analyse large systems, the power of analytical tools is in many cases lost due to lack of knowledge of parameters or the complexity of the network.

Various techniques have been developed to reduce the complexity of models yet preserve specific features of the behaviour of the system. For Boolean networks it is possible to remove specific components while restructuring other parts of the network to conserve the logical structure, but this requires quite specific network topologies [9]. In the case of mass action reactions, one approach to simplification is to generate intermediate species in a systematic way such as to make different systems comparable. This then allows the identification of motifs that can be simplified. This approach comes with the natural limitation that the species involved cannot be chosen a priori and must conform to a particular network structure [10]. Other methods to reduce complexity involve removing nodes or perturbing the system and analysing a posteriori the effect that this has on the dynamics of the network. This can be helpful to develop intuition, but generally requires the investigation of a large number of combinations of perturbations [11, 12]. A further approach, although applicable only in rather specific cases, is to enforce a timescale separation between network nodes that is chosen not to perturb the dynamics significantly [13]. One can find “morphisms” that relate network structure to function in order to then simplify a set of nodes into potentially simpler motifs; this can reveal an interpretation of how certain motifs work but does not provide a direct conversion from the original network [14]. Less formal methods for replacing well understood network motifs with simpler motifs producing similar behaviour have also been extensively explored [15].

Along similar lines, entire parts of a network can be substituted with effective nodes with more complicated dynamics that produce similar output for a set of chosen species, although the feasiblity of this depends on network structure and the complexity of the original dynamics [16].

We focus in this paper on reducing model complexity by tracking a subnetwork that is embedded within the remainder—or “bulk”—of a larger system. Analysing the behaviour of the subnetwork and its interaction with the bulk can help reveal and rationalise the properties of the entire network. With any coarse graining approach, including those discussed above, a compromise must be reached between the precision of the method, i.e. how well it captures the dynamics of the full network, and the simplicity and interpretability of the resulting description. For example, the introduction of intermediate species can produce model reductions with high precision but these may not be interpretable [10], while morphism approaches provide interpretable insights into the structure but capture the full model dynamics only for specific initial conditions and system parameters [14]. We consider important for interpretability the ability to choose a subnetwork, guided by biological relevance or the availability of data; none of the methods described above allow this and instead identify a subnetwork based on their internal criteria.

One class of methods that offers the potential to balance accuracy and interpretability are Zwanzig-Mori projections [17, 18]. Originally developed to allow the extraction of macroscopic equations from a microscopic description of the dynamics—a brief overview can be found in e.g. [19]—these methods have since been used to separate a network into an arbitarily chosen subnetwork and bulk in ways that preserve substantial features of the original temporal dynamics [20, 21]. Used in this way, the approach describes the concentration of components in the subnetwork in detail, while the activities of the species in the bulk are replaced with so-called ‘memory functions’. These memory functions are derived from the detailed kinetic description of the remainder of the network and summarise how, by acting via the bulk, the past behaviour of the subnetwork influences its current state (Fig 1). The resulting memory functions are functions of time difference, describing the amplitude of a signal returning from the bulk some specified time after the original signal left the subnetwork. Additionally, it is possible to separate each memory function in order to analyse through which bulk species the signal is flowing, thus providing a means to gain a more systematic understanding of a system.

Evolving from its original applications to critical dynamics and supercooled liquids near the glass transition [22], this approach has been applied to biochemical networks [20], where it has been used to analyse the dynamics of signaling in the EGFR (epidermal growth factor receptor) network. It has since been developed further to include also enzymatic Michaelis-Menten reactions [21]. Here we apply the Zwanzig-Mori approach to the analysis of GRNs. We derive a general method that is applicable to any thermodynamic state ensemble representation of protein dynamics including transcriptional networks. These thermodynamic equations mechanistically represent interactions between proteins via transcriptional regulation [23] and can be expanded into mass action reactions to which the Zwanzig-Mori projection formalism can be applied. This approach not only justifies a heuristic method where the thermodynamic equations are expanded to quadratic order around a steady state to obtain a set of effective unary and binary reactions, but also opens up the possibility of developing more advanced projection approximations. We test the method on a cross repressive motif, where we show that we are able to qualitatively reproduce the dynamics of the system and obtain a simple intuitive interpretation of the memory functions. Finally, we use the approach on a GRN from the vertebrate neural tube, a well characterised developmental GRN composed of four interacting transcription factors [24]. This reveals the importance of specific links within the network and identifies features that appear to be present primarily to increase the robustness of the network to initial conditions, rather than maintain the steady states.

Taken together, the study provides evidence that Zwanzig-Mori projections are a powerful and efficient method to simplify and explore the behaviour of gene regulatory networks, with memory functions providing new tools for probing the dynamics near steady state.

Several classes of model have been developed to describe gene regulation. One such class consists of thermodynamic state ensemble models in which equations are constructed that relate interactions between key components of the regulatory mechanism to gene expression [23, 25, 26]. In these models, all possible states of a regulated gene are enumerated. Each state consists of a specific set of transcription factors bound to a DNA cis-regulatory element, weighted by the affinity of interactions. The rate at which a protein is produced is then represented as a ratio of the weight of the subset of states that promote gene expression, to the total weight of all states. We define the vector x such that it contains the protein concentrations for each gene. We can then generically define the time evolution for any protein concentration, indexed by j, as

Here βj is the decay rate of protein j. We use the summation over n to represent all possible binding conformations for the DNA that produces protein j. The term α(j,n) represents the production rate of protein j from a particular DNA conformation n. The factor w(j,n) is the specific affinity for individual proteins to the DNA that produces species xj. Finally x(j,n) is the concentration of DNA producing protein j and in conformation n. This is assigned in (2) as the weight of the conformation, dependent on the protein concentrations in a mass-action form, normalised by the total weight of all conformations. This ensures that x(j,n) lies in the range [0, 1], so it is not an absolute concentration; any overall DNA concentration scale is therefore to be understood as incorporated into the protein production rates α(j,n). Accordingly one can also think of x(j,n) as the probability of finding a certain conformation n of the DNA producing protein j. As made explicit in (2), the DNA conformation label n is a collection of integers ni counting how many copies of protein i are bound to the DNA.

The Zwanzig-Mori method in general starts from the choice of a set of observables for which dynamical equations are to be obtained. Initially we select, as the simplest possible observables, the deviations from steady state of the concentrations of the chosen subnetwork species. For these observables the Zwanzig-Mori projection leads to a set of dynamical equations of the form

Here δxi are the subnetwork concentration deviations as defined explicitly below and Ns indicates the number of subnetwork species. The Ωji define a rate matrix Ω that represents subnetwork interactions. The ri(t) are so-called random forces that arise from the uncertainty about the initial (time t = 0) state of the bulk concentrations. (For the linear dynamics they can be expressed in closed form [20] but this is not useful in our context as discussed in Sec. Nonlinear projected equations below.)

The key quantities in the eq (3) are the memory functions Mji(Δt), which describe how past concentration fluctuations δxj(t′) influence the current time evolution. However, while the projected equations are formally exact, it is in general impossible to evaluate these memory functions explicitly. One scenario where this can be done is a network of unary reactions, where the full reaction equations are linear in the concentrations [20]. This suggests a heuristic approach to calculating memory functions in the current context, which is to linearize the thermodynamic equations. One can think of this as generating a set of effective unary reactions (at least loosely; see Sec. Justifying the heuristics: Network expansion). Explicitly, one linearly expands the time evolution eqs (1) with (2) inserted. The expansion is performed around a steady state with concentrations yi, i.e. in terms of the deviations δxi = xi − yi. The r.h.s. of (1) then contains only linear terms in the δxi since constant terms cancel out because of the steady state condition. One can therefore write the time evolution equations in matrix form [20]

From the matrix L one can calculate the terms in the (linear) projected dynamics of the subnetwork concentrations from (3) as shown in [20]. We summarize the method here. One separates the set of all (linear in concentration) observables into subnetwork (s) and bulk (b), and assumes as above that these concentration deviations δxi are numbered so that the first Ns are the subnetwork observables while the rest are the bulk observables. The matrix L for the linearised dynamics then separates into blocks according to the different observables:

where labels S and B represent respectively the set of species in the subnetwork and in the bulk.

We use the labels S and B here to represent subnetwork and bulk species respectively, this will be useful once more subnetwork and bulk observables arise (see Sec. Mathematical results). The rate matrix Ω with entries Ωji in the projected eq (3) is then simply the subnetwork block of L,

The memory functions similarly form the elements of a memory matrix, which can be expressed in terms of the blocks of L as [20]

This form allows straightforward evaluation of the memory functions starting from any given matrix L describing the linearised dynamics.

To capture the leading nonlinear corrections to the subnetwork dynamics, it is natural to enlarge the set of subnetwork observables by adding quadratic observables, i.e. products of concentration deviations. The projected equations then become

This form is analogous to (3) but now includes the quadratic observables that we have retained, both in the subnetwork interactions (rate matrix terms) and in the memory terms. Once again ri(t) represents the random forces but in this nonlinear case further characterization would require knowledge of the statistics of the initial bulk fluctuations, which is not generally available. We therefore disregard these terms, making (3) a closed system of equations for the time evolution of the subnetwork concentrations. Note that neglecting the random forces is equivalent to assuming that the bulk concentrations are at their steady state values at the initial time. Formally these initial bulk concentrations should be random rather than deterministic to make the Zwanzig-Mori projection well defined; we follow the strategy of [20] here and use initial Poisson distributions in the limit of vanishing variance. This limit has already been taken in the results above [20]. How large protein copy numbers have to be to allow stochastic effects to be safely neglected is something that could be studied in future work; the answer will depend, among other things, on the degree of nonlinearity of the GRN time evolution equations [27].

Moving on to the nonlinear memory functions, we are again faced with the difficulty that these cannot in general be evaluated explicitly. However, previous work [20] has developed a systematic approximation technique, which is applicable to the case of reaction networks with unary and binary reactions described by mass action kinetics. In such networks the time evolution equations for linear observables have linear and quadratic terms on the r.h.s. Inserting these into the equations for quadratic observables using the product rule gives e.g.

and one sees that cubic terms arise. The approximation in [20] neglects these terms in the spirit of an expansion to second order in the changes from steady state δx. These changes are therefore implicitly assumed to be small. Defining a vector z that collects all δx variables and their products such as δx12, δx1 δx2, one can now again write a matrix form of the time evolution equations:

The rate matrix entries and memory functions can then be calculated [20] from the general formulae (6 and 7), applied to the expanded L matrix constructed as explained above. Because the subnetwork block S now contains linear {s} and quadratic observables {ss}, the matrices have the corresponding block structures. The coefficients Ωji are collected in the block Ωs,s, for example, and the Ω(jk)i in the block Ωss,s of the overall rate matrix Ω. The memory functions Mji(Δt) and M(jk)i(Δt) are similarly contained in blocks Ms,s and Mss,s of the memory matrix M.

To apply the above technique in order to derive nonlinear projected subnetwork equations for GRN dynamics, the obvious heuristic route is again to expand the dynamical eq (1), this time to second order in the δxi, to obtain an effective set of unary and binary reactions. This set defines the matrix L, from which the rate matrix and memory functions can then be found as explained above.

The heuristic method set out above for deriving linear or nonlinear subnetwork equations for GRN dynamics is a useful practical recipe, but its mathematical basis is not obvious. One difficulty is that the thermodynamic GRN eq (1) involve rational functions of the concentrations such as 1/(1 + x1), and while these can be Taylor expanded they have only finite radius of convergence. Not only does the heuristic method then have to throw away an infinite number of higher order terms when truncating the Taylor expansion after linear or quadratic order, but it may implicitly be using the resulting approximation outside the regime where the expansion is even convergent. In addition, while we have loosely talked about the expanded time evolution equations as representing networks of effective reactions, such an interpretation is in no way assured as there is no underlying mechanistic model and hence no constraints that would ensure an appropriate stoichiometry or even positivity of rate constants.

Our mathematical contribution in this paper is to demonstrate that the above limitations can be resolved, and the heuristic approach justified, by initially expanding the thermodynamic eqs (1 and 2) into a network of well-defined unary and binary reactions. The way to achieve this is already suggested by (2): we represent each possible existing DNA conformation as an additional species whose concentration has its own dynamical evolution, such that in an appropriate limit of fast DNA binding and unbinding, it reaches the quasi-steady state (QSS) values of the DNA conformation concentrations (2). The prefix “quasi” refers to the fact that this steady state is for given protein concentrations, which themselves change slowly in time. We can then apply the Zwanzig-Mori method to this network to extract memory functions for the chosen subnetwork. In doing this we can justify removing higher order terms as in [20] as the system is composed of only first and second order reactions. Finally the limit of fast DNA binding and unbinding has to be taken. Our result is that this procedure results in exactly the same projected equations as the heuristic approach described above, thus putting the method on a firm mathematical footing. The conclusion applies not only in the context of GRNs but in all networks with the appropriate timescale separation, provided that all fast observables are treated as part of the bulk. More sophisticated projection methods can then be derived by retaining some fast observables (in our case, subnetwork DNA) in the subnetwork. For the case of Michaelis-Menten dynamics, previous work has shown that such an approach allows one to derive projected equations that retain all subnetwork nonlinearities [21].

More generally, by embedding GRN models based on thermodynamic state ensembles into the broad class of mass action reaction networks, our network expansion method opens up the possibility of incorporating stochastic effects and of applying a broad range of other approximation and model reduction techniques [28].

In this section we first provide our mathematical results, where we use a network expansion technique to provide a justification for the heuristic approach to finding memory functions for GRNs that was described in Secs. Projected equations for linearised dynamics and Nonlinear projected equations above. We then demonstrate applications of the approach with case studies. We first illustrate the memory function method on a simple example system. Finally, the results from the biological application to neural tube patterning are given in Sec. Application to neural tube network onwards.

Our network expansion approach necessitates defining new reaction rates for binding and unbinding of protein to/from DNA that are consistent with the affinities w from the thermodynamic description. To be explicit, we define k(j,n)p+ as the rate constant for protein p to bind to DNA coding for protein j and in binding conformation n, with a similar meaning for k(j,n)p-. These rates will thus describe the interaction each DNA species has with the protein species that bind to it. With these rate constants we write down an equation for the time evolution of each DNA concentration x(j,n):