Activity sequences have long been recognized as a fundamental constituent of neural processing. Lorente de No suggested that reverberatory activity sequences in small networks could sustain activity (Lorente de No, 1938). Hebb's idea of cell assemblies suggested that ensembles of cells encoded a particular neuronal concept, but also that there was sequential activation within the group of cells forming the assembly (Hebb, 1949).

Many sensory systems process sequential stimuli, and these are typically mapped to ensembles of sequentially active neurons (Bouchard and Brainard, 2016; Broome et al., 2006; Carrillo-Reid et al., 2015). Deeper in the brain, hippocampal place cells represent a higher-order cognitive map of space, yet here too sequences occur when the animal moves through spatial locations represented by the place cells (Wilson and McNaughton, 1994). The hippocampus exhibits other forms of sequential activity in the form of fast replay events (Jadhav et al., 2012; Wilson and McNaughton, 1994) and stimulus-bridging activity that emerges during associative learning (MacDonald et al., 2011) and trace conditioning (Modi et al., 2014). Thus, there are numerous neural correlates both of sequences, and of processing steps that recognize them.

The idea of synfire chains examines conditions for self-sustaining sequential activity to occur in multiple layers of a network (Abeles, 1982). This is non-trivial, as excessive activation can lead to runaway epileptiform activity, whereas insufficient activation causes decay of the activity wave (Kumar et al., 2008; Mehring et al., 2003). Neural networks that can recognize such sequences are well-established. For example, events may be run through a delay line, so that the oldest event is delayed more, the next event less, and so on, so that they all arrive at the neural network at the same time. Thus, the time-sequence is flattened in time and classical attractor networks can recognize patterns in the sequence (Tank and Hopfield, 1987). Time-varying attractor networks have also been shown to be able to implement and recognize sequences (Lee, 2002). Time-invariance is a desirable feature of sequence recognition circuits, since the same order of events may take place at different speeds. Recursive networks using supervised learning, and short-term synaptic plasticity have been shown to be implement time-invariance (Barak and Tsodyks, 2006; Goudar and Buonomano, 2015; Laje and Buonomano, 2013).

While networks can carry out sequence recognition, they make limited use of the rich dynamics of biological neurons. The theory of Hierarchical Temporal Memory builds on the idea of sequence-recognizing neurons and networks as a way to perform complex temporal computation (George and Hawkins, 2009; Hawkins and Ahmad, 2016). Even the ability to recognize simultaneous closely-localized synaptic input has interesting computational implications (Hawkins and Ahmad, 2016), The current study extends this to patterns of input both in time and space.

One of the first biophysical proposals for subcellular sequence recognition was made by Rall, who showed theoretically that synaptic events propagate down the dendritic tree with a small delay. By timing synaptic inputs to coincide with this delay, Rall predicted that ordered input should yield a larger response than reversed input (Rall, 1964). A stronger version of this mechanism was experimentally demonstrated by Branco et al. (2010) who used glutamate uncaging to provide sequential input along a pyramidal neuron dendrite. They showed that NMDA receptor amplification of 'inward' sequences (distal to proximal) gave rise to about 40% larger somatic depolarization than 'outward' sequences, for rapid (~40 ms) sequences. However, many interesting neural sequences occur at slower time-scales.

At least three attributes should converge for single neurons to achieve and report sequence recognition in the noisy context of neural activity. First, the neurons should recognize inputs coming in the correct order in space and time. Second, the selectivity for the correct input over scrambled input and background noise should be strong enough for there to be reliable discrimination. Third, recognized sequences should trigger activity changes either by way of changed firing, or by way of plasticity. In the current study, we use theory and reaction-diffusion modeling to show how the first condition may be achieved, and develop detailed multi-scale models to address the second and third.

As the setting for this study, we considered a network in which ensembles of neurons are active in sequence, for example, place cells in the hippocampus as an animal moves along a linear track (Figure 1). We assumed that a single neuron from each of these ensembles projects onto a given postsynaptic cell, which is the focus of our analysis. The projections are ordered such that they converge onto a succession of spines located on a short stretch of dendrite on this cell, in the same spatial and temporal order as the activation of the ensembles. We ignored all other network context. We first treated this system in the abstract, as sequentially ordered activation of segments of a one-dimensional reaction-diffusion system with abstract chemistry. We asked if the reaction system could discriminate sequential from scrambled input. We then mapped the discrimination mechanism to a signaling pathway modeled as mass-action reaction-diffusion chemistry in a similar cylindrical geometry, but ornamented with dendritic spines. We then tested discrimination when we embedded mass-action chemistry into a morphologically and electrically detailed neuronal model receiving sequential input on a series of synapses, through ligand-gated ion channels. Finally, we asked if ion channel modulation by sequence discrimination chemistry in small dendritic zones, could affect neuronal firing.

Figure 1

Download asset Open asset

Abstract reaction-diffusion systems support sequence recognition

We first analyzed the requirements for chemical reaction-diffusion systems to achieve selectivity for spatially and temporally ordered inputs. In doing these calculations, we utilized 'spherical cow' models of chemistry: highly reduced formulations with just two diffusive state variables A and B, interpreted as molecules undergoing 1-dimensional diffusion (Figure 2A,B). The reaction system received input from a third molecule, designated as Ca (Figure 2A,B). In all cases, molecule B inhibited A.

Figure 2

Download asset Open asset

The abstract models were initially formulated to implement reaction-diffusion wave propagation at a speed that matched the arrival of sequential input. We implemented two models with this property: the FitzHugh-Nagumo form, which is a known stimulus-triggered, oscillating, and propagating wave system (Fitzhugh, 1961; Nagumo et al., 1962); and a bistable reaction-diffusion system, which again is a known substrate for propagating waves (Keener and Sneyd, 1998). To the bistable switch system, we added inhibitory feedback to restore the switch to baseline after a delay. In addition to these models with complex dynamics, we implemented a system with feedback inhibition, and a system with feedforward inhibition. These simpler models were designed to test if inhibitory feedback or feedforward systems, coupled with diffusion, could achieve selectivity in the absence of an amplifying positive feedback process.

Our criterion for sequence recognition was that ordered inputs should elicit high total activity of a signaling molecule A, whereas scrambled inputs should elicit low levels of activity. Total activity here means the sum of activity of A over the duration and spatial extent of the stimulus.

To investigate the temporal responses of these 'reaction' systems, we first delivered impulse (red dot) and step-function (green line) input to non-spatial versions of the models (Figure 2C). Each model exhibited a large impulse response. Two of the models (feedforward inhibition and state-switching) exhibited only a transient response to the sustained input. The negative feedback model, as expected, had a transient strong response followed by a shallow sustained response. The FitzHugh-Nagumo model, again as expected, oscillated.

We then implemented 1-dimensional reaction-diffusion versions of these models (Figure 2, methods). We delivered Ca²⁺ stimuli at five equally spaced points on the reaction system. Each Ca²⁺ stimulus followed a Gaussian time-profile:

(1) ${\displaystyle Ca=exp(-((phase-t{)}^2)/(2widt{h}^2))}$

Here, phase defines the timing of each input, t is time, and width is the duration of the Ca²⁺ pulse.

The only difference between the sequential and scrambled stimuli was their order ([0,1,2,3,4] and [2,1,4,3,0], respectively (methods)).

The results of these simulations are compared in Figure 2D–G. First, we found that the simple feedback inhibition model showed no input sequence-dependent difference in total activity. This can be seen in terms of time-response at the five stimulus points (Figure 2D,E), and also in terms of spatial profile of A at the end of the stimulus (Figure 2F,G, Video 1).

Video 1

Download asset

Then, we observed that the feedforward inhibition model showed a small amount of selectivity (Figure 2D–G). This arose because the A response was diminished at two of the input points when the input was scrambled (Video 2).

Video 2

Download asset

The FitzHugh-Nagumo (FHN) model was strongly selective, with a large buildup of response only when the input was sequential (Figure 2 third column, Video 3). We interpret this as arising when the propagating wave from the FHN equation arrived at successive input points just when the input was also present.

Video 3

Download asset

The switching model was the most selective (Figure 2 last column, Video 4). This was because it built up to a large, sustained response lasting several seconds, provided the successive inputs were in the same position as the diffusively propagating activity of A.

Video 4

Download asset

Thus, this set of simulations showed that several reaction-diffusion like systems were capable of sequence selectivity. Strong selectivity emerged from a supralinear buildup of responses when diffusively propagating activity was aligned with successive inputs, and suppression due to inhibition when the alignment was off.

Reaction systems select for distinct speeds and length-scales of sequential input

We next asked how the different reaction systems were tuned to different spatial and temporal intervals. We devised a scalar metric, Q, of the degree of sequential ordering of stimuli, as a first step. We plotted the ordinal position of the input against the ordinal value of its arrival time. Thus, a perfect sequence would arrive at positions [0,1,2,3,4] at times [0,1,2,3,4]. We used the quantity

(2) ${\displaystyle Q=m{R}^2}$

as the measure of how ordered the sequence was. Here, m = slope and R = regression coefficient of best-fit line. R² provides a measure of linearity of the sequence, and m assigns it a magnitude and sign. With this measure, the sequence [0,1,2,3,4] has Q = 1; [4,3,2,1,0] has Q = −1, and [4,0,2,1,3] has Q = −0.001. Examples of regression plots used to generate Q for different sequences are presented in Figure 3A,B.

Figure 3

Download asset Open asset

We then compared our reaction readout A with the metric Q, to see how well each of our four reaction systems could discriminate sequence order (Figure 3C,E,G,I). Here, we ran simulations with each of the possible permuted input sequences. In each run, we summed the level of A over all five stimulus points, and from the time of the first stimulus till the end of the run, in order to obtain a cumulative estimate of the response. This total A value (Atot) was plotted against the Q metric for each of these permuted sequences (Figure 3C,E,G,I). If the chemical system was highly selective, we expected that only one or two points in the scatter plot should have a high value of Atot, and these should be the points with Q around +1 or −1. Chemical systems that did not exhibit sequence discrimination should have similar values of Atot for all sequences.

The results were as expected from the individual runs from Figure 2. Specifically, the inhibitory feedback system had poor tuning, there was a small amount of selectivity for the feedforward inhibition case, and strong selectivity for the FitzHugh-Nagumo and switching cases. Each of these outcomes had been observed in Figure 2 for just the sequential and a single scrambled stimulus, and in Figure 3 we found that it generalized to the entire set of sequence permutations. Note that the modeled reaction-diffusion system was spatially symmetric and did not distinguish between forward and backward sequences. An analysis of symmetry breaking arising from incorporation of biological detail is out of the scope of the current paper.

We next asked how selective each reaction system became, when the inputs were delivered at different time and space intervals. The time and space intervals were specified as parameters for each simulation run. We devised a metric for reaction network sequence selectivity, based on a comparison of responses to sequential vs. scrambled inputs:

(3) ${\displaystyle Selectivity=(Asequential-mean(Atot))/max(Atot)}$

Provided the system responds more strongly to sequential input than to any other order, selectivity should be in the range from zero (unselective) to one (highly selective). Asequential is Atot for the sequential stimulus, mean(Atot) is the mean of Atot over all permutations of the stimulus order, and max(Atot) is the maximum value of Atot for all permutations of stimulus order.

To compute Selectivity, we carried out simulations of the model for all possible non-repeating patterns of input, and obtained Atot in each case (Figure 3, first column). Selectivity was then estimated as per Equation 3. We repeated these calculations for each point in the matrix of timing and spacing values, to obtain a grid of selectivity values (Figure 3D,F,H,J). As already seen from Figure 2, the feedback inhibition model had low selectivity in all cases. The feedforward inhibition model showed a diagonal band of weak selectivity such that rapid stimuli in close proximity as well as slower stimuli at greater distances were preferred. The FitzHugh-Nagumo case was much more strongly selective, but also had a similar diagonal band. In both these cases we interpret this as there being an intrinsic propagation speed of the chemical activity wave, and if the inputs were timed and spaced accordingly, the response built up. The switching model was strongly selective, and had a large diffuse zone of strong selectivity centred around (2 s, 4 μm).

Thus, each of the three sequence-selective reaction systems had a preferred range of stimulus timing and spacing, but the FHN and switch forms had much higher selectivity.

Sequence speed selectivity scales with reaction rates and diffusion constants

As a further analysis on our abstract models, we asked how sequence speed selectivity scaled with rates and diffusion parameters. This scaling is important because it determines the spatial extent of the sequence selective zones that we propose to exist in the dendritic tree. It is also important as it defines the time-scales of sequential events that may be recognized by such reaction-diffusion mechanisms.

We first tested the most straightforward assumption, that sequence propagation speed was proportional to chemical rate and diffusion constants. To do this, we scaled all rates and the diffusion constants by the same factor, and reduced the stimulus width by the same factor. We ran the same grid of stimulus space and time intervals as above. As expected, scaling all the chemical and diffusion rates also scaled the sequence speed to which the network had the strongest response. As the rates were increased, the best tuning was at shorter time intervals but roughly the same spatial intervals. (Figure 4A–E). However, different stimulus strengths were needed to obtain strong tuning in these runs. We therefore also examined dependence of tuning on stimulus strength, and observed that the tuning zone also depended on stimulus strength (Figure 4F,G). This is because the selectivity occurs when the stimulus is strong enough that the ordered sequence leads to pathway activation, but not so strong that the stimulus overrides the inhibitory reactions and produces indiscriminate activation. To better understand the range of sequence selectivity, we varied each of the parameters of the model one at a time with respect to the reference model. We carried out this 1-dimensional parameter sensitivity sweep for the FHN and switching models (Figure 4—figure supplements 1 and 2). In each case we computed Asequential (blue traces), and compared it to mean(Atot) averaged over 12 input patterns, 11 of which were scrambled (green traces, methods). We also computed sequence selectivity as per Equation 3 (red traces). We found that selectivity was robust to some parameters, such as diffusion constants, but fragile (range ~20%) with respect to some rate constants and stimulus amplitude. These models were intended to be illustrative and to explore the properties of sequence-selective systems. Hence, we did not optimize the models for robustness.

Figure 4 with 2 supplements see all

Download asset Open asset

Thus, by varying rates for reactions and diffusion, we were able to achieve good sequence selectivity over a broad range of time-intervals (0.5–4 s between stimuli) and spatial intervals (2–6 μm between stimuli).

Mass-action reaction-diffusion system based on MAPK feedback discriminates sequences

We next asked if above design principles for sequence selectivity could be applied to a mass-action reaction-diffusion system. We based our mass-action model on a published, reduced model of MAPK feedback and bistability (Bhalla, 2011). The model already exhibited the switch-like turnon behavior of our abstract 'switch' model. The other key aspect of the abstract model was delayed turnoff of activity. There are several known inhibitory feedback turnoff mechanisms for the MAPK pathway (Lake et al., 2016). We therefore added a MAPK-activated protein phosphatase to cause delayed turnoff of the kinase. The other changes we made to the published MAPK model were (a) to remove the synaptic signaling leading to receptor turnover in the spine and (b) to add diffusible CaM as a Ca²⁺ buffer. With the exception of active PKC, Raf, and the ion channels, all molecules were diffusible. The model schematic is shown in Figure 5A, and the full model specification is presented in supplementary material. This model formulation is a semi-quantitative map to the far more detailed and tightly constrained models of this pathway in the literature (Bhalla and Iyengar, 1999; Resat et al., 2003). Ca²⁺ stimuli were delivered in the PSD and diffused into the dendrite (Methods). Model responses to brief (1 s) and square step function (100 s) Ca²⁺ inputs are shown in Figure 5B,C, both measured in the dendrite. As expected, there is a strong switch-like turnon to Ca²⁺ stimuli, followed by delayed inhibition. We ran this model in a cylindrical geometry ornamented with spines at ~1.1 μm intervals (Materials and methods, Figure 5D). We provided sequential and scrambled Ca²⁺ input to 5 of the spines separated by ~3 μm each, and observed good sequence selectivity (Figure 5E,F).

Figure 5

Download asset Open asset

Thus, we showed that a mass-action reaction-diffusion system exhibited sequence selectivity when it had switch-like turnon with delayed turnoff by negative feedback. This configuration had been predicted by the abstract models.

Sequence selectivity works in a biologically detailed multiscale neuronal model with noisy input

We then asked whether biochemical sequence recognition would 'work' in the more complex context of active neurons receiving background activity in a network. We brought in biological detail at the following levels: (a) Model neuron morphology based on anatomical reconstructions (>3000 segments), including >5700 spines spaced at ~1 μm. (b) Voltage-gated ion channels distributed throughout the cell, based on published models. (c) Background glutamatergic Poisson synaptic input at 0.1 Hz on all spines. (d) Background GABAergic synaptic input with an 8 Hz theta modulation on proximal dendritic compartments. Due to the background synaptic input, the model exhibited theta-modulated sub-threshold oscillations, with spiking activity at ~1 Hz (Figure 6B). (e) Stimuli in the form of spike trains arriving on AMPA and NMDA receptors on subsets of spines. (f) Ca²⁺ dynamics following synaptic ion flux, including diffusion between spine and dendrite, and buffering by calmodulin, in all spines and dendrites (Video 5). (g) The reduced MAPK model described above, was distributed throughout the dendritic tree in ~6000 diffusive compartments (Figure 6A). (h) Chemical calculations in the spines were carried out using a stochastic method (Gillespie Stochastic Systems Algorithm, methods) for the runs used to calculate selectivity.

Figure 6 with 1 supplement see all

Download asset Open asset

Video 5

Download asset

We first ran a direct comparison of a sequential input train [0,1,2,3,4] compared with a scrambled train with sequence [4,0,3,1,2]. We picked time and space intervals of (3 s, 4 μm) based on preliminary calculations for tuning. We delivered the input sequence on a basal dendritic branch (Figure 6A). We selected a set of 5 spines, spaced at ~3 μm. Each spine was stimulated with a Poisson spike burst lasting 2 s at 20 Hz on the AMPA and NMDA receptors. With these parameters, the [Ca²⁺] reached ~40 μM in the PSD, ~35 μM in the spine head, and ~1 μM in the region of dendrite immediately below the stimulated spine (e.g., Figure 6C). Similar Ca²⁺ levels were obtained when stimuli were delivered in other regions of the dendrite. The exception was the primary apical dendrite, in which dendritic Ca²⁺ only reached ~0.5 μM. This is an expected outcome of diluting out the Ca²⁺ arriving from the spine, since the diameter of the primary dendrite was ~2 μm as compared to ~1.0 μm in the other stimulus regions. With these parameters, we observed strong selectivity in MAPK activity for sequential stimuli as compared to scrambled in the basal dendrite (Figure 6D, solid lines for sequential vs. dashed for scrambled).

We then asked how this multiscale model responded to a range of time and space intervals in each of the stimulated zones. Because the calculations were expensive, we used a modified version of the selectivity metric (Equation 3) in which we computed outcomes of only 12 stimulus patterns rather than the exhaustive 120 possible permutations (Materials and methods). In each case, we used MAPK-P in place of the molecule ‘A’ for calculating selectivity (Equation 3). For efficiency, we delivered the stimulus patterns simultaneously in the 4 zones of the cell (Figure 6A). Similar tuning was observed in test simulations where stimuli were delivered only in a single zone.

We observed strong sequence selectivity in the basal dendrite zone, for a restricted range of space and time intervals, as expected from the earlier calculations (Figure 6E–G). We repeated the calculations using different random seeds to generate different background and stimulus spiking input to the model neuron. We found that the responses of individual runs, and hence selectivity, were somewhat noisy (Figure 6E–G).

We next investigated whether the response selectivity was manifested in different regions of the cell. We found that there were indeed strongly sequence selective responses in oblique dendrites and distal apical dendritic regions, but these required stronger stimuli (Figure 6H,I; Videos 6 and 7). In these runs, we had to increase the diameters of spine heads and spine shafts by 40% and 20% respectively to obtain sufficient dendritic Ca²⁺ influx, and hence selectivity. The reference spine models had a shaft of 1 µm length x 0.2 µm diameter, and a head of length and diameter 0.5 µm. We did not observe sequence recognition in the primary apical dendrite (Figure 6J). Even the strongest stimuli applied in this location elicited only small Ca²⁺ elevations, with weak downstream MAPK responses, and no pattern selectivity.

Video 6

Download asset

Video 7

Download asset

We performed a limited parameter sensitivity analysis for the detailed cell model, focusing on biophysical and biochemical parameters involved in triggering or responding to Ca²⁺. We chose these parameters since stimulus range had been one of the most sensitive parameters for the abstract models. We individually varied stimulus strength, L-type Ca²⁺ channel density, Ca²⁺ diffusion, and the activation of Raf by Ca²⁺. Each of these runs was performed with distinct random seeds and stochasticity in spike-trains and in chemical kinetics in spines. Even though we had not fine-tuned the model for robustness, we found that the model retained selectivity (at least 50% of best selectivity) over a factor of ~2 or more for all parameters except the Raf affinity for Ca²⁺ (Figure 6—figure supplement 1). The latter parameter could be varied by ±15% about its reference. While this list is not exhaustive, it does suggest that the sequence selectivity of the detailed model is fairly robust, and somewhat more so than the abstract models.

Thus, we showed that intracellular signaling pathways exhibited selectivity for ordered synaptic input, even when considerable morphological and electrical detail and noisy input were incorporated. Selectivity became noisier with this additional detail, and different zones of the cell required different synaptic efficacies to achieve selectivity. With the current parameters the selectivity was on the time-scale of seconds, with successive inputs spaced apart by a few microns.

Sequence-signaling triggered channel modulation may influence cellular firing

There are multiple possible outcomes of local sequential activity-triggered chemical signaling. These include synaptic plasticity, dendrite remodeling, and electrical changes through channel modulation. For the purposes of the current study, we focused on the third, namely channel modulation leading to changes in electrical activity, since this defines how the neuron may report the presence of a sequence to other cells in the neuronal network. Similar brief, strong synaptic input on small sections of distal dendrites has been reported to cause strong changes in activity of striatal spiny neurons by triggering a transition to an ‘up’ state (Plotkin et al., 2011).

We utilized the same morphologically detailed pyramidal neuron model as above, and asked if local channel modulations over the length scale of about 20 μm could alter somatic spiking. We assumed that MAPK activity, or other signaling events downstream of sequence selective chemistry, could perform channel modulation through phosphorylation or triggering of channel insertion. For example, MAPK phosphorylates and modulates KA (Yuan et al., 2002) Kir6.2 (Lin and Chai, 2008), and Na(v) 1.7channels (Stamboulian et al., 2010). MAPK is also implicated in control of trafficking of AMPA receptors (Keifer et al., 2007). We did not simulate these modulation events explicitly since the computations were very lengthy and we wished to perform many repeats to estimate the distribution of the firing rates following channel modulation. Instead we ran just the electrical cell model for an initial settling period of 1 s, and then used the simulation script to modify the selected channel conductance to represent its modulation. Following this we ran the model for another 3 s in each case. We applied these channel modulations to the same four regions of the cell which we had tested for sequence discrimination. We examined modulation of KA, Na, nonselective leak, NMDAR and AMPAR channels.

We estimated the statistics of firing following these modulatory changes through repeated simulations of the control condition (100 repeats) and each modulation condition (40 repeats), each with different random number seeds. We found that modulation of the Na channel, and increase in nonselective leak conductance in the apical and primary dendrites could increase firing rates by large factors (Figure 7A,B).

Figure 7

Download asset Open asset

In initial runs, modulation of NMDA and AMPA receptors had no effect. This was because the random background synaptic input to these receptors was low, at 0.1 Hz. Since the receptors were activated at this low rate, their effect on cell firing was minimal. However, in these calculations, we had omitted the synaptic input that triggered the chemical cascades. Based on the time-course of build-up of MAPK-P (Figure 5E, Figure 6D), the MAPK-P activity was already high by the time of the last synaptic burst in the sequence of inputs. We therefore repeated the calculations, and additionally incorporated a single burst of synaptic input to the NMDA and AMPA receptors for the same duration and frequency as in the full multiscale cell model (Figure 6).

With this change, we found that a 10x increase in AMPAR conductance in the 20 μm zone in the apical dendrite led to a doubling in cell firing (Figure 7C). A similar manipulation in the primary dendrite led to a shallower increase, about 1.5x.

Thus, chemical signaling, even in rather small regions of the apical dendrite, may lead to a rapid change in cell firing by local channel modulation. However, large channel modulations were needed to elicit sufficiently large (e.g., 50%) changes in in firing rate that could be detected over the noisy background.

We postulate that recursive sequence recognition is a fundamental computational operation in the brain, and that this operation may be implemented by neurons able to recognize spatially and temporally sequential input on behavioral time-scales (seconds). In this study, we have examined three key parts of such a hypothesis. (1) At the chemical level, we have examined a set of abstract rate equations coupled by linear diffusion, and identified key motifs that support sequence recognition. (2) At the dendrite level, we have shown that these abstract principles of reaction-diffusion mediated sequence recognition carry over to much more detailed and biologically motivated neuronal models. These models are multiscale composite models, and explore the length-scales and physiological detail relevant to the current hypothesis. (3) At the cell level, we have shown that local, 20 μm channel modulation can elicit changes in somatic spiking. 20 μm is within the range of our predicted sequence-activated signaling that could lead to channel modulation by phosphorylation or insertion.

We envision a network configuration where at the input level, there are ensembles of neurons that are active in a defined sequence in some sensory, motor or other context. We focus our analysis on a specific postsynaptic neuron, which receives a projection from each of the preceding ensembles, such that each synapse is located within a few μm of the next.

All modeling was carried out using MOOSE, the Multiscale Object-Oriented Simulation Environment (Ray and Bhalla, 2008). MOOSE is freely available, open source, and licensed under the GNU Public License version 3. It can be downloaded from moose.ncbs.res.in and GitHub https://github.com/BhallaLab/moose (Bhalla et al., 2016). MOOSE utilizes the GNU Scientific Library (GSL) Runge-Kutta-Fehlberg fifth order method for chemical computations, and a custom-written branching 1-D diffusion solver using the backward Euler (implicit) method. Stochastic chemical calculations were carried out using a custom-written optimized version of the Gillespie Stochastic Systems Algorithm. Pseudo-random numbers were generated by the Mersenne twister (Matsumoto and Nishimura, 1998). Electrical computations utilized a custom version of the branched nerve equation solution methods described by Hines (Hines, 1984). This has been validated by comparison with other simulators (Gleeson et al., 2010). Interfaces between chemical and electrical signaling components of the model utilized adaptor classes in MOOSE, which average over spatial and temporal discretization differences between the two methods. For example, electrical calculations (yielding Ca²⁺ values) utilize a much smaller timestep (~50 μs; 5 μm) but typically employ a larger spatial step than chemical calculations (~1 ms, 1 μm). The adaptors synchronized the chemical and electrical models every millisecond using first-order corrections to each model system. Electrical model morphologies were either geometrical (cylindrical) or derived from published cell reconstructions available on NeuroMorpho.org (Ascoli et al., 2007; Dougherty et al., 2012). Spines were positioned along dendrites using Poisson statistics with a specified mean spacing between spines (1 μm in the full cell model). Spines were modeled as a cylindrical head (0.5 μm length and diameter) on a cylindrical shaft (1 μm length, 0.2 μm diameter). In some simulations the spine dimensions were scaled up by 20% or 40%. Channel kinetics were derived from Traub et al. (1991). Analysis and plotting was done using Python, NumPy, and MatPlotLib. Simulations were carried out on a variety of Linux workstations, and large calculations were carried out on Linux clusters. Figure generation code for Figures 2, 5 and 6 is available as supplementary material.

Abstract reaction-diffusion calculations (Figures 2–4) were carried out in a 1-dimensional geometry with spatial discretization of 1 μm. In order to represent reactions occurring in a spatially restricted region in or under dendritic spines, the reaction systems were active only in five equally spaced 1 μm patches, but diffusion took place throughout the 100 μm length of the model.

Selectivity calculations were performed according to Equation 3. In this equation, Atot is computed for each permutation of the sequence. This is expensive as there are 120 permutations for a sequence of length 5. Therefore exhaustive permutations were only performed for Figures 3 and 4. In order to generate sequence selectivity matrices for Figure 6 and for all the sensitivity analyses, we took every tenth permutation (including the first one, [0,1,2,3,4]), for a total of 12 permutations. The order of permutations was as generated by the Python library call itertools.permutations. The resultant sequences were: [0,1,2,3,4], [0,2,4,1,3], [0,4,2,1,3], [1,2,0,3,4 , 1,3,4,0,2 , 2,0,3,1,4 , 2,3,0,1,4 , 2,4,3,0,1 , 3,1,2,0,4 , 3,4,0,1,2 , 4,0,3,1,2 , 4,2,1,0,3]. Here [0,1,2,3,4] means that synapse 0 fires first, then synapse 1, then synapse two and so on. In Figure 6—figure supplement 1 we modified Equation 3 slightly by computing Asequential as the mean of 5 runs with the same, sequential stimuli. This was done because the stochastic spiking and chemical calculations introduced considerable variability in the estimate of Asequential. In this figure the mean(Atot) term in Equation 3 took this mean value of Asequential, and Atot for the other 11 permutations, as the set over which the mean was computed.

Chemical models are available in the supplementary material in tabular form and as GENESIS/kkit files, and will be hosted on http://doqcs.ncbs.res.in and ModelDB. Electrical models are available in supplementary material in tabular form and as MOOSE scripts. The morphology file is from NeuroMorpho.org (Ascoli et al., 2007; Dougherty et al., 2012) and attached in the supplementary material as the file VHC-neuron.CNG.swc.

1. 1. Ascoli GA
   2. Donohue DE
   3. Halavi M

   (2007) NeuroMorpho.org

   Publicly available at NeuroMorpho.org (accession no: NMO_09573).

   http://neuromorpho.org/neuron_info.jsp?neuron_name=VHC-neuron

1. Book

   1. Abeles M

   (1982)

   Local Cortical Circuits: An Electrophysiological Study

   Berlin: Springer.

   • Google Scholar
2. 1. Ascoli GA
   2. Donohue DE
   3. Halavi M

   (2007) NeuroMorpho.Org: a central resource for neuronal morphologies

   Journal of Neuroscience 27:9247–9251.

   https://doi.org/10.1523/JNEUROSCI.2055-07.2007

   • PubMed
   • Google Scholar
3. 1. Barak O
   2. Tsodyks M

   (2006) Recognition by variance: learning rules for spatiotemporal patterns

   Neural Computation 18:2343–2358.

   https://doi.org/10.1162/neco.2006.18.10.2343

   • PubMed
   • Google Scholar
4. 1. Bhalla US
   2. Iyengar R

   (1999) Emergent properties of networks of biological signaling pathways

   Science 283:381–387.

   https://doi.org/10.1126/science.283.5400.381

   • PubMed
   • Google Scholar
5. 1. Bhalla US

   (2002) Mechanisms for temporal tuning and filtering by postsynaptic signaling pathways

   Biophysical Journal 83:740–752.

   https://doi.org/10.1016/S0006-3495(02)75205-3

   • PubMed
   • Google Scholar
6. 1. Bhalla US

   (2011) Multiscale interactions between chemical and electric signaling in LTP induction, LTP reversal and dendritic excitability

   Neural Networks 24:943–949.

   https://doi.org/10.1016/j.neunet.2011.05.001

   • PubMed
   • Google Scholar
7. Software

   1. Bhalla US
   2. Ray S
   3. HarshaRani GV
   4. Singh D
   5. Dudani N

   (2016) MOOSE simulator, version e5ac4a0

   GitHub.

   https://github.com/BhallaLab/moose
8. 1. Bouchard KE
   2. Brainard MS

   (2016) Auditory-induced neural dynamics in sensory-motor circuitry predict learned temporal and sequential statistics of birdsong

   PNAS 113:9641–9646.

   https://doi.org/10.1073/pnas.1606725113

   • PubMed
   • Google Scholar
9. 1. Branco T
   2. Clark BA
   3. Häusser M

   (2010) Dendritic discrimination of temporal input sequences in cortical neurons

   Science 329:1671–1675.

   https://doi.org/10.1126/science.1189664

   • PubMed
   • Google Scholar
10. 1. Broome BM
    2. Jayaraman V
    3. Laurent G

    (2006) Encoding and decoding of overlapping odor sequences

    Neuron 51:467–482.

    https://doi.org/10.1016/j.neuron.2006.07.018

    • PubMed
    • Google Scholar
11. 1. Carrillo-Reid L
    2. Miller JE
    3. Hamm JP
    4. Jackson J
    5. Yuste R

    (2015) Endogenous sequential cortical activity evoked by visual stimuli

    Journal of Neuroscience 35:8813–8828.

    https://doi.org/10.1523/JNEUROSCI.5214-14.2015

    • PubMed
    • Google Scholar
12. 1. Dougherty KA
    2. Islam T
    3. Johnston D

    (2012) Intrinsic excitability of CA1 pyramidal neurones from the rat dorsal and ventral Hippocampus

    The Journal of Physiology 590:5707–5722.

    https://doi.org/10.1113/jphysiol.2012.242693

    • PubMed
    • Google Scholar
13. 1. Fitzhugh R

    (1961) Impulses and physiological States in Theoretical models of Nerve membrane

    Biophysical Journal 1:445–466.

    https://doi.org/10.1016/S0006-3495(61)86902-6

    • PubMed
    • Google Scholar
14. 1. Frey U
    2. Morris RG

    (1997) Synaptic tagging and long-term potentiation

    Nature 385:533–536.

    https://doi.org/10.1038/385533a0

    • PubMed
    • Google Scholar
15. 1. George D
    2. Hawkins J

    (2009) Towards a mathematical theory of cortical micro-circuits

    PLoS Computational Biology 5:e1000532.

    https://doi.org/10.1371/journal.pcbi.1000532

    • PubMed
    • Google Scholar
16. 1. Gleeson P
    2. Crook S
    3. Cannon RC
    4. Hines ML
    5. Billings GO
    6. Farinella M
    7. Morse TM
    8. Davison AP
    9. Ray S
    10. Bhalla US
    11. Barnes SR
    12. Dimitrova YD
    13. Silver RA

    (2010) NeuroML: a language for describing data driven models of neurons and networks with a high degree of biological detail

    PLoS Computational Biology 6:e1000815.

    https://doi.org/10.1371/journal.pcbi.1000815

    • PubMed
    • Google Scholar
17. 1. Goudar V
    2. Buonomano DV

    (2015) A model of order-selectivity based on dynamic changes in the balance of excitation and inhibition produced by short-term synaptic plasticity

    Journal of Neurophysiology 113:509–523.

    https://doi.org/10.1152/jn.00568.2014

    • PubMed
    • Google Scholar
18. 1. Hagenston AM
    2. Fitzpatrick JS
    3. Yeckel MF

    (2008) MGluR-mediated calcium waves that invade the soma regulate firing in layer V medial prefrontal cortical pyramidal neurons

    Cerebral Cortex 18:407–423.

    https://doi.org/10.1093/cercor/bhm075

    • PubMed
    • Google Scholar
19. 1. Hawkins J
    2. Ahmad S

    (2016) Why neurons have thousands of Synapses, a Theory of sequence memory in Neocortex

    Frontiers in Neural Circuits 10:.

    https://doi.org/10.3389/fncir.2016.00023

    • PubMed
    • Google Scholar
20. Book

    1. Hebb DO

    (1949)

    The Organization of Behavior

    Wiley.

    • Google Scholar
21. 1. Hines M

    (1984) Efficient computation of branched nerve equations

    International Journal of Bio-Medical Computing 15:69–76.

    https://doi.org/10.1016/0020-7101(84)90008-4

    • PubMed
    • Google Scholar
22. 1. Jadhav SP
    2. Kemere C
    3. German PW
    4. Frank LM

    (2012) Awake hippocampal sharp-wave ripples support spatial memory

    Science 336:1454–1458.

    https://doi.org/10.1126/science.1217230

    • PubMed
    • Google Scholar
23. 1. Kapur A
    2. Yeckel M
    3. Johnston D

    (2001) Hippocampal mossy fiber activity evokes Ca2+ release in CA3 pyramidal neurons via a metabotropic glutamate receptor pathway

    Neuroscience 107:59–69.

    https://doi.org/10.1016/S0306-4522(01)00293-7

    • PubMed
    • Google Scholar
24. Book

    1. Keener J
    2. Sneyd J

    (1998)

    Mathematical Physiology

    New York: Springer.

    • Google Scholar
25. 1. Keifer J
    2. Zheng ZQ
    3. Zhu D

    (2007) MAPK signaling pathways mediate AMPA receptor trafficking in an in vitro model of classical conditioning

    Journal of Neurophysiology 97:2067–2074.

    https://doi.org/10.1152/jn.01154.2006

    • PubMed
    • Google Scholar
26. 1. Kim M
    2. Park AJ
    3. Havekes R
    4. Chay A
    5. Guercio LA
    6. Oliveira RF
    7. Abel T
    8. Blackwell KT

    (2011) Colocalization of protein kinase A with adenylyl cyclase enhances protein kinase A activity during induction of long-lasting long-term-potentiation

    PLoS Computational Biology 7:e1002084.

    https://doi.org/10.1371/journal.pcbi.1002084

    • PubMed
    • Google Scholar
27. 1. Kumar A
    2. Rotter S
    3. Aertsen A

    (2008) Conditions for propagating synchronous spiking and asynchronous firing rates in a cortical network model

    Journal of Neuroscience 28:5268–5280.

    https://doi.org/10.1523/JNEUROSCI.2542-07.2008

    • PubMed
    • Google Scholar
28. 1. Laje R
    2. Buonomano DV

    (2013) Robust timing and motor patterns by taming Chaos in recurrent neural networks

    Nature Neuroscience 16:925–933.

    https://doi.org/10.1038/nn.3405

    • PubMed
    • Google Scholar
29. 1. Lake D
    2. Corrêa SA
    3. Müller J

    (2016) Negative feedback regulation of the ERK1/2 MAPK pathway

    Cellular and Molecular Life Sciences 73:4397–4413.

    https://doi.org/10.1007/s00018-016-2297-8

    • PubMed
    • Google Scholar
30. 1. Larkum ME
    2. Watanabe S
    3. Nakamura T
    4. Lasser-Ross N
    5. Ross WN

    (2003) Synaptically activated Ca2+ waves in layer 2/3 and layer 5 rat neocortical pyramidal neurons

    The Journal of Physiology 549:471–488.

    https://doi.org/10.1113/jphysiol.2002.037614

    • PubMed
    • Google Scholar
31. 1. Lee DL

    (2002) Pattern sequence recognition using a time-varying Hopfield network

    IEEE Transactions on Neural Networks 13:330–342.

    https://doi.org/10.1109/72.991419

    • PubMed
    • Google Scholar
32. 1. Lee KF
    2. Soares C
    3. Thivierge JP
    4. Béïque JC

    (2016) Correlated synaptic inputs drive dendritic calcium amplification and Cooperative plasticity during Clustered Synapse Development

    Neuron 89:784–799.

    https://doi.org/10.1016/j.neuron.2016.01.012

    • PubMed
    • Google Scholar
33. 1. Lin YF
    2. Chai Y

    (2008) Functional modulation of the ATP-sensitive potassium channel by extracellular signal-regulated kinase-mediated phosphorylation

    Neuroscience 152:371–380.

    https://doi.org/10.1016/j.neuroscience.2008.01.003

    • PubMed
    • Google Scholar
34. 1. Lisman JE
    2. Zhabotinsky AM

    (2001)

    A model of synaptic memory: a CaMKII/PP1 switch that potentiates transmission by organizing an AMPA receptor anchoring assembly

    Neuron 31:191–201.

    • PubMed
    • Google Scholar
35. 1. Loewenstein Y
    2. Sompolinsky H

    (2003) Temporal integration by calcium dynamics in a model neuron

    Nature Neuroscience 6:961–967.

    https://doi.org/10.1038/nn1109

    • PubMed
    • Google Scholar
36. 1. Lorente de No R

    (1938)

    Analysis of the activity of the chains of internuncial neurons

    Journal of Neurophysiology 1:195–206.

    • Google Scholar
37. 1. MacDonald CJ
    2. Lepage KQ
    3. Eden UT
    4. Eichenbaum H

    (2011) Hippocampal "time cells" bridge the gap in memory for discontiguous events

    Neuron 71:737–749.

    https://doi.org/10.1016/j.neuron.2011.07.012

    • PubMed
    • Google Scholar
38. 1. Matsumoto M
    2. Nishimura T

    (1998) Mersenne twister: a 623-dimensionally equidistributed uniform pseudo-random number generator

    ACM Transactions on Modeling and Computer Simulation 8:3–30.

    https://doi.org/10.1145/272991.272995

    • Google Scholar
39. 1. Mehring C
    2. Hehl U
    3. Kubo M
    4. Diesmann M
    5. Aertsen A

    (2003) Activity dynamics and propagation of synchronous spiking in locally connected random networks

    Biological Cybernetics 88:395–408.

    https://doi.org/10.1007/s00422-002-0384-4

    • PubMed
    • Google Scholar
40. 1. Modi MN
    2. Dhawale AK
    3. Bhalla US

    (2014) CA1 cell activity sequences emerge after reorganization of network correlation structure during associative learning

    eLife 3:e01982.

    https://doi.org/10.7554/eLife.01982

    • PubMed
    • Google Scholar
41. 1. Nagumo J
    2. Arimoto S
    3. Yoshizawa S

    (1962) An active pulse transmission line simulating nerve axon

    Proceedings of the IRE 50:2061–2070.

    https://doi.org/10.1109/JRPROC.1962.288235

    • Google Scholar
42. 1. Narayanan R
    2. Johnston D

    (2010) The h current is a candidate mechanism for regulating the sliding modification threshold in a BCM-like synaptic learning rule

    Journal of Neurophysiology 104:1020–1033.

    https://doi.org/10.1152/jn.01129.2009

    • PubMed
    • Google Scholar
43. 1. Plotkin JL
    2. Day M
    3. Surmeier DJ

    (2011) Synaptically driven state transitions in distal dendrites of striatal spiny neurons

    Nature Neuroscience 14:881–888.

    https://doi.org/10.1038/nn.2848

    • PubMed
    • Google Scholar
44. 1. Plotkin JL
    2. Shen W
    3. Rafalovich I
    4. Sebel LE
    5. Day M
    6. Chan CS
    7. Surmeier DJ

    (2013) Regulation of dendritic calcium release in striatal spiny projection neurons

    Journal of Neurophysiology 110:2325–2336.

    https://doi.org/10.1152/jn.00422.2013

    • PubMed
    • Google Scholar
45. Book

    1. Rall W

    (1964)

    Theoretical significance of dendritic trees for neuronal input-output relations

    In: R.Reiss, editors. Neural Theory and Modeling. Stanford, CA: Stanford Univ. Press. pp. 73–97.

    • Google Scholar
46. 1. Ray S
    2. Bhalla US

    (2008) PyMOOSE: interoperable scripting in Python for MOOSE

    Frontiers in Neuroinformatics 2:6.

    https://doi.org/10.3389/neuro.11.006.2008

    • PubMed
    • Google Scholar
47. 1. Resat H
    2. Ewald JA
    3. Dixon DA
    4. Wiley HS

    (2003) An integrated model of epidermal growth factor receptor trafficking and signal transduction

    Biophysical Journal 85:730–743.

    https://doi.org/10.1016/S0006-3495(03)74516-0

    • PubMed
    • Google Scholar
48. 1. Rosenblum K
    2. Futter M
    3. Voss K
    4. Erent M
    5. Skehel PA
    6. French P
    7. Obosi L
    8. Jones MW
    9. Bliss TV

    (2002)

    The role of extracellular regulated kinases I/II in late-phase long-term potentiation

    Journal of Neuroscience 22:5432–5441.

    • PubMed
    • Google Scholar
49. 1. Ross WN

    (2012) Understanding calcium waves and sparks in central neurons

    Nature Reviews Neuroscience 13:157–168.

    https://doi.org/10.1038/nrn3168

    • PubMed
    • Google Scholar
50. 1. Smolen P
    2. Baxter DA
    3. Byrne JH

    (2006) A model of the roles of essential kinases in the induction and expression of late long-term potentiation

    Biophysical Journal 90:2760–2775.

    https://doi.org/10.1529/biophysj.105.072470

    • PubMed
    • Google Scholar
51. 1. Stamboulian S
    2. Choi JS
    3. Ahn HS
    4. Chang YW
    5. Tyrrell L
    6. Black JA
    7. Waxman SG
    8. Dib-Hajj SD

    (2010) ERK1/2 mitogen-activated protein kinase phosphorylates sodium channel na(v)1.7 and alters its gating properties

    Journal of Neuroscience 30:1637–1647.

    https://doi.org/10.1523/JNEUROSCI.4872-09.2010

    • PubMed
    • Google Scholar
52. 1. Tank DW
    2. Hopfield JJ

    (1987) Neural computation by concentrating information in time

    PNAS 84:1896–1900.

    https://doi.org/10.1073/pnas.84.7.1896

    • PubMed
    • Google Scholar
53. 1. Traub RD
    2. Wong RK
    3. Miles R
    4. Michelson H

    (1991)

    A model of a CA3 hippocampal pyramidal neuron incorporating voltage-clamp data on intrinsic conductances

    Journal of Neurophysiology 66:635–650.

    • PubMed
    • Google Scholar
54. 1. Tsokas P
    2. Ma T
    3. Iyengar R
    4. Landau EM
    5. Blitzer RD

    (2007) Mitogen-activated protein kinase upregulates the dendritic translation machinery in long-term potentiation by controlling the mammalian target of rapamycin pathway

    Journal of Neuroscience 27:5885–5894.

    https://doi.org/10.1523/JNEUROSCI.4548-06.2007

    • PubMed
    • Google Scholar
55. 1. Tyler WJ
    2. Alonso M
    3. Bramham CR
    4. Pozzo-Miller LD

    (2002) From acquisition to consolidation: on the role of brain-derived neurotrophic factor signaling in hippocampal-dependent learning

    Learning & Memory 9:224–237.

    https://doi.org/10.1101/lm.51202

    • PubMed
    • Google Scholar
56. 1. Wilson MA
    2. McNaughton BL

    (1994) Reactivation of hippocampal ensemble memories during sleep

    Science 265:676–679.

    https://doi.org/10.1126/science.8036517

    • PubMed
    • Google Scholar
57. 1. Yuan LL
    2. Adams JP
    3. Swank M
    4. Sweatt JD
    5. Johnston D

    (2002)

    Protein kinase modulation of dendritic K+ channels in Hippocampus involves a mitogen-activated protein kinase pathway

    Journal of Neuroscience 22:4860–4868.

    • PubMed
    • Google Scholar

Author details

1. Upinder Singh Bhalla

   National Centre for Biological Sciences, Tata Institute of Fundamental Research, Bangalore, India

   Contribution

   USB, Conceptualization, Resources, Data curation, Software, Formal analysis, Supervision, Funding acquisition, Validation, Investigation, Visualization, Methodology, Writing—original draft, Project administration, Writing—review and editing

   For correspondence

   bhalla@ncbs.res.in

   Competing interests

   USB: Reviewing editor, eLife

   "This ORCID iD identifies the author of this article:" 0000-0003-1722-5188

• 2,339

  Page views

• 427

  Downloads

• 19

  Citations

Article citation count generated by polling the highest count across the following sources: Crossref, PubMed Central, Scopus.