extfieldsofstudy:
Psychology
Medicine
Computer Science
Biology

text:
Learning Mid-Level Auditory Codes from Natural Sound Statistics

Interaction with the world requires an organism to transform sensory signals into representations in which behaviorally meaningful properties of the environment are made explicit. These representations are derived through cascades of neuronal processing stages in which neurons at each stage recode the output of preceding stages. Explanations of sensory coding may thus involve understanding how low-level patterns are combined into more complex structures. Although models exist in the visual domain to explain how mid-level features such as junctions and curves might be derived from oriented filters in early visual cortex, little is known about analogous grouping principles for mid-level auditory representations. We propose a hierarchical generative model of natural sounds that learns combinations of spectrotemporal features from natural stimulus statistics. In the first layer the model forms a sparse convolutional code of spectrograms using a dictionary of learned spectrotemporal kernels. To generalize from specific kernel activation patterns, the second layer encodes patterns of time-varying magnitude of multiple first layer coefficients. Because second-layer features are sensitive to combinations of spectrotemporal features, the representation they support encodes more complex acoustic patterns than the first layer. When trained on corpora of speech and environmental sounds, some second-layer units learned to group spectrotemporal features that occur together in natural sounds. Others instantiate opponency between dissimilar sets of spectrotemporal features. Such groupings might be instantiated by neurons in the auditory cortex, providing a hypothesis for mid-level neuronal computation.

Introduction
The ability to interact with the environment requires an organism to infer characteristics of the world from sensory signals. One challenge is that the environmental properties an organism must recognize are usually not explicit in the sensory input. A primary function of sensory systems is to transform raw sensory signals into representations in which behaviorally important features are more easily recovered. In doing so the brain must generalize across irrelevant stimulus variation while maintaining selectivity to the variables that matter for behavior. The nature of sensory codes and the mechanisms by which they achieve appropriate selectivity and invariance are thus a primary target of sensory system research.
The auditory system is believed to instantiate such representations through a sequence of processing stages extending from the cochlea into the auditory cortex. Existing functional evidence suggests that neurons in progressively higher stages of the auditory pathway respond to increasingly complex and abstract properties of sound [Chechik and Nelken, 2012, Atencio et al., 2012, Carruthers et al., 2015, Bizley et al., 2005, Elie and Theunissen, 2015, Russ et al., 2008, Mesgarani et al., 2008, Obleser et al., 2010, Bendor and Wang, 2005, Overath* et al., 2015, Norman-Haignere et al., 2015. Yet our understanding of the underlying transformations remains limited, particularly when compared to the visual system.
Feature selectivity throughout the auditory system has traditionally been described using linear receptive fields [Aertsen and Johannesma, 1981, Theunissen et al., 2000, Depireux et al., 2001. The most common instantiation is the spectrotemporal receptive field (STRF), which typically characterizes neural activity with a one-dimensional linear projection of the sound spectrogram transformed with a nonlinearity [Sharpee et al., 2011]. As a neural data analysis technique, STRFs are widespread in auditory neuroscience and have generated insight in domains ranging from plasticity to speech coding (e.g. [Fritz et al., 2003, Woolley et al., 2005).
Despite their utility, it is clear that STRFs are at best an incomplete description of auditory codes, especially in the cortex [Machens et al., 2004, Sahani and Linden, 2003, Williamson et al., 2016. Experimental evidence suggests that auditory neural responses are strongly nonlinear. As a consequence, auditory receptive fields estimated with natural sounds differ substantially from estimates obtained with artificial stimuli [Theunissen et al., 2000]. STRF descriptions also fail to capture the dimensionality expansion of higher representational stages. In contrast to the brainstem, neurons in the auditory cortex seem to be sensitive to multiple stimulus features at the same time [Atencio et al., 2012,Kozlov and Gentner, 2016, Harper et al., 2016. The presence of strongly non-linear behavior and multiplexing necessitates signal models more sophisticated than one-dimensional, linear features of the spectrogram such as STRFs.
An additional challenge to characterizing mid-level features of sound is that humans lack strong intuitions about abstract auditory structure. In specific signal domains such as speech, progress has been made by cataloging phonemes and other frequently occurring structures, but it is not obvious how to generalize this approach to broader corpora of natural sounds.
An alternative approach to understanding sensory representations that is less reliant on domain-specific intuition is that of efficient coding [Barlow, 1961, Attneave, 1954. The efficient coding hypothesis holds that neural codes should exploit the statistical structure of natural signals, allowing such signals to be represented with a minimum of resources. Numerous studies have demonstrated that tuning properties of neurons in early stages of the visual and auditory systems are predicted by statistical models of natural images or sounds [Srinivasan et al., 1982, Olshausen and Field, 1997, Bell and Sejnowski, 1997, van Hateren and van der Schaaf, 1998, Lewicki, 2002, Smith and Lewicki, 2006,Carlson et al., 2012,Terashima and Okada, 2012,Młynarski, 2015,Mlynarski, 2014,Carlin and Elhilali, 2013,Lee et al., 2008,Hosoya and Hyvärinen, 2015. Although the early successes of this approach engendered optimism, applications have largely been limited to learning a single stage of representation, and extensions to multiple levels of sensory processing have proven difficult. The underlying challenge is that there are many possible forms of high-order statistical dependencies in signals, and the particular dependencies that occur in natural stimuli are typically not obvious. The formulation of models capable of capturing these dependencies requires careful analysis and design [Hyvärinen and Hoyer, 2000, Karklin and Lewicki, 2005, Karklin and Lewicki, 2009, Cadieu and Olshausen, 2008, and perhaps good fortune, and is additionally constrained by what is tractable to implement. In the auditory system in particular, it remains to be seen whether modeling statistical signal regularities can reveal the complex acoustic structures and invariances that are believed to be represented in higher stages of the auditory system.
The primary goal of the present work was to discover such high-order structure in natural sounds and generate hypotheses about not-yet-observed intermediate-level neural representations. To this end we developed a probabilistic generative model of natural sounds designed to learn a novel stimulus representation -a population code of naturally occurring combinations of basic spectrotemporal patterns, analogous to STRFs.
The representations learned by the model from corpora of natural sounds suggest grouping principles in the auditory system. Model units learned to pool sets of similar spectrotemporal features, presumably because they occur together in natural sounds. In addition, some model units encode opponency between different sets of features.
Although not yet described in the auditory system, such tuning patterns may be analogous to phenomena such as end-stopping or cross-orientation suppression in the visual system. The representations learned by our model also resemble some recently reported properties of auditory cortical neurons, providing further evidence that natural-scene statistics can predict neural representations in higher sensory areas.

Methods and Models
Overview of the hierarchical model To learn mid-level auditory representations, we constructed a hierarchical, statistical model of natural sounds. The model structure is depicted in Fig 1. The model consisted of a stimulus layer and two latent layers that were adapted to efficiently represent a corpus of audio signals. Because our goal was to learn mid-level auditory codes, we did not model the raw sound waveform. Instead, we assumed an initial stage of frequency analysis, modeled after that of the mammalian cochlea. This frequency analysis results in a spectrogram-like input representation of sound, which we term a 'cochleagram,' that provides a coarse model of the auditory nerve input to the brain (Fig 1 A,  The first layer of the model (Fig 1, middle row) was intended to learn basic acoustic features, analogous to STRFs of early auditory neurons. Through the remainder of the paper, we refer to the features learned by the first layer as spectrotemporal kernels (STKs), in order to differentiate them from neurally derived STRFs. The second layer, depicted in the top row of Fig 1 A, was intended to learn patterns of STK co-activations that frequently occur in natural sounds.
The model specifies a probability distribution over the space of natural sounds, and its parameters can be understood as random variables whose dependency structure is depicted in Fig 1 B. The spectrogram x is represented with a set of spectrotemporal features φ convolved with the latent activation time-courses s. The second latent layer encodes the magnitudes of s with the basis functions B convolved with their activation time-courses v.
In the following sections, we present the details of each layer. First layer of model -convolutional, non-negative sparse coding of spectrograms The first layer of the model was designed to learn basic spectrotemporal features. One previous attempt to learn sparse spectrotemporal representations of natural sounds [Carlson et al., 2012] produced structures reminiscent of receptive fields of neurons in the auditory midbrain, thalamus, and cortex. Here, we extended this approach by learning a spectrogram representation which is sparse and convolutional. Figure 2: Explanation of first layer of the model. A cochleagram x of arbitrary length is represented as a sum of spectrotemporal kernels φ convolved with time-courses of corresponding coefficients s. Kernels span the full extent of the frequency dimension, and convolution occurred only in time. Coefficients s are non-negative and have sparse distributions (i.e. remain close to zero most of the time).
A schematic of the first layer is depicted in Fig 2. We modeled the cochleagram (x t,f ) as a linear combination of spectrotemporal kernels φ convolved in time with their activation time-courses s and distorted by additive Gaussian noise ξ t,f with variance σ 2 .: Kernel activations s i,t are assumed to be independent, i.e. their joint distribution is equal to the product of marginals: We assumed that each spectrotemporal kernel remains inactive for most of the time, i.e. that the distribution of its activations is sparse. Moreover, we imposed a non-negativity constraint on the coefficients s. This facilitates interpretations in terms of neural activity and improves the interpretability of the learned representation. These constraints were embodied in an exponential prior on the coefficients s: where λ i is the scale parameter. The scale of the exponential distribution determines its "spread" along the real line. Small values of λ i yield distributions tightly peaked at 0.
Large λ i generate broader distributions allowing s i,t to attain large values with higher probability. When training the first layer we assumed all λ i to be constant and equal to 1. As we will describe in the following sections, the second layer of the model relaxes that assumption and learns a representation of time-varying scale parameters λ i . The first layer of the model specifies the following negative log-posterior probability of the data: This negative log-probability can be viewed as a cost function to be minimized when inferring the value of coefficients s: while maintaining low-reconstruction error (first term on the right-hand side), the sparsity of representation should be maximized (the second term).
Dependencies between coefficients -a signature of mid-level structure Although the sparse coding strategy outlined above learns features that are approximately independent across the training set, residual dependencies nonetheless remain. In part this is because not all dependencies can be modeled with a single layer of convolutional sparse coding. However, due to the non-stationary nature of natural audio, coefficient dependencies can also be present locally despite not being evident across a large corpus as a whole. Empirically, the learned features exhibit dependencies specific to particular sounds [Karklin and Lewicki, 2009], and thus deviate locally from their (approximately independent) marginal distribution. For example, a spoken vowel with a fluctuating pitch contour would require many harmonic STKs to become activated, and their activations would become strongly correlated on a local time scale. Such local correlations reflect higher-order structure of particular natural sounds. Statistically speaking, this is an example of marginally independent random variables exhibiting conditional dependence (in this case conditioned on a particular point in time or a type of sound).
Fig 3 depicts such dependencies via cross-correlation functions of selected STK activations for a white noise sample and for three different natural sounds. Coefficient correlations (black curves in each subplot) vary from sound to sound, but in all cases deviate from those obtained with noise (gray bars within subplots), revealing dependence. These dependencies are indicative of "mid-level" auditory features, perhaps analogous to the correlations between oriented Gabor filters induced by an elongated edge. In this work, we exploited the fact that intermediate level representations can be learned by modeling dependencies among first-layer features [Karklin and Lewicki, 2009, Karklin and Lewicki, 2005, Cadieu and Olshausen, 2008.

Variability of STK activations
A second phenomenon evident in STK activations is that particular patterns of co-activation occur with some variability. This is visible in Fig 4, which depicts activations of selected STKs when encoding multiple exemplars of the same sound -the word "one" spoken twice by the same speaker ( Fig 4A) and two exemplars of water being poured into a cup ( Fig  4B). Figure 3: Dependencies among spectrotemporal feature activations. A) When encoding white noise (cochleagram depicted on top), sparse coefficients are uncorrelated on short time-scales. This is visible in the shape of the coefficient cross-correlation functions (black lines) corresponding to the four most strongly activated STKs (bottom, right column). Due to the lack of temporal structure, the cross-correlation between STKs is flat, and the autocorrelation of individual STKs is a Dirac delta function at 0. Here and in other panels, the mean was substracted from the (non-negative) coefficient trajectories prior to computing the cross-correlation. B) When encoding structured stimuli such as an owl vocalization, STK activations reveal strong local correlations -the cross-correlations deviate from those for white noise (thick gray lines). C) Same as B, for a hammer hit. D) Same as B, for a speech excerpt.  It is apparent that the STK coefficient trajectories for the two exemplars in each case (black and gray lines) reflect the same global pattern even though they differ somewhat from exemplar to exemplar. The similarity suggests that the trajectories could be modeled as different samples from a single time-varying distribution parameterized by a nonstationary coefficient magnitude. When the magnitude increases, the probability of a strong STK activation increases. Retaining the (inferred) time-varying magnitude instead of precise values of STK coefficients would yield a representation more invariant to lowlevel signal variation, potentially enabling the representation of abstract regularities in the data. Such a representation bears an abstract similarity to the magnitude operation used to compute a spectrogram, in which each frequency channel retains the time-varying energy in different parts of the spectrum. Here we are instead estimating a scale parameter of the underlying distribution, but the process similarly discards aspects of the fine detail of the signal.

Second layer of model -encoding of STK combinations
The second layer of the model was intended to exploit the two statistical phenomena detailed in the previous section: conditional dependencies between STKs and their variation across exemplars. Similarly to the first layer, the second layer representation is formed by a population of sparsely activated basis functions. These basis functions capture local dependencies among STK magnitudes by encoding the joint distribution of STK activations rather than exact values of STK coefficients. The resulting representation is thus more specific than the first-layer code -instead of encoding single features independently, it signals the presence of particular STK combinations. It is also more invariant, generalizing over specific coefficient values.
Because the proposed representation is a population code of a distribution parameter, it bears conceptual similarity to previously proposed hierarchical models of natural stimuli that encoded patterns of variance [Karklin andLewicki, 2005, Bumbacher andMing, 2012], covariance [Karklin and Lewicki, 2009] or complex amplitude [Cadieu andOlshausen, 2008, Hyvärinen andHoyer, 2000]. The novelty of our model structure lies in being convolutional (i.e., it can encode stimuli of arbitrary length using the same representation) and in parameterizing distributions of non-negative STK coefficients, increasing the inter-pretability of the learned spectrogram features. The novelty of the model's application is to learn hierarchical representations of sound (previous such efforts have largely been restricted to modeling images; though see [Lee et al., 2009]). Figure 5: Explanation of second layer of the model. A) An array of STK activations s (a "STK-gram") serves as an input to the second layer. Rows correspond to first layer features φ i and columns to time points. B) The second layer uses a population of features B to encode the logarithm of STK activation magnitudes. C) Coefficient trajectories s (thin grey lines) and their magnitudes λ (thick red and black lines) for two example STKs. D) Distributions of s at time points t 1 and t 2 are depicted in the right column.
The second layer of the model is depicted schematically in Fig 5A. We assume that STK activations s are samples from a non-stationary exponential distribution with timevarying scale parameter λ i,t , relaxing the assumption of stationary λ i made in learning the first layer: When λ i,t is high, the distribution of s i,t becomes broader ( Fig 5C, black line in the top row, red line in the bottom). This allows the coefficient s i to attain large values. For small values of the scale parameter, the probability density is concentrated close to 0 ( Fig  5D), red line in the top row, black line in the bottom), and coefficients s i,t become small. To model the magnitudes λ (which are non-negative, akin to variances), we took their logarithm, mapping their values onto the entire real line so that they could be represented by a sum of real-valued basis functions.
Patterns of STK magnitudes are represented in the second layer by a population of features B convolved with coefficients v: where ρ is a bias vector. Each second-layer basis function B represents a particular temporal pattern of co-activation of first-layer STKs. Their corresponding coefficients v are assumed to be sparse and independent: where α controls the degree of sparsity.
As with the first layer, learning and inference are performed by gradient descent on the negative log posterior. Because the second layer units encode combinations of sparse first-layer coefficients, we placed a sparse prior on the L 1 norm of the basis functions B. The overall cost function to be minimized during learning in the second layer is then: where β controls the strength of the sparse prior on B, and T b is the temporal extent of each second-layer basis function. As in the first layer cost function (Eq 5), the first term on the right hand side of Eq 10 enforces a match of the representation to data (the magnitudes λ are pushed away from zero towards the observed coefficients s), while the second and third terms promote sparsity of second-layer coefficients and basis functions, respectively. Sparsity of second layer units is a reasonable assumption given the sparsity of first-layer coefficients, and we found that using a sparse prior on second-layer units (i.e. setting β to be larger than 0) was necessary to achieve convergence. Others have found that the addition of sparse priors can enable faster learning with less data without substantially altering the obtained solution [Hyvärinen and Raju, 2002].

Learning procedure
The two layers of the model were trained separately, i.e. the training of the second layer occurred after the first layer training was completed. In the first layer, training was performed with an EM-like procedure that iteratively alternated between inferring STK coefficients and updating STK features [Olshausen andField, 1997, Cadieu andOlshausen, 2008]. Spectrotemporal features φ were initialized with Gaussian white noise. For each excerpt in the training set, coefficients s i,t were inferred via gradient descent on the energy function (5). Because inference of all coefficients s is computationally expensive, we adopted an approximate inference scheme [Blumensath and Davies, 2006]. Instead of inferring values of all coefficients for each excerpt, we selected only a subset of them to be minimized. This was done by computing the cross-correlation between a sound excerpt and features φ i and selecting a fixed number of the largest coefficients s i,t . The inference step adjusted only this subset of coefficients while setting the rest to 0. Given the inferred coefficients, a gradient step on the spectrotemporal features φ was performed.
Each learning iteration therefore consisted of the following steps: 1. Draw a random sound excerpt from the training data set. Excerpts were 403 ms in length (129 time samples of the spectrogram, sampled at 320 Hz).
2. Compute the cross-correlation of all basis functions φ i with the sound excerpt. Select the 1024 pairs of coefficient indices and time-points (i, t) that yield the highest correlation values.
3. Infer the values of the selected coefficients by minimizing Eq. 5 with respect to s i,t via gradient descent. Set the rest of coefficients to 0.
4. Compute the gradient step on the basis functions as the derivative of Eq. 5 with respect to basis functions φ using inferred coefficient valuesŝ. Update basis functions according to the gradient step.
5. Normalize all basis functions to unit norm.
This procedure was terminated after 200000 iterations. The second layer was then learned via the same procedure used for the first layer. In each iteration a 528 ms long (169 samples at 320 Hz) randomly drawn sound excerpt was encoded by the first layer, and the resulting matrix of coefficients s served as an input to the second layer. A subset of coefficients v was selected for approximate inference by computing the cross correlation between features B i and the logarithm of the first-layer coefficients s (analogous to step 2 in the procedure described above for the first layer). The energy function E 2 was first minimized with respect to coefficients v followed by a gradient update to the basis functions B (analogous to steps 3 and 4 for the first layer). Entries in the bias vector ρ corresponding to each coefficient s i were set to the expectation of the coefficient across the entire training set: , the estimate of the marginal scale parameter λ i for the corresponding STK). Learning was again terminated after 200000 iterations. To validate the learning algorithm, we ran it on a toy data set generated using known second-layer units. The algorithm successfully recovered the parameters of the generating model, as desired (see Appendix C for these results).

Training data and spectrogram parameters
We trained the model on two different sound corpora. The first corpus was the TIMIT speech database [Garofolo et al., 1993]. The second corpus combined a set of environmental sounds (the Pitt sound database [Lewicki, 2002]) and a number of animal vocalizations downloaded from freesound.org. The environmental sounds included both transient (breaking twigs, steps, etc.) and ambient (flowing water, wind, etc.) sounds; the animal vocalizations were mostly harmonic.
We computed cochleagrams by filtering sounds with a set of 65 bandpass filters intended to mimic cochlear frequency analysis. Filters were equally spaced on an equivalent rectangular bandwidth (ERB) scale [Glasberg and Moore, 1990]), with parameters similar to that from a previous publication [McDermott and Simoncelli, 2011]. Center frequencies ranged from 200 Hz to 8 kHz. We computed the Hilbert envelope of the output of each filter and raised it to the power of 0.3, emulating cochlear amplitude compression [Robles and Ruggero, 2001]. To reduce dimensionality, each envelope was downsampled to 320 Hz.
We set the number of features in the first layer to 128 and in the second layer to 100. Pilot experiments yielded qualitatively similar results for alternative feature dimensionalities. Each first layer feature encoded a 203 ms interval (65 time samples of the spectrogram).

First layer: Basic spectrotemporal features of natural sounds
The first-layer features learned from each of the two sound corpora are shown in Fig 6  A and B. These features could be considered as the model analogues of neural STRFs.
The vast majority of features are well localized within the time-frequency plane, encoding relatively brief acoustic events. The STKs learned from speech included single harmonics and harmonic "stacks" (Fig 6 A -features numbered 1 and 2), frequency sweeps (feature 3), and broadband clicks (feature 4). The features learned from environmental sounds also included single harmonics and clicks (Fig 6 B -features 1, 2 and 4). In contrast to the results obtained with speech, however, harmonic stacks were absent, and a number of high-frequency hisses and noise-like features were present instead (feature 3).   Figure 6: Spectrotemporal kernels learned by the first layer. A) A population of STKs learned from a speech database. Representative STKs for each corpus are magnified and numbered from 1-4 for ease of reference. B) Population of STKs learned from environmental sounds. C) Speech-trained STK population plotted on time-frequency (left) and spectral-temporal modulation (right) planes. Each dot corresponds to a single STK, and its color encodes the log of its mean coefficient (averaged over entire training dataset). D) The STK population trained on environmental sounds, represented as in (C).

B) A)
The properties of the learned dictionaries are also reflected in distributions of feature locations in the time-frequency and spectrotemporal modulation planes, as visible in Fig  6C and D. Each dot position denotes the center of mass of a single STK, while its color signals the feature's average coefficient value over the stimulus set. For both speech and environmental sounds, the learned STKs uniformly the audio frequency spectrum (Fig 6 C and D, left panels). Due to the convolutional nature of the code, the energy of each feature is concentrated near the middle of the time axis. The modulation spectra of the learned STKs (Fig 6 C and D, right panels) are somewhat specific to the sound-corpus. Features trained on a speech corpus were more strongly modulated in frequency, while environmental sounds yielded STKs with faster temporal modulations. STKs learned from both datasets exhibit a spectrotemporal modulation tradeoff: if a STK is strongly temporally modulated its spectral modulation tends to be weaker. This tradeoff is an inevitable consequence of time-frequency conjugacy [Singh and Theunissen, 2003], and is also found in the STRFs of the mammalian and avian auditory systems [Miller et al., 2002, Woolley et al., 2005.

Second layer: Combinations of spectrotemporal features
STKs captured by the first layer of the model reflect elementary features of natural sounds. By contrast, the features learned by the second layer capture how activations of different STKs cluster together in natural sounds, and thus reflect more complex acoustic regularities. We first present several ways of visualizing the multi-dimensional nature of the second-layer representation, then make some connections to existing neurophysiological data, and then derive some neurophysiological predictions from the model.

Visualizing second-layer features
The second-layer features encode temporal combinations of STK log-magnitudes. Each feature can be represented as a N × T b dimensional matrix, where rows correspond to STKs and columns to time-points. An example second-layer unit is shown in Fig 7A  (left panel). A positive value in the i-th row and t-th column of a feature B j encodes a local increase in the magnitude of the i-th STK. A negative value encodes a decrease in magnitude.
An alternative visualization is to examine the spectrotemporal structure of the STKs that have large weights in the second-layer feature. The same feature B j is depicted in this way in the center panel of Fig 7A, which displays the four STKs with highest average absolute weights for this particular feature. To the right of each STK are its weights (i.e., the corresponding row of the B j matrix). It is apparent that the weights increase and decrease in a coordinated fashion, and thus likely encode particular dependencies between the STKs.
To summarize the full distribution of STKs contributing to a second-layer unit, we adopted the visualization scheme illustrated in the right panel of Fig 7A. We plot the center of mass of each STK in the modulation (top row) and time-frequency (middle row) planes, as in Fig 6C and D. The dot for a first-layer STK is colored red or blue, depending on the sign of their time-averaged weight, with the average absolute value of the weight signaled by the intensity of the color. The bottom row of the panel depicts the temporal pattern of STK magnitudes -line colors correspond to dots in the top and middle rows of the panel. Although the weights of most STKs maintain the same sign over the temporal support of the second-layer unit, they were not directly constrained in this regard, and in some cases the weight trajectories cross zero. Representative examples of second-layer basis functions are depicted in this way in Fig  7B. We separated them into two broad classes -"excitatory" units (columns 1 and 2 in Fig 7), which pool STKs using weights of the same sign, and therefore encode a pattern of coordinated increase in their magnitudes, and "excitatory-inhibitory" units (columns 3 and 4 in Fig 7) which pool some STKs with positive average weights and others with negative average weights. We note that excitatory-only and inhibitory-only units are functionally interchangeable in the model, because the encoding is unaffected if the sign of both the STK weights and coefficients are reversed. From inspection of Fig 7B it is evident that some second-layer units pool only a few STKs (e.g. A1, C1, D1) while others are more global and influence activations of many first-layer units (e.g. A2, A3, C3).

Second layer features encode patterns of STK dependencies
To better understand the structure captured by the second-layer units, we considered their relationship to the two kinds of dependencies in STK activations that initially motivated the model. In Methods and Models section we observed that STK activations to natural sounds exhibit strong local cross-correlations as well as variations of particular temporal activation patterns. The second layer of the model was formulated in order to capture and encode these redundancies.   To first test whether the second layer captures the sorts of residual correlations evident in the first layer output, we measured cross-correlations between STK activations conditioned on the activation of particular second-layer units. Fig. 8A-D depicts four example second-layer units (top row, left column) along with an example stimulus that produced a strong positive response in the unit (selected from the TIMIT corpus). The bottom section of each panel depicts cross-correlation functions of activations of four STKs, averaged over 25 stimulus epochs that produced a strong positive response of the second-layer unit. The cross-correlation functions deviate substantially from 0, as they do when conditioned on excepts of natural sounds (Fig. 3). These correlations reflect the temporal pattern of STK coefficients that the second-layer unit responds to.   Figure 9: Second-layer units generalize across variation in STK coefficients.Panels A-D correspond to second-layer units depicted in Fig. 8 Each panel shows four STKs pooled with strongest weights (left column) by the corresponding second-layer unit. Next to each STK are their activation patterns (right column) for the 5 stimuli eliciting the strongest response in the second-layer unit. Despite some variability they share a global trend.
We next examined whether the second-layer units respond to STK activations fluctuating around particular global patterns (as depicted in Fig.4). Because the second layer of the model represents the magnitude, rather than the precise values of first-layer coefficients, it should be capable of generalizing over minor STK coefficient variation. Fig.9 plots STK coefficient trajectories for stimuli eliciting a strong response in the second-layer features from Fig. 8. The STK activation traces reveal variability in each case, but nonetheless exhibit a degree of global consistency, as we saw earlier for natural sound exemplars (Fig.4). These results provide evidence that the second layer is capturing the dependencies it was intended to model.

Experimental predictions
Pooling of similiar spectrotemporal patterns What do the model results suggest about mid-level auditory structure? Our visualizations of second layer features in Fig 7B revealed that many represent the concurrent activation of many STKs. Examination of Fig 7B (columns 1 and 2) suggests that the first layer STKs that such units pool are typically similar either in their spectrotemporal or modulation properties. For example, the unit depicted in panel A2 encodes joint increases in the magnitude of high-frequency STKs, while that in panel C2 encodes joint increases in lowfrequency STKs.
In order to quantitatively substantiate these observations, we measured the spread of the center-of-mass of each of the STKs pooled by "excitatory-only" second-layer units (those whose STK weights exceeding 5% of the maximum absolute weight were all of the same sign). Specifically, we computed the standard deviation of the center-of-mass on each dimension of the time-frequency and modulation planes, for all STKs with weights exceeding 5% of the maximum weight for the unit, and compared this to the standard deviation of centers of mass for random samples of STKs (matched in size to the mean number of STKs pooled by the excitatory-only units). These standard deviations measure the spread of STKs in spectrotemporal and modulation domains; if the pooled STKs are similar, the spread should be small.
As shown in Fig.10, the STKs pooled with the same sign were more similar in each of the four analyzed dimensions than random groups of STKs, for both sound corpora (red bar vs. gray bar, t-tests yielded p<.05 in all cases). This result suggests the hypothesis that STRFs that are pooled by downstream auditory cortical neurons with excitatory weights [Atencio et al., 2012, Kozlov andGentner, 2016] should also tend to be more similar than expected by chance. Fig 7B also shows examples of a distinct set of second-layer units in which increased activation of one group of STKs (again, typically similar to each other) is associated with a decrease in activity of another group of STKs. For example, the unit in panel B3 encodes an increase in magnitude of temporally modulated features (clicks) together with a simultaneous decrease of activation of a strongly spectrally modulated (harmonic) STK. Such "opponency" was evident in a large subset of second layer units, and represents the main novel phenomenon evident in our model.

Opponency patterns in mid-level audition
To our knowledge no such opponent tuning has been identified in auditory neuroscience, but qualitatively similar opponent behavior is evident in visual neurons exhibiting endstopping or cross-orientation inhibition. The results raise the possibility that coordinated excitation and inhibition could be a feature of central auditory processing, and we thus examined this model property in detail.
To quantitatively substantiate the observation of opponency, we again examined the variation in STKs pooled by the second-layer units, in this case those that pooled STKs with both positive and negative weights. We analyzed all units that had at least one STK weight of each sign whose absolute value exceeded 5% of the maximum weight. For each such second-layer unit we identified centers of mass of the STKs pooled with weights exceeding 5% of the maximum and again computed their standard deviation along timefrequency and modulation coordinates. We compared the standard deviation for STKs pooled with the same sign (separately computed for positive and negative signs) to that of STKs pooled with both signs. If opponent units contrast different types of features we expect the standard deviation of STKs pooled with both signs to be higher than that of STKs pooled with a single sign.
As shown in Fig.10, even for STKs pooled with the same sign, the spread was typically higher for opponent units than for excitatory-only units, for both speech and environmental sounds (blue vs. red bars in Fig.10A and B respectively). However, the average spread of STKs pooled with the both signs by opponent units (Fig.10, dark blue bars) was significantly higher than the spread of STKs pooled with the same sign (Fig.10, light blue bars). This finding held for all dimensions and both corpora (p < .05 in all cases, via t-test). ). Moreover, the spread of STKs pooled with the same sign by opponent units were also significantly lower than the spread of random groups of STKs (light blue bars vs. gray bars; p < .05 in all cases, via t-test). These results provide quantitative evidence that second-layer units pool features in a structured way and that opponent units assign opposite signs to to groups of STKs that are similar within groups but differ across groups. max std.dev.
[ms]  One might imagine that examination of the STKs pooled by each unit would be sufficient to determine the sort of stimuli eliciting strong positive or negative responses. But because the activation of each unit is the result of a non-linear inference process in which units compete to explain the stimulus pattern [Olshausen and Field, 1997], it is often not obvious what a set of STKs will capture. Thus to understand which stimuli 'excite' or 'inhibit' second-level features, we inferred coefficients v by encoding the entire training dataset and selected two sets of 25 sound epochs that elicited the strongest positive and strongest negative responses, respectively, in each unit. Examples of opponent stimulus patterns encoded by second-layer features are visualized in Fig 11. These positive and negative stimuli are depicted in the second and third columns from the left, respectively.
In the last two columns, the center of mass of each of these stimuli is plotted on the timefrequency plane (fourth column) and modulation plane (fifth column). Although the center of mass of each stimulus is admittedly a crude summary, the simplicity of the representation facilitates visualization and analysis of the clustering of positive and negative stimuli. We note also that because second-layer coefficients and weights can be sign-reversed without changing the representation, the designation of stimuli (and STK weights) as 'excitatory' or 'inhibitory' is arbitrary.
In some cases, some natural function can be ascribed to the unit, particularly upon listening to stimuli eliciting positive and negative responses. For instance, positive activations of the unit shown in Fig 11B appear to encode onsets of voiced speech, while its negative activations encode voicing offsets. By contrast, the unit shown in Fig 11H appears to code speaker gender. For this unit, representations of positive and negative stimuli were mixed on the time-frequency and modulation planes, but the excitatory stimuli exhibited somewhat coarser spectral modulation (red circles in the right column are more concentrated in the lower part of the plane than blue circles). Listening to the stimuli revealed a clear difference in pitch/gender (23 out of 25 of the negative stimuli were vowels uttered by a female speaker, while all of the positive stimuli were vowels uttered by a male speaker, as confirmed by the first author upon listening to the stimuli). This observation underscores the fact that the time-frequency and modulation planes are but two representational spaces in which to assess stimuli, and are used here primarily because they are standard. Moreover, the centroids in these planes are but one way to summarize the STK. Although opponent sets of STKs and positive/negative stimuli often exhibited separation in one of the two planes when analyzed in terms of centroids, they need not do so in order to be distinct.
Taken together, our results demonstrate opponent units that encode STK activation patterns which are mutually exclusive and presumably do not co-occur in natural signals. The phenomenon is a natural one to investigate in the auditory cortex.
[cyc/oct] Figure 11: Opponent stimulus selectivity in second layer units. Each row corresponds to a particular second-layer unit. The leftmost column plots the center of mass of each first-layer STK in the modulation and time-frequency planes, along with the time courses of their weights in the second-layer unit (as in Fig 7B). The second and third columns from the left depict 25 stimuli eliciting strong positive and negative responses in the corresponding second-layer unit. In the fourth and fifth columns, positive and negative stimuli are visualized as red and blue circles, respectively, in time-frequency and modulation planes (the circle is located at the center of mass of the stimulus). Large circles correspond to centroids of positive and negative stimulus clusters.

Comparison with neurophysiological data
Although our primary goal was to generate predictions of not-yet observed neural representations of sound, we also sought to test whether our model would reproduce known findings from auditory neuroscience. The first layer STKs replicated some fairly standard findings in the STRF literature, as discussed earlier. To compare the results from the second layer to experimental data, we examined their receptive field structure and the specificity of their responses, and compared each to published neurophysiology data. Figure 12: Comparison with neurophysiological data. A) Spectrotemporal receptive fields estimated for second layer units trained on speech and environmental sounds. B) Comparison of spectral and temporal modulation tuning of first and second layer receptive fields (panels in left column) to experimental measurements in auditory thalamus and cortex of the cat (panels in right column were provided by Lee Miller [Miller et al., 2002]). Higher processing stages both in the model and the auditory system exhibit tuning for coarser modulations in both frequency and time. Sub-panels with experimental data reprinted with permission of original author. C) Example activation trajectories of first (middle row) and second layer (bottom row) units to an excerpt of speech (top row). The activation of the first layer unit is tightly locked to presence of a preferred stimulus. Responses of the second layer unit are less specifically locked to particular spectrotemporal structures. D) Comparison of tuning specificity, as measured with the Feature Selectivity Index (FSI), for first and second layers of the model trained either on speech or environmental sounds (gray and red bars, respectively), and for different layers of the auditory cortex of the cat (right panel replotted from [Atencio et al., 2009]).

Differences in modulation tuning between first and second model layers
First, we estimated spectrotemporal receptive fields for units in both layers of the model. The receptive fields of the first-layer units are simply the STK of the unit. To estimate receptive fields of a second-layer unit, we drew inspiration from the spike-triggered average, generating a number of cochleagram samples from each basis function and averaging them. To generate samples from the j−th basis function, we set a coefficient v j to 1 with all other v i =j set to zero. We then sampled STK activation trajectories from the distribution dictated by the second-layer feature's coefficient and weights, convolved them with the corresponding STKs and summed the results. We then averaged multiple such samples together. Although we could have computed something more directly analogous to a spike-triggered average, the average sample (which we can compute only because we have the underlying generative model, unlike when conducting a neurophysiology experiment) has the advantage of alleviating the influence of stimulus correlations on the signature of the receptive field. See Appendix B for an illustration of this analysis. "Receptive fields" obtained in this way are depicted in Fig 12A. To compare the model units to neurophysiological data, we generated histograms of average spectral and temporal modulation frequency (center of mass in the modulation plane) of first-and second-layer receptive fields and plot them next to distributions of preferred modulation frequencies of neurons in the auditory thalamus and cortex of the cat [Miller et al., 2002] (Fig 12B). The same trend is evident in the model and the auditory system: the second-layer prefers features with slower/coarser spectral and temporal modulations relative to the first-layer, mirroring the difference seen between the cortex and thalamus. This analysis used features trained on speech, but environmental sounds yielded qualitatively similar results. In the model as well as the brain, lower modulation frequencies may result in part from combining multiple distinct STKs in downstream units.

Model units reflect neural hierarchical trends of response specificity
Neuronal tuning in early and late stages of the auditory system also tends to differ in specificity [Chechik andNelken, 2012, Kozlov andGentner, 2016]. Compared to the auditory brainstem, cortical neurons are less selective and respond to multiple features of sound [Kozlov andGentner, 2016, Atencio et al., 2012], consistent with an increase in abstraction of the representation [Chechik and Nelken, 2012]. Suggestions of similar behavior in our model are apparent in the activations of first-and second-layer units to sound, an example of which is shown in Fig 12C. The first layer feature (middle row) becomes activated only when it is strongly correlated with the stimulus. In contrast, activations of a typical second layer feature (bottom row) deviate from zero during many, seemingly different parts of the stimulus.
We quantified the specificity of tuning with the feature selectivity index (FSI), a measure introduced previously to quantify how correlated a stimulus has to be with a neuron's spike-triggered average to evoke a response [Miller et al., 2001] (see Appendix B for definition). An FSI equal to 1 implies that a neuron spikes only when a stimulus is precisely aligned with its STRF (defined as the spike-triggered average), whereas an FSI equal to 0 means that the stimuli triggering neural firing are uncorrelated with the STRF. We computed the FSI using the 25 cochleagram excerpts that most strongly activated each of the first-and second-layer units. The average FSI of each model layer is plotted in Fig 12D. Second-layer features are substantially less specific than first-layer features. A similar effect occurs across different cortical layers [Atencio et al., 2009] (Fig 12D, right). Analogous differences seem likely to occur between thalamus and cortex as well, although we are not aware of an explicit prior comparison. The decrease in response specificity in our model can be explained by the fact that second layer units can become activated when any of the pooled first-layer features (or their combination) appears in the stimulus.

Discussion
Natural sounds are highly structured. The details of this structure and the mechanisms by which it is encoded by the nervous system remain poorly understood. Progress on both fronts is arguably limited by the shortage of signal models capable of explicitly representing natural acoustic structure. We have proposed a novel statistical model that captures an unexplored type of high-order dependency in natural sounds -correlations between the activations of basic spectrotemporal features. Our model consists of two layers. The first layer learns a set of elementary spectrotemporal kernels and uses them to encode sound cochleagrams. The second layer learns a representation of co-occurence patterns of the first layer features.
We adopted a generative modelling approach, inspired by its previous successes in the domain of natural image statistics. Previous hierarchical, probabilistic models of natural images were able to learn high-order statistical regularities in natural images [Karklin and Lewicki, 2009, Karklin and Lewicki, 2005, Cadieu and Olshausen, 2008, Hyvärinen et al., 2009, Hoyer and Hyvärinen, 2002, Lee et al., 2008, Hosoya and Hyvärinen, 2015, Garrigues and Olshausen, 2008, Berkes et al., 2009. In many cases the representations learned by these models exhibit similarity to empirically observed neural codes in the visual system.
The hierarchical representations learned by our model provide predictions about neural representations of mid-level sound structure. Moreover, the model reproduces certain aspects of the representational transformations found through the thalamus and cortex. The results suggest that principles of efficient coding could shape mid-level processing in the auditory cortex in addition to the auditory periphery [Lewicki, 2002].

Statistical dependencies in natural sounds
Our model exploits statistical dependencies between first-layer spectrotemporal kernels. These dependencies have two likely causes. First, dependencies almost surely remain from limitations of the convolutional sparse coding model, in that first-layer coefficients are never completely marginally independent even after learning (due to insufficient expressive power of the code). Second, even if the first-layer coefficients were fully marginally independent they could exhibit local dependencies when conditioned on particular sound excerpts. Observations of such conditional dependence have been made previously, mostly in the context of modelling natural image statistics [Karklin and Lewicki, 2009,Karklin and Lewicki, 2005,Hyvärinen et al., 2009,Cadieu and Olshausen, 2008,Karklin et al., 2012. We similarly observed strong non-zero cross-correlations between STKs for particular natural sound excerpts, and found that the second-layer units captured these sorts of dependencies (Fig. 3,8).
We also observed that different instances of the same acoustic event (such as the same word uttered twice by the same speaker) yield similiar STK coefficient trajectories (Fig. 4), which can be thought of as samples from a nonstationary distribution with a time-varying magnitude parameter. By modeling this magnitude the model learned representations with some degree of invariance (Fig. 9).

Model results -Convolutional spectrotemporal features
The model first learned a sparse convolutional decomposition of cochleagrams. Convolutional representations are less redundant than "patch-based" codes and can represent signals of arbitrary length, and have been succesfully applied to spectrogram analysis in engineering (e.g. [Blumensath andDavies, 2006,Grosse et al., 2012]). The features learned with convolutional sparse coding improved sound classification accuracy, musical note extraction and source separation in prior work [Blumensath andDavies, 2006, Grosse et al., 2012]. By contrast, the model presented here is, to our knowledge, the first use of convolutional sparse coding of audio in a neuroscience context. Accordingly, we provide an in-depth analysis of feature shapes learned by the first layer and relate them to known spectrotemporal tuning properties of auditory neurons.
The learned spectrotemporal kernels span clicks, harmonics, combinations of harmonics, bandpass noise, frequency sweeps, onsets, and offsets. Some of these structures are present in previous learned codes of spectrograms, but due to the convolutional nature of the learned code, the model learns only a single version of each feature rather than replicating it at different time points. Although difficult to quantify, our impression is that the code is more diverse than that obtained with previous patch-based approaches [Carlson et al., 2012. As in previous neurophysiological measurements of cortical STRFs, the model STKs tile the spectrotemporal modulation plane, subject to the contraints of the tradeoff between time and frequency [Singh and Theunissen, 2003].
As with features learned from sound waveforms [Lewicki, 2002, Smith andLewicki, 2006], we found the first-layer spectrogram features to depend on the training corpus (speech or a set of environmental sounds). Although certain spectrotemporal patterns (e.g. single harmonics or clicks) appeared in both features sets, others were corpus-specific. Corpus-specific structure was also evident in the distributions of features in the modulation plane. This corpus dependence contrasts with the relatively consistent occurrence of Gabor-like features in sparse codes of images. These results raise the possibility that auditory cortical neuronal tuning might exhibit considerable heterogeneity, given that the full range of natural audio must be encoded by cortical neurons.

Model results -Second-layer features
The second layer of the model encodes a probability distribution of spectrotemporal feature activations. By jointly modelling the magnitudes of the first layer responses, the second layer basis functions capture patterns of spectrotemporal kernel covariation. Magnitude modeling was essential to learning additional structure -we found in pilot experiments that simply applying a second layer of convolutional sparse coding did not produce comparable results. This observation is consistent with previous findings that nonlinear transformations of sparse codes often help to learn residual dependencies [Karklin and Lewicki, 2009, Shan et al., 2007, Hyvärinen et al., 2009.
The second layer of our model learned a representation of spectrotemporal feature co-activations that frequently occur in natural sounds. Typically, the spectrotemporal features pooled by a second-layer unit shared some property, often evident in at least one of the time-frequency or modulation planes -e.g. frequency content, temporal pattern or modulation characteristics (Fig. 7), quantified by the similarity in these planes. Our model also identified 'opponent' patterns -sets of spectrotemporal features that are pooled with opposite sign, presumably because they are rarely active simultaneously in natural sounds (Fig. 11). The opponent features for a second-layer unit were often at least partially separated in at least one of the spectrotemporal modulation or time-frequency planes (10).

Relation to auditory neuroscience
The structure learned by the model from natural sounds replicates some known properties of the auditory system. The modulation frequencies preferred by units dropped from the first to the second layer (Fig. 12B), as has been observed between the thalamus and cortex. Unit "tuning" specificity (i.e. the similarity between stimuli eliciting a strong response) also decreased from the first to the second layer (Fig. 12C,D). A similar specificity decrease has been observed between granular and supra-and infra-granular layers of the cortex [Atencio et al., 2009]. Although we do not suggest a detailed correspondence between layers of our model and particular anatomical structures, these similarities indicate that some of the principles underlying hierarchical organization of the auditory pathways may derive from the natural sound statistics and architectural choices that constrain our model.
The combinations of spectrotemporal features learned by our model provide hypotheses for neural tuning that might be present in the auditory system. Some units combined only STKs with similar acoustic properties (Fig. 7, 10). Others encode activations of 'opponent' sets of first-layer kernels. Such opponent kernel sets are presumably those that do not typically become active simultaneously in the training corpus. The results are somewhat analogous to excitation-inhibition phenomena in visual neurophysiology (such as end-stopping, length and width suppression or cross-orientation inhibition) emerging in models of natural images [Karklin and Lewicki, 2009, Coen-Cagli et al., 2012, Zhu and Rozell, 2013, Hosoya and Hyvärinen, 2015. The second-layer units of our model provide candidates of potentially analogous phenomena in the auditory cortex.
Our results are also relevant to recent evidence that central auditory neurons in some cases are driven by more than one stimulus feature [Atencio et al., 2009, Sharpee et al., 2011, Kozlov and Gentner, 2016. Our model reproduces that trend, pooling up to dozens of single layer features with strong weights, and predicts that distinct dimensions may be combined in opponent fashion.

Relation to other modeling approaches
The model described here represents one of several approaches to construction of hierarchical signal representations. The representations in our model are learned from natural sounds, and thus contrast with hand-engineered models that seek to replicate known or hypothesized features of sensory coding [Chi et al., 2005,McDermott andSimoncelli, 2011]. One advantage of learning models from natural signals is that any similarities with known neural phenomena provide candidate normative explanations for these phenomena, e.g. that they arise from the demands of efficient coding or some other optimality constraint imposed by the model [Olshausen andField, 1997, Lewicki, 2002]. Another potential advantage is that one might hope to learn structures that have not yet been observed neurally, but that could provide hypotheses for future experiments. The latter was the main motivation for our modeling approach.
The model is also unsupervised, and generative, specifying a joint probability distribution over the data and coefficients in both latent layers. An alternative approach to learning hierarchical representations of data is to use discriminative models, which optimize a representation for performance of a particular task. Prominent recent examples of such discriminative learning come from the field of deep neural networks [Yamins et al., 2014]. Models of this class are comparatively straightforward to optimize and can yield high performance on classification tasks, but typically require large numbers of labeled training examples. The requirement of labeled training data is sometimes a practical limitation, and raises questions about the extent to which the learning procedure could relate to the brain. Generative models currently have the disadvantage of requiring custom optimization procedures. However, their generative nature facilitates certain applications (e.g. denoising) and allows samples to be straightforwardly generated, potentially for use in experiments. Moreover, they do not require labeled data, and learn representations that are independent of any particular task, much like the unsupervised learning believed to occur in sensory systems.

Conclusion
We have presented a hierarchical model of natural sounds. When trained on natural sound corpora, the model learned a representation of spectrotemporal feature combinations. The properties of the model layers resemble aspects of hierarchical transformations previously observed in the brain, suggesting that efficient coding could shape such transformations throughout the auditory system. The learned mid-level features provide hypotheses for auditory cortical tuning as well as a means to parameterize stimuli with which to probe mid-level audition.

Appendix B
In this appendix we describe the methods used to estimate and analyze receptive fields of model units.   Figure B1: Estimation of spectrotemporal receptive fields of model units. A) Top: an example first layer STK (leftmost column) and stimuli which strongly activate it (second column from the left). The "spike triggered average" i.e., the mean of strongly activating stimuli (third column from the left) differs strongly from the true underlying kernel (rightmost column) due to stimulus correlations. Bottom: same as top row for an example second layer unit. Second column from the left depicts strongly activating positive stimuli. Spectrotemporal pattern in rightmost column is derived from the procedure shown in (B). B) Generative estimation of second layer spectrotemporal receptive fields. To avoid the influence of stimulus correlations, we generate a number of STK activation trajectories (second column from the left; gray and blue lines plot two different sets of sampled trajectories) from a distribution encoded by a second-layer unit of interest (leftmost column). For each set of sampled trajectories, we convolve the trajectories with the corresponding STK and sum the results to yield a sampled spectrotemporal pattern (2 example patterns are shown in the third column from the left, resulting from the gray and blue samples, respectively). The set of such spectrotemporal patterns are then averaged to estimate the "true" receptive field.

Receptive field estimation
The "true" spectrotemporal pattern encoded by a STK is simply the STK itself. To identify average spectrotemporal pattern encoded by a second-layer unit we take advantage of the fact that our model is generative (Fig. B1). For each unit we generated a number of STK coefficient trajectories (second panel from the left, Fig. B1 B) from the distribution encoded by that particular unit (e.g. leftmost panel, Fig. B1 B). We then generated spectrotemporal patterns by convolving the trajectories with the corresponding STKs and summing the results together (see Fig.2). Due to the stochasticity of the generated STK trajectories, some variability is visible among the generated spectrotemporal patterns (third panel from the left, Fig.B1 B). We then averaged these patterns to form an estimate of a spectrotemporal receptive field of a second layer unit (right panel, Fig.  B1 B). Because the estimation process does not involve averaging stimuli, we avoid biasing the estimate of the average preferred spectrotemporal pattern with stimulus correlations.
It is apparent that "spike-triggered averages" (third column, Fig. B1 A and B) differ from the true stimulus representation in the model for both layers (as embodied by the STK or the average spectrotemporal pattern generated by a second-layer unit). The discrepancy is due to presence of strong correlations in natural stimuli.

Feature Selectivity Index
For comparison with neural data, we computed the feature selectivity index (FSI) proposed in [Miller et al., 2001]. The FSI is a number lying in the [0, 1] interval. FSI values close to 1 imply that stimuli eliciting a response of a neuron (or, in our case, a model unit) are similar to each other (specifically, the stimuli are close to the mean stimulus eliciting a response). When the FSI value is close to 0, the corresponding neuron spikes at random i.e. stimuli preceding spikes are uncorrelated with the spike-triggered average. The relevance of the FSI for our purposes is that a neuron or model unit that exhibits invariance to some type of stimulus variation should have an FSI less than 1. By comparing the FSI across layers we hoped to quantify differences in the degree of representational abstraction.
The FSI computation procedure are described in detail in [Miller et al., 2001]. The only discrepancy between the use of the FSI here and its prior use in neurophysiological studies is that units in our model are continuously active, and do not discretely spike. We emulated the selection of stimuli eliciting spikes by selecting the stimulus excerpts yielding the highest activation of model layer units. In the first layer, we selected the 25 stimuli yielding highest positive activation. In the second layer, we separately computed FSI indices for stimuli eliciting positive and negative responses, using 25 stimuli per unit in each case. We then averaged the indices obtained for the two sets of stimuli for each unit.
Computing FSI for an i-th unit consists of the following steps: 1. Compute average of strongly activating stimuli (analogous to spike-triggered average -STA) -separately for positive and negative stimuli in the case of second-layer units.
2. Compute correlations c S,i,n between each strongly activating stimulus n and its respective STA.
3. Compute correlations c R,i,n between the STA and a randomly selected subset of stimuli n.

Appendix C
In order to verify the correctness of the learning algorithm, we generated a toy dataset using 9 second-layer kernels, out of which 3 were "opponent" (Fig.C1 A, left panel). We generated sample STK activations by randomly superimposing the second-layer kernels to create a "variance map" (Fig.C1 B), from which we sampled STK trajectories (Fig.C1 C). Each iteration of the learning procedure made a gradient step on the model parameters using a single sample of STK trajectories (from a single sampled variance map). After 1500 iterations (Fig.C1 D), the model converged on the solution visualized in the right panel of Fig.C1 A. The generative kernels were recovered up to permutation, sign changes and in some cases, a temporal shift. These latter differences are expected due to the inherent arbitrariness of positive and negative signs in the model and the shift-invariance inherent to a convolutional code. These results demonstrate that the learning algorithm is capable of converging to the correct solution when the data is well-described by the model. -+ 0 max Figure C1: The model recovers ground truth second-layer kernels from a training dataset. A) Left panel -nine ground truth second-layer kernels used to generate the toy dataset. The three kernels in the top row are "opponent". Right panel -the generative second-layer kernels recovered by the model. They match the ground truth kernels up to permutation, sign change and time shifts. B) Example variance map generated by the model from the second-layer kernels. Such maps result from superimposing second-layer kernels with coefficients sampled from their sparse prior. This example was generated from three second-layer kernels from panel (A). The variance map encodes temporal changes in the distribution of STK activations. The toy training dataset was generated from 1500 such sampled maps (one per iteration of the learning algorithm). C) STK activations sampled from the variance map in panel B. In the actual training procedure, a single sample was used for each sampled variance map, but here we show 9 examples to illustrate the stochastic nature of the procedure. D) The convergence of the learning algorithm over 1500 learning epochs. The gray line depicts values of the second-layer energy function (see Eq.10). The black line is a 10-epoch moving average.


extfieldsofstudy: Medicine Biology

text:
Stabilization of HIF-2α impacts pancreas growth

Hypoxia inducible factors (HIFs) are critical regulators of the response to oxygen deficiency by activating target genes involved in a variety of biological functions. HIFs have been implicated in the pathophysiology of numerous pathologies including cancer. Patients with mutations in the von Hippel-Lindau (VHL) gene, an essential regulator of HIF activity, develop tumors in several organs including the pancreas. Previous functional studies of HIF activation in the pancreas have used Vhlh (the murine homolog of VHL) deficient mice. However, the role of each specific HIF transcription factors in the pancreas has not been thoroughly examined. We derived mice that constitutively express a normoxia-stable form of HIF2α in the pancreas. Activation of HIF2α in the pancreas severely impairs postnatal exocrine pancreas. Mice with pancreas-specific activation of HIF2α develop histological features reminiscent of pancreatitis including loss of acinar cells, ductal dilation and fibrosis. Moreover, we provide evidence that signaling pathways important for acinar cell homeostasis are altered in HIF2α-overexpressing pancreata.

the activation of HIF pathway on pancreas formation. To this end, we derived mice that constitutively express HIF2α in the pancreas. We focused on HIF2α given the prominent role of this HIF factor on proliferation, apoptosis and tumor formation 23 . Also, HIF2α appears to play a critical role in pancreas embryonic development 24 . We show that activation of HIF2α in the pancreas severely impairs postnatal exocrine homeostasis. Mice with pancreas-specific activation of HIF2α display histological features reminiscent of pancreatitis including loss of acinar cells, ductal dilation and fibrosis. Moreover, we provide evidence that signaling pathways important for acinar cell homeostasis are altered in HIF2α-overexpressing pancreas.

Results
Specific activation of HIF2α in the pancreas. To investigate the role of HIF2α in pancreas formation, we specifically activated HIF2α in the pancreas by crossing mice expressing Cre recombinase under control of the pancreatic and duodenal homeobox 1 (Pdx1) promoter (Pdx1-Cre mice) 25 with HIF2dPA mice (Pdx1-Cre;HIF2dPA mice). In HIF2dPA mice 26 , a form of human HIF2α that escapes recognition by the von Hippel-Lindau tumor suppressor protein (i.e. it is not degraded under normal oxygen conditions) is expressed upon Cre recombination. Pdx1-Cre;HIF2dPA mice were born at the expected Mendelian ratio and reached adulthood without any sign of compromised health. No differences in total body weight were found between Pdx1-Cre;HIF2dPA mice and littermate control mice at both 2 and 8 weeks of age (Fig. 1A). However, a substantial reduction in total pancreatic mass of Pdx1-Cre;HIF2dPA mice was found as early as 2 weeks of age, even when pancreatic weight was normalized to body weight (Fig. 1A). Western Blot (Fig. 1B) and immunofluorescence analysis (Fig. 1D) confirmed the efficient accumulation of HIF2α in Pdx1-Cre;HIF2dPA pancreas. In control pancreas, endogenous HIF2α expression was found in islets but not in the exocrine compartment (Fig. 1C), as previously reported 20 . Exocrine cell atrophy and ductal dilation in postnatal Pdx1-Cre;HIF2dPA pancreata. Gross morphology inspection and analysis of Haematoxylin and Eosin-stained tissue of newborn (P0) Pdx1-Cre;HIF2dPA pancreata did not reveal any apparent defects (Fig. 1E,F). However, loss of acinar cells and increased dilation of duct-like structures was observed at 2 weeks of age (Fig. 1I,J). These lesions rapidly progressed, and by 8 weeks of age most of the Pdx1-Cre;HIF2dPA acinar tissue had been replaced by duct-like tubular structures, stroma and adipose tissue (Fig. 1M,N). To further characterize the pancreatic exocrine defects of Pdx1-Cre;HIF2dPA mice we performed double immunohistochemistry for a ductal marker, duct-specific cytokeratin 19 (KRT19), and an acinar cell marker, amylase. The Pdx1-Cre;HIF2dPA pancreata appeared normal at P0 (Fig. 1G,H) but beginning at 2 weeks of age amylase expression was substantially reduced concomitant with expansion of duct-like tubular structures (Fig. 1K,L). The acinar architecture appeared disorganized in two-week-old Pdx1-Cre;HIF2dPA mice, marked by the loss of the typical round morphology of the acini as well as the dilation of the intracinar lumen (Fig. 1L) and concomitant increase of KRT19-positive epithelial cells. This acinar-ductal metaplasia, broadly defined as replacement of acinar cells by ductal-like cells, did not seem to involve transdifferentiation of acinar cells into ductal-like cells since intermediate cells co-expressing acinar and ductal markers were not observed ( Fig. 1L and Supplementary Fig. 1). By 8 weeks of age, the Pdx1-Cre;HIF2dPA pancreata appeared severely atrophic with extensive loss of acinar cells and dilated ducts; only residual clusters of acinar cells could be observed (Fig. 1O,P). Islets of Pdx1-Cre;HIF2dPA mice displayed a morphology similar to those observed in control mice ( Supplementary Fig. 1). Indeed, blood glucose levels in Pdx1-Cre;HIF2dPA mice were normal compared to control mice ( Supplementary Fig. 1). These results indicate that ectopic HIF2α stabilization in the pancreas leads to extensive loss of acinar cells and expansion of duct-like cells.
HIF2α activation leads to pancreatic abnormalities reminiscent of pancreatitis. Loss of acinar cells and ductal expansion are histological hallmarks of pancreatitis. We examined whether other classical features of pancreatitis were also present in Pdx1-Cre;HIF2dPA postnatal pancreata. Haematoxylin and Eosin staining showed increased stroma in 8-week-old Pdx1-Cre;HIF2dPA pancreata ( Fig. 1M,N). Gomori trichome staining revealed an marked increase in connective tissue in Pdx1-Cre;HIF2dPA mice ( Fig. 2A,B). Pdx1-Cre;HIF2dPA pancreatic tissue exhibited a substantial increase in cells expressing smooth muscle actin (SMA), a fibroblast activation marker (Fig. 2C,D) and vimentin, a mesenchymal marker (Fig. 2E-G), thus confirming the formation of pancreatic fibrosis. An increase in inflammatory cells was also apparent in 8-week-old Pdx1-Cre;HIF2dPA mice, as evidenced by CD11b immunohistochemistry, a marker of leukocytes (macrophages, neutrophils and granulocytes) (Fig. 2H-J) and F4/80 immunohistochemistry, a specific marker of macrophages ( Fig. 2K-M). Thus, HIF2α activation results in fibrosis and increased inflammatory cells in the pancreas, histological features of pancreatitis.
Increased acinar cell apoptosis and proliferation in Pdx1-Cre;HIF2dPA mice. To elucidate the underlying cause of the massive acinar loss in Pdx1-Cre;HIF2dPA mice, we assessed acinar cell proliferation, apoptosis and differentiation in two-week-old mice, an age in which Pdx1-Cre;HIF2dPA pancreata exhibit areas with preserved acinar tissue. Apoptosis was measured using TUNEL assay. Pdx1-Cre;HIF2dPA pancreata showed a marked increase in apoptotic cells compared to control pancreata ( Fig. 3A-C). Immunohistochemistry for cleaved caspase-3 confirmed the increase in apoptotic cells of Pdx1-Cre;HIF2dPA pancreata ( Fig. 3D,E). Cell proliferation was assessed by Ki-67 immunohistochemistry. Control mice displayed substantial pancreatic cell proliferation, as expected at this early postnatal age. However, both acinar (Fig. 3F,G,J) and ductal cell (Fig. 3H,I,J) proliferation was significantly higher in Pdx1-Cre;HIF2dPA pancreata. Thus, despite extensive acinar cell loss, pancreatic cell proliferation was prominent in Pdx1-Cre;HIF2dPA mice.
Remodeling of pancreatic acini upon HIF2α activation. Although the Pdx1-Cre;HIF2dPA pancreata appeared largely unaffected at birth, a substantial replacement of acinar cells by cells expressing ductal markers was observed as early as two weeks of age. We analyzed whether acinar cell dedifferentiation could also play a role in this postnatal exocrine pancreatic degeneration. Close examination of acini in two-week-old Pdx1-Cre;HIF2dPA mice revealed that dilation of intracinar ducts was associated with reduction of amylase immunoreactivity ( Fig. 1L and Supplementary Fig. 1). A key characteristic of the mature exocrine pancreas is the architecture of acini consisting of individual acinar cells that exhibit apical-basal cellular polarity and secrete digestive enzymes into a central lumen. Disruption of acinar cell organization and apical-basal polarity has been associated with postnatal pancreatic exocrine degeneration in various mutant mouse strains [27][28][29] . Thus, we decided to analyze the accumulation of several markers of acinar cell organization in Pdx1-Cre;HIF2dPA mice at two weeks of age when substantial acinar tissue was still preserved. The localization of the basal marker laminin (Fig. 4A,B) as well as basolateral markers β-catenin (Fig. 4C,D) and E-cadherin ( Fig. 4E-L) was not largely affected in Pdx1-Cre;HIF2dPA acini. In control acinar cells, carboxypeptidase A1 (CPA1) is localized close to the lumen of the acini (Fig. 4E). However, CPA1 did not show a clear apical localization in Pdx1-Cre;HIF2dPA acinar cells (Fig. 4F). Immunofluorescence for apical markers mucin-1 (Fig. G,H) and PKCε (Fig. 4I,J) confirmed the mislocalization of apical markers in Pdx1-Cre;HIF2dPA mice. Mucin-1 and PKCε were localized exclusively at the apical surface of control acinar cells (Fig. 4G,I). While these proteins still displayed an apical localization in Pdx1-Cre;HIF2dPA acinar cells, a marked cytoplasmic accumulation was observed (Fig. 4H,J). Immunofluorescence for mucin-1 also revealed the dilation of acini lumens in Pdx1-Cre;HIF2dPA mice ( Fig. 4H) compared to the small lumens found in control acini (Fig. 4G). The increase in lumen size of Pdx1-Cre;HIF2dPA acini was confirmed by immunofluorescence analysis of the tight junction protein ZO-1 (Fig. 4K,L). Therefore, HIF2α stabilization leads to extensive acini disorganization including intracinar lumen dilation and mislocalization of apical markers while basal markers are appropriately localized.  VHL patients commonly exhibit pancreatic lesions that include cysts and microcystic adenomas 10,11 . These lesions display prominent fibrous stroma and endothelial cells 11 . Also, epithelial cells lack the characteristic papillary structures and mucin found in pancreatic mucinous cystic neoplasms 11 . Similar to the pancreatic lesions of VHL patients, the duct-like tubular structures of Pdx1-Cre;HIF2dPA pancreata did not exhibit mucin-producing papillary structures, as revealed by alcian blue (Fig. 5A,B) and PAS staining (Fig. 5C,D). Indeed, Pdx1-Cre;HIF2dPA duct-like epithelial cells expressed markers of ductal cells including KRT19 (Fig. 1L,P) and duct binding lectin Dolichos biflorus agglutinin (Fig. 2C,D). Also, cilia, which are normally present in control ductal cells 30,31 but not in pancreatic intraepithelial neoplasia (PanIN) lesions and pancreatic cancer cells 32 , were found in pancreatic duct-like epithelial cells of Pdx1-Cre;HIF2dPA mice ( Fig. 5E,F). Pdx1-Cre;HIF2dPA pancreata displayed increased vascularization (Fig. 5G-K). MECA32-positive endothelial cells were closely intermixed with mucin-1-positive epithelial cells (Fig. 5I,J), as described in pancreatic cysts of VHL patients 11 . In agreement with the increased vascularization, an upregulation of the expression of the HIF-target gene VEGF was found in Pdx1-Cre;HIF2dPA pancreata ( Fig. 5K). Thus, HIF2α stabilization in the pancreas leads to the formation of pancreatic lesions similar to those found in the pancreas of VHL patients.

Ductal metaplasia and increased vascularization in
Altered signaling pathways in exocrine pancreas upon HIF2α activation. The expression of ductal transcription factors such as Sox9 is increased during acinar-ductal metaplasia in models of pancreatic injury [33][34][35] . In agreement with this, we observed prominent Sox9 expression in both disorganized acini and well-defined duct-like tubular structures of Pdx1-Cre;HIF2dPA pancreata (Fig. 5L,M). The MAPK signaling pathway plays an important role in acinar cell homeostasis in response to injury. During pancreatitis the MAPK pathway is up-regulated 36,37 and it has been shown to be necessary for acinar-ductal metaplasia 38,39 . We evaluated activation of the MAPK pathway, as assessed by immunofluorescence of phosphorylated ERK (pERK). In control pancreata, only a few scattered cells were positive for pERK (Fig. 6A, top panels). In stark contrast, strong nuclear pERK signal was observed in Pdx1-Cre;HIF2dPA acinar cells. This activation of MAPK pathway was largely limited to disorganized acini (Fig. 6A, bottom panels), the preserved acini in Pdx1-Cre;HIF2dPA mice did not shown substantial pERK signal (Fig. 6A, bottom panels). Of note, prominent MAPK activity was also observed in the stromal cells of Pdx1-Cre;HIF2dPA pancreata (Fig. 6A, bottom panels). The AKT/mTOR pathway has also been implicated in acinar-ductal metaplasia [40][41][42] . Thus, we evaluated the phosphorylation status of the ribosomal protein S6 (RPS6), a downstream target of the AKT/mTOR signaling pathway. As previously shown 42,43 , RPS6 was phosphorylated in acinar (Fig. 6B, top panels), but not ductal cells of control pancreata (Fig. 6C, top panels). In the pancreas of Pdx1-Cre;HIF2dPA mice, phosphorylated RPS6 (pRPS6) was easily detected in well-preserved acini (Fig. 6B, bottom panels). However, a dramatic decrease in pRPS6 was observed in areas with disorganized acini (Fig. 6B, bottom panels). Interestingly, we observed an induction of the stress protein clusterin in disorganized acini (Fig. 6D, bottom panels). No expression of clusterin was detected in normal acini of Pdx1-Cre;HIF2dPA mice (Fig. 6D, bottom panels) as well as in acini of control mice (Fig. 6D, top panels). Thus, HIF2α activation in the pancreas affects signaling pathways important for exocrine pancreas homeostasis.

Discussion
Here, we show that specific activation of HIF2α in the pancreas severely disrupts homeostasis of the exocrine pancreas. We derived transgenic mice that specifically expressed a normoxia-stable form of HIF2α in the pancreas. Continuous activation of HIF2α in pancreatic cells resulted in severe abnormalities in the postnatal pancreas including loss of acinar cells, acinar-ductal metaplasia, increased vascularization, and fibrosis.
To specifically activate HIF2α in the pancreas, we used the Pdx1-Cre transgenic line. Pdx1 is expressed during early stages of the developing pancreatic epithelium and HIF2α is expressed in pancreatic progenitor cells 24 . Indeed, HIF2α inactivation causes pancreatic hypoplasia indicating an important role for HIF2α in pancreas development 24 . However, newborn Pdx1-Cre;HIF2dPA pancreata did not display apparent pancreatic abnormalities thus indicating the HIF2α overactivation do not overtly affect embryonic development of the pancreas. After birth, a gradual loss of acinar cells was observed in Pdx1-Cre;HIF2dPA pancreata. Indeed, Pdx1-Cre;HIF2dPA mice develop features reminiscent of chronic pancreatitis including acinar cell loss, acinar-ductal metaplasia, inflammatory cell infiltration, fibrosis and lipomatosis. Despite the extensive acinar cell loss, Pdx1-Cre;HIF2dPA mice did not display loss of body weight after birth. However, it should be noted that the loss of pancreatic enzyme synthesis needs to be almost complete to cause overt physiological defects 44,45 .
Concomitant with the acinar cell loss, a marked increase in duct-like tubular structures was observed in Pdx1-Cre;HIF2dPA pancreata. Acinar-ductal metaplasia is commonly observed in chemical models of pancreatic injury as well as in several mutant mice deficient in genes critical for exocrine cell homeostasis and it is considered an adaptive response to acinar cell injury 46,47 . Although increased acinar cell proliferation was prominent in Pdx1-Cre;HIF2dPA mice, this did not result in tissue regeneration. Indeed, two-month-old Pdx1-Cre;HIF2dPA mice displayed a massive loss of acinar cells indicating that sustained HIF2α was detrimental to acinar cell regeneration. Lineage tracing experiments have demonstrated transdifferentiation of acinar cells into ductal-like cells in several animal models of acinar cell injury (reviewed in 47,48 ). In our study, lineage tracing could not be utilized to determine the acinar origin of ductal-like cells in Pdx1-Cre;HIF2dPA pancreata because the Pdx1-Cre transgene is expressed in common pancreatic progenitors. We did not detect intermediate cells co-expressing acinar and ductal markers in Pdx1-Cre;HIF2dPA pancreata that could suggest acinar-duct transdifferentiation but this is just circumstantial evidence. We did observe increased cell apoptosis in postnatal Pdx1-Cre;HIF2dPA pancreata. Thus, apoptosis could account for the severe loss of acinar cells although we cannot rule out an additional contribution from transdifferentiation of acinar cells into ductal-like cells. Although HIF2α usually promote cell survival and proliferation, a role of HIF2α in inducing cell apoptosis has also been reported. For example, HIF2α overexpression in hepatocellular carcinoma 49 and glioblastoma cells 50 increases tumor cell apoptosis. The proapoptotic effect of HIF2α overexpression on acinar cells could also be indirect as a response to acinar cell injury, as observed in pancreatitis. In agreement with this notion, we observed expression of the injury marker clusterin in Pdx1-Cre;HIF2dPA pancreata.
The duct-like tubular structures of Pdx1-Cre;HIF2dPA mice are reminiscent of the cysts and microcystic adenomas commonly found in VHL patients 10,11 , which are characterized by increased vascularization and prominent fibrous stroma. However, most pancreatic lesions in VHL patients are asymptomatic and no extensive acinar cell loss is usually found. Only rare cases of pancreatitis have been reported and these appear to be secondary to obstruction of the main pancreatic duct of cysts. A potential explanation for the difference in pancreatic phenotypes between VHL patients and Pdx1-Cre;HIF2dPA mice is that HIFα activation (by VHL inactivation) occurs locally in VHL patients (second hit hypothesis of hereditary cancer) while broad HIFα activation in the pancreas is achieved in Pdx1-Cre;HIF2dPA mice. In agreement with this notion, pancreas-specific Vhlh (the murine homolog of VHL) inactivation in mice results in the formation pancreatic cysts, loss of exocrine pancreas and fat replacement 21 . Interestingly, these pancreatic lesions were observed only in older Vhlh mutant mice, not in young mice 21 , as we observed in Pdx1-Cre;HIF2dPA mice. This apparent discrepancy might be related to the specificity and/or levels of HIF2α activation obtained in our mouse genetic model. Nonetheless, it is important to note that Vhlh mutant mice display severe perinatal lethality 15,21 and thus, the analysis of pancreata in older Vhlh mutant mice was limited to few surviving mice 21 , which could have not achieved full HIF activation on exocrine pancreas. In this regard, it would be interesting to determine which specific HIF transcription factors are activated in the pancreatic lesions of VHL patients. Vhlh mutant mice exhibit impaired glucose homeostasis including severe hypoglycemia and defective insulin and glucagon secretion. Although we did not perform a complete characterization of glucose homeostasis, Pdx1-Cre;HIF2dPA mice showed normal fed glucose levels. This apparent discrepancy might be due to differences in the specific HIFα overactivated between Pdx1-Cre;HIF2dPA and Vhlh mutant mice. In support of this hypothesis, the abnormal glucose homeostasis of ß-cell-specific Vhlh mutant mice is HIF1α-dependent since deletion of HIF1α restores impaired glucose homeostasis 12,16 . Thus, our results suggest that HIF2α overactivation in the pancreas do not markedly alter glucose homeostasis.
The pancreatic defects of Pdx1-Cre;HIF2dPA mice are similar to those observed in mice with defects in primary cilia formation and it has been shown that pVHL contributes to cilia maintenance and stability in kidney cells [51][52][53] . Indeed, cilia are lost in kidney cysts from VHL patients 53 . To the best of our knowledge, it has not been reported whether pancreatic cysts in VHL patients lack primary cilia 54 . However, primary cilia were clearly present on the ductal cells of Pdx1-Cre;HIF2dPA pancreata indicating that HIF2α overactivation does not overtly compromise primary cilia integrity. Interestingly, it has been reported that only HIF1α (not HIF2α) activity is important for cilia formation 53 .
The loss of acini organization in Pdx1-Cre;HIF2dPA mice was associated with alteration of the MAPK and AKT/mTOR signaling pathways. These pathways play a pivotal role in acinar cell homeostasis in response to injury. Thus, MAPK signaling is required for the dedifferentiation of acinar cells into duct-like cells 39 . In concordance with this notion, upregulation of pERK signal was observed in acinar cells located in disorganized acini of Pdx1-Cre;HIF2dPA mice. Interestingly, we also observed prominent MAPK activity in the stromal cells of Pdx1-Cre;HIF2dPA pancreata, a phenomenon also reported in fibrosis associated to pancreas cancer formation 55 . Disorganization of acini structure was also accompanied by decreased levels of phosphorylated RPS6, a downstream target of the AKT/mTOR signaling pathway that has been implicated in acinar-ductal metaplasia [40][41][42] . Crosstalk between these signaling pathways and HIFs have been reported 23,56 , thus providing a potential causal link between HIF2α overactivation and impairment of acinar cell homeostasis.
Depletion of Ptf1a, a master regulator of acinar differentiation and function in adult mice also results in decreased levels of phosphorylated RPS6 in acinar cells concomitant with loss of acinar cell identity 43 . Indeed, the postnatal pancreatic exocrine degeneration observed in Pdx1-Cre;HIF2dPA mice is similar to defects observed in mutant mice deficient in transcription factors critical for acinar cell formation and identity such as MIST1, NR5A2, GATA6 and PTF1A 28,29,43,57 . Pdx1-Cre;HIF2dPA pancreata displayed features compatible with loss of acinar cell identity including decreased expression of amylase, disruption of acini organization, mislocalization of acinar cell apical markers and increased expression of ductal markers. However, these phenotypes could be secondary effects of acinar cell injury, as reported in animal models of pancreatitis. Temporal, synchronous activation of HIF2α in adult acinar cells would be necessary to elucidate the direct role of HIF2α in acinar cell homeostasis.
SCIENtIFIC REPORTS | (2018) 8:13713 | DOI:10.1038/s41598-018-32054-5 A potential role for HIF2α in pancreatic cancer formation has been recently described. HIF2α is highly expressed in preneoplastic lesions but is downregulated during malignant progression of pancreatic cancer 20 . HIF2α inactivation impairs the progression of cancerous lesions in the Kras-mediated mouse model of pancreatic cancer 20 . We did not observe tumor formation in Pdx1-Cre;HIF2dPA mice but our analysis was limited to relatively young mice (two-month-old mice). Future studies in genetic models of pancreatic cancer will help to determine the effect of HIF2α activation in pancreatic cancer.

Methods
Mice. Mice were housed in a pathogen-free facility at the Institute of Biomedicine of Seville (IBiS). All procedures involving experimental animals were performed in accordance with European and local animal welfare laws, guidelines and policies. All animal studies were approved by the IBiS-Virgen del Rocio Ethics Committee. Pdx1-Cre and Hif2dPA (Rosa26Sor Tm4(Hif2A*)Kael ) mice carrying a floxed allele of human HIF2a have been previously described 25,26 . Hif2dPA mice were purchased from Jackson Labs.
Histological and immunohistochemical analyses. Histological (including alcian blue, periodic acid-Schiff and Gomori's one-step trichrome staining) as well as immunohistochemical analyses were performed as previously described 58,59 . Fluorescent and bright-field images were captured using a BX-61 microscope (Olympus) and LTC SP2 confocal microscope (Leica). All photomicrographs shown are representative of at least 3 independent samples of the indicated genotype. Primary antibodies were used at the To quantify cell proliferation, following double immunofluorescence for the proliferation marker Ki-67 (RM-9061; Thermo Scientific) and amylase or KRT19, positive cells were then counted from 10 non-overlapping fields at 200x magnification from 3 control and 3 Pdx1-Cre;Hif2dPA mice. Cell apoptosis was detected using In Situ Cell Death Detection Kit (Roche) according to the manufacturer's protocol in 10 non-overlapping fields at 200x magnification from 4 control and 5 Pdx1-Cre;Hif2dPA mice. Pancreatic sections treated with DNAse I (Sigma-Aldrich) were used as positive control. Acinar cell apoptosis was quantified by averaging the number of TUNEL-positive cells per total number of amylase-positive cells. Morphometric analysis of vimentin, CD11b and F4/80 were performed by using ImageJ software. For CD11b and F4/80, the number of positive cells per field was calculated. The extent of vimentin accumulation was measured as the percentage of vimentin-positive area per total pancreatic area. For each group, 10 randomly selected low-power fields from pancreatic sections at least 60 µm apart in three mice were analyzed. Western Blot analysis. Tissue samples were lysed and homogenized in lysis solution (50 mM Tris-HCl [pH 7.4], 100 mM NaCl, 5 mM MgCl 2 , 0,1% SDS, 1% Triton X-100) and Complete Ultra protease inhibitors (Roche). Protein concentration was determined using the Pierce BCA Protein Assay Kit (Thermo Fisher Scientific). After SDS-PAGE separation, samples were transferred to PVDF membranes (BioRad), processed for immunoblotting against anti-HA (1:500, 1167423001; Roche) and anti ß-actin (1:1000, ab8226; Abcam) antibodies, developed using HRP-conjugated secondary antibodies (Jackson) and Western ECL Substrate (Bio-Rad). ImageQuant LAS 4000 Mini Gold (GE Healthcare) was used for imaging the blots.

RNA isolation and Quantitative Real Time PCR. RNA isolation was performed using RNeasy Kit
(Qiagen) with deoxyribonuclease treatment, following manufacturer's instruction. Total RNA purity and concentration was assessed using Nanodrop 2000 spectrophotometer (Thermo Scientific, Wilmington, NC, USA), and subsequently reverse-transcribed using random hexamer primers and the cDNA Omniscript reverse transcriptase (Qiagen). Quantitative PCR was performed using with SYBR Green PCR Master Mix (Applied Biosystems) using a 7900HT real-time PCR system (Applied Biosystems). Relative quantification of RNA levels was calculated using the ΔΔCt method. Cyclophilin A (peptidylprolyl isomerase A-Mouse Genome Informatics) was used as housekeeping gene 59 . The results were expressed as fold relative to levels in control pancreata (value of 1). Primer sequences for cyclophylin A and VEGF were as follows: CycloA-forward 5′-TCACAGAATTATTCCAGGATTCATG-3′; CycloA-reverse 5′-TGCCGCCAGTGCCATT-3′; VEGF-forward 5′-CCAGCACATAGGAGAGATGAGCTT-3′; VEGF-reverse 5′-TCTGTCTTTCTTTGGTCTGCATTC-3′.

Statistical analysis.
Significance was determined using two-tailed Student's t-test, and one-way ANOVA (post-hoc Tukey HSD test). P < 0.05 was considered significant. Data are presented as mean ± SD.