Memory functions reveal structural properties of gene regulatory networks

We first examined the properties of the linear memory functions. These are the most substantial contributions, close to the specific steady states, of the species in the bulk (Pax6 and Irx3) on the temporal changes in the activity of Nkx2.2 and Olig2. We note that the system is multistable (Fig 3E) and different combinations of steady states are available at different positions along the dorsal ventral axis (Shh gradient). We therefore analysed each possible steady state along the neural tube. In this case, as the bulk species are mutual repressors with the subnetwork species, they form a positive feedback loop with the factors in the subnetwork. As a consequence, each of the linear memory functions is positive and exponentially decaying on the time-scale of protein degradation of the bulk species. Thus the relative contribution of each of the memory functions can be assessed by comparing their amplitude, defined as the memory function value at zero time difference.

We calculated the amplitudes of the memory functions and determined the memory effects that Nkx2.2 and Olig2 receive over time. Furthermore, we used the method developed in [20] to decompose the memory into contributions from memory signals passing through the two different bulk species (Fig 4A–4D). This allowed us to determine the importance of specific regulatory interactions between transcription factors (TFs) at every position of the neural tube model, and assess their respective contributions. The results indicate that, for the most part, there is only a small memory amplitude of Nkx2.2 to the past of Olig2 passing through Pax6 (Fig 4B; note the logarithmic y-axis): the memory of (past) Olig2 on Nkx2.2 is largely dominated by memory through Irx3. The only exception is for steady states that are not reached during normal neural tube patterning (lower pair of curves in the neural tube position range ≈ [0.2 … 0.45] in Fig 4B). These steady states exist because the system is multistable, but the initial conditions that lead to them are incompatible with physiological conditions. We additionally observe in Fig 4A that the memory of Nkx2.2 to itself is for the most part shared between Irx3 and Pax6, with no particular TF being dominant. There are again exceptions from this, but only in steady states that are biologically unreasonable. In the case of the memory of Olig2, it can only be influenced by Irx3 as is clear from the structure of the network in Fig 3C. Channel decomposition is therefore unnecessary. We observe that the memory effects of Olig2 to itself are strongest in the p3 domain, and the opposite is true (memory of Nkx2.2 on Olig2 is stronger) in the pMN domain.

We observe in Fig 3 that the order of magnitude of the memory amplitudes changes substantially with neural tube position. This is a consequence of the system becoming more or less sensitive to concentration fluctuations of a given species within the network. An example of this is Nkx2.2 memory to itself via Pax6 in the pMN domain (Fig 4A), where dorsally the lower levels of Shh mean that the binding sites of Nkx2.2 are less occupied by the active form of Gli (Gli is the transcriptional effector downstream of Shh that binds to Nkx2.2 & Olig2, see [24] for details on how this is implemented). Active Gli sets the rate of production of Nkx2.2 and thus has a direct effect on the amplitude of the memory functions. In addition to this, as one moves dorsally in the pMN domain, the steady state concentration of Pax6 is increasing (Fig 3E, also observed experimentally). This increase causes the system to become insensitive to fluctuations of Pax6 proteins as Pax6 binding sites are typically already saturated. This phenomenon is at work also in the other memory function amplitudes that change across neural tube position as seen in Fig 4A–4D.

We next examined the second order memory functions. These provide corrections to the linear memory terms and encode information about how the bulk nodes drive the system towards a steady state. Quantitative accuracy in capturing the dynamics of the full nonlinear network will also be improved to the extent possible within an expansion around a steady state, but such accuracy is not our main focus—it is the additional insights gained from the memory functions that we are after.

It is apparent that once memory terms are included, different trajectories in the subnetwork plane cross. This is visually most obvious with the full nonlinear memory terms (S1C and S1G Fig) and consistent with the full dynamics (S1D and S1H Fig). The crossings arise from the fact that, in the presence of memory, systems in the same state can follow different trajectories towards a steady state depending on their past. In the full dynamics, the additional information as to which path a system will take is contained in the bulk species, which effectively hold the information on the past of the system. In the case of the memory functions, the memory is not contained in network nodes outside of the subnetwork, but is represented explicitly via the memory terms, with the time difference dependence of the memory functions indicating, for example, the timescale of the memory. This offers a useful abstraction: while it is not easy to visualize the four-dimensional concentration space for the full dynamics, the projection approach reduces the system to two components with memory of their own past.

In the nonlinear memory functions, we can assess the importance of the link from Olig2 to Pax6 by considering the squared deviations of Olig2: the Pax6 channel contribution to this memory necessarily involves signals being propagated via this Olig2→Pax6 link. We find that the Pax6 channel makes a very small contribution to the nonlinear memory amplitude in comparison to the Irx3 channel. This is shown in Fig 4E for the p3 domain but is true also for all other states. We additionally checked the time-dependence of the memory function: from Fig 4F one sees that the relative contribution of the Pax6 channel to the nonlinear memory of (past) Olig2 squared on Nkx2.2 remains small at all time-differences. Together with our results above for linear memory functions, we thus conclude that the repression of Pax6 by Olig2 is dispensable for maintaining an established dorsal-ventral pattern. We plot an example of nonlinear memory functions (Fig 4G and 4H) to illustrate that their dependence on time difference can include non-monotonicities and sign changes that mean memory amplitudes do not tell the full story of the effect of nonlinear memory. The effects of including the memory terms can be confirmed and further visualised by considering the effective drift vector, which contains the rates of change of the concentrations of Nkx2.2 and Olig2. We describe this approach in S2 Text, where we also discuss further how specific contributions of the nonlinear memory affect the dynamics, for example in the approach to the state in the p2 domain.

The analysis of both linear and nonlinear memory functions in the previous section suggested that the repression of Pax6 by Olig2 is not critical for the dynamics of the steady states observed in neural tube patterning. We based this on the observation that the amplitude of memory that is transmitted through this link is relatively small compared to the other memory functions, for all biologically relevant steady states.

To test this prediction of the memory function analysis, we removed the repressive link from Olig2 to Pax6 in the full model. Consistent with our prediction of the relative insignificance of the link from Olig2 to Pax6, the resulting steady state concentrations of the simulations are similar to those in the original system (Fig 5A). The only qualitative change is the increased level of Pax6 between neural tube position 0.2 and 0.7, which results from the loss of the repressive link. But the spatial aspects, such as the positional sequence of genes and qualitative concentrations, of the domains are otherwise unaltered.

Since the repressive link between Olig2 and Pax6 has been experimentally documented [33–35] but from our analysis is not required for the dynamics around the steady state, we sought to understand what purpose it might serve during development. We performed simulations from varying initial conditions for both WT (Fig 5B) and the system lacking the link (Fig 5C). These show that the p3 domain is reached from only a small range of initial conditions following the removal of the repressive link between Olig2 and Pax6. To rationalise this, we performed a bifurcation analysis (Fig 5D). This indicated that the removal of the Olig2 inhibition of Pax6 markedly increased the range of positions at which the Olig2 steady state is present, thus facilitating its invasion into the Nkx2.2 domain. In particular, under initial conditions of high Pax6 and Irx3, the pMN domain is induced ventrally. In the unperturbed network it is not possible to reach this steady state in the ventralmost part because the system is monostable, allowing only the p3 domain to arise. On the other hand, in the absence of Olig2 inhibition of Pax6, the ventral region is bistable. Hence, depending on the initial conditions a pMN fate can be reached (Fig 5D). Taken together, this analysis suggests that the repressive link from Olig2 to Pax6 contributes to the robustness of the system for patterning the neural tube.

A potential explanation for the seemingly dispensable regulatory interaction between Olig2 and Pax6 arises from considering the normal development of the neural tube. During embryonic neural tube development, neural progenitors are generated by the process of neural induction and this initiates the expression of neural genes including Pax6 and Irx3. At the same time, cells of the notochord (which underlies the ventral midline of the neural tube) begin to secrete Shh. Shh spreads into the neural tube and neural progenitors respond to the signal. The consequence is that in the ventral neural tube neural progenitors begin responding to Shh at intermediate levels of Pax6 and Irx3. A more appropriate initial condition, representative of the in vivo situation, is for Pax6 and Irx3 to be somewhere between 0 and maximum. Moreover, heterogeneity between cells, as well as variation in the exact timing of developmental events along the rostral-caudal axis of the embryo, means that the precise initial conditions—the levels of Pax6 and Irx3 when cells respond to Shh signalling—will vary. Hence the robustness provided by the Olig2-Pax6 link may play an important role in ensuring reliable pattern formation in the developing neural tube.

We have developed a generally applicable method for applying Zwanzig-Mori projections to transcriptional networks that allows the analysis of subnetworks extracted from a larger network. We demonstrated the approach on a simple genetic cross repressive motif, to illustrate the minimal example of a memory function, and on a more complicated, but well-characterised, transcriptional network operating in the ventral part of the vertebrate neural tube. This showed how the method allows the function and importance of specific links within a network to be defined in an intuitive manner. We used these insights to identify structural features of the neural tube network that appear primarily to increase the robustness of the network to initial conditions, rather than maintain its steady state.

The projection technique is straightforward to implement using the methodology described, and can be applied to any GRN following the thermodynamic formalism (including activation, repression, competitive and cooperative binding). Further protein-protein interactions could readily be incorporated as these would be represented by simple first order (unary) or second order (binary) interactions. Thus, in addition to gene regulation, nonlinear protein-protein interactions mediated by enzymes such as those found in signal transduction pathways could be included in the analysis [21]. It would also be possible to include Hill functions as long as these possess integer exponents [27]. Note that our treatment relies on mass action kinetics, which implicitly assumes that the different molecular species are well mixed by diffusion processes that are faster than any of the reaction kinetics [36]. The size of the matrices involved in constructing the memory functions scales as the square of the number of bulk species, which in our experience means that networks with up to 200 nodes can be investigated without difficulty. Thus, the generality of the method and the ability to implement it algorithmically provide a comprehensive mathematical toolkit to simplify and analyse dynamical systems describing a range of cellular and molecular processes.

In our mathematical treatment we took advantage of an argument of timescale separation for molecular mechanisms operating in transcriptional networks, in which fast processes (DNA-protein interactions) are assumed to be in QSS compared to the slow processes (changes in protein concentration). We also assume the effective concentration of DNA species to be small in comparison to that of protein species. We consider these to be reasonable assumptions because DNA-protein interactions occur at a much faster rate than the changes in protein concentrations produced by transcription and translation [37]. Moreover, there are usually at most four DNA copies of each gene per cell, whereas the number of individual protein molecules is normally considerably higher. The production of individual proteins per DNA copy number, i.e. the relevant rate constant, can then be large while the resulting relative changes to protein concentration are slow. These assumptions allowed us to construct an extended reaction network involving fast molecular (DNA) species and at most binary reactions. This construction is sufficiently general for it to be applicable to other systems with a timescale separation, as long as they satisfy the above conditions for the relevant concentrations and rate constants.

An important component of this Zwanzig-Mori projection method is that the memory functions describe all the behaviours of a system in the vicinity of a steady state. Even though this comes at the cost of predictions becoming increasingly inaccurate the further the system is from the fixed point, the information around the steady state is sufficient to distil relevant properties of the network. We therefore employ the method as a tool for probing the dynamics near a chosen steady state, rather than for tracking dynamics. The benefit is that the resulting memory functions are substantially simpler than the original system. There is no need to perform simulations to analyse memory functions or to test different positions in phase space as a steady state is approached. Moreover, the method can also be applied to an arbitrarily chosen subnetwork of a generic network whereas other model reduction approaches—as discussed in the Introduction—do not allow this, and are also more restrictive as they require specific network motifs or functional forms of the dynamics.

The memory functions produced by the Zwanzig-Mori approach allow exploration of the amplitudes and timescales of the interactions between the subnetwork and bulk. This can reveal information that is not otherwise easy to discern from the parameters of the full dynamics or from the functional form of the equations. Each memory function represents the total of the contributions from a given subnetwork component that feeds back to the subnetwork over time after passing through the bulk. This total memory can be decomposed into channels that describe the flow of signals through specific components in the bulk. In this way, the channels that dominate a memory effect can be identified, providing insight into the dynamical mechanisms responsible for achieving and maintaining a steady state. The memory functions thus help to extract features of the dynamics that would not be readily detectable in direct simulations of the time evolution equations.

An example of the type of insight provided by the projection method came from our analysis of the neural tube network. We focused on a subnetwork comprising Nkx2.2 and Olig2. Decomposition of the memory functions suggested that an experimentally documented repression of Pax6 by Olig2 [34] plays only a minor role in sustaining the steady state pattern. Consistent with this prediction, simulations in which this regulatory link had been eliminated in the full model showed qualitatively unchanged steady states. We note that removal of the repressive link between Olig2 and Pax6 is distinct from the removal of Olig2 itself [24, 30, 33]. In the Olig2 mutant, the repressive effect of Olig2 on both Nkx2.2 and Pax6 are eliminated whereas the removal of the Pax6-Olig2 link leaves the repressive effect of Olig2 on Nkx2.2 intact. The memory function analysis thus raised the question of what purpose the Olig2-Pax6 regulation might serve and prompted us to explore the transient system dynamics, specifically the effect of the initial conditions, which are unknown in the in vivo system. Simulations in the absence of the Olig2-Pax6 repression, compared to the full model, revealed a marked increase in sensitivity to these initial conditions. Comparison of the bifurcation diagrams of the systems, with and without the Pax6-Olig2 link, indicated that the well-defined Nkx2.2 monostable region (p3 domain)—which appears at high levels of Shh signalling in the full system—was replaced by a region of bistability for Nkx2.2 and Olig2 in the absence of the Olig2-Pax6 regulation. In this latter case the steady state induced by high levels of signal therefore depends on the initial concentrations of Pax6 and Irx3.

This suggests that, while not necessary to maintain steady states, the Olig2-Pax6 regulation ensures access to the appropriate steady states irrespective of the initial levels of Pax6 and Irx3. As these may differ between cells and at different locations along the rostral-caudal axis of the neural tube, the Olig2-Pax6 regulatory link might make the system less sensitive to these variations. This implies that the primary purpose of Olig2 repression of Pax6 is to increase the robustness of pattern formation.

Alternative implementations of the Zwanzig-Mori projection [38–40] have been developed that do not rely on expanding around a fixed steady state and we are currently adapting these for transcriptional networks. While these alternative methods may provide a solution to the specific limitations of our current Zwanzig-Mori projection method, the nonlinear memory functions they produce are much richer. Further work will be needed to fully understand the more complex information they encode. Projection methods might also provide a way to analyse systems for which only partial information is available. Given a known subnetwork that is not able to fully reproduce experimental observations, the incorporation of memory functions provides a means to explore the possible effects of unknown factors. This could then be used to identify plausible network structures that generate these memory functions. Taken together therefore, the projection approach provides tools to simplify, visualise and explore the behaviour of large networks that would otherwise be difficult to analyse in their entirety.

EHD and PS acknowledge the stimulating research environment provided by the EPSRC Centre for Doctoral Training in Cross-Disciplinary Approaches to Non-Equilibrium Systems (CANES, EP/L015854/1) (https://www.epsrc.ac.uk/skills/students/centres/profiles/crossdisciplinaryapproachestononequilibriumsystems/). RPC was supported by the Wellcome Trust (WT098325MA) EHD and JB are supported by the Francis Crick Institute which receives its core funding from Cancer Research UK (FC001051) (https://www.cancerresearchuk.org/), UK Medical Research Council (FC001051) (https://www.mrc.ac.uk/), Wellcome Trust (FC001051 and WT098326MA) (https://wellcome.ac.uk/). The funders had no role in study design, data collection and analysis, decision to publish, or preparation of the manuscript.