Introduction

Sensory environments are complex – rich with objects and events that specific brain circuits must decipher to ensure their accurate perception and to guide appropriate behavior. To efficiently perceive and navigate sensory environments, these brain circuits combine information from multiple sensory modalities, including vision, somatosensation, or hearing, a process known as multisensory integration. Several brain regions have been shown to generate multimodal representations, typically by aligning connections from the sensory periphery where specific sensory features are first encoded to downstream projection areas. In the rodent, whiskers transmit tactile sensations which are combined with visual information to guide many behaviors including goal-directed behavior [1], gap crossing and food source detection [2], navigation [3], prey-capture behavior [4], and recognizing object features and shapes [5]. Multisensory integration has been reported in subcortical areas such as the superior colliculus [6, 7] and the striatum [8]. Higher-order areas (HOAs) situated between primary sensory areas have also been found to play a critical role in multisensory integration [1, 3, 9]. HOAs receive excitatory feedforward input via long-range axonal projections [10] and inherit and match cortical maps from primary sensory areas [1, 3, 11]. Moreover, connectivity is organized in a topographic manner across the modalities, whereby sensory neurons converge onto higher-order cortex while preserving space [3, 12]. Specifically, visual and tactile stimuli have overlapping representations in a higher order association area, area RL, that lies between the primary visual (V1) and the primary somatosensory (S1) cortices. There is significant spatial alignment between V1 and S1 in area RL; e.g. the lower whisker pad aligns to the lower visual field of view and can elicit responses in the same neurons, called bimodal neurons, in RL [3]. Additionally, this multisensory integration allows for cross-modal generalization. For example, rodents can adapt to changing sensory environments and can generalize sensorimotor tasks to a different modality, i.e. visual, after learning a task from only tactile stimulus, due to aligned abstract representations of these modalities in RL [1]. Finally, evidence shows that rats combine and integrate visual and tactile signals in a supralinear manner, significantly enhancing sensory processing efficiency [5]. Although full functional maturation of visual HOAs occurs only after eye opening, basic properties such as responsiveness to stimuli and topographic organization between V1 (as well as S1) and RL are already present at eye opening [13, 14].

Before the onset of sensory experience, primary sensory cortices generate spontaneous activity, which plays a pivotal role in their proper maturation in terms of cellular properties and activity-dependent connectivity refinements [15–18]. Across the sensory areas, a general developmental profile of spontaneous activity can be observed: early global and dense spontaneous events transition to more localized and sparse patterns, typical of adult networks [15, 17, 19–21]. By promoting the refinement of appropriate connectivity [22–24], these spontaneous activity patterns prepare the developing cortex for future sensory input [25–27], but at the same time are affected by changing circuitry [28]. Spontaneous activity has been characterized and shown to guide the establishment of hierarchical architectures in HOAs in the visual system before eye opening [29], but much less is known about its role in the development and topographic co-alignment of multiple sensory modalities in multisensory HOAs.

Whereas activity-independent mechanisms are known to set up the first cues instructing axons where to go [30, 31], activity-dependent synaptic plasticity mechanisms instructed by either spontaneous activity or sensory experience, are crucial in refining this initial imprecise connectivity [22, 32, 33]. Indeed, various forms of Hebbian synaptic plasticity, which promotes the potentiation of synaptic strength as a result of correlated pre- and postsynaptic activity, have been proposed to guide the activity-dependent refinement of circuit connectivity in development [28, 34–38]. Some experimental work has also measured Hebbianlike learning rules specific to development [39]. Activity-dependent mechanisms likely also contribute to the establishment and alignment of maps from different sensory modalities to downstream circuits, which has been mainly investigated in the rodent superior colliculus and avian optic tectum [32]. For example, how visual experience aligns the auditory and visual maps during development has been mainly studied by manipulations that cause the two maps to misalign [40–44]. This suggests that these same principles might apply not only to realign and maintain maps after manipulation, but also to the activity-dependent emergence of aligned representations in developing circuits.

Here, we combine computational modeling and calcium imaging experiments to investigate the role of spontaneous activity in the emergence of aligned topographic maps from two primary sensory cortices (visual, V1 and somatosensory, S1) to RL, an associative cortical area rostro-lateral to V1, and the emergence of bimodal responses in RL. We find that spontaneous activity in V1, S1, and RL sparsifies during development, but S1 exhibits more mature activity patterns compared to V1, which desynchronize as early as postnatal day eight. Using these experimental insights, we expand a computational model of activity-dependent connectivity refinement between primary sensory areas (V1 and S1) and RL using biologically plausible Hebbian synaptic plasticity which was previously used to study topographic map-formation and receptive field refinement between the thalamus and V1 [28]. Our model predicts that aligned topographic connectivity refinement and the emergence of bimodal neurons requires an intermediate amount of correlation between the primary sensory cortices (V1 and S1), which we investigate in our data. This connectivity refinement through spontaneous activity between the primary sensory areas and RL is sufficient to explain the developmental sparsification in RL whereby RL neurons become responsive to either or both primary sensory cortices in a ratio that ensures optimal information processing.

Results

Spontaneous activity in V1, S1 and RL during development

To investigate how spontaneous activity matures differentially across the mouse’s developing primary visual (V1), somatosensory (S1), and higher-order areas, specifically area RL, we performed two-photon imaging with one field of view per developing cortex in unanesthetized mouse pups at postnatal days (PN) 8-16 (Figure 1A,B). Spontaneous activity was observed across all three cortical areas (Figure 1C). We used a linear mixed model analysis [45] to estimate the differential effect of age on different aspects of spontaneous activity. As described in previous experiments [17, 21], we found that spontaneous activity in V1, sparsifies throughout development, i.e. the average amplitude, duration, and participation rate of spontaneous activity decrease with age, as reflected in the significant negative effect of age on these properties (Figure 1D-F, H, see Methods). In contrast, the frequency of events increases with age, as reflected in the significant positive relationship between age and event rate (Figure 1G, H). Similar trends of activity sparsification were also observed in S1 and RL. There were no differences between participation rate across areas, with all areas following a similar decreasing trend over time (Figure 1F,H). Furthermore, comparing area RL to V1 also showed no significant differences in duration of events and event rate between these two areas (Figure 1E,G,H), although amplitudes were slightly higher in RL compared with V1. However, we observed some significant differences, especially in S1. Specifically, spontaneous events in S1 have significantly lower amplitudes, lower durations, and higher rates compared to the other regions, at younger ages (Figure 1D,E,G,H). Together, this suggests that spontaneous activity in S1 has more mature amplitudes, duration, and event rates earlier in development, whereas RL and V1 develop more similarly during the second postnatal week. RL was found to slightly lag behind V1 and S1.

Figure 1.

In summary, spontaneous activity in developing primary cortices and RL sparsifies with different temporal profiles, where the somatosensory cortex leads before the visual cortex and RL.

The primary somatosensory and visual cortex exhibit correlated spontaneous activity

Because RL needs to integrate activity from both V1 and S1 simultaneously, we hypothesized that there could be correlated activity between the primary sensory cortices to guide the emergence of this multisensory integration via some activity-dependent mechanism. To this end, we reanalyzed our published wide-field calcium imaging recordings from both V1 and S1 [46] to assess the temporal relationship between spontaneously occurring activity in these two primary sensory cortices (Figure 2A). We found that activity between the visual and the somatosensory cortex is often, but not always, temporally synchronized (Figure 2B). Indeed, analyzing the average activity between V1 and S1 of several animals at postnatal days PN9-PN12 revealed a range of Pearson correlation coefficients, with a mean of ∼0.5 (Figure 2C).

Figure 2.

Previous work has shown that connectivity between V1 and RL as well as S1 and RL is topographically organized in mature animals to preserve space [3, 11]. Furthermore, this organization is also aligned across the modalities, i.e. stimuli in the same region of sensory space tend to converge onto the same cells in RL [3, 12, 29]. Hence, in addition to temporal correlations, we evaluated whether spontaneous activity also carries information about spatial topography between the two primary sensory cortices and RL. We computed functional correlation maps for three distinct locations in V1 and RL, as well as three other locations in S1 and RL, showing that spontaneous activity in RL is coordinated topographically relative to both V1 and S1 (Figure 2D,E). Thus, spontaneous activity between V1 and RL and between S1 and RL is also spatially correlated and carries topographic information. This suggests that some activity-dependent mechanism might be using these temporal and spatial correlations during development to establish topographic connectivity maps between the primary sensory areas and RL, and hence enable multisensory integration in the adult animal.

Correlated spontaneous activity can generate topographic connectivity between the primary sensory cortex and RL and bimodal cells in RL

To test the hypothesis that an activity-dependent mechanism may use the spatio-temporally structured spontaneous activity in the primary sensory cortices to set up topographic connectivity from each area into RL and generate bimodal cells, we built and analyzed a computational model with three populations of neurons for V1, S1, and RL (Figure 3A) [28]. We modeled the activity-dependent synaptic plasticity of feedforward connections between V1 and RL, and S1 and RL, with each area subject to spontaneous activity as characterized before (Figure 1). While we initialized the connectivity as random, we included a slight bias for neighboring cells in the primary sensory cortices to project onto neighboring cells in RL (Figure 3B). This bias incorporates activity-independent effects like molecular gradients that are known to set up a coarse topography before the onset of activity-dependent plasticity, at least for connectivity between the sensory periphery and subcortical or primary cortical areas [38, 47–49]. Activity-dependent plasticity was implemented through a synaptic plasticity rule involving Hebbian and heterosynaptic terms. The Hebbian term potentiates synaptic weights based on coincident pre- and postsynaptic activity. The additional heterosynaptic term was implemented to depress synaptic weights, but only in the presence of presynaptic activity (see Methods)

Figure 3.

This rule was previously derived from interacting molecular neurotrophins [39, 50, 51] and shown to produce selectivity to correlated inputs from random initial connectivity [37, 52]. We then implemented spatiotemporally correlated spontaneous activity in areas V1 and S1 as measured experimentally (Figure 2C) which this synaptic plasticity rule used to modify the synaptic weights in the network. Specifically, to implement temporal correlations, S1 and V1 neurons were driven by a combination of events, independent and shared in each V1 and S1. To implement spatial correlations within a given area (V1 or S1), each independent event activated a randomly chosen, contiguous set of neurons, which we have previously shown can establish topographic connectivity between two areas [28]. To implement spatial correlations between V1 and S1, the same set (in terms of topographic identity) of neurons in both S1 and V1 were activated. The common inputs giving rise to shared events might reflect shared inputs from the subcortical areas such as the thalamus or superior colliculus into the primary sensory cortices [53–55] or correlation induced by direct connections between V1 and S1 [56, 57].

Using this spatiotemporally correlated spontaneous activity in the primary sensory cortices, the Hebbian plasticity rule (Eq. 1) refined and stabilized connectivity between V1/S1 and RL (Figure 3D), as shown before in a two-layer network [28]. We investigated the emergent connectivity properties for a range of V1 and S1 temporal correlations (Figure 3E). To quantify the topographic organization of connectivity between V1 and RL, and S1 and RL, and the alignment of maps from V1/S1 to RL, we defined two mathematical measures, topography and alignment (Figure 3F, see Methods). Weak temporal correlations (where independent events dominate over shared events) resulted in strong competition between the primary sensory cortices with dominant synaptic weight depression. The strong competition generated sparse connectivity matrices to RL, with many RL neurons decoupling from the primary sensory cortices (Figure 3E, left). This resulted in poor topography and map alignment of the maps between V1/S1 and RL (Figure 3G,H). RL had almost no bimodal cells, i.e. neurons that receive input from both primary sensory cortices, in contrast to what is found in mature RL experimentally [3] (Figure 3I). Increasing temporal correlations between the primary sensory cortices increased alignment of the V1/S1 to RL maps (Figure 3H). However, very strong correlations were detrimental to V1/RL and S1/RL topography (Figure 3E,G). Strong temporal correlation between V1 and S1 rendered many RL neurons selective to visual and somatosensory inputs at the same topographic location as reflected in the identical connectivity matrices (Figure 3E, right), overriding the initial topographic bias (Figure 3B). Only a medium-strength temporal correlation between V1 and S1, as found experimentally (Figure 2C), generated a connectivity profile which prevented RL cells to decouple from V1 and S1, while preserving topography (Figure 3E, middle).

We quantified the percentage of bimodal neurons in our model as a function of the temporal correlation between V1 and S1 (Figure 3I). By performing a steady-state analysis of our model (see Methods), we derived the critical amount of temporal correlation, c∗, between S1 and V1 required for the emergence of bimodal neurons. This parameter depends solely on two parameters of the model, i.e. the average firing rate of cells in the sensory cortices, ⟨pre⟩, and the heterosynaptic offset from Eq. 1,

With these connectivity refinements, the generated activity in RL became sparser in terms of amplitude and participation rate (Figure 3J), consistent with the sparsification observed in the experimental data (Figure 1C,D,F).

In summary, our proposed computational model can account for the developmental sparsification of spontaneous activity in higher-order cortical area RL through the refinement of feedforward connectivity between two primary sensory cortices, V1 and S1, and RL. Furthermore, the model predicts that only medium-strength correlations between primary cortices produce topographically organized and aligned connectivity in V1/RL and S1/RL, and a mixture of bimodal and unimodal neurons in RL consistent with experimental data.

More mature spontaneous events in S1 can guide map alignment between V1 and RL

Having established that the spatio-temporal correlations between V1 and S1 can set up topographic and aligned connectivity to RL as well as generate bimodal cells in RL, we next turned to investigate a possible role of the differentially more mature spontaneous activity in S1 compared to V1 (Figure 1H). To examine this, we systematically increased the bias of the topographic connectivity between S1 and RL (Figure 4A). While unsurprisingly a stronger bias improved the topographic connectivity between S1 and RL (Figure 4B, right), it also improved the topographic organization of connectivity between V1 and RL (Figure 4B, left, and C). Therefore, spontaneous activity in S1, which matures more rapidly than V1, can guide the refinement of connectivity between V1 and RL.

Figure 4.

A mixture of bimodal and unimodal neurons in RL performs optimal decoding of activity

RL is involved in integrating information from different sensory modalities to obtain a more accurate representation of statistically rich sensory environments. In particular, an area that integrates visual and tactile stimuli is able to more accurately infer such information; e.g. when rats use both tactile and visual information, their performance in a task involving determining object orientation is better than the sum of the performance of each modality individually [5]. To investigate how the topographically aligned cortical connectivity between V1 and RL, and S1 and RL, might facilitate the integration of information from different sensory modalities, we quantified how much information about the activity of each primary sensory cortex can be decoded from the activity in RL. We quantified the amount of decoded information as the coefficient of determination, R², which can be interpreted as the fraction of variance in V1 or S1 that can be reconstructed from activity in RL. To do this, we fitted a linear regression model using RL activity and both V1 and S1 activity, and calculated R² from the error between our fitted linear model and each sensory area individually, which we refer to as the ‘V1 or S1 reconstruction’ (Figure 5A). [58, 59] (See Methods). R² is naturally affected by how much activity in V1 or S1 is propagated into RL. In our model (Figure 3) and in biology [11, 60], some regions of the primary sensory cortices are not connected to RL; for example, RL contains only a partial representation of the visual field, hence V1 cells outside of that region do not project to RL [10,61]. Consistently, we found that the fraction of variance reconstructed, R², for either V1 or S1 independently increased with the percentage of connected neurons from the respective cortex into RL (Figure 5B, see Methods). We showed analytically that a locally connected network, where every primary sensory cortical neuron makes at least one connection to RL and every RL neuron receives at least one connection from the primary sensory cortices, can maximize the total fraction of variance reconstructed for both V1 and S1 if the activity in V1 or S1 is perfectly temporally correlated (correlation = 1) (see Methods). These conditions translate to the scenario where all neurons in RL are bimodal. However, when activity across V1 and S1 is moderately correlated, having some unimodal RL neurons can achieve a higher total maximum fraction of variance for both V1 and S1 compared to the purely bimodal case (Figure 5C, see Methods). Having unimodal neurons increases variance explained when activity is not perfectly correlated because incomplete information from one modality can be completed with information from the other modality. Furthermore, the activity of unimodal neurons provides unambiguous information about whether only one or both of the primary sensory cortices is active, whereas in the case with all bimodal neurons this is indistinguishable. This suggests that the mixture of unimodal and bimodal neurons in RL [3] is the optimal solution for encoding the maximum amount of information from the activity in V1 and S1 [62].

Figure 5.

Discussion

We investigated the role of spontaneous activity in the maturation of connectivity from two primary sensory cortices, V1 and S1, to RL, a higher-order association area, by combining computational modeling with experiments. We found that while spontaneous activity sparsifies during development in all investigated areas, S1 exhibits relatively mature, desynchronized activity patterns as early as postnatal day eight. Using these experimental insights, we constructed a computational model of activity-dependent connectivity refinement from each V1 and S1 to RL. To generate topographic connectivity from V1 and S1 to RL, aligned connectivity maps for the two sensory projections and a mixture of unimodal and bimodal neurons in RL, our model predicts that spontaneous activity between V1 and S1 should exhibit moderate temporal correlations, which we also found experimentally. Furthermore, we demonstrated how a more mature somatosensory cortex can guide the refinement of visual cortex connectivity to RL, suggesting that inter-cortical influence plays a central role in the establishment of refined and topographically aligned connectivity during development. Finally, we found that moderate correlations between V1 and S1 result in a mixture of bimodal and unimodal neurons in RL which can optimize information processing in RL when either one modality or both are active. Thus, by combining modeling with experiments, we found that activity-dependent plasticity can instruct the formation of refined and topographically aligned connectivity between primary sensory and high-order cortices in development to produce efficient information integration from multiple modalities.

Origin of correlated events in the primary sensory cortices

We find that early spontaneous activity with medium correlations can lead to the optimal amount of bimodal neurons and alignment of these sensory maps before sensory experience, which is crucial for effective multisensory integration in higherorder cortical areas. Such medium correlations indeed exist in the experimental data (Figure 2C). This result makes the prediction that alterations in the correlations in early spontaneous activity might affect later multisensory integration, which can be empirically tested. For instance, spontaneous activity in V1 has been shown to have different correlation properties in a mouse model of Fragile X [63] and in cases of disrupted inhibition [19, 46]; our model therefore suggests that multisensory integration in RL may be perturbed in these cases.

What might be the origin of the correlated activity we observed in the visual and somatosensory cortex (Figure 2)? One possibility is that the source is external to the animal [64], i.e. a single stimulus with both visual and somatosensory aspects. Indeed, while a substantial portion of activity in the developing cortex occurs independently of sensory stimulation [23], both the immature visual [65] and somatosensory [66] cortices exhibit some degree of sensory-driven activity. Although it is unlikely that this external stimuli gives any topographic information before eye opening, it could help to correlate the activity between sensory cortices. Alternatively, subcortical areas like the (primary or higher-order) thalamus are in an ideal position for coordinating activation across distinct regions of the cortex [29, 54, 67]. The developing thalamus exhibits spontaneous propagating waves, which have been shown to guide the formation of sensory maps to the cortex and support circuit plasticity [54,68,69]. Spontaneous calcium waves propagate among the different sensory-modality thalamic nuclei up to the cortex and hence provide a means of communication among sensory systems [69]. Finally, correlated events in V1 and S1 could also originate from lateral long-range connections between the visual and somatosensory cortex. These types of connections emerge during the first two postnatal weeks in mice [70] and potentially directly propagate activation from one area to another [56, 57]. More specifically, it has been shown that postsynaptic excitatory neurons in a subregion near the anterior border of V1, situated in close proximity to S1, receive input from local fast spiking inhibitory neurons that have direct cortico-cortical input from the posterior area of S1 [71]. This cross modal suppression may contribute to these medium level correlations during development.

Precise timing of map formation

We have shown that the somatosensory cortex’s early maturation is able to drive the aligned topographic development of the visual cortex Figure 4. This agrees with previous work showing earlier development of somatosensory cortex. For example, layer 4 (L4) to L2/3 excitatory synapses in somatosensory cortex were shown to mature 2 weeks earlier than those in visual cortex [72]. A similar template for the maturation of one map provided by another is observed in the refinement of cortico-collicular connections. Projections from S1 to the superior colliculus arrive as early as postnatal day 2 while projections from V1 enter the colliculus around postnatal day 6 [73], giving more time for S1 maps to refine. The precise timing of connectivity refinement between S1 and the superior colliculus likely guides the later developmental connectivity refinement between V1 and the colliculus and helps align the two maps [73, 74]. Moreover, this process of visual and somatosensory map alignment in the superior colliculus seems to be activity-dependent, with spontaneous activity in the retina playing an important role by driving the alignment of cortico-collicular connections from V1 to a genetically-induced duplicated retino-collicular map [74]. This modality difference in map refinement with the somatosensory leading the visual map seems to be present also in the refinement of connections between S1 vs. V1 and RL as we have shown earlier. Although the activity-dependent plasticity mechanisms behind this are experimentally unknown, inspired by earlier works [28, 50], our models suggest that similar activity-dependent plasticity involving Hebbian and heterosynaptic terms are sufficient.

Role of recurrent connectivity

Neural dynamics in the adult cortex are strongly affected by recurrent connections [75, 76]. In our model, we have chosen to focus on feedforward connections between the primary sensory and RL for several reasons. First, recent experimental evidence from the developing somatosensory cortex suggests that feedforward connectivity might be established first [77], and second, the strength and functional specificity of recurrent connections in the visual cortex increases only later in development [78]. Furthermore, abstraction is important for deriving robust biological insights from models [79] and here we are focused on how these feedforward connections arise and refine themselves over development building on our previous work of refinements between the thalamus and V1 [28].

Uni- and bimodal neurons in other species

Neurons responding to only a single or to multiple modalities are common in other animals apart from mice, like macaques [80–82], cats [83–85], or ferrets [86, 87]. The surprising presence of unimodal neurons in association cortices has inspired a number of studies attempting to explain their functional significance [85, 88, 89]. Close examination of the receptive fields of neurons in the associative cortex of the cat revealed that stimulus location is likely represented in a distributed population code and that unimodal neurons could contribute by improving the accuracy of the encoded stimulus location [85]. Furthermore, the response patterns of unimodal neurons in the cat extrastriate cortex revealed that while these neurons do not respond to their non-preferred modality, their activation is nonetheless modulated when non-preferred and preferred stimulus occur at the same time [88], further suggesting that unimodal neurons contribute to a population code. Varying the proportion of uni- and bimodal neurons might thus contribute to a continuum of multisensory responses across different cortical areas [90]. Thus, our proposed explanation for the mixture of uni- and bimodal neurons in terms of optimal decoding of input origin integrates well into the existing literature.

Methods

Experimental procedure Animals

All experimental procedures were approved by the institutional animal care and use committee of the Royal Netherlands Academy of Sciences. C57BL/6J mice of both sexes were used. All animals were aged between postnatal days (PN) 8-16. Mice of this strain open their eyes at P14.

In utero electroporation

For wide-field imaging (Figure 2) pyramidal neurons in layer 2/3 of the visual cortex were transfected with GCaMP6s (2 mg/ml) and DsRed (2 mg/ml) at embryonic day (E) 16.5 using in utero electroporation [91]. Pregnant mice were anesthetized with isoflurane and a small incision (1.5–2 cm) was made in the abdominal wall. The uterine horns were carefully removed from the abdomen, and DNA was injected into the lateral ventricle of embryos using a sharp glass electrode. Voltage pulses (five square wave pulses, 30 V, 50-ms duration, 950-ms interval, custom-built electroporator) were delivered across the brain with tweezer electrodes covered in conductive gel. Embryos were rinsed with warm saline solution and returned to the abdomen, after which time the muscle and skin were sutured.

Surgery

Animals were anesthetized with isoflurane (3 percent in 1 l/min O₂). After anesthesia had become effective, lidocaine was used for local analgesia and a head bar with an opening (Ø 4 mm) above the visual cortex (0.5–2.5 mm rostral from lambda and 1–3 mm lateral from the midline) was attached to the skull with superglue and dental cement. For calcium imaging, a small craniotomy above the visual cortex was performed. The exposed cortical surface was kept moist with cortex buffer (125 mM NaCl, 5 mM KCl, 10 mM glucose, 10 mM HEPES, 2 mM MgSO₄ and 2 mM CaCl₂ [pH 7.4]).

Calcium indicator application

For the imaging experiments shown in Figure 1, the calcium-sensitive dye Cal-590 (Aat Bioquest) was dissolved in 4 µl pluronic F-127, 20 percent solution in DMSO (Invitrogen) and further diluted (1:10) in dye buffer (150 mM NaCl, 2.5 mM KCl, and 10 mM HEPES) to yield a final concentration of 0.5 mM. The dye was then pressure-ejected into S1, RL and V1 at 10–12 psi for 8–13 min with a micropipette (3–5 MΩ) attached to a picospritzer (Toohey).

Image acquisition

One hour before imaging, anesthesia was stopped, and animals were recorded in the absence of anesthesia. In vivo calcium imaging with cellular resolution (Figure 1) was performed on a Nikon multi-photon microscope (A1R-MP) with a 0.8/16x water-immersion objective and a Ti:Sapphire laser (Chameleon II, Coherent) Nikon software. We recorded the movement signal of the scan mirrors to synchronize calcium imaging and electrophysiology. Pixel size was 650 nm and images of 330 by 330 µm were recorded at 5-10 Hz.

In utero electroporated pups were used for wide-field calcium imaging of visual and somatosensory cortex (Figure 2). Calcium events were recorded with a Movable Objective Microscope (MOM, Sutter Instrument). Time-lapse recordings were acquired with a 4x objective (0.8 NA, Olympus) and blue light excitation from a Xenon Arc lamp (Lambda LS, Sutter Instrument Company). A CCD camera (Evolution QEi, QImaging) was controlled by custom-made LabVIEW (National Instruments) based software and images were acquired at a frame rate of 20 Hz.

Image processing

Images were processed as described previously [46]. For two-photon image processing, drift and movement artifacts were removed from each recording using NoRMCorre [92]. Each recording was aligned to the first recording in the series to remove any movements between recording sessions. Delta F stacks were made using the mean fluorescence per pixel as baseline, F₀. ROIs were hand-drawn using ImageJ (NIH). Automated transient detection and further data processing was performed using custom-made MATLAB software (MathWorks).

Epifluorescence image processing: Delta F stacks were made using the average fluorescence per pixel as baseline. V1 was identified based on activity coordinates and shape after this method of identification was confirmed through immunohistochemistry for vGluT2 [19].

Statistical analysis

In two-photon recordings (Figure 1), activations within each cell were detected as in [93] as a peak with both an absolute height and relative prominence of at least 5% ΔF/F₀. When more than 20% of recorded cells were active simultaneously, the activation was considered an event. In wide-field calcium imaging (we re-analyzed data from [46] (Figure 1)), events were detected as calcium peaks with an amplitude above 2x the standard deviation.

Functional correlation maps

Functional correlation maps were generated as in [29] as follows: Three seed areas were chosen (10x10 pixels each) within V1 or S1 and their time courses across a recording were determined by averaging the pixel values in each seed area for each frame of the recording. Next, the Pearson coefficient of correlation between the time courses of each seed area and those of all pixels in the field of view were determined and represented in one image for each seed area. Finally, the three images of each seed area were combined into one RGB image, such that each channel represents one functional correlation map. The actual color of each pixel represents the highest coefficient of correlation value across the three channels. The final RGB map reflects the topographic organization of the recorded spontaneous activity patterns.

Linear mixed model

To correctly account for the nested data structure commonly encountered in neuroscience research [45], we fit a linear mixed-effects model that treats animal and sample identity as random factors. We included interaction effects between age and region within our fixed effects. No transformations of the data were performed. Estimated effect sizes are shown in (Figure 1). Effect sizes on cortical region affect the intercept, but age and age-region interaction estimates affect the slope of the regression line. One-hot encoding was used because the regions are categorical variables. This means that one region, i.e. V1, is chosen as baseline and every other region is compared to this baseline. Any differences due to region are added/subtracted from the intercept for V1. Age and age-region interaction effects change the slope and should therefore be added or subtracted from the estimated slope for region V1. Since PN8 was the earliest age used in our data, 8 was subtracted from all ages so that PN8 is the estimated intercept. Our linear model may not be valid outside of the ages used (PN8-PN16) and no extrapolation should be done without experimental confirmation. This multilevel analysis allowed us to examine the effects of age and cortical region on four response variables: amplitude, event rate, participation rate, and duration. To obtain p-values, we used the Satterthwaite’s degrees of freedom method, as implemented in the lmerTest package [94] in R. The linear mixed model formula is:

where y is one of the response variables, i.e. amplitude, event rate, participation rate, or duration and age is postnatal day. β₀ is the intercept for V1 at age PN8, and β₁ and β₂ are the differences in intercept for cortical region (S1 and RL are dummy variables set to 0 if not predicting that region or 1 if predicting that region). β₃ is the slope for V1 over age and β₄ and β₅ are the differences in slope (age-region interaction effects) from V1 over age for the cortical region of interest. ɛ_ij is the error from both random effects (mouse ID and sample ID) and residual error.

Rate-based network model

Extending the model introduced by [28], we studied a feedforward, rate-based network with two onedimensional input layers, v and s, and one output layer r with N ^v, N ^s, and N ^r neurons each. We used periodic boundary conditions in each layer to avoid edge effects. The initial connectivity between v and r, W_v, and s and r, W_s, is initially all-to-all with uniformly distributed weights in the range w_ini ∈ [a, b]. In addition, a topographic bias was introduced by modifying the initially random connectivity matrix to have the strongest connections between neurons at the matched topographic location. These biased weights decay with a Gaussian profile with increasing distance (Figure 3)), with amplitude a and spread s. During the evolution of the weights, soft bounds were applied on the interval [0, W_max].

The input populations v and s receive spontaneous activation from a homogeneous Poisson process with an inter-event interval of 150s, implemented using a Binomial random variable with sufficiently small time steps. A spontaneous event increased the firing rate of a fraction of neighboring neurons to 1 for L_dur seconds. The activity of the output layer r is given as a leaky integrator of the input populations,

We studied the weight evolution under activity-dependent learning rules. In particular, following [28] we chose a rule that combines Hebbian and heterosynaptic effects that can produce selectivity,

and

where θ is a heterosynaptic offset constant that defines a distance-dependent heterosynaptic depression. Heterosynaptic plasticity occurs when a nearby unstimulated synapse stabilizes or depresses due to the plasticity in a stimulated one. This helps to stabilize the Hebbian plasticity [50]. Table 1 lists all parameters. For simulations with varying correlation values, the % of correlated events is chose from a uniform distribution between 0 and 1, U[0,1]. For simulations with varying S1 bias, a Gaussian spread of 4 and a random bias from a uniform distribution between 0.05 to 0.6 was chosen. Simulation time of 500,000 ms was used.

Table 1.

Expected weight dynamics and bimodality threshold

Since the dynamics of the output layer are much faster than the change in weight, r << w, we made a steady state assumption and computed

This allowed us to write the expected change in weight between (e.g.) the v and r populations as a function of the input statistics of the v and s population,

This formulation allowed us to investigate the critical required amount of correlation between cells in v and s that allows for the formation of stable bimodal cells in r. In particular, we investigated the expected change of weight in the configuration where all connections between v and r are fully depressed and only one connection between s and r is potentiated:

and

The expected weight change in this configuration is

where Cov is the covariance and µ_v, µ_s are the average firing rates. The critical amount of correlation, where potentiation and depression balance each other out, is characterized by the condition

which we solved for the correlation, γ_sv = Cov(s_i, v_i)/(σ_vσ_s), as

Assuming that rates and variances in v and s are equal and that the activity is Poisson, µ = σ², the condition for the critical correlation further simplifies to

Optimal weight matrix for correlated input populations

We considered a concatenation of the two “input” populations, X = [v; s] of size M × 1, that drive activation in the “output” population, r = ϕ(WX), where ϕ is a (possibly non-) linear transfer function, r is size N ^r × 1 and M = N ^v + N ^s. The decoding problem introduced in Figure 5 can be phrased as inferring activity of the input populations from the output population, = Aϕ(WX), where A is the pseudoinverse of W, A = (W^T W)⁻¹W^T. We then searched for the optimal weight matrix, W_opt = [W^v, W^s] of size N ^r × M, where W^v represents the weights from v to r, and W^s represents the weights from s to r, that minimizes the squared error between X and , E = (X − )² and therefore maximizes R².

As the dimension of r is substantially smaller than the dimension of X, dim(r) ≪ dim(X), the optimal weight matrix needs to preserve the largest possible amount of variability of the input while projecting it into a lower dimensional space. Using insights from principal component analysis theory [95], we can demonstrate that the optimal weight matrix is composed of eigenvectors of the covariance matrix of X with the largest eigenvalues,

To better understand the structure of the eigenvectors, we can exploit the fact that in our set-up the covariance matrix, C is circular, i.e. every row of X^T X can be obtained by a cyclic shift of the generating vector, c = [c₀, c₁, …, c_M−1]:

Eigenvectors and eigenvalues of circulant matrices have well-known closed form solutions [96]. The eigenvectors are the Fourier modes:

where ω = e^2πi/M, a primitive M -th root of unity. The corresponding eigenvalues can be computed as

Since the covariance matrix is symmetric, i.e. c_j = c_M−j for all j and all values are real in the covariance matrix we can rewrite the eigenvalues in terms of their real parts:

For a real-symmetric circulant matrix, the real and imaginary parts of the eigenvectors are themselves eigenvectors:

The first M/2 eigenvectors with nonzero eigenvalues from v and w will comprise W_opt which maximizes R². Only the first M/2 are used since the dimension of r is M/2 in our case where N ^v = N ^s and can have a maximum of M/2 eigenvectors.

Our experimental set-up with correlated input populations translates into a generating vector c that has nonzero terms around the diagonal and offset from the diagonal by M/2. If the inputs between V1 and S1 are fully correlated, then c_k = 1 for k ∈ {0, . . . , W/2, (M − 1)/2 − W/2, . . . , (M − 1)/2 + W/2, (M − 1) − W/2, . . . , M − 1} and c_k = 0 else, where W is the extent of spatial correlation of the neurons. Taking this into consideration and removing any terms that simplify to zero due to the periodicity of the cosine function gives us the eigenvalues:

Hence, the eigenvectors with nonzero eigenvalue of the covariance matrix are exactly the first half of the eigenvectors of the sin and cos functions with an even-numbered period. Due to this simplification of the odd j’s leading to λ_j = 0 and even numbered periods, this results in eigenvectors that are identical on the indices {0, . . . , M/2 − 1} and {M/2, . . . , M } which is the definition of a bimodal cell and therefore all cells must be bimodal to maximize the R² when input correlation = 1 between V1 and S1.

Topography

Network topography measures the degree to which the initial topographical bias is preserved in the final receptive field. To calculate topography, we calculated the Pearson’s correlation coefficient, r, of the final receptive field with an ideal receptive field template (Figure 3). The ideal receptive field template that was used has weights of 1 with a Gaussian spread of 4 from center.

Alignment

Alignment is a measure of similarity between the final receptive fields of V1 and S1. To calculate alignment, we calculated the Pearson’s correlation coefficient, r, between the final receptive field of S1 and the final receptive field of V1.

Reconstruction

To calculate the reconstruction of V1 and S1 activity from RL in Figure 5, we fitted a linear model between all 3 layers’ activity patterns for 1000 sec using the final weight matrix of RL (independent variable) and V1 and S1 combined (dependent variable). Duplicates of RL activity pattern were removed from all layers activity and 80% of the data was used to avoid over-fitting. Then the coefficient of determination was calculated as

where SS_res is the sum of the residuals squared and SS_tot is the total sum of squares, using only one sensory area’s activity at a time, e.g. RL vs. V1 or RL vs. S1 activity compared with the linear fit of all 3 layers.

Acknowledgements

We thank Mark Hü bener for feedback on the manuscript and all members of the Gjorgjieva and Lohmann labs for comments and discussions. This work was supported by the Max Planck Society (J.H.K., M.E.W., D.K. and J.G.), the Technical University of Munich (J.M.D. and J.G.), the European Research Council (ERC) under the European Union’s Horizon 2020 research and innovation program (Grant agreement No. 804824 to J.G.) and Nederlandse Organisatie voor Wetenschappelijk Onderzoek (grant numbers ALWOP.216, OCENW.KLEIN.535, 865.12.001, OCENW.M.22.310) to C.L., ZonMw (9126021) to C.L. and Stichting Vrienden van het Herseninstituut (805254845) to C.L.

Author contributions

Conceptualization: J.G. and C.L. Experiments: N.Z., A.R., G.J.H. and P.P.M. Modeling and analysis: J.M.D., J.H.K., M.E.W., D.K. and J.G.; Writing: J.M.D., J.H.K., P.P.M., C.L. and J.G.

Data and code availability

All code for generating the model and figures can be found at: https://github.com/comp-neural-circuits/MultisensoryIntegration.git