Dynamical properties of gene regulatory networks involved in long-term potentiation

The long-lasting enhancement of synaptic effectiveness known as long-term potentiation (LTP) is considered to be the cellular basis of long-term memory. LTP elicits changes at the cellular and molecular level, including temporally specific alterations in gene networks. LTP can be seen as a biological process in which a transient signal sets a new homeostatic state that is “remembered” by cellular regulatory systems. Previously, we have shown that early growth response (Egr) transcription factors are of fundamental importance to gene networks recruited early after LTP induction. From a systems perspective, we hypothesized that these networks will show less stable architecture, while networks recruited later will exhibit increased stability, being more directly related to LTP consolidation. Using random Boolean network (RBN) simulations we found that the network derived at 24 h was markedly more stable than those derived at 20 min or 5 h post-LTP. This temporal effect on the vulnerability of the networks is mirrored by what is known about the vulnerability of LTP and memory itself. Differential gene co-expression analysis further highlighted the importance of the Egr family and found a rapid enrichment in connectivity at 20 min, followed by a systematic decrease, providing a potential explanation for the down-regulation of gene expression at 24 h documented in our preceding studies. We also found that the architecture exhibited by a control and the 24 h LTP co-expression networks fit well to a scale-free distribution, known to be robust against perturbations. By contrast the 20 min and 5 h networks showed more truncated distributions. These results suggest that a new homeostatic state is achieved 24 h post-LTP. Together, these data present an integrated view of the genomic response following LTP induction by which the stability of the networks regulated at different times parallel the properties observed at the synapse.

1. Introduction

Living cells are equipped with a robust and yet plastic analog system, which allows them to respond to environmental inputs. Information sensed from the changing environment is integrated and processed by networks of interacting elements to generate an adequate response. Despite the complexity of the underlying mechanisms, some aspects of cellular behavior are of apparent multistable nature, leading to discrete changes that can last for long periods of time. Cell cycle checkpoints, cell differentiation, and apoptosis represent phase transition of multi-component switches that generate robust (and potentially irreversible) transitions (Siegal-Gaskins et al., 2011). From a genetic control perspective, it has been postulated that attractors in the gene and protein expression dynamics, which are the more stable position to which systems tends to evolve, define the cell's character (Kauffman, 1969). Namely, if individual expression values get close enough to the attractor values, these will remain constant even if disturbed. The characterization of the structure of genetic networks from a dynamical perspective using different theoretical methods (Mestl et al., 1997) predicts two broad regimes. The ordered regime, robust against perturbations, and a chaotic regime, sensitive to perturbations. While robustness is a hallmark of homeostasis, it is reasonable to expect that transitions between cellular states require an enhanced sensitivity. In such scenario, a compromise between robustness and sensitivity could potentially be attained by a rewiring of the network or by the recruitment of different networks.

The change in synaptic efficiency known as long-term potentiation (LTP) represents the cellular correlate for long-term memory in the mammalian brain. From a systems perspective, LTP offers an attractive model of a cellular switch, whereby activation promotes movement to a new cellular state. Just as any other biological switch, LTP accommodates a compromise between robustness to genetic and environmental noise and sensitivity to discriminate meaningful signals. This characteristic is likely to be distributed at different levels of biological organization. For example, LTP requires activation and trafficking of glutamate receptors to the postsynaptic membrane, in addition to protein synthesis and de novo gene expression (Abraham and Williams, 2003). Indeed, specific patterns of gene expression have indeed shown to be regulated at different times following LTP induction, and are crucially involved in the maintenance of LTP (Park et al., 2006; Håvik et al., 2007; Ryan et al., 2011, 2012). Specifically, we reported that the networks derived 20 min post-LTP induction, comprised many transcription factors (TFs), including all members of the early growth response (Egr) family, and were associated with functions such as development, proliferation, and neurogenesis. By contrast the networks derived at 5 h contained molecules associated with calcium dynamics, dendritogenesis and neurogenesis and in the networks derived at 24 h neurotrophin-NFKB driven pathways of neuronal growth were identified. Our analysis also revealed several mechanisms controlling the temporal shifts in gene expression such as regulation of specific microRNA and histone deacetylases. Thus, the variety of functions held by the genes offers a glimpse of the potential complexity underlying the genomic response to LTP.

The involvement of complex gene networks in the maintenance of LTP suggests that a multistable switch is also present in the structure of the gene regulatory networks recruited at different times following LTP induction. From a systems perspective, we expect that the transition between the pre- and post-LTP homeostatic state is paralleled by the genomic response. We propose that while gene networks rapidly and transiently activated following LTP-inducing stimulation will show a less stable architecture, networks recruited later will exhibit an increased stability, being more directly related to LTP consolidation and post-LTP homeostatic state. Using random Boolean networks (RBN) simulations, we addressed this hypothesis by studying the dynamical stability of networks previously identified (Ryan et al., 2011, 2012). We also use the gene expression profiles provided in these studies to explore the overall co-expression network architecture. Following identification of tightly co-expressed modules, we used functional analysis to investigate the intramodular differential connectivity at different times post-LTP induction. Our results offer an integrated picture of the genomic response following LTP and support the conclusion that a new homeostatic state is achieved 24 h post-LTP.

2. Materials and Methods

2.1. Network Topologies

LTP-related gene expression profiles investigated in this study were taken from Ryan et al. (2011, 2012). Briefly, 20 min, 5 h, and 24 h following induction of LTP at perforant path synapses in the dentate gyrus in awake freely moving rats, gene expression profiling was carried out using Affymetrix RAT230.2 microarrays and the functional relationships of the differentially expressed genes (±1.15-fold change; p < 0.05; moderated paired t-test between the stimulated and control hemispheres) was explored using the network analysis tool Ingenuity pathway analysis, version 7 (IPA) (Ingenuity Systems, USA; https://www.analysis.ingenuity.com). In the present study we analyzed the highest scoring network from each time point (N = 35) alongside the yeast transcriptional network (N = 30) as a benchmark for RBN modeling (Lee et al., 2002). The yeast transcriptional network represents potential pathways that yeast cells can use to regulate global gene expression. It provides a useful comparison for our analysis for a number of reasons. First, it was constructed by determining experimentally the binding sites of most of the known yeast TFs. In addition to being comprehensive, the yeast network is of similar size to the LTP networks identified by IPA. Finally, it has been used previously in the literature for RBN models (Kauffman et al., 2003; Karlsson and Hörnquist, 2007; Tuğrul and Kabakçıoğlu, 2010). For a node (gene), the number of incoming connections (edges) is called the in-degree of the node and the number of outgoing connections (edges) is its out-degree. The analyses using RBNs were also applied to two different sets of null-hypothesis random networks. First, an ensemble of 100 random networks was generated for each of the 4 biological networks studied (20 min, 5 h, 24 h, and yeast) by preserving the same number of nodes and edges. In order to construct these random networks, pairs of genes are connected randomly with equal probability from the list of N = 35 or N = 30 genes until the total number of edges of the biological network has been set.

A more stringent control consisted of 4 ensembles of 100 rewired networks constructed such that each of the genes had the same in- and out-degree as the biological network. These networks are constructed by randomly choosing two edges of the biological network and swapping them so that A → B and C → D, become A → D and C → B. Swapping is prevented if either A → D or C → B exist already (Kannan et al., 1997). Hence, these sets of rewired networks preserve not only the number of genes and the total number of interactions, but also the original degree sequence and hence their local connectivity is identical to the biological network. These controls allowed us to discriminate the effects of the network's local structure from the effects of the general topology on the robustness.

2.2. Random Boolean Networks

RBNs represent one of the simplest models for gene regulatory networks (Kauffman, 1969). We assumed the expression of each gene in a network to be only “on” or “off” values (corresponding to transcriptionally active and transcriptionally inactive). If we define the state of the network as the set of values of expression of all the genes, for a network of size N, a bit vector G(x₁, …, x_N) with x = 0, 1 suffices to fully describe it. We implemented the connectivity of a given network as a matrix W of size N × N where w_ij = 1 if gene j acts upon gene i, and w_ij = 0 otherwise.

To model the gene expression dynamics of the network, we allowed the vector G to evolve across discrete time steps. We considered the value x_i of the gene i at a certain time step t to be dependent only on the values at t − 1 of the k_i genes that act on it. For each gene in the network, hence, a fixed Boolean function b_i mapped every possible combination of values of the k_i inputs of i to an updated value for the gene i.

On each interaction, we updated all the values of the state vector G synchronously. This means that given a set of values for the state vector G(t) and a set of Boolean functions B, the dynamics of the system are deterministic and the number of possible states for G is discrete. Furthermore, it is known that using synchronous RBN dynamics leads to the expression values of the genes to converge into a number of recurring states or attractors. Crucially, these states can be regarded as different homeostatic cell states (Kauffman, 1993). The size and number of these attractors characterize the dynamical stability of a network. Note that while we allowed the values of the genes in the network to evolve over time, we kept the connectivity W and the set of Boolean functions fixed.

The choice of the synchronous Boolean approximation was guided by an optimal compromise between both conceptual simplicity and computational feasibility, while still holding the capacity to approximate a general stability characterization of a biological network (Kauffman et al., 2003). Indeed, since the gene networks studied were built based on gene expression profiles that represent averages of gene product, a choice of an asynchronous stochastic schema is difficult to justify.

The networks studied using RBN in the present work (see Section 2.1) do not convey any information on the rules for the interactions that could be used to constrain the sample space of all possible Boolean functions. However, as the present study deals with the comparison of the stability conferred by the architecture of specific networks, it suffices to construct the Boolean functions from a flat distribution. Namely, we opted not to introduce a bias in the outputs of the Boolean functions. This means that for a given Boolean function b_i the probability of the output being one of a particular combination of input values is P(b_i = 1) = 0.5.

2.3. Dynamical Robustness

In general, we consider a stable network one in which small perturbations are not amplified in time but rather converge into the same attractor. In this scenario, arbitrarily changing the values of few genes would not have a dramatic effect—after a few iterations, the perturbed genes would return to their original values. The opposite behavior is also possible. We consider an unstable (or chaotic) network one in which the perturbation of few genes results in a generalized change in expression values after few iterations. To simulate the effect of a perturbation in a network, we first chose a set of initial conditions G_A(t = 0). We generated an additional instance of the network, with the same initial conditions but shifting some of them [a total of H(t = 0) genes] to the opposite Boolean value, so that the two instances G_A(t = 0) and G_B(t = 0) are, in terms of gene expression, H(t = 0) far apart. This measure is known as Hamming distance, and is equivalent to the size of the perturbation. We then ran RBN dynamics in parallel for both G_A and G_B, with the same Boolean functions for τ time steps. We calculated the new updated Hamming distance H(t = τ) between G_A(t = τ) and G_B(t = τ) to establish whether the perturbation converged into the same values, H(τ) < H(0), or propagated across the network, H(τ) > H(0). In the study of dynamical RBN systems, these distinct outcomes are said to fall either on the ordered or chaotic regime, respectively (Fox and Hill, 2001).

Different effects of the perturbations are expected with differing network architectures, Boolean functions or specific initial conditions. As we wanted to characterize the contribution of the topology to the stability of the network, we characterized a network by recording the effects of a perturbation over a large number of initial distances H(0), initial conditions G_A(t = 0) and G_B(t = 0), and random Boolean functions.

In order to visually represent the average effect of perturbations in a particular network, we plotted H(0) (size of the perturbation at the beginning of the parallel runs, x-axis) against H(τ) (y-axis), where τ is a small number of discrete steps. Plotting different values of H(0) (sampling different perturbation sizes) results in a Derrida plot, a popular tool used in RBNs (Derrida and Weisbuch, 1986). While some network architectures tend to absorb small perturbations and the final Hamming distances H(τ) are on average smaller than the initial perturbation, some topologies tended to amplify them—few genes with different values of expression lead to dramatically different network states. Crucially, in the Derrida plots these different behaviors fall in the opposite halves of the plot, with robust architectures represented by curves underneath the diagonal, H(0) > H(τ), and sensitive architectures above, H(0) < H(τ) (Fox and Hill, 2001), and a network whose Derrida mapping appears tangent to the diagonal is said to exhibit criticality. Note that choosing small integer values for τ (shorter dynamics) captures the effects of the network's local geometry, while larger values reflect the general characteristics of the structure of the network, since information has more time to spread across the network (Aldana et al., 2003). In practice, the slope for the small H(0) region reveals the average outcome of a small perturbation. Curves below the diagonal indicate a tendency toward stability (ordered regime) whereas curves above imply instability (chaotic regime). The diagonal H(0) = H(τ) represents the transition from order to chaos.

To construct the Derrida plots, for a given network we assigned 1000 random sets of Boolean rules and for each set of randomly generated rules we ran 100 parallel simulations with random initial values, uniformly sampling different values of H(0). We did not imply any structure in the Boolean rules and these were randomly generated from a flat distribution without any explicit bias.

2.4. Gene Co-expression Network Construction

Using the gene expression profiling data obtained from Ryan et al. (2011) differentially expressed genes were ranked according to the significance of their p-value. For computational reasons and to enhance the signal in the data, we used only the top 1700 genes of each of the temporal contrasts. The contribution to the final list of each of the different time groups was equivalent in terms of number of genes, and resulted in a set of differentially expressed genes across early and late LTP of 4804 genes with p-values ranging from 7.7 × 10⁻⁶ to 6.7 × 10⁻². The total number of genes used for network construction analysis is in the same order of magnitude of other co-expression studies (e.g., Ghazalpour et al., 2006).

For the formation of a gene co-expression network for each of these time points we followed the weighted gene co-expression network analysis (WGCNA) protocols (Zhang and Horvath, 2005) as implemented in the WGCNA package of R software (Langfelder and Horvath, 2008). Briefly, for every pair of genes i, j the Pearson correlation is calculated and transformed into an adjacency measure with a power function, which serves to further separate the highly co-expressed pairs from the weakly co-expressed pairs. We used the scale-free topology criterion to choose the soft threshold p = 5 for the adjacency measure calculation (Zhang and Horvath, 2005). As a measure of connectedness, the topological overlap (TO, Ravasz et al., 2002) was used to compute the similarity between genes, resulting in four undirected weighted networks of the same size (N = 4804) but with varying connectivity values of TO—control (unstimulated hemispheres), 20 min, 5 h, and 24 h. TO can be understood as a measure of “agreement” between the nearest-neighbors of two genes. It has been shown to be one of the most biologically meaningful similarity measures used in gene co-expression analysis (Ravasz et al., 2002).

2.5. Identification and Characterization of Gene Co-expression Modules

For each of the four co-expression networks, we constructed specific gene co-expression modules, clusters of densely interconnected genes. This analysis provided a summary of the networks by reducing their complexity to a small number of modules uncovering potential biological associations. For each of the four co-expression networks, we performed a hierarchical clustering using TO as a similarity measure. The branches of the resulting dendrogram were cut using the default parameters implemented in the WGCNA R package (Langfelder and Horvath, 2008). Modules were considered for further analysis if they contained at least 50 genes. Note that the intersections between the modules are not always empty since module detection was performed independently on each temporal expression dataset. As larger modules in the control and 20 min samples appear to segregate into different sets of genes, we chose to keep all the modules and not to merge overlapping modules (see Supplementary Figure S1). Intramodular functional enrichment was calculated using the topGO R package (Alexa et al., 2006) with with a significance criteria of p < 0.01 (Fisher's exact test).

2.6. Co-expression Network Reconfiguration

For each of the 58 modules obtained by WGCNA, we quantified the modular differential connectivity (MDC, Zhang et al., 2013), which corresponds to the ratio of the average connectivity for any pair of module-sharing genes at time T₁ compared to that of the same genes at time T₂ where w_ij is the TO between two genes in a given network. The statistical significance was assessed by a false discovery rate (FDR) based on permutation of the gene labels (Zhang et al., 2013). MDC > 1 indicates enhanced co-regulation between genes, whereas MDC < 1 indicates reduced co-regulation. The MDC was calculated for each of the 58 networks for each of the three temporal transitions (control → 20 min, 20 min → 5 h, 5 h → 24 h).

3. Results

3.1. Dynamic Stability of Temporal LTP Networks

3.1.1. Gene Networks Recruited Earlier Following LTP have a More Unstable Architecture

To test the hypothesis that the gene networks induced more rapidly following LTP in vivo show a less stable architecture when compared to the network induced later, we have drawn on data from our previously published microarray data studies Ryan et al. (2011, 2012). To understand the complexity of the gene networks regulated following the induction of LTP, we used Affymetrix DNA microarrays to identify genes differentially expressed at 20 min, 5 h, and 24 h post-LTP induction in vivo Ryan et al. (2011, 2012). Analysis of the gene regulatory networks derived using IPA suggested that these networks made an important contribution to the stabilization of LTP. Furthermore, not only were subsets of genes confirmed to be differentially expressed by quantitative qPCR, but also specific microRNA predicted to act as key regulatory hubs within these networks were shown to be differentially expressed in the hours following LTP induction (Ryan et al., 2013; Joilin et al., 2014). Here, we have used RBN modeling to assess the stability of the architecture of the three highest scoring networks as identified by IPA at each time point. These networks were analyzed alongside the yeast transcriptional network as a benchmark for RBN modeling (Lee et al., 2002; Kauffman et al., 2003; Karlsson and Hörnquist, 2007; Tuğrul and Kabakçıoğlu, 2010).

Consistent with our hypothesis, the output of the RBN analysis (Figure 1) demonstrates that the network identified 24 h following LTP induction is considerably more ordered than either of the earlier networks (20 min and 5 h) or the RBN benchmark, the yeast transcriptional network. The curve corresponding to the late (24 h) network lies underneath the others, which means that the average outcomes of perturbations to the gene expression levels do not spread across the network to the same extent (Figure 1A; τ = 1). This observation is even more apparent if the simulations are evolved for more iterations, allowing the new values for the gene expression to be used as inputs for next iteration (Figure 1A; τ = 5) before plotting the Hamming distances H(0) vs H(τ). This amplification of the differences between the temporal networks with longer dynamics indicates that both local motifs and long-distance interactions contribute to the differential stability observed between the temporal networks.

FIGURE 1

Figure 1. Results of the RBN dynamical stability analysis. (A) Derrida plots for the LTP networks previously identified by Ryan et al. (2011, 2012) at different times post-LTP induction (20 min, 5 h, 24 h in green, orange, and blue, respectively) and the yeast transcriptional network (black). The left plot corresponds to the initial Hamming distance H(0) plotted against the Hamming distance after 1 iteration, H(1). Hence, only the nearest-neighbor interactions (local motifs) affect the dynamics. The right plot depicts H(0) vs H(5), where long-distance indirect influences between genes have an effect on the dynamics. Longer dynamics allow to reveal the influence of the overall network structure on its stability. The earlier networks (20 min and 5 h) are more unstable than the 24 h network, and the curve corresponding to the latter lies near the diagonal of the plot, which represents the border between the chaotic (white background) and the ordered regime (gray background). (B) Average time evolution of perturbed fixed points starting from Hamming distance H(0) = 1. This small difference tends to be amplified in these biological networks. The latest network recruited following LTP induction (24 h, blue), shows a less pronounced tendency to amplify the perturbation. Furthermore, from t = 2 to t = 5 the Hamming distance shows a slight decrease. The yeast transcriptional network (black dashed line) lies between the earlier LTP networks and the 24 h. (C) Derrida plots for each LTP network and to ensembles of random networks. The same stability profile shown in (A) for H(0) vs H(5) is shown separately for each of the temporal networks. The range of stability exhibited by the two ensembles of random networks (red shade: same number of nodes randomly connected by the same number of edges; blue shade: same number of nodes, same number of edges, and same in- and out-degree). These contrasts allow to isolate the effect of the specific degree sequence from the effect of the average degree (blue shade vs. red shade). In addition, it shows that if an evolutionary constraint were to act on the degree sequence, the real networks choose the less unstable option among all the possible network architectures with the same degree sequence (namely identical local motifs, blue shade). Each point in the plots is the average over 1000 random rule assignments for 100 random initial conditions (increasing these numbers has no effect on the results). Shades for random networks (red) and rewired networks (blue) correspond to the ranges observed using 100 topologies for each. Hamming distances are normalized by the number of nodes.

The same conclusion can be drawn from the panel depicted in Figure 1B, where the average Hamming distance of changing one random gene [H(0) = 1] is plotted at each consecutive time step. In other words, the two network states differing in only one position are independently evolved over 5 time steps and the distances are monitored at each time step. The amplification of the perturbation is clearly less pronounced than the one observed for the other networks, the yeast network lies between the earlier networks (20 min and 5 h) and the more stable 24 h network.

3.1.2. Stability of Random and Rewired Versions of the Real Networks

Using a similar RBN model, Kauffman et al. (2003) compared the stability of the yeast transcriptional network with networks of the same number of nodes and edges that also preserved the degree sequence (the same sequence of in- and out-going edges); so-called “rewired” networks. The study demonstrated that the yeast network was more stable than these rewired networks, which suggested that an evolutionary pressure may be acting on the network geometry. To assess if that was the case for the temporally specific LTP-related gene networks, we conducted the same analysis by studying the stability of rewired versions of the real networks. We found that the LTP networks lean toward a more ordered regime than their rewired counterparts (Figure 1C), in a manner similar to the yeast transcriptional network analysis. Furthermore, we also analyzed the stability of a less constrained set of random networks that only preserve the number of nodes and the number of edges of the real network (see Section 2). This contrast isolates the effect of the specific degree sequence from the effect of the average degree. The plots indicate that the real networks lie in the less unstable margin of the possible network architectures with the same degree sequence (see Figure 1C). This observation is particularly marked in the cases of the 5 h and 24 h networks, and supports the idea that the stability is not only dependent on the local structural motifs but rather is distributed across the global architecture of the network (Wagner, 2005).

3.2. Co-expression Analysis

Although analysis of the IPA generated networks has provided validated and biologically meaningful data (Ryan et al., 2011, 2012; Joilin et al., 2014), some limitations are inherent to the methodology. First, potential key interactions may be excluded as the interactions of only 35 genes per network have been considered. Secondly, the architecture of each network is directly dependent upon the information contained within the IPA Knowledge base, a manually curated database, which makes these networks susceptible to false negatives. Thirdly, genes that are modestly but consistently regulated at each time point will be excluded if they do not reach the inclusion criteria at any time. Finally, the analysis is incompatible with the identification of a control network, allowing the characterization of a pre-LTP homeostatic state, as networks are based on differentially expressed genes. Thus, to rigorously test our findings that not only are biologically relevant groups of genes regulated following LTP, but that the resultant networks have specific architectural properties and become more stable with time, we used the WGCNA methodology (Zhang and Horvath, 2005) to construct weighted gene co-expression networks based on each of the sample classes (unstimulated hemispheres, 20 min, 5 h, and 24 h stimulated hemispheres). We next assessed the connectivity distributions within these networks.

3.2.1. Identification of Co-expression Modules

Co-expression matrices were formed using TO as a similarity measure, which represents the degree of “connectedness” between two genes (Figure 2). The top hubs in the co-expression networks according to their degree (TO with the other genes in the network are shown in Table 1). Within the TO matrices we identified a total of 58 densely connected modules through hierarchical clustering (Control: 9 modules containing between 69 and 2327 genes; 20 min: 9 modules with 64–1184 genes; 5 h: 24 modules with 31–535 genes; 24 h: 16 modules with 42–1164 genes) (see Table 1 and Supplementary Figure S2). To explore the functional relationships of the genes within these modules, modules were tested for Gene Ontology (GO) term enrichment using the topGO R package. We present a brief summary here and a more comprehensive list in Supplementary Material.

FIGURE 2

Figure 2. (A) The co-expression networks (N = 4804) corresponding to the different temporal microarray datasets (control, 20 min, 5 h, 24 h) represented as heatmaps. The darker shade of red represents higher TO values between a pair of genes at that particular time. The TO measure represents the degree of “connectedness” between two genes, and it is based on the adjacency measure calculated from the Pearson correlation (see Section 2). (B) Scale-free fit index as a function of the soft-thresholding power p for the co-expression networks constructed using the time-course microarray data. The R² fit to a scale-free distribution of the unstimulated control and the 24 h networks (black and blue curves, respectively) are both higher and saturate at lower values of p than the earlier networks (20 min and 5 h co-expression networks, green and orange, respectively). The latter reaches a saturation only of around R² = 0.5. Scale-free networks have been shown to be more robust against small random perturbations, while at the same time are sensitive to specific directed perturbations, which confers them a high degree of sensitivity to meaningful signals. These results are in agreement with the notion drawn from the results using RBNs on the IPA networks. (C) Mean connectivity as a function of the soft-thresholding power p. While the temporal co-expression networks fall into two different categories according to their scale-free distribution fit, the average connectivity does not show a clear temporal-specific pattern. The dashed lines in the plots indicate the value of p = 5, chosen to conduct the module identification. (D) Fraction of nodes with degree k in the above co-expression networks using p = 5 and p = 20. These degree distributions are log-transformed both in the x- and y-axes. The black lines represent the linear model fit with the values of R².

TABLE 1

Table 1. Top hubs in the co-expression networks according to their degree (TO with the other genes in the network).

Consistent with our previous analysis, WGCNA identified a number of transcriptional modulators as hubs in the 20 min co-expression modules. Our results stress the importance of the Egr family, previously reported to be expressed following LTP induction (Cole et al., 1989; Richardson et al., 1992). In particular, Egr1 appears among the top 10 hubs of the overall 5 h co-expression network (see Table 1). In addition, Wt1, a member of the same family, appears as a key regulator in one of the co-expression modules (see Supplementary Figure S2, module brown_20). Its role as a repressor of the other Egr family members (Haber et al., 1991) suggests that it may play an important role in regulating Egr gene expression after LTP. Similarly, the Homer family of TFs has been implicated in LTP (Kato et al., 1997) and Homer2 appears as a hub in a co-expression module activated at 20 min and 24 h (see Supplementary Figure S2, module turquoise_24).

A representative GO term for each module is shown in Table 2. Functions overrepresented in the modules identified at 20 min show “positive regulation of endocytosis,” “neuron part and cytoplasmic microtubule,” “axogenesis,” “calmodulin-dependent kinase activity,” and “transcription from RNApolI promoter” among others. At 5 h, “CNS neuron axonogenesis,” “anion homeostasis and synapse assembly” are salient overrepresented functions. Regulation of gene expression is represented by the GO terms “regulation of DNA methylation” and “chromatin DNA binding,” “histone H3-K27 methylation.” Finally, modules identified at 24 h show “neuron projection membrane,” “histone demethylation,” and “response to calcium” among others.

TABLE 2

Table 2. Summary of the modules identified by WCGNA with at least 50 genes.

3.2.2. Scale-free Distributions are Distinctive in the Pre-LTP and 24 h LTP Co-expression Networks

It is widely accepted that biological networks tend to have connectivity distributions that approximate scale-free distributions (Jeong et al., 2000), which define networks with few highly connected nodes and many sparsely connected nodes (Barabási and Albert, 1999). This structural property may be selected for in biological networks due to its robustness against random perturbations while retaining a high sensitivity to directed signals (Albert et al., 2000). From a perspective in which LTP is considered as a high-level switch, these characteristics are to be expected of the genetic networks associated with the pre- and post-LTP homeostatic states. On the contrary, it is plausible that transient topological rearrangements taking place during the transition exhibit architectures that depart from a scale-free architecture.

The WGCNA methodology transforms the correlation between the profiles of expression of two genes into an adjacency measure using a power function. The degree distribution of the network is hence dependent on the choice of the parameter p for the power function. To examine this dependency and characterize the temporal networks in terms of their resemblance to scale-free networks, we plotted the scale-free fit index of each temporal network (control, 20 min, 5 h, 24 h) as a function of the parameter p (Figure 2B). The results demonstrate that the networks reach an asymptotic maximum fit at different values of p, with the unstimulated control and the 24 h networks saturating earlier (p ≈ 5 − 10) than the 20 min and 5 h networks (p = 20). In addition, the saturation value of the scale-free fit coefficient is higher for the control and 24 h than for the 20 min and 5 h networks. In summary, the control and the 24 h networks fit to a greater degree to a scale-free distribution, independently of the average connectivity (Figure 2C). Furthermore, a closer inspection to the scale-free model fit for the values of p = 5 and p = 20 reveals that the earlier networks following LTP stimulation (20 min, 5 h), have few nodes with a low degree (truncated left side of the distribution) (Figure 2D).

3.2.3. Changes in Connectivity Parallel Gene Up-regulation Following LTP

Interestingly, the co-expression networks corresponding to the later times (5 h and 24ḣ) appear to be more dissociated than the earlier networks (control and 20 min), splitting up into more modules for the same power threshold (see Table 2).

The modular differential connectivity (MDC, Zhang et al., 2013) corresponds to the ratio of the average connectivity for any pair of module-sharing genes at time T₁ compared to that of the same genes at time T₂ where w_ij is the TO between two genes in a given network. We calculated the MDC in the 58 co-expression modules across each temporal transition—control → 20 min, 20 min → 5 h, and 5 h → 24 h. The findings are summarized in Figure 3. While the fraction of modules with a significant increase in MDC is similar along the time samples, the fraction of modules with a significant loss of connectivity increases with time (stacked bars plot, Figure 3D).

FIGURE 3

Figure 3. Module reconfiguration following LTP. (A) Modules identified using the 20 min co-expression data all share a similar trend—they become tightly co-regulated at 20 min to lose the connectivity at 5 h. (B) Modules identified at 5 h exhibit a significant decrease in TO at 24 h. (C) Modules identified using the 24 h data together with two modules identified with the control data show a significant gain in intramodular connectivity at 24 h. Green lines represent significant gain of intramodular connectivity (MDC > 1, FDR < 10%) and red lines represent significant loss of modular connectivity (MDC > 1, FDR < 10%). Gray lines represent no statistical significance for the MDC. (D) Distribution of modules with gain, loss or conservation of intramodular connectivity.

The rapid response elicited after LTP induction at the gene expression level corresponds to the transition observed between the control and the 20 min dataset in terms of MDC. Out of the 58 modules, 57 show a gain of connectivity (MDC > 1), of which 20 are significant with FDR < 10% (a more strict FDR < 1% gives still a total of 10 significant modules). This marked increase in the average TO connectivity takes place in combination with the documented enrichment in up-regulated genes within the set of differentially expressed genes (Ryan et al., 2011).

Our results, however, highlight a global loss of connectivity between the 20 min and 5 h networks (0.187–0.089 average TO, respectively). Yet, the proportion of modules exhibiting gain/loss/conservation of connectivity is very similar (17/17/24, respectively). More interestingly, the fraction of the modules with a significant loss of MDC corresponds to modules that showed an increase in MDC in the earliest transition. These modules subsequently undertake a loss of connectivity in the final transition to the 24 h phase.

4. Discussion

Complex biological networks are expected to maintain a certain level of stability against environmental perturbations. Whereas some classically studied mechanisms such as gene redundancy and epistasis suggest that the dynamical properties of biological networks are restricted to small sets of genes (Sanjuán et al., 2004; Moore, 2005), other authors point in the direction of a “distributed robustness” scenario suggesting that all the regulatory interactions among genes play a role in the dynamical characterization of the network (Shmulevich et al., 2005; Wagner, 2005). These studies suggest that stable, robust attractors for genetic networks may underlie homeostatic cell states. Likewise, the capacity to generate a new phenotype by sensing specific environmental signals represents a crucial property of living organisms. The apparent antagonistic relationship between robustness and sensitivity could potentially be eased by rewiring the network architecture or by bringing other networks into operation. Network architectures recruited transiently during phenotypic transitions do not uphold a homeostatic state and conceivably they are not under a selective pressure to exhibit a stable architecture.

While some biological systems are designed to elicit a graded response to an input, other systems show behaviors that resemble multistable transitions. These systems can range from relatively simple switch-like responses to more complex multistable switches (e.g., Ferrell and Machleder, 1998; Pomerening et al., 2003; Yao et al., 2011). During LTP in particular, the cell state needs to transit from a stable point to another stable point. If this assumption holds for the gene expression profile, LTP induction can be seen as the perturbation needed to shift the gene expression equilibrium to the post-LTP attractor. The neuron has to be able to discriminate the changes eliciting LTP from the environmental noise. From this perspective, LTP represents a cellular mechanism, which operates as a switch. The essentiality of the genomic component for the maintenance of late-LTP suggests that the gene regulatory tier of information processing may have characteristics of a high-level switch. Ultimately, other tiers of regulation acting at different levels of organization (networks of neurons and nerve fiber projections between brain areas) act jointly to create, maintain, and retrieve memories.

4.1. Dynamical Stability of Temporal LTP-related Gene Networks

Previous work shows that there is a critical temporal window after LTP induction in which a rapid nuclear response takes place (Nguyen et al., 1994). This early phase is in the order of minutes, and is characterized by a rapid up-regulation of gene expression, which persists for many hours. The network identified 20 min after induction represents these early response genes. The set of genes identified 5 h post-LTP induction are not closely related to the 20 min early responding genes, as demonstrated by the expression profiles (Ryan et al., 2012). The nature of this rapid transcriptional response following LTP induction suggests that the underlying mechanisms are facilitating a switch-like response. In this line, Saha et al. (2011) documented recently the presence of stalled RNA polymerase II in LTP immediate early genes, which they interpreted as a mechanism for the rapid neuronal induction observed. However, other mechanisms may be acting jointly at different levels to complement the gene expression trigger.

The 24 h post-LTP induction represents a temporal and functionally different data set, as indicated both by the lack of overlap in gene expression (Ryan et al., 2012) as well as by the fact that mRNA-synthesis inhibitors are only effective in blocking LTP when delivered 4–6 h after stimulation (Vickers et al., 2005). We hypothesized that if the 24 h network was representative of a new homeostatic state brought about by LTP induction, its architecture should display an enhanced stability. While there seems to be an early critical time window of transcription for the induction of late-LTP (Nguyen et al., 1994; Vickers et al., 2005), the functional significance of gene expression at 24 h may be coupled to the activation and coordination of pathways related to growth and/or neurogenesis (Ryan et al., 2012). The predominant downregulation of gene expression observed at 24 h may be partly responsible for reorganizing the transcriptional layout toward homeostasis.

Using RBN modeling, we found that the networks derived in the early time points (20 min, 5 h) by the IPA software were more labile, while the most significant network derived at 24 h was markedly more stable (see Figure 1). Furthermore, the WGCNA showed that the degree of co-expression at the different times evidenced a contrasting distribution of the connectivity. The unstimulated control and the 24 h networks fit to a scale-free distribution fairly well while on the contrary, the co-expression networks corresponding to the 20 min and 5 h datasets displayed a truncated distribution. As higher robustness is expected of scale free distributions (Albert et al., 2000), this observation is consistent with the presence of two different homeostatic states before and after LTP induction, whereas transient topological rearrangements are characteristic of intermediate networks.

Interestingly, this temporal effect on the vulnerability of the networks is mirrored by what is known about the vulnerability of LTP and memory itself. Previous studies have shown that LTP can be reversed within hours of induction, but then becomes resistant to reversal (e.g., Xu et al., 1998; Manahan-Vaughan et al., 2000; Woo and Nguyen, 2002). It is of particular relevance to our studies that this resistance to reversal is dependent on new protein synthesis. Thus, our new data support the conclusion that the LTP-related gene networks contribute to the stabilization of LTP.

These results reinforce the view by which the architecture of the networks is under a selective pressure. Yet the contribution of the structural properties to the overall robustness of biological circuits remains to be further clarified—this tendency toward the stable regime represents only one mechanism yielding robust behavior and does not rule out other genetic mechanisms (Wagner, 2005).

4.2. Functional Analysis of WGCNA Modules

While LTP is considered the gold standard model for the cellular mechanism underlying long-term memories, it is becoming clear that LTP encompasses a family of different processes by which neurons integrate and process the information to change their synaptic weights. For example, it has been argued recently that in the Schaffer collateral-commissural pathway at least three mechanistically different forms of synaptic plasticity co-exist, all N-methyl-D-aspartate (NMDA) receptor-dependent (Park et al., 2014). These forms can overlap partially in time, and the combination of these processes can increase their functional utility. In turn, a specific form of LTP consists of a number of mechanistically distinct phases that operate at different levels of organization and time scales. This should come as no surprise if we recognize that neurons capable of modifying the synaptic efficacy for long periods of time are to overcome the limitations imposed by protein and mRNA half-life, typically in the scale of minutes or hours. In contrast to the TFs, which act in the nucleus, some effector genes exert their functions in the distal axonal and dendritic extensions. The translocation of mRNA granules to the synaptic terminals, for example, carries a significant temporal lag between transcription and translation (Knowles et al., 1996; Steward and Schuman, 2003). The control of gene expression must act in coordination with the different time constraints posed by the different subcellular destinations. It is reasonable to assume that the genes transcribed following LTP are effectively being translated at different times and in different subcellular loci.

While understanding the genomic component underlying LTP has been the focus of recent research using differential expression analysis (Lee et al., 2005; Park et al., 2006; Ryan et al., 2011, 2012), these methods are more error-prone for genes with a large expression variation than co-expression analysis. Potentially, genes which are not detected by differential expression can be detected by co-expression if they activate other genes which change enough to be detected. Arguably, a pair of genes with a high value of co-expression are likely to be forming complexes, pathways, or participate in the same cellular circuits (Eisen et al., 1998). Using WGCNA we found that the rapid increase in gene expression observed by the differential expression analysis is complemented by the increase in intramodular connectivity. Out of the 58 modules identified by WGCNA, only one exhibits a significant decrease in average TO from the control to the 20 min time point, while a total of 20 show a significant increase. This suggests that the rapid genomic response that follows LTP induction does not only involve a marked up-regulation of gene expression, but also a tight coordination of the components that ultimately allow the transition to a new homeostatic cellular state. The transitional 5 h dataset shows a loss of intramodular connectivity that parallels the onset of a general down-regulation of gene expression, similarly to the 24 h co-expression network. We believe that the early phase following stimulation is critical in the onset of the genomic changes that are known to be essential for late-LTP. A fundamental fraction of the genes that are transcribed rapidly after LTP induction may be of crucial importance at later times (Nguyen et al., 1994).

The functional analysis confirms the central role of the Egr and Homer families in LTP consolidation and maintenance. Changes in transcription are also likely to be driven by the constitutive transcription factor NFKB in a transcription-independent manner. In fact, our results are consistent with its peaks in activity observed in learning paradigms. In agreement with previous studies, we found that the control of gene expression following LTP is, at least to a certain extent, driven by epigenetic changes. Furthermore, it appears that epigenetic control does not work in isolation, but rather in conjunction with other mechanisms (Lubin et al., 2011). Our study identifies Akt (protein kinase B) in the 20 min dataset even though its expression does not change significantly. The PI3K-Akt-mTOR is regulated via lipid signaling and its role in LTP may have been overlooked in previous studies (although see Sanna et al., 2002). Finally, while regulators of membrane composition are common across all the datasets, neuronal morphological changes are and the amplification of the ubiquitin-proteasome pathway are characteristic of later stages.

4.3. Conclusions

We have presented a view of LTP as a biological process in which a transient signal sets a new homeostatic state that is “remembered” by the cellular systems. Central to this process is the regulation by gene expression, in which the central role played by the Egr TFs early after LTP induction was highlighted by differential expression and co-expression analyses. In addition, we found a rapid enrichment in connectivity at 20 min followed by a systematic decrease. This observation provides a potential explanation for the down-regulation of gene expression at 24 h documented by previous studies. From a systems perspective, we have provided evidence that these networks will show less stable architecture, while networks recruited later will exhibit increased stability, consistent with the fact that are more directly related to LTP consolidation. The architecture exhibited by a control and the 24 h LTP co-expression networks fit well to a scale-free distribution, known to be robust against perturbations, whereas the earlier 20 min and 5 h networks showed truncated distributions. Moreover, using the RBN paradigm we have shown that the network derived at 24 h exhibited an enhanced stability when compared to those derived at earlier times post-LTP. This temporal effect on the vulnerability of the networks is mirrored by what is known about the vulnerability of LTP and memory. Taken together, these results suggest that a new homeostatic state is achieved 24 h post-LTP, and defines an integrated view of the genomic response following LTP induction by which the stability of the networks regulated at different times parallel the properties observed at the synapse.

Conflict of Interest Statement

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Supplementary Material

The Supplementary Material for this article can be found online at: http://journal.frontiersin.org/article/10.3389/fnmol.2015.00042

References

Abraham, W. C., and Williams, J. M. (2003). Properties and mechanisms of ltp maintenance. Neuroscientist 9, 463. doi: 10.1177/1073858403259119

PubMed Abstract | CrossRef Full Text | Google Scholar

Albert, R., Jeong, H., and Barabási, A.-L. (2000). Error and attack tolerance of complex networks. Nature 406, 378–382. doi: 10.1038/35019019

PubMed Abstract | CrossRef Full Text | Google Scholar

Aldana, M., Coppersmith, S., and Kadanoff, L. P. (2003). “Boolean dynamics with random couplings,” in Perspectives and Problems in Nolinear Science, eds E. Kaplan, J. E. Marsden, and K. R. Sreenivasan (New York, NY: Springer), 23–89.

Alexa, A., Rahnenführer, J., and Lengauer, T. (2006). Improved scoring of functional groups from gene expression data by decorrelating GO graph structure. Bioinformatics 22, 1600–1607. doi: 10.1093/bioinformatics/btl140

PubMed Abstract | CrossRef Full Text | Google Scholar

Barabási, A.-L., and Albert, R. (1999). Emergence of scaling in random networks. Science 286, 509–512.

PubMed Abstract | Google Scholar

Cole, A. J., Saffen, D. W., Baraban, J. M., and Worley, P. F. (1989). Rapid increase of an immediate early gene messenger rna in hippocampal neurons by synaptic nmda receptor activation. Nature 340, 474–476.

PubMed Abstract | Google Scholar

Derrida, B., and Weisbuch, G. (1986). Evolution of overlaps between configurations in random boolean networks. J. Phys. 47, 1297–1303.

Google Scholar

Eisen, M. B., Spellman, P. T., Brown, P. O., and Botstein, D. (1998). Cluster analysis and display of genome-wide expression patterns. Proc. Natl. Acad. Sci. U.S.A. 95, 14863–14868.

PubMed Abstract | Google Scholar

Ferrell, J. E. Jr., and Machleder, E. M. (1998). The biochemical basis of an all-or-none cell fate switch in Xenopus oocytes. Science 280, 895–898.

PubMed Abstract | Google Scholar

Fox, J. J., and Hill, C. C. (2001). From topology to dynamics in biochemical networks. Chaos 11, 809–815. doi: 10.1063/1.1414882

PubMed Abstract | CrossRef Full Text | Google Scholar

Ghazalpour, A., Doss, S., Zhang, B., Wang, S., Plaisier, C., Castellanos, R., et al. (2006). Integrating genetic and network analysis to characterize genes related to mouse weight. PLoS Genet. 2:e130. doi: 10.1371/journal.pgen.0020130

PubMed Abstract | CrossRef Full Text | Google Scholar

Haber, D. A., Sohn, R. L., Buckler, A. J., Pelletier, J., Call, K. M., and Housman, D. E. (1991). Alternative splicing and genomic structure of the wilms tumor gene wt1. Proc. Natl. Acad. Sci. U.S.A. 88, 9618–9622.

PubMed Abstract | Google Scholar

Håvik, B., Røkke, H., Dagyte, G., Stavrum, A., Bramham, C., and Steen, V. (2007). Synaptic activity-induced global gene expression patterns in the dentate gyrus of adult behaving rats: induction of immunity-linked genes. Neuroscience 148, 925–936. doi: 10.1016/j.neuroscience.2007.07.024

PubMed Abstract | CrossRef Full Text | Google Scholar

Jeong, H., Tombor, B., Albert, R., Oltvai, Z. N., and Barabási, A.-L. (2000). The large-scale organization of metabolic networks. Nature 407, 651–654. doi: 10.1038/35036627

PubMed Abstract | CrossRef Full Text | Google Scholar

Joilin, G., Guévremont, D., Ryan, B., Claudianos, C., Cristino, A. S., Abraham, W. C., et al. (2014). Rapid regulation of microrna following induction of long-term potentiation in vivo. Front. Mol. Neurosci. 7:98. doi: 10.3389/fnmol.2014.00098

PubMed Abstract | CrossRef Full Text | Google Scholar

Kannan, R., Tetali, P., and Vempala, S. (1997). “Simple markov-chain algorithms for generating bipartite graphs and tournaments,” in Proceedings of the Eighth Annual ACM-SIAM Symposium on Discrete Algorithms (Philadelphia, PA: Society for Industrial and Applied Mathematics), 193–200.

Karlsson, F., and Hörnquist, M. (2007). Order or chaos in boolean gene networks depends on the mean fraction of canalizing functions. Physica A 384, 747–757. doi: 10.1016/j.physa.2007.05.050

CrossRef Full Text | Google Scholar

Kato, A., Ozawa, F., Saitoh, Y., Hirai, K., and Inokuchi, K. (1997). Vesl, a gene encoding vasp/ena family related protein, is upregulated during seizure, long-term potentiation and synaptogenesis. FEBS Lett. 412, 183–189.

PubMed Abstract | Google Scholar

Kauffman, S. (1993). The Origins of Order: Self-organization and Selection in Evolution. New York, NY: Oxford University Press.

Kauffman, S., Peterson, C., Samuelsson, B., and Troein, C. (2003). Random boolean network models and the yeast transcriptional network. Proc. Natl. Acad. Sci. U.S.A. 100, 14796. doi: 10.1073/pnas.2036429100

PubMed Abstract | CrossRef Full Text | Google Scholar

Kauffman, S. A. (1969). Metabolic stability and epigenesis in randomly constructed genetic nets. J. Theor. Biol. 22, 437–467.

PubMed Abstract | Google Scholar

Knowles, R. B., Sabry, J. H., Martone, M. E., Deerinck, T. J., Ellisman, M. H., Bassell, G. J., et al. (1996). Translocation of rna granules in living neurons. J. Neurosci. 16, 7812–7820.

PubMed Abstract | Google Scholar

Langfelder, P., and Horvath, S. (2008). WGCNA: an R package for weighted correlation network analysis. BMC Bioinformatics 9:559. doi: 10.1186/1471-2105-9-559

PubMed Abstract | CrossRef Full Text | Google Scholar

Lee, P. R., Cohen, J. E., Becker, K. G., and Fields, R. R. (2005). Gene expression in the conversion of early-phase to late-phase long-term potentiation. Ann. N. Y. Acad. Sci. 1048, 259–271. doi: 10.1196/annals.1342.023

PubMed Abstract | CrossRef Full Text | Google Scholar

Lee, T., Rinaldi, N., Robert, F., Odom, D., Bar-Joseph, Z., Gerber, G., et al. (2002). Transcriptional regulatory networks in Saccharomyces cerevisiae. Sci. Signal. 298, 799. doi: 10.1126/science.1075090

PubMed Abstract | CrossRef Full Text | Google Scholar

Lubin, F. D., Gupta, S., Parrish, R. R., Grissom, N. M., and Davis, R. L. (2011). Epigenetic mechanisms critical contributors to long-term memory formation. Neuroscientist 17, 616–632. doi: 10.1177/1073858410386967

PubMed Abstract | CrossRef Full Text | Google Scholar

Manahan-Vaughan, D., Kulla, A., and Frey, J. U. (2000). Requirement of translation but not transcription for the maintenance of long-term depression in the ca1 region of freely moving rats. J. Neurosci. 20, 8572–8576.

PubMed Abstract | Google Scholar

Mestl, T., Bagley, R., and Glass, L. (1997). Common chaos in arbitrarily complex feedback networks. Phys. Rev. Lett. 79:653.

Google Scholar

Moore, J. H. (2005). A global view of epistasis. Nat. Genet. 37, 13–14. doi: 10.1038/ng0105-13

PubMed Abstract | CrossRef Full Text | Google Scholar

Nguyen, P. V., Abel, T., and Kandel, E. R. (1994). Requirement of a critical period of transcription for induction of a late phase of ltp. Science 265:1104.

PubMed Abstract | Google Scholar

Park, C. S., Gong, R., Stuart, J., and Tang, S. J. (2006). Molecular network and chromosomal clustering of genes involved in synaptic plasticity in the hippocampus. J. Biol. Chem. 281, 30195. doi: 10.1074/jbc.M605876200

PubMed Abstract | CrossRef Full Text | Google Scholar

Park, P., Volianskis, A., Sanderson, T. M., Bortolotto, Z. A., Jane, D. E., Zhuo, M., et al. (2014). Nmda receptor-dependent long-term potentiation comprises a family of temporally overlapping forms of synaptic plasticity that are induced by different patterns of stimulation. Philos. Trans. R. Soc. Lond. B Biol. Sci. 369, 20130131. doi: 10.1098/rstb.2013.0131

PubMed Abstract | CrossRef Full Text | Google Scholar

Pomerening, J. R., Sontag, E. D., and Ferrell, J. E. (2003). Building a cell cycle oscillator: hysteresis and bistability in the activation of cdc2. Nat. Cell Biol. 5, 346–351. doi: 10.1038/ncb954

PubMed Abstract | CrossRef Full Text | Google Scholar

Ravasz, E., Somera, A. L., Mongru, D. A., Oltvai, Z. N., and Barabási, A.-L. (2002). Hierarchical organization of modularity in metabolic networks. Science 297, 1551–1555. doi: 10.1126/science.1073374

PubMed Abstract | CrossRef Full Text | Google Scholar

Richardson, C., Tate, W., Mason, S., Lawlor, P., Dragunow, M., and Abraham, W. (1992). Correlation between the induction of an immediate early gene, zif/268, and long-term potentiation in the dentate gyrus. Brain Res. 580, 147–154.

PubMed Abstract | Google Scholar

Ryan, M. M., Mason-Parker, S. E., Tate, W. P., Abraham, W. C., and Williams, J. M. (2011). Rapidly induced gene networks following induction of long-term potentiation at perforant path synapses in vivo. Hippocampus 21, 541–553. doi: 10.1002/hipo.20770

PubMed Abstract | CrossRef Full Text | Google Scholar

Ryan, M. M., Morris, G. P., Mockett, B. G., Bourne, K., Abraham, W. C., Tate, W. P., et al. (2013). Time-dependent changes in gene expression induced by secreted amyloid precursor protein-alpha in the rat hippocampus. BMC Genomics 14:376. doi: 10.1186/1471-2164-14-376

PubMed Abstract | CrossRef Full Text | Google Scholar

Ryan, M. M., Ryan, B., Kyrke-Smith, M., Logan, B., Tate, W. P., Abraham, W. C., et al. (2012). Temporal profiling of gene networks associated with the late phase of long-term potentiation in vivo. PLoS ONE 7:e40538. doi: 10.1371/journal.pone.0040538

PubMed Abstract | CrossRef Full Text | Google Scholar

Saha, R., Wissink, E., Bailey, E., Zhao, M., Fargo, D., Hwang, J., et al. (2011). Rapid activity-induced transcription of arc and other iegs relies on poised rna polymerase ii. Nat. Neurosci. 14, 848–856. doi: 10.1038/nn.2839

PubMed Abstract | CrossRef Full Text | Google Scholar

Sanjuán, R., Moya, A., and Elena, S. F. (2004). The contribution of epistasis to the architecture of fitness in an rna virus. Proc. Natl. Acad. Sci. U.S.A. 101, 15376–15379. doi: 10.1073/pnas.0404125101

PubMed Abstract | CrossRef Full Text | Google Scholar

Sanna, P. P., Cammalleri, M., Berton, F., Simpson, C., Lutjens, R., Bloom, F. E., et al. (2002). Phosphatidylinositol 3-kinase is required for the expression but not for the induction or the maintenance of long-term potentiation in the hippocampal ca1 region. J. Neurosci. 22, 3359–3365.

PubMed Abstract | Google Scholar

Shmulevich, I., Kauffman, S. A., and Aldana, M. (2005). Eukaryotic cells are dynamically ordered or critical but not chaotic. Proc. Natl. Acad. Sci. U.S.A. 102, 13439–13444. doi: 10.1073/pnas.0506771102

PubMed Abstract | CrossRef Full Text | Google Scholar

Siegal-Gaskins, D., Mejia-Guerra, M. K., Smith, G. D., and Grotewold, E. (2011). Emergence of switch-like behavior in a large family of simple biochemical networks. PLoS Comput. Biol. 7:e1002039. doi: 10.1371/journal.pcbi.1002039

PubMed Abstract | CrossRef Full Text | Google Scholar

Steward, O., and Schuman, E. M. (2003). Compartmentalized synthesis and degradation of proteins in neurons. Neuron 40, 347–359. doi: 10.1016/S0896-6273(03)00635-4

CrossRef Full Text | Google Scholar

Tuğrul, M., and Kabakçıoğlu, A. (2010). Anomalies in the transcriptional regulatory network of the yeast Saccharomyces cerevisiae. J. Theor. Biol. 263, 328–336. doi: 10.1016/j.jtbi.2009.12.008

CrossRef Full Text | Google Scholar

Vickers, C. A., Dickson, K. S., and Wyllie, D. J. (2005). Induction and maintenance of late-phase long-term potentiation in isolated dendrites of rat hippocampal ca1 pyramidal neurones. J. Physiol. 568, 803–813. doi: 10.1113/jphysiol.2005.092924

PubMed Abstract | CrossRef Full Text | Google Scholar

Wagner, A. (2005). Distributed robustness versus redundancy as causes of mutational robustness. Bioessays 27, 176–188. doi: 10.1002/bies.20170

PubMed Abstract | CrossRef Full Text | Google Scholar

Woo, N. H., and Nguyen, P. V. (2002). Silent metaplasticity of the late phase of long-term potentiation requires protein phosphatases. Learn. Mem. 9, 202–213. doi: 10.1101/lm.498402

PubMed Abstract | CrossRef Full Text | Google Scholar

Xu, L., Anwyl, R., and Rowan, M. J. (1998). Spatial exploration induces a persistent reversal of long-term potentiation in rat hippocampus. Nature 394, 891–894.

PubMed Abstract | Google Scholar

Yao, G., Tan, C., West, M., Nevins, J. R., and You, L. (2011). Origin of bistability underlying mammalian cell cycle entry. Mol. Syst. Biol. 7:485. doi: 10.1038/msb.2011.19

PubMed Abstract | CrossRef Full Text | Google Scholar

Zhang, B., Gaiteri, C., Bodea, L.-G., Wang, Z., McElwee, J., Podtelezhnikov, A. A., et al. (2013). Integrated systems approach identifies genetic nodes and networks in late-onset Alzheimer's disease. Cell 153, 707–720. doi: 10.1016/j.cell.2013.03.030

PubMed Abstract | CrossRef Full Text | Google Scholar

Zhang, B., and Horvath, S. (2005). A general framework for weighted gene co-expression network analysis. Stat. Appl. Genet. Mol. Biol. 4:1128. doi: 10.2202/1544-6115.1128

PubMed Abstract | CrossRef Full Text | Google Scholar