Calcium regulation of HCN channels supports persistent activity in a multiscale model of neocortex

Neuronal persistent activity has been primarily assessed in terms of electrical mechanisms, without attention to the complex array of molecular events that also control cell excitability. We developed a multiscale neocortical model proceeding from the molecular to the network level to assess the contributions of calcium (Ca(2+)) regulation of hyperpolarization-activated cyclic nucleotide-gated (HCN) channels in providing additional and complementary support of continuing activation in the network. The network contained 776 compartmental neurons arranged in the cortical layers, connected using synapses containing AMPA/NMDA/GABAA/GABAB receptors. Metabotropic glutamate receptors (mGluR) produced inositol triphosphate (IP3) which caused the release of Ca(2+) from endoplasmic reticulum (ER) stores, with reuptake by sarco/ER Ca(2+)-ATP-ase pumps (SERCA), and influence on HCN channels. Stimulus-induced depolarization led to Ca(2+) influx via NMDA and voltage-gated Ca(2+) channels (VGCCs). After a delay, mGluR activation led to ER Ca(2+) release via IP3 receptors. These factors increased HCN channel conductance and produced firing lasting for ∼1min. The model displayed inter-scale synergies among synaptic weights, excitation/inhibition balance, firing rates, membrane depolarization, Ca(2+) levels, regulation of HCN channels, and induction of persistent activity. The interaction between inhibition and Ca(2+) at the HCN channel nexus determined a limited range of inhibition strengths for which intracellular Ca(2+) could prepare population-specific persistent activity. Interactions between metabotropic and ionotropic inputs to the neuron demonstrated how multiple pathways could contribute in a complementary manner to persistent activity. Such redundancy and complementarity via multiple pathways is a critical feature of biological systems. Mediation of activation at different time scales, and through different pathways, would be expected to protect against disruption, in this case providing stability for persistent activity.


INTRODUCTION
Computational models of functional neural activity patterns in neuronal networks have traditionally focused primarily on the role of electrical activity in shaping these patterns, neglecting the rich chemical complexity that complements electrical signaling in neurons. Persistent neuronal activity, lasting several seconds, has been proposed to underlie several functions in the central nervous system including short-term working memory (Goldman-Rakic, 1995;Braver et al., 1997;Kane and Engle, 2002) and motor preparatory set (Ames et al., 2014). Additional functions probably also depend on similar mechanisms subserved by the ''UP state," identified originally in sleep and in slice but now also demonstrated in the visual cortex (Cossart et al., 2003) and other brain areas (Major and Tank, 2004;Poskanzer and Yuste, 2011;Oikonomou et al., 2014;Zhou et al., 2015). Computational models of network persistent activity that have been proposed largely rely on continued interactions of neurons maintaining activity in one another through http://dx.doi.org/10.1016/j.neuroscience.2015.12.043 0306-4522/Ó 2015 The Authors. Published by Elsevier Ltd. on behalf of IBRO. This is an open access article under the CC BY license (http://creativecommons.org/licenses/by/4.0/). mutual synaptic activation (Lisman et al., 1998;Wang, 1999;Wang, 2001;Goldman, 2013, 2014).
A separate class of studies has focused on modulation of the single neuron via chemical signaling sequences (Egorov et al., 2002;Loewenstein and Sompolinsky, 2003;Fall et al., 2005;Fall and Rinzel, 2006;Franse´n et al., 2006;Ramakrishnan and Bhalla, 2008;Sidiropoulou and Poirazi, 2012;Honnuraiah and Narayanan, 2013;Ashhad et al., 2015;Tiganj et al., 2015). Many of these have demonstrated that calcium (Ca 2+ ) signaling pathways underlie a large repertoire of cell dynamics that complements electrical signalinginteractions of chemophysiology and electrophysiology. Ca 2+ pathways can enable persistent activity at the single-cell level through effects on hyperpolarizationactivated cyclic nucleotide-gated (HCN) channels (with I h current) (Destexhe et al., 1996;Winograd et al., 2008). This molecular/cellular mechanism has also been proposed as one underpinning for working memory (Thuault et al., 2013). HCN channels are an important control point (Burdakov, 2005) in other respects as well, regulating somatic and dendritic responsivity and having effects on network oscillations (Kocsis and Li, 2004;. These network and molecular views of persistent activity and memory are entirely complementary. Indeed, it would be remarkable if multiple scales and multiple mechanisms were not involved in such a basic phenomenology . In this paper, we combine these observations through a multiscale computer model of the neocortex that ranges in scale from intracellular Ca 2þ dynamics, up through cellular electrochemical coupling and on to network activity. Intracellular species simulated include Ca 2+ , inositol triphosphate (IP 3 ), Ca 2+ buffers, and cAMP.
Both single cells and the network produced persistent activity following excitatory stimulation. Metabotropically, the stimulus triggered an intracellular signaling cascade a b  ) and second-(and third-etc.) messenger modulation (red back-beginning arrows). We distinguish membrane-associated ionotropic and metabotropic receptors and ion channels involved in reaction schemes in red (in reality, it is likely that almost every membrane-bound protein is modulated). External events are represented by yellow lightning bolts -there is no extracellular diffusion; the only extracellular reaction is glutamate binding, unbinding and degradation on mGluR1 after an event. Ca 2+ is shown redundantly in blue -note that there is only one Ca 2+ pool for extracellular, 1 pool for cytoplasmic, and 1 pool for ER. (PLC: phospholipase C, DAg: diacyl-glycerol, cAMP: cyclic adenosine monophosphate; PIP 2 : phosphatidylinositol 4,5-bisphosphate). (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.) which caused a persistent increased firing rate in the individual pyramidal cell. Placing this cell into a network with realistic connectivity produced effects augmenting and altering the persistent activity patterns through the combination of metabotropic and ionotropic recurrent activations, with non-stimulated cells suppressed via feedback inhibition. Network firing rates differed across cortical layers, and oscillation frequency shifted. The network model also demonstrated multiple roles of excitation/inhibition (E/I) balance in maintaining both optimal network activation, and in controlling membrane mechanisms mediated through the voltage sensitivity of I h . Multiple Ca 2+ sources -extracellular and intracellular -played synergistic roles in the induction of persistent activity. Ca 2+ buffer in the model was critical for regulating cytosolic free Ca 2+ , and therefore the buffer level had an inverse relationship with persistent activity. By bridging multiple interacting scales and signaling modalities, electrical and chemical, our model demonstrates novel inter-scale interactions providing complementary routes for induction and maintenance of persistent activity.

EXPERIMENTAL PROCEDURES
Network simulations consisted of 776 reduced compartmental cell models with single compartments for inhibitory cells (I cells) and five compartments for pyramidal cells, arrayed by layer with connectivity taken from experimental results on the motor cortex (Weiler et al., 2008) (Fig. 1a). Parallel-conductance electrophysiological simulation in the pyramidal cells was complemented by chemophysiological simulation focused on Ca 2+ handling, based on a prior study (Neymotin et al., 2015) (Fig. 1b). Simulations were run in the NEURON (version 7.4) simulation environment (Carnevale and Hines, 2006) utilizing the reaction-diffusion (RxD) Python module (McDougal et al., 2013a;McDougal et al., 2013b) and NMODL (Hines and Carnevale, 2000). 16 s of simulation time took $6 min using 82 nodes on a Linux cluster with parallel NEURON, run with a fixed time-step of 0.1 ms. The full model is available on ModelDB (#185858).
We describe the scales of the multiscale model sequentially from smaller to larger in the following sections. 1. The molecular scale of 10-100 microns: molecular interactions in the model primarily revolve around Ca 2+ dynamics and the effects of Ca 2+ on various intracellular properties, most notably on the HCN (I h ) channel. 2. The cellular scale of up to 1 mm: the description considers both cell morphology (anatomy) and electrophysiological properties. Although the primary electrophysiological mediators are ion channels at the molecular scale, they are encapsulated at the higher scale as electronic components. This encapsulation and linking across scales occurs as part of multiphysics simulation -the electrophysiological simulation based on electrical circuit theory is coupled to reaction-diffusion formalisms. 3. The local network scale of up to 1 cm. Description of the network is based on the identification of synapse types. The overall simulation is multialgorithmic as well as multiphysical: event-driven algorithms for synapse and network interact with the ordinary and partial differential equation levels of electrophysiological and reaction-diffusion simulation.

Intracellular molecular scale
Our Ca 2þ dynamics (Fig. 1b), similar to those used in Neymotin et al. (2015), were derived from Wagner et al. (2004), itself a spatial variant of Li and Rinzel (1994). Parameters are shown in Table 1. We modeled a onedimensional RxD system of intracellular neuronal Ca 2+ signaling in all compartments of neocortical pyramidal (PYR) neurons. Within each compartment, we modeled cytosolic and endoplasmic reticulum (ER) subcompartments by using a fractional volume for each: for a given cell volume, f ER denotes the fraction occupied by the ER, and f cyt denotes the fraction occupied by the cytosol: f cyt þ f ER ¼ 1.
IP 3 was produced through a reaction sequence initiated by glutamate binding to the metabotropic glutamate receptor (mGluR), based on a reaction scheme developed by Ashhad and Narayanan (2013) (ModelDB #150551). IP 3 diffused outward from the synapse location and decayed following first-order kinetics with s IP 3 of 1 s. Baseline mGluR synaptic weight was normalized to represent the increase in the amount of glutamate bound to mGluR. Extracellular glutamate did not diffuse but was represented by a local Glu value that was incremented in response to an event delivered due to a presynaptic spike. Glu showed bind/unbind kinetics on mGluR and was eliminated by first-order degradation (lower left of Fig. 1b).
The ER Ca 2+ model involves IP 3 receptors (IP 3 Rs), sarco/ER Ca 2+ -ATP-ase pumps (SERCA), and a Ca 2+ leak. IP 3 R dynamics involved a slow Ca 2+ inactivation binding site state (De Young and Keizer, 1992;Li and Rinzel, 1994). This IP 3 R bell-shaped dependence on Ca 2+ is evident in all mammalian isoforms of IP 3 Rs (Tu et al., 2005a,b;Taylor and Tovey, 2010), and provides  Hartsfield, 2005;Fall and Rinzel, 2006;Peercy, 2008). We denote by J IP3R ; J SERCA , and J leak the mass flux per unit volume due to the IP 3 receptor, SERCA pump, and leak channels, respectively. Dividing the mass flux by the volume fraction gives the change in concentration.
Cytosolic Ca 2+ concentration (Ca 2þ cyt ), ER Ca 2+ concentration (Ca 2þ ER ), and IP 3 concentration (IP 3 ) followed: J memb denotes the net flux into the cytosol due to channels on the plasma membrane. Fluxes from the IP 3 R (J IP3R ), SERCA pump (J SERCA ), and leak channels (J leak ) follow: Here, N ¼ N A =10 18 % 6:02 Á 10 5 , the number of molecules in a cubic micron at a concentration of 1 mM, where N A is Avogadro's constant. J IP3R has three gating parameters. Activation gates m Ca and m IP 3 are represented by their steady state values m Ca1 and m IP 3 1 due to their rapid activation rates (Li and Rinzel, 1994). Inactivation via h IP3R is a slow Ca 2+ -dependent process with s IP3R of 400 ms.
The SERCA pump is a pump rather than a channel and so is modeled with Hill-type dynamics. Initial conditions for Ca 2þ cyt ; Ca 2þ ER ; IP 3 , and h IP3R were set to 0.0001 mM, 1.25 mM, 0.0 mM, and 0.8, respectively. Ca 2þ ER was set to 1.25 mM to reflect primed ER Ca 2+ stores. This 1.25 mM is consistent with experimentally-obtained values (Bygrave and Benedetti, 1996;Neymotin et al., 2015).
Calcium buffering followed where B is a diffusible buffer with diffusion coefficients D = 0.043 lm 2 /ms for both B and CaB, about half the rate of Ca 2+ diffusion (Table 1) (Anwar et al., 2012). Ca 2+ extrusion across the plasma membrane was modeled by firstorder decay with s ex ¼ 5 ms except where varied in Fig. 13.
Chemical dynamics are typically slower than electrical dynamics but vary more quickly with position (less stiff in time but more stiff in space). We therefore increased chemical discretization to make sure that our baseline compartmentalization was adequate. Testing Ca 2+ dynamics using 1 lm subcompartments demonstrated that the Ca 2+ values were well-approximated (Fig. 2).
Steady-state x 1 and time constant s x are either related to channel opening aðVÞ and closing kinetics bðVÞ as x 1 ¼ a=ða þ bÞ; s x ¼ 1=ða þ bÞ, or are directly parameterized: x 1 ðVÞ; s x ðVÞ. Kinetics for channels were scaled by Q 10 from an experimental temperature (where available) to a simulation temperature of 37°Celsius. Q 10 ¼ 3 was used when no experimental value was available. All cells contained leak current, transient sodium current I Na , and delayed rectifier current I KÀDR , to allow for action potential generation. Additionally, PYR cells contained in all compartments I KÀA ; I KÀM providing firing-rate adaptation (McCormick et al., 1993). Pyramidal cells also had I h , voltage-gated Ca 2+ channels (VGCCs) in all compartments: I L ; I T ; I N (Kay and Wong, 1987;McCormick and Huguenard, 1992a;Safiulina et al., 2010;Neymotin et al., 2015), and SK and BK Ca 2+ -activated potassium currents (I KCa ). LTS cells contained I L , non-Ca 2+ -dependent I h , SK, and Ca 2+ decay. LTS I h parameterization used Hodgkin-Huxley parameterization: g ¼ 0:00015 S=cm 2 ; E rev ¼ À40 mV;  Fig. 3), modified to provide the faster voltage-sensitivity time constants found in the cortex (Harnett et al., 2015). The mechanism for Ca 2+ regulation of HCN channels in PYR cells in Winograd et al. (2008) is modeled empirically in order to produce the relationship between cytosolic Ca 2+ levels and I h activation without simulating the details of Ca 2+ effects on adenyl cyclase (see schematic for HCN chan in Fig. 1b). This empiric model provides combined control of conductance by cAMP (represented as a normalized level 'p1') and voltage using the following:  from soma with 325 lm space constant, hence e-fold augmented at 325 microns as described by Kole et al. (2006). Apical dendrite I KÀDR ; I KÀA ; I KÀM density also increased with the same length constant, based on data showing HCN and Kv channel colocalization (Harnett et al., 2015;Harnett et al., 2013).
Connection probabilities p ij in the network were dependent on pre-and post-synaptic type and inversely scaled based on distance (Table 2) x, z plane perpendicular to the y-direction of layering, based on data from rodent motor cortex mapping (Weiler et al., 2008;Kiritani et al., 2012a). Individual neurons were placed randomly with uniform distribution. Weights from excitatory cells (E cells) displayed in Table 2 are for the AMPA synapses, with colocalized NMDA weights at 10% of AMPA weights. This provided for AMPA/NMDA current ratio of approximately 1.2, consistent with experiment (Wang et al., 2008). Since connectivity was random, specific divergence could vary. Synaptic delays were randomized between 1.8 ms and 5 ms with additional delay based on distance. Other parameters were based on the literature where available, as well as on previous models (Neymotin et al., 2011a(Neymotin et al., ,b, 2013Song et al., 2013;Chadderdon et al., 2014).
Background activity was simulated by excitatory and inhibitory synaptic inputs following a Poisson process, sent to all cells (Table 3) from other cortical areas and other inputs (Destexhe et al., 2003a). High background rates were required to maintain activity in the model, due to the small network size (Neymotin et al., 2011a). The synaptic stimulus used to assess persistent activity for both single-cell and network simulations was provided to apical dendrites. For most simulations, a single period of stimulation was applied to 50% of PYR cells across all cortical layers. In Fig. 7, stimulation was provided selectively to a fixed number of cells in individual layers. Our current neocortical model, based on motor cortex, does not have cells in L4 so that L4 inputs only synapse onto apical dendrites of L5 cells, an important part of L4 projections in sensory cortex (Shipp, 2005;Hooks et al., 2013;Suter and Shepherd, 2015).
To quantify the extent to which population firing rates changed during and after the stimulus, we formed multiunit activity (MUA) time-series, which count the number of spikes in each bin for a given population. We used several bin sizes (10 ms, 100 ms, 1000 ms) to display temporal activity at multiple resolutions. To quantify the efficacy of the network's response to a stimulus we used firing-rate distinction (FRD): the average ratio of firing rates of stimulated to nonstimulated E cells using 10 ms bin size over the 5 s following the end of stimulation (measured from 0.5 s after the end of input events, to avoid ongoing effects from NMDAR activation). Average measures of Ca 2+ and HCN channel-related variables were also taken from the same period. For single-cell simulations, FRD was calculated as the ratio of the firing rate after the stimulus to the firing rate before the stimulus. To quantify long-term ensemble dynamics, we used population firing rate vectors formed by setting each entry in the vector as the corresponding cell's spikes in a 1000 ms interval. To quantify similarity of these firingrate vectors over time, we took the pairwise Pearson correlation, forming a firing-rate vector similarity matrix (Kelemen and Fenton, 2010;Neymotin et al., 2010).

RESULTS
Simulation overview Over 1600 simulations, with simulation durations of 10-245 s, were analyzed. Simulations of individual pyramidal neurons and simulations of entire networks of neurons were run. A typical simulation of 16 s of network time took $6 min using 82 nodes on a Linux cluster with parallel NEURON. We ran simulations with initial Ca 2+ concentration in the ER set to 1.25 mM (Bygrave and Benedetti, 1996), to mimic the effects of ER Ca 2+ priming via prior excitatory synaptic stimulation (Ross et al., 2005;Hong and Ross, 2007;Fitzpatrick et al., 2009;Neymotin et al., 2015). This allowed the ER Ca 2+ to contribute to HCN channel regulation after mGluR stimulation.

Persistent activity in a single pyramidal neuron
A single layer 5 pyramidal neuron (E5) receiving random background synaptic activation followed by strong excitatory (AMPAR/NMDAR/mGluR) stimulation displayed a persistent increase in firing rate, via the effects of a set of molecular signaling cascades (Fig. 4). Prior to the stimulus, the E5 cell displayed a low firing rate of 1.6 Hz due to approximately balanced excitatory and inhibitory synaptic currents. After stimulation, the cytosolic Ca 2+ concentration increased due to VGCC and NMDAR Ca 2+ influx ( Fig. 4b). At the same time, the mGluR-activated cascade produced IP 3 (Fig. 4c), causing Ca 2+ efflux from ER to cytosol (Fig. 4d). These multiple sources of cytosolic Ca 2+ augmented cAMP, an HCN channel-regulating second-messenger, over tens of seconds (Fig. 4e). Prior to the stimulation, HCN channel conductance (g h ) in the soma and dendrite were at similar levels ( Fig. 4e), but diverged after stimulation. During the upstroke of individual action potentials, g h in the soma followed spiking rapidly while g h in dendrite did not vary as much due to more moderate alterations in dendritic voltage but showed a persistent activation due to the second messenger effect  (Fig. 4g) briefly decreased (note negative inward current) due to membrane depolarization. Afterward, dendritic I h increased, following g h despite some negative feedback due to depolarization effects on both activation and driving force. The associated dendritic depolarization drove spiking in the soma (Fig. 4h). Elevated dendritic voltage lasted >65 s with somatic action potential firing increasing to 52 Hz, a 32Â increase over baseline, during the 10 s after stimulation. The single cell's firing remained elevated above 10 Hz for $62 s after stimulation. The continued activation of I h depended on an ongoing balance between excitation and inhibition. In the absence of inhibition, I h turned off due to its voltage dependence (hyperpolarization activates; depolarization deactivates).

Persistent activity in network
Adding cell-based I h -dependent persistent activity to a network produced a model with synergistic network-and cell/molecular-effects. Individual cells in the network model demonstrated the same pattern of persistent activity after stimulation seen in the single neuron model (Fig. 5). Activation provided augmented inhibition as well as increased excitation. This ongoing inhibition helped maintain I h effect by preventing excessive depolarization that would deactivate I h . The relative isolation of I h in the dendrites also helped prevent this deactivation in the setting of high ongoing activity. During baseline conditions, the average firing rates of excitatory (E) and inhibitory (I) neurons ranged between 2.5-5.4 and 2.3-9.4 Hz, respectively, with rhythmic population oscillations at a high alpha/low beta frequency of $14 Hz. A brief, single stimulation period applied to a randomly selected 50% of the excitatory neurons across layers was used to represent the activation pattern. The immediate effects were depolarization of the stimulated cells and increase in cytosolic IP 3 . During stimulation, average membrane potential was increased, causing a reduction in AMPAR/ NMDAR driving-force and a brief fall-off of firing rate.
After stimulation was turned off, the firing rate of stimulated E cells remained elevated for $60 s, dropping below 10 Hz after 35 s. Rates of the stimulated pyramidal cells increased from 3.7 to 20.0 Hz while rates of the non-stimulated pyramidal cells decreased from 3.6 to 0.8 Hz. The population rhythm increased from 14 Hz to $20 Hz. This resultant oscillation was at the firing rates of the stimulated cells, and was of lower amplitude, being carried primarily by these cells.
The network connectivity enabled several emergent effects produced via synaptic interactions. The higher activity of stimulated E cells produced increased drive from E to I cells: I rates transiently increased to $20 Hz during stimulation (Fig. 5b gray), and then remained elevated throughout the persistent activity. This dynamic increase in inhibition throughout the network was the cause of depressed firing rates of the non-stimulated E cells (>75% reduction) which allowed the stimulated cells to dominate network dynamics. The increased activation of I cells also shortened the network's response duration, compared to the isolated single cell.
Although the stimuli were identical across cortical layers, there was a larger post-stimulus response in layers 2/3 and 6, compared to the two E populations in layer 5 (L5A and L5B). These differences in response correlated with different ratios of convergent inhibitory inputs to convergent excitatory inputs for the E cells in the different layers (L2: 0.4, L5A: 0.2, L5B: 0.1, L6: 0.4). Higher I?E connectivity would translate into greater hyperpolarization of E cells in a layer, greater activation of I h , and greater activation despite the countervailing direct effect of inhibition.

Molecular cascades within cells within network
The same sequence of molecular activations that led to persistent activity in the simulation of the individual E5 cell ( Fig. 4) was also seen in the E5 cells within the full network (Fig. 6). Stimulation-induced depolarization caused further opening of NMDA channels, as well as the activation of VGCCs, which admitted additional Ca 2+ into the cytosol. On average, E5 cells in the stimulated population showed a transient elevation in the cytosolic Ca 2+ of $0.01 mM (Fig. 6a). mGluR stimulation occurred at the same time as AMPA/NMDA synapse activation and produced IP 3 , with subsequent IP 3 binding to ER IP 3 Rs. Together with the available Ca 2+ , this caused the ER to release additional Ca 2+ from its stores into the cytosol after a short interval, further contributing to the population Ca 2+ increase (Fig. 6b). Upon stimulus-induced Ca 2+ influx, the level of Ca 2+ -buffering proteins (B) bound to Ca 2+ (CaB) increased rapidly for the stimulated population (Fig. 6c), preventing the high levels of free cytosolic Ca 2+ which can cause excitotoxicity Taxin et al., 2014;Seidenstein et al., 2015). After stimulation, the Ca 2+ -bound buffer levels gradually decayed as Ca 2+ was taken up into the ER and expelled via the Ca 2+ extrusion pumps. CaB for the non-stimulated population decayed due to hyperpolarization in nonstimulated cells reducing influx and decreasing intracellular Ca 2+ . Upon stimulation, the Ca 2+ influx into the cytosol provided the preparatory steps for persistent activity, causing cAMP to increase and remain high in the stimulated population (Fig. 6d). During the stimulation itself, HCN channel conductance (g h ) declined briefly due to the immediate effect of stimulus-induced depolarization (Fig. 6e). Thereafter, increased inhibition from feedback interneurons reduced stimulated-cell potential, thereby permitting rapid conductance increase, due to the quick responsivity of conductance to voltage change (Harnett et al., 2015), contrasted with the slower second messenger effects. In this way, HCN channel activation and its role in persistent dynamics was critically dependent on network effects, namely ongoing inhibition from the highly activated E?I pathway, as well as on the intrinsic chemophysiological effects of increased second messenger. Similarly, the non-stimulated population showed a slight increase in g h , caused by hyperpolarization due to augmented inhibitory input. However, in the absence of second messenger, this had little effect on the non-stimulated population.

Targeted stimulation of cortical layers
Cortical layers receive layer-specific inputs from other parts of cortex and from other brain areas (Hooks et al., 2013;Suter and Shepherd, 2015). We tested the effects  of inputs into individual layers, stimulating the same number of cells in each case (Fig. 7). In the baseline simulation, activity in the stimulated layer remained confined to that layer. Duration of activity, measured by time from stimulus until the activated neuron rates dropped below 10 Hz, was about 40 s with either individual or all-layer stimulation (Fig. 7a). This differed considerably from the $60 s seen with stimulation of an isolated L5 pyramidal cell (not shown), due to feedback inhibition in the network. Compared with the other 2 stimulation locations, L5 stimulation provided a different activation pattern with lower rate but with an initial plateau period of 15 s before activity decay. This appeared to be a network consequence of the higher density of interconnectivity among E5 cells compared to other layers. The only prominent cross-layer effect was seen with stimulation of L2/3 ( Fig. 7a left column) which provided strong feedforward inhibition in L5. Activation across all layers largely overcame this downward inhibition (Fig. 7 right column). As with L5 focal stimulation, all-layer activation produced an initial plateau and abbreviated activity duration in L5. In order to assess the influence of excitatory/inhibitory balance (Table 2) (Weiler et al., 2008), we also assessed the network with higher internal excitatory gain (Fig. 7b). This gain increase generally produced longer response durations, except that L5 stimulation paradoxically produced a weaker and shorter response (32 s) despite the increase in excitatory strength compared to the baseline network. Instead of interlayer feedforward inhibition from L2/3 to L5, this network showed feedforward excitation with interlayer activity propagation. Stimulation of all layers now produced more prolonged activation with some spread from stimulated to non-stimulated cells.

Prolonged weak stimulation also activated network
Prolonged weak (10% strength) stimulation also produced increases in the firing rate of the stimulated population with further augmentation with greater duration of stimulation (Fig. 8). Weak stimulation produced a less sharp increase in firing rates than was seen with the stronger stimulation and peak rates were reduced, even with 20 s of stimulation. In this low-strength case, the stimulated cells' Ca 2+ transients were broader due to lengthier stimulation, but reached a lower peak. Poststimulus firing-rates of stimulated neurons increased from 7.8 to 11.2 Hz with stimulus duration. Increased excitatory firing promoted higher firing of interneurons, which also increased slightly with stimulus duration (average from 6.0 to 6.8 Hz). As a result of increasing inhibition, the non-stimulated E neurons decreased their post-stimulus rates with increased duration (average 2.6 to 1.8 Hz). We quantified the degree of activation of the stimulated population by comparing it to the firing rate of the non-stimulated cells. We defined a firing-rate distinction (FRD) measure as the average ratio of firing rates for stimulated versus non-stimulated E cells measured by multiunit activity (MUA) over the first 5 s following stimulation (Fig. 5b). FRD increased nearly linearly over stimulus durations from 1 to 14 s, and then saturated with longer duration stimulation (Fig. 8c). As with the brief, strong stimulation, high firing rates persisted after the stimulation was turned off, albeit at a lower level, and produced suppression of the nonstimulated E neurons (purple line).

Neuronal ensemble dynamics during persistent activity
Stimulation of persistent activity was repeatable with similar overall activity as well as specific activity pattern (Fig. 9). The baseline state ensemble activity could be distinguished from that of the stimulated state by noting the change in activity vectors (Fig. 9b), with different elements and layers showing somewhat different degrees of persistence -note that the strands of activation are of different lengths in different layers. The activated state provided a reproducible firing pattern (Fig. 9c). The strong correlations in the block starting at ðx; yÞ coordinates (185 s, 60 s) in Fig. 9c demonstrates that the activity after the second stimulation starting at 185 s is very similar to the activity after the first  Table 4. Low similarity across pre-stimulation and poststimulation states can be seen in the blue-green blocks -e.g. that between the first persistent-activity state (after 60 s) and the first baseline state (starting at time 0 s) is 0.45 as seen in the block from (60,0) -(80,60).

E/I balance modulates persistent activity
Classically E/I balance is considered as a fundamental attribute of network dynamics, balancing activity between the pathological extremes of epileptiform activity (too much E) and inactivity (too much I). This balance also plays this primary role in the current simulation. However, in addition, inhibition was permissive for stimulationinduced persistent activity due to its role in preventing deactivation of I h through excessive depolarization due to synaptic excitation and due to depolarization from I h itself (E h ) at À30 mV with m 1h inflection at about the same voltage depending on ½Ca þþ i (Fig. 3).
Changing the static synaptic weights from E?I neurons, and from I?E neurons, set the overall dynamic level of inhibition in the network. The dynamics of this feedback inhibition was approximately proportional to the product of E?I and I?E synaptic strengths, due to the E?I increasing I cell firing which then produced negative feedback via the I?E weights. Due to this relationship, inhibition increases from bottom-left to top-right in Fig. 10. The dynamic level of inhibition was reflected in neuronal firing rates (Fig. 10b) and in the overall activation of g h (Fig. 10e). Increasing the synaptic weights along the E?I pathway tended to cause the I firing rates to increase (Fig. 10a). As the E cells activated I cells more strongly, the extra feedback inhibition decreased E firing rates (Fig. 10b). Increasing synaptic weights from I?E cells hyperpolarized the E cells and dampened their activity. When weights from E? I cells were set to zero, the I firing rates were sparse.
The modulation of firing rates and the overall depolarization level of E cells influenced the level of stimulus-induced Ca 2+ transients via VGCCs and voltage-sensitive NMDA Ca 2þ flux (Fig. 10c). Ca 2+ transients in E cells were largest when there was minimal inhibition in the network. Elevated Ca 2+ , as a second messenger, then produced additional intracellular molecular effects. The availability of cytosolic Ca 2+ after a stimulus translated directly into the level of cAMP (Fig. 10d), with values increasing in parallel with cytosolic Ca 2+ . Once cAMP increased, it contributed to regulating HCN channel activation (g h ; Fig. 10e). However, in addition to Ca 2þ control, g h is also voltage-sensitive, being deactivated by depolarization. Activation of g h (Fig. 10e) therefore did not follow the patterns of cAMP (Fig. 10d) but instead nearly mirrored it, with large values at high levels of inhibition.
The ratio of g h for stimulated vs. non-stimulated cells supported stimulus-specific representations due to these molecular/cellular mechanisms (Fig. 10f). This ratio was maximally activated over a relatively narrow crescentshaped region of E?I and I?E gains with moderate overall inhibition. At values above and to the right of the crescent, high inhibition controlled HCN channels through voltage effects so that the Ca 2+ transients no longer provided primary regulation of the HCN channel dynamics of the stimulated population.  FRD followed a similar crescent-shaped pattern, with maximal activation at moderate levels of inhibition (Fig. 10g). At these levels of inhibition, there was sufficient cytosolic Ca 2+ to regulate HCN channels in a population-specific manner, and adequate inhibition to activate HCN channels (Fig. 10f, h). At the lowest E?I synaptic weights interneuron firing rates were sparse (left-most column in Fig. 10a), and when I?E weights were also low, there was not sufficient inhibition/ hyperpolarization to activate HCN channels and produce substantial FRD. Increasing I?E weights at these low E?I levels caused enhanced inhibition within the network, enhanced activation of I h , and increased FRD (rising portion of FRD in Fig. 10h). FRD then peaked at intermediate levels of inhibition. At higher E?I weights, increasing I?E weights caused high feedback-inhibition within the network, dampening the FRD response (decaying portion of Fig. 10h). Fig. 10h shows the discrepancy between molecular and network effects, demonstrating that the molecular effects seen in the single-cell model do not fully explain the persistent activity differences in the network model. Relative activation of g h , decreasing with increasing inhibition in the network (symbol color in Fig. 10h), only partially explained FRD, shown on the y-axis of Fig. 10h. Peak FRD did not occur with lowest inhibition both because low inhibition did not adequately suppress the non-stimulated population and because excessive depolarization in the stimulated cells deactivates HCN (Fig. 10e).
Stimulation of both E cells and I cells rather than E cells alone produced results qualitatively similar to Fig. 10, but with reduced FRD (peaking at 31 vs. 39). Stimulating the I cells directly produced additional inhibition in the network which reduced activation of the stimulated cells and reduced Ca 2+ entry, consequently reducing I h effects.

Internal vs. external Ca 2+ stores
Induction of persistent activity in our model relies on Ca 2+ , which comes to the cytosol from both extracellular stores via NMDARs and VGCCs, and from ER via IP 3 R. We investigated the effects of modulating the relative contribution of external and internal Ca 2+ sources by regulating the strength of stimuli influencing each store (Fig. 11) and by adjusting properties of the neurons (VGCC density, Ca 2+ in ER) that lead to different relative strengths of each Ca 2+ store (Fig. 12).
Varying the strength of ionotropic (AMPAR/NMDAR) vs. metabotropic (mGluR) stimulation regulated FRD (Fig. 11). The AMPAR/NMDAR stimulus strength both directly regulated depolarization/firing produced by stimulation, and also admitted Ca 2+ via both NMDARs and VGCCs, and produced augmentation of both [Ca 2þ ] and activity. Peak [Ca 2þ ] was reduced in the absence of AMPAR/NMDAR stimulus, since Ca 2þ was then sourced from the limited ER stores. With increased AMPAR/ NMDAR stimulus strength Ca 2+ continued to increase, whereas it plateaued with increased mGluR strength (Fig. 11a). The differences in Ca 2+ transients from differential AMPAR/NMDAR activation translated into differences in the FRD (Fig. 11b). There was a nearly linear [Ca 2þ ] increase with AMPAR/NMDAR translated into a sigmoidal FRD while the saturating [Ca 2þ ] response from mGluR produced a more rapid increase in FRD, with both plateauing at an FRD value of $25. In the regime of 0-1 gain of AMPAR/NMDAR, it appeared that stimulus strength could be encoded by the FRD, the differential activation of the two populations of neurons.
Simultaneously increasing the AMPAR/NMDAR and mGluR stimulus strengths had only a small additional impact on Ca 2+ transients, with increases adding sublinearly. Increases in Ca 2+ transients at high stimulus strengths did not translate into substantial increased FRD, due to Ca 2þ buffering and extrusion pumps quickly clearing cytosolic Ca 2þ , and preventing it from impacting HCN channel dynamics.
With low Ca 2+ availability from both sources (VGCCs and ER; Fig. 12a lower left and Fig. 12d), the FRD was low, because of insufficient Ca 2+ to trigger a poststimulus response. High entry from external stores only (high VGCC density) produced an increase in prestimulus firing, a small increase in post-stimulus firing, and weak suppression of non-stimulated neuron activity, resulting in moderate FRD (Fig. 12a lower right and  Fig. 12e). High entry from both stores (Fig. 12a upper  right and Fig. 12c) caused a substantial increase in FRD a b Fig. 11. Persistent activity can be induced via either AMPAR/NMDAR (external stores) or mGLUR (internal store) activation. AMPAR/NMDAR and mGLUR stimulus strength were each varied independently from 0 to 3Â baseline (n ¼ 31), while holding the other stimulus strength constant at 1Â baseline. (a) Average cytosolic calcium in the 5 s after stimulus as a function of stimulus strength of AMPAR/NMDAR circles) and mGLUR (triangles). (b) Ratio of stimulated to non-stimulated population firing rates (FRD) as a function of stimulus strength of AMPAR/NMDAR (circles) and mGLUR (triangles).
due to higher post-stimulus firing rates which led to more suppression of non-stimulated neurons. The highest FRD was seen with low flux from external and low to mid flux from internal stores (Fig. 12a middle left and upper left, Fig. 12b), which resulted in the largest increase in stimulated neuron firing, and the largest suppression of non-stimulated neurons. Overall, FRD increased along the y-axis of Fig. 12a, in proportion to initial ER Ca 2+ levels, because ER Ca 2þ only contributes to cytosolic Ca 2+ due to stimulation. In contrast, FRD tended to have an inverted-U relationship with VGCC density (around y=1.0) because increasing VGCC density allowed post-stimulus Ca 2þ influx, but at high levels produced increases in baseline firing of both populations both before and after the stimulus, which reduced FRD. We also modified the FRD measure for use in single-cell simulations by taking the ratio of post-stimulus firing rate to baseline firing rate; across the same set of variations in VGCC density and initial ER Ca 2+ , qualitatively similar results were found (data not shown).

Free Ca 2+ regulates persistent activity
Modulating the plasma membrane Ca 2+ extrusion pump time constant (s ex ) altered Ca 2+ levels, with downstream effects on the network's ability to display Ca 2+ -dependent persistent activity and stimulus-specific representations (Fig. 13). Both stimulation-induced Ca 2+ transients and post-stimulus Ca 2+ levels increased as the Ca 2+ extrusion pump s ex increased (Fig. 13a). This pattern was evident for both the stimulated (light-blue) and non-stimulated (purple) population, although at any s ex , the stimulated population had a higher Ca 2+ concentration due to stimulation it received. At a fast s ex of 0.5 ms the extrusion pump cleared cytosolic Ca 2+ quickly, with prestimulus Ca 2+ levels for both populations being reduced to low levels. The fast action of the extrusion pump also reduced the stimulation-induced Ca 2+ transient, which was considerably smaller than at the baseline s ex of 5 ms (2.5 vs. 13 lM). With this fast s ex , the Ca 2+ was more quickly cleared after stimulation, leaving only a small amount of Ca 2+ available to impact HCN channels, and producing a relatively low FRD (1.2). With slow Ca 2+ extrusion (s ex ¼ 20 ms), pre-stimulus Ca 2+ levels were rising for both populations prior to stimulation, since the Ca 2+ clearance rate was less than the rate of cytosolic Ca 2+ entry from NMDA/ VGCC/ER leak Ca 2+ channels. Stimulation also produced a large Ca 2+ transient of 50 lM, which decayed more gradually. As a result of the slow clearance rate of Ca 2+ , the Ca 2+ remained elevated in the stimulated population longer, contributing to positive feedback and increasing modulation of the stimulated population's g h . These factors allowed the stimulated population to dominate the network dynamics, creating a large overall FRD (20). As s ex increased, the ratio of Ca 2+ in the two populations showed an initial sharp rise followed by a slower rise (Fig. 13a black points). FRD, which depends on Ca 2+ -regulation of HCN channels, also showed a similar relationship with the Ca 2+ extrusion pump s ex (Fig. 13b).
IP 3 R density increase allowed release of more Ca 2+ from ER after mGluR activation (not shown). Changes between 0 and 3Â baseline modulated Ca 2+ transients (0.003-0.014 mM) and increased FRD over a substantial range (12-28). The value with low IP 3 R density reflected a situation where primary Ca 2þ came in from extracellular sources, less effective in the absence of any internal release as shown in Fig. 12. Increasing IP 3 R demonstrated a synergy between Ca 2þ sources.
Modulating Ca 2+ buffer efficacy (FRateÂ buffer concentration) altered the ability of free Ca 2+ to contribute to persistent activity (Fig. 14). Buffering affects the motility and amount of free Ca 2+ . Increased buffer prevented Ca 2+ from being taken up and stored in the ER. Overall, the amount of Ca 2+ released from the ER upon stimulation was inversely proportional to the buffer efficacy, with 1.1 vs. 1.2 mM Ca 2+ released at high and low buffering efficacy, respectively. However, with low buffering efficacy, the ER Ca 2+ levels were only slightly larger than low buffering efficacy, since ER stores were initialized to 1.25 mM and were only allowed 5 s to refill. Increased buffering efficacy meant more Ca 2+ bound prior to stimulation and less available to be taken into the ER, while at the same time reducing baseline Ca 2+ effects on HCN channels. Increases in buffering efficacy also decreased the Ca 2+ transient with stimulation (Fig. 14) -free Ca 2+ availability was primarily a function of buffer concentration with an effect of FRate primarily at the left margin. HCN channel activation showed a similar response profile (Fig. 14b) but FRD alterations were more patchy since FRD is based on effects on the nonstimulated as well as on the stimulated cells (Fig. 14c).

DISCUSSION
We developed a multiscale neocortical model from molecular to network level by evaluating a model of HCN modulation via a second-messenger signaling cascade at cell and network scales. This model displayed the persistent activity that has been hypothesized to contribute to cognitive, mnemonic and behavioral functions that require the maintenance of state over a period of seconds (Goldman-Rakic, 1995;Cossart et al., 2003;Major and Tank, 2004;Ames et al., 2014). The model confirmed prior simulations that showed how the molecular dynamics of Ca 2+ regulation a b of HCN channels (Fig. 6) could produce persistent electrical activity at the single-cell level (Fig. 4). Embedding these molecular models into cells within a network produced emergent effects based on interactions between excitatory and inhibitory neurons (Figs. 5 and 10), and the specific details of neocortical architecture (Fig. 7).

Scale overlaps in neural systems
To a greater extent than in other bodily organ systems, the nervous system features an interplay of scales which makes encapsulation of one scale for use in a higher scale difficult. A weakness of the widely-used point-neuron model is the encapsulation of activity in a single input-output dynamic which not only omits obvious non-point interactions such as dendrodendritic synapses (Rall and Shepherd, 1968), but also subtle interactions such as that between signal integration in apical dendrite and signal integration in local microcircuit. By contrast, the multiscale modeling approach used here builds interacting scales into a coherent model which allows us to explore these interactions across levels of organization (Sejnowski et al., 1988). In the present simulations, we focused on Ca 2+ and on multiple roles for inhibition, bridging from the molecular to cellular and network scales (Lytton and Sejnowski, 1991;Neymotin et al., 2011c). Synaptic connectivity and Ca 2+ handling intersected at multiple locations, including the interplay between extracellular and intracellular sources of Ca 2+ . Due to these second-messenger roles, inhibition plays a major role in activity maintenance through its influence on preventing depolarization with deactivation of I h , while at the same time contributing to activity suppression. Hence the activity maintenance that is classically attributed to excitatory/inhibitory balance is here also provided within the context of inhibitory tone (Fig. 10). We defined dynamic inhibition as a network balance which itself depends on balance provided by E?I and I?E strengths across multiple cell subpopulations in the neocortex, as well as on the excitability of the I cells. Changes in V memb alters Ca 2+ flux through VGCCs based on both V memb and [Ca 2þ cyt ] via the Goldman-Hodgkin-Katz equation, as well as through NMDARs based on voltage sensitivity. As has been noted in other models and in tissue, synaptic bombardment greatly altered both membrane potential and input impedance with multiple subsequent effects (Bernander et al., 1991;Destexhe et al., 2003b). Overly high inhibition dampened neuronal firing while moderate inhibition was consistent with persistent activity over a wide range of parameter settings.
Calcium plays a pivotal role as a second messenger in many neural signaling pathways (Blackwell, 2013;Evans and Blackwell, 2015). Our model evaluated various scenarios regarding Ca 2þ sourcing, Ca 2þ removal, Ca 2+ binding, and time constants of interaction with target molecules. Calcium in cytosol comes from the extracellular space, modeled here as an inexhaustible source, and from organelles which are exhaustible. Compared to mitochondria, ER is likely to be a more accessible Ca 2+ source across multiple locations including spines -it spreads diffusely throughout the neuron, foliating at key points to provide additional surface area (Berridge, 1998;Shemer et al., 2008).
As part of being exhaustible, intracellular sources need to be replenished. Therefore our model implicitly required priming in order to be set up to produce Ca 2+ -mediated persistent activity. We predict that interventions that reduce intracellular stores or prevent their replenishment would have an adverse effect on these mechanisms. Ca 2+ handling is strongly affected by binding, made more complex by the existence of different binding proteins with different affinities and mobilities, many specific to particular cell types (Schwaller, 2010). In our highly simplified model, we considered only a single diffusible binding protein. Increased binding reduced Ca 2+ availability which led to reduced priming for persistent activity as well as less availability for its direct effects after stimulation.

Predictions
Our model makes the following experimentally testable predictions.

Relative activation of E and I neurons has indirect
effects on excitability via Ca 2+ -mediated secondmessenger cascades as well as direct effects. Testable with stimulation of different subsets of neurons (i.e. E vs. I) using optogenetics (Tonnesen et al., 2009;Suter et al., 2014) while imaging Ca 2+ (and ideally cAMP as well). 2. Stimulus-specific persistent activity is only viable in the presence of moderate inhibition. Excessive inhibition dampens activity and prevents Ca 2+ influx from enabling persistent activity. Disinhibition produces high cytosolic Ca 2+ levels in all neurons, causing all neurons to fire at high rates, and lose stimulus-specific information. Testable in vitro by regulating excitability of localized populations of neurons via injectrode application of muscimol (GABA A R agonist) to increase inhibition or bicuculline/gabazine (GABA A R antagonists) to decrease inhibition. 3. In contrast to attractor network models that display a permanent up-state (Wang, 1999(Wang, , 2001, our model predicts that persistent activity will have a limited duration with an upper bound dictated by the multiple time constants associated with Ca 2+ and cAMP regulation of HCN channels. Testable in vitro (Winograd et al., 2008) with intracellular recordings while providing tonic low-grade electrical stimulation followed by heightened stimulation to induce persistent activity, and measuring the duration of increased firing rates.

Compensation, complementarity and synergy
Nervous systems, and biological systems in general, utilize multiple pathways for critical functions. This multiplicity is the key to the adaptability and robustness of these complex systems. This seeming extravagance of mechanism can be arrayed along a spectrum from full redundancy at the one extreme, to the case where a separate system B with an entirely non-overlapping function can nonetheless be called into service to compensate when a system A is compromised or lost. Full redundancy provides a full backup copy -such redundancy may be present in the case of critical messenger peptides or proteins that might be produced without relevant variation from different genes to forestall disastrous mutation. Our simulations showed some degree of redundancy: activation of persistent activity depended on two stores of Ca 2+ . The internal store is potentially exhaustible; in case of failure, additional Ca 2+ can be drawn from extracellular stores. Systems that feature both compensation and complementarity rather than full redundancy are likely to be more common. Complementarity implies some degree of redundancy, but with major differences that give additional features and advantages. The simulations evinced complementarity in the ability of stronger excitatory mechanisms (Fig. 7b) to provide transfer of activation along the network circuitry, while large internal Ca 2þ stores provide the substrate for activating the I h mechanism. Additionally, the different sources of cytoplasmic Ca 2+ , intra-and extracellular, each introduced different time courses of activation that could be utilized for nuanced learning algorithms based on relative degrees of activation of ionotropic and mGluRs.
Moving further along this spectrum, we can consider to what extent these systems are synergistic, working together in a way that produces advantages that neither system would provide individually. We can infer synergy from the many points of interaction between the molecular and network mechanisms which demonstrate persistent activity. For example, excitation has its effects at two very different time scales associated with the mechanisms at different spatial scales. It has immediate effects on membrane potential to activate sets of cells. This activity can then be prolonged due to the relatively long voltage-dependent currents provided by NMDARs. Still slower and longer activations can then be provided by the molecular effects via IP 3 and interaction with internal stores. These activations work together synergistically, boosting the duration of persistent activity. This is still only the tip of this iceberg since the triggering of these multiple intracellular signaling pathways will produce subsequent interactions with the nucleus and genome.

Network effects
The network showed emergent effects not predictable from the single-neuron simulations, through recurrent network loops that provided both excitation and inhibition. Inhibition dominated in our baseline network, strongly reducing activity in non-stimulated cells. This feedback-driven inhibition provided our simulated network with an off switch to terminate persistent activity earlier than in the single cell. Truncating the activation would also ready a network to more quickly switch internal representations, and also provide energysaving. Differences in wiring among different layers of cortex produced longer persistent activity in L2/3 and L6 compared to L5. These differences are likely to have functional consequences due to the different sources/targets of each layer's inputs/outputs (Hooks et al., 2011). Temporally precise interactions between excitatory and inhibitory neurons also produced alpha (14 Hz) oscillations (Lytton and Sejnowski, 1991;White et al., 2000;Neymotin et al., 2011c), which shifted to beta (20 Hz) upon induction of persistent activity, due to the increased firing rates of stimulated excitatory neurons.
In addition to providing a memory mechanism, persistent activation in the network could enhance input/ output for a column or for a layer in a column in a number of ways. Neurons in an activated layer are ''prepared cells," with enhanced responsiveness to subsequent synaptic inputs due to having depolarized potentials nearer to threshold and having fast time constants that enable them to respond quickly to inputs (Bernander et al., 1991). We found that individual cortical layers could also be turned on into a persistent firing state (Fig. 7). An activated layer would better respond to inputs in order to route information to its distinct cortical and subcortical targets (Hunter et al., 2006;Fukushima et al., 2012;Schroeder and Lakatos, 2012;Kumar et al., 2013).
Similarly the major interlaminar route from L2/3 to L5 could be primed and set up to permit signal flowthrough. In this respect our simulated network displayed suppressive and facilitating output modes (Fig. 7 left column). The switch between modes could be effected by neuromodulation by acetylcholine which can increase excitatory tone in neocortex (Hsieh et al., 2000;Hasselmo and McGaughy, 2004;Giocomo and Hasselmo, 2007;Carcea and Froemke, 2013). Excitatory strength determined the magnitude of L5 activation, the major cortical output layer for both intratelencephalic and descending outputs (Anderson et al., 2010;Kiritani et al., 2012b;Yamawaki and Shepherd, 2015). In suppressive mode (Fig. 7a left), targeted stimulation of L2/3 suppressed L5 via strong top-down inhibition from activated E2?I5L?E5 cells. In this mode, the cortical circuit might monitor signals arriving at L2/3, a major input layer (Petersen and Crochet, 2013) without passing that information on. Increasing E?E synaptic weights shifted the network into a facilitating mode (Fig. 7b left), which activated the strong top-down excitatory projections from L2/3?L5 (Fig. 7b left) (Weiler et al., 2008). The network model also showed that stimulation of all cortical layers together with E?E augmentation, produced an additional degree of facilitation (Fig. 7b right).
The simulated network could also sense stimulus magnitude (Fig. 11) and duration (Fig. 8) in progressive levels of persistent activation providing further gradations of encoding. During tonic low-strength stimulation, the stimulated neurons displayed a gradual rise in firing rates, which may allow the neurons to slowly accumulate stimulus-related evidence. As stimulus-related evidence gathers in competing neuronal assemblies, the assembly with the highest firing rate could dominate cognitive representations (winner-take-all) or multiple assemblies could remain active to varying degrees.
The model demonstrated repeatable persistent activity patterns which could contribute to identifiable features of neural representation and memory function in vivo. Stimulation caused activation of an ensemble, which then represented the pattern in a reproducible ensemble activation pattern across the network (Fig. 9). This dynamic representation could be utilized as a form of short-term memory via ensemble encoding for a particular stimulation pattern (Goldman-Rakic, 1995;Durstewitz et al., 2000;Barbieri et al., 2005;Quiroga and Panzeri, 2009;Ghitza, 2011).

Inter-scale interactions produce persistent activity
Many prior models of persistent activity have focused only on cellular, or only on network, dynamics. Our model is unusual in that we develop a multiscale composite of a subcellular-scale model of calcium-induced calcium release (CICR) (Neymotin et al., 2015), a cellular-scale model of hyperpolarization-activated graded persistent activity based on I h modulation (Winograd et al., 2008), and a network-scale neocortical model (Neymotin et al., 2011a,c;Chadderdon et al., 2014). The multiscale model demonstrates many points of inter-scale interaction that would not be fully captured by a series of simplifications, encapsulations and embeddings, such as the classical approach of the point-neuron approximation with its simple input/output function. Due to the inter-scale interactions, our cellular model showed additional dynamics based on the ER's role in Ca 2+ handling. Inter-scale interactions are perhaps most clearly seen where network effects meet Ca 2+ effects at the nexus of I h : Ca control of HCN channels is preparatory for I h activation through relative hyperpolarization from inhibitory network effects. The network model largely followed the pattern of activity set by the cellular model, but showed additional features due to the complexities of multiple types of inputs, to the competing roles of excitation and inhibition, and to the complexity of inter-laminar wiring.
I h appears to have a major role to play in the dynamics of both cells and networks. In addition to the role in providing persistent activity, it also has large effects in maintaining the resting potential, augmenting somatic sub-threshold resonance, and determining dendritic subthreshold responses (Chen et al., 2001;Accili et al., 2002;Poolos et al., 2002;Santoro and Baram, 2003;Dyhrfjeld-Johnsen et al., 2008Zemankovics et al., 2010). Further complexity at the network levels, not explored in this study, arises from inhomogeneous distributions of HCN isoforms and HCN density in different cell types, differential dendritic expression in pyramidal neurons and colocalization with other ion channels (Bender et al., 2001;Accili et al., 2002;Santoro and Baram, 2003;Aponte et al., 2006). Additionally I h is under subtle second messenger control, likely through more than the one effect used in this model. Multiple functions, multiple types, and multiple routes for modulation make HCN channels potent control points in the circuit.
Our model is consistent with a variety of experiments which have demonstrated that I h contributes to neuronal persistent activity: Ca 2+ regulates HCN channels and firing dynamics in cat sensorimotor cortex neurons (Schwindt et al., 1992); blocking I h pharmacologically with ZD7288 abolishes persistent activity in supra-and infra-granular pyramidal neurons from prefrontal cortex of ferrets and guinea pigs (Winograd et al., 2008), and in mouse prefrontal cortex pyramidal neurons (Thuault et al., 2013); HCN channels dynamically regulate resting membrane potential and cellular excitability Thuault et al., 2013) after neuromodulators increase I h (Wang et al., 2007;Heys and Hasselmo, 2012;Tsuno et al., 2013). However, cAMP reduced persistent activity in primates suggesting that it also plays some countervailing roles (Wang et al., 2007;Arnsten et al., 2010). cAMP controls a number of other channels and other cytosolic processes that were not included in our model, including augmentation of K þ flux (e.g. KCNQ) that would reduce cellular excitability (Arnsten and Jin, 2012). HCN channels and certain voltage-sensitive K þ channels (Kv) are colocalized along the apical dendrites (Harnett et al., 2013) which would allow HCN and K þ channels to interact competitively, with HCN channels contributing to a neuronal excitability dampened by Kv (Poolos et al., 2002;Migliore and Migliore, 2012). Overall effects would then depend on relative sensitivity to cAMP, which could be effected through a level of metamodulation via both intrinsic factors and external neuromodulators.
In the current study, we see how a persistent activity attractor landscape can be formed by complementing the effects of mutual excitation with the multiple complementary effects of excitatory and inhibitory interactions across cortical layers (Marder et al., 1996;Nelson et al., 2003;Papoutsi et al., 2013;Papoutsi et al., 2014). By contrast, Wang's ''bump attractor" network model was based on localized connectivity in a single layer of integrate-and-fire neurons (Compte et al., 2000;Wang, 2001). That model produced persistent activity through reverberant NMDAR excitatory activation. In that model NMDAR activation outlasted the effects of fast GABA A activation, with no slower effects from slower dendritic GABA A or from GABA B inputs. Although adding many features, our model is generally consistent with the Wang model, since it demonstrates the importance of excitatory tone in shaping the network's response to targeted stimulation and eliciting interlayer activity propagation. However, a major difference from the bump attractor model is that our network wiring has cortical layering (y-direction) as the major geometrical feature in addition to minor lateral fall-off in connectivity (x-z plane). Our network therefore did not produce enhanced activation of a highly localized subset of strongly interconnected cells -the localized bump.
We attempted to include major neurobiological features that would allow us to assess the many interscale interactions that make the dynamics of the local neocortical network so complex and so fascinating. We necessarily left out many important features, both due to computational limitations and due to the fact that many of these are not understood in sufficient detail for us to adequately portray. Additionally, we of course left out many unimportant features, a key role of modeling being to eventually determine which are the important and which the unimportant parameters related to a given functional outcome. Additional features that would be valuable to more fully understand persistent activity include I CAN cationic currents (Egorov et al., 2002;Tiganj et al., 2015), persistent sodium channels (Oikonomou et al., 2014;Zhou et al., 2015), and neuromodulatory synaptic receptors such as acetylcholine and dopamine (Sawaguchi and Goldman-Rakic, 1991;Sidiropoulou et al., 2009). A further omitted feature that is believed to be important is the aforementioned distribution of HCN isoforms by cell type and location. An additional unexplored consideration is our limited approximation of Ca 2+ handling due to the use of homogeneous diffusion that does not incorporate the important roles of localized subdomains due to molecular barriers, fixed buffers, and colocalizations of targets. Incorporation of some of these considerations will bring us near to the domain of molecular dynamics, a further scale of exploration that remains to be incorporated at some, perhaps remote, future date.