Memory storage was modeled with an Amari-Hopfield network consisting of 10 binary neurons (± states). Small networks were studied to enable full enumeration of the state space (1,024 states for a 10-neuron vector). The network weights W were constructed using a batch Hebbian learning rule. Given a set of p stored patterns {x₁, x₂, …, x_p}, each of length 10, the weight matrix was computed as:

Where X is the p × 10 matrix whose μth row is the pattern x_μ. The diagonal elements were set to zero (i.e., W_ii = 0) to eliminate self-feedback.

In the experimental paradigm, the network was first “pretrained” with a diverse set of 50 random stimuli. To amplify the influence of these stimuli (which would later gradually be overcome by repeated presentation of the target stimulus), each of these 50 patterns was multiplied by 10, then combined into a 50 × 10 matrix which was used to produce initial network weights via the Hebbian learning rule. This pretraining phase established an initial attractor landscape with a rich set of stored memories with which the target stimulus would compete to attract states to its attractor basin.

The target stimulus was then generated randomly. In the target presentation phase, the network was presented repeatedly with the target stimulus, with the number of target presentations being incremented one at a time from 0 to 1,000. For each target presentation count j ∈ {1, 2, …,1,000}, the complete memory matrix was X was constructed by concatenating the pretraining matrix (of size 50 × 10) and j copies of the target stimulus. The weight matrix was then recomputed using the Hebbian learning rule.

The network dynamics were simulated using synchronous updates:

Where sign (·) returns 1 for nonnegative inputs and −1 for negative inputs. Convergence was determined when the state no longer changed between updates.

At each step during the target presentation phase, the attractor basin size was computed corresponding to the target pattern—that is, the number of initial states (sampled exhaustively) that converged to the target attractor. These basin sizes were stored and used to compute increases (or jumps) in attractor basin size during the target presentation phase.

Simulations were implemented in R (R Core Team, 2013). The experiment was repeated over 100 runs with new random pretraining stimuli and target stimulus being generated for each run.

The basin size trajectories for the target stimulus were aggregated across runs to analyze the distribution of positive jump sizes in the target attractor basin. Jump size for stimulus presentation j was defined as the number of initial states that converged on the target state after stimulus presentation j minus the number of initial states that converged on the target state after stimulus presentation j − 1. Discrete jumps were defined as positive changes in basin size given that most stimulus presentations did not lead to a change in basin size. Excess kurtosis for the empirical distribution of positive jump sizes was estimated using the e1071 package in R. Several candidate distributions were fit to the observed positive jump sizes using maximum likelihood estimation implemented in the R package fitdistrplus. The candidate models included a lognormal distribution, an exponential distribution, a half-normal distribution (centered at 1), and a power law distribution (with lower bound x_min = 1 and with the scaling exponent α estimated via the R package poweRlaw). Akaike information criterion (AIC) was computed for each distribution fit and cumulative distribution function (CDF) on a linear scale and complementary cumulative distribution function (CCDF) on a log scale were plotted for the empirical data and the fitted distributions. Parameter stability estimates for distributions were assessed using nonparametric bootstrapping (resampling 1,000 bootstrap replicates from observed data with replacement).

In order to investigate the structure of memory representations after a basin size jump, the state transition graph immediately after each increase in basin size was examined. After an increase in basin size, the set of new states was defined as states belonging to the target basin after the jump but not belonging to the target basin before the jump. Transitions were based on single updates from one state to the next. Transition graphs were then constructed where each node represented a state of the network and a directed edge connected each state to the state obtained after one synchronized update. This yielded a branching tree structure in which each state’s “parent” was defined as the next state obtained after a synchronous update, and the root node was the target attractor state. To quantify the graph structure of the new states added after an increase in basin size, these new states were sorted into branches, and the size of each branch was recorded. To accomplish this, a set of branch heads was first identified within the new states, defined as new states whose immediate parent was not a new state (i.e., which had belonged to the attractor basin prior to the jump). For each of these branch head states, the number of its descendants were computed (while also counting the branch head state itself), thus quantifying the size of each branch of new states which had been added to the target basin. This procedure partitioned all new states into branches, with the sum of branch sizes being equal to the total number of new states after an increase in basin size. Branch sizes were defined as the number of states that eventually flow into a node in the state-transition graph. They were measured by exhaustively determining the descendent state of each possible network state.”

A network with 10 neurons was pretrained with 50 random patterns, each being multiplied by 10, and then a randomly selected target stimulus was presented 1,000 times, as in the previous experiments. Transition graphs were examined after each increase in target basin size, and branch sizes of new states were recorded. This process was repeated over 100 runs, with new random pretraining stimuli and target stimulus being generated for each run. The distribution of branch sizes for new states was then examined and fit with candidate distributions including a lognormal distribution, an exponential distribution, a half-normal distribution, and a power law distribution, as in the previous experiment.

For comparison purposes, the distribution of branch sizes was also examined for all states across a network (i.e., the total number of descendant states, also counting the state itself, for all states across a network). This was performed for a network with 10 neurons that had been pretrained with 50 random patterns, each being multiplied by 10. For visualization, the graph structure of a target attractor basin immediately before and after a basin size jump was plotted for a smaller network of 8 neurons.

In order to examine basin growth in larger networks, networks with 100 and 1,000 neurons were examined. Because it is not feasible to enumerate all possible starting states in networks of this size, starting states were randomly sampled from subsets of states with differing Hamming distances (i.e., number of flipped bits) k from the target stimulus. Distances k of 1, 2, 3, 5, 10, and 20 were sampled for the 100-neuron network and distances 1, 2, 3, 5, 10, 20, 50, and 100 were sampled for the 1,000-neuron network. For each distance k, 100 states of distance k from the target stimulus were randomly sampled. As in previous experiments, the networks were pretrained with random stimuli (one stimulus multiplied by weight 30 for the 100-neuron network and one stimulus multiplied by a weight of 100 for the 1,000-neuron network). As before, the target stimulus was repeatedly presented to the network (the target was presented 500 times). For each presentation of the target, the proportion of sampled initial states at each distance k which converged on the target (i.e., were in the target stimulus basin) was recorded. The same sampled starting states were used across stimulus presentations to avoid noise caused by resampling starting states.

Simulations were performed to investigate the effect of varying levels of memory interference (pretraining stimulus presentations) and stimulus degradation (distance k of the probe stimulus from the target stimulus) on the distribution of sizes of jumps in sampled initial states which are in the target stimulus basin. Simulations were performed on a network with 100 neurons. In order to examine the effect of pretraining stimulus presentations, k was held constant at 20 and simulations were performed with 1, 2, 3, 4, and 5 pretraining stimuli (each multiplied by a weight of 20). In order to examine the effect of k, the number of pretraining stimuli was held constant at 1 (multiplied by a weight of 20), and simulations were run with k of 5, 10, 15, 20, and 25. For each run, the target stimulus was presented 300 times. For each condition, 100 runs were performed. All jump sizes (with jumps defined as any increase in number of sampled starting states which converged on the target) were recorded, and the coefficient of variation of the distribution of positive jump sizes was calculated for each run. Welch’s ANOVA tests were performed to determine whether stimulus degradation (k) or memory interference (number of pretraining memories) significantly affected the coefficient of variation of jump sizes.

Examination of plots of target attractor basin size with repeated target stimulus presentations revealed that attractor basin size increases in intermittent, abrupt jumps, rather than gradually and smoothly increasing (see Figure 1a for an example plot for one experimental run). The proportion of all stimulus presentations which caused an increase in attractor basin size was 0.0126. See Figure 1b for a plot of average increase in target basin size across experimental runs, which shows a relatively smooth increase with repeated stimulus presentations.

Estimated excess kurtosis for the distribution of all positive jump sizes was 11.9, greater than the expected excess kurtosis for a half-normal distribution (~0.87) and for an exponential distribution (Smith and Kounios, 1996). This suggests that the distribution of positive jump sizes is heavy-tailed.

Based on AIC, the best-fitting distribution for positive jump sizes was the lognormal distribution (mean = 2.58, 95% CI [2.49, 2.67]; sd = 1.56, 95% CI [1.51, 1.61]; AIC = 11176.16), which was a superior fit to the half-normal (sigma = 73.21, 95% CI [65.11, 81.14]; AIC = 12626.58), exponential (rate = 0.026, 95% CI [0.024, 0.029]; AIC = 11,653.53), and power law (alpha = 1.39, 95% CI [1.38, 1.40]; AIC = 11,382.23) distributions (see Figures 1c,d). The lognormal distribution was preferred based on AIC in 999 out of 1,000 bootstrap replicates.

Sizes of branch structures of new states belonging to target attractor basins were computed after increases in basin size. The proportion of branches consisting of a single state was 0.69. Of all branches consisting of two or more states, the estimated excess kurtosis was 50.3, greater than the expected excess kurtosis for a half-normal distribution (~0.87) and for an exponential distribution (Smith and Kounios, 1996). The best fitting distribution was a power law distribution, (AIC = 10501.03), which was a superior fit to the half-normal (AIC = 19,054.53), exponential (AIC = 15,130.96), and lognormal (AIC = 12,508.08) distributions (see Figures 2a,b). The power law fit is further supported by the linear appearance of the empirical distribution on the log–log CCDF plot (with the exception of decreased frequency for the largest branch sizes, which may be explained by truncation of the tail due to the finite size of the network).

For the distribution of branch sizes consisting of two or more states across the entire network, estimated excess kurtosis was 74.0. The best fitting distribution was a power law distribution, (alpha = 1.67, 95% CI [1.67, 1.68]; AIC = 1,389,856), which was a superior fit to the half-normal (sigma = 28.37, 95% CI [27.87, 28.84]; AIC = 1,963,401), exponential (rate = 0.088, 95% CI [0.087, 0.088]; AIC = 1,657,266), and lognormal (mean = 1.49, 95% CI [1.48, 1.49]; sd = 1.25, 95% CI [1.25, 1.26]; AIC = 1,510,063) distributions (see Figures 2c,d). The power law distribution was preferred based on AIC in 1000 out of 1,000 bootstrap replicates.

A plot of the graph structure of a target attractor basin immediately before and after a jump in basin size reveals a branching structure of new states (i.e., states which belong to the target basin after the jump but did not belong to the target basin prior to the jump; see Figure 3).

Similarly to the small networks, large networks displayed evidence of discontinuous jumps in target basin with repeat stimulus presentations. The proportion of sampled initial states at various distances k from the target stimulus which fell within the target basin increased in intermittent, abrupt jumps, rather than gradually and smoothly increasing, in networks with 100 and 1,000 neurons (Figures 4A,B).

A Welch’s ANOVA test revealed that the coefficient of variation of accuracy jump sizes was negatively associated with distance from the target stimulus (i.e., stimulus degradation) (F = 184.14, p < 0.001; Figure 4C). Similarly, the coefficient of variation of accuracy jump sizes was negatively associated with pretraining stimuli (i.e., memory interference) (F = 195.86, p < 0.001; Figure 4D).