The mixed selectivity and overlapping representations for different movements in PPC make it an ideal site for mediating complex behaviors like eye-hand coordination. Indeed, growing evidence suggests that coordinated behaviors, such as eye-hand movements, rely on parietal circuits.  The mixing of various types of signals in single PPC neurons and sectors could be interpreted as a manifestation of coordinated activity. For example, most LIP neurons fire stronger when a combined reach and saccade is planned compared to a saccade alone (Hagan et al., 2012). Neural correlates for single and combined eye- and arm-related movements were reported in several PPC fields (Battaglia-Mayer et al., 2001; Calton et al., 2002; Dickinson et al., 2003), with activity being usually weaker for the non-preferred movement. Moreover, neural responses are modulated by static eye and arm position in PEc, V6A and the caudal inferior parietal lobe (Battaglia-Mayer et al., 2000, 2007; Breveglieri et al., 2012, 2014; Piserchia et al., 2017).

The mixing of signals within PPC may result from the short-range intrinsic connections between adjacent parietal areas (Caminiti et al., 2017). In order to understand the mixed selectivity and how it relates to complex behaviors, simultaneous recordings from multiple PPC areas are necessary. However, very few works have employed this method in PPC (e.g., Cui and Andersen, 2007; Dean et al., 2012). By comparison, increasingly interactions between areas of the frontal and parietal cortex are being studied. Multi-area recordings in primates allow for correlations between the activity across areas to be studied and have complemented non-invasive work using fMRI and MEG.

In electrophysiological studies, the local field potential (LFP) has been instrumental in understanding the relationship in neural activity across brain areas. The LFP is composed of synaptic and spiking activity in the vicinity of the recording electrode (Mitzdorf, 1985), and gives an estimate of the population activity. Like spiking-activity, the LFP power is tuned to saccade and reach direction in LIP and PRR, respectively (Pesaran et al., 2002; Scherberger et al., 2005). Synchrony, or coherence, between the firing rates of individual neurons and the LFP at different frequencies may reflect the processing of different types of information (Fries, 2005). During coordinated eye-hand movements, the beta-band (∼15–30 Hz) LFP activity decreases around movement initiation in both LIP and PRR, and correlates with the reaction times for coordinated reach and saccades (but not for saccades made alone, Dean et al., 2012). Furthermore, LIP neurons with reduced activity during eye-hand movements, compared to saccades, tend to be coherent with the beta-band LFP (Hagan et al., 2012) and their firing rate predicts the reaction times of coordinated eye-hand movements. This suggests that these neurons participate in a neural circuit that orchestrates coordinated eye-hand movements (Dean et al., 2012). Coherent activity across areas may also contribute to the processing of cognitive signals such as decision-making (Hawellek et al., 2016; Wong et al., 2016) and visual attention in PPC (Buschman and Miller, 2007; Saalmann et al., 2007; Gregoriou et al., 2009).

The studies of LFP-firing coherence are limited in their ability to provide causal evidence of the role of the PPC in eye-hand coordination. In this regard, a number of inactivation studies in PPC have provided more direct evidence, with two works reporting effects on limb (but not eye) movements (Hwang et al., 2012; Yttri et al., 2014), whereas another one found disrupted eye-hand correlations after bilateral inactivation (Battaglia-Mayer et al., 2013). Furthermore, unilateral inactivation of LIP combined with fMRI resulted in rapid spatial reorganization in the active hemisphere (Wilke et al., 2012), suggesting that the functions of PPC are likely spread over a wider network that extends over both hemispheres. This could also explain recent evidence showing no effect of unilateral LIP inactivation on decision-making (Katz et al., 2016). Similarly, inactivation of VIP had no effect on behavior in a heading discrimination task (Chen et al., 2016). In humans, fMRI-guided transcranial magnetic stimulation demonstrated a causal role of the anterior portion of the IPS to reaching (Reichenbach et al., 2011). Overall, inactivation evidence should be treated cautiously. More sensitive activity manipulations could be useful to determine how PPC nodes contribute to motor behaviors. The use of sophisticated tools such as optogenetics in primates (Jazayeri et al., 2012; Watakabe et al., 2016; El-Shamayleh et al., 2017) could help overcome current limitations.

## IMPLICATIONS OF MIXED SELECTIVITY IN THE PPC FOR MEDICAL INTERVENTIONS

The diversity of signals within the PPC has sparked great interest to the neuroprosthetic community. For patients suffering from loss of function due to paralysis or amputation of a limb, there can be great difficulty in interacting with people or everyday objects. Brain machine interfaces (BMIs) offer some hope in helping remedy these difficulties. A BMI is a device that can record neural activity from the brain while subjects think about a certain task, and then via a decoder, extract the subject's intentions. These decoded intentions are used to control external devices that can vary from a cursor on a monitor, to an anthromorphic robotic arm and hand, to a functional electrical stimulator to activate paralyzed muscles.

Most commonly, electrodes are implanted in the primary motor and premotor areas while patients use motor imagery to provide the necessary input to these BMIs (Markowitz et al., 2011; Hochberg et al., 2012; Collinger et al., 2013). Devices implanted in the motor areas typically decode the trajectory of an effector. Early studies showed that PPC neurons could be used in conjunction with frontal motor areas to control closed loop BMIs, however it was unclear to what extent the PPC neurons contributed to the efficacy of these devices (Wessberg et al., 2000). In a study that compared offline decoding of hand position and velocity in non-human primates, decoding with PPC neurons was inferior to the decoding performance achieved with primary motor and dorsal premotor cortex (Carmena et al., 2003), possibly indicating that the PPC neurons were not contributing much to the overall control.

However, Musallam et al. (2004) went on to demonstrate that high level movement goal information as well as expected reward values of different targets could be decoded from signals in PRR to control a cursor on the screen during a BMI task. These control signals could be generated in the absence of an actual movement. The goal signals allow an abstraction away from the low-level commands necessary to achieve the wanted action as well as the device that actually enacts the action. These low-level commands can be generated through external optimal control algorithms. Goals for multiple sequential movements are planned in PRR (Baldauf et al., 2008) but not in the superior parietal convexity (Li and Cui, 2013) providing a rich mix of signals.

However, soon after this, trajectory information was successfully decoded from the medial bank of the IPS as well as the dorsal convexity to allow control of a 2-dimensional (2D; Mulliken et al., 2008a,b) as well as 3D (Hauschild et al., 2012) cursor on a screen. Decoding algorithms to incorporate the cognitive neural signals and the trajectory information will also provide increased performance compared to each type of signal alone (Shanechi et al., 2013a,b). These studies primarily focused on decoding of spiking activity, but similar information could be extracted from the LFP (Andersen et al., 2004; Scherberger et al., 2005).

The clinical relevance of the PPC to neural prosthetics was demonstrated in the first human trial of a BMI that utilized neural signals from the PPC (Aflalo et al., 2015). In this study, a tetraplegic patient was implanted with electrode arrays in putative areas 5d/PE and AIP and could successfully control 2D and 3D cursors as well as a robotic limb. Therefore, exploiting the richness of information in the PPC may be an advantageous strategy for developing more efficient BMIs.

# CONCLUDING REMARKS

Despite decades of research, a definitive understanding of how individual brains areas are defined, perform distinct computations, and interact with other brain areas remains elusive. The PPC has proved an ideal test bed for understanding how the underlying neural architecture supports a range of sensory, motor and cognitive functions. Anatomy and physiology provide distinct lines of evidence for characterizing the brain areas of the PPC less as a cluster of finite regions and more as a network of integrated areas that may flexibly form the neural basis for diverse functions. The future of systems neuroscience is in understanding how these brain areas work in concert with one another and how the neural dynamics can be used for powering the next generation of prosthetic devices.

# AUTHOR CONTRIBUTIONS

KH, SB, YW and MH contributed to the preparation, writing and revising of this text.

### FUNDING

We acknowledge the Australian Research Council (DE120102883, DE180100344), National Health and Medical

### REFERENCES

Research Council (1020839, 1082144), H2020-MSCA-734227- PLATYPUS and EU Fellowship FP7-PEOPLE-2011-IOF 300452 (SB) for financial support.

## ACKNOWLEDGMENTS

We thank Marcello Rosa and Patrizia Fattori for helpful discussions and support.

input node to the eye/hand coordination system. J. Neurosci. 31, 1790–1801. doi: 10.1523/jneurosci.4784-10.2011

of cortical neurons in primates. Nature 408, 361–365. doi: 10.1038/350 42582

Zhang, C. Y., Aflalo, T., Revechkis, B., Rosario, E. R., Ouellette, D., Pouratian, N., et al. (2017). Partially mixed selectivity in human posterior parietal association cortex. Neuron 95, 697.e4–708.e4. doi: 10.1016/j.neuron.2017.06.040

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

Copyright © 2019 Hadjidimitrakis, Bakola, Wong and Hagan. This is an open-access article distributed under the terms of the Creative Commons Attribution License (CC BY). The use, distribution or reproduction in other forums is permitted, provided the original author(s) and the copyright owner(s) are credited and that the original publication in this journal is cited, in accordance with accepted academic practice. No use, distribution or reproduction is permitted which does not comply with these terms.

# Toward a Biologically Plausible Model of LGN-V1 Pathways Based on Efficient Coding

Yanbo Lian1 , David B. Grayden1,2, Tatiana Kameneva1,3, Hamish Meffin4,5 \* † and Anthony N. Burkitt 1 \* †

*1 Department of Biomedical Engineering, The University of Melbourne, Melbourne, VIC, Australia, 2 Centre for Neural Engineering, The University of Melbourne, Melbourne, VIC, Australia, 3 Faculty of Science, Engineering and Technology, Swinburne University, Melbourne, VIC, Australia, 4 Department of Optometry and Visual Science, The University of Melbourne, Melbourne, VIC, Australia, 5 National Vision Research Institute, The Australian College of Optometry, Melbourne, VIC, Australia*

Increasing evidence supports the hypothesis that the visual system employs a sparse code to represent visual stimuli, where information is encoded in an efficient way by a small population of cells that respond to sensory input at a given time. This includes simple cells in primary visual cortex (V1), which are defined by their linear spatial integration of visual stimuli. Various models of sparse coding have been proposed to explain physiological phenomena observed in simple cells. However, these models have usually made the simplifying assumption that inputs to simple cells already incorporate linear spatial summation. This overlooks the fact that these inputs are known to have strong non-linearities such the separation of ON and OFF pathways, or separation of excitatory and inhibitory neurons. Consequently these models ignore a range of important experimental phenomena that are related to the emergence of linear spatial summation from non-linear inputs, such as segregation of ON and OFF sub-regions of simple cell receptive fields, the push-pull effect of excitation and inhibition, and phase-reversed cortico-thalamic feedback. Here, we demonstrate that a two-layer model of the visual pathway from the lateral geniculate nucleus to V1 that incorporates these biological constraints on the neural circuits and is based on sparse coding can account for the emergence of these experimental phenomena, diverse shapes of receptive fields and contrast invariance of orientation tuning of simple cells when the model is trained on natural images. The model suggests that sparse coding can be implemented by the V1 simple cells using neural circuits with a simple biologically plausible architecture.

Keywords: efficient coding, LGN-V1 pathways, biological plausibility, separated ON and OFF sub-regions, push-pull effect, phase-reversed feedback, receptive fields, contrast invariance

# 1. INTRODUCTION

In early experimental studies of cat striate cortex, Hubel and Wiesel found two main types of cells: simple cells and complex cells (Hubel and Wiesel, 1959, 1962). Simple cells exhibit linear spatial summation of visual stimuli, while complex cells have significant non-linear behavior. This difference is reflected in receptive field (RF) structures of the two types of cells. Receptive fields (RFs) describe spatial patterns of light and dark regions in the visual field that in combination are effective at driving neural response. They are frequently modeled as linear spatial filters. Simple

### Edited by:

*Greg Stuart, Australian National University, Australia*

### Reviewed by:

*Marco Mainardi, Scuola Normale Superiore di Pisa, Italy C. Daniel Meliza, University of Virginia, United States*

### \*Correspondence:

*Hamish Meffin<EMAIL_ADDRESS>Anthony N. Burkitt<EMAIL_ADDRESS>
*†These authors share joint senior authorship*

Received: *31 October 2018* Accepted: *19 February 2019* Published: *14 March 2019*

### Citation:

*Lian Y, Grayden DB, Kameneva T, Meffin H and Burkitt AN (2019) Toward a Biologically Plausible Model of LGN-V1 Pathways Based on Efficient Coding. Front. Neural Circuits 13:13. doi: 10.3389/fncir.2019.00013* cells have a single RF filter, reflecting the linear spatial summation properties, while complex cells pool the output for two or more RF filters in a non-linear fashion.

Over the past decades, some important characteristics of simple cell RF have been observed experimentally (with emphasis on cat and primates, but also ferrets). First, simple cells show a range of selectivity for the orientation of visual stimuli, from highly oriented RFs, which are selective to an optimal orientation, to non-oriented RFs, which are insensitive to orientation. Many RFs of simple cells in V1 are oriented, localized, and bandpass (Hubel and Wiesel, 1962, 1968), while non-orientated RFs are seen in all layers of V1 (Hawken et al., 1988; Chapman and Stryker, 1993). Second, RFs of orientation tuned simple cells can be well-described by two-dimensional Gabor functions (Jones and Palmer, 1987a; Ringach, 2002). In addition, both these studies found some blob-like RFs, which are broadly tuned in orientation. Third, RFs of simple cells have spatially segregated ON and OFF sub-regions (Hubel and Wiesel, 1962; Martinez et al., 2005); i.e., the spatial region that excites the simple cell in response to bright (ON) stimuli is separated from the region that excites the cell in response to dark (OFF) stimuli (left column of **Figure 1**). Fourth, simple cells show push-pull responses; i.e., if one stimulus excites a simple cell, the stimulus with opposite contrast, but same location, will inhibit the simple cell (Jones and Palmer, 1987b; Ferster, 1988; Hirsch et al., 1998; Martinez et al., 2005). One example of the push-pull effect can be seen on the left of **Figure 1** where a simple cell is excited by input from a cell in the lateral geniculate nucleus (LGN) responding to dark spots (an OFF LGN cell) but is effectively inhibited by LGN cells responding a bright spot in the same location (an ON LGN cell). Fifth, feedback from simple cells to LGN cells frequently has a phase-reversed influence compared to the feedforward input (Wang et al., 2006); i.e., where the RF of an ON (OFF) LGN cell is overlapped with the ON (OFF) sub-region of the RF of a simple cell, i.e., feedforward excitation, feedback from the simple cell to the LGN cell is suppressive; where an ON (OFF) LGN cell coincides with the OFF (ON) sub-region of a simple cell RF, i.e., effective feedforward suppression, the feedback is facilitatory. This effect of phase-reversed feedback is also illustrated in **Figure 1**, where the influence from a simple cell to LGN cells is opposite to the influence from LGN cells to the same simple cell. Lastly, the orientation tuning property of simple cells are contrast invariant; i.e., the shape and width of orientation tuning curves remain the same for different stimulus contrasts (Sclar and Freeman, 1982; Skottun et al., 1987; Finn et al., 2007; Priebe, 2016).

On the other hand, insights from computational modeling of V1 cells have also been used to explain experimental data. Sparse coding has been proposed by many researchers as a principle employed by the brain to process sensory information. Olshausen and Field reproduced localized, oriented and spatially bandpass RFs of simple cells based on a sparse coding model that aimed to reconstruct the input with minimal average activity of neurons (Olshausen and Field, 1996, 1997). However, the original model failed to generate non-oriented RFs observed in experiments (Ringach, 2002). Subsequently, Olshausen and colleagues found that the sparse coding model can produce RFs that lack strong

orientation selectivity by having many more model neurons than the number of input image pixels (Olshausen et al., 2009). Rehn and Sommer introduced hard sparseness to classical sparse coding, which minimizes the number of active neurons rather than the average activity of neurons in the original model, and demonstrated that the modified sparse coding model can generate diverse shapes of simple cell RFs (Rehn and Sommer, 2007). Zhu and Rozell showed that many visual non-classical RF effects of V1 such as end-stopping, contrast invariance of orientation tuning can emerge from a dynamical system based on sparse coding (Zhu and Rozell, 2013).

These studies were important in explaining the RF structure, but made a number of simplifying assumptions that overlooked many details of biological reality, include some or all of the following. First, the responses of neurons (e.g., firing rates) should be non-negative. Second, the learning rule of synaptic connections should be local where the changes of synaptic efficacy depend only on pre-synaptic and post-synaptic responses. Third, the learning rule should not violate Dale's Law, namely that neurons release the same type of transmitter at all their synapses, and consequently, the synapses are either all excitatory or all inhibitory (Strata and Harvey, 1999). Fourth, the computation of the response of any neuron should be local, such that only neurons synaptically connected to this target neuron can be involved. In addition, a biologically plausible model should also be consistent with important experimental evidence. For LGN-V1 visual pathways, experimental evidence includes the existence of a large amount of cortico-thalamic feedback (Swadlow, 1983; Sherman and Guillery, 1996), longrange excitatory but not inhibitory connections between LGN and V1, and separated ON and OFF channels for LGN input (Hubel and Wiesel, 1962; Ferster et al., 1996; Jin et al., 2008, 2011). The original sparse coding model neglects many of the biological constraints described above.

Several recent studies addressed the issue of biological plausibility by incorporating some of these constraints, while continuing to neglect others. For example, Zylberberg and colleagues designed a spiking network (based on sparse coding) that can account for diverse shapes of simple cell RFs using lateral inhibition (Zylberberg et al., 2011). The local learning rule and the use of spiking neurons bring some degree of biological plausibility to the model, but the model employs connections that can change sign during learning, which violates Dale's law, and there are not separate channels for ON and OFF LGN input. Additionally, the effect of sparse coding is achieved by competition between units via lateral inhibition, but a recent study suggested that dominant lateral interactions are excitatory in the mouse cortex (Lee et al., 2016). In another modeling work of simple cell RFs, Wiltschut and Hamker designed an efficient coding model with separated ON and OFF LGN cells, and, feedforward, feedback, and lateral connections that can generate various types of simple cell RFs (Wiltschut and Hamker, 2009), but their model does not incorporate Dale's law.

As with earlier studies (Olshausen and Field, 1996, 1997; Rehn and Sommer, 2007; Olshausen et al., 2009), these more recent studies (Wiltschut and Hamker, 2009; Zylberberg et al., 2011), incorporating biological constraints, have continued to focus on the RF structure of simple cells, while largely neglecting the experimental phenomena shown in **Figure 1**. This is because they have typically not separated inputs from ON and OFF LGN cells, which is a key distinction underlying all the phenomena listed in **Figure 1**. One important question in this regard is how these non-linear (half-wave rectified) LGN inputs are combined to give linear RFs for simple cells and whether this causes the experimental phenomena listed in **Figure 1**. To our knowledge, Jehee and Ballard are the only researchers that have explicitly explained the effect of phase-reversed feedback using a model based on predictive coding (Jehee and Ballard, 2009). However, the RFs generated by their model do not match well with those observed in experiments and the push-pull effect for simple cells has not been explained. In addition, the formula for calculating responses of model neurons (Jehee and Ballard, 2009, Equation 7) is not local and the learning rule neglects Dale's law.

In this paper, we propose a two-layer model of LGN-V1 visual pathways that can account for experimental phenomena:

Our model is biologically plausible by incorporating:

The first layer consists of ON and OFF LGN cells and the second layer consists of simple cells. The connections from the first layer to the second layer (feedforward connections) and from the second layer to the first layer (feedback connections) consist of separate excitatory and inhibitory connections. Even though the inhibitory connections between LGN and V1 should be implemented via intermediate populations of inhibitory interneurons, we use neurons that have both excitatory and inhibitory connections to simplify the circuit. This aspect of the model is not biologically plausible, but possible biologically plausible neural circuits for implementing inhibitory connections are proposed in the Discussion section. The model presented here is relevant to visual cortices both with and without an orientation columnar organization.

The novelty of the model proposed here is that it models LGN-V1 pathways using segregated ON and OFF LGN channels and separate excitatory and inhibitory connections to investigate the structure of connections between LGN and simple cells to explain a wide range of experimental phenomena. In addition, it can generate a wide variety of experimentally observed RFs of simple cells. Also, the model is biologically plausible by respecting many biological constraints and important experimental evidence. Finally, the experimental phenomena explained in this paper are all caused by the structure of learned connections between LGN and V1 after the model is trained on natural image data.

# 2. MATERIALS AND METHODS

## 2.1. Sparse Coding

The original sparse coding model (Olshausen and Field, 1996) proposed that simple cells represent their sensory input in such a way that their spiking rates in response to natural images tend to be statistically independent and rarely attain large values (near the top of the cells' dynamic range). Mathematically this means that the joint distribution of spike rates over natural images is the product of the distributions for individual cells, and that each of these individual distributions has a long tail (i.e., high kurtosis). Additionally it was proposed that the representation should allow the reconstruction of the sensory input through a simple weighted sum of visual features with minimal error. This can be formulated as an optimization problem of minimizing the cost function,

$$E(\mathbf{A}, \mathbf{s}) = \frac{1}{2} \|\mathbf{x} - \mathbf{A}\mathbf{s}\|\_2^2 + \lambda \sum\_{i} Q(s\_i), \tag{1}$$

where **x** represents the input, columns of the matrix **A** represent basis vectors that are universal visual features from which any image can be constructed from a weighted sum, **s** is the vector of responses, si , of model units that represent the corresponding coefficients for all basis vectors, Q(·) represents a penalty function that favors low activity of model units, and λ is a parameter that scales the penalty function (Olshausen and Field, 1996, 1997). The term **As** in Equation (1) is the reconstruction of the input from the model, so the first term on the right-hand-side of Equation (1) represents the sum of squared difference between the input and model reconstruction. The second term on the right-hand-side of Equation (1) tends to push **s** to small values. Therefore, by solving this minimization problem, the model finds a sparse representation for the input. By taking the partial derivatives of Equation (1) in terms of the elements of **A** and **s**, and applying gradient descent, the dynamic equations and the learning rule are given by

$$\begin{aligned} \dot{\mathbf{s}} &= \mathbf{A}^T \mathbf{r} - \lambda Q'(\mathbf{s}) \\ \Delta \mathbf{A} &\propto \langle \mathbf{r} \mathbf{s}^T \rangle, \end{aligned} \tag{2}$$

where **r** = **x** − **As**, h·i is the average operation, the dot notation represents differentiation with regard to time, and Q ′ (·) represents the derivative of Q(·).

Based on Equation (2), a network implementation of sparse coding, shown in **Figure 2**, was proposed by Olshausen and Field (1997) who suggested that a feedforward-feedback loop can implement sparse coding. The input to the model was natural images that had been whitened using a filter that resembles the center-surround structure of retinal ganglion RFs. However, the original sparse coding model was not biologically plausible in several aspects, such as the possibility of negative spiking rates and the violation of Dale's law. In addition, the input the the model was not split into separate ON and OFF channels. Finally, this network imposed feedback synaptic connections that were anti-symmetric to the corresponding feedforward connections (i.e., equal but opposite in sign) and it was unclear how such symmetry could be achieved using biologically plausible mechanisms.

### 2.2. Structure of Our Model

We propose a two-layer network with rate-based neurons that models the activities of LGN cells (first layer), and simple cells (second layer), respectively (**Figure 3**). The model is based on a locally competitive algorithm that efficiently implements sparse coding with neural dynamics with non-negative spiking rates (Rozell et al., 2008).

We first define the parameters of the model that will be used throughout the paper. A summary of all symbols defined below is shown in **Table 1**. There are 2N LGN cells in the first layer, with N ON LGN cells and N OFF LGN cells, and M simple cells in the second layer. Denote **x** = [x1, · · · , x2N] T as the vector of input stimuli to the first layer. Denote **x**ON as the input to ON LGN cells (the first N elements of **x**) and **x**OFF as the input to OFF LGN cells (the last N elements of **x**), i.e., **x** = [**x** T ON, **x** T OFF] T .

Denote **v** L and **s** L as 2N × 1 vectors that represent membrane potentials and firing rates of LGN cells in the first layer. Denote **v** L ON, **s** L ON, **v** L OFF, and **s** L OFF as N × 1 vectors that represent the membrane potentials and firing rates of ON and OFF LGN cells, i.e., **v** L = [**v** L ON T , **v** L OFF T ] T and **s** L = [**s** L ON T ,**s** L OFF T ] T . Similarly, **v** C and **s** C are M ×1 vectors that represent membrane potentials and firing rates of M cortical simple cells in the second layer.

In our model, there are several important connections: feedforward (up) excitatory and inhibitory connections from LGN cells to simple cells, feedback (down) excitatory and inhibitory connections from simple cells to LGN cells, and self-excitatory connections of simple cells that represent selfexcitation. Definitions of connections are described below. One aspect of the model that lacks biological plausibility is existence of inhibitory connections between thalamus and cortex, but we propose biologically plausible neural circuits of implementing this aspect of the model in the Discussion section.

Denote **A** u,+ ON as an N × M matrix with non-negative elements that represents the feedforward excitatory connections from ON LGN cells to simple cells. Each column of **A** u,+ ON represents connections from N ON LGN cells to a simple cell. Similarly, denote **A** u,+ OFF as an N × M matrix with non-negative elements that represents the feedforward excitatory connections from OFF LGN cells to simple cells. Denote **A** u,− ON and **A** u,− OFF as N × M matrices with non-positive elements that represent inhibitory connections from ON and OFF LGN cells to simple cells, respectively. Denote **A** u,+ and **A** u,− as 2N × M matrices that

inhibitory connections, respectively. Upward and downward arrows are for feedforward and feedback pathways. Notation defined in the main text.

### TABLE 1 | Model symbols.

represents all excitatory and inhibitory connections from LGN to V1; then we have **A** u,+ = [**A** u,+ ON **A** u,+ OFF] and **A** u,− = [**A** u,− ON **A** u,− OFF].

For the feedback pathway, similar notation is used except superscript "d" represents feedback connections from simple cells to LGN cells. Therefore, we have **A** d,+ = [**A** d,+ ON **A** d,+ OFF] and **A** d,− = [**A** d,− ON **A** d,− OFF].

Using the notation defined above, the dynamics of ON and OFF LGN cells located in the first layer are given by

$$\begin{aligned} \tau\_{\rm L} \dot{\mathbf{v}}\_{\rm ON}^{\rm L} &= -\mathbf{v}\_{\rm ON}^{\rm L} + \mathbf{x}\_{\rm ON} + \mathbf{A}\_{\rm ON}^{\rm d,+} \mathbf{s}^{\rm C} + \mathbf{A}\_{\rm ON}^{\rm d,-} \mathbf{s}^{\rm C} + s\_{\rm b} \\ \mathbf{s}\_{\rm ON}^{\rm L} &= \max(\mathbf{v}\_{\rm ON}^{\rm L}, 0) \end{aligned} \tag{3}$$

and

$$\begin{aligned} \mathbf{r}\_{\rm L} \dot{\mathbf{v}}\_{\rm OFF}^{\rm L} &= -\mathbf{v}\_{\rm OFF}^{\rm L} + \mathbf{x}\_{\rm OFF} + \mathbf{A}\_{\rm OFF}^{\rm d,+} \mathbf{s}^{\rm C} + \mathbf{A}\_{\rm OFF}^{\rm d,-} \mathbf{s}^{\rm C} + \mathbf{s}\_{\rm b}, \\ \mathbf{s}\_{\rm OFF}^{\rm L} &= \max(\mathbf{v}\_{\rm OFF}^{\rm L}, \mathbf{0}), \end{aligned} \tag{4}$$

where τL is the time constant of the membrane potentials of LGN cells, sb is a constant that represents the instantaneous firing rate of the background input (i.e., from neurons outside the network), and the max operation represents the firing dynamics such that a cell only fires when the membrane potential is above a threshold.

Therefore, using the combined notation for ON and OFF LGN cells, the dynamics of LGN cells can be written as

$$\begin{aligned} \tau\_{\mathbf{L}} \dot{\mathbf{v}}^{\mathcal{L}} &= -\mathbf{v}^{\mathcal{L}} + \mathbf{x} + (\mathbf{A}^{\mathcal{d},+} + \mathbf{A}^{\mathcal{d},-}) \mathbf{s}^{\mathcal{C}} + s\_{\mathbf{b}} \\ \mathbf{s}^{\mathcal{L}} &= \max(\mathbf{v}^{\mathcal{L}}, \mathbf{0}). \end{aligned} \tag{5}$$

The dynamics of simple cells located in the second layer is given by

$$\begin{split} \tau\_{\text{C}} \dot{\mathbf{v}}^{\text{C}} &= -\left(\mathbf{v}^{\text{C}} - \mathbf{v}\_{\text{leak}}^{\text{C}}\right) + \mathbf{A}\_{\text{ON}}^{\text{u}, + \text{T}} \mathbf{s}\_{\text{ON}}^{\text{L}} + \mathbf{A}\_{\text{ON}}^{\text{u}, - \text{T}} \mathbf{s}\_{\text{ON}}^{\text{L}} \\ &+ \mathbf{A}\_{\text{OFF}}^{\text{u}, + \text{T}} \mathbf{s}\_{\text{OFF}}^{\text{L}} + \mathbf{A}\_{\text{OFF}}^{\text{u}, - \text{T}} \mathbf{s}\_{\text{OFF}}^{\text{L}} + \mathbf{s}^{\text{C}}, \end{split} \tag{6}$$

which can be reformulated as

$$\begin{aligned} \tau\_{\mathbf{C}} \dot{\mathbf{v}}^{\mathbf{C}} &= -\mathbf{v}^{\mathbf{C}} + \mathbf{v}\_{\text{leak}}^{\mathbf{C}} + (\mathbf{A}^{\mathbf{u},+} + \mathbf{A}^{\mathbf{u},-})^{T} \mathbf{s}^{\mathbf{L}} + \mathbf{s}^{\mathbf{C}} \\ \mathbf{s}^{\mathbf{C}} &= \max(\mathbf{v}^{\mathbf{C}} - \lambda, \mathbf{0}), \end{aligned} \tag{7}$$

where τC is the time constant of the membranes of simple cells and λ is the threshold of the rectifying function of firing rates. In addition, λ is a positive constant that introduces sparseness into the model, **s** C represents the self-excitation of simple cells, which comes from reformulating the model equations of the locally competitive algorithm (Rozell et al., 2008), and **v** C leak represents the change of membrane potential caused by leakage currents. The leakage currents drive the membrane potentials of simple cells to their resting potentials when there is no external input, i.e., **v** C is zero. Therefore, the steady states of the model dynamics are **v** L = sb, **s** L = sb, **v** C = 0, and **s** C = 0, which implies that **v** C leak = −(**A** u,+ + **A** u,−) T **s**b, where **s**b is a vector whose elements are all equal to sb. Equations 5 and 7 are solved simultaneously by iteration to obtain values of membrane potentials and firing rates.

The codes to run the model are available from ModelDB (http://modeldb.yale.edu/247970).

### 2.3. Learning Rule

The learning process of the model is based on a Hebbian or anti-Hebbian rule, namely that the change of synaptic strength is related only to local pre-synaptic and post-synaptic activities.

The learning rules are given by

$$\begin{aligned} \Delta \mathbf{A}^{\mathbf{u},+} &= \eta \langle (\mathbf{s}^{\mathcal{L}} - s\_{\mathbf{b}}) \mathbf{s}^{\mathcal{C}^{T}} \rangle \\ \Delta \mathbf{A}^{\mathbf{u},-} &= \eta \langle (\mathbf{s}^{\mathcal{L}} - s\_{\mathbf{b}}) \mathbf{s}^{\mathcal{C}^{T}} \rangle \\ \Delta \mathbf{A}^{\mathbf{d},+} &= -\eta \langle (\mathbf{s}^{\mathcal{L}} - s\_{\mathbf{b}}) \mathbf{s}^{\mathcal{C}^{T}} \rangle \\ \Delta \mathbf{A}^{\mathbf{d},-} &= -\eta \langle (\mathbf{s}^{\mathcal{L}} - s\_{\mathbf{b}}) \mathbf{s}^{\mathcal{C}^{T}} \rangle, \end{aligned} \tag{8}$$

where η is the learning rate, h·i is the ensemble average operation over some samples, **s** L −sb is the vector such that each element of vector **s** L is subtracted by scalar sb, and (**s** L − sb)**s** CT is the matrix given by the outer product of vectors **s** L − sb and **s** C.

The change of synaptic strength depends only on the presynaptic activity (**s** L ) and post-synaptic activity (**s** C). Therefore, this learning rule is local and thus biophysically realistic. In obedience to Dale's law, all the weights of **A** u,+ and **A** d,+ are kept non-negative and all weights of **A** u,− and **A** d,− are kept non-positive during learning. If any synaptic weight changes sign after applying Equation (8), the synaptic weight is set to zero. In addition, after each learning iteration, synaptic weights are multiplicatively normalized to ensure that Hebbian learning is stable. Specifically, each column of **A** u,+ and **A** d,− is normalized to norm l1 and each column of **A** u,− and **A** d,+ is normalized to norm l2. The multiplicative normalization of synaptic weights may be achieved by homeostatic mechanisms (Turrigiano, 2011), but these are not implemented here as they are not the focus of this paper.

### 2.4. Input

The data set used in our simulation consists of 10 pre-whitened 512 × 512 pixel images of natural scenes provided by Olshausen and Field Olshausen and Field (1996). Some previous studies of sparse coding (efficient coding) also used this data set (Olshausen and Field, 1996; Wiltschut and Hamker, 2009; Zylberberg et al., 2011; Zhu and Rozell, 2013). The input stimuli to the model are chosen to be 16 × 16 pixel image patches sampled from these 10 pre-whitened 512 × 512 pixel images, similar to previous studies (Zylberberg et al., 2011; Zhu and Rozell, 2013).

The pre-whitening process mimics the spatial filtering of retinal processing up to a cut-off frequency determined by the limits of visual acuity (Atick and Redlich, 1992). This process is realized by passing the original natural images through a zerophase whitening filter with root-mean-square power spectrum,

$$R(f) = f e^{-(f/f\_\epsilon)^4},\tag{9}$$

where fc = 200 cycles/picture (Olshausen and Field, 1997). **Figure 4** shows the spatial and frequency profiles of the prewhitening filter. The spatial profile of the filter (**Figure 4C**), obtained by taking the 2D inverse Fourier transform of the filter in the 2D frequency domain, approximates center-surround RFs of LGN cells in a pixel image. The pre-whitening filter described in Equation (9) is widely used in computational studies (Jehee et al., 2006; Jehee and Ballard, 2009; Wiltschut and Hamker, 2009; Zhu and Rozell, 2013).

The pre-whitened images are then scaled to variance 0.2 similar to Olshausen and Field (1997).  If a pixel intensity in a pre-whitened image patch is negative, we assign the absolute value of the pixel intensity to the input of the OFF LGN cell and set the input of the corresponding ON LGN cell to zero; if the pixel intensity is positive, we set the input of the ON LGN cell to the pixel intensity and set the input to the OFF LGN cell to zero.

### 2.5. Training

Since we use 16 × 16 pixel images as the input to our model, 256 ON and 256 OFF LGN cells (N = 256) are required in the first layer. We use 256 simple cells (M = 256) in the second layer. The first-order Euler method is implemented to solve the dynamical system described by Equation 5 and 7. We choose

a time scale in which the passive membrane time constant is τL = τC = 12 ms, within the physiological range (Dayan et al., 2001), and sparsity level λ = 0.6 similarto Zhu and Rozell (2013). The spontaneous firing rate, sb, is chosen as sb = 2 Hz, the median of spontaneous firing rates of the mouse LGN cells in the experimental study of Tang et al. (2016). There are 30 integration time steps, with an integration time step of 3ms, for calculating neuronal responses per stimulus with the assumption that neural responses will converge after 30 iterations.

Learning rules in Equation (8) are used to update the synaptic weights. For the normalization step after each learning iteration, each column of **A** u,+ and **A** d,− is normalized to have norm l1 and each column of **A** u,− and **A** d,+ is normalized to have norm l2. Elements of **A** u,+ and **A** d,+ are non-negative and initialized randomly using an exponential distribution with mean 0.5. **A** u,− and **A** d,− are initialized randomly with non-positive elements that are sampled from an exponential distribution with mean −0.5. Then, synaptic weights are normalized before the learning process starts. Results shown in this paper are from simulations with l1 = l2 = 1 (unit norm), as used in the previous study by Rozell et al. (2008). The learning rule based on the average activities of a mini-batch is applied; i.e., in every epoch, a minibatch that consists of 100 randomly selected 16 × 16 pixel images sampled from the data set is used. Before the training process of natural image patches, the model is pre-trained on white noise for 10, 000 epochs to mimic the process of pre-development of the visual system; the learning rate is 0.5 in pre-training. To ensure that the weights converge after learning on natural image patches, we use 30, 000 epochs in the training process, where the learning rate is 0.5 for the first 10, 000 epochs, 0.2 for the second 10, 000 epochs and 0.1 for the third 10, 000 epochs. Learning rates were chosen to ensure stable convergence of the weights in a reasonable time; but the results are not sensitive to moderate changes.

### 2.6. Recovering Receptive Fields of Model Simple Cells Using White Noise

In order to estimate the RFs of model simple cells in a systematic way, we use the method of spike-triggered averaging to find the pattern that each simple cell responds to on average (Schwartz et al., 2006). Using K 16 × 16 white noise stimuli **n**1, · · · , **n**K, we present pre-processed stimuli to the model, record the firing rates of a simple cell, s1, · · · ,sK, and then estimate the RF, **F**, of the simple cell as the weighted average,

$$\mathbf{F} = \frac{s\_1 \mathbf{n}\_1 + \dots + s\_K \mathbf{n}\_K}{s\_1 + \dots + s\_K}. \tag{10}$$

We used 70, 000 white noise stimuli, i.e., K = 70, 000.

In our simulations, we have two versions of estimated RFs using the two different methods of pre-processing the white noise stimuli: the same pre-whitening filter for natural scenes (Equation 9) and a low-pass filter described by

$$L(f) = e^{-\left(f/f\_0\right)^4}.\tag{11}$$

# 2.7. Fitting Receptive Fields to Gabor Functions

The RFs of visual cortical cells are often modeled using a 2D Gabor function G(x, y) of the form

$$\begin{split} &G(\mathbf{x}, \mathbf{y}; \mathbf{x}\_0, \mathbf{y}\_0, \sigma\_{\mathbf{x}}, \sigma\_{\mathbf{y}}, \hat{f}\_s, \boldsymbol{\beta}, \boldsymbol{\theta}, \boldsymbol{\phi}, \boldsymbol{\phi}) \\ &= \beta \cos(2\pi f\_s \mathbf{x}' + \boldsymbol{\phi}) e^{-(\frac{\mathbf{x}'}{\sqrt{2}\sigma\_{\mathbf{x}}})^2 - (\frac{\mathbf{y}'}{\sqrt{2}\sigma\_{\mathbf{y}}})^2} \end{split} \tag{12}$$

with

$$\begin{aligned} \mathbf{x}' &= (\mathbf{x} - \mathbf{x}\_0)\cos\theta + (\mathbf{y} - \mathbf{y}\_0)\sin\theta\\ \mathbf{y}' &= -(\mathbf{x} - \mathbf{x}\_0)\sin\theta + (\mathbf{y} - \mathbf{y}\_0)\cos\theta,\end{aligned} \tag{13}$$

where β is the amplitude, (x0, y0) is the center, σx and σy are standard deviations of the Gaussian envelope, θ is the orientation, fs is the spatial frequency, and φ is the phase of the sinusoid wave (Ringach, 2002). These parameters are fitted using the built-in MATLAB (version R2016b, MathWorks, MA, USA) function, lsqcurvefit, that efficiently solves non-linear curvefitting problems using a least-squares method. The fitting error is defined as the square of the ratio between the fitting residual and RF.

To ensure that results were only reported for RFs that were well-fitted to Gabor functions, we excluded RFs for which either (1) the synaptic fields had fitting error larger than 40% or (2) the center of the fitted Gabor functions lay either outside the block, or within one standard deviation of the Gaussian envelope of the block edge (Zylberberg et al., 2011). After applying these two quality control measures, 140 out of 256 model cells remained for subsequent analysis.

# 2.8. Measuring the Overlap Index Between ON and OFF Sub-regions

To investigate the extent of overlap between ON and OFF subregions, we used an overlap index that was used in experimental studies (Schiller et al., 1976; Martinez et al., 2005). Similar to the method used in Martinez et al. (2005), each ON and OFF excitatory sub-region was fitted by an elliptical Gaussian function:

$$h(\mathbf{x}, \mathbf{y}; \mathbf{x}\_0, \mathbf{y}\_0, a, b, \theta, \nu) = \frac{\mathcal{Y}}{2\pi \, ab} e^{-\frac{\mathbf{x}^2}{2a^2} - \frac{\mathbf{y}^2}{2b^2}} \tag{14}$$

where γ is the amplitude, a and b are half axes of the ellipse, and x ′ and y ′ are the transformed coordinates given by Equation (13). If there are more than one ON (or OFF) sub-regions for the simple cell, only the most significant sub-region was fitted by the elliptic Gaussian. If either the ON or OFF sub-region of a simple cell has fitting error larger than 40% or has the half axis, a, larger than 3 pixels, this simple cell is excluded. 92 simple cells remained for the analysis of overlap index.

The overlap index, Io, is then defined as

$$I\_{\rm o} = \frac{W\_{\rm ON} + W\_{\rm OFF} - d}{W\_{\rm ON} + W\_{\rm OFF} + d}, \ (-1 < I\_{\rm o} \le 1) \tag{15}$$

where WON and WOFF are the half width measured at the point where the response is 30% of the maximal response, and d is the distance between the centers of ON and OFF sub-regions. Smaller values of Io indicate more segregation between ON and OFF sub-regions.

### 2.9. Measuring the Push-Pull Index

The push-pull effect of the model was measured by a pushpull index (Martinez et al., 2005). First, for each simple cell, we recorded the membrane potential, P, when the preferred input (the synaptic field) was presented to the model. Next, we recorded the membrane potential, N, while presenting the opposite of preferred input to the model. To make the measurement independent of the relative strength of different simple cells, P and N were normalized by

$$P = \frac{P}{\max(|P|, |N|)} \text{ and } N = \frac{N}{\max(|P|, |N|)}. \tag{16}$$

The Push-pull index, Ip, is then defined as

$$I\_{\mathfrak{P}} = |P + N|, \text{ ( $0 \le I\_{\mathfrak{P}} \le 2$ )}.\tag{17}$$

Smaller values of Ip indicate stronger push-pull effect.

### 2.10. Measuring Contrast Invariance of Orientation Tuning

The method in (Zhu and Rozell, 2013) was used to investigate contrast invariance of orientation tuning and the procedure is as follows. First, an exhaustive search was performed to find the preferred circular sinusoidal grating in the parameter space of the following ranges: radius of the grating was between 1 pixel and 2.5 min(σx, σy) (smaller than 8 pixels which is the maximum radius for a 16 × 16 image patch) with the stepsize of 1 pixel ; spatial frequency was between 0.05 and 0.3 cycles/pixel with the stepsize of 0.05 cycles/pixel; orientation was between 0 and 180 degrees with the stepsize of 5 degrees; phase was between 0 and 360 degrees with the stepsize of 30 degrees. Next, we measured the mean response to the drifting grating with orientations between 0 and 180 degrees with the stepsize of 5 degrees while varying the contrast of the stimuli from 0.2 to 1 in increments of 0.2, where contrast is defined as the amplitude of the sinusoidal grating. The orientation tuning curve for each contrast level was then fit to the Gaussian function and the half-height bandwidth of the Gaussian fit was calculated. The slope of the linear fit to half-height bandwidth vs. contrast for the cell was used to plot the population statistics of contrast invariance (Alitto and Usrey, 2004). Here, only 68 model simple cells that have oriented RFs located well within the 16 × 16 image patch were selected for the analysis.

### 3. RESULTS

After learning, synaptic weights between LGN and V1 display spatial structures similar to those observed in recordings of neurons in V1, such as oriented Gabor-like filters and nonoriented blobs. Since both excitatory and inhibitory connections from ON and OFF LGN cells contribute to the responses of simple cells, we use the synaptic field (**S**f ) defined as

$$\mathbf{S}\_{\rm f} = (\mathbf{A}\_{\rm ON}^{\rm u,+} + \mathbf{A}\_{\rm ON}^{\rm u,-}) - (\mathbf{A}\_{\rm OFF}^{\rm u,+} + \mathbf{A}\_{\rm OFF}^{\rm u,-}) \tag{18}$$

to visualize the overall synaptic weights from ON and OFF LGN cells. The synaptic fields of 140 model simple cells that meet the two quality control measures (see the Materials and Methods section) are shown in **Figure 5**, where each block represents the overall effect of the feedforward connections from ON and OFF LGN cells to a simple cell. Note that although **Figure 5** displays spatial patterns that are similar to experimental RFs, strictly they represent the synaptic weights from LGN cells to simple cells, which ignores the early visual processing before LGN. However, the RFs of the model are systematically investigated in the following sections.

In the remaining results, we show that the synaptic weights exhibit several properties that have been observed experimentally, including segregation of ON and OFF subregions, push-pull effect, phase-reversed feedback, diverse shapes of simple cell RFs, and contrast invariance of orientation tuning.

# 3.1.  However, it is impossible for a model that combines ON and OFF LGN input into a single linear input to explain this important phenomenon. Our model separates ON and OFF LGN input and shows that the learned feedforward excitatory connections from ON and OFF LGN cells to simple cells can

account for the segregation of ON and OFF sub-regions of simple cells.

ON and OFF excitatory regions of some example simple cells are displayed in **Figure 6A**. In our model, there are 256 ON LGN and 256 OFF LGN cells located evenly on a 16 × 16 image, so each block in **Figure 6A** represents 256 excitatory weights from ON or OFF LGN cells to a simple cell. **Figure 6A** shows that these excitatory connections form spatial patterns such as bars and blobs. Furthermore, a careful examination of the patterns shows that excitatory connections from ON LGN cells are normally adjacent to patterns of excitatory connections from OFF LGN cells, but the two patterns do not overlap, as can be seen when contour plots for the ON and OFF excitatory regions are overlaid in **Figure 6B**.

We quantified the segregation of ON and OFF sub-regions using the overlap index (defined in the Materials and Methods section). The histogram of the overlap index for simple cells in an experimental study (Martinez et al., 2005) is re-plotted in **Figure 6C**. Consistent with the experimental data, 88 out of 92 model simple cells had an overlap index smaller than 0.1 (**Figure 6D**), which indicates that the ON and OFF sub-regions are well-separated in a large majority of the population. The synaptic fields of simple cells whose overlap indices are larger than 0.1 are shown in **Figure 6E**, revealing that most of them have low spatial frequencies.

### 3.2. Push-Pull Effect

Simple cells are also found to have push-pull responses; i.e., if one contrast polarity excites a cell, the opposite contrast polarity tends to inhibit it (Jones and Palmer, 1987b; Ferster, 1988; Hirsch et al., 1998; Martinez et al., 2005). Even though this effect has been

connections from ON or OFF LGN cells to a simple cell. The color magenta represents excitatory connections. (B) Red and blue contours represent excitatory connections from ON and OFF LGN cells, respectively. Connections that are smaller than 20% of the maximal connection were removed to only show the substantial weights. The number in each block indicates the overlap index. (C) Histogram of the overlap index for simple cells in cat V1. Re-plotted from Figure 3C in Martinez et al. (2005). (D) Histogram of the overlap index for model simple cells. (E) Synaptic fields of the four simple cells with overlap index larger than 0.1.

observed in many experimental studies, to our knowledge there has not been a learning model proposed that can explain how this effect emerges. Again, a model that separates ON and OFF LGN input is necessary to investigate the emergence of the push-pull effect. In this section, we show that the push-pull effect for simple cells naturally emerges as a result of neural learning.

Some examples of ON excitatory and OFF inhibitory synaptic weights (**A** u,+ ON and **A** u,− OFF, respectively) are shown in **Figure 7A**. The patterns of **A** u,+ ON are similar to the ones of **A** u,− OFF and they are located at similar locations, as can be seen from the highly overlapped contours in **Figure 7B**. However, the degree of overlap is different between the examples.

Analogous results to the above also hold for learned excitatory connections from OFF LGN cells, **A** u,+ OFF, and inhibitory connections from ON LGN cells, **A** u,− ON (data not shown).

We then quantified the push-pull effect using push-pull index (defined in the Materials and Methods section). Both the histograms of push-pull index for experimental data (**Figure 7C**) and model simple cells (**Figure 7D**) peaked near zero and showed an decreasing trend. Model simple cells showed even stronger push-pull index with more simple cells having push-pull index close to zero. The synaptic fields of simple cells with push-pull indices larger than 0.2 are shown in **Figure 7E**, showing that most of them have low spatial frequencies.

# 3.3. Phase-Reversed Feedback

The experimental study of Wang and colleagues suggests that the synaptic feedback from V1 to LGN is phase-reversed with respect to the feedforward connections (Wang et al., 2006). For example, the connection from a simple cell to an ON-center LGN cell will be excitatory if the ON-center is aligned in visual space to the OFF sub-field of simple cell (i.e., phase-reversed). Conversely, if the ON-center is aligned to the ON sub-field of the simple cell, the connection will be inhibitory. Our learning model with separate ON and OFF LGN cells enables us to investigate the feedback effect from simple cells to LGN cells. In this section, we show that phase-reversed feedback arises in the structures of learned connections.

connections from OFF LGN cells (A u,− OFF), respectively. Connections that are smaller than 20% of the maximal connection were removed to only show substantial weights. The number in each block indicates the push-pull index. (C) Histogram of the push-pull index for simple cells in cat V1. Re-plotted from Figure 4B in Martinez et al. (2005). (D) Histogram of the push-pull index for model simple cells. (E) Synaptic fields of the six simple cells with push-pull index larger than 0.2.

Feedback from simple cells to LGN cells occurs via both excitatory connections, **A** d,+ x , and inhibitory connections, **A** d,− x , with the overall effect characterized by **A** d x = **A** d,+ x + **A** d,− x , where x = ON or OFF depending on the type of LGN cell. Therefore, the overall feedback to ON LGN cells, denoted as **A** d ON, can be represented by **A** d ON = **A** d,+ ON + **A** d,− ON.

The ON and OFF sub-fields of simple cells receptive fields are characterized by the positive and negative regions of the synaptic field defined in Equation (18). The scatter plots in **Figure 8** show that relationship expected from phase-reversed feedback. **S**f is highly positively correlated with **A** d OFF (correlation coefficient r = 0.90), while **S**f is highly anti-correlated with **A** d ON (correlation coefficient r = −0.92). According to the figure, wherever **S**f is positive, indicating the ON sub-field, the feedback to ON LGN cells, **A** d ON, is very likely to be negative and the feedback to OFF LGN cells, **A** d OFF, tends to be positive; however, wherever **S**f is negative, indicating the OFF-field, the converse is true: the feedback to ON LGN cells, **A** d ON, is very likely to be positive and the feedback to OFF LGN cells, **A** d OFF, tends to be negative. This corresponds to a phase-reversed feedback from V1 to LGN.

This phase-reversed feedback from V1 to LGN can be explained by the learning dynamics of LGN and simple cells described in Equation 8. The learning rule shows that **A** u,+ and **A** d,− are updated with the same magnitude of synaptic change but opposite in sign (and are normalized with the same norm l1). Similarly, **A** u,− and **A** d,+ are updated with the same magnitude of synaptic change but opposite in sign (and are normalized with the same norm l2). These anti-symmetries are a consequence of having Hebbian learning for the forward weights and anti-Hebbian learning for the feedback weights. In both cases the magnitude of weight change is proportion to the production of pre- and post-synaptic spike rates, but the sign of the change is opposite. The anti-symmetry arises because roles of pre- and post-synaptic rates are interchanged in forward vs. feedback directions, in combination with the sign change. Simulation results show that **A** u,+ converges to −**A** d,−

FIGURE 8 | Synaptic fields, Sf (defined in Equation 18), vs. feedback to ON and OFF LGN cells, A d ON and A d OFF. Sf is highly positively correlated with A d OFF (correlation coefficient *r* = 0.90) and Sf is highly anti-correlated with A d ON (correlation coefficient *r* = −0.92).  On the contrary, A d OFF tends to be smaller than zero and A d ON tends to be greater than zero if Sf is negative.

and **A** u,− converges to −**A** d,+ even during pre-development when white noise is used as the input to the model, as illustrated in **Figure 9**.

# 3.4. The Diversity of Model Receptive Fields Resembles That Observed Experimentally for Simple Cells

In this section, we show that the range of spatial structures of RFs of our model have a close resemblance to experimental data.

RFs were calculated from the model by simulating experiments in which Gaussian white noise is presented as a visual stimulus, and the spike triggered average is used to estimate RFs. As the presentation of white noise may cause adaptive effects in the early stages visual system relative to natural images, we considered two versions of the model, one with the standard pre-whitening filter (Equation 9) modeling center-surround processing, and a second without pre-whitening in which the filter is replaced by a low-pass filter (Equation 11) with the same upper cut-off frequency as pre-whitening filter. We use pre-whitened RFs and low-pass RFs to represent of simple cell RFs estimated using the pre-whitening filter and low-pass filter.

Some examples of pre-whitened RFs, low-pass RFs and synaptic fields are shown in **Figure 10**, which shows that prewhitened RFs and low-pass RFs are similar to synaptic fields. However, pre-whitened RFs tend to have more and thinner stripes, which indicates a narrower tuning to a somewhat higher spatial frequency. For a simple cell tuned to very low spatial frequencies (bottom right blocks), the RF recovered with prewhitening was a poor match to the original synaptic field, but for RF recovered with low-pass filtering it was fair.

Early studies show that RFs of simple cells can be welldescribed by 2D Gabor functions described in Equation (12) (Jones and Palmer, 1987a; Ringach, 2002). For our model, most RFs could be well-fitted by Gabor functions with suitable choices of parameters with small fitting errors, as shown in **Figure 11A**. Note that although the fitting error of blob-like RFs might be low, the parameter choices are not necessarily reasonable, in that they are poorly constrained and the process of Gabor fitting imposes an a priori hypothesis that the RF is a 2D-Gabor function even though it is clearly not Gabor-like. The pre-whitened RFs with fitting errors larger than 40% (**Figure 11B**) are cells whose synaptic fields have low spatial frequencies (**Figure 11C**), because pre-whitened RFs of these cells matched poorly to the original synaptic fields (**Figure 10B**). Low-pass RFs of all 140 selected model cells have fitting errors smaller than 40% with 132 of them having fitting errors smaller than 20% (data not shown).

Using fitted parameters of Gabor functions, Ringach constructed a scatter plot of nx = σxfs vs. ny = σyfs to analyze the spatial structures of RFs in V1 over the population (Ringach, 2002). Such plots have subsequently been used by many researchers to investigate the distributions of model simple cell RFs (Rehn and Sommer, 2007; Wiltschut and Hamker, 2009; Zylberberg et al., 2011). nx and ny are the width and length of the Gabor function measured in the number of cycles of the spatial frequency (i.e., across and along the stripes). Ringach noted that blob-like RFs are mapped to points near the origin, while RFs with elongated sub-regions are mapped to points away from the origin (Ringach, 2002). In addition, nx and ny are directly related with the half-magnitude spatial frequency bandwidth 1f and

example model cells. The low-pass filter described in Equation (11) was used to filter white noise stimuli.

orientation bandwidth 1θ of the fitted Gabor function,

$$\begin{aligned} \Delta f &:= h(n\_{\mathbf{x}}) = \log\_2 \left( \frac{1 + \frac{\sqrt{2\ln 2}}{2\pi n\_{\mathbf{x}}}}{1 - \frac{\sqrt{2\ln 2}}{2\pi n\_{\mathbf{x}}}} \right) \text{ in cotaves} \\ \Delta \theta &:= g(n\_{\mathcal{V}}) = 2 \arctan \left( \frac{\sqrt{2\ln 2}}{2\pi n\_{\mathcal{V}}} \right) \text{ in degrees.} \end{aligned} \tag{19}$$

Both h(nx) and g(ny) are monotonically decreasing functions; i.e., the larger nx and ny, the smaller 1f and 1θ. Note that h(nx) is not well-defined when nx < √ 2 ln 2/2π (≈ 0.13), i.e., when the lower half-magnitude frequency do not exist. This corresponds to the region in which Gabor fitting gives ambiguous fits for parameters like spatial frequency and orientation, because oriented RFs with low spatial frequency might lie in this region as well.

We plot nx vs. ny and 1f vs. 1θ for RFs obtained from both the model and experimental studies in **Figure 12**. However, the different pre-processing filters for white noise stimuli have a dramatic influence on the distributions of nx vs. ny, shifting the distribution for low-pass RFs to the left of pre-whitened RFs, in closer agreement to the experimental data. As mentioned earlier, pre-whitened RFs tend to have more stripes relative to the lowpass RFs, so they are mapped to points away from the origin compared to low-pass RFs. In addition, the distribution of lowpass RFs is continuous from the origin, while there is a gap between points near the origin and points away from the origin for pre-whitened RFs. The inset sub-plots of **Figure 12** show that data points near the origin might be orientated RFs with low spatial frequencies and blob-like RFs might not be necessarily mapped to points near the origin.

In general, oriented RFs are well-described by Gabor functions and low-pass RFs better resemble the distribution of experimental data compared with pre-whitened RFs.

## 3.5. Contrast Invariance of Orientation Tuning

Another important property of simple cells is contrast invariance of orientation tuning; i.e., the width of the orientation tuning curve is maintained when the contrast of the stimulus changes, as demonstrated in **Figure 13A**. The orientation tuning curves with various stimulus contrasts for a model simple cell are shown in **Figure 13B**, where the bandwidths of each curve remain the same while the responses become larger when the stimulus contrast increases. For a study of contrast invariance of V1 population in ferret, the histogram of the slope of the linear fit of half-width bandwidth vs. contrast (**Figure 13C**) showed that most cells were contrast invariant with the slope close to zero (Alitto and Usrey,

2004). **Figure 13D** shows that most model cells have the slope around zero, which is consistent with experimental data.

# 4. DISCUSSION

## 4.1. Relationship With Sparse Coding

Sparse coding has been successful in modeling simple cell receptive fields (RFs) and has been used by many researchers over the past years. Our model is based on an algorithm that efficiently implements sparse coding (Rozell et al., 2008), and is therefore closely related to the original concept of sparse coding (Olshausen and Field, 1996).

If we define **A** as a 2N × M matrix that represents the overall effect caused by excitatory and inhibitory connections from 2N LGN cells to M simple cells, we have **A** = **A** u,+ + **A** u,−. The dynamics of simple cells described in Equation (7) can be rewritten as

$$\mathbf{r}\_{\rm C}\dot{\mathbf{v}}^{\rm C} = -\mathbf{v}^{\rm C} + \mathbf{A}^{T}(\mathbf{s}^{\rm L} - s\_{\rm b}) + \mathbf{s}^{\rm C}.\tag{20}$$

As illustrated in **Figure 9**, **A** u,+ → −**A** d,− and **A** u,− → −**A** d,+ during learning. Therefore, we have **A** d,− + **A** d,+ = −**A** u,+ − **A** u,− = −**A**. The dynamics of LGN cells described in Equation (5) can be rewritten as

$$
\tau\_{\mathbf{L}} \dot{\mathbf{v}}^{\mathcal{L}} = -\mathbf{v}^{\mathcal{L}} + \mathbf{x} - \mathbf{A} \mathbf{s}^{\mathcal{C}} + s\_{\mathbf{b}}.\tag{21}
$$

If the columns of **A** are seen as the basis vectors of a generative model, **As**C can be seen as the linear reconstruction of the input using learned basis vectors and thus **x** − **As**C represents the residual error, which is similar to **r** of the sparse coding formulation given in Equation (2). Therefore, the residual error used to update the basis vectors of the original sparse coding model is represented by the responses of LGN cells in our model.

To incorporate Dale's law, non-negative connections, **A** u,+, and non-positive connections, **A** u,−, are employed in our model to represent the positive and negative elements of **A**. **A** u,+ and **A** u,− are not co-active in general, which suggests that **A** u,+ ≈ [**A**]+ and **A** u,− ≈ [**A**]−, where [ · ]+ preserves the positive elements and sets negative elements to zero and [ · ]− preserves the negative elements and sets positive elements to zero.

In other words, our model is essentially a variant of sparse coding that employs separate connections to learn the positive and negative part of the overall connections.

## 4.2. Relationship With Predictive Coding

Our model is a hierarchical model with feedforward and feedback connections based on a locally competitive algorithm (Rozell et al., 2008). The structure of our model is essentially very similar to that of predictive coding models. To be more specific, the feedback from the second-layer neurons reconstruct the input. The residual error is computed at the first layer and then propagated to the second layer via feedforward connections.

Although our model presented here and the predictive coding model of Jehee and Ballard (2009) can explain phase-reversed feedback, the models differ in several respects. First, sparse

coding in our model is simply realized by the threshold of the rectifying function of firing rates for simple cells and this simple mechanism leads to simple neural circuits. Second, compared to the mechanism for determining simple cell responses one by one in their model, our model computes the responses in parallel. Third, our model generates diverse types of RFs that correspond well to experimental data. Finally, the phase-reversed effect is simply accounted for by the special pattern of learned connections, which also explains the segregation of ON/OFF sub-regions and push-pull effect for simple cells.

## 4.3. The Function of Spontaneous Activity

In the model proposed here, the dynamics of LGN cells described in Equation (5) has the background firing rate, sb, as part of the input to LGN cells. This spontaneous firing rate introduces a shift of the operating point for LGN cells. Given the responses of simple cells, **s** C, **x** − **As**C in Equation (21) represents the reconstruction residual error between the input and reconstruction. The residual error gives the difference between the real input and the representation produced by the model and it can be either positive or negative. To code for the signed quantities (residual error), Ballard and Jehee carried out a case-by-case study, leading to very complicated neural circuits (Ballard and Jehee, 2012). However, our model has a straightforward method for the implementation of solving signed quantities. The background firing rate, sb, in Equation (5) increases the residual errors by sb. Therefore, the membrane potential of LGN cell, **v** L , represents the residual error shifted up by sb. The threshold function in Equation (5) gives the firing rate of the LGN cell and it preserves the residual error in the interval of [−sb, ∞], which preserves the information of whether the model under-estimates or over-estimates the input stimuli and forces the connections to evolve through learning in the correct direction. In Equation (7), which describes simple cell dynamics, the effect of the spontaneous firing rate, sb, is removed by **v** C leak, a homeostatic mechanism employed by simple cells to maintain resting membrane potentials when there is no external input. The local learning rule described by Equation (8) also eliminates the effect of the spontaneous firing rate by subtracting it. The use of spontaneous firing rate makes the model much simpler and offers a new approach for solving the problem of signed quantities (residual errors). Experimental evidence shows that thalamocortical neurons can fire with bursts of action potentials without any synaptic input (Kandel et al., 2013), which suggests that the spontaneous firing activities might be used to encode the difference between input and feedback information.

# 4.4. Pre-processing of the Early Visual System

Atick and Redlich suggest that the retinal goal is to whiten the visual input up to a transition frequency such that input noise can also be suppressed (Atick and Redlich, 1992). The pre-whitening filter (Equation 9) approximately whitens the natural scenes up to the cut-off frequency.

However, for pre-processing white noise stimuli, two hypotheses are considered here. First, the filtering process of the early visual system can be described by the pre-whitening filter (Equation 9) whether or not the visual stimuli are natural scenes. Second, the early visual system is adaptive such that the visual stimuli are whitened up to a cut-off frequency. In this case, a low-pass filter (Equation 9) should be used, because white noise stimuli are already whitened across all frequencies. Our results suggest that estimated RFs using low-passed white noise match the experimental data much better than estimated RFs using prewhitened white noise. Further investigation of how visual stimuli are processed before they are fed to the visual cortex is needed to better understand the properties of simple cells.

# 4.5. The Role of l1 and l2

Each column of **A** u,+ and **A** d,− is normalized to norm l1 and each column of **A** u,− and **A** d,+ is normalized to norm l2. In other words, l1 represents the overall strength of feedforward excitatory connections and feedback inhibitory connections while l2 represents the overall strength of feedforward inhibitory connections and feedback excitatory connections. The results shown in this paper are based on l1 = 1 and l2 = 1; i.e., the strength of feedforward excitatory connections is equivalent to feedforward inhibitory connections, which leads to a strong push-pull effect in **Figure 7D**. If l2 is smaller than l1, the pushpull effect will be weaker and the distribution of the push-pull index will shift to the right. In addition, reducing l2 results in more blob-like receptive fields (data not shown).

# 4.6. Neural Circuits

Biologically realistic neural models can provide deeper insights into how real neural circuits function. The model proposed here contains a number of features that correspond to those in its biological counterpart, namely in terms of ON and OFF channels for LGN cells, positive neuronal responses, local computation, local learning rule, existence of feedback, and obedience to Dale's law.

In addition, our model incorporates inhibitory effects between LGN cells and cortical simple cells. As pointed out in the Materials and Methods section, for simplicity, inhibitory effects are implemented by direct inhibitory connections between two layers. However, in reality, long-range inhibitory effects should be implemented via interneurons that have inhibitory synapses. In this section, we will discuss several neural circuits of implementing inhibitory connections of our model.

Possible neural circuits that may be used to implement long-range inhibition are displayed in **Figure 14**. Assume that the overall inhibitory effects from LGN cells (with activity **s** L ) to cortical simple cells (with activity **s** C) can be represented by inhibitory connections, **A** −, between populations. We also assume that the learning rule of **A** − is local, i.e., that only depends on the responses of two populations (**s** L and **s** C). Long-range inhibition in our model is implemented via direct inhibitory connections, which is not biologically realistic (**Figure 14A**).

The circuit in **Figure 14B** implements inhibitory connections, **A** − (with non-positive weights), via a population of interneurons that have inhibitory connections, **A** −, with cortical simple cells. LGN cells are connected to interneurons via long-range identical excitatory connections, **I**; i.e., the interneurons copy the responses of LGN cells. For this structure, long-range excitatory connections, **I**, are fixed while **A** − are learned using the same learning rule in **Figure 14A**. In this case, the learning rule of **A** − is still local because the responses of interneurons are just **s** L and the model is still biologically plausible in terms of the local learning rule. Furthermore, the RFs of interneurons in the same layer as cortical simple cells should be LGN-like. Though V1 cortical cells with blob-like RFs were found in different species (Kretz et al., 1986; Jones and Palmer, 1987a; Hawken et al., 1988; Muly and Fitzpatrick, 1992; Chapman and Stryker, 1993;

Ringach, 2002), we are not sure whether this neural circuit is the most likely candidate because the fixed identical connection between LGN cells and the interneurons seems artificial unless they can be learned.

**Figure 14C** shows another possible neural circuit for implementing **A** −. LGN Cells are connected to interneurons via long-range excitatory connections, −**A** −. There is a one-to-one mapping between interneurons and cortical simple cells. In this case, the overall effect from LGN cells to simple cells is equivalent to **A** −. In addition, the RFs of inhibitory interneurons should resemble simple cells and show orientation tuning since the learned **A** − has spatial structures such as oriented bars, which is consistent with the smooth simple cells found in cat V1 of the experimental study (Hirsch et al., 2003). The positive connections −**A** − can be learned by Hebbian learning and the identical connections between interneurons and cortical simple cells can be learned by anti-Hebbian learning. Therefore, this neural circuit is more feasible than than the circuit in **Figure 14B**.

# 4.7. Discrepancies Between Model and Experimental Data

Our model can capture the most significant features of experimental phenomena such as the segregation of ON and OFF sub-regions, push-pull effect and contrast invariance of orientation tuning. However, there are also discrepancies between the distributions of model and experimental data. In general, the histograms of experimental data (**Figures 6C**, **7C**, **13C**) are wider than model data (**Figures 6D**, **7D**, **13D**), which shows that experimental data is more diverse. One possible explanation is that model cells in this paper are only a subset of the rich repository of real cortical cells. Furthermore, choices of free parameters in the model might also lead to different results.

## REFERENCES

# 5. CONCLUSION

In this paper, we presented a biologically plausible model of LGN-V1 pathways to account for many experimental phenomena of V1. We found that the segregation of ON/OFF sub-regions of simple cells, push-pull effect, and phase-reversed corticothalamic feedback can all be explained by the structure of learning connections when the model incorporates ON and OFF LGN cells and is trained using natural images. Furthermore, the model can produce diverse shapes of receptive fields and contrast invariance of orientation tuning of simple cells, consistent with experimental observations.

# DATA AVAILABILITY

The datasets generated for this study are available on request to the corresponding author.

# AUTHOR CONTRIBUTIONS

All authors listed have made a substantial, direct and intellectual contribution to the work, and approved it for publication.

### FUNDING

This research was funded by the Australian Research Council Discovery Projects scheme (Project DP140102947). HM acknowledges funding from the ARC Centre of Excellence for Integrative Brain Function (CE140100007).

# ACKNOWLEDGMENTS

We would like to thank Michael Ibbotson and Ali Almasi for helpful discussion and comments.

cortex of the old world monkey. J. Neurosci. 8, 3541–3548. doi: 10.1523/JNEUROSCI.08-10-03541.1988

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

Copyright © 2019 Lian, Grayden, Kameneva, Meffin and Burkitt. This is an openaccess article distributed under the terms of the Creative Commons Attribution License (CC BY). The use, distribution or reproduction in other forums is permitted, provided the original author(s) and the copyright owner(s) are credited and that the original publication in this journal is cited, in accordance with accepted academic practice. No use, distribution or reproduction is permitted which does not comply with these terms.

# rSK1 in Rat Neurons: A Controller of Membrane rSK2?

Eleonora Autuori, Petra Sedlak, Li Xu, Margreet C. Ridder, Angelo Tedoldi and Pankaj Sah\*

Queensland Brain Institute, The University of Queensland, Brisbane, QLD, Australia

In mammalian neurons, small conductance calcium-activated potassium channels (SK channels) are activated by calcium influx and contribute to the afterhyperpolarization (AHP) that follows action potentials. Three types of SK channel, SK1, SK2 and SK3 are recognized and encoded by separate genes that are widely expressed in overlapping distributions in the mammalian brain. Expression of the rat genes, rSK2 and rSK3 generates functional ion channels that traffic to the membrane as homomeric and heteromeric complexes. However, rSK1 is not trafficked to the plasma membrane, appears not to form functional channels, and the role of rSK1 in neurons is not clear. Here, we show that rSK1 co-assembles with rSK2. rSK1 is not trafficked to the membrane but is retained in a cytoplasmic compartment. When rSK2 is present, heteromeric rSK1-rSK2 channels are also retained in the cytosolic compartment, reducing the total SK channel content on the plasma membrane. Thus, rSK1 appears to act as chaperone for rSK2 channels and expression of rSK1 may control the level of functional SK current in rat neurons.

### Edited by:

George Augustine, Nanyang Technological University, Singapore

### Reviewed by:

Robert Brenner, The University of Texas Health Science Center at San Antonio, United States Paul F. Kramer, National Institute of Neurological Disorders and Stroke (NINDS), United States

### \*Correspondence:

Pankaj Sah<EMAIL_ADDRESS>
Received: 09 July 2018 Accepted: 11 March 2019 Published: 03 April 2019

### Citation:

Autuori E, Sedlak P, Xu L, Ridder MC, Tedoldi A and Sah P (2019) rSK1 in Rat Neurons: A Controller of Membrane rSK2? Front. Neural Circuits 13:21. doi: 10.3389/fncir.2019.00021 Keywords: spike frequency adaptation, potassium channel, afterhyperpolarization, excitability, calcium activated K+ channels (KCa1–KCa5)

### INTRODUCTION

Potassium channels are widely expressed in the central nervous system (CNS) where they play an important role in regulating the intrinsic excitability of neurons. A subset of potassium channels expressed in central neurons are Ca2+-activated K+ channels that regulate cellular excitability and spike frequency adaptation (Coetzee et al., 1999; Adelman et al., 2012). These are divided into three families that comprise the large conductance (BK) channels (KCa1.1), small conductance (SK) channels KCa2.1, KCa2.2, KCa2.3 (SK1, SK2, and SK3), and the intermediate conductance (IK) channels KCa3.1 (Vergara et al., 1998). In neurons, these channels are generally driven by calcium influx during action potentials, and activation of BK currents contributes to spike repolarization, while SK channel activity is slower, contributing to the afterhyperpolarization (AHP) that follows (Sah, 1996). These channels can also be activated by calcium release from intracellular stores again leading to inhibition of neural activity (Fiorillo and Williams, 1998), and are also expressed at glutamatergic synapses, where they are activated by calcium influx during synaptic activity and play a role in tuning synaptic plasticity (Faber et al., 2005; Ngo-Anh et al., 2005; Lin et al., 2008).

SK channels are encoded by three genes, SK1, SK2 and SK3 (Köhler et al., 1996; Stocker and Pedarzani, 2000; Sailer et al., 2002), which are expressed throughout the mammalian brain in distinct but partially overlapping distributions. In the rodent brain, SK1 and SK2 are generally co-expressed while SK3 channels are present in a complementary distribution (Stocker and Pedarzani, 2000). Functional studies in heterologous expression systems have shown that rat SK2 (rSK2) and SK3 (rSK3), and human SK1 (hSK1) channels form functional homomeric channels

**92**

that are voltage-insensitive and gated by the binding of Ca2+ to calmodulin (CaM; Xia et al., 1998; Lee and MacKinnon, 2018), which is covalently linked at their cytosolic, carboxyterminal region (Lee and MacKinnon, 2018). In contrast, the rat SK1 (rSK1), does not generate functional homomeric channels (Köhler et al., 1996; D'Hoedt et al., 2004), but does appear to co-assemble with rSK2 (Ishii et al., 1997; Benton et al., 2003; Church et al., 2015). This difference is due to differences in the sequence identity of the two channels. Thus, rSK2 and rSK3 are highly homologous to the human SK2 (97.6%) and SK3 (94.4%), and have similar pharmacology and functional profiles (Köhler et al., 1996; Joiner et al., 1997; Desai et al., 2000), while rSK1 is only 84% homologous to hSK1 (D'Hoedt et al., 2004). Indeed, replacing the carboxy terminal of rSK1 with that of rSK2 restores functional expression of rSK1, and swapping C- and N-termini of hSK1 with those from rSK1 prevents expression of functional hSK1 channels on the cell membrane of HEK293T cells (D'Hoedt et al., 2004).

In summary, while rodent rSK1, rSK2 and rSK3 channels are widely expressed in neurons of the rodent nervous system, rSK1, unlike rSK2 and rSK3, does not form functional homomeric channels. Homomers and dimers consisting of SK2 and SK3 form functional channels and contribute to the AHP and are also involved in synaptic plasticity (Adelman et al., 2012). It has been reported (Benton et al., 2003) that co-expression of rSK1 and rSK2 in HEK293 cells results in larger calcium-activated currents with altered pharmacology suggesting that rSK1 can form functional heteromeric assemblies with rSK2, though biochemical evidence for this co-assembly is lacking. Thus, while rSK1 is widely co-expressed with rSK2, the functional role of rodent SK1 remains unclear. In this study, we test the function of rSK1 channels in a heterologous system and cultured rat hippocampal pyramidal neurons. Firstly, we tested if rSK1 and rSK2 could interact in a heterologous system. We then looked at the effect of either overexpressing or downregulating rSK1 channels and tested the level of rSK2 channels expressed on the cells membranes of both a heterologous system and in hippocampal pyramidal neurons. Finally, we tested if changing rSK1 levels was related to a change in medium AHP (mAHP) and its underlying current (IAHP) in hippocampal pyramidal neurons that are generated by the activation of SK2 channels (Stocker et al., 1999, 2004).

## MATERIALS AND METHODS

### Neuronal Cultures

Primary E18 Wistar rat hippocampal cultures were prepared as previously described (Delaney et al., 2013). Briefly, hippocampi were digested for 20 min at 37◦C in papain (12 U/ml, Worthington, suspension 28.4 U/mg) made up in dissection medium (1× HBSS, 1% penicillin/streptomycin, 1% pyruvate, 10 mM HEPES, 30 mM glucose) and supplemented with 1% DNase (Sigma). The digested material was washed three times in plating medium [Neurobasal medium (Life Technologies) supplemented with 5% heat-inactivated fetal bovine serum (FBS), 1% penicillin/streptomycin, 1% Glutamax (Life Technologies) and 2% B27 (Life Technologies)] and the resulting pellet was triturated in plating medium using a fire-polished glass Pasteur pipette. Cells were plated at a density of 1 × 105 cells/ml for immunocytochemistry and the biotinylation assays and at 4 × 104 cells/ml for the electrophysiology experiments, on glass coverslips precoated with poly-D-lysine in plating medium and maintained at 37◦C in 5% CO2.

### Biotinylation and Western Blot Analysis

Hippocampal neurons were cooled on ice and were then rinsed three times with phosphate buffered saline (PBS) prior to exposure with EZ-Link Sulfo-NHS-SS-Biotin (Pierce), dissolved in PBS to a concentration of 1.22 mg/ml. The biotinylation reaction was undertaken on ice for 30 min, followed by biotin removal and washing twice with 100 mM glycine in PBS and then once with PBS. Cells were lysed in 200 µl lysis buffer (20 mM sodium phosphate, pH 7.5, 150 mM NaCl, 0.5% NP40, 0.5% sodium deoxycholate, 0.1% SDS, 1× complete protease inhibitors, Roche) and left on ice for 30–45 min. The lysates were then centrifuged at 6,500 g for 5 min. Twenty microliters of the resulting supernatant was then removed and this was classified as the total lysate (LYS). Biotinylated proteins were captured from the supernatant by the addition of 25 µl Streptavidin Agarose Resin (Pierce) at 4◦C for 30 min. The beads were then centrifuged at 10,000 g for 1 min and the resulting supernatant was classified as the cytoplasmic phase (CYTO). The beads were washed a further three times in PBS. This was classed as the membrane fraction (MEMB). All samples were then boiled for 5–10 min in 1× SDS sample buffer containing 100 mM DTT and then fractionated on 4%–12% Bis-Tris gradient gels (Life Technologies) and transferred to polyvinylidene difluoride membrane (Immobilon-P, Merck Millipore) at 150 V in 1× MOPS transfer buffer (Life Technologies). Blots were blocked in Tris-buffered saline (TBS) containing 5% skim milk powder, probed with primary antibody [mouse α Myc (9B11; 1/5,000, Cell Signalling Technology), rabbit α HA (1/500, Cell Signalling Technology), mouse α β-actin (AC-15; 1/20,000, Cell Signalling Technology), rabbit α EGFR (1/2,000, Cell Signalling Technology), mouse α Na+/K+ ATPase (alpha 1; clone C464.6; 1/4,000, Merck Millipore), rabbit α SK2 (c-39; 1/2,000, gift of J. Adelman)] followed by incubation with horseradish peroxidase-conjugated goat anti-rabbit or mouse IgG (1/20,000, Biorad) and detection by SuperSignal West Pico or FEMTO chemiluminescent substrate (Pierce). Blots were scanned using the Odyssey Infrared Imaging System (Li-cor) and densitometry analysis was carried out using Image Studio Lite software.

All data are expressed as mean ± standard error of the mean (SEM). Statistical analysis was performed using GraphPad Prism (GraphPad Software). All the data were tested for normal distribution using a normality test. If data were normally distributed, a Student's unpaired t-test was used. If the data did not pass the normality test, the Mann-Whitney test was used. Significance was determined at p < 0.05.

### Co-immunoprecipitation

Supernatants were precleared for 1 h at 4◦C with 75 µl of 50% Sepharose bead slurry (Amersham) and 1.5 µg of speciesspecific IgG (Sigma, 1 mg/ml). Protein G sepharose bead slurry and mouse IgG were used for the samples immunoprecipitated with mouse α Myc (9B11; Cell Signalling Technology); protein A sepharose bead slurry and rabbit IgG for the samples immunoprecipitated with the rabbit α HA (Cell Signalling Technology) antibody. Samples were then spun at 13,000 rpm for 5 min and the resulting supernatants were incubated overnight at 4◦C respectively with mouse α Myc antibody (2 µl/1.5 mg proteins) or rabbit α HA (1:50) antibodies. The following day, the α Myc samples were incubated with 75 µl of 50% protein G sepharose beads and the α HA samples with 75 µl of 50% protein A sepharose beads for a further 3 h at 4◦C. Samples were centrifuged at 3,000 rpm for 5 min at 4◦C and the recovered pellets were washed three times in RIPA buffer [150 mM NaCl, 1% Triton X-100, 0.5% sodium deoxycholate, 0.1% SDS, 50 mM Tris pH = 8.0 and protease inhibitors (Roche)]. Samples were eluted with sample buffer containing DTT and denatured by boiling for 5 min. Co-immunoprecipitation samples (IP) were run together with the input (IN) and the supernatant (S) resulting from the immunoprecipitation, transferred as described above and blotted with mouse α Myc and rabbit α HA as described above.

### Lentiviral Constructs

Several knock-down hairpin constructs were produced using annealed phosphorylated oligonucleotides, which were ligated into the pll3.7 dsRed vector digested with HpaI and XhoI. The pll3.7 dsRed vector has the EGFP of pll3.7 vector replaced with dsRed. HEK293T cells were plated at a density of 5 × 105 cells per 35 mm well and were transfected 2 h later using Lipofectamine 2,000 (Life Technologies) as per the manufacturer's instructions with constructs expressing rat SK1-YFP (gift of G. Moss) and three different plasmids carrying SK1 knock down hairpins (rSK1-KD 2.4, 3.6 and 4.5) at a ratio of 1:1. Forty-eight hours later, the cells were checked for YFP and dsRed expression, washed with PBS and lysed in 600 µl 1× sample buffer. Thirty microliters of this was run on an SDS-PAGE gel, transferred and western blots carried out using rabbit α GFP (1/500, Merck Millipore) primary antibody as per methods above except that x-ray film was used for detection. Mouse β-actin (AC-15; 1/20,000, Cell Signalling Technology) was used as a loading control. Out of the three plasmids, rSK1-KD 3.6 was the most successful at reducing SK1-YFP expression (**Figure 1**). The primers used to produce this knockdown construct are as follows–

To enable us to determine the virally infected neurons, we replaced the CMV promoter (expressing the dsRed) in pll3.7 dsRed with the hSyn promoter using the Not and Nhe restriction sites and the following primers–

A scrambled (SCRAM) control plasmid was also produced using a random combination of the hairpin sequence of the successful knockdown construct. The primers used to produce the SCRAM construct are as follows–

The rat SK2 overexpression plasmid pJPA5.rSK2.Myc3 (gift of J. Adelman) has a triple Myc tag located extracellularly between the S3 and S4 domains. The Myc-tagged SK2 was further cloned into the lentiviral FUGW plasmid (gift of P. Osten) using the AgeI and EcoRI restriction sites and the following primers–

The rat SK1 overexpression plasmid pcDNA3.1-rSK1\_HA (gift of G. Moss) has a HA tag located intracellularly between the S2 and S3 domains. The HA-tagged SK1 was further cloned into the lentiviral FUGW plasmid (gift of P. Osten).

For lentivirus production, all plasmids (pMDG, pMDL g/p RRE, pRSV Rev and the pll3.7 transfer vector) were prepared using the Qiagen Endofree maxiprep kit and lentivirus was prepared via calcium phosphate transfection of 80% confluent triple T175 flasks of HEK293T cells. Cells were grown in DMEM (Life Technologies) plus 10% FBS and 1% penicillin/streptomycin (Life Technologies). Forty microgram of transfer vector together with 20 µg of the other three plasmids were transfected per triple T175 flask. Plasmids were diluted to 0.5 µg/ml in TE buffer and made up to a volume of 3 ml with sterile distilled water. Three-hundred microliter of 2.5 M CaCl2 was added together with 3 ml of 2× BBS (50 mM BES, 280 mM NaCl, 1.5 mM Na2HPO4, pH = 6.95). This was incubated for 20 min at room temperature before adding the mix dropwise to the cells. The cells were incubated for 3–4 h at 37◦C plus 5% CO2 before the removal of the calcium phosphate mix and addition of new media. Viral supernatants were removed 48 h post-transfection, filtered using a 0.45 µm filter unit and then centrifuged for 4 h at 20,000 g over a 20% sucrose cushion. The resulting pellet was resuspended overnight in PBS containing 1% BSA and then ultracentrifuged at 44,000 g for 90 min to concentrate the virus. The virus pellet was resuspended in PBS containing 1% BSA.

Neuronal cultures were infected at 4–5 DIV and the virus was left on for 4–12 h. The virus was then removed and conditioned media added to the cultures. Biotinylation experiments and quantitative PCR (qPCR) were undertaken 4–5 days post infection. Knockdown and overexpression lentiviral constructs were also tested for specificity and expression using qPCR. Briefly, RNAs were purified using a RNeasy kit (Qiagen), DNase 1 (Qiagen) treated and reverse transcribed to produce complementary DNA using random hexamers and the SuperScript III First Strand Synthesis System (Life Technologies), as per the manufacturer's instructions. Standard qPCR was carried out on a Rotorgene RG-3000 Thermocycler using the Platinum SYBR Green qPCR UDG SuperMix (Life Technologies) and 0.2 µM of the following primers:

Reactions were performed in either duplicate or triplicate. The PCR program was as follows: 50◦C 2 min, 95◦C 10 min, 35 cycles −95◦C 10 s, 54◦C 15 s, 72◦C 20 s. Relative expression of rSK1 and rSK2 mRNA was obtained using the comparative Ct method.

### Cosm6 Transfections

For the Cosm6 transfections- pJPA5.rSK2.Myc3 and pcDNA3.1 rSK1\_HA along with pcDNA3 for the single construct transfections were transfected using the same protocol as described above but using a ratio of 1:2:1.

### Immunohistochemistry

Hippocampal cultures were washed in PBS and fixed in 4% PFA containing 4% sucrose. After three washes in PBS, the cells were first blocked in a 0.3% BSA/0.1% Triton X-100 solution, then incubated in primary antibody overnight at 4◦C, washed in PBS and left in secondary antibody (Life Technologies) for a further 1 h at room temperature. Stained cultures were analyzed on a confocal fluorescence microscope (Zeiss LSM510).

### Electrophysiology

Rat hippocampal neuronal cultures were recorded at least 2 weeks after viral infection (DIV 15–17). Coverslips were transferred from the incubator to the recording chamber, located on a stage of a fluorescence microscope (Zeiss 710). Ringer solution (140 mM NaCl, 5 mM KCl, 2 mM CaCl2, 1 mM MgCl2 10 mM glucose and 10 mM HEPES) was constantly perfused at 1–1.5 ml/min at room temperature. Whole-cell recordings were obtained from dsRed-positive neurons (SCRAM or SK1 KD) or neurons (rSK1-HA and rSK2-myc) using a K-methyl-sulfate based internal solution containing (mM): 135 KMeSO4, 5 NaCl, 10 Hepes, 2 Mg2-ATP, 0.3 Na3-GTP, 0.1 spermine and 7 phosphocreatine (pH = 7.3, ∼290 mOsmol).

Tetrodotoxin (TTX, 1 µM, Sigma), D-(−)-2-Amino-5 phosphonopentanoic acid (D-APV, 20 µM, Tocris) and 6-cyano-7-nitroquinoxaline-2,3-dione disodium salt (CNQX, 20 µM, Tocris) were added to the ringer solution to reduce spontaneous activity and isolate mAHP and IAHP. For a recording of mAHP, cells were injected with a 400 pA/100 ms step either at resting membrane potential or by clamping the cells at −50 mV. For recordings of IAHP, cells were held at −50 mV and a depolarizing step (+60 mV) to +10 mV for 100 ms was used to increase Ca2+ influx through voltage-gated Ca2+channels, which are necessary to activate IAHP (Stocker et al., 2004).

Data were collected using Axograph X software and a Multiclamp 700B amplifier (Molecular Devices). Signals were filtered at 10 kHz and digitized at 50 kHz using an ITC-16 A/D converter (InstruTech).

### RESULTS

### rSK1 and rSK2 Co-assemble in Cosm6 Cells

To test expression of SK channels, rSK1 and rSK2 were epitopetagged with HA and Myc, respectively (see ''Materials and Methods'' section). As shown previously (D'Hoedt et al., 2004; Church et al., 2015), transfection in Cosm6 cells shows that rSK2 was expressed throughout the cell, while rSK1-HA was restricted to the cytosolic somatic compartment (**Figure 2**), likely the endoplasmic reticulum and Golgi (Church et al., 2015). As shown by others (Strassmaier et al., 2005), western blot analysis using Myc and HA antibodies shows bands of several sizes that correspond to monomeric protein (∼60 kDa, band 3), dimers (∼120 kDa, band 2) and a large molecular band (∼200 kDa, band 1), which is likely to be a tetramer (**Figure 2B**). Separation of the membrane fraction (MEMB on blots) using surface biotinylation (see ''Materials and Methods'' section) shows that unlike rSK2, rSK1 channels are not detectable on the plasma membrane. When rSK1 and rSK2 were co-expressed, again rSK2 was detectable in the membrane fraction but rSK1 was not (**Figure 2B**), but the level of rSK2 protein in the membrane fraction was reduced (**Figures 2B,C**). There was a clear reduction

FIGURE 2 | rSK1 regulates the membrane expression of rSK2 channels. (A) Immunocytochemistry of Cosm6 cells transfected with rSK2-Myc or rSK1-HA. rSK2-Myc (pink), rSK1-HA (green), DAPI (blue). (B) Immunoblots from Cosm6 cells transfected with rSK1-HA, rSK2-Myc and hSK1 in different combinations. Shown are the total lysate (LYS) and membrane (MEMB) fraction isolated using surface biotinylation probed with either Myc (rSK2; left) or HA (rSK1; right) antibody. EGFR was used as the loading control. rSK2 transfection results in three size bands that represent monomeric protein (band 3, ∼60 kDa), dimers (band 2, ∼120 kDa) and tetramers (band 1, >200 kDa). Co-transfection of rSK2 and rSK1 reduced the amount of rSK2 expressed in the cell membrane. (C) Bar graphs show normalized relative density of SK2, as described in the "Materials and Methods" section, in the membrane fraction from Cosm6 cells transfected with rSK2 alone or co-transfected with rSK1 or hSK1. Co-transfection with rSK1-HA significantly reduced the 60 kDa band (p = 0.002). Co-transfection with hSK1 increased expression of rSK2. (D) Co-immunoprecipitation (co-IP) assay of transfected Cosm6 cells shows co-assembly of rSK1-rSK2. Solubilized protein from Cosm6 cells transfected with rSK1 and/or rSK2 was immunoprecipitated with Myc (left) or HA (right) and immunoblotted with either Myc or HA as indicated.

in monomeric rSK2 (n = 3, p = 0.002), but other higher molecular weight species were also greatly reduced (**Figure 2C**). In contrast to rSK1, hSK1 is known to express as functional channels (D'Hoedt et al., 2004), and forms heteromultimeric channels (Ishii et al., 1997; Benton et al., 2003; Church et al., 2015). In agreement, co-expression of hSK1 with the rSK2 increased the amount of rSK2-Myc on the cell membranes (**Figures 2B,C**).

What explains the reduction in membrane rSK2 when it is co-expressed with rSK1? As rSK1 is retained in cytosolic compartments and not trafficked to the membrane (Benton et al., 2003; Church et al., 2015), it is possible that rSK2 co-assembles with rSK1 with the resulting heteromultimers being retained, leading to a reduction of membrane rSK2. We, therefore, tested for interactions between these two channels using co-immunoprecipitation (co-IP). rSK1-HA and rSK2-Myc were expressed in Cosm6 cells, and Myc or HA antibodies used to isolate rSK2 or rSK1 using sepharose A/G beads followed by immunoblot analysis of SK protein. IP of Cosm6 cell homogenates using anti-Myc antibody revealed the presence of rSK2 in both supernatant, input and co-immunoprecipitated samples, from single transfected (rSK2-Myc) and co-transfected cells (rSK2-Myc+rSK1-HA; **Figure 2D**). Probing with HA antibody revealed the presence of rSK1-HA in the multimeric complex (**Figure 2D**, left). Complimentary experiments using the anti-HA antibody to precipitate rSK1 protein lysates and probing with anti-Myc also revealed the presence of rSK2 in the multimeric high molecular weight band (**Figure 2D**, right). It was notable that following IP with myc and probe with HA or the reverse, while very high molecular weight (>200 kD) complexes are present, there is a dearth of complexes at ∼110 kD (**Figure 2D**, upper panels). We interpret this to show the presence of tetrameric heteromeric protein.

# Overexpression of rSK1 in Neurons Decreases Plasma Membrane rSK2

These results show that when expressed in Cosm6 cells rSK1 and rSK2 can co-assemble. Moreover, the interaction between the two channels affects the total rSK2 present in the cell membrane. Thus, rSK1 channels appear to ''trap'' rSK2 in the cytoplasmic compartment, therefore reducing the total amount of rSK2 trafficked to the cell membrane.

To test if rSK1 acts similarly in neurons, we turned to cultured rat hippocampal neurons. rSK1-HA and rSK2-Myc were delivered using lentivirus (see ''Materials and Methods'' section). We first tested the effect of overexpression of rSK1 or rSK2 on mRNA levels, using qPCR. Increasing the volume of rSK1 or rSK2 expressing virus added to neuronal cultures increased the relative expression of mRNA of the specific gene delivered without affecting the other: infection of 1 × 105 rat hippocampal cultures with 0.2 µl and 0.5 µl rSK1-HA lentivirus increased rSK1 relative expression by 60 and 200-fold, respectively, while infection with 0.2 µl and 0.5 µl of rSK2-Myc lentivirus increased rSK2 relative expression by 10 and 40-fold, respectively (data not shown). Next, expression of SK2 protein was tested using immunocytochemistry (**Figure 3**). Transduction of rSK2-Myc (**Figure 3A**) or rSK1-HA (**Figure 3B**) showed clear expression of Myc and HA in cells positive for the neuronal marker NeuN. Co-expression of rSK1 and rSK2 show that the two are at least partially co-localized in neurons (**Figure 3C**), consistent with heterodimer formation of rSK1 and rSK2. It can also be seen that in neurons when rSK1-HA is expressed alone (**Figure 3B**), its distribution is not confined to the soma, but spreads throughout the cell, consistent with the expression of endogenous SK2 in these neurons (**Figure 3D**).

To test for the presence of rSK2 in the plasma membrane we again turned to the surface biotinylation assay to separate the membrane compartment. rSK1-HA and rSK2-Myc were transduced using lentivirus. As in Cosm6 cells (**Figure 2**), we found that co-expression of rSK2 with rSK1 decreased the amount of rSK2 detected in the membrane fraction of infected neurons (**Figures 4A,B**). This reduction was particularly significant in the high molecular weight band, again suggesting delivery of a heteromultimeric protein (n = 3, p = 0.04; **Figure 4B**). To account for possible effects of overexpressing rSK2, we tested the impact of rSK1 expression on endogenous rSK2 using the specific SK2 antibody (c-39). When rSK1-HA was transduced in hippocampal neurons, the endogenous levels of rSK2 protein detected in the cell membrane was reduced (**Figure 4C**, upper panel), while the cytoplasmic fraction increased (**Figure 4C**, bottom panel). Quantification of the relative density of endogenous rSK2 (c-39) normalized to β-actin shows that the density of the high molecular band was significantly reduced in the membrane fraction in comparison to uninfected neurons (n = 5, p = 0.001; **Figure 4D**), but was increased in the cytosolic fraction (n = 4, p = 0.006; **Figure 4D**, bottom panel). This increase in SK2 in the ''cytoplasmic'' fraction is likely due to the fact that rSK1-rSK2 heteromultimers are not trafficked to the plasma membrane, but remain trapped in the cytoplasmic compartment, possibly the endoplasmic reticulum.

### rSK1 Modulates Functional SK Channels

In mammalian neurons, SK channels contribute to the medium duration AHP that follows action potentials (IAHP), and in some types of neurons they are also present at excitatory synapses where they modulate synaptic strength (Faber et al., 2005; Ngo-Anh et al., 2005; Adelman et al., 2012). We have shown that overexpression of rSK1 seems to reduce the amount of membrane rSK2 protein. We next tested if altering endogenous rSK1 levels would have an impact on functional membrane rSK channels in neurons. To disrupt endogenous rSK1, we used short hairpin RNA interference (RNAi) constructs to reduce rSK1 expression in hippocampal neurons using lentivirus (see ''Materials and Methods'' section). rSK2 protein levels in the membrane were assessed using the biotinylation assay and western blot analysis (**Figures 5A,B**). Knockdown of rSK1 (rSK1-KD) significantly reduced the amount of rSK2 monomers (60 kDa) expressed in the cell membrane in comparison to neurons infected with scrambled virus (n = 4, p = 0.008), but the total membrane rSK2 tetramers (>200 kDa) appeared to be higher than the scrambled control (**Figure 5A2**). In the cytosolic compartment, high molecular weight (>200 kDa) rSK2 protein was significantly reduced as compared to scrambled control cells (n = 4, p = 0.03; **Figure 5B2**).

To test for functional expression of rSK channels, whole-cell recordings were obtained from infected hippocampal neurons in culture (**Figures 5C,D**). The passive membrane properties of neurons in infected with SK2 RNAi hairpins and scrambled controls were not different, and are given in **Supplementary Table S1**. Reducing rSK1 with RNAi increased the amplitude of IAHP in comparison to that in scrambled infected neurons (p = 0.03; **Figure 5C**). As expected, switching to current clamp also revealed an increase in the amplitude of the medium AHP (p = 0.03; **Figure 5D**), and this enhanced IAHP current was blocked by apamin (**Figure 5E**). Together, these results support a role for rSK1 channels in regulating the expression of rSK2 on the cell membrane.

### DISCUSSION

SK1, SK2 and SK3 encode calcium-activated potassium channels that are widely expressed in the mammalian brain. In

FIGURE 4 | rSK1 regulates membrane expression of rSK2 channels in cultured neurons. Cultured hippocampal neurons were transduced to express rSK1-HA or rSK2-Myc or both. Surface protein was biotinylated, and solubilized protein from neurons separated into total lysate (LYS) and membrane fraction (MEMB) as detailed in the "Materials and Methods" section. Immunoblots were probed using the Myc antibody to test for rSK2 expression. Na/K ATPase or β-actin was used as loading controls. (A) Co-expression of rSK1 and rSK2 reduces membrane expression of rSK2. (B) Graph shows normalized relative density of SK2, in the membrane fraction from neurons transduced with rSK2 alone or rSK2 and rSK1. Co-infection significantly reduced the band at 200 kDa (p = 0.04). (C) Expression of rSK1 in neurons reduced surface expression of endogenous rSK2. Cultured neurons were transduced to express rSK1 using lentivirus. Blots show solubilized protein from the membrane fraction (upper blots) and cytosolic fraction (lower blots) from control neurons (left) and neurons transduced with rSK1(right). Blots were probed with anti SK2 antibody (c-39). (D) Quantification of SK2 protein expression in the membrane (upper graph) and cytosolic (lower graph) fraction from control neurons and neurons transduced with rSK1. Neurons infected with rSK1 had a significantly lower amount of endogenous high molecular weight rSK2 in the membrane fraction and this was balanced by an increase in the total cytoplasmic fraction.

rodents, rSK1 and rSK2 are expressed in strongly overlapping distributions. However, while expression of rSK2 channels produces functional calcium-activated potassium channels, rSK1 channels are made but not trafficked to the plasma membrane (D'Hoedt et al., 2004; Church et al., 2015), and the functional role of these channels is not clear. In this study, we show that rSK1 channels co-assemble with rSK2, and regulate the plasma membrane levels of rSK2. Thus, expression of rSK1 in

and membrane fraction (MEMB) as detailed in the "Materials and Methods" section. Blots were immunoprobed using the rSK2 antibody (c-39). (A1) Knockdown of SK1 protein results in an overall reduction of total endogenous monomeric SK2 protein, but enhancement of high molecular weight SK channels in the plasma membrane. Na+/K+ ATPase was sued as the loading control. (A2) Normalized relative density of endogenous rSK2 in the membrane fraction [NORM (MEMB)] from neurons infected with rSK1-KD. rSK1-KD significantly decreased the amount of endogenous monomeric rSK2, but increased high molecular weight rSK2. Na+/K+ ATPase (∼110 kDa) was used to normalize all the loaded samples for the MEMB fraction. (B1) Knockdown of rSK1 results in a reduction in high molecular weight SK2 in the cytoplasmic fraction. (B2) Quantification of blot shown in B1. β-actin (∼40 kDa) was used to normalize all the loaded samples for the CYTO fraction. (C) Knockdown of rSK1 in hippocampal neurons (red traces) increased the amplitude of the IAHP in voltage clamp (left) and the medium AHP in current clamp (right) in comparison to SCRAM infected neurons (traces in gray). (D) Quantification of the impact of the increase in the IAHP current shown as pA/pF and the medium AHP following knockdown of rSK1. (E) The enhanced IAHP current following knockdown of rSK1 is apamin (100 nM) sensitive consistent with SK2 channels.

neurons leads to the formation of heteromultimers containing rSK2 that are trapped in the cytoplasmic compartment and a reduction in the total rSK2 in the plasma membrane. In support of this, knockdown of endogenous rSK1 in cultured neurons results in an increase in membrane rSK2. These results suggest that rSK1 channels, rather than functioning as ion channels, are involved in regulating the expression of rSK2 channels in the plasma membrane.

While SK channels were cloned more than 20 years ago, antibodies for immunoprecipitation and immunohistochemistry for all three members are not easily available, thus we have used epitope-tagged channels to study their expression and trafficking. Here using rSK1 tagged with HA (rSK1-HA) and rSK2 tagged with Myc (rSK2-Myc), in agreement with previous studies (D'Hoedt et al., 2004; Church et al., 2015), we show that when expressed in Cosm6 cells, rSK2 channels are present in the plasma membrane, rSK1 channels are made but not trafficked from the cytoplasmic compartment, with little or no protein detectable in the membrane compartment. When rSK1 and rSK2 are co-expressed in Cosm6 cells, we show using co-immunoprecipitation that these channels co-assemble and there is a significant reduction in the amount of rSK2 in the plasma membrane. It is well known that hSK1 channels can form functional channels that are expressed on the cell surface (D'Hoedt et al., 2004; Church et al., 2015). In agreement with this, co-expression of rSK2 with hSK1 had an opposite effect to that of rSK1 with an increase in membrane rSK2.

These results in Cosm6 cells were confirmed in rat hippocampal neurons. Co-expression of rSK1-HA and rSK2-Myc resulted in an overall reduction in the amount of rSK2-Myc trafficked to the membrane. Importantly, expression of exogenous rSK1-HA in cultured neurons also resulted in a reduction in the total endogenous rSK2, showing that expression of rSK1 can modify the levels of membrane rSK2 channel. Finally, we show that by reducing the amount of rSK1 channels using RNAi, increases the amount of endogenous rSK2 channels expressed in the cell membrane, and is consistent with the finding that transgenic mice lacking SK1 show no overall effects on hippocampal neuronal electrophysiology (Bond et al., 2004). We show that reducing the rSK1 content of hippocampal neurons reduces the overall amount of rSK2 monomeric protein in the membrane but increases the total large molecular weight fraction (**Figure 5**). This is perhaps expected as we have not modified the expression level of rSK2 protein, and as rSK1 levels are lower, more channels are trafficked to the membrane where they co-assemble as high molecular weight homomeric multimers. As a result, total monomeric protein levels are lower. This increase in membrane rSK2 levels has a functional impact as it also increases the amplitude of the SK-mediated IAHP in infected rat neurons. Thus, by increasing or reducing the amount of rSK1 channels expressed in a heterologous system or rat neurons, we can decrease or increase rSK2 channels expressed on the cell membrane.

Previous studies have shown that when expressed in HEK293 cells, the interaction of rSK1 and rSK2 channels resulted in an overall larger SK-channels mediated current and a change in their pharmacology (Benton et al., 2003; Church et al., 2015). This appears to be in contrast with our results in both Cosm6 cells and rat neurons where overexpression of rSK1-HA reduced the amount of rSK2 channels expressed in the cell membrane. However, it remains possible that these channels assemble as heteromultimers with the total fraction of SK1-SK2 heteromeric protein in the membrane being very low.

Our results suggest that in the rodent brain rSK1 channels, rather than acting as independent membrane ion channels, are involved in trafficking of rSK2 channels to the plasma membrane. hSK1, which is ∼84% homologous to the rodent channels (D'Hoedt et al., 2004), behaves entirely differently, being translated as a functional ion channel (Church et al., 2015). The reason for this difference is not clear. It is possible, that the rSK1-rSK2 interaction and the regulation of rSK2 channels expressed on the cell membrane is important during development. rSK1 and rSK2 channel transcripts within the hippocampus (CA1, CA3 and dentate gyrus) start being expressed at embryonic day 19 (E19) in rodents (Gymnopoulos et al., 2014), with their expression patterns being very similar, showing colocalization (Stocker and Pedarzani, 2000). These similar expression patterns could indicate a developmental change within the expression of rSK2 channels on the cell membrane that could affect neuronal physiology. Interestingly, rSK2 was found to be present mainly in the ER of CA1 pyramidal neurons at P5. By P30, however, most rSK2 was present at spines and dendrites changing the cellular physiology (Ballesteros-Merino et al., 2012). It is, therefore, possible that rSK1 channels modulate the amount of rSK2 channels expressed on the cell membrane. While rSK2 channels are trapped in the ER by rSK1 channels, neurons can receive more inputs and increased the excitation in the early developmental stages, which will increase the memory acquisition. Once later developmental stages are reached, however, rSK1 channels may have a different role in regulating the expression of rSK2 channels and the neuronal physiology.

# ETHICS STATEMENT

This was approved by the University of QLD Ethics Committee.

# AUTHOR CONTRIBUTIONS

EA discussed experiments, did experiments, and wrote the manuscript. PSe did experiments and wrote the manuscript. LX made virus. MR did experiments and made figures. AT did experiments, made figures and wrote the manuscript. PSa designed experiments, made figures and wrote the manuscript.

# FUNDING

This work was supported by grants from the National Health and Medical Research Council of Australia and the Australian Research Council (CE140100007) to PSa.

### ACKNOWLEDGMENTS

We thank Victor Anggono for comments on the manuscript. We thank John Adelman and Guy Moss for supplying constructs and antibodies as indicated in the ''Materials and Methods'' section.

### REFERENCES

### SUPPLEMENTARY MATERIAL

The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fncir. 2019.00021/full#supplementary-material

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

Copyright © 2019 Autuori, Sedlak, Xu, Ridder, Tedoldi and Sah. This is an open-access article distributed under the terms of the Creative Commons Attribution License (CC BY). The use, distribution or reproduction in other forums is permitted, provided the original author(s) and the copyright owner(s) are credited and that the original publication in this journal is cited, in accordance with accepted academic practice. No use, distribution or reproduction is permitted which does not comply with these terms.

# Hierarchical and Nonlinear Dynamics in Prefrontal Cortex Regulate the Precision of Perceptual Beliefs

Leonardo L. Gollo1,2 \* † , Muhsin Karim3,4† , Justin A. Harris 5 , John W. Morley 6 and Michael Breakspear 1,2,3,4,7,8 \*

1QIMR Berghofer Medical Research Institute, Brisbane, QLD, Australia, 2Centre of Excellence for Integrative Brain Function, QIMR Berghofer Medical Research Institute, Brisbane, QLD, Australia, 3School of Psychiatry, Faculty of Medicine, University of New South Wales, Sydney, NSW, Australia, 4The Black Dog Institute, Sydney, NSW, Australia, 5School of Psychology, The University of Sydney, Sydney, NSW, Australia, 6School of Medicine, Western Sydney University, Sydney, NSW, Australia, 7Metro North Mental Health Service, Brisbane, QLD, Australia, 8Hunter Medical Research Institute, University of Newcastle, New Lambton Heights, NSW, Australia

### Edited by:

Gary F. Egan, Monash University, Australia

### Reviewed by:

Dirk Ostwald, Freie Universität Berlin, Germany Joachim Lange, Heinrich Heine Universität Düsseldorf, Germany

### \*Correspondence:

Leonardo L. Gollo<EMAIL_ADDRESS>Michael Breakspear michael.breakspear@ newcastle.edu.au

†These authors have contributed equally to this work

> Received: 31 July 2018 Accepted: 29 March 2019 Published: 24 April 2019

### Citation:

Gollo LL, Karim M, Harris JA, Morley JW and Breakspear M (2019) Hierarchical and Nonlinear Dynamics in Prefrontal Cortex Regulate the Precision of Perceptual Beliefs. Front. Neural Circuits 13:27. doi: 10.3389/fncir.2019.00027 Actions are shaped not only by the content of our percepts but also by our confidence in them. To study the cortical representation of perceptual precision in decision making, we acquired functional imaging data whilst participants performed two vibrotactile forcedchoice discrimination tasks: a fast-slow judgment, and a same-different judgment. The first task requires a comparison of the perceived vibrotactile frequencies to decide which one is faster. However, the second task requires that the estimated difference between those frequencies is weighed against the precision of each percept—if both stimuli are very precisely perceived, then any slight difference is more likely to be identified than if the percepts are uncertain. We additionally presented either pure sinusoidal or temporally degraded "noisy" stimuli, whose frequency/period differed slightly from cycle to cycle. In this way, we were able to manipulate the perceptual precision. We report a constellation of cortical regions in the rostral prefrontal cortex (PFC), dorsolateral PFC (DLPFC) and superior frontal gyrus (SFG) associated with the perception of stimulus difference, the presence of stimulus noise and the interaction between these factors. Dynamic causal modeling (DCM) of these data suggested a nonlinear, hierarchical model, whereby activity in the rostral PFC (evoked by the presence of stimulus noise) mutually interacts with activity in the DLPFC (evoked by stimulus differences). This model of effective connectivity outperformed competing models with serial and parallel interactions, hence providing a unique insight into the hierarchical architecture underlying the representation and appraisal of perceptual belief and precision in the PFC.

Keywords: decision making, dynamic causal modeling, fMRI, prefrontal cortex, vibrotactile

# INTRODUCTION

Percepts underpin all our interactions with the world. Perceptual precision, the confidence with which we hold those percepts, informs this interaction, such as when a decision is biased toward a precisely represented percept (Ernst and Banks, 2002). Although high perceptual precision may be advantageous in some contexts, such as when driving a car, there exist other situations where a degree of imprecision is crucial: if percepts were held with infinite precision then it would be impossible to recognize any object encountered for a second time. For example, the texture of a surface would feel unique and surprising on every touch. Whereas the neurobiology of perception has been a long-studied subject, research into the basis of perceptual precision and its impact on decision making has been more recent (Knill and Pouget, 2004; Moran et al., 2013; Pouget et al., 2013; Navajas et al., 2017).

The neural basis of perceptual decision-making has been extensively studied using two-alternative forced-choice tasks in the somatosensory (Romo and Salinas, 2003) and visual domain (Britten et al., 1992). These prototypical experiments consist in presenting two sequential stimuli that are followed by a forced response between two choices involving a comparison between the properties of these two stimuli (see **Figure 1**). In the somatosensory modality, a wealth of neurophysiological research using vibrotactile stimuli has established the crucial role of the prefrontal cortex (PFC) during the performance of such tasks (Gold and Shadlen, 2007; Hegner et al., 2007; Heekeren et al., 2008; Wang, 2012). While the primary somatosensory cortex is clearly involved in stimulus representation (Hernández et al., 2000; Harris et al., 2002; Sörös et al., 2007), the PFC holds the representation of the first stimulus in working memory for subsequent comparison against representation of the second stimulus (Preuschhof et al., 2006; Wang, 2008), as well as the final decision process (Miller et al., 2003; Pleger et al., 2006; Heekeren et al., 2008; Wang, 2008; Barak et al., 2010). With very few exceptions (Engel and Wang, 2011), decisions in these forced-choice experiments are only dependent on magnitude comparisons of the perceived frequencies. A sensory percept can be viewed probabilistically (as a probability distribution) and to first order can hence be decomposed into its magnitude (here, the perceived frequency) and its precision (the inverse of the variance of the probability distribution; see **Figure 2**). Whilst perceptual precision—classically captured by the signal-to-noise ratio—impacts upon the performance accuracy of a fasterslower comparison, the decision itself does not explicitly require representing and acting on the precision of those perceptions. This is because the final decision only rests upon deciding whether the second stimulus is faster or slower than the first and does not depend upon the subjective confidence in that judgment. That is, a faster-slower decision can be made by a simple subtraction and does not crucially depend upon the precision of either percept.

The anterior cingulate and ventromedial PFC appear to play critical roles in assessing the value of current information in an environment of uncertain outcome and reward (Daw et al., 2005; Kennerley et al., 2006; Behrens et al., 2007). These regions also represent changes in this value (that is, when the link between stimulus, outcome and reward is volatile; Rushworth and Behrens, 2008). Whilst the value of the percept to an external reward is uncertain in these studies (Fiorillo et al., 2003; Yu and Dayan, 2005; Hsu et al., 2005; Huettel et al., 2006; Behrens et al., 2007; Tobler et al., 2007), the percept itself is not ambiguous. Hence, it is not clear from these studies whether these regions are also involved in representing the intrinsic precision of the percept itself, or whether other regions are recruited when the stimulus is noisy but the task contingencies are fixed (Kayser et al., 2010; Bach and Dolan, 2012).

Here, we sought to disentangle the representation of stimulus properties from the precision of those representations in the PFC. Functional neuroimaging data were acquired whilst paired vibrotactile flutter stimuli (10–50 Hz) were sequentially applied to the index finger. In separate tasks, participants were requested to decide if the second stimulus was faster than the first, or if the second stimulus was different from the first. As rehearsed above, the ''faster-slower'' task can be performed by simply encoding and subtracting an estimate of each stimulus frequency—that is, decisions only explicitly depend on comparing the likely value of each of the flutter frequencies. In the ''same-different'' task, the magnitude of this subtraction must be weighed against the precision of the perceptual beliefs, such that a difference that is perceived as small may be inferred as significant if each percept is held precisely (and conversely for imprecise representations). The precision of a percept is the composite of the roughness of the stimulus and the perceptual imprecision due to stochastic effects in perceptual systems: to manipulate stimulus precision, noise was introduced to the vibrotactile oscillatory frequency as an additional experimental factor (Harris, 2006; Harris et al., 2006; Karim et al., 2012). Note that we refer to precision in the statistical sense of the inverse of the noise variance (**Figure 2**).