Interhemispheric synchrony of spontaneous cortical states at the cortical column level Kazunori O'Hashi*1, Tomer Fekete*1,2, Thomas Deneux1, Rina Hildesheim1, Cees van Leeuwen2, Amiram Grinvald1

Keywords: Neural mass; Orientation; Spontaneous activity; Visual cortex; Voltage sensitive dye;

1

Department of Neurobiology, the Weizmann institute of science, Rehovot 7610001, Israel

2

Laboratory for perceptual dynamics, KU Leuven, Leuven3000, Belgium

*

equal contribution

Corresponding author: Tomer Fekete - [email protected]

1

ABSTRACT In cat early visual cortex, neural activity patterns resembling evoked orientation maps emerge spontaneously under anesthesia. To test if such patterns are synchronized between hemispheres, we performed bilateral imaging in anesthetized cats using a new improved voltage-sensitive dye. We observed map-like activity patterns spanning early visual cortex in both hemispheres simultaneously. Patterns virtually identical to maps associated with the cardinal and oblique orientations emerged as leading principal components of the spontaneous fluctuations, and the strength of transient orientation states was correlated with their duration, providing evidence that these maps are transiently attracting states. A neural mass model we developed reproduced the dynamics of both smooth and abrupt orientation state transitions observed experimentally. The model suggests that map-like activity arises from slow modulations in spontaneous firing in conjunction with interplay between excitation and inhibition. Our results highlight the efficiency and functional precision of interhemispheric connectivity.

2

INTRODUCTION Spontaneous cortical activity exhibits a rich and highly non-random spatiotemporal variability, even under anesthesia (Arieli et al. 1996; Tsodyks et al. 1999;Kenet et al. 2003). In early visual cortex (area17/18) of anesthetized cat spontaneous activity exhibits instantaneous patterns that resemble stimulus evoked orientation maps (Kenet et al. 2003; Altavini et al. 2016). Since also in the awake state, evoked and spontaneous activity resemble each other at least to some extent (Fiser et al. 2004; Omer et al. 2008; Luczak et al. 2009; Fukushima et al. 2012; Mohajerani et al. 2013), activity patterns appear to be similarly constrained by the underlying neuronal circuitry (Gong and van Leeuwen 2009; Luczaket al. 2009; Berkes et al. 2011; Mohajerani et al. 2013; Xie et al. 2016). In explaining these spontaneous patterns, recurrent processing takes a prominent role. Recurrent processing is mainly viewed as a matter of intrahemispheric lateral connectivity, (e.g. Goldberg et al. 2004; Blumenfeld et al. 2006;Galán 2008).The question could, however, be raised: do spontaneous maps within each hemisphere arise independently of each other, or do they involve both hemispheres? If the latter is true, activity patterns will span retinotopically corresponding loci (at least near the transition zone between areas 17 and 18 in the case of the cat visual cortex; Rochefort et al. 2007). Previous studies have shown bilateral synchrony between corresponding cortical areas during spontaneous activity in human and animals. These studies used both electrophysiology, voltage-sensitive dyes imaging (VSDI) and fMRI (Delucchi et al. 1962; Raichle et al. 2001; Lu et al. 2007; Nir et al. 2008; Mohajerani et al. 2010). However, the important question whether synchrony exists also at the fundamental columnar level remained unresolved. Yet, one might expect the symmetry to occur at that level: Fast monosynaptic connections exist between corresponding areas in area 17 as well as in area 18, mainly along the 3

midline (Innocenti et al. 1995; Olavarria 2001), and severing callosal connectivity abolishes stimulus induced interhemispheric synchrony in the cat early visual cortex (Engel et al. 1991). Still, it remains unclear if these connections are specific and powerful enough to induce synchrony in spontaneous orientation states across both hemispheres. To answer this question we chose to apply voltage sensitive imaging (VSDI) using a new improved voltage sensitive dye, to the early visual cortex of anesthetized cat bilaterally under two conditions: 1) in absence of stimulation with the animals' eyes covered and 2) while presenting full-field oriented gratings. Imaging the visual cortex bilaterally allowed us to observe to what extent callosal connectivity is able to synchronize information flow at the columnar level. To better understand the mechanisms underlying the dynamics we observed, we first characterized the precise dynamics of the emergence of spontaneous map-like activity patterns and the switching from one orientation state to another. To address the question of how such dynamics arise, we propose a neural mass model that can reproduce them, and allows to benchmark how efficient the callosal connectivity must be to achieve the experimentally observed interhemispheric synchrony.

MATERIALS AND METHODS Animals All surgical and experimental procedures were approved by the Institutional Animal Care and Use Committee at the Weizmann Institute of Science (No. 03140610-3). All experiments were performed in strict accordance with NIH guidelines. Recordings were performed in four adult cats, initially anaesthetized with medetomidine hydrochloride (0.1 mg/kg, i.m) followed by ketamine hydrochloride (15mg/kg, i.m) sedation. Animals were placed in a stereotaxic frame immediately after venous canulation and tracheal intubation. They were artificially ventilated 4

with a 60:40% mixture of N2O and O2 containing 1.0-1.5% halothane during surgery and 0.81.0% during recording. ECG, EEG, end-tidal CO2 and rectal temperature were continuously monitored. Paralysis was induced with pancuronium bromide (0.2 mg/kg/h) during imaging. Both eyes were protected from drying by contact lenses with appropriate curvatures. A rectangular craniotomy (12 mm × 17 mm) was performed and the dura mater was removed. Through the cranial window, areas 17 and 18 in both hemispheres were accessible. The chamber was filled with 2% agar and sealed with a cover glass during imaging to suppress the cortical pulsation originating from respiratory and cardiovascular movements.

Voltage-sensitive dye (VSD) imaging A new improved voltage sensitive dye was used here. The cortex was stained for 2-3 h with this dye; RH2115.The full chemical name is: 4-(5-hydroxy-4-{5-[3-allyl-4-oxo-2thiothiazilidin-5-ylidene]-penta-1,3-dienyl}-3-ethyl-pyrazol-1-yl)benzenesulfonicacid,triethylammonium salt). A MiCAM ULTIMA (BrainVision Inc., Tokyo) fast large well VSD camera was used for data acquisition. The cortical surface inside a part of the rectangular cranial window was imaged with a tandem-lens macroscope arrangement, comprising two photographic lenses (the 50-mm lens having a 1.2f number), coupled front to front. Each pixel of the CMOS camera detected light from a small rectangular cortical domain of 102  102 m2. Each VSD imaging session (for both spontaneous and evoked activity) lasted 20.48 s, and 4096 frames were obtained at 200Hz. Visual stimuli, produced by a stimulus generator (VSG, Cambridge Research Systems, Rochester, UK), consisted of full-field oriented (0, 45, 90 and 135) gratings (0.15 cycles/degree). To obtain orientation maps, visual stimulation consisted of full-field gratings (presented for 50ms) alternated with an inter-stimulus interval (gray screen, presented for 150 ms

5

to both eyes). The gratings were cycled in a fixed order (0° → 45° → 90° → 135° → 0 ...). The screen background was kept isoluminant for the entire trial period, including the period of blank screen prior to stimulus onset.

Data analysis Numerous papers and reviews, many of them by our group, discussed the elaborate analyses and image processing of VSDI data in great detail since 1975. Here we outline it only briefly. For more details see a recent review and references cited therein (Grinvald et al. 2015). Bleaching artifacts were removed by fitting data with a sum of two exponentials (Arieli et al. 1995; Tsodyks et al. 1999; Kenet et al. 2003; Deneux and Grinvald 2016). Next, heart pulsations were removed by removing the ECG triggered average of each time series (Arieli et al. 1995; Tsodyks et al. 1999; Kenet et al. 2003). Breathing pulsations were removed following the method described in (Fekete et al. 2011). If time series analysis was carried out data were transformed into fractional units (∆𝐹/𝐹) by dividing each optical movie by the average frame, and subtracting 1. To carry out pattern analysis, data were temporally smoothed to the 0.2-8Hz band, and spatially high pass filtered with a 510m radius filter to correct for uneven illumination (Tsodyks et al. 1999; Kenet et al. 2003; Chavane et al. 2011). Next, principal component analysis (PCA) was carried out separately for each experimental condition. Due to the vast number of frames PCA was carried out sequentially, similar to batch ICA analyses (Calhoun et al. 2001): PCA was carried out on each consecutive 1024 optic frames. The 100 leading components from each batch were concatenated, and finally PCA was carried out again, resulting in a basis which we refer to as signal. The same procedure was applied to the difference between the raw data and the preprocessed data resulting in an orthonormal basis we refer to as noise.

6

The resulting bases were used for denoising the data: data were reconstructed using the signal basis, and to reduce shot noise, only the components explaining 80% of the variance were used. Signal components were projected out of the data in an semi-automatic procedure informed by noise components: The correlation coefficient between the signal and noise components was computed, and signal components which were spatially similar to noise components (r>0.3) were discarded after visual inspection to verify that they did not contain functional patterns. In one instance of evoked activity it was necessary to manually adjust threshold to 0.5. An example of the effect of the procedure is shown is Fig. S8

Derivation of functional maps Maps were computed for each trial using a general linear model (GLM): box car regressors were convolved with a custom tailored impulse response constructed via the SPM8 software (http://www.fil.ion.ucl.ac.uk/spm/software/spm8/). The resulting coefficient maps were then averaged across trials. Next, orthogonal orientation maps were subtracted to obtain differential maps

Correlation analysis Optical frames were compared to the evoked maps using the correlation coefficient as a measure of similarity. To enhance the contribution of functional signals to the resulting scores images were (i) masked to include only the regions that gave rise to strong mapping signals (ii) smoothed with a Gaussian isotropic kernel of HWHH of a radius of 51 m.

Approximating the orientation domain Although when collecting data only a small subset of possible orientations can be presented to the animals, the resulting functional maps can be used to interpolate over the entire range of possible orientation states: Let VH denote the differential map obtained by subtracting 7

the Horizontal orientation map from the vertical one, and Obl, the oblique differential pattern. In the ideal scenario, both patterns are orthonormal, therefore can serve as a basis for orientationlike functional patterns. To ensure that both patterns will be orthonormal a two step procedure can be carried out. First, an orthogonal pattern for each map can be derived using the Graham Schmidt method, i.e. VH  VH  VH , Obl Obl and similarly for the oblique map. Next the

original and modified patterns can be averaged to obtain two corrected orthogonal basis patterns, i.e. 𝑉𝐻 = 𝑉𝐻 + 𝑘𝑉𝐻 ⊥ and 𝑂𝑏𝑙 = 𝑂𝑏𝑙 + 𝑘𝑂𝑏𝑙 ⊥ where k is given by : k   1

1  VH , Obl

2

.

These patterns can be normalized, and then the orientation of an instantaneous pattern I at time t is given by: 𝐴𝑁𝐺 𝑡 = 𝑎𝑡𝑎𝑛2 𝐼, 𝑂𝑏𝑙 , 𝐼, 𝑉𝐻 . The Similarity Index, 𝑆𝐼 𝑡 , is given by the correlation coefficient of each optic frame to the best fitting orientation template, i.e. 𝑆𝐼 𝑡 = 𝑐𝑜𝑟𝑟 𝐼, 𝑉𝐻 𝑉𝐻 + 𝐼, 𝑂𝑏𝑙 𝑂𝑏𝑙, 𝐼

A neural mass model for spontaneous activity The general layout of the model is presented in Fig. 5, and in the RESULTS section entitled: A neural mass model reproduces spontaneous dynamics. The dynamics of the excitatory neurons is determined according to: 𝑗

𝑗

(1) 𝜏𝑒−1 𝑣𝑒 = −𝑣𝑒 + 𝑕𝑒 + 𝑔𝑒𝑒

𝑁 k 𝑘=1 𝑤𝑗𝑘 fe

j

− 𝑔𝑖𝑒 fi + 𝜁𝑗

𝑗

where 𝑣𝑒 is the membrane potential of the excitatory population within the jth column, 𝜏𝑒 is the excitatory membrane time constant, 𝑕𝑒 the resting potential. 𝑔𝑒𝑒 , 𝑔𝑖𝑒 are synaptic gain terms for the interactions within and between populations, 𝑤𝑗𝑘 the lateral connection between columns j j

j

and k, fe and fi are the instantaneous firing rates of the jth excitatory and inhibitory populations respectively and 𝜁𝑗 a noise term representing random firing in other areas (c.f. the thalamus). The dynamics of the inhibitory neurons is determined by: 8

𝑗

j

𝑗

(2) 𝜏𝑖−1 𝑣𝑖 = −𝑣𝑖 + 𝑕𝑖 + 𝑔𝑒𝑖 fe The membrane potential to firing transfer functions are: j

𝑗

j

𝑗

(3) fe = Θ 𝑣𝑒 , fi = Θ(𝑣𝑖 ) where Θ is the sigmoidal function, i.e.: 𝑓

𝑚𝑎𝑥 (4) Θ 𝑥 = 1+𝑒 −𝑟𝑥

where 𝑓𝑚𝑎𝑥 is the maximal firing rate a column can produce and 𝑟 is the slope of the function. Following (Blumenfeld et al. 2006) the synaptic coupling strength between two columns was determined by the respective preferred orientation (given in radians) and orientation selectivity of each column, and also by the distance between two columns. Both parameters derive from the polar map that encodes the angle preference and selectivity at each cortical site (Bonhoeffer and Grinvald 1993): If the animal is presented with several evenly spaced oriented gratings spanning the orientation domain (0-), the information from all the experimental conditions can be used to derive the following map: (5) 𝑃𝑥 ≝ rx eiθ x = K2

j K j=1 sx

e−iφ j

j

where sx is the amplitude of the optical response to an oriented grating of angle φj at cortical location x(a given pixel in the optical frames), θx ∈ [0° 180°] is the preferred orientation of the column at cortical location x and rx its selectivity index. The map can then be used to determine the coupling between each pair of columns, i.e.: (6) 𝑤𝑥𝑦 = rx ry cos θx − θy Gσ ( x − y ) where Gσ x = 2π1σ 2 𝑒 −𝑥

2 /2𝜎 2

and 𝜎 is given in mm. and determines the spatial extent of the

interactions. Polar maps were normalized by setting 𝑟𝑥2 = 1 (which in turn normalizes the magnitude 9

of synaptic weights). Approximations of individual orientation maps for any desired angle φ could be derived from the polar map using: 𝜑

(7) 𝑀𝑥 = rx cos⁡ (θx − φ)(Blumenfeld et al. 2006). Our choice of parameters for the model - and the logic behind them - was inspired by (Jansen and Rit 1995; Wendling et al. 2000) .

Unbiasing and biasing polar maps To begin with, the model is based on Polar Maps (PM) that are strictly isotropic (both in terms of orientation selectivity and angle preference), obtained by slightly modifying the empirical maps as follows: First the intensity values of differential orientation maps in both hemispheres were forced to lie on the same symmetric Cauchy distribution (this was achieved by curve fitting followed by bin-wise adjustment of values to transform the original distribution to the best matching one, and was needed in particular to homogeneously distribute large values). Then the polar map derived from these differential maps was corrected to be isotropic as described in(Blumenfeld et al. 2006), both within each hemisphere and globally (Fig. 5 A). Such isotropic maps could later be biased in terms of: 1) bias in the area of the map allocated to, for example, the cardinal orientations 2) increased selectivity of columns with preference to e.g. the cardinal orientations or 3) both. Bias in the area allocated to an orientation could be achieved by multiplying the orientation selectivity map rx according to preferred angle multiplying rx by 1 + 𝛼 ⋅ cos⁡ (4(𝜃𝑥 − 𝜃0 )), where 𝜃0 could be e.g., 0 to bias for the cardinal orientations. Biasing both the preferred angle and selectivity distributions could be achieved by multiplying e.g. the 0 and 90 orientation maps approximated from the isotropic map (eq. 7) by 1 + 𝛼. Finally, bias in preferred 10

angle while maintaining an equal distribution of selectivity could be achieved by applying the second method followed by the first, with 𝛼 matched to offset the bias in area.

RESULTS Spontaneous orientation map patterns appear synchronously in both hemispheres In general, as reported previously (Delucchi et al. 1962; Raichle et al. 2001; Lu et al. 2007; Nir et al. 2008; Mohajerani et al. 2010), spontaneous fluctuations of optical signals were highly synchronous across hemispheres. In addition we observed spontaneous activity patterns resembling orientation maps (Kenet et al. 2003) spanning both hemispheres in all four animals used in this study (Fig. 1A, Fig.S1,Fig.S2). To quantify the prevalence and significance of the occurrence of bilateral spontaneous orientation states we conducted frame by frame correlation analysis employing the evoked differential orientation maps (the 0°/90° and 45°/135° differential maps obtained by the subtraction of the average single condition maps) as templates. As a control, the same procedure was repeated after the template maps were flipped horizontally. Note that negative correlations denote correlation to the orthogonal templates (e.g. r=-0.5 for the0°/90° map is by definition r=0.5 for the 90°/0° map).The results are shown in Fig.1B. We found that 34% of the time the maximal correlation to either of the templates (the maximum of the absolute of the two correlation values) exceeded the (p<0.01) threshold derived from thecontrol distributions. Repeating the analysis using different control patterns yielded similar results - see Fig. S3. In order to compare within and between hemisphere synchrony, we also conducted similar correlation analysis for each hemisphere separately: indeed, if spontaneous orientation states are not synchronous between hemispheres, correlation values are expected to be higher for single hemispheres compared to both taken together. However the distributions of correlations 11

Figure 1: Spontaneously emerging orientation states across hemispheres. (A) Spontaneously emerging orientation states span both hemispheres. Left - the differential orientation map (0/90) derived using a general linear model (GLM). Middle - a single frame captured while animal's eyes were covered in a dark environment. Right, a single evoked frame. (B, C, D, F data pooled from all animals) (B) Distributions of correlations between spontaneous states and orientation maps across both hemispheres. Red: the correlations to the maps pooled across four animals (8 hemispheres). Blue: correlation to flipped functional patterns. (C) Distributions of correlations between spontaneous states and orientation maps in both hemispheres simultaneously (blue), or in the left or the right hemisphere alone (red and yellow). Correlations were pooled from data originating from 4 animals (8 hemispheres). Each distribution was fit to a Gaussian. (D) Difference between the distributions of correlations when both hemispheres are considered, or only one. The slightly higher proportion of low correlation values in the twohemisphere distribution could be explained by the particulars of the computation: Indeed it was reproduced when using synthetic patterns in which synchrony across hemispheres was the ground truth (see text for details) (E) Top: differential functional maps evoked through repetitive stimulation with oriented moving gratings. Bottom: principal components derived from spontaneous data. Emergence of leading principal components in the spontaneous data indicates the prominence of orientation like states in the ongoing dynamics. Moreover, the correlation to

12

the evoked pattern spans both hemispheres indicating the tendency for inter-hemisphere synchrony in the dynamics: the correlation of the components to the evoked patterns across both hemispheres was 0.81 (0°/90°) and 0.82. The correlation of the evoked and spontaneous patterns across the left hemisphere was 0.86 (0°/90°) and 0.85 while the correlations across the right one were 0.77 (0°/90°) and 0.8 respectively. (F) Distribution of instantaneous orientation states across both hemispheres pooled across four animals. The states corresponding with horizontal and vertical orientations were the most common. for single and double hemispheres turned to be nearly identical (Fig. 1C). While a small difference existed, it could be fully explained by the particulars of the analysis: namely the fact that chance values for high correlations tend to be larger when the number of points decreases (on average the number of pixels spanning both hemispheres was double that in single hemisphere). Indeed, we carried out the same correlation analysis on a synthetic dataset to check the effect of pixel number, where the spontaneous frames were replaced by the functional maps to which we added time varying spatially smoothed normally distributed noise (Gaussian 204m radius). Therefore we exploited the fact that the ground truth that in the synthetic dataset patterns had the same amount of synchrony within and across hemispheres. We then computed the correlations to the templates for each hemisphere on its own vs. both taken together. The difference between distributions of correlation in both or single hemispheres obtained on this synthetic dataset reproduced that of the original data (Fig.1D), strongly suggesting that synchrony was as high between hemispheres as within the single hemispheres. We also carried out principal component analysis (PCA) on our spontaneous data. In all four data sets components nearly identical to the functional maps emerged as leading principal components (Fig. 1E, Fig.S4).To quantify the extent to which the principal components (PCs) of the spontaneous data resembled the functional maps, we computed correlation coefficients between the PCs and the differential maps (the 0°/90° and 45°/135° maps). The values for the cardinal orientation ranged between 0.72-0.82 (mean 0.78), and between 0.60-0.82 (mean 0.71)

13

for the oblique map. Furthermore, in all four animals, the components associated with the cardinal orientations appeared before (that is explained more variance than) the component associated with the oblique map (45°/135°). The ratio between the eigen-values associated with the cardinal maps and the oblique ones was 1.09, 1.47, 1.23 and 1.12 for the four cats.

Orientation synchrony between hemispheres is as strong as within single hemispheres Next, each frame was projected on an ideal ring structure derived by correcting the cardinal and oblique maps to be orthonormal (see methods). This resulted in an instantaneous orientation state (between 0°-180°) and similarity index (SI; the correlation to the combination of the corrected maps according to the projection weighs) for each frame. As was reported for the single hemisphere (Kenet et al. 2003), the distribution of momentary orientations showed a strong bias toward cardinal orientations, whereas oblique states were the least common (Fig. 1F). The correlation of SIs computed for each hemisphere was 0.46 (𝑝 = 0). These computations of instantaneous orientation were used to compare the difference in orientation between and within the hemispheres during spontaneous activity (Fig.2A). We divided each hemisphere into posterior and anterior parts using masks of roughly the same size. We then derived the instantaneous orientation state and similarity index for each of the quadrants of the imaged cortex. The within and between hemisphere distributions of difference in quadrant orientation state were compared using the Smirnov Kolmogorov test, and were not statistically different. To rule out the possibility that underlying differences were masked by carrying the analysis on all spontaneous frames - which include numerous frames in which prominent orientation map-like states were absent - we repeated the analysis after first selecting frames according to quadrant SI. This was done for three threshold values - 𝑝 < 0.05, 𝑝 < 0.01 and

14

Figure 2: Both hemispheres are functionally coupled: (A,B) data pooled from all animals. (A) The two hemispheres were divided into four quadrants (half hemisphere). For each quadrant the spontaneous distributions of similarity indices and angle states were derived. The resulting four distributions were compared either within hemisphere (anterior left and right to posterior left and right respectively) or across hemisphere (anterior left to right, and posterior left to right). (B) Histogram of angle differences within and across hemisphere, restricted to optical frames for which at least one of the quadrants expressed a clear orientation map-like activity pattern (SI of p<0.01). The resulting distributions were not statistically different as determined by the Smirnov-Kolmogorov test. 𝑝 < 0.001, and again no significant differences were found (Fig. 2B shows the 𝑝 < 0.01 case). Similarly, the correlation between the within and between angle difference across all frames and animals was 0.54 ( 𝑝 = 0) . Altogether, no difference was found between intra- and interhemispheric synchrony of orientation representations.

Trajectories of instantaneous orientation state “crawl” and “hop”, and oversample cardinal orientations Next we characterized the dynamical aspects of orientation state switching by analyzing the time series of instantaneous SIs and orientations. We observed two characteristic types of salient events (Fig.3A): "crawling" that is prolonged stretches in which orientation states changed gradually, and "hopping", rather abrupt changes in orientation state in which the cortex tended to switch to an orthogonal orientation state (e.g. 45°→135°). Spontaneous orientation states tended to be transiently attracting (Fig.3B, C, D), that is 15

Figure 3: Ongoing dynamics of global orientation state as motion on a ring of transiently attracting states. (B-E) data pooled from all animals. (A) The dynamics of transient orientation state across both hemispheres during a single trial. As can be seen in the example the orientation state transitions exhibited two salient features: "crawling" - prolonged stretches of moderate change in orientation state - see e.g. the period marked by the grey bar, and "hopping" - sudden transitions to orthogonal states - marked by the red bar. (B) The rate of change in orientation state per optical frame as a function of time (all data analyzed from 4 animals concatenated). Similarity Index (SI, the correlation to the best matching state on the idealized ring) is color coded. As can be seen high SI values are associated with near zero change in orientation. (C) The distribution of orientation change rates shows that counter and clockwise trajectories on the ring are as likely. ~91.2% of changes in orientation state are smaller than 5 per frame (i.e. per 5ms). (D) The average instantaneous similarity index for different orientation change rates (error bars: s.e.). As in B it can be seen that stable states are associated with higher similarity to ring states. (E) Similarity to ring states as a function of instantaneous orientation. Cardinal orientations are associated with stronger SIs as compared to oblique patterns. Taken together with the frequency of instantaneous orientations (Fig. 1F), this suggests stronger attracting properties of cardinal orientations.

16

the stronger the state in is terms of its correlation to the best matching map, the more stable it tended to be: the instantaneous SI of a state was inversely proportional to the rate of change in instantaneous orientation (Fig. 3B, expressed in °per optical frame). Clock-wise and counter clock-wise changes in angle were just as likely (Fig.3C). Attraction strength was strongest for cardinal orientations and weakest for the oblique angles (Fig.3E). The of rate of angle change between hemispheres was highly correlated during the expression of spontaneous orientation states - see Fig. S5. The bias for cardinal orientations was not readily explained by the spatial layout of the maps in our sample: We found no correlation between the fraction of pixels whose preferred orientation lies in each of eight 22.5° wide orientation bins and either the fraction of the time spent in each angle state (Fig. 1F; r=-0.11, p>0.6, df = 22), or similarity index (Fig. 3E; r=-0.32, p>0.1, df = 22).

Orientation patterns emerge during depolarization events To explore the link between depolarization and emergence of orientation states, we used the local maxima in the SI time series that were in the top 10 percentile in terms of SI magnitude as triggers to average the optical signal around these events. As can be seen in Fig. 4A,B, on average, strong orientation states were preceded by slow depolarization spanning the entire imaged area (Fig. 4B, inset) and could persist provided that depolarization was sufficiently large. The model of Blumenfeld and colleagues (Blumenfeld et al. 2006) predicts that the emergence of mosaic states - mixed orientation patterns, e.g. a pattern composed of the 0° map in the anterior and the 45° map posteriorly (see Fig. S6A) -is contingent on the spatial extent of effective connectivity, that is, such states are ruled out by wide-reaching effective connectivity.

17

Figure 4: Emergence of spontaneous orientation states. (A-D) data were pooled form all animals. (A-B) Local maxima in the SIs of spontaneous orientation states were used as triggers for averaging (A) the resulting average SI suggests that stable orientation states tended to last around 200 msec. (B) The mean optical trace averaged around the triggers. It can be seen that emergence of strong spontaneous orientation states was associated with brief depolarizing events. As on average these events seem to span longer durations than the corresponding activated orientation maps this suggests that global depolarization is an enabling condition for spontaneous orientation dynamics. Inset: the average spatial spread of depolarization presented in z-scores computed relative to each time series (C-D) Local analysis: SIs and orientation states were derived for each pixel using10 pixel radius circular windows (see text for details). Using a varying SI threshold the local similarity maps were analyzed for clusters. Local maxima in cluster size was chosen for each threshold value such that there was no overlap between cluster of lower threshold and higher threshold (i.e. clusters thresholded by 0.35 strictly contained values smaller than 0.45). (C) The optical trace was averaged across each chosen cluster for a given threshold value and averaged. As can be seen, depolarization was positively correlated with SI. This means that for spontaneous orientation states to emerge locally, sufficient depolarization was necessary. Inset: the spatial spread of depolarization for the 0.55 threshold presented in z-scores (D) The circular standard deviation of orientation states in each cluster was computed. As SI increases the variance in orientation states within each cluster decreases showing that high SI values indicate coherent rather than mosaic states. Taken together (A-D) suggest that orientation states "ride" on widespread depolarizing events presumably either standing or slowly propagating waves.

To explore to what extent mosaic states were prevalent in our data, orientation state and SI were computed for each pixel using a circular neighborhood of 10 pixels. As can be seen is Fig.S6, mosaic states, like global orientation states, are expected to result in widespread large local SI values, but unlike global states also in large variance of the local orientation state over different locations due to coexistence of multiple orientation states simultaneously (e.g. bimodal partition 18

in the example in Fig. S6C). Next, in each optical frame we sought out clusters (connected components of pixels) with SI values in the [0.35 0.45) range, the [0.45 0.55) range, or above 0.55., choosing temporal local maxima in cluster size as triggers for averaging the optical frames. As can be seen in Fig. 4C, as was the case for global analysis, emergence of orientation states was triggered by depolarization that was proportional to SI magnitude. As was the case for global analysis, the spatial extent of strongly expressed orientation states (i.e. strong SI) spanned the entire imaged area (Fig. 4C, inset). We also computed the circular standard variation of the orientation states in each such cluster. As can be seen in Fig. 4D the S.D. decreased as SI increased, indicating that only globally coherent orientation states were associated with high SI values (i.e. orientation states were mostly coherent), which implies a wide effective connectivity. Taken together these analyses suggest that spontaneous orientation states "ride" on widespread depolarizing events - presumably either standing or slowly propagating waves(Benucci et al. 2007) and tend to span the entire early visual cortex in the two hemispheres. Several studies reported "replay" of evoked patterns observed in spontaneous activity following exposure to stimuli during sleep (Ji and Wilson 2006) and anesthesia (Han et al. 2008). This raises the potential concern that our results represent such replay rather than bona fide spontaneous dynamics. However, half of our ongoing data were collected prior to evoked data, and we did not observe differences in the reported properties of the spontaneous data before and after exposure to stimuli (see Fig. S7)

A neural mass model reproduces spontaneous dynamics Theoretical studies suggest that spontaneously emerging map like activity patterns could result from two mutually exclusive scenarios (Goldberg et al. 2004): 1) a "background state" scenario, in which in the absence of external stimulation the network hovers around a baseline

19

state, and 2) wandering around multiple attractor states (MASs), namely cortical maps which encode various visual attributes (e.g., orientation, ocular dominance, spatial frequency ...). The scenarios differ in their predictions in two critical regards: First, under the background state regime, the width of the distribution of correlations to the evoked maps is inversely proportional to the number of imaged columns K (~

1 𝐾

), i.e. correlation is expected to diminish as the size of

imaged cortex increases. In contrast, the MASs scenario predicts that the width of the distribution will remain invariant regardless of the size of the patch of imaged cortex (as long as it is within early visual cortex), given that during the expression of a map like activity pattern the correlation to the map will be similar for different sized parts of the image. In a similar vein, only under the multiple states regime, it is predicted that the principal components of the spontaneous data will highly resemble the encoded feature maps, such as orientation maps. This is due to the fact that under MASs, the network will dynamically switch between map-like attractor states, and therefore a large fraction of the variance in activity will reside in a low dimensional subspace associated with the maps. Our data are consistent with the MASs regime - the width of the distribution of correlations of spontaneous activity patterns to the functional maps shown here is similar to that reported by (Kenet et al. 2003), even though the area imaged in our study is more than 3 times larger. Further still, the width remains the same when comparing the single hemisphere case, to the distribution arising from both hemispheres taken together - a cortical area roughly double in size. Secondly, map like patterns emerged as principal components of spontaneous fluctuations. While the model of Goldberg et al. (Goldberg et al. 2004) fits our data along these regards, it predicts that in the MASs regime the network will exhibit a random walk on the manifold of attracting states. However, apart from "crawling" - extended stretches of wandering

20

through orientation states - the dynamics also exhibited "hopping", abrupt transitions in orientation state usually to the orthogonal state (e.g. 0°→90°), calling for a more comprehensive modeling approach. We chose to employ a neural mass approach to model the observed dynamics (for different applications of neural mass models to VSDI, see Markounikau et al. 2010; Trong et al. 2013). In our model, basic units represent "columns" within early visual cortex, each corresponding to one pixel in the imaged cortex (therefore capturing the activity of a few thousands of cells; note that here the term “column” refers to a processing unit that is smaller than an orientation column). Each such column comprises two populations - excitatory units and inhibitory ones. Excitatory and inhibitory pools are connected bilaterally within each individual column, and only the excitatory pools of different columns are connected laterally (see Fig. 5 for a schematic representation of the model). The strength of the lateral connections between columns is determined by the Polar Map (PM, Fig. 5. and see methods), which encodes the orientation preference (an angle between 0180°), and selectivity (a number between 0-1) of each recorded cortical site (pixel). The connection strength decreases monotonously with the difference in preferred angle, and is weighted by the orientation selectivity of both columns (e.g. will be larger for selective columns). In addition, the lateral connections are spatially attenuated according to distance (i.e. will tend to 0 for remote connections). The empirical orientation maps shown in Fig. 1 were used to derive the polar map utilized in our simulations. However they were standardized as follows to eliminate possible confounds that could result from uneven illumination and the partial sampling of the maps (relative to the entire map of each animal): First in each hemisphere, an equal region of the map,

21

Figure 5: A neural mass model for cat V1. (A) An experimentally derived Polar Map (PM) corresponding to the orientation maps shown in Fig. 1A. At each recorded cortical position x the map encodes the preferred orientation 𝜃𝑥 (color) and orientation selectivity 𝑟𝑥 (brightness). In the zoom in below two cortical sites (i.e. pixels) are marked, one with low orientation selectivity (x, a pinwheel), the second with high orientation selectivity (y). (B) Each cortical site is modeled as consisting of a local pool of inhibitory neurons (interneurons) mutually connected to an excitatory pool of (pyramidal) neurons (top left- represented as two neural mass equations in the model - see methods for details). Different cortical sites are connected laterally by their excitatory populations. The strength of these lateral connections is determined by the difference in preferred orientation between the columns, weighted by their respective orientation selectivity (the more selective they are, the stronger the connection), and modulated by their distance (the extent of lateral connectivity is determined by a Gaussian term whose width is controlled by the parameter ).

approximately 3x6mm in which the functional pattern was clearly observed was chosen. Second, orientation preference and selectivity values were slightly modified to yield equally distributed preferred orientations (which corresponds to equal area in the map dedicated to each angle) and selectivity indices, within each hemisphere and globally (Fig. 5 A, see Methods for details). To simulate spontaneous activity the network was driven with i.i.d. noise equally

22

distributed in a defined interval, i.e.: 𝜁(𝑡) ∈ 0, 𝑁 , and the equations were evaluated using the Euler method with ∆𝑡 = 5 𝑚𝑠𝑒𝑐. Following (Markounikau et al. 2010) we defined the VSD 𝑗

𝑗

signal in our network as: 𝑉𝑆𝐷 = 𝛼 ∙ 𝑣𝑒 + (1 − 𝛼) ∙ 𝑣𝑖 . This reflects that VSDI signals originate from membrane potential, and most of the membrane area in cortex belongs to cell bodies and dendritic trees. For simplicity's sake we employed the value 𝛼 = 0.5. Under noise stimulation for values of N smaller than 75mV the network exhibited Noise Oscillations, in which spontaneous orientation map like activity patterns were not observed. From that point onward, map like patterns begin to emerge with increasing probability. At relatively moderate values of 𝑁 "hopping" like transitions dominate the dynamics - rather abrupt switching from a map to its (near) orthogonal counterpart (see Fig.6A - in which N was 80mV−and compare to Fig.3A). As 𝑁 is increased further, short epochs of "crawling" - a gradual progression along the ring of orientation states either clockwise or counter-clockwise - emerge with increasing probability, until the dynamics are dominated by them (see Fig.6B in which 𝑁 was 120mV, and compare to 3A). As the empirical data exhibited both hopping and crawling, we reasoned that stimulating the network with time varying noise could lead to similar "mixed" dynamics. We therefore set 𝑁 to be time dependent and vary between 60 and 110mV. We explored three noise regimes, in which 𝑁was taken to vary in time within the following frequency bands: 0-0.3 Hz, 0-1Hz and 04Hz. Note that under such time-varying parameter N noise inputs to different columns become positively correlated (instead of uncorrelated as before), e.g. large values of N lead to large average noise values over all columns. Indeed, this regime gave rise to mixed dynamics, and modulations in the 0-1HZ range were found to best resemble the empirical transitions - see

23

Fig.6C (in which the connectivity was also biased towards cardinal orientations: see below as well and Fig.7).

Figure 6: simulation of spontaneous activity. The model described in the text was driven with noise representing spontaneous firing in other areas, to generate spontaneous activity. (A)The network was fed with i.i.d. noise lying in [0 80 mV] leading to dynamics dominated by hopping - abrupt switching from a map to its (near) orthogonal counterpart. (B) Driving the network with i.i.d equally distributed noise lying in [0 120 mV] led the dynamics to be dominated by crawling - gradual progression along the ring of orientation states either clockwise or counter-clockwise. (C) When N, the limit of the noise, was set to be a time dependent global parameter oscillating in the 0-1Hz range between 60 and 110 mV – representing slow wave activity across cortex and the thalamus – mixed dynamics qualitatively similar to the empirical dynamics emerged. In this example the connectivity was biased towards cardinal orientations (see text, Fig.7).

As noted above, spontaneous activity under anesthesia exhibited clear bias towards cardinal orientations (Fig. 2). We hypothesized that this could result from 3 scenarios: 1) bias in the area of the map allocated to the cardinal orientations 2) increased selectivity of columns with preference to cardinal orientations or 3) both. Starting from the isotropic map described above, we systematically biased the map to conform to one of these scenarios. We then proceeded to simulate the network with the slowly modulating noise described above. We observed that while 1) and 2) could systematically bias the resulting distribution of 24

spontaneous angle states (e.g. Fig. 1F), the resulting distribution of SIs (see Fig.3E) did not show any systematic effect. In contrast, when both preference and selectivity were biased to cardinal orientations (Fig.7 A, B) both SIs and spontaneous angle states reflected the chosen bias. Figure 7: biasing the emergence of spontaneously emerging maps towards cardinal orientations. An empirical polar map was biased in both angle preference and selectivity using a factor of 0.04 (see Methods for details and figures 1-2 for comparison to the empirical results), thus modifying the connectivity between units and the dynamics of simulated spontaneous activity: (A) Similarity Indices were larger on average for cardinal orientations.(B) Spontaneous VH states were more frequent than oblique states. (C) Biased functional maps (top row, the bias toward cardinal orientations appears in the stronger values of the 0°/90° map), leading principal components (PCs) computed from the resulting simulated activity (middle) and leading PCs of the biased connectivity matrix used in the simulation (bottom). As predicted by(Galán 2008), PCs of the activity and connectivity matrix are highly similar, and exhibit a similar ratio of the corresponding eigenvalues, in the range found for the empirical data.

We also computed the principal components (PCs) of this simulated activity. The leading components (typically the 1st and 2nd) were highly correlated to the orientation maps (Fig.7C), and the ratios between the respective eigen values (i.e. 𝜆𝑣𝑕 𝜆𝑜𝑏𝑙 ) were in the range found empirically (1.1-1.4), reflecting the imposed bias towards cardinal orientations. A previous modeling study by Galán (Galán 2008) suggested that in scenarios where spontaneous activity is driven by noise, it can be shown through a linear approximation that the PCs of the resulting data, and those of the connectivity matrix between units in the network are expected to be the same, and that their eigenvalues determine the probability of the emergence of the associated patterns. 25

While our model is more complex than the one described there, and the regime we focused on violates the assumptions enabling the linearization of the system at the basis of Galán's approximation, we nevertheless derived the PCs of the connectivity matrix. Indeed we found that the leading PCs were highly similar both to the spontaneous PCs and the approximated orientation maps. Moreover, the ratio between the eigenvalues associated with the leading map-like components were similar, suggesting that the theory developed by Galán is true for our model as well. More broadly, this suggests that cortical feature maps should be thought of as the PCs of the underlying patterns of connectivity, and as a consequence, because PCs are orthogonal by definition, that different cortical feature maps tend to be orthogonal to each other as a generic organizational principle (Kremkow et al. 2016; Lee et al. 2016).

Bilaterally emerging orientation states necessitate precise functional specificity of callosal connections Our correlation analysis suggests that, functionally speaking, both hemispheres are coupled as if they were a single hemisphere. This raises the question of the conditions the callosal connections must satisfy to enable such coupling. We could address this question by testing different modifications of our model: We manipulated selectively the L-R connections, generated spontaneous data accordingly, and studied how the angle distributions derived for each quadrant were modified and compared either within or across hemispheres. First, we asked to what extent L-R connections could differ in magnitude from lateral connections to maintain L-R coupling. Therefore, we systematically varied only the magnitude of inter-hemispheric connections in our model by multiplying them by a positive value smaller than one. This was followed by simulating 20 trials of 4096 frames for each gain factor, and comparing the half-within and half-across angle distributions using the Kolmogorov Smirnov

26

test after binning the angle data to a 250 bins. As can be seen in Fig.8A, even perturbations as small as 10% were sufficient to significantly decouple the hemispheres. This lends credibility to the notion that the strength of callosal connections is of the same order of magnitude as within hemispheric (lateral) connections in early visual cortex.

Figure 8: Interhemispheric (L-R) connections are strong, orientation specific and fast .Our model was used to examine the extent to which the emergence of bilateral spontaneous orientation states depends on the specific attributes of callosal connections. (A) Decreasing the magnitude of interhemispheric synapses by as little as 10% led to noticeable decoupling of the hemispheres suggesting that these connections must be of similar magnitude to lateral connections in early visual cortex. (B) Scrambling the targets of the LR connections completely decoupled the hemispheres, supporting the notion that callosal connections are orientation specific at the columnar level. (C) Introducing small delays selectively only to L-R connections resulted in noticeable decoupling even for small delays of 5 msec., suggesting that callosal fibers must be faster than those in early visual cortex to avoid introducing signaling delays.

By the same token, one would expect that the functional efficacy of the callosal connections would depend on their orientation specificity. We therefore conducted similar analysis after first randomizing increasing fractions of the inter-hemispheric connections in the model via permutation (leaving the total distribution identical), in increments of 25%. As can be seen in Fig.8B, this effectively decouples the hemispheres, suggesting that L-R connections of bilateral early visual cortex are orientation specific at the columnar level. Finally, one would expect that for such functional coupling to manifest, callosal axons must be faster than lateral connections within early visual cortex, due to the increased distance they must transverse. While our model does not directly incorporate synaptic delays, we tested this idea by introducing simple uniform delays only to the interhemispheric (L-R) synapses in

27

our model. As can be seen in Fig.8C even the small delays we explored - 5-20 msec. (corresponding to 1-4 time steps in the simulated dynamics) - could produce a significant effect, and while the effect was not as strong as the previous two, still we found a monotonic dependence of the KS statistic on the magnitude of the delay, lending support to the finding that callosal fibers are fast (Innocenti et al. 1995; Olavarria 2001).

DISCUSSION In this study we applied VSDI to early visual cortex (areas 17/18) bilaterally in anesthetized cat. Previous studies have shown bilateral synchrony during spontaneous activity in human and animal using electrophysiology, fMRI and VSDI (Delucchi et al. 1962; Raichle et al. 2001; Lu et al. 2007; Nir et al. 2008; Mohajerani et al. 2010).This is the first study that explored the bilateral synchrony at the columnar level. We found spontaneously arising activation patterns highly similar to evoked orientation maps spanning both hemispheres, demonstrating spontaneous, interhemispheric synchrony at the columnar level. Furthermore, inter-hemisphere synchrony was as high as intra-hemisphere synchrony, implying that interhemispheric connections are as functionally efficient and fast as intra-area connections. The dynamics of bilaterally synchronized activity: The precise characterization of spontaneous pattern dynamics informed us as to the nature of the underlying network: The strength of the transient orientation states was inversely proportional to the rate of change in angle, providing direct evidence that these maps are transiently attracting states. These dynamics also exhibited two prominent features: extended stretches in which orientation state changed gradually, or "crawling", but also "hopping", abrupt transitions in orientation state usually to the orthogonal state (e.g. 0°→90°), which we could attribute to time-varying levels of input noise. Relationship to previous theoretical models: Our results partially support the 28

theoretical outlook suggested in (Kenet et al. 2003)with regards to the observed spontaneous dynamics, which is referred to as wandering through multiple attractor states (Goldberg et al. 2004) (MASs): first as expected by this model, correlations of spontaneous activity patterns to orientation maps remained similar (the width of their distribution did not decrease) when doubling the area used for comparison. Second, as expected from MASs, the principal components of spontaneous fluctuations coincided with the evoked orientation maps: patterns virtually identical to maps associated with the cardinal (0°/90°) and oblique (45°/135°) orientations emerged as leading principal components of the spontaneous fluctuations. However, the model of (Goldberg et al. 2004) predicts that in the MASs regime the network should exhibit a random walk on the manifold of attracting states: Although we indeed observed extended stretches of wandering through orientation states ("crawling"), the dynamics also exhibited "hopping". Therefore, to explain these dynamics we took a more biophysically realistic modeling approach, employing a neural mass model with both inhibitory and excitatory pools of neurons, as well as incorporating a more realistic wiring scheme, taking into account not only columnwise orientation preference and selectivity, but also the spatially restricted nature of connections (Blumenfeld et al. 2006). Mechanistic insight from the present model: When driving our model with slowly modulating noise (<1Hz), it reproduced the salient features of the data: the emergence of spontaneous orientation states forming a ring of transiently attracting states, emergence of maplike patterns as principal components, as well as both crawling and hopping state transitions. Even

though

previous

studies

hypothesized

that

synaptic

mechanisms

leading

to

depression/facilitation of lateral connections are likely at play (Blumenfeld et al. 2006; Romani et al. 2012), our model demonstrates that the interplay between excitation and inhibition is

29

sufficient for such dynamics. Further still, our model allowed us to specifically understand which properties the callosal connections need to have, in order to support the bilateral similarity in orientation state. Our model in conjunction with the data suggests that callosal connections must be fast (to cover the added distance compared to within hemispheric lateral connections; Innocenti et al. 1995; Olavarria 2001), functionally precise (orientation specific; Rochefort et al. 2007) and similar in synaptic efficacy compared to within hemispheric connectivity. A previous imaging study in humans showed that spontaneous bilateral synchrony at the area level persists in callosal agenenis (Tyszka et al. 2011). However, this finding is not relevant to orientation specific synchrony in cats or monkeys, given that orientation specificity is achieved in cortex rather than lower brain regions (Kremkow et al. 2016; Lee et al. 2016). Recent studies also lend strong support to our interpretation that the callosal connections are involved in establishing the bilateral columnar synchrony (Altavini et al. 2016). Future experiments can prove this point by repeating the described experiments after severing callosal connections. During anesthesia (and slow wave sleep) cortex is induced into slow wave states, a phenomenon known to occur at slow rates (<1Hz and delta). These oscillations originate from intrinsic inter cellular dynamics (up and down states) induced by neuromodulators, and are likely synchronized though thalamo-cortical loops (Steriade et al. 1993). Introduction of time varying noise at 1Hz in effect mimics this phenomenon. Therefore, our results suggest that slow oscillations are a necessary (but not sufficient) condition for the emergence of the spontaneous map dynamics found in visual cortex. We believe that the fact that the simulated dynamics preserve this salient feature of the intrinsic dynamics adds to the ecological validity of the conclusions drawn from our simulation regarding the requisite functional properties of callosal

30

connections. The above hypothesis can be checked by imaging spontaneous activity employing either different levels of anesthesia, or perhaps different anesthetic protocols, which would differ in the extent (depth of modulation) of slow waves. Such future experiments are warranted. We predict that regimes that will induce only shallow modulations will exhibit only hopping transitions. Bias towards the cardinal orientations states: As previously reported, our data exhibited bias towards the cardinal orientations, as compared to oblique states (Kenet et al. 2003), both in dwell times and in the strength of the spontaneous patterns. In our model, we could recreate this bias only if we biased both the area in the simulated cortical sheet allocated to representing the cardinal orientations, and the selectivity of the columns according to their orientation preference, making them more orientation selective the closer their preference was to cardinal orientations. In contrast, changing either of these factors on its own was not sufficient to produce the bias towards cardinal orientations. Although the empirical maps in our study did not provide conclusive evidence of such bias, the introduced factors are in line with a previous study in single cell electrophysiology that found more cells preferring cardinal orientations, as well as heightened orientation selectivity in cells with preference to cardinal orientations in cat early visual cortex (Li et al. 2003), but see (Altavini et al. 2016). Spontaneous maps and the functional architecture of primary visual cortex: As was predicted by previous theoretical work (Galán 2008),under our model the leading principal components(PCs) of the spontaneous activity, the orientation maps and the leading PCs of the connectivity matrix coupling the columns in the cortical sheet are virtually identical, and the loads of these PCs determine the probability of their emergence during spontaneous activity. This suggests that visual cortical feature maps should be thought of as the PCs of the underlying

31

patterns of connectivity, further implying that the features processed by the visual cortex are naturally ordered according to their saliency for upstream processing. Other types of functional maps: Apart for orientation maps, ocular dominance, direction, and spatial frequency maps are well known functional maps in cat early visual cortex (Hübener et al. 1997; Shoham et al. 1997). Therefore, according to the above reasoning, we would expect to find other such functional patterns expressed in the spontaneous PCs. However, in our study we only recorded orientation evoked activity, and hence lacked the necessary templates to directly verify the correspondence between our spontaneous data and other maps. In contrast, in an independent set of experiments in the primate, in which early visual cortex was imaged for both orientation and ocular dominance maps, we could determine that spontaneous data collected in visual cortex exhibited both spontaneous ocular dominance and orientation states (Grinvald et al. 2015). These results are expected from the well known strength of the various types of cortical maps in different species. The amplitude of evoked orientation maps in the adult cat is much larger than that of ocular dominance (OD) maps, direction selective maps, and spatial frequency maps (Shmuel and Grinvald 1996). In contrast, the opposite was found in monkeys (Shtoyerman et al. 2000; Slovin et al. 2002), for which OD maps are stronger than orientation maps. If indeed feature maps in visual cortex are PCs of the lateral connectivity matrix, and if as our results from both monkey and cat indicate, this is also true for parts of such maps, this would suggest that the wiring principles behind the emergence of the overlapping maps in early visual cortex are such that they ensure that the maps will be orthogonal throughout. This could be satisfied through maintaining invariant geometric relationships between the various maps. Therefore our results are in line with recent suggestions regarding the functional origins of

32

known feature maps in early visual cortex(Kremkow et al. 2016; Lee et al. 2016), and moreover indirectly support the notion of uniform feature coverage in early visual cortex (Hubel and Wiesel 1977 ; Swindale et al. 2000). Concluding remarks: Our results highlight the utility of the study of spontaneous activity under anesthesia: in this state, the impact of external perturbations is largely diminished and long range interaction strength appears to be dampened down, in effect highlighting lateral connectivity (Ferrarelli et al. 2010). Accordingly, local network structure and function can be studied in near independence. This enables to glean the computational principles resulting from the architecture of a given brain network, affording indispensable insight into the principles at the root of neuronal processing of information.

REFERENCES Altavini TS, Conde-Ocazionez SA, Eriksson D, Wunderle T, Schmidt KE. 2016. Selective interhemispheric circuits account for a cardinal bias in spontaneous activity within early visual areas. NeuroImage. Arieli A, Shoham D, Hildesheim R, Grinvald A. 1995. Coherent spatiotemporal patterns of ongoing activity revealed by real-time optical imaging coupled with single-unit recording in the cat visual cortex. Journal of Neurophysiology 73:2072-2093. Arieli A, Sterkin A, Grinvald A, Aertsen A. 1996. Dynamics of ongoing activity: explanation of the large variability in evoked cortical responses. Science 273:1868-1871. Benucci A, Frazor RA, Carandini M. 2007. Standing waves and traveling waves distinguish two circuits in visual cortex. Neuron 55:103-117. Berkes P, Orbán G, Lengyel M, Fiser J. 2011. Spontaneous cortical activity reveals hallmarks of an optimal internal model of the environment. Science 331:83-87. Blumenfeld B, Bibitchkov D, Tsodyks M. 2006. Neural network model of the primary visual cortex: From functional architecture to lateral connectivity and back. Journal of computational neuroscience 20:219-

33

241. Bonhoeffer T, Grinvald A. 1993. The layout of iso-orientation domains in area 18 of cat visual cortex: optical imaging reveals a pinwheel-like organization. The Journal of neuroscience 13:4157-4180. Calhoun V, Adali T, Pearlson G, Pekar J. 2001. A method for making group inferences from functional MRI data using independent component analysis. Human brain mapping 14:140-151. Chavane F, Sharon D, Jancke D, Marre O, Frégnac Y, Grinvald A. 2011. Lateral spread of orientation selectivity in V1 is controlled by intracortical cooperativity. Frontiers in Systems Neuroscience 5:4. Delucchi M, Garoutte B, Aird RB. 1962. The scalp as an electroencephalographic averager. Electroencephalography and Clinical Neurophysiology 14:191-196. Deneux T, Grinvald A. 2016. Milliseconds of Sensory Input Abruptly Modulate the Dynamics of Cortical States for Seconds. Cerebral Cortex. Engel AK, König P, Kreiter AK, Singer W. 1991. Interhemispheric synchronization of oscillatory neuronal responses in cat visual cortex. Science 252:1177-1179. Fekete T, Rubin D, Carlson JM, Mujica-Parodi LR. 2011. The NIRS Analysis Package: Noise Reduction and Statistical Inference. PloS one 6:e24322. Ferrarelli F, Massimini M, Sarasso S, Casali A, Riedner BA, Angelini G, Tononi G, Pearce RA. 2010. Breakdown in cortical effective connectivity during midazolam-induced loss of consciousness. Proceedings of the National Academy of Sciences 107:2681-2686. Fiser J, Chiu C, Weliky M. 2004. Small modulation of ongoing cortical dynamics by sensory input during natural vision. Nature 431:573-578. Fukushima M, Saunders RC, Leopold DA, Mishkin M, Averbeck BB. 2012. Spontaneous high-gamma band activity reflects functional organization of auditory cortex in the awake macaque. Neuron 74:899910. Galán RF. 2008. On how network architecture determines the dominant patterns of spontaneous neural activity. PloS one 3:e2148. Goldberg JA, Rokni U, Sompolinsky H. 2004. Patterns of ongoing activity and the functional architecture

34

of the primary visual cortex. Neuron 42:489-500. Gong P, van Leeuwen C. 2009. Distributed dynamical computation in neural circuits with propagating coherent activity patterns. PLoS computational biology 5:e1000611. Grinvald A, Omer D, Naaman S, Sharon D. 2015. Imaging the Dynamics of Mammalian Neocortical Population Activity In-Vivo. In. Membrane Potential Imaging in the Nervous System and Heart Springer p 243-271. Han F, Caporale N, Dan Y. 2008. Reverberation of recent visual experience in spontaneous cortical waves. Neuron 60:321-327. Hubel DH, Wiesel TN. 1977. Ferrier lecture: Functional architecture of macaque monkey visual cortex. Proceedings of the Royal Society of London B: Biological Sciences 198:1-59. Hübener M, Shoham D, Grinvald A, Bonhoeffer T. 1997. Spatial relationships among three columnar systems in cat area 17. The journal of neuroscience 17:9270-9284. Innocenti GM, Aggoun-Zouaoui D, Lehmann P. 1995. Cellular aspects of callosal connections and their development. Neuropsychologia 33:961-987. Jansen BH, Rit VG. 1995. Electroencephalogram and visual evoked potential generation in a mathematical model of coupled cortical columns. Biol Cybern 73:357-366. Ji D, Wilson MA. 2006. Coordinated memory replay in the visual cortex and hippocampus during sleep. Nature neuroscience 10:100-107. Kenet T, Bibitchkov D, Tsodyks M, Grinvald A, Arieli A. 2003. Spontaneously emerging cortical representations of visual attributes. Nature 425:954-956. Kremkow J, Jin J, Wang Y, Alonso JM. 2016. Principles underlying sensory map topography in primary visual cortex. Nature 533:52-57. Lee K-S, Huang X, Fitzpatrick D. 2016. Topology of ON and OFF inputs in visual cortex enables an invariant columnar architecture. Nature 533:90-94. Li B, Peterson MR, Freeman RD. 2003. Oblique effect: a neural basis in the visual cortex. Journal of neurophysiology 90:204-217.

35

Lu H, Zuo Y, Gu H, Waltz JA, Zhan W, Scholl CA, Rea W, Yang Y, Stein EA. 2007. Synchronized delta oscillations correlate with the resting-state functional MRI signal. Proceedings of the National Academy of Sciences 104:18265-18269. Luczak A, Barthó P, Harris KD. 2009. Spontaneous events outline the realm of possible sensory responses in neocortical populations. Neuron 62:413-425. Markounikau V, Igel C, Grinvald A, Jancke D, Sporns O. 2010. A dynamic neural field model of mesoscopic cortical activity captured with voltage-sensitive dye imaging. PLoS Comput Biol 6:e1000919. Mohajerani MH, Chan AW, Mohsenvand M, LeDue J, Liu R, McVea DA, Boyd JD, Wang YT, Reimers M, Murphy TH. 2013. Spontaneous cortical activity alternates between motifs defined by regional axonal projections. Nature neuroscience 16:1426-1435. Mohajerani MH, McVea DA, Fingas M, Murphy TH. 2010. Mirrored bilateral slow-wave cortical activity within local circuits revealed by fast bihemispheric voltage-sensitive dye imaging in anesthetized and awake mice. The Journal of Neuroscience 30:3745-3751. Nir Y, Mukamel R, Dinstein I, Privman E, Harel M, Fisch L, Gelbard-Sagiv H, Kipervasser S, Andelman F, Neufeld MY. 2008. Interhemispheric correlations of slow spontaneous neuronal fluctuations revealed in human sensory cortex. Nature neuroscience 11:1100-1108. Olavarria JF. 2001. Callosal connections correlate preferentially with ipsilateral cortical domains in cat areas 17 and 18, and with contralateral domains in the 17/18 transition zone. Journal of Comparative Neurology 433:441-457. Omer DB, Rom L, Ultchin U, Grinvald Aeditors. Year Published|. Title|, Conference Name|; Year of Conference Date|; Conference Location| Place Published|:Publisher|. Pages p|. Raichle ME, MacLeod AM, Snyder AZ, Powers WJ, Gusnard DA, Shulman GL. 2001. A default mode of brain function. Proceedings of the National Academy of Sciences 98:676-682. Rochefort N, Buzas P, Kisvarday Z, Eysel U, Milleret C. 2007. Layout of transcallosal activity in cat visual cortex revealed by optical imaging. Neuroimage 36:804-821. Romani S, Ahmadian Y, Grinvald A, Tsodyks M, Miller KDeditors. Year Published|. Title|, Conference Name|; Year of Conference Date|; Conference Location| Place Published|:Publisher|. Pages p|.

36

Shmuel A, Grinvald A. 1996. Functional organization for direction of motion and its relationship to orientation maps in cat area 18. The Journal of Neuroscience 16:6945-6964. Shoham D, Hübener M, Schulze S, Grinvald A, Bonhoeffer T. 1997. Spatio–temporal frequency domains and their relation to cytochrome oxidase staining in cat visual cortex. Nature 385:529-533. Shtoyerman E, Arieli A, Slovin H, Vanzetta I, Grinvald A. 2000. Long-term optical imaging and spectroscopy reveal mechanisms underlying the intrinsic signal and stability of cortical maps in V1 of behaving monkeys. The Journal of Neuroscience 20:8111-8121. Slovin H, Arieli A, Hildesheim R, Grinvald A. 2002. Long-term voltage-sensitive dye imaging reveals cortical dynamics in behaving monkeys. Journal of neurophysiology 88:3421-3438. Steriade M, McCormick DA, Sejnowski TJ. 1993. Thalamocortical oscillations in the sleeping and aroused brain. Swindale NV, Shoham D, Grinvald A, Bonhoeffer T, Hübener M. 2000. Visual cortex maps are optimized for uniform coverage. Nature neuroscience 3:822-826. Trong MN, Bojak I, Knösche TR. 2013. Associating spontaneous with evoked activity in a neural mass model of visual cortex. NeuroImage 66:80-87. Tsodyks M, Kenet T, Grinvald A, Arieli A. 1999. Linking spontaneous activity of single cortical neurons and the underlying functional architecture. Science 286:1943-1946. Tyszka JM, Kennedy DP, Adolphs R, Paul LK. 2011. Intact bilateral resting-state networks in the absence of the corpus callosum. The Journal of Neuroscience 31:15154-15162. Wendling F, Bellanger J-J, Bartolomei F, Chauvel P. 2000. Relevance of nonlinear lumped-parameter models in the analysis of depth-EEG epileptic signals. Biological cybernetics 83:367-378. Xie Y, Chan AW, McGirr A, Xue S, Xiao D, Zeng H, Murphy TH. 2016. Resolution of high-frequency mesoscale intracortical maps using the genetically encoded glutamate sensor iGluSnFR. The Journal of Neuroscience 36:1261-1272.

37

𝑓𝑚𝑎𝑥 = 50 maximal firing rate

𝑟 = 0.25 𝑚𝑉 −1 slope rate transfer function

𝑔𝑒𝑒 = 25 𝑚𝑉excit. to excit. synaptic gain factor

𝑔𝑖𝑒 = 10 𝑚𝑉inhib. to excit. synaptic gain factor

𝑔𝑒𝑖 = 5 𝑚𝑉excit. to inhib. synaptic gain factor

𝜏𝑒 = 8 𝑠𝑒𝑐 −1 excitatory pool time constant

𝜏𝑖 = 16𝑠𝑒𝑐 −1 inhibitory pool time constant

𝑕𝑒 = 𝑕𝑖 = −60 𝑚𝑉 resting membrane potential

𝜎 = 1.2𝑚𝑚 (12 𝑝𝑖𝑥𝑒𝑙𝑠) spatial connectivity extent

Table 1: model parameters and values

ACKNOWLEDGEMENTS We thank S. Kaufman for technical support.

38

Supplementary figures

Supplementary Figure 1: Interhemispheric synchrony of spontaneous oscillations in the    band. (A) the average VSD signal for each hemisphere, and both hemispheres taken together during a single trial.(B) The average VSD signal across the midline taken from both hemispheres during the same trial.

39

Supplementary figure 2: Spontaneously emerging orientation maps. Same as figure 1A, Examples taken from three of the animals used in this study.

40

Supplamentary figure 3: 6 different pairs of control patterns (T1,T2) were generated (see table). C1 is the control pattern used in Kenet et al. 2003, in which the maps were flipped horizontally (Hf). Vf denotes a vertical flip, and GT the ground truth (the corrected 0/90 and 45/90 maps). For spontaneous frames, the correlation to T1 and T2 was computed and the maximum of the absolute values was stored. Top: the differential maps were used to generate templates representing the expected patterns associated with the orientation domain (0 to 180 in steps of 1). The max correlation is shown for C1-C6 as well for the ground truth (GT) . Middle: table listing the exact composition of C1-C6. Bottom: the empirical distributions of max correlation pooled across animals. The GT distribution was significantly larger than all of the control distribution. As expected from A, the C4 and C6 control maps resulted in distributions that were significantly larger than the control reported in Fig. 1 due to the high resemblance of the templates to projected orientation states (e.g. C6 was highly correlated (r>.7) to the 65 map)

41

Supplementary figure 4: Correspondence between the principal components derived from spontaneous fluctuations and orientation maps. For each animal the top row depicts the functional maps obtained via a GLM applied to evoked data while the bottom row displays the principal components derived from spontaneous data which were the most similar to the corresponding maps.

42

Supplementary figure5: correlation of rate of orientation change (/frame - see figure 3) between hemispheres as a function of SI across all animals. Clearly expressed orientation states were associated with increasingly correlated rate of change in orientation between hemispheres.

43

Supplementary Figure 6: Detection of mosaic orientation states. (A) The anterior half of the 0° orientation map and posterior half of the 45° orientation map were combined to create an artificial frame displaying a mosaic orientation state (Blumenfeld et al. 2006). (B-D) Assignment of local angle and similarity index to every pixel of single optical frames is achieved by applying the same calculations used for entire frames (see text for details) in 10 pixel radius windows around each pixel. (B) the pixel-wise best matching template matched to (A): for each pixel the best fitting template for the surrounding 10 pixel neighborhood was computed, and the value for that pixel location in the resulting template was stored. This was repeated for all pixels.(C) Local angle calculation carried out on (A) for each pixel using the pixels within a 10 pixel radius from it (D) The image in (B) was compared to (A) by computing the local SI for each pixel. This was done by correlating corresponding circular regions of radius 10 in both images centered around each pixel in turn. As can be seen, mosaic states are expected to result in large patches of high SI (D) showing greater angle (circular) variance as compared to non-mosaic states (C).

44

Supplementary figure 7 - spontaneous maps are not replay: The pooled distribution across animals of frame by frame maximal correlation. The average max correlation (the largest value of the correlations to the 0,45 90 and 135 maps in each optical frame) before and after exposure to the stimuli was 0.1287 and 0.1234 respectively, the variances were 0.0058 and 0.0057 respectively, and the maximally correlated frame values were 0.62 0.61 respectively. A one sided Wilcoxon rank sum test indicated that the Median of the distribution after exposure was not larger (p=1).

45

Supplementary figure 8 - PCA denoising: (A-B) Single frame evoked by presenting a fullfield oriented moving horizontal grating before (A) and after (B) our PCA based denoising (see Methods). (C-D) Single frame with spontaneously occurring activity pattern before (C) and after (D) PCA denoising. (E) The difference C-D. (F) The GLM map corresponding to a 0 stimulus. Each frame was scaled to +/- 1.8 standard deviations to enhance contrast.

46