Subject sex and partner sex modulate social touch responses across multiple cortical areas

Touch is a fundamental aspect of mammalian social, parental and sexual behavior. Human affective touch is critical for healthy child development and shows great promise as a novel therapeutic strategy for mental disorders characterized by social dysfunction, such as anxiety, depression and autism spectrum disorder. However, despite our detailed knowledge about cortical processing of non-social touch, we still know very little about how social touch modulates cortical circuits. We investigated the activity patterns of single neurons (N = 1156) across five sensory and frontal cortical areas in both male and female rats (N = 28) engaging in naturalistic social facial touch with male and female conspecifics. We found that information about social touch is widely available across cortex. Besides touch, the sex of the interaction partner (a biologically significant feature) is a major determinant of single neuron activity, and across cortex we observed 25.7% ‘touch’ and 11.9% ‘sex-touch’ responses. Although all areas investigated had access to social touch and partner sex information, social touch modulated different cortical areas in different ways. Finally, we found that network activity patterns during social touch depend on both subject sex and partner sex. Interestingly, these sex-differences in network activity patterns were differences in response magnitude and would not be evident without single cell resolution. Our observations suggest that socio-sexual characteristics of touch (subject and partner sex) widely modulate cortical activity and need to be investigated with cellular resolution.


Introduction
Social touch is a powerful emotional stimulus [1][2][3][4] . Harnessing the ability of social touch to modulate emotion for example by caress or massage has emerged as a protective and therapeutic strategy for various mental health conditions, such as anxiety, depression and autism spectrum disorder [5][6][7] . However, despite our detailed knowledge about cortical processing of discriminative, non-social touch, we still know very little about how social touch modulates cortical circuits.
Human social touch is different from non-social touch, in part due to activation of a specialized 'bottom-up' pathway. Light touch and caress activates a particular class of Ctactile mechanosensory afferents in hairy skin 3,8 . Activation of these fibers leads to decreased levels of the stress hormone cortisol, mediates release of the neuropeptide oxytocin and modulates the activity of insular cortex, a key region for emotional regulation [9][10][11][12] . However, two major lines of evidence also point to a major role of 'top-down' modulation in the processing of social touch.
First, unlike discriminative touch, the psychophysics of social touch are extremely dependent on individual, social and sexual context. For example, even though the actual tactile input is essentially identical, the touch of a loved one can feel comforting and pleasant while the touch of someone repulsive can feel intensely aversive 4,8,13 . Second, activation of the C-tactile fibers is not required to make social touch 'special'. Interpersonal touch delivered by the glabrous skin of the palm, thus not activating the C-tactile fibers, elicits a social softness illusion in the giver 14 and has a different neural activation pattern from non-social touch [15][16][17] .
Pioneering functional imaging studies in humans have investigated which cortical regions contribute to the topdown modulation of social touch processing by identifying brain regions, where activity patterns depend on the social 'meaning' of the touch, rather than the actual haptic input. For example, although in fact an identical pattern of touch was always given by the same experimenter, activity in anterior cingulate, orbitofrontal, somatosensory and insular cortices is different when subjects believe they are being touched by a man or a woman 18,19 or by a partner or a stranger 20 . This social-context-dependent modulation of cortical touch responses is negatively correlated with autisticlike traits 19 and is increased by intranasal oxytocin administration 20 . In line with a role of these regions in top-down modulation of social touch processing, somatosensory cortex, for example, is activated in a social context before any actual touch input 16,17 , when observing others being touched 21 , and when simply imagining pleasant or sexual touch 22 .
In this paper, we apply techniques with cellular resolution in the rat, which has been pioneered in the primate by Romo and colleagues [23][24][25][26] , namely to investigate and compare how single cortical neurons across multiple cortical areas respond during the same somatosensory stimulus. We focus on rat social facial touch, during which rats align their snouts and palpate each other's faces with their whiskers 27 . Social facial touch has attractive behavioral characteristics.
First, the behavior is ecologically valid. The animals are untrained, social interactions are jointly initiated by the animals themselves and the animals are freely moving, thus their cortical activity presumably closely resembles activity in a natural social setting. Second, by letting animals interact with partners of both sexes, we can manipulate the social context of the touch, while keeping the actual haptic input identical or similar. The socio-sexual 'meaning' of touching male and female conspecifics is different, but during social facial interactions, male rats whisk with equal power onto conspecifics of both sexes. Females whisk onto females like a male rat would, but whisk with lower whisking amplitudes onto males 27 . Similar to humans, previous work has shown that even though whisking amplitude is lower during social facial interactions than when investigating objects, population firing rate changes 28 and membrane potential modulations 29 in rat somatosensory cortex are larger during social touch than object touch and do not correlate with whisker movements. Also similar to humans, rat somatosensory activity is modulated in a social context before actual social facial touch 29 and modulation depends on socio-sexual context, such as estrus state 28,30 and emotional state 31 .
We present a flexible regression approach, which allows us to ask how social touch impacts activity, despite the large behavioral variability in social interactions. We ask the following questions: ( ) How is information about touch and social context (partner sex) represented at the level of single neurons? ( ) How widely is this information available across five different cortical areas? ( ) How does social context impact population dynamics during touch? ( ) How does the population response structure depend on cortical area, partner sex and sex of the subject animal itself?

Data
In this study, we analyze 1,156 single neurons, recorded from five cortical areas, over 7,408 episodes of social facial touch in 28 rats (58,591 unique cell-touch pairs, averaging 51 touch episodes per cell).

Diverse responses to social facial touch
To investigate how social touch modulates cortical processing across frontal and sensory cortices, we implanted tetrodes to record single-unit responses from freely moving, socially interacting male and female rats in the 'social gap The social gap paradigm 27. Rats separated on two platforms will reach across the gap to engage in social facial touch. In social facial touch, rats align their snouts and whisk to palpate each other's face with their mystacial vibrissae. We recorded the behavior of male and female rats by videography from above in visual darkness under infrared illumination 28,32,33 . (B) Anatomical location of vibrissa motor cortex ('VMC' -the primary motor representation of the mystacial vibrissae), cingulate cortex ('ACC' -a putative homolog of human/primate anterior cingulate cortex), prelimbic cortex ('PrL' -a putative homolog of human/primate medial prefrontal cortex), barrel cortex ('S1' -the primary somatosensory representation of the mystacial vibrissae) and auditory cortex ('A1'). (C) Example peri-stimulus time histograms (PSTHs) of touch-activated single neurons (Wilcoxon signed-rank test), aligned to the first whisker-to-whisker touch in each social touch episode. Raster plot above shows single trials, black dots indicate single spikes and vertical line indicates beginning of social facial touch. Black line below indicates mean firing rate, smoothed with a Gaussian kernel (s = 100 ms), shaded area indicates s.e.m. (D) Example PSTHs of touch-suppressed single neurons (Wilcoxon signed-rank test), aligned to the first whisker-to-whisker touch in each social touch episode. Same plotting convention as in (C). (E) Top: Histograms of touch duration (red) and inter-touch-time (blue) of social facial touch episodes. The time axes of both plots are clipped at 10 s. Below: Probability distributions of touch duration (red) and inter-touch-interval (blue) of social touch episodes, plotted on a logarithmic time scale. The two distributions span multiple orders of magnitude and strongly overlap. paradigm' (Figure 1A). In this paradigm, rats reach across a gap between elevated platforms to engage in social facial touch; an un-trained, self-initiated social behavior where the animals align their snouts and palpate each other's faces with their whiskers 27 .
We recorded the activity of neurons throughout the cortical column in five cortical areas: Two sensory areas, barrel cortex ('S1', the primary somatosensory representation of the mystacial vibrissae) and auditory cortex ('A1'), and three frontal areas, vibrissa motor cortex ('VMC', the primary motor representation of the mystacial vibrissae), cingulate cortex ('ACC', a putative homolog of human anterior cingulate cortex) and prelimbic cortex ('PrL', a putative homolog of human medial prefrontal cortex) (Figure 1B). To investigate how cortical processing of social facial touch depends on the subject sex and the sex of the social interacting partner, we recorded from both male and female rats, interacting with multiple male and female conspecifics (see Methods). A portion of the data analyzed here has been presented in previous studies where we investigated other questions (ref. 28,32,33 , see Methods).
We plotted the activity of single neurons, across all social interactions, aligned to the onset of whisker-towhisker touch (peri-stimulus time histograms, 'PSTHs'). From inspecting these plots, we made two preliminary observations. First, we noticed that cortical responses to social facial touch are widespread and diverse: In all brain areas, we both found neurons, which significantly increased ( Figure 1C) and significantly decreased ( Figure 1D) their activity at the onset of social facial touch. Second, we noticed that the PSTHs were highly variable from trial to trial. This variability could reflect a low correlation between firing patters and behavior. However, if the neural responses are tightly locked to behavior, but the behavior itself is highly variable, we would expect the same pattern. To discern these possibilities, we investigated the temporal statistics of the social facial touch episodes.
The median duration of a social facial touch was 1.33 s (IQR: 0.77-2.53 s, Figure 1E, top) and the median inter-touchinterval (ITI) was 4.68 s (IQR: 1.24-23.78 s, Figure 1E, middle), but both were Raster plot of example touch (activated S1 L5b neuron and suppressed VMC L5b neuron) and sex-touch neurons (activated VMC L2/3 neuron and suppressed VMC L6 neuron). Raster plots show spike times (black dots) aligned to the first whisker-to-whisker touch in each social touch episode. Social touch episodes are sorted by partner sex (female: pink, male: blue) and by duration (indicated by length of colored bar). Many touch episodes happen close together in time and there is a large variability in the touch duration. (B) Peri-stimulus time histograms of the example neurons shown in (A), separated by partner sex. Black line indicates mean firing rate, smoothed with a Gaussian kernel ( = 100 ms), shaded area indicates s.e.m, pink/blue color indicates female/male partner animals. (C) Peri-stimulus time histograms of the example neurons shown in (A), calculated from the fitted regression model, shown for comparison (plot conventions as in (B)). (D) Estimating touchmodulation: Log-likelihood values of models fitted to the neurons shown in (A). The log-likelihood of models depending on touch is indicated by the green arrow, the log-likelihood of models without touch is indicated by the grey arrow, and the log-likelihood distribution of shuffled touch-models is indicated by green bars. All neurons are significant at p < 0.05 (the green arrow is outside shuffled distribution). (E) Estimating sex-touch-modulation: Log-likelihood values of models fitted to the neurons in (A). The log-likelihood of models depending on both partner sex and touch is indicated by brown arrow, the log-likelihood of models without sex is indicated by green arrow and the log-likelihood distribution of shuffled sex-touch-models is indicated by brown bars. The two touch neurons are not significantly modulated by sex (the brown arrows are inside shuffled distribution), both sex-touch neurons are significant at p < 0.05. (F) Number of neurons that modulated by touch ('touch neurons', green color) and neurons that are modulated by touch, but respond differently to male and female conspecifics ('sex-touch neurons', pink/blue striped color) across cortical areas. (G) Mosaic plot of the distribution of touch neurons, sex-touch neurons and non-significant neurons across cortical areas (p-value indicates χ2-test of independence). Colors indicate significantly increased (blue) and decreased (red) proportions (standardized Pearson residuals at p < 0.05). highly variable. The duration of touch episodes varied across two orders of magnitude, the ITI varied across four orders of magnitude and the distributions of touch durations and inter-touch-intervals were highly overlapping ( Figure 1E, bottom). In other words, social facial touch episodes happen in bouts where the animals engage in several short touches, separated by inter-touch-intervals, which are often on the same order of magnitude of the touch durations themselves. This behavioral observation suggests that a PSTH-based analysis ( Figure 1C-D) will under-estimate the magnitude and over-estimate the trial-to-trial variability of neuronal responses to social facial touch. The 'baseline' period (-2000 to 0 ms) is not 'clean', but contaminated with touch episodes, and the post-stimulus period (0 to 500 ms) contains a mixture of trials with a wide range of touch durations.

Touch and sex-touch responses across cortex
To overcome the challenges presented by the temporal statistics of naturalistic social interactions, we used a generalized linear regression approach. Briefly, we modeled the spiking activity of the neurons as a Poisson process with history dependence and baseline fluctuations 34,35 Figure S1A). We used a maximum likelihood approach combined with a non-parametric shuffling procedure to investigate the effect of social touch and the sex of the social interaction partner, while maintaining the information about duration and inter touch interval of every single touch episode (see Methods).
Using this approach, we identified neurons that were modulated by touch ('touch neurons'), and neurons that were modulated by touch, but responded differently when touching male and female conspecifics ('sex-touch neurons'). Touch and sex-touch neurons had very diverse response patterns. We both found touch and sex-touch neurons that increased and decreased their firing rates during social touch (Figure 2A-E, left), and sex-touch neurons that responded more strongly to either female or male conspecifics ( As we expected based on the temporal patterns of the behavior ( Figure 1E), a lot of the inter-trial variability in the baseline period was due to variations in behavior. For example, many of the spikes in the baseline period of the example activated touch neuron (a layer 5b neuron in S1) and many of the pauses in firing in the example suppressed touch neuron (a layer 5b neuron in VMC) coincided with social touch episodes ( Figure 2A). Similarly, much of the inter-trial variability in the post-stimulus period was due to variations in the length of the social touch ( Figure 2A).
Across all investigated areas, we found that a large proportion of neurons were modulated during social facial touch episodes. On average, across all areas, we found 25.7% touch neurons, and 11.9% sex-touch neurons ( Figure 2F). We found that the proportions of touch and sex-touch neurons depended on the brain area (p < 0.001, 2 -test of independence, Figure 2G) and a mosaic plot 36,37 revealed that this was driven by the fact that AC had more sex-touch neurons than other brain areas and PrL had less sex-touch neurons, less touch-neurons and more non-significant neurons than other brain areas (all: |standardized Pearson residual| > 2, all: p < 0.05, Figure 2G). This suggests that, although there are differences in the proportions, information about touch and sex of the interaction partner is available across all investigated brain areas.

Cortical neurons carry significant sex information
In a typical experiment, a subject animal would interact with at least two male and two female conspecifics (see Methods). In our analysis so far, we have naïvely grouped female and male interaction partner rats together and identified sextouch neurons, which respond differently depending on the sex of the interaction partner. We wondered if this observation was really due to the sex of the partner animal, or if the neurons perhaps respond in an individual-specific manner. If the neurons did not encode the sex of the stimulus animals at all, but had individual-specific responses, we might also expect to identify artefactual 'sex-touch' neurons, simply because we are comparing two groups of animals with individual responses.
To investigate if the neurons really encode sex, we used an information theoretical approach. We calculated the information per spike 38 and used a shuffling procedure to identify putatively individual-specific neurons (see Methods). These neurons are a subset of our sex-touch neurons, most likely to carry individual-specific information. To ask if the responses of these neurons are individual-specific or sexspecific is difficult at the single-cell level. One way to frame the question is to ask if grouping the neurons by the real sex yields a higher information content than other groupings of the data. However, if we presented two male and two female partners, we can only perform two shuffled partitions of the data per single cell ( Figure S3A, top). If we presented more animals, we can perform more shuffles ( Figure S3A, bottom), which creates dependence between data points across the dataset. To overcome this imbalance, we used a mixedeffects modeling approach 39 (see Methods) and found that putatively individual-specific neurons were significantly more informative about sex than all other possible partitions of the data (mean DInfo = 0.007 bits/spike, p = 0.0027, mixed-effects model, Figure S3B). This analysis does not exclude the possibility that some neurons with individual-specific responses could be present in these brain areas. It shows that the sex of the partner animal is indeed a major overall determinant of the firing patterns, and that partner-sex-specific modulation is not an artifact better explained by individual-specific effects.

Prototypical response patterns vary by cortical area
Based on the diversity of touch and sex-touch responses, we wondered if the different cortical regions had different population response patterns. We plotted the estimated touch regression coefficient ( ℎ ) of all neurons by area and found that -as a population -S1 neurons significantly increased their firing rate during social touch (mean ℎ = 0.17, p = < 0.001, Wilcoxon signed-rank test, Figure 3A). VMC, ACC and A1 neurons significantly decreased their firing rate (VMC/ACC/A1: mean ℎ = -0.084/-0.077/-0.073, all p < 0.01, Wilcoxon signed-rank test, Figure 3A) and PrL neurons did not show any modulation at the population level (mean ℎ = -0.036, p = 0.20, Wilcoxon signed-rank test, Figure 3A).
When we only looked at the number of neurons where the touch coefficient was significant at the single cell level, we found that these significant neurons generally responded in the same direction as the population. In S1, more neurons were significantly increased by touch (38 v. 89 neurons, p < 0.001, binomial test, Figure 3B, top) and the significantly increasing neurons were the most strongly modulated neurons (p = 0.041, Mann-Whitney U test, Figure 3B, bottom). In VMC, we found the opposite picture. More neurons were decreasing (73 v. 37 neurons, p < 0.001, binomial test, Figure  3B) and the decreasing neurons were the most strongly modulated (p = 0.035, Mann-Whitney U test, Figure 3B). Similarly to VMC, ACC and A1 had more significantly decreasing neurons (ACC: 18 v. 7 neurons, A1: 53 v. 25, both p < 0.05, binomial test, Figure 3B). Also like for VMC, there was a tendency for the decreasing neurons in ACC and A1 to be more strongly modulated, but this difference was not significant in our dataset (ACC: p = 0.11, A1: p = 0.29, Mann-Whitney U-test, Figure 3B). In line with the lack of a population response in PrL, we found that both the number of increasing and decreasing PrL neurons (14 v. 10 neurons, p = 0.12, binomial test) and their modulation strengths were similar (p = 0.75, Mann-Whitney U test, Figure 3B).
From the profile of the peri-stimulus time histogram, we estimated the onset time of the modulation of the firing rate of single cells ( Figure  S4A, see Methods). In line with previous investigations 26,40,41 , we found that neurons in somatosensory cortex generally responded with a shorter latency than neurons in other cortical areas (p = 0.017, Kruskal-Wallis test of independence, S1 v. all p < 0.05, except S1 v. PrL: p = 0.083, Mann-Whitney U-tests, Figure S4B), which may suggest a role of cortico-cortical connections from somatosensory cortex in driving the modulation of neurons in the other cortical areas during social touch.
From inspecting the PSTHs, we noticed that some responses appeared qualitatively more transient (responding at the onset of social touch) whereas some appeared more tonic (responding constantly throughout the social touch episode) (Figure 1-2). Thus, we wondered if there might be systematic differences in the temporal profile of the touch modulation between areas or between increasing and decreasing neurons. For example, neurons in some cortical areas might have mainly transient neural responses while neurons in other areas might be mainly tonically modulated. The temporal profile could also be different between increasing and decreasing neurons. We used a principle component analysis of the temporal profile of the peri-stimulus time histograms ( Figure S5A, see Methods) and found that the response profiles did not cluster into distinct clusters of either tonically or transiently modulated neurons. Rather, we found that the temporal profiles of both increasing and

Population dynamics depend on subject sex and partner sex
To investigate if the cerebral cortices of male and female rats process social touch differently, we first plotted all touch and sex-touch neurons recorded in the respective areas and separated neurons recorded in males and females ('male neurons' and 'female neurons'). Across areas, we found that male and female neurons respond similarly to touch and mirror the population responses described in Figure 3. In S1, firing rates of both male and female neurons were increased by social touch and there was no difference between the sexes ( , < 0.001, In VMC, ACC and AC we found that both male and female neurons were decreased by touch, and we found no differences between the sexes there either (all: , < 0.05, male v. female p > 0.05, Figure 4A).
In PrL, neither male nor female neurons were significantly modulated as a population, and there was no difference between the sexes ( , > 0.05, male v. female p > 0.05, Figure 4A).
Next, we investigated if population responses depend on the sex of the interaction partner. In the posterior parietal cortex (PPC), there is evidence that modulation by sensory stimuli and movement features is distributed randomly across neurons with no or little population structure 42 . We plotted the modulation of neurons during interactions with male and female conspecifics. Colored dots indicate neurons recorded in female (red) and male (blue) animals. In S1, firing rates of both male and female neurons were increased by social touch (median male βtouch = 0.32, pmale = 0.00011, median female βtouch = 0.30, pfemale = 0.00051, male v. female p = 0.38, N = 165). In VMC, ACC and AC, firing rates of both male and female neurons were decreased by touch (VMC: median male βtouch = -0.14, pmale = 0.0039, median female βtouch = -0.25, pfemale = 0.017, male v. female p = 0.61, N = 125, ACC: median male βtouch = -0.26, pmale = 0.046, median female #$%&' = -0.14, pfemale = 0.048, male v. female p = 0.19, N = 29, A1: median male βtouch = -0.18, pmale = 0.025, median female βtouch = -0.10, pfemale = 0.0042, male v. female p = 0.23, N = 115). In PrL, neither male nor female neurons were significantly modulated as a population, and there was no difference between the sexes (median male βtouch = -0.019, pmale = 0.67, median female βtouch = -0.091, pfemale = 0.21, male v. female p = 0.76, N = 34). All: t-tests (pfemale, pfemale) and unpaired t-test with unequal variance (male v. female) if normal by a Lilliefors test, else Wilcoxon signed rank tests and Mann-Whitney U-test, respectively. * indicates p < 0.05, ** indicates p < 0.01, *** indicates p < 0.001 (B) Modulation of activity (in fold change) during social touch with male and female conspecifics is highly correlated (VMC/S1/A1: Kendall's = 0.40/0.38/0.32, all p < 10 -23 /10 -26 /10 -13 ). Touch neurons are indicated by green dots, female/male preferring sex-touch neurons are indicated by pink/blue dots and non-significant neurons are indicated by grey dots, Kendall's and p-value above. (C) Top: A biased response to one partner sex corresponds to a shift of the regression line. For example, the red line corresponds to a situation where neurons always fire more spikes when touching males than females. Bottom: A potentiated response to one partner sex corresponds to a change in slope of the regression lin. For example, the red line corresponds to a situation where neuronal responses to female conspecifics are always larger in magnitude than responses to male conspecifics. (D) Population response pattern depends on subject sex and partner sex. Dots indicate neurons recorded in female (red) and male (blue) subject animals, lines indicate maximum-likelihood fit of regressing modulation with males as a function of modulation with females (red = female, blue = male, shaded area indicates 95% C.I.) * indicates slope different from unity (outside 95% C.I. for both males and females), ** indicates p < 0.01, *** indicates p < 0.001 (full model specification in Suppl. Note 1).
Similar to the PPC, we found that responses were very diverse and that neurons populated all quadrants (Figure 4b,  S6). Sex-touch neurons were both in the first and third quadrants (corresponding to differences in magnitude of modulation) and second and fourth quadrants (corresponding to different directions of modulation, e.g. increasing during interactions with males and decreasing during interactions with females). As expected from the definition, touch neurons were generally near the diagonal and non-significant neurons were near the origin ( Figure 4B). First, we focused on the three brain areas where we had a substantial dataset (VMC, S1 and A1) and found that responses to male and female interaction partners were highly correlated (VMC/S1/A1: Kendall's = 0.40/0.38/ 0.32, all p < 0.001, Figure 4B).
The high correlation values suggest that the population response to male and female stimuli is nonrandom. But what is the structure of the population response, and does it depend on the sex of the subject animal? High correlation values could reflect that there is no difference between responses to males and females, which would correspond to having all neurons along the diagonal (with some noise). If responses to males and females are different, the responses could be biased towards one sex, responses when touching one sex could be potentiated, or it could be a combination of both. A biased response would for example correspond to having all neurons spike more when interacting with males than females ( Figure 4C, top). A potentiated response would for example mean that responses when touching females are always stronger in magnitude than responses when touching males ( Figure 4C, bottom). Finally, these response patterns might depend on the sex of the subject animal. To identify the population pattern, we performed mixed-effect generalized linear modeling (see Methods). We did not find any evidence of a bias in responses (all: > 0.05) and no subject-sex-dependent bias (all: _ > 0.05, full model specification in Suppl. Note 1). However, in all three areas, we found evidence of potentiation. There was a highly significant dependence of responses to males on the response to females (all: _ < 0.001) and the estimate of the slope was significantly smaller than unity, suggesting that responses to females are stronger in magnitude than responses to males (VMC/S1/A1: Figure 4D). In S1 and VMC, there was also a significant interaction between the subject sex and this potentiation (both: _ * _ < 0.01, Suppl. Note 1). Interestingly, the interaction had opposite signs (VMC/S1: Figure 4D).
We also found significant correlation values between responses to male and female conspecifics in the areas where we only had a limited sample (ACC/PrL: Kendall's = 0.22/0.17, both p < 0.01, Figure S7A). In these smaller datasets, there was no statistically significant subject-sex-dependent or stimulus-sex-dependent potentiation, but -although not significant -the maximal-likelihood fit also estimated _ to be less than unity for both ACC and PrL ( Figure S7B, Suppl. Note 1). However, in ACC we found a significant bias, indicating that ACC neurons fire less with male than female conspecifics, in both male and female subjects ( < 0.05, Suppl. Note 1). In summary, for both male and female subjects, cortical responses were larger in magnitude with female than male interaction partners (statistically significant in VMC, S1 and A1, same direction, but not significant in ACC or PrL). In S1 the potentiation of responses to females was stronger in female subjects than male subjects. In VMC the potentiation was stronger in male subjects and in A1 the potentiation did not depend on the subject sex.

Summary
We developed a flexible regression approach that allowed us to investigate how social context (partner sex and subject sex) impacts cortical processing during naturalistic social touch. Social context modulates the firing patterns of individual cortical neurons in a highly structured manner. Information about social touch was available across all cortical regions investigated and partner sex was a major determinant of touch responses. Although every area investigated has access to touch information, touch modulates different brain structures in different ways.

Population patterns in cortical networks
A previous study found that firing rate modulation by sensory stimuli and movement features in the posterior parietal cortex was distributed randomly across neurons with no or little population structure 42 . In contrast to this observation, we found that social touch responses were highly structured, and that the network activity pattern signaled the sex of the interaction partner by partner-sex-dependent potentiation of responses. Moreover, in somatosensory and vibrissa motor cortex we found that the sex of the subject animal itself was encoded by a subject-sex-dependent magnitude of this potentiation. In line with the encoding of partner sex observed, there is previous evidence of population coding of touch stimuli in the somatosensory cortex 43,44 . Previous investigations in the rodent whisker system have found that putatively 'simpler' stimuli presented in highly controlled conditions, such as object location 45,46 and texture roughness 47 , are encoded with higher fidelity by precise spike timing than by firing rates. However, in line with our observations, previous investigations also found that firing rates carry extra information in addition to the information conveyed by spike timing 48-50 . The relevant time scale for temporal coding by spike timing in the somatosensory cortex seems to be within tens of milliseconds from the stimulus, and the analysis requires extremely precise information about stimulus onset time 46,51,52 , which is not available in our naturalistic experimental paradigm. More generally, structured population dynamics in vibrissa motor cortex aligns well with other observations of population coding in motor cortex, where the population activity vector of both preparatory 53 and movementrelated activity 54,55 correlates with movement features [56][57][58] .

Social touch responses in frontal cortex
In line with our observations of touch-and sex-touch neurons in rat cingulate cortex, previous human studies have found that social context strongly modulates responses in cingulate and orbitofrontal cortex 18-20 . We found that -when examined at single cell resolution -the majority of cingulate neurons decrease in activity during social touch. In light of the tight anatomical integration of prelimbic cortex with brain structures responsible for social behavior and emotions 59,60 and the putative homology with human and primate prefrontal cortex, a key structure in social cognition 61,62 , we found it curious that prelimbic cortex was so weakly modulated by social touch. We did not see any overall changes in firing rate, only weakly modulated touch and sex-touch neurons at both edges of a zero-centered distribution. While previous rodent studies did not investigate social touch as such, in line with our observations, it has been shown that during social approach (to a caged conspecific) prelimbic responses are weak and diverse 63,64 and that social behaviors are encoded by sparse groups of neurons either increasing or decreasing in activity during behavior ('on' and 'off' ensembles 65 ).
In line with what we described for other whisking behaviors 33 , we found that despite active whisker movement, neurons in vibrissa motor cortex mainly decrease their activity during social touch. We also found a large proportion of sextouch neurons in vibrissa motor cortex. We already know that neurons in rodent motor cortex can have low-latency responses to touch 40, [66][67][68][69][70] , that motor cortex is an important node in a distributed network for touch processing [71][72][73][74][75] and that motor cortical activity modulates somatosensory processing in both somatosensory cortex [76][77][78] and thalamus 79 . Vibrissa motor cortex has been implicated in diverse aspects of sensorimotor cognition (for reviews, see 80,81 ) and our study adds a potential role in sex-dependent processing of social cues to this complex picture.

Sex differences in cortical processing of social touch
In somatosensory and motor cortex, we found that social touch leads to different network activity patterns in male and female subjects. Generally, while there clearly are some systematic sex differences in both anatomy [82][83][84][85][86][87] and functional connectivity 88 , male and female cortices are overall remarkably similar. For example, even though male and female genital anatomy is extremely different, the layout of the somatosensory body map is essentially identical in both sexes 89,90 and has similar projection targets 91 . In recent years, several studies have identified sex differences in subcortical circuits involved in regulating social behavior. For example, circuits involved in the control of aggression in the ventromedial hypothalamus are wired differently in males and females 92 , medial amygdala has striking sex differences in olfactory responses 93 , and galanin-positive neurons in the medial preoptic area 94 and their subcortical input nuclei 95 are activated during parental behavior in a sex-and reproductivestate-specific manner. In contrast to these subcortical examples, sex-differences in cortical processing are rare and subtle 62,96,97 .
Our data invites two mutually not exclusive interpretations. On one hand, the differences in processing may be a signature of a universal sex difference in how male and female cortices process the same sensory stimuli. Thus, these sex differences might generalize across other cortical regions and sensorimotor modalities. On the other hand, male and female cortices may use identical computational strategies, and our observed differences simply reflect the fact that for male and female subjects, the same partner sex presents a different social situation with different cognitive and behavioral challenges. For example, a male rat interacting with another male rat might be mainly assessing dominance and aggression, while a female rat meeting a male rat might be assessing aggression as well as reproduction 27, 98 .
We still know very little about sex differences in human cortical processing of social touch stimuli. As noted in the introduction, our investigation of partner-sex dependent differences in touch responses was based on pioneering studies with human subjects 18-20 . These studies have not reported sex differences in the processing of social touch stimuli, but then again, previous studies have not compared experimental subjects of both sexes 13,[16][17][18][19][20][21][99][100][101] . It would be most interesting to know if our observed network activity differences in the rat are paralleled in primates and humans. Our findings highlight the importance of complementing human work with the single cell resolution offered by animal studies.

Cellular mechanisms underlying differences in male and female responses
During social facial interactions, the actual haptic input from male and female partners is very similar 27 . However, our naturalistic paradigm is inherently multisensory, and even though the whisking and touch input is similar, male and female partners convey very dissimilar olfactory cues. The vomeronasal organ is important for determining the sex of the conspecific in both male-male 102 and male-female interactions 103 , and might be an important 'bottom-up' pathway for the modulation of cortical responses during social touch observed here.
We found that both partner sex and subject sex was encoded in a potentiation of social touch responses. This pattern suggests that a potential mechanism underlying the sexdependent modulation of cortical responses could be a change in inhibitory drive. Such a change could explain why we observed that both increasing and decreasing neurons changed the magnitude of their responses in tandem. For example, increased responses in parvalbumin-positive interneurons would presumably be accompanied by increased suppression of somatically inhibited principal cells. This idea aligns well with other observations showing that inhibitory neuron subtypes control context-dependent modulation of sensory responses in visual 104,105 and auditory cortex 106 (for review see 107 ). A likely candidate mechanism for the regulation of cortical interneuron activity during social touch is by neuromodulatory hormones, such as estrogen and oxytocin. We recently found that deep layer parvalbumin-positive interneurons in rat somatosensory cortex express estrogen receptor β and that fast-spiking interneuron activity changes with the estrus-cycle and estradiol injection 30 . Oxytocin action on neural circuits plays a major role in cognitive and emotional aspects of sociosexual behavior, such as social bonding, anxiety and trust 108,109 and oxytocin receptors are expressed in cortex mainly by interneurons 110,111 . Oxytocin acts on interneurons in auditory cortex to enable pup-retrieval by mother mice 112 and modulates interneuron activity in prefrontal cortex 111 and olfactory cortex 113 to enable social recognition. Gentle touch and stroking activates oxytocinergic neurons in the paraventricular hypothalamus 114 which project directly to sensory cortices 110 . This provides a potential circuit, by which touch-related oxytocin release can impact cortical interneuron activity to modify cortical responses and network activity patterns during social facial touch. Moreover, since the pattern of cortical oxytocin receptor expression is sex-dependent 110,115 , this provides a potential basis for our observed sex differences in network activity patterns.

Conclusion
Across multiple cortical regions, social context (partner sex) is a major determinant of single cell responses and population dynamics. Moreover, activity patterns depend on the sex of the subject animal. Identifying potential differences in how male and females cortices process social stimuli could help shed light on the biological basis of sex differences in the etiology, prevalence and symptoms of autism, depression and anxiety disorders, all characterized by social dysfunction [116][117][118] .

Competing interest statement
The authors declare no competing interests.

Animal welfare
All experimental procedures were performed according to German animal welfare law under the supervision of local ethics committees. Rats (Wistar) were purchased from Janvier Labs (Le Genest-Saint-Isle, France). Rats presented as partner animals were housed socially in same-sex cages, and postsurgery implanted animals were housed in single animal cages. All animals were kept on a 12h:12h reversed light/dark cycle and all experiments were performed in the animals' dark phase. Rats had ad libitum access to food and water.

Social gap paradigm
Behavioral experiments were done using the social gap paradigm 27 ( Figure  1A). The experimental paradigm consists of two elevated platforms, 30 cm long and 25 cm wide surrounded by walls on 3 sides, positioned approximately 20 cm apart. The distance between platforms was varied slightly depending on the size of the rats. The platforms and platform walls were covered with soft black foam mats to provide a dark and non-reflective background and to reduce mechanical artifacts in tetrode recordings. All experiments were performed in darkness or in dim light, and behavior was recorded from above under infrared light. The implanted rat was placed on one platform, and on the other platform we either presented various objects or other rats. The implanted rats were not trained, just habituated to the setup and room, and spontaneously engaged in social facial interactions. The rat behavior was recorded at low video speed from above with a 25 Hz or 30 Hz digital camera, synchronized to the electrophysiological data acquisition using TTL pulses. Typically, recording sessions were performed in four to eight 10-15 min blocks, where we would present either one or two conspecifics (of both sexes) in each block, randomly, see 28,32,33 . Video frames were labeled blind to the spike data.

Tetrode recordings and histology
In tetrode recording experiments, we used ~p60 Wistar rats, which were handled for 2-3 days before being implanted with a tetrode microdrive. Surgery was done as previously described 28 . The implanted microdrive had eight separately movable tetrodes driven by screw microdrives (Harlan 8drive; Neuralynx, Bozeman, MT, USA). The tetrodes were twisted from 12.5 µm diameter nichrome wire coated with polymide (California Fine Wire Company), cut and examined for quality using light microscopy and goldplated to a resistance of ca. 300 kOhm in the gold-plating solution using an automatic plating protocol ("nanoZ", Neuralynx). For tetrode recordings targeting VMC, ACC and PrL, a craniotomy of 1x2 mm was made 0.75-2.75 mm anterior and 1-2 mm lateral to bregma 119,120 . ACC includes data from the area 'Cg1' and dorsal 'Cg2' in the atlas 119 . A1 includes data recorded in the rat atlas regions 'Au1', 'AuD' and 'AuV' 119 . Steel screws for stability and two gold-plated screws for grounding the headstage were drilled and inserted into the skull, and the gold-plated screws were soldered and connected to the headstage PCB using silver wire. After fixation of all screws, the dura was removed, the implant fixated above the craniotomy, the craniotomy sealed with 0.5% agarose and the tetrode drive fixed in place with dental cement (Heraeus). The tetrodes were arranged in a 2-by-4 grid with a spacing of ~500 µm. Neural signals were recorded through a unity-gain headstage preamplifier and transmitted via a soft tether cable to a digital amplifier and A/D converter (Digital Lynx SX; Neuralynx) at 32 kHz. We filtered the signal between 600 Hz and 6 kHz and detected spikes by crossing of a threshold (typically ~50 µV) and saved each spike (23 samples, 250 µs before voltage peak to 750 µs after voltage peak). At the end of the experiment, animals were again anaesthetized with a mix of ketamine and xylazine, and the single tetrode tracks were labelled using small electrolytic lesions made by injecting current through the tetrode wire (10 µA for 10 s, tip-negative DC). After lesioning, animals were perfused with phosphate buffer followed by a 4% paraformaldehyde solution (PFA). Brains were stored overnight in 4% PFA before preparing 150 µm coronal sections. Sections were stained for cytochrome oxidase to reveal the areal and laminar location of tetrode recording sites, which could be calculated from the location of tetrode tracks and lesions. We only analyzed data from recording sites where the lesion pattern could unanimously identify the tetrode and the recording sites.
All spike analysis was done in Matlab (MathWorks, Natick, MA, USA). Spikes were preclustered off-line on the basis of their amplitude and principal components by means of a semiautomatic clustering algorithm ('KlustaKwik', https://github.com/klusta-team/klustakwik , ref. 121 ). After preclustering, the cluster quality was assessed and the clustering refined manually using MClust (http://redishlab.neuroscience.umn.edu/MClust/-MClust.html, A. D. Redish, University of Minnesota). The spike features used for clustering were energy and the first principle component of the waveform. To be included in the analysis as a single unit, clusters had to fulfill the following criteria: first, the L-ratio, a measure of distance between clusters, was below 0.5. Second, the histogram of inter-spike intervals (ISIs) had to have a shape indicating the presence of single units, e.g. a refractory time of 1-2 ms, or the appearance of a bursty cell (many short ISIs). Multiunit clusters were not included in the analysis.
In order to directly compare cortical areas with different cytoarchitectonics, we did not separate neurons into putative interneurons and putative excitatory neurons based on their extracellular spike shape. In some brain regions, for example the hippocampus, separation of cell type based on extracellular spike shape appears to work well and reliably identify thinspiked neurons which suppress simultaneously recorded neurons with short latency 122,123 . However, in other brain regions the use is more controversial 124 . For example, motor cortical pyramidal projection neurons have extremely thin spikes 125 and some cortical interneurons have very wide spikes 126,127 . Further, since the spike width and shape depends strongly on the location of the soma to the electrode 128 , spike shape recorded with tetrodes is likely less reliable with the lack of cytoarchitectonic stereotypy in the agranular rat frontal cortex compared to for example the hippocampus 120,129 . In line with this, we previously did not find a clear bimodal distribution of spike width allowing us to confidently identify putative GABAergic neurons in deep layers of vibrissa motor cortex 33 . In our previous analysis of the tetrode recordings from somatosensory cortex, on the other hand, we did see such a bimodal distribution which suggested that social touch responses of regular-spiking neurons in somatosensory cortex change with the estrus cycle 28 . However, in a follow-up study using juxtacellular recordings (which capture the spike shape with much higher fidelity than tetrodes), we were unable to reproduce that finding 30 . Clarifying if and how interneuronal activity shapes cortical responses is probably better tackled by future studies, using genetically encoded tools to image or tag neuronal subpopulations.

Previous use of data
Part of the data presented in this study has already been presented in other studies. The barrel cortex data, recorded throughout all cortical layers, has previously been published in ref. 28 . The auditory cortex data, recorded throughout all cortical layers, has previously been published in ref. 32 . The recordings from deep layers of vibrissa motor cortex has previously been presented in ref. 33 . The previous study on barrel cortex 28 investigated how touch-evoked activity depended on male and female subject and partner animals in a PSTH-based way (like our Figure 1C-D), the other two previous studies did not investigate any sex differences 32,33 . Neuons were recorded throughout the cortical column, with a majority of neurons in the deep layers. The laminar distribution was (S1/VMC/PrL/ACC/A1): Layer 6: 32/16/39/31/46, layer 5b: 88/153/70/17/100, layer 5a: 42/37/33/47/63, layer 4: 56/0/0/0/19, layer 2/3: 37/90/0/0/4, layer uncertain: 129/0/0/0/7. All data from the superficial layers of vibrissa motor cortex, data from all layers of cingulate cortex and data from all layers of prelimbic cortex have not previously been published.

Statistical modeling -PSTHs
To identify significantly increasing and decreasing neurons, as shown in Figure 1C-D, we calculated the mean firing rate in the 'baseline' period (-2500 to 0 ms) and the post-stimulus period (0 -500 ms) and determined significant changes using Wilcoxon signed-rank tests.

Statistical modeling -Touch and sex-touch neurons
We used a generalized linear regression approach 34,35 to identify touch and sex-touch responses. We discretize the spike train in 1-ms bins, reduced the data amount by removing baseline periods, which were more than five seconds from any social interaction and model the firing rate as a Poisson process. If we assume that the discharge of spikes within each time bin is generated by a homogenous Poisson point process, then the probability of observing y spikes in a single time bin is where ∆ = 1 ms is the width of the time bin and λ > 0 s-1 is the expected discharge rate of the cell. If we assume that each time bin is independent, the probability of the entire spike train, ̅ is where , is the observed number of spikes and the expected discharge rate in the i'th time bin, respectively. If we model the expected discharge rate, ̅ , so that it depends on the parameters, ̅ , we have the log-likelihood function We model ̅ so that it depends linearly on spike history, experimental recording, touch and partner sex and -since the expected firing rate cannot be negative -we model where P is a predictor matrix and ̅ is a vector of regression coefficients 34 . The predictor matrix has the following columns: a constant baseline rate; five 1-ms spike history bins; six 25-ms history bins; ( -1) one-hot columns to model a change in baseline between the recordings; a one-hot column indicating all social touch episodes; a one column indicating the sex of the stimulus animal (0 = female, 1 = male). Due to the refractory period of the cell, it is not correct to assume that all time bins are statistically independent, so following 35 , we include 11 spike history parameters, h1…h11, to model the inter-spike-interval distribution of the cell. The spike history is binned to 11 successive bins, five 1-ms bins (vectors with no. of spikes in the previous 0-1 ms, 1-2 ms, 2-3 ms, 3-4 ms, 4-5 ms) and six 25-ms bins (vectors with no. of spikes in the previous 5-30 ms, 30-55 ms, 55-80 ms, 80-105 ms, 105-130 ms, 130-155 ms). We also include constant bias terms to allow for variations in baseline firing rate between each recording. Thus, ] .
An example predictor matrix is shown in Figure S1A. We used the function package 'neuroGLM' (https://github.com/pillowlab/neuroGLM , ref. 130 ) to calculate and numerically fit the models using Matlab (MathWorks, Natick, MA, USA). To assess the statistical significance of touch and sex-touch responses, we used a non-parametric, shuffling-based model selection approach. To assign statistical significance to , we compared the log-likelihood of the model including touch and sex as predictors with the distribution of log-likelihood values from the same models fitted to predictor matrices, where we randomly shuffled the label (male/female) of the partner animal (N = 100, Figure 2D). We do not assume that random effects have a Gaussian distribution. Rather, to assign statistical significance to ℎ , we compared the log-likelihood of the model including touch (and not sex) as predictors with the distribution of log-likelihood values from the same models fitted to predictor matrices, where we circularly permutated the column indicating where the social touch episodes had happened (N = 100, Figure  2E). We defined 'sex-touch neurons' as neurons with < 0.05, 'touch neurons' as neurons with > 0.05 and ℎ < 0.05, and 'non-significant' neurons as neurons with both > 0.05 and > 0.05. We compared the proportions of non-significant, touch and sex-touch neurons across cortical areas by calculating standardized Pearson residuals and visualized the contingency table as a mosaic plot 36,37 . Standardized residuals and mosaic plots were generated using the 'vcd' package 131 for R 132 .

Information theory
To avoid assuming anything about the direction of firing rate changes, we calculated the information per spike, where 〈 〉 is the average spike rate, is the spike rate during the i'th stimulus, and and is the probability of the i'th stimulus 38 . To identify neurons, which could potentially have individual-specific or sex-specific response patterns we treated as 'stimuli' the social touch episodes with all individual partner animals. Significance was assessed by a shuffling procedure, where we circularly shifted the timing of the social touch episodes over the recordings (N = 200, p < 0.05). To investigate if these putatively individual-specific neurons carry more information about the real sex of the partner animal than randomly assigned sex labels, we also used a shuffling procedure. For all neurons, we calculated the mutual information per spike with three situations, baseline, social touch with male partners and social touch with female partners, and calculated the same value, where we shuffled the sex of the partner animals ( Figure S3A). Since the number of possible shuffles depend on the number of partner animals, the dataset had strong dependency per neuron 39 , so we fitted a mixed effects model, where Δ is the real minus the shuffled value of the information per spike between the three 'stimuli'.

Estimation of response onset time
To estimate the response onset time, we estimated the average firing rate in 1-ms bins (convolved with a Gaussian with = 100 ms) in neurons, where there was a significant response in the PSTH (assessed as in Figure 1). We only considered social touch episodes, where there was no social touch in the preceding 1000 ms. Response onset time was defined as the first time bin after -100 ms where the firing rate escaped the 1% -99% confidence interval of the firing rate in the second just before touch, as estimated by fitting a Gaussian distribution. This method follows ref. 26 . It should be noted, however, that since we have a much more limited number of trials, our estimates of onset times are much more noisy. For example, since the width of the confidence interval depends on the variance in the 'baseline' period just before touch, and we have a limited number of trials, we may systematically estimate the onset latency as much larger than they 'really' are. Thus, while a relative ranking in onset time between areas is possible, the temporal estimates of the actual onset are probably over-estimates.

Principal component analysis of PSTH profiles
To perform the principal component analysis, we again only considered social touch episodes, where there was no social touch in the preceding 1000 ms, and selected neurons that had a significant touch response. This was assessed by a Wilcoxon signed-rank test comparing the firing rate in the 1000 ms before touch with the first 500 ms after touch at p < 0.05. We compiled a matrix of the mean firing rate from -1000 ms to +1000 ms around the onset of touch, down-sampled the traces to 20 ms bins, and whitened the matrix to have zero mean and unit variance per cell. We calculated the principal component coefficients, scores and explained variance using the built-in 'pca' function in Matlab. We estimated the number of clusters by fitting 2d Gaussian mixture models (with different numbers of components, and no constraints on the covariance) to the coefficients of the first two principal components and calculated the Bayesian information criterion (BIC) using the built-in 'fitgmdist' functions in Matlab.

Statistical modeling -Magnitude of responses
To estimate the average depth of modulation of increasing and decreasing neurons across both partner sexes (Figure 3), we calculated the fold modulation by touch as the ratio where ℎ , are the firing rates and ℎ is the fitted regression coefficient of the GLM models including touch (and not sex) as a predictor. To plot and compare the magnitude of increases and decreases, we calculated the numerical value of the base-2 logarithm and plotted the data as fold increases/decreases (e.g. log 2 (ratio) = 1 corresponds to a two-fold increase (double the firing rate), log 2 (ratio) = -1 corresponds to a two-fold decrease (half the firing rate), etc.). To estimate the modulation when touching male and female conspecifics (Figure 4), we calculated the fold modulation from the fitted GLM models as where ℎ , ℎ , are the firing rates and ℎ , are the fitted regression coefficients of the full GLM models including touch and sex (i.e. this is a different fitted value of ℎ than above, Touch mod. ≠ Female mod.). To determine the population response pattern of the modulation, we fitted the mixed-effects regression male_mod~1 + female_mod + subject_sex + female_mod * subject_sex + (1|subject), where male_mod and female_mod are the modulation ratios, subject_sex is a one-hot vector indicating the sex of the subject animal (0 = female, 1 = male) and subject is a categorical variable indicating the subject animal. To avoid biasing the regression fit by the few cells with extremely low firing rates (indicated by circles in Figure S6A), we removed neurons with a more than 32-fold increase/decrease from the fit. The (1|subject) is a constant error term per subject animal, which we add due to the dependence introduced by the fact that we have unequal numbers of neurons recorded from the different subject animals 39 . Figure S1: Spike train regression model. (A) Example predictor matrix for the generalized linear regression approach. We discretize the spike train in 1-ms bins, down-sample the data to five seconds around every social interaction and model the firing rate as a Poisson process 34,35 (see Methods).
The predictor matrix has the following columns (left to right): a constant baseline rate; five 1-ms spike history bins; six 25-ms history bins; three onehot columns to model a change in baseline between the (in this example case) four recordings; a one-hot column indicating all social touch episodes; a one-hot column indicating the sex of the stimulus animal (in this example case, the interaction partner animals were male in recording one and recording three). The vector indicates regression coefficients, which we fit by likelihood maximization. The log-likelihood of models depending on touch is indicated by a green arrow, the log-likelihood of models without touch is indicated by a grey arrow and the log-likelihood distribution of shuffled touch-models is indicated by green bars. Some sex-touch neurons would not be significant, if the partner sex was not considered in the model. For example, the ACC L5a neuron is suppressed by males, but shows (almost) no response with females, so the green arrow is not in the 0.05 fraction of the shuffled distribution if all touches are pooled, without parsing out the partner sex. (E) Estimating sex-touch-modulation: Loglikelihood values of models fitted to the neurons in (a). The log-likelihood of models depending on partner sex and touch is indicated by a brown arrow, the log-likelihood of model without sex is indicated by a green arrow and the log-likelihood distribution of shuffled sex-touch-models is indicated by brown bars. All neurons are significant at p < 0.05 (brown arrow outside the shuffled distribution).  . This method follows ref. 26 , except that since we have a much more limited number of trials, the estimates of onset times are more noisy (both larger and smaller than they 'really' are). (B) Distribution of onset times (green dots) for all neurons with significant modulation in the PSTH. Despite the noisy estimates, in agreement with 26,40,41 , we find that somatosensory neurons generally respond earlier during touch than other cortical areas (all comparisons are significant, except for prelimbic cortex which is close to significant: S1 vs. VMC/ACC/A1: all p < 0.05, S1 vs. PrL: p = 0.083, Mann-Whitney U-tests, pindependence = 0.017, Kruskal-Wallis test). Ring and whisker plots indicate median and 95% confidence interval of the median.

Figure S5: Principal component analysis identifies two segregated clusters of modulated neurons. (A)
Two observations motivated this analysis. From inspecting the PSTHs, we noticed that some responses appeared qualitatively more transient (responding at the onset of social touch), whereas some appeared more tonic (responding constantly throughout the social touch episode) (Figure 1-2 Figure 4D, but showing data from ACC and PrL. Although not significant, the maximal-likelihood fit also estimated βsubject_sex to be less than unity for both ACC and PrL. This pattern is in line with the pattern in S1, VMC and AC ( Figure 4D). Red/blue dots indicate neurons recorded in female/male subject animals, red/blue lines indicate maximum-likelihood fit of regressing modulation with males as a function of modulation with females, for female/male subjects (see Methods and Suppl. Note 1 for model specification), shaded area indicates 95% C.I.