# Dopaminergic mechanisms of dynamical social specialization in mouse microsocieties

C. Solié, A. Nicolson, R. Justo, Y. Layadi, C. Morin, B. Batifol, Lm Reynolds, T. Le Borgne, S.L. Fayad, A. Gulmez, et al.

► To cite this version:

C. Solié, A. Nicolson, R. Justo, Y. Layadi, C. Morin, et al. Dopaminergic mechanisms of dynamical social specialization in mouse microsocieties. 2025. hal-05346504

HAL Id: hal-05346504

https://hal.science/hal-05346504v1

Preprint submitted on 4 Nov 2025

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
1 Dopaminergic mechanisms of dynamical social specialization in mouse microsocieties

2 C. Solié¹, A. Nicolson¹, R. Justo¹, Y. Layadi¹, B. Morin¹, C. Batifol¹, LM Reynolds¹, T. Le Borgne¹, SL

3 Fayad¹, A. Gulmez¹, J. Allegret--Vautrot¹, F. de Chaumont⁴, C. Viollet¹, S. Didienne¹, N. Debray¹, J-P.

4 Hardelin¹, B. Girard², A. Mourot¹, J. Naudé³, F. Marti¹, B. Delord²*, P. Faure¹*

5

6 1. Brain Plasticity Laboratory, CNRS UMR 8249, ESPCI Paris, PSL Research University, Paris, France.

7 2. Institut des Systèmes Intelligents et de Robotique (ISIR), Sorbonne Université, CNRS, Paris, France.

8 3. Montpellier Inserm, CNRS, Université de Montpellier ; Institut de Génomique Fonctionnelle, Montpellier, France.

9

10 4. Human Genetics and Cognitive Functions, CNRS UMR UMR3571, Pasteur Institute, Paris, France.

11

12 *Corresponding authors. Emails: phfaure@gmail.com ; Bruno Delord: bruno.delord@gmail.com

13
14 **Summary:**

15 Social organization and division of labor are fundamental to animal societies, but how do these structures
16 emerge from individual interactions, and what role does neuromodulation play in shaping them? Using
17 behavioral tracking in a semi-natural environment, neural recordings, and computational models
18 integrating reinforcement and social learning, we show that groups of three isogenic mice spontaneously
19 develop specialized roles while solving a foraging task requiring individual decisions under social
20 constraints. Moreover, these roles are shaped by dopaminergic activity in the ventral tegmental area.

21 Strikingly, despite minor sex-differences in behavior when mice were tested alone, male triads formed
22 stable worker-scrounger relationships driven by competition, whereas female triads adopted uniform,
23 cooperative strategies. Model analysis revealed how intra- and inter-sex parameter differences in
24 resource exploitation, combined with contingent and dynamic social interactions, drive behavioral
25 specialization and labor division. Most notably, it highlighted how contingency, amplified by competition,
26 magnifies individual differences and shapes social profiles. The plastic, adaptive nature of social
27 organization within triads was confirmed by manipulating dopaminergic cell activity, which reshaped social
28 roles and altered group structure. Our findings support a feedback loop where social context shapes
29 neural states, which in turn reinforce behavioral specialization and stabilize social structures.

30

31
Social animals interact and coordinate their behaviors to form social organizations: structured patterns of relationships and interactions that stabilize into collective social forms (e.g., division of labor, norms, institutions) and specialized individual roles¹⁻⁵. A key challenge in behavioral science is deciphering how individual differences in behavior emerge, how they contribute to collective actions, and how they are linked to neural activity. While the broad effects of social organization on individual behavior and specialization are well documented ⁶⁻¹², the underlying cognitive and neurophysiological mechanisms remain largely unexplored. In particular, it is unclear how individual traits shape collective behaviors and social structures. This gap persists largely due to the difficulty of studying such mechanisms in controlled, artificial social settings that also allow for neurobiological investigation.

In social groups, the production of and access to shared resources drives strategies such as competition or cooperation⁷,¹³,¹⁴. This dynamic is exemplified by the producer-scrounger game, where individuals either produce resources (e.g., press a lever to obtain food) or scrounge from others who have already produced them¹³,¹⁵⁻¹⁹. Foraging strategy thus refers to the behavioral patterns and decision-making processes individuals adopt to obtain food, balancing effort, risk, and reward¹⁴. These strategies can range from independent food acquisition to social exploitation of others' efforts, shaped by ecological constraints and social dynamics. Evolutionary game theory²⁰ predicts stable outcomes under such interactions, such as the equilibrium between producer and scrounger strategies, depending on resource availability. However, this theory typically assumes behaviors are predetermined and stable²¹, overlooking the internal mechanisms that guide an individual's behavior, including the neural and cognitive processes which grant adaptation and learning-based flexible strategies²¹⁻²³. Here, we hypothesize that social foraging strategies, including producer-scrounger dynamics, emerge flexibly through socially-constrained reinforcement learning, with dopaminergic systems playing a central role in behavioral specialization. As social organization forms and norms develop, individuals are expected to adjust their actions based on expected payoffs, which affect both their neural activity and activity-dependent learning processes. In
turn, these learning-induced changes in neural activity and in cognitive decisions would shape how
individuals behave in social context²⁴.

Advances in tracking mice in semi-natural environments²⁵⁻³⁰ provide a unique opportunity to explore how social organization emerges from learning. By continuously monitoring individual behavior and interactions, we investigated how foraging specialization develops in small groups (n=3) housed in a controlled semi-natural setting. Males rapidly adopted distinct roles with divisions of labor, while females maintained uniform strategies. Modeling showed that learning and cognitive processes drive specialization, with exploration-exploitation trade-offs shaping social roles. Dopaminergic manipulations confirmed dopamine's key role in structuring social organization and reinforcing sex-specific adaptations.

In a lone context, mice exhibit two foraging strategies.

We first assessed mice foraging behavior and underlying neural mechanisms in a lone context. Female
or male mice (N=62) were placed alone in a 50 × 50 cm multi-compartment environment and tracked
continuously for 5 days and 4 nights using the Live Mouse Tracker (LMT²⁵) system (Figure 1a). This
setup includes a lever on one side and a food dispenser with a beam on the opposite side, enabling mice
to obtain pellets by pressing the lever and consume them at any time (Figure 1a). Mice displayed
nocturnal foraging activity, with lever presses (LP) and nose-pokes in the food dispenser occurring
primarily during the dark cycle, as expected (Figure 1b, left). Over time, animals increased LPs and
reduced nose-pokes, indicating learned association between LP and food retrieval (Figure 1b, right and
Figure S1a-c). Accordingly, trajectory durations between LP and food dispenser showed a skewed
distribution, with 54% below 6 seconds, which we defined as a complete sequence (Figure 1c, left). The
percentage of complete sequences (number of complete sequences / numbers of LPs) increased over
sessions, stabilizing at 35%, suggesting instrumental use of LPs (Figure 1c, right).

To explore the neural bases of this association, we used *in vivo* fiber photometry in DAT-Cre mice
expressing GCaMP7c in ventral tegmental area (VTA) dopaminergic (DA) neurons (Figure 1d and Figure
S1d). A phasic increase in VTA DA neuron activity was observed at the time of food delivery, regardless of the completeness of sequence types or of learning progression, consistent with reward signaling (Figure 1d). In contrast, we found a specific phasic increase in VTA DA neuron activity at LP time, but only after the learning phase and for complete sequences, consistent with prediction of food reward upon LP³¹. These findings further suggests that LPs associated with complete sequences acquired an instrumental status during learning, reflected in reward prediction-like VTA DA activity.

We next analyzed individual decision-making strategies. Clustering analysis, based on the number of LPs and the percentage of complete sequences, identified two distinct groups of mice: Achievers (N=31) and Storers (N=31) (Figure 1e). While these groups displayed a similar number of LPs, their percentages of complete sequences were different (Figure 1f, Figure Supp 1d). Specifically, Achievers made a high proportion of complete sequences, whereas Storers often made several LPs before retrieving the food, indicating different foraging strategies. Moreover, their distribution differed between sexes, with a higher proportion of Storers in female mice than in males (Figure 1g). Spontaneous VTA DA neuron activity – recorded post-experiment in anesthetized males– showed higher firing frequency in Achievers than in Storers (Figure 1h, Figure Supp 1e), suggesting a link between VTA DA activity and foraging strategies²⁶,³²,³³.

## 97 Distinct sex-specific strategies emerge within groups

We next explored the foraging strategies of WT mice (N=84) in a social context, by housing triads of either all-male or all-female task-naive mice for 8 days and 7 nights (Figure 2a, top). Overall, LP activity mirrored that in the lone context, with higher nocturnal activity (Figure 2a, bottom). Clustering, based on the percentage of LPs and complete sequences, identified four distinct foraging strategies (Figure 2b-c): i) Independents, who pressed the lever sufficiently for self-feeding (i.e., 150 pellets per day) and thus showed minimal social impact, ii) Workers, who pressed the lever extensively, producing more food than they consumed, effectively subsidizing other group members, iii) Scroungers, who rarely pressed the lever
and mainly relied on food produced by others, and (iv) *Storers*, who, as in the lone context, frequently pressed the lever but, contrarily to *Workers*, *Independents*, and *Scroungers*, performed few complete sequences. These roles or profiles rapidly emerged, with distinct behaviors observable from day 1 (Figure 1a) and maintained over time (Figure 1b). Indeed, within hours on day 1, *Workers* exhibited frequent LPs and complete sequences, while *Scroungers* displayed few LPs and complete sequences.

Similarly, *Storers* began pressing early (within the first five hours) and frequently, while performing few complete sequences. On the first day, *Storers* experienced frequent food “theft” (i.e., the food they produced was consumed by others), but also engaged in theft themselves (Figure 2c-d). These results indicate that role specialization and consistent strategies develop rapidly and are stable.

Sex-based differences were evident in both role types and their distribution within each cage. While males pressed the lever as often as females on average, they showed greater variability and more extreme behaviors (Figure 2d, Supp 2e). Females also performed significantly fewer complete sequences (Figure 2d, Supp 2e). Male triads (n = 14) typically comprised one or two *Scroungers* coexisting with *Workers* and *Independents*, but no *Storer* (Figure 2e-f), reflecting pronounced specialization. In contrast, female triads displayed remarkable uniformity in foraging strategies, with approximately 90% *Storers* (Figure 2e-f). This pronounced divergence between male and female triads suggests potential sex-based differences in neural, cognitive, and behavioral properties and/or social dynamics. Furthermore, the non-uniform distribution of profiles within triads (Figure 2f)—with certain triad compositions being absent, rare, or frequent—indicates that individual profiles are stabilized through dynamic social interactions. For instance, the emergence of three *Scroungers* appears unsustainable, suggesting that intrinsic social constraints govern the feasibility of specific role distributions.

To investigate the neural underpinnings of these behavioral strategies, we recorded VTA DA activity during the task, using fiber photometry in DAT-Cre mice (N=19) (Figure 2g). All groups showed increased VTA DA activity at food delivery, but *Workers* and *Independents* displayed a sharp rise during LP for complete sequences (Figure 2g), consistent with the use of LP as an instrumental action predicting
reward. In contrast, *Scroungers and Storers* lacked this response, with *Scroungers* instead displaying heightened DA activity when conspecifics pressed the lever (*Figure 2g*), indicating that for scroungers, LPs of their peers, rather than their own, acquired the status of a reward predictor. This suggests a social alertness mechanism in *Scroungers* for tracking food produced by others. These distinct dopaminergic signatures cross-validates the behavioral classification (*Figure 2b*) and strengthen the link between observable strategies and their neural substrates. Finally, post-experiment recordings of VTA DA neurons revealed sex- and strategy-dependent differences (*Figure 2h*, *Figure Supp 3*): males in social contexts exhibited increased VTA DA bursting activity compared to unexposed males, whereas females showed the opposite. *Workers and Independents* had higher DA bursting activity than *Scroungers and Storers*. These results underscore the dynamic and sex-dependent interplay between DA signaling and behavioral roles, suggesting plastic neural adaptations to social contexts.

## Sex-specific decision making in social context accounts for behavioral specialization

To assess the mechanisms of behavioral specialization, we developed a reinforcement Q-learning model (*Full model*, see methods). In this model, one or three e-mice learned to find food in an environment (*Figure 3a*) with four spatial states (1-4) and two action states — lever-press (L) and food dispenser (D)— where lever presses increase (and consumption decreases) reward availability by one pellet. Over time, state transition values formed a spatial gradient leading to the rewarded location (D; **Figure 3a**). Specialization could arise from intra- or inter-sex individual parameter differences, social constraints on behavior (i.e., lever press, food availability), or both. We assessed how model parameter variations could account for behavior specialization in lone and social conditions. Learning rate ($\alpha$) and temporal discount factor ($\gamma$) parameters had moderate effects, whereas inverse temperature ($\beta$), which controls the exploitation-exploration tradeoff (high-valued / random decisions), was critical (**Figure Supp 4a-c, e-g**). In the lone context, varying $\beta$ (*Figure 3a*) recreated the overall behavioral pattern observed in lone mice (*Figure 3b*). Lower $\beta$ favored exploration, leading to shallower value gradients (*Figure Supp 4c*), and
less constrained trajectories, resulting in lower percentage of complete sequences (Figure Supp 4b), resembling Storers (Figure 1f). By contrast, higher $\beta$ promoted steeper value gradients (Figure 3a) and higher percentage of complete sequences (Figure Supp 4b), akin to Achievers (Figure 1f). Crucially, introducing intra- and inter-sex-specific $\beta$ variability — where males exhibited higher $\beta$ values than females (Figure 3c) — allowed the model to replicate the observed sex differences (Figure 3b).

Strikingly, even with identical $\beta$, e-triads (e-mice triads) with high $\beta$s (favoring exploitation; $\beta \approx 1$), displayed distinct gradients and behaviors (Figure 3d), whereas e-triads with low $\beta$s (favoring exploration; $\beta \lesssim 1$) displayed no specialization (Figure Supp 4d). A narrow high-$\beta$ distribution produced large variability in percentage of complete sequences, with clustering unravelling Worker, Independent and Scrounger types, whereas (Figure 3e) a narrow low $\beta$ distribution produced Storer e-triads. Thus, sex-specific $\beta$ distributions produced patterns that closely matched those observed in male and female triads (Figure 3f). Moreover, these $\beta$ distributions accounted for the distribution of behavioral types and triads (Figure 3g) found in mice.

In contrast, contiguous low and high $\beta$ distributions with a larger variability led to specialized profiles without Storers (Figure 3h), as well as specific types and e-triads distributions (Figure Supp 4h, i), with high $\beta$ individuals as Scroungers or Independents and low $\beta$ as Independents or Workers (Figure 3h).

## Contingence and competition constrained behavioral specialization

We refined the causal mechanisms of social behavioral differentiation in a mathematically tractable reduced model (see methods) with two e-mice and two locations (L and D). Social differentiation arose through symmetry breaking of individual behaviors at a supercritical pitchfork bifurcation at $\beta=\beta_{Bifurc}$ value (Figure 4a). Below $\beta_{Bifurc}$ (e-females), the single fixed-point for Q-values toward D ($q_D$) and residence probability at D ($p_D$), lever press ($p_{LP}$) and complete sequence ($p_{PS}$) probabilities corresponded to Storer (Figure 4a and Figure Supp 5a-c). Above $\beta_{Bifurc}$ (e-males), the two stable fixed points branches, correspond to Workers and Scroungers (bistability; Figure 4a and Figure Supp 5a-c).
Qualitative analysis³⁴ also revealed the causal origin underlying symmetry breaking and behavioral bistability for e-male (Figure 4b-c and Figure Supp 5d-e). Initially, contingency (i.e. individual random choices) determined which e-mouse first accessed the dispenser (D), increased its Q-Value ($q_D$) and selection probability ($p_D$), which both were then amplified through the positive feedback loop between Q-learning and softmax decision processes. This feedback led the e-mouse to predominantly occupy the dispenser area and become a *Scrounger*. The reduced food access for the other e-mouse of the dyad implied low reward and learning at D, and consequently, a higher presence probability at L (i.e., more lever presses). The other e-mouse thus became a *Worker*. Thereby, competition for shared resources is the causal factor driving behavioral specialization in male e-triads. Higher $\beta$ intensified specialization, indicating that increased competition due to greater resource exploitation promoted differentiation in males. In contrast, e-female maintained uniform behaviors, attributed to lower levels of exploitation and resource competition.

Behavioral specialization thus arises from several intricate mechanisms. Sex-specific difference in $\beta$ determines whether e-mice become *Storers* (females) or specialize into distinct behavioral profiles (males). In males, specialization arises from both random contingency and higher $\beta$ values, with larger $\beta$ e-mice being more likely to become *Scroungers*. Therefore, the model highlights $\beta$ as a key determinant of behavioral specialization and social organization, and differences in exploitation between males and females as the critical factor underlying sex-specific social structures. Additionally, the model predicts that pairing high-$\beta$ with low-$\beta$ individuals should lead to a *Scrounger-Worker* dynamics, where the more exploitative individual becomes *Scrounger*, while the less exploitative one becomes the *Worker* (Figure 4d and 3h). Three experiments can empirically test these predictions. First, introducing a single high-$\beta$ male with a low-$\beta$ female pair should induce specialization in females, shifting them toward *Workers* or *Independents* (Figure 4e), and damper the competition between males, resulting in fewer *Scroungers*. Second, lowering $\beta$ in one male within a triad should increase its likelihood of becoming a *Worker* (Figure
205 **4f**). Finally, reducing *β* in all three males should produce a *Storer*-like dynamics, as seen in female triads
206 (**Figure 4g**). We next experimentally tested these hypotheses to validate our model's predictions.

207

208 **Testing model predictions with mixed-sex triads: The influence of low-β females on specialization**
209 We experimentally tested the first of these predictions (**Figure 4e**) by forming mixed-sex triads composed
210 of one male and two females (**Figure 5a**, Methods). Males performed fewer LPs than females, but unlike
211 in sex-matched triads, there was no difference in the percentage of complete sequences between sexes
212 (**Figure 2d**). Clustering analysis identified three profiles — *Scroungers*, *Workers* and *Independents*
213 (**Figure 5b**) — mirroring those seen in male-only triads. Females neither adopted the *Storer* profile (due
214 to increased competition) nor the *Scrounger* profile (due to their lower *βs*). Furthermore, in nearly half of
215 the triads (5/12 triads), no *Scrounger* mouse emerged, validating our first prediction: incorporating low-β
216 individuals (females) reduced the competition driving specialization, instead fostering more cooperative
217 dynamics. The emergence of an individual's role is thus not solely driven by its own traits, but by the
218 dynamics induced by trait distribution within the group.

219

220 **Artificially increasing VTA DA neuron activity influences role specialization in male microsocieties**
221 VTA DA activity has been linked to exploration-exploitation parameters, further supporting its involvement
222 in shaping role-specific behaviors ²⁴,³³,³⁵,³⁶. Accordingly, we observed variations in tonic VTA DA neuron
223 activity at the end of the experiment across behavioral roles and sexes. While our model does not explicitly
224 predict a direct relationship between β and VTA DA firing rate, we tested whether modulating DA activity
225 could impact specialization dynamics (Figure 4f, g). We selectively increased basal VTA DA neuron
226 activity in male DAT-Cre mice³⁷ using excitatory channelrhodopsin (ChR2) (Figure 5d, e, and Figure
227 **Supp 6a-d**). Mice received 15 minutes of photostimulation 24 hours before, and again 30 minutes before
228 the start of the microsociety experiment (Figure 5f). Electrophysiological recordings confirmed that this
229 stimulation reliably increased VTA DA neuron activity for at least 6 hours post-stimulation (Figure **Supp**
230

**6a-d).** This allowed us to assess whether altering one (Figure 4f) or all (Figure 4g) individual's DA activity affected role distribution.

231

232

We thus compared three types of male triads: (i) three YFP control mice (n=13 triads), (ii) one ChR2-stimulated mouse and two YFP controls (n=13 triads), and (iii) three ChR2-stimulated mice (n=12 triads).

233

234

Behavioral clustering identified four strategies – *Workers*, *Independents*, *Scroungers*, and *Storers* (Figure 5h) – whose behavior was consistent with our previous observations (Figure 5i-j; Figure 2b,c).

235

236

Surprisingly, *Storer* strategy emerged in photostimulated male triads, shifting the overall role distribution compared to YFP controls (Figure 5k; chi-squared test, p = 0.051). Importantly, the *Storer* was not necessarily the stimulated individual, emphasizing the role of group dynamics in shaping behavioral outcomes (Figure Supp 6e). This suggests that enhancing VTA DA activity prior to the social experiment leads to a lower β, facilitating the emergence of *Storers* by diminishing exploitative strategies. This effect was further confirmed in triads where all three mice were ChR2-stimulated, in which *Storer* frequency was increased (Figure 5k; chi-squared test, p = 0.0003). These findings demonstrate that while increasing VTA DA activity did not disrupt existing roles, it facilitated a female-like *Storer* profile in male triads. This highlights the complex interplay between VTA DA neurophysiology and the social environment in shaping individual strategies, providing a neural basis for the adaptive specialization observed within these microsocieties.

237

238

239

240

241

242

243

244

245

246

247

248

## Discussion

249

250

In this study, we highlight the importance of exploring social organization through the lens of division of labor to gain deeper insights into behavioral specialization in animal microsocieties. Traditional accounts of behavioral specialization depict fixed social roles as evolutionarily stable strategies²⁰,³⁸, but our findings challenge this static view by instead revealing a dynamic, learning-driven processes ²¹,²³. Our results demonstrate that social norms, conceptualized as implicit rules governing group behavior, emerge rapidly from animals' coexistence and dynamical structure, guiding both individual roles and collective dynamics.

251

252

253

254
A key contribution here is the demonstration that an individual cognitive property, i.e. the exploration-exploitation tradeoff, constraints the emergent social structure. Using reinforcement-learning models and experimental data, we show that variability in neural mechanisms, particularly dopaminergic signaling in the VTA, shapes how behavioral roles evolve. This neurophysiological foundation emphasizes that behavioral specialization is not merely an adaptive response to external social demands, but arises from a permanent causal interaction between intrinsic neural properties and social interactions, through a feedback loop: social contexts constraints individual's neural dopaminergic levels, which in turn reinforce specific behavioral strategies, eventually stabilizing social structures over time.

A striking aspect of our findings is the pronounced sex-specific divergence in social organization. Male groups exhibited distinct and stable divisions of labor, including roles such as *Workers*, *Scroungers* and *Independents*, whereas female groups displayed remarkable behavioral uniformity, predominantly behaving as *Storers*. These differences highlight how sex-based cognitive and neurophysiological differences interact with social environments to shape collective dynamics. For example, male-specific dopaminergic activity patterns, coupled with higher exploitation behavior (high $\beta$), promote competition and role specialization. Conversely, females' explorations tendency (low $\beta$) appears minimize competition and favor behavioral uniformity. These findings resonate with observations that males and females exhibit distinct responses to social competition and resource access <sup>39</sup>. Males are more prone to develop divergent roles under competitive pressures, while females tend to maintain consistent strategies, highlighting sex-specific adaptations in response to social dynamics.

The role of early-life contingency in shaping behavioral and social outcomes warrants discussion. It has been proposed that early-life experiences and chance events can magnify small initial differences between individuals, leading to diverging trajectories<sup>24,29,39-41</sup>. In competitive social environments, such as those often encountered by male mice, these feedback loops amplify disparities and drive the emergence of distinct behavioral roles. By contrast, female mice may experience reduced competitive pressures, allowing for more equitable outcomes within their groups. In the model, $\beta$ distributions required
280 to reproduce experimental behaviors were compatible with the hypothesis that the presence of other
281 males or other females in the unisex triads increases or decreases exploitation, respectively, possibly due
282 to social effects on the internal state of the individual 3,42,43. This contingency-based framework
283 emphasizes the importance of considering both stochastic aspects (e.g., unpredictable group composition
284 or initial random reward discovery) and deterministic factors (e.g., one animal becoming a *Scrounger*
285 favors others to become *Workers*) in shaping sex-specific behaviors and social structures.

286 Our findings demonstrate that dopaminergic signaling plays a central role in behavioral specialization,
287 with *Workers* and *Independents* showing VTA activity linked to reward prediction, while *Scroungers*
288 respond more to conspecific lever presses, suggesting social alertness. These distinct dopaminergic
289 profiles reinforce the idea that neural plasticity underlies the emergence of specialized roles within social
290 groups. Furthermore, variations in tonic VTA DA activity, shaped by both intrinsic predispositions and
291 social experiences, likely contribute to the stabilization of individual strategies over time. Future research
292 should determine whether these neural adaptations persist across successive different social settings or
293 revert when individuals are removed from their original groups.

294 In conclusion, this work provides a comprehensive analysis of how individual and social factors interact
295 to shape behavioral specialization in animal microsocieties. By bridging neurophysiology, behavioral
296 science, and computational modeling, our findings offer a nuanced understanding of the dynamical
297 interplay between individual cognition and collective behavior, with broad implications for the study of
298 social systems across taxa.

299

300 References

301

302 1. Wilson, E. O. Sociobiology: The New Synthesis. (Belknap Press, 1975).

303 2. Beshers, S. N. & Fewell, J. H. Models of division of labor in social insects. Entomology 46, 413-440

304 (2001).

305 3. Bergmüller, R. & Taborsky, M. Animal personality due to social niche specialisation. Trends in

306 ecology & evolution 25, 504-511 (2010).

307 4. Gartland, L. A., Firth, J. A., Laskowski, K. L., Jeanson, R. & Ioannou, C. C. Sociability as a personality

308 trait in animals: methods, causes and consequences. Biol Rev (2021) doi:10.1111/brv.12823.
309 5. Sumpter, D. J. T. The principles of collective animal behaviour. Philos. Trans. R. Soc. B: Biol. Sci.

310 361, 5–22 (2006).

311 6. Haun, D. B. M., Leeuwen, E. J. C. van & Edelson, M. G. Majority influence in children and other

312 animals. Dev. Cogn. Neurosci. 3, 61–71 (2013).

313 7. Ioannou, C. C. & Laskowski, K. L. Conformity and differentiation are two sides of the same coin.

314 Trends Ecol. Evol. 38, 545–553 (2023).

315 8. Cooper, G. A. & West, S. A. Division of labour and the evolution of extreme specialization. Nat Ecol

316 Evol 2, 1161–1167 (2018).

317 9. Webster, M. M. & Ward, A. J. W. Personality and social context. Biol Rev 86, 759–773 (2011).

318 10. Loftus, J. C., Perez, A. A. & Sih, A. Task syndromes: linking personality and task allocation in social

319 animal groups. Behav Ecol 32, araa083 (2021).

320 11. Spiegel, O., Leu, S. T., Bull, C. M. & Sih, A. What’s your move? Movement as a link between

321 personality and spatial dynamics in animal populations. Ecol Lett 20, 3–18 (2017).

322 12. Jolles, J. W., King, A. J. & Killen, S. S. The Role of Individual Heterogeneity in Collective Animal

323 Behaviour. Trends in ecology & evolution 35, 1–14 (2019).

324 13. Barnard, C. J. & Sibly, R. Producers and scroungers: a general model and its application to captive

325 flocks of house sparrows. Animal Behaviour 29, 543–550 (1981).

326 14. Giraldeau, L.-A. & Dubois, F. Chapter 2 Social Foraging and the Study of Exploitative Behavior. Adv

327 Stud Behav 38, 59–104 (2008).

328 15. Oldfield-Box, H. Social Organization of Rats in a “Social Problem” Situation. Nature 213, 533–534

329 (1967).

330 16. Ahn, S. et al. A role of anterior cingulate cortex in the emergence of worker–parasite relationship.

331 Proc National Acad Sci 118, e2111145118 (2021).

332 17. Colin, C. & Desor, D. Differenciations comportementales dans des groupes de rats soumis a une

333 difficile d’acces a la nourriture. Behavioural processes (1986).

334 18. Grasmuck, V. & Desor, D. Behavioural differentiation of rats confronted to a complex diving-for-food

335 situation. Behav Process 58, 67–77 (2002).

336 19. Alfaro, L. & Cabrera, R. Effect of group size on producer–scrounger strategies of Wistar rats. Behav

337 Process 182, 104280 (2021).

338 20. Smith, J. M. Evolution and the Theory of Games. (Cambridge University Press, 1982).

339 21. Fawcett, T. W., Hamblin, S. & Giraldeau, L.-A. Exposing the behavioral gambit: the evolution of

340 learning and decision rules. Behav. Ecol. 24, 2–11 (2013).

341 22. Laland, K. N., Sterelny, K., Odling-Smee, J., Hoppitt, W. & Uller, T. Cause and Effect in Biology

342 Revisited: Is Mayr’s Proximate-Ultimate Dichotomy Still Useful? Science 334, 1512–1516 (2011).

343 23. Dridi, S. & Lehmann, L. On learning dynamics underlying the evolution of learning rules. Theor.

344 Popul. Biol. 91, 20–36 (2014).

345 24. Faure, P., Fayad, S. L., Solié, C. & Reynolds, L. M. Social Determinants of Inter-Individual Variability

346 and Vulnerability: The Role of Dopamine. Front Behav Neurosci 16, 836343 (2022).

347 25. Chaumont, F. de et al. Real-time analysis of the behaviour of groups of mice via a depth-sensing

348 camera and machine learning. Nature biomedical engineering 7, 462–942 (2019).

349 26. Torquet, N. et al. Social interactions impact on the dopaminergic system and drive individuality.

350 Nature Communications 9, 3081 (2018).

351 27. Weissbrod, A. et al. Automated long-term tracking and social behavioural phenotyping of animal

352 colonies within a semi-natural environment. Nature Communications 4, (2013).
353 28. Zilkha, N. et al. Sex-dependent control of pheromones on social organization within groups of wild house mice. Curr Biol (2023) doi:10.1016/j.cub.2023.02.039.

354 29. Freund, J. et al. Emergence of Individuality in Genetically Identical Mice. Science (New York, N.Y.) 340, 756–759 (2013).

355 30. Shemesh, Y. et al. High-order social interactions in groups of mice. eLife 2, e00759–e00759 (2013).

356 31. Schultz, W., Dayan, P. & Montague, P. R. A Neural Substrate of Prediction and Reward. Science 275, 1593–1599 (1997).

360 32. Fayad, S. L. et al. Individualistic reward-seeking strategies that predict response to nicotine emerge among isogenic male mice living in a micro-society. PLOS Biol. 22, e3002850 (2024).

362 33. Dongelmans, M. et al. Chronic nicotine increases midbrain dopamine neuron activity and biases individual strategies towards reduced exploration in mice. Nat Commun 12, 6945 (2021).

364 34. Strogatz, S. H. Nonlinear Dynamics and Chaos : With Applications to Physics, Biology, Chemistry, and Engineering. (1994).

366 35. Beeler, J. A., Daw, N., Frazier, C. R. M. & Zhuang, X. Tonic dopamine modulates exploitation of reward learning. Frontiers in Behavioral Neuroscience 4, 170 (2010).

368 36. Cinotti, F. et al. Dopamine blockade impairs the exploration-exploitation trade-off in rats. Scientific Reports 9, 6770 (2019).

370 37. Bariselli, S. et al. Role of VTA dopamine neurons and neuroligin 3 in sociability traits related to nonfamiliar conspecific interaction. Nat. Commun. 9, 3173 (2018).

372 38. Laland, K. N. & O'Brien, M. J. Cultural Niche Construction: An Introduction. Biol. Theory 6, 191–202 (2011).

374 39. Zipple, M. N. et al. Competitive social feedback amplifies the role of early life contingency in male mice. Science 387, 81–85 (2025).

376 40. Honegger, K. & Bivort, B. de. Stochasticity, individuality and behavior. Current biology : CB 28, R8–R12 (2018).

378 41. Wolf, M., Bierbach, D. & Laskowski, K. L. Behavioural individuality in clonal fish arises despite near-identical rearing conditions. Nature Communications 8, 1–7 (2017).

380 42. Battivelli, D. et al. Dopamine neuron activity and stress signalling as links between social hierarchy and psychopathology vulnerability. Biol. Psychiatry (2023) doi:10.1016/j.biopsych.2023.08.029.

382 43. Wei, A. et al. Sexually dimorphic dopaminergic circuits determine sex preference. Science 387, eadq7001 (2025).

384 44. Turiault, M. et al. Analysis of dopamine transporter gene expression pattern -- generation of DAT-iCre transgenic mice. Febs J 274, 3568 3577 (2007).

386 45. Ungless, M. A. & Grace, A. A. Are you or aren't you? Challenges associated with physiologically identifying dopamine neurons. Trends Neurosci 35, 422 430 (2012).

388 46. Pinault, D. A novel single-cell staining procedure performed in vivo under electrophysiological control: morpho-functional features of juxtacellularly labeled thalamic cells and other central neurons with biocytin or Neurobiotin. J Neurosci Meth 65, 113–136 (1996).

391

392
<figure><img src="image_2.png" /><figcaption>b</figcaption></figure>

<figure><img src="image_3.png" /><figcaption>c</figcaption></figure>

<figure><img src="image_4.png" /><figcaption>d</figcaption></figure>

<figure><img src="image_5.png" /></figure>

<figure><img src="image_6.png" /><figcaption>e</figcaption></figure>

<figure><img src="image_7.png" /><figcaption>f</figcaption></figure>

<figure><img src="image_8.png" /><figcaption>g</figcaption></figure>

<figure><img src="image_9.png" /><figcaption>h</figcaption></figure>

<figure><img src="image_10.png" /></figure>

Figure 1
393 **Figure 1: Different reward-seeking strategies emerged in a non-social semi-natural environment.**

394 **(a)** Mice were housed individually for 5 days in habitats where food was delivered via a lever-press (LP),
395 positioned opposite to the food dispenser. **(b)** (Left) Cumulative polar LP histogram over 24 hours,
396 showing peak activity during the dark cycle. (Right) Time course of LP and beam breaks over 4 nights.

397 **(c)** (Left) Distribution of trajectory durations between LP and dispenser, with a 6 sec threshold
398 distinguishing complete and incomplete (over 6 sec) sequences. (Right) % of complete sequences over
399 time. **(d)** (Top) Schema of injections and implantations. (Bottom) VTA DA activity (fiber photometry) differs
400 before and after learning (after day 3) and based on sequence completion, with activation at food delivery.

401 **(e)** (Top) Clustering distinguished two food seeking strategies: “Achievers” and “Storers”. (Bottom)
402 Examples of 6 sec trajectory after a LP, from two different mice. **(f)** Percentage of complete sequences
403 over time by strategy. **(g)** Strategy distribution by sex. **(h)** (Left) In-vivo electrophysiology setup. (Right)
404 VTA DA activity in anesthetized males at experiment end, depending on strategy.

405

406
bioRxiv preprint doi: https://doi.org/10.1101/2025.02.10.637303; this version posted February 11, 2025. The copyright holder for this preprint (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made available under a CC-BY 4.0 International license.

<caption>a</caption>

<figure><img src="image_11.png" /><figcaption>Food Dispenser</figcaption></figure>

Independent (N=19)

<figure><img src="image_12.png" /><figcaption>b</figcaption></figure>

Worker (N=9)

<figure><img src="image_13.png" /></figure>

Scrounger (N=18)

<figure><img src="image_14.png" /></figure>

Storer (N=38)

<figure><img src="image_15.png" /></figure>

<figure><img src="image_16.png" /></figure>

<figure><img src="image_17.png" /></figure>

<figure><img src="image_18.png" /><figcaption>d</figcaption></figure>

N=42; N=42

<figure><img src="image_19.png" /><figcaption>e</figcaption></figure>

<figure><img src="image_20.png" /><figcaption>f</figcaption></figure>

<figure><img src="image_21.png" /><figcaption>g</figcaption></figure>

1 AAV1-Syn-FLEX-GCamp7c

2 Optic Fiber

<figure><img src="image_22.png" /></figure>

1 2

<figure><img src="image_23.png" /></figure>

<figure><img src="image_24.png" /><figcaption>h</figcaption></figure>

<figure><img src="image_25.png" /></figure>

<figure><img src="image_26.png" /></figure>

In vivo anesthetized recording
VTA DA neurons (REC)

Naive ♂: 100 neurons | 6 mice

Post-Task ♂: 335 n|22 m

Naive ♀: 75 n|7 m

Post-Task ♀: 131 n|9 m

Independents: 191 n|13 m

Workers: 79 n|6 m

Scroungers: 96 n|6 m

Storers: 100 n|7 m

<figure><img src="image_27.png" /></figure>

<figure><img src="image_28.png" /></figure>

<figure><img src="image_29.png" /><figcaption>Figure 2</figcaption></figure>
**Figure 2: Distinct sex-specific strategies emerge within groups (a) (Top) Triads of mice in the same environment as figure 1. (Bottom) Polar histogram of cumulative LP over 24 hours. (b) (Left) Clustering analysis identified four strategies Independents, Workers, Scroungers, and Storers. (Right) Example of a 6 sec trajectory after LP for each strategy, except for Scroungers, for which trajectory after Worker LPs is shown. (c) Time course of LPs (top left) and % of complete sequences (bottom left), with strategy-averaged values for the last two days (top right, bottom right). (d) Mean number of LPs and % of complete sequences by sex. (e) Strategy distribution by sex. (f) Triad composition (g) (Top) Schema of injection and implantation for fiber-photometry recording and example image. (Bottom) VTA DA activity varies at Food, LP, and conspecific LP depending on strategy. (h) (Top) In vivo electrophysiology setup and protocol. (Bottom) Firing patterns of VTA DA cells (percentage of spikes within burst - %SWB) in anaesthetized mice: comparison between naïve and experienced mice (left) and across strategies (right).**
<figure><img src="image_30.png" /><figcaption>b</figcaption></figure>

<figure><img src="image_31.png" /><figcaption>c</figcaption></figure>

<figure><img src="image_32.png" /><figcaption>d</figcaption></figure>

<figure><img src="image_33.png" /></figure>

Three e-mice
in a cage

<figure><img src="image_34.png" /><figcaption>e</figcaption></figure>

<figure><img src="image_35.png" /><figcaption>f</figcaption></figure>

<figure><img src="image_36.png" /><figcaption>g</figcaption></figure>

<figure><img src="image_37.png" /><figcaption>h</figcaption></figure>

<figure><img src="image_38.png" /><figcaption>Figure 3</figcaption></figure>
**Figure 3 Reinforcement learning model and emergent strategies.** (a) (Left) RL model (full model, methods) with one e-mouse, four states and two actions: lever press (L) and eat at dispenser (D). Q-values represent agent learning (e-mice). $\beta$=2. (Right) LP count vs. % of complete sequences for 200 e-mice with $\beta$ uniformly distributed in [0.5, 2.5]. (b) Experimental LP vs. % of complete sequences for 62 mice, with males (empty squares) and females (filled circles). (c) Mean proportions of *Achievers* and *Storers* over repeated clustering (see methods) (d) Simulation as in (a), with three e-mice. (e) LP vs. % of complete sequences for 200 e-triads, clustered into *Independents* (red), *Workers* (black), *Scroungers* (blue), and *Storers* (yellow). (f) Experimental data for 84 mice, with color code for profile as in (e). (g) Strategy distribution and triad composition for 200 single-sex e-triads. (h) LP vs. % of complete sequences for 200 mixed e-triads (two females, one male) clustered into *Independents* (red), *Workers* (black), and *Scroungers* (blue).

430

431
<figure><img src="image_39.png" /><figcaption>a</figcaption></figure>

<figure><img src="image_40.png" /><figcaption>b</figcaption></figure>

<figure><img src="image_41.png" /><figcaption>d</figcaption></figure>

<figure><img src="image_42.png" /><figcaption>e</figcaption></figure>

<figure><img src="image_43.png" /><figcaption>g</figcaption></figure>

<figure><img src="image_44.png" /></figure>

<figure><img src="image_45.png" /><figcaption>Figure 4</figcaption></figure>
**Figure 4: Behavioral specialization from competition-induced bifurcation.** (a) Reduced model with two states (L: lever, D: dispenser) and two mice (1 and 2). (Right) Bifurcation diagram of the Q-value at D as a function of $\beta$, with two characteristic values, $\beta_{bifurc}$ a bifurcation point and $\beta_{NoFood}$ (see methods). (b) Similar low $\beta$s e-mice results in a single stable fixed point (uniform profiles) in the velocity landscape. Black curve: simulated dynamics example. (c) Similar high $\beta$ values yields two stable fixed points corresponding to a *Worker* (here *e-mouse #1*) and a *Scrounger*. (d) One high-$\beta$ e-male and one low-$\beta$ e-female yields a stable fixed point, with a male *Scrounger*. (e-g) Experimental design to test model predictions. (e) A triad with a high-$\beta$ male and a low-$\beta$ female pair to assess role specialization. (f) Lowering $\beta$ in a male to test the likelihood of adopting *Worker* roles. (g) Reducing $\beta$ in a male triad to evaluate the emergence of *Storer*-like dynamics analogous to female triads.

442

443
<figure><img src="image_46.png" /><figcaption>b</figcaption></figure>

<figure><img src="image_47.png" /></figure>

<figure><img src="image_48.png" /><figcaption>d DAT-Cre</figcaption></figure>

<figure><img src="image_49.png" /></figure>

1 AAV9-Ef1α-DIO-ChR2-YFP
or AAV5-Ef1α-DIO-YFP

2 Optic fiber

<figure><img src="image_50.png" /></figure>

<figure><img src="image_51.png" /></figure>

<figure><img src="image_52.png" /><figcaption>h p<10⁻⁴</figcaption></figure>

<figure><img src="image_53.png" /></figure>

<figure><img src="image_54.png" /></figure>

<figure><img src="image_55.png" /><figcaption>Figure 5</figcaption></figure>
Figure 5: Testing model predictions and exploring the role of VTA dopaminergic activity in behavioral specialization. (a) (Top) Mixed-sex triads (1 male and 2 females). (Down) Mean LP and % of complete sequences for males (green N=12) and females (blue N=24). (b) (Left) Clustering identified three strategies: *Scroungers*, *Workers*, and *Independents*. (Right) Number of LPs and % of complete sequences for the last two days per strategy. (c) (Top) Strategy distribution by sex. (Bottom) Triads composition. (d) Experimental design for VTA DA manipulation using ChR2 in male DAT-Cre mice. (e) Example image of viral infection and fiber implantation. (f) Experimental setup: VTA DA stimulation applied 24 hours and 30 minutes before microsociety experiment. (g) Triad compositions: (i) all YFP controls, (ii) two YFP with one ChR2-stimulated mouse, (iii) all ChR2-stimulated. (h) Clustering in optogenetically manipulated male triads identified four strategies: *Workers*, *Independents*, *Scroungers*, and *Storers*. (i) Number of LPs (top) and % of complete sequences (bottom) over time by strategy. (j) Mean LP (top) and % of complete sequences (bottom) for the last 2 days. (k) Strategy distributions across triad types: (i) all YFP control mice, (ii) one ChR2-stimulated mouse with two YFP controls, and (iii) all ChR2-stimulated mice. (l) Triad composition.

458

459
460 Methods

461

462 Animals

463 Males and females wild-type (WT) C57BL/6J mice (8-12 weeks old, Janvier Labs, France) or DAT-iCre mice (C57BL/6J background, BAC-DAT-iCre, from François Tronche⁴⁴) were housed in groups in an animal facility and maintained on a 12-hour light-dark cycle (7:00 a.m.–7:00 p.m.), at an average temperature of 22°C (22 ± 2°C) and ~50% humidity. All the animals lived in triads in the animal facility before entering in the task for 8 days with the same triad for the same-sex experiments. In lone condition, the mice were housed individually in the cage for 5 days only, to prevent side effect of long-term social isolation.

470 All experiments and procedures were performed in accordance with European Commission directives 219/1990, 220/1990 and 2010/63, and approved by the ESPCI and the ethical committee #059 under APAFIS #34335-2021121318085835.

473

474 Behavior

475 Mice were placed either individually or in groups of three (triads) in a 50 × 50 cm large environment and continuously tracked using the LMT system²⁵, allowing to keep the identity of each animal over extended periods. The mice in groups (3 males, 3 females or 2 females / 1 male) stayed for 8 days and 7 nights in the environment, while mice placed alone in the cage stayed only 5 days and 4 nights.

479 The cage was composed of different zones freely accessible by all the mice at any moment. A lever was placed on one side of the environment, while a food dispenser with a magazine were located at the opposite side. Every time the lever is pressed, a small food pellet of 20 mg (TestDiet 5TUL/1811142 purified, Bio-Concept) is delivered at the magazine and the lever becomes inactive for 5 seconds. A nose-poke in the magazine leads to a beam break that is a proxy of the animal's consumption of the pellet. Outside the 5 seconds delay between each lever-press, all the mice can press at any moment and
consume the pellet at any time following its release. The mice were checked and weighted every day to
keep track of the good health of each individual. No mouse was removed from the task either for weight
loss or for aggression.

488

489 **Live Mouse-Tracker system**

490 This setup, in our paradigm, provided a controlled context for examining how individual mice learn and
491 adapt to the task, either in isolation or within a standard social setting, allowing comparisons between
492 lone and group behaviors. We employ comprehensive 24/7 data collection within the experimental
493 environment, a method that provides rich longitudinal datasets. The data collected fall into two main
494 categories: discrete events and continuous video streams. Discrete events are automatically captured by
495 sensors integrated into the environment, including infrared beams at the lever and nose poke, and food
496 dispenser activations (MedAssociates). Each event is tagged with the corresponding sensor number,
497 timestamp, and the identity of the animal involved. These events are logged into a structured database,
498 enabling detailed analysis of individual and group-level behavioral patterns over time. In addition to
499 discrete events, we record continuous 24/7 video streams, which provide a complete visual record of
500 animal's behavior within the environment. To manage the substantial storage demands of such video
501 data, we implement real-time processing techniques to extract essential features. This includes tracking
502 each animal's center of mass and body contours, allowing us to quantify locomotion, spatial organization,
503 and interactions. These extracted metrics are stored in a compressed format for subsequent analysis,
504 ensuring that critical information is retained without overwhelming storage capacity. This dual approach
505 to data collection provides a nuanced understanding of behavior, capturing both fine-grained event data
506 and broader spatial and social dynamics.

507 Every day the LMT and MedAssociates softwares were relaunched to prevent huge file weight.

508 The number of lever-presses, the number of nosepokes and the number of complete sequences were
509 extracted directly from the database using PyCharm (Python) and MySQL. We defined two specific zones
around the lever and the magazine to be able to know which animal was located there during a lever press or a nose poke. We then determined which TTL was associated with which animal in the zone. For each day of analysis, it was not possible to associate a little percentage of trials (< 5%) with specific mice since it happened that several mice were present at the same time either at the lever or at the magazine.

Once the association made between the position of the animal at a specific timestamp and the time of the TTL, we obtained the number of lever-presses and nose pokes for each animal. For complete sequences, we quantified the number of times that the mouse which had pressed the lever was the first arrived at the subsequent nose poke in less than 6 seconds, otherwise the sequence was not counted (typically if a mouse was already at the magazine waiting for the food for example).

## Behavioral analyses

After extraction of the data from the database using PyCharm and MySQL, MATLAB was used to perform more detailed analysis. The percentage (%) of complete sequences was calculated as the number of complete sequences divided by the number of lever presses.

The clustering analysis was made using "linkage" followed by "dendrogram" functions on MATLAB. We normalized all the data for all the mice as follows:

- For the lone context: we used % of complete sequences and the ratio lever-press on nose poke for each mouse. We then pooled all the mice of all the cages, males and females, and we performed the clustering on the last two days of the experiments.

- For the social context: we used % of complete sequences and % of lever presses in the triad (number of lever presses of one mouse divided by the total number of lever presses of the triad). We pooled all the mice together, males and females for figure 2, and all the virus conditions for figure 5, independently of the triad. We performed the clustering analysis on the last five days.
535 *In vivo electrophysiology*

536 Mice were deeply anaesthetized with isoflurane delivered continuously (3% induction, 1-2% maintenance;

537 TEMSega). The scalp was opened and a hole was drilled in the skull above the location of the VTA.

538 Intravenous administration of saline or nicotine (30µg/kg) was carried out through a catheter (30G needle

539 connected to polyethylene tubing PE10) connected to a Hamilton syringe, into the saphenous vein of the

540 animal. Extracellular recording electrodes were constructed from 1.5 mm outer diameter / 1.17 mm inner

541 diameter borosilicate glass tubing (Harvard Apparatus) using a vertical electrode puller (Narishige). The

542 tip was broken straight and cleaned under microscopic control to obtain a diameter of about 1 µm. The

543 electrodes were filled with a 0.5% NaCl solution containing 1.5% of neurobiotin® tracer (VECTOR

544 laboratories) yielding impedances of 6-9 MΩ. Electrical signals were amplified by a high-impedance

545 amplifier (Axon Instruments) and monitored audibly through an audio monitor (A.M. Systems Inc.). The

546 signal was digitized, sampled at 25 kHz, and recorded on a computer using Spike2 software (Cambridge

547 Electronic Design) for later analysis. The electrophysiological activity was sampled in the central region

548 of the VTA (coordinates: between 3.1 to 4 mm posterior to bregma, 0.3 to 0.7 mm lateral to midline, and

549 4 to 4.8 mm below brain surface). Individual electrode tracks were separated from one another by at least

550 0.1 mm in the horizontal plane. Spontaneously active DA neurons were identified based on previously

551 established electrophysiological criteria⁴⁵. After recording, a subset of nicotine-responsive cells were

552 labelled by electroporation of their membrane: successive currents squares were applied until the

553 membrane breakage, to fill the cell soma with neurobiotin contained into the glass pipet⁴⁶. To be able to

554 establish correspondence between neurons responses and their location in the VTA, we labeled one type

555 of response per mouse: solely activated neurons or solely inhibited neurons, with a limited number of cells

556 per brain (1 to 4 neurons maximum, 2 per hemisphere), always with the same concern of localization of

557 neurons in the VTA.

558
559 Stereotaxic surgeries

560 For fiber photometry experiments: injection of AAV1-Syn-FLEX-GCamp7c (Titer ≥ 9.4x10¹² vg.mL⁻¹, Addgene, Watertown, MA, USA) was performed in 6-8 week-old DAT-Cre mice. Mice were anesthetized with a mixture of oxygen (1L/min) and isoflurane 2% (TEMSega, Pessac, France). After shaving the skin, mice were placed in a stereotaxic frame (David Kopf, Tujunga, USA), the skin was disinfected and locally anesthetized with 70 µL of lurocaine (0.5 %). Craniotomies were then performed, and unilateral injections were performed over the VTA at AP: -3.20 mm ML: ± 0.5 mm from the bregma and DV: -4.20 ± 0.05 mm from the brain surface. Injections were carried out with a glass micropipette (Drummond Scientific Company, Broomall, PA), at a 100 nL/min rate for a total volume of 300 nL in the VTA. At the end of the surgery, 100 µL of buprenorphine (0.1mg./kg) was injected subcutaneously. Unilateral optic fiber (200 µm core, NA= 0.39, ThorLabs) implantations were performed 2 – 3 weeks later above VTA with the same coordinates, and dental acrylic (SuperBond, Sun medical) was used for fixing the ferrule to the skull.

571

572 For optogenetic experiments: AAV5-Ef1α-DIO-ChR2(H134R)-EYFP (Titer ≥ 4.7x10¹² vg.mL⁻¹, VVF Zürich, Switzerland) or AAV5-Ef1α-DIO-EYFP (Titer ≥ 5.4x10¹² vg.mL⁻¹, VVF Zürich, Switzerland) were injected in 6-8 week-old DAT-Cre mice. Mice were anesthetized and disinfected as previously described.
573
574
575 The animals were placed in a stereotaxic frame and bilateral injections were performed in the VTA at the same coordinates than above, 300 nL per side, using a glass micropipette. After 2 – 3 weeks, one optic fiber (200 µm core, NA= 0.39, ThorLabs) was implanted unilaterally above the VTA with 10° angle (AP: -3.20 mm ML: ±0.9 mm from the bregma and DV: -3.95 mm from the brain surface), and fixed to the skull using dental acrylic.
576
577
578
579

580 After all the surgeries, mice were placed alone in a heated cage for one hour recovery, and were subsequently checked every day. The behavioral experiments started at least 1 week after the surgery.
581
582
583
Injections and implantations were all confirmed with post-hoc immunochemistry.
584 **Immunohistochemistry**

585 Following euthanasia, brains were rapidly extracted and fixated in a solution of 4% paraformaldehyde for
586 a minimum of three days at 4°C. Subsequently, serial sections of 60 micrometers in thickness were
587 obtained from the midbrain using a vibratome.

588 For immunostaining experiments (TH-staining), free-floating brain sections were initially incubated for one
589 hour at 4°C in a phosphate-buffered saline (PBS) blocking solution containing 3% bovine serum albumin
590 (BSA) by volume/volume ratio (Sigma; A4503) and supplemented with 0.2% Triton X-100. This was
591 followed by overnight incubation at 4°C with a primary antibody solution containing a mouse monoclonal
592 anti-tyrosine hydroxylase antibody (anti-TH, Sigma, T1299) diluted 1:500 in PBS supplemented with 1.5%
593 BSA and 0.2% Triton X-100.

594 The following day, sections were thoroughly rinsed with PBS and then incubated for three hours at room
595 temperature (22-25°C) with the secondary antibody Cy3-conjugated goat anti-mouse (Jackson
596 ImmunoResearch, 715-165-150) diluted at 1:500 in a solution of 1.5% BSA in PBS. After undergoing three
597 additional rinses with PBS, the slices were mounted using Prolong Gold Antifade Reagent with DAPI
598 (Invitrogen, P36930). Microscopic analysis was performed using a Zeiss epi-fluorescent microscope, and
599 image capture was performed by a camera. To preserve the raw data, images were captured in grayscale
600 format using the Zeiss software (ZEN microscopy software). Following acquisition, false coloring was
601 applied using ImageJ to enhance visualization of the immunofluorescence signal.

602

603 **Fiber photometry**

604 Mice injected with AAV1-Syn-FLEX-GCamp7c and implanted in the VTA underwent fiber photometry
605 experiments in the lone or social context. A Doric Lenses fiber photometry system was employed to record
606 fluorescence signals reflecting dopamine neuron activity in the VTA. The system utilized a fiber
607 photometry console linked to an LED driver, controlling a connectorized blue light-emitting LED (465 nm)
608 in Lock-in mode at 220.537 Hz. This signal was then transmitted through an optic patch cord to the Mini
609 Cube (FMC4_AE(405)_E(460-490)_F(500-550)_S). An optical fiber, connected to the Mini Cube's sample
610 port and the animal's implanted fiber via a zirconia sleeve, served as the conduit for both light stimulation
611 and recorded fluorescence. The received light signal was converted to electrical signals by a
612 photoreceiver using the AC Low setting, before being transmitted through another optic patch cord to the
613 Mini Cube via a dedicated fiber optic adaptor. Finally, the Doric Neuroscience Studio software acquired
614 the signal at a sampling rate of 12.0 kHz and applied a low-pass filter with a 12.0 Hz cutoff frequency.

615 For lone condition, the animals were recorded between one and two hours at the beginning of the dark
616 cycle, when the mice were active, during the first day and the last two last of the experiment.

617 For social condition, the animals were also recorded at the beginning of the dark-cycle, one after the
618 other, between one and two hours each.

619
620 Analysis of fiber photometry recordings

621 We first subtracted a biexponential detrend to our fluorescence signal in order to eliminate the slow decay
622 of the signal due to sensor photobleaching, and then added on offset equal to the signal mean before
623 detrending. The ΔF/F ratio was calculated as follow:

$$
624 \frac{\text{fluo signal} - \text{baseline fluo signal}}{\text{baseline fluo signal}}
$$

625 The analysis of VTA DA calcium variations and its correlation with behavior relied on the construction of
626 peri-event time histograms, which were obtained by aligning and centering specific behavioral events.

627 The identification of these events was done with TTLs received on the fiber photometry console at lever
628 press or nose-poke events. The temporal resolution of the recordings was determined by the bin size
629 chosen for analysis. We used a bin size of 100 ms to achieve a more precise representation of the calcium
630 variations over time. Subsequent analyses were conducted on normalized activity data. To obtain
631 normalized calcium fluorescence variations, a specific formula was applied: Zscore ΔF/F =
632 $\frac{(m\Delta F/F - \mu\Delta F/F)}{\sigma\Delta F/F}$, where Zscore ΔF/F represents the normalized calcium activity for a specific time bin,
mΔF/F represents the average calcium signal variations within that specific time bin, μΔF/F represents the mean calcium signal variations across the baseline activity a few seconds before each events (5 seconds), and σΔF/F represents the standard deviation of the calcium signal variations across the baseline activity a few seconds before each events.

Following normalization, a smoothing technique was employed on the data. This involved convolution with a Gaussian kernel using a sliding window encompassing 100 bins. The "gausswin" function within the MATLAB software was utilized to achieve this smoothing process.

## Optogenetic experiments

In males DAT-Cre mice injected with AAV5-Ef1α-DIO-ChR2(H134R)-EYFP or AAV5-Ef1α-DIO-EYFP, optical stimulations were performed with an ultra-high-power LED (470 nm, Prizmatix) coupled to a patchcord (500 µm core, NA = 0.5, Prizmatix) with an output intensity of 5 – 10 mW. We used a 20 Hz stimulation protocol of 5 ms light pulse for 15 minutes 24 hours and 30 minutes prior to the social task or in vivo electrophysiological recording.

## Modeling

### Building a behavioral model of e-mouse lone and social behaviors

The environment of experiments was modeled as 6 states (Figure 3a, d, rooms 1-4, and lever and dispenser positions). The number and sex of agents ("e-mice") present in the environment was varied, with either 1) one (male or female) e-mouse in the lone experiment, 2) three male or female e-mice in social experiments, or 3) 1 male and 2 females in the mixed-box experiment. State transition occurred at each time step, with transitions probabilities determined by a softmax based on Q values of all accessible k states from the current state $s_i$:

$$p(s_i, s_j) = \frac{e^{\beta Q(s_i, s_j)}}{\sum_{s_k} e^{\beta Q(s_i, s_k)}}$$
with $\beta$, the inverse temperature parameter, favoring either exploration (i.e., random choices; lower $\beta$s) or exploitation of previous learning (higher $\beta$s). $Q$ learning occurred upon each transition from a departure state $s_i$ to an accessible arriving state $s_j$, as

$$Q(s_i, s_j) \leftarrow (1 - \alpha) Q(s_i, s_j) + \alpha (r w d_{s_j} + r w d_{s_j}^{\text{future}}),$$

with $\alpha$ the learning rate, $r w d_{s_j}$ the reward obtained at the arriving state $s_j$ and $r w d_{s_j}^{\text{future}}$ the maximal expected reward obtainable in states accessible from the arriving state $s_j$.

$$r w d_{s_j}^{\text{future}} = \gamma \arg \max_{\forall k} (Q(s_j, s_k)),$$

$\gamma$ being a temporal discount factor. Furthermore, e-mice encountered satiety, which scaled action probabilities, the learning rate, and fatigue, which affected both pressing and eating. Satiety and fatigue were used to scale an action pace in the simulation that was consistent with the experimental measures.

The lever could not be pressed for 3-time steps after each press (i.e., to mimic the 5 sec lever unavailability in experiments). The number of lever presses and the percentage of complete sequences depicted in Figure 3a, e and h were used to categorize e-mice and e-mice triads, as in experiments (Figure 3c, e, g, h, i). Complete modeling information regarding the full model is given in Supplementary Methods.

Reduced model of social interactions.

We built a reduced theoretical model (hereafter termed the reduced model) to assess, within a mathematically tractable framework, the causal mechanisms whereby specialized behaviors emerge under social interactions. To do so, we derived the reduced model from the reinforcement-learning one, based on a continuous time version of $Q$ dynamics (i.e., ordinary differential equations, ODEs). In these ODEs, learning and behavioral dynamics operated at a slower time scale – compared to that of individual choices in the full model – so that actions (i.e., state transitions, pressings, eatings) were described probabilistically. In this framework, we performed qualitative analysis of ODEs to determine the number and stability of fixed points of learning and behavioral (state) variables as a function of parameters (with as focus on $\beta$, which is essential in setting social interactions in the full model, see Results). Moreover, to
reduce dimensionality for better tractability, we considered a simpler setup, where the environment contains only two positions for e-mice (the lever, L, and the food dispenser, D), only two e-mice (which allowed us to assess social interactions), and we did not consider fatigue or satiety. Reduced model ODEs could be expressed under a tractable form:

$$ \left\{ \begin{aligned} \frac{dq_{D1}}{dt} &= \alpha \left( \frac{1+\gamma}{4} \left( 1 + \frac{\beta q_{D1}}{2(1+\gamma)} \right) \left[ 1 - \frac{\beta(q_{D1} + 2q_{D2})}{2(1+\gamma)} \right]^+ r - q_{D1} \right) \\ \frac{dq_{D2}}{dt} &= \alpha \left( \frac{1+\gamma}{4} \left( 1 + \frac{\beta q_{D2}}{2(1+\gamma)} \right) \left[ 1 - \frac{\beta(q_{D2} + 2q_{D1})}{2(1+\gamma)} \right]^+ r - q_{D2} \right) \end{aligned} \right. . $$

with $q_{XYi}$ being the q value of X → Y transition and r the unit reward per eating. This system admitted one or several fixed points ($q_{D1}^*, q_{D2}^*$), which number and nature depended on parameters (see the bifurcation schema and phase portraits, Figure 4a). We focused on the dependence of $q_{Di}^*$ steady-states of individual e-mice on $\beta$, i.e., $q_{Di}^*(\beta)$, denoted branches. As in the full model, the value of $\beta$ separated distinct regimes of social interaction between e-mice. We found that for $\beta$ values below a critical value $\beta_{Bifurc}^+$, only one stable fixed point, $q_{Mono}^{FP}$ existed, with both e-mice admitting a similar q steady-state. Above $\beta_{Bifurc}^+$, $q_{Mono}^{FP}$ became a saddle node. Also, two other symmetric stable fixed points appeared, $q_{Bist}^{FP-A} = (q_{Bist}^{Up*}, q_{Bist}^{Down*})$ and $q_{Bist}^{FP-B} = (q_{Bist}^{Down*}, q_{Bist}^{Up*})$. The fixed point $q_{Bist}^{FP-A}$ corresponded to the case in which e-mouse #1 displayed a higher q value at dispenser ($q_{Bist}^{Up*}$), and, as a consequence, a higher residence probability ($p_{Bist}^{Up*}$) at dispenser, and lower probabilities of lever presses ($p_{Bist}^{LP*}$) and direct sequences ($p_{Bist}^{CS*}$), compared to e-mouse #2 (see below). In other words, $q_{Bist}^{FP-A}$ was a fixed point at which e-mouse #1 behaved as a scrounger and e-mouse #2 as a worker, while behaviors were inverted at $q_{Bist}^{FP-B}$. In this context, the $q_{Mono}^{FP}$ saddle node branch could be interpreted in this reduced model as corresponding to independent e-mice found in both social male boxes in experiments and in the full model (starting from (0,0) with small levels of noise yields animals to this saddle branch, see below). For $\beta$ values above a second critical value $\beta_{NoFood} > \beta_{Bifurc}^+$, the probability of available food at the dispenser became null for the worker e-mouse (see detail of the analysis in Supplementary Methods), the system
704 wrote differently and stable branches were continued by two branches $q_{NoFood}^{Down*}$ and $q_{NoFood}^{Up*}$. The main results of the qualitative analysis of the system are recapitulated in Table 1. Full information on the reduced model derivation and analysis in Supplementary Methods.

707

708 **Statistics**

709 A priori power analyses were not conducted to predetermine the exact number of animals utilized in the experiments. Randomization was implemented at the point of viral infection or behavioral testing to assign animals to their respective groups. Statistical analyses were performed using MATLAB. The normality of data distributions was assessed using the Shapiro-Wilk test. For data sets failing the normality test, non-parametric statistical tests were employed. Normally distributed data were analyzed using independent or paired t-tests depending on the study design. In cases where data normality was violated, the Mann-Whitney test or the Wilcoxon matched-pairs signed-rank test were implemented. For experiments involving multiple comparisons, repeated-measures (RM) ANOVA was used. RM two-way ANOVA was employed for analyses involving two factors, assuming data normality and applying the Bonferroni-Holm correction post-hoc. To compare proportions and ratios between groups, the chi-square ($\chi^2$) test was applied. All statistical tests were two-tailed. Data are presented as mean ± s.e.m., with statistical significance established at a p-value threshold of $P < 0.05$.

721

722 **Data and materials availability:** A table summary of the experimental data is in the Supplementary Material. Source data are provided with this paper. In addition, some of the data generated in this study will be deposited in the Zenodo database. Full behavioral data generated during and/or analyzed during the current study are available from the corresponding author on reasonable request.

726

727 **Code availability:** Code for the model and for the data analyses will be made available in the Zenodo database.

728
729 Acknowledgements: We are grateful for support from Otilia de Oliveira, Noemie Karakaplan-Dherbe and Emilie Tubeuf at ESPCI animal facilities. Clément Solié would like to thank "La Fondation des Treilles" for supporting this work.

732
733 Funding: Centre national de la recherche scientifique (CNRS UMR 8249), ESPCI Paris - PSL, Fondation pour la recherche Médicale (FRM DEQ2013326488 to PF, FRM SPF202005011922 to CS), French National Cancer Institute Grant (SPAV1-21-002 and SPAV1-23-005 to P.F.), French state funds managed by Agence Nationale de la Recherche (ANR-23 VarSeek to PF and BD). The funders had no role in study design, data collection and analysis, decision to publish, or preparation of the manuscript.

738
739 Author Contributions

740 Conceptualization: CS and PF

741 Experimental Methodology: CS AN RJ FdC CV AM SD ND

742 Investigation: CS AN BM CB LR TleB SIF AG JAV

743 Experimental data Analysis: CS BM FM RJ PF

744 Numerical modelling: BD JN BG YL PF

745 Reduced model: BD

746 Funding acquisition: PF BD CS

747 Writing – original draft: CS BD PF

748 Writing – review & editing: CS AN YL CV JPH AM JN

749
750 Competing interest

751 The authors declare no competing financial interests.

752
753
754 Supplementary Figure legends

755

756
<figure><img src="image_56.png" /><figcaption>Supplementary Figure 1</figcaption></figure>
**Supplementary Figure 1:** (a) Number of lever-presses (LPs) through the different days in lone condition for each animal. (b) Time course of the number of beam breaks (BB) for each animal. (c) Time course of the percentage of complete sequences for each animal in lone condition. (d) In-vivo electrophysiology in anesthetized mice. (Left) Spontaneous frequency in hertz (Hz) in function of the percentage of spikes within burst (%SWB) for each neuron recorded, depending on the profiles of the mice. (Right) Average of the %SWB between *Achievers* and *Storers* male mice. (e) Cumulative probability of the number of LPs (left) and of the percentage of complete sequences (right) between *Achievers* and *Storers* mice for the last two days of experiment.

765

766
a

Day1

<figure><img src="image_57.png" /></figure>

Number of LP/hour

<figure><img src="image_58.png" /></figure>

W

<figure><img src="image_59.png" /></figure>

S

<figure><img src="image_60.png" /></figure>

Day1

<figure><img src="image_61.png" /></figure>

Number of Complete Sequence/hour

<figure><img src="image_62.png" /></figure>

W

<figure><img src="image_63.png" /></figure>

S

<figure><img src="image_64.png" /></figure>

b

b=64.5, a=0.5, R²=0.39, p<0.001

<figure><img src="image_65.png" /></figure>

b=46.5, a=0.64, R²=0.47, p<0.001

<figure><img src="image_66.png" /></figure>

b=25.4, a=0.85, R²=0.7, p<0.001

<figure><img src="image_67.png" /></figure>

c

Day1

<figure><img src="image_68.png" /></figure>

Stolen

p<10⁻⁴

<figure><img src="image_69.png" /></figure>

d

Day1

<figure><img src="image_70.png" /></figure>

Thief

p<10⁻⁴

<figure><img src="image_71.png" /></figure>

e

<figure><img src="image_72.png" /></figure>

Male

p=0.014

<figure><img src="image_73.png" /></figure>

Female

p<0.001
**Supplementary Figure 2:** (a) (Left) Time course of the number of LPs during the first day for the different profiles of mice in the social condition. (Middle-Light) Average number of LPs per hours the first day between the different profiles. (Middle-Right / Right) Same analysis than described previously for the percentage of complete sequences. Seq. (b) Correlations analysis for the number of LPs between the last day and day2 (left), day4 (middle) and day6 (right), for each mouse. (c) Number of food pellet stolen to a mouse (Stolen) after an LP during the first day depending on the profile. (d) Number of food pellets stolen by a mouse (Thief) after an LP during the first day depending on the profile. (e) Cumulative probability of the number of LPs (left) and in the percentage of complete sequences (right) between male and female mice in social condition.

776

777

778
<figure><img src="image_74.png" /><figcaption>b</figcaption></figure>

<figure><img src="image_75.png" /><figcaption>c</figcaption></figure>

<figure><img src="image_76.png" /><figcaption>Supplementary Figure 3</figcaption></figure>
**Supplementary Figure 3: (a)** In-vivo electrophysiology in anesthetized mice. (Left) Spontaneous frequency in hertz (Hz) in function of the percentage of spikes within burst (%SWB) for each neuron recorded in female mice in naïve condition and after the social condition. (Right) Average of the VTA DA firing rate between naive and after the task in female mice. (b) Same as (a) for male mice. (c) In-vivo electrophysiology in anesthetized mice. (Left) Spontaneous frequency in function of %SWB for each neuron recorded depending on the behavioral profiles after the social condition. (Right) Average of the VTA DA firing rate between the different profiles after the task.

786

787
a

low $\beta$

<figure><img src="image_77.png" /></figure>

b

<figure><img src="image_78.png" /></figure>

d

<figure><img src="image_79.png" /></figure>

Three e-mice
in a cage

<figure><img src="image_80.png" /></figure>

e

<figure><img src="image_81.png" /></figure>

f

g

<figure><img src="image_82.png" /></figure>

Three e-mice
in a cage

<figure><img src="image_83.png" /></figure>

h

<figure><img src="image_84.png" /></figure>

i

<figure><img src="image_85.png" /></figure>
**Supplementary Figure 4: (a) Q-learning simulation of one e-mouse with « low » β (β=0.5; compare to**
788
Fig. 3a, where β is high). **(b-c) Dependence of LP count (b) and in the percentage of complete sequences**
789
**(c) on the learning rate (α), inverse temperature (β) and temporal discount factor (γ) in the full model and**
790
lone condition. In each panel, the considered parameter is varied, with other parameters set at their
791
standard value (methods). **(d) Same as (a), with three e-mice with low β (β=0.5; compare to Fig. 3d, where**
792
β is high). **(e-f) Dependence of LP count (e) and in the percentage of complete sequences (f) on**
793
parameters, as in (b-c), but in the social condition, with three e-mice having the same parameter values.**
794
**795 (g) Same as (d), with two e-mice with a low β and one e-mouse with a high β (compare to (d), with three**
796
low β e-mice and Fig. 3d, with three high β e-mice). **(h-i) Mean proportions of the three profiles (h) and of**
797
mean probabilities of box types (i) of Fig. 3h.

798

799
a

<figure><img src="image_86.png" /></figure>

b

<figure><img src="image_87.png" /></figure>

<figure><img src="image_88.png" /></figure>

d

<figure><img src="image_89.png" /></figure>

<figure><img src="image_90.png" /></figure>

<figure><img src="image_91.png" /></figure>

e

<figure><img src="image_92.png" /></figure>

<figure><img src="image_93.png" /></figure>

<figure><img src="image_94.png" /></figure>

f

<figure><img src="image_95.png" /></figure>

<figure><img src="image_96.png" /></figure>

<figure><img src="image_97.png" /></figure>

g

<figure><img src="image_98.png" /></figure>

h

<figure><img src="image_99.png" /></figure>

i

<figure><img src="image_100.png" /></figure>

j

<figure><img src="image_101.png" /></figure>
**Supplementary Figure 5: (a-c) Bifurcation diagrams corresponding to Fig. 4a, as a function of the residence probability at the dispenser (a), lever press (b) and complete sequence (c) probabilities in the reduced model. (d-f) Model dynamics corresponding to Fig. 4b-c in the residence probability at Dispenser (left), lever press (center) and complete sequence (right) probabilities planes, for females (d), males (e) and in the mixed box condition. (g-j) Same as (a-c) for the reduced model when no linearization of the softmax function is done. The bifurcation diagrams are qualitatively and quantitatively similar to those of the reduced model with linearization (compare to Fig. 4a, right and (a-c)).**

807

808
<figure><img src="image_102.png" /><figcaption>Supplementary Figure 6</figcaption></figure>
809 **Supplementary Figure 6: (a) Time course of YFP or ChR2 transfection in the VTA DA neurons of DAT-**
810 Cre mice and optic fiber implantation, followed by stimulation protocol 24 hours and 2 hours before in-
811 vivo anesthetized electrophysiology. **(b) In-vivo electrophysiology in anesthetized mice.** Spontaneous
812 frequency in hertz (Hz) in function of the percentage of spikes within burst (%SWB) for each neuron
813 recorded in YFP and ChR2 mice after the stimulation protocol. **(c) Average of the VTA DA firing rate**
814 between YFP and ChR2 mice after the stimulation protocol. **(d) Average of the VTA DA %SWB between**
815 YFP and ChR2 mice after the stimulation protocol. **(e) Comparison of the repartition of the different**
816 behavioral profiles between the group of mice expressing YFP and the group of mice expressing ChR2,
817 in triads containing one ChR2 and two YFP mice.

818

819
## Modeling (Supplementary Methods)

### Building a behavioral model of e-mice lone and social behaviors

The environment of experiments was modeled as 6 states (Fig. 3 and Supp 4; rooms 1-4, and lever and dispenser positions). The number and sex of e-mice present in the environment was varied, with either 1) one (male or female) e-mouse in the lone experiment, 2) three males or females e-mice in social experiments, or 3) 1 male and 2 females in the mixed-box experiment. In the model, we considered e-mice satiety, to account for the convergence of pellets' consumption, $n_{Eat}$, toward a maximum, $n_{Sat}$. Motivation was considered to decrease with satiety:

$$p_{Motiv} = 1 - p_{Sat} = 1 - \exp\left(-\alpha_{Sat}\left[1 - \frac{n_{Eat}}{n_{Sat}}\right]^{+}\right) \quad (1)$$

and scaled both action probabilities and the learning rate (see below). Moreover, e-mice were considered to encounter both pressing and eating fatigue in the model, to get paces consistent with those observed in real mice (without fatigues, paces were unrealistically high in the model). Both fatigue processes were described as action probabilities decreasing upon each individual action and restoring back to their maximal value with first-order dynamics :

$$\left\{ \begin{array}{l} \tau_{PressFat} \frac{dp_{PressFat}}{dt} = p_{PressFat}^{max} - p_{PressFat} \\ p_{PressFat} \leftarrow p_{PressFat} - \Delta p_{PressFat} \delta(t - t_{Press}) \end{array} \right. \quad (2)$$

and

$$\left\{ \begin{array}{l} \tau_{EatFat} \frac{dp_{EatFat}}{dt} = p_{EatFat}^{max} - p_{EatFat} \\ p_{EatFat} \leftarrow p_{EatFat} - \Delta p_{EatFat} \delta(t - t_{Eat}) \end{array} \right. \quad (3),$$

with $\delta$ being the Dirac delta function. Both fatigue processes and satiety concurred to decrease action probabilities in the model. When at the lever, e-mice pressed the lever with probability

$$p_{Press} = p_{PressFat}p_{Motiv} \quad (4).$$

The lever could not be pressed for 3 time steps after each press (i.e., to mimic the 5s lever unavailability in experiments). When at the dispenser, e-mice ate with probability:

$$p_{Eat} = p_{EatFat}p_{Motiv} \quad (5),$$

when pellets were present. State transition occurred at each time step, with softmax probabilities based on Q values of all accessible k states from the current state $s_i$:

$$p(s_i, s_j) = \frac{e^{\beta Q(s_i, s_j)}}{\sum_{\forall k} e^{\beta Q(s_i, s_k)}} \quad (6),$$
with $\beta$ the inverse temperature parameter favoring either exploration (i.e., random choices; lower $\beta$s) or exploitation of previous learning (higher $\beta$s). $Q$ learning occurred upon each transition from a departure state $s_i$ to an accessible arriving state $s_j$, as

$$Q(s_i, s_j) \leftarrow (1 - \alpha) Q(s_i, s_j) + \alpha (r w d_{s_j} + r w d_{s_j}^{\text{future}}) \quad (7),$$

with $\alpha$ the learning rate, $r w d_{s_j}$ the reward obtained at the arriving state $s_j$ and $r w d_{s_j}^{\text{future}}$ the maximal expected reward obtainable in states accessible from the arriving state $s_j$:

$$r w d_{s_j}^{\text{future}} = \gamma \arg \max_{\forall k} \left( Q(s_j, s_k) \right) \quad (8),$$

$\gamma$ being a temporal discount factor. The learning rate scaled with motivation, which was considered to decrease with satiety,

$$\alpha = \alpha_0 (1 - p_{\text{Sat}}) \quad (9).$$

As a consequence, learning stopped when e-mice were satiated. The $Q$ matrix was initialized with value 1 at all elements corresponding to possible transitions (and 0 elsewhere). Consistent with observed trajectories in real mice, stops were forbidden at rooms 3 and 4, i.e., on trajectories back and forth between lever and dispenser states.

### Observables

The number of lever presses was computed – for each e-mice – from the beginning of the simulation up to its satiety. The percentage of complete sequences was computed as the number of direct sequences from lever to dispenser (within 3 times steps, i.e. 6 seconds, as in experiments), normalized by the number of lever presses:

$$P_{\text{CompleteSequence}} = 100 \frac{n_{\text{CompleteSequence}}}{n_{\text{Presses}}} \quad (10).$$

These two observables are depicted in Fig. 3a, e and h and were used to categorize e-mice and e-mice triads, as done in experiments (Fig. 3c, e, g, h, i). In each experimental condition (alone, social, mixed), a hierarchical clustering of e-mice was achieved with Matlab linkage method, based on the number of lever presses and the percentage of complete sequences, after removal of outliers on both variables (5% bilaterally). This procedure was essential as outliers commonly led to unbalanced clustering. This occurred in cases of greatest differences in random realizations of $\beta$ values between e-mice, leading to extreme behaviors (i.e., very « worker » and « scourger » behaviors for very low and high $\beta$s, respectively), which induced irrelevant clusters including one or a few outliers. Also, both variables were normalized before clustering:
$$ \left\{ \begin{aligned} P_{CompleteSequence}^{Normalized} &= \frac{P_{CompleteSequence}}{\langle P_{CompleteSequence} \rangle} \\ n_{Presses}^{Normalized} &= \frac{n_{Presses}}{\langle n_{Presses} \rangle} \end{aligned} \right. \qquad (11), $$

so that the clustering was influenced equally by both variables.

## Reduced model of social interactions.

We built a reduced theoretical model (hereafter termed the reduced model) to assess, within a mathematically tractable framework, the causal mechanisms whereby specialized behaviors emerge under social interactions. To do so, we derived the reduced model from the full one, based on a continuous time version of Q dynamics (i.e., ordinary differential equations, ODEs). In these EDOs, learning and behavioral dynamics operated at a slower time scale – compared to that of individual choices in the full model – so that actions (i.e., state transitions, pressings, eatings) were described probabilistically in the reduced model. In this framework, we performed qualitative analysis of ODEs to determine the number and stability of fixed points of learning and behavioral (state) variables as a function of parameter (with as focus on $\beta$, which is essential in setting social interactions in the full model). Moreover, to reduce dimensionality for better tractability, we considered a simpler setup, where the environment contains only two positions for e-mice (the lever, $L$, and the food dispenser, $D$), only two e-mice (which allowed to assess social interactions) and we did not consider fatigue or satiety. Here, $q$ dynamics were derived from equation (7), to account for their convergence to the expectation of present and future rewards (as for the dynamic of their discrete time equivalent $Q$ values in the full model). Hence,

$$ \left\{ \begin{aligned} \frac{dq_{DDi}}{dt} &= \alpha (rwd_{Di} + rwd_{Di}^{Future} - q_{DDi}) \\ \frac{dq_{LDi}}{dt} &= \alpha (rwd_{Di} + rwd_{Di}^{Future} - q_{LDi}) \\ \frac{dq_{DLi}}{dt} &= \alpha (rwd_{Li} + rwd_{Li}^{Future} - q_{DLi}) \\ \frac{dq_{LLi}}{dt} &= \alpha (rwd_{Li} + rwd_{Li}^{Future} - q_{LLi}) \end{aligned} \right. \qquad (12), $$

with $q_{XYi}$ being the $q$ value of $X \to Y$ transition, and $rwd_{Xi}$ and $rwd_{Xi}^{Future}$ the present and future expected reward at position $X$ in e-mouse $i$. Note that, despite all the simplifications made, the dynamical system under consideration is still 8-dimensional. Fortunately, it is straightforward to note that dynamics of $q_{DDi}$ and $q_{LDi}$ are identical (as that of $q_{LDi}$ and $q_{LLi}$), so that, starting from identical null initial conditions
$$q_{LLi}(t=0) = q_{DLi}(t=0) = q_{LDi}(t=0) = q_{DDi}(t=0) = 0 \quad (13),$$

one has

$$\begin{cases} q_{LLi}(t) = q_{DLi}(t) \\ q_{DDi}(t) = q_{LDi}(t) \end{cases} \quad (14),$$

so that only one q value at each position is sufficient to describe dynamics:

$$\begin{cases} q_{Li} \equiv q_{LLi} = q_{DLi} \\ q_{Di} \equiv q_{DDi} = q_{LDi} \end{cases} \quad (15),$$

reducing the theory to a 4-dimensional problem. For each e-mouse, one could express q values evolved according to :

$$\left\{ \begin{aligned} \frac{dq_{Li}}{dt} &= \alpha (r w d_{Li} + r w d_{Li}^{\text{Future}} - q_{Li}) \\ \frac{dq_{Di}}{dt} &= \alpha (r w d_{Di} + r w d_{Di}^{\text{Future}} - q_{Di}) \end{aligned} \right. \quad (16).$$

The mean expectation of current reward at L was null ($r w d_{Li} = 0$), while that at D scaled with the probability of being at the dispenser, the probability of eating, the probability of available food at the dispenser, and the reward value of eating, which writes, for e-mouse #1, e.g.,

$$r w d_{D1} = p_{Eat} p_{D1} \left[ \underbrace{\frac{p_{Press}(p_{L1} + p_{L2})}{p_{common \ production}}}_{p_{food at Dispenser}} - \underbrace{\frac{p_{Eat} p_{D2}}{p_{consumption e-mouse 2}}} \right]^{+} r \quad (17),$$

with $p_{Xi}$ the residence probability of being at position X for mouse i, $p_{Eat}$ and $p_{Press}$ the probabilities of eating and pressing, and r the unit reward per eating. The positive part function, $[]^+$, is necessary to account for a positive probability of available food at the dispenser. Note also that because $p_{Di} > 0.5$ and $p_{Li} = 1 - p_{Di} < 0.5$ (as there is a reward at D and not at L), the probability of available food at the dispenser is inferior to 1. Note that equation (17) can be written for mouse #2's by swapping indices #1 and #2, as can be done for other equations in the following. Whenever possible, we keep the generic index i for e-mice. By else, at both positions, the expectation of future rewards expressed as

$$r w d_{Di}^{\text{Future}} = \gamma r w d_{Di} \quad (18),$$

as getting a reward at D is the maximal expected possible reward at the next step from both positions L and D. Therefore, q derivatives can be re-expressed as
$$ 
\left\{ 
\begin{aligned}
\frac{dq_{L1}}{dt} &= \alpha \left( \gamma p_{Eat} p_{D1} \left[ p_{Press} (p_{L1} + p_{L2}) - p_{Eat} p_{D2} \right]^+ r - q_{L1} \right) \\
\frac{dq_{D1}}{dt} &= \alpha \left( (1 + \gamma) p_{Eat} p_{D1} \left[ p_{Press} (p_{L1} + p_{L2}) - p_{Eat} p_{D2} \right]^+ r - q_{D1} \right)
\end{aligned}
\right.
\quad (19). $$

In the following, we considered, for the sake of simplicity that $p_{Press} = p_{Eat} = 1$. Also, because $p_{Li} = 1 - p_{Di}$ for both e-mice i, q derivatives simply expressed as

$$ 
\left\{ 
\begin{aligned}
\frac{dq_{L1}}{dt} &= \alpha \left( \gamma p_{D1} \left[ 2 - p_{D1} - 2p_{D2} \right]^+ r - q_{L1} \right) \\
\frac{dq_{D1}}{dt} &= \alpha \left( (1 + \gamma) p_{D1} \left[ 2 - p_{D1} - 2p_{D2} \right]^+ r - q_{D1} \right)
\end{aligned}
\right.
\quad (20). $$

Note that the first equation is actually the second one scaled by a factor $\frac{\gamma}{1 + \gamma}$. As both $q_{D1}$ and $q_{L1}$ both start from null initial conditions (equation (13)), their dynamics are homothetic, so that at any point in time

$$ q_{Li} = \frac{\gamma}{1 + \gamma} q_{Di} \quad (21). $$

In the following, we proceed with $q_{D1}$ equation, which, together with that for $q_{D2}$, describes the system with two-dimensional dynamics:

$$ 
\left\{ 
\begin{aligned}
\frac{dq_{D1}}{dt} &= \alpha \left( (1 + \gamma) p_{D1} \left[ 2 - p_{D1} - 2p_{D2} \right]^+ r - q_{D1} \right) \\
\frac{dq_{D2}}{dt} &= \alpha \left( (1 + \gamma) p_{D2} \left[ 2 - p_{D2} - 2p_{D1} \right]^+ r - q_{D2} \right)
\end{aligned}
\right.
\quad (22). $$

In such a setup, with two states, residence probabilities of e-mouse i, i.e., $p_{Li}$ and $p_{Di}$ can be directly related to transition probabilities. Indeed, assuming first-order state transitions, $p_{Di}$ evolves as:

$$ \frac{dp_{Di}}{dt} = p_{Li}p_{LDi} - p_{Di}p_{DLi} \quad (23), $$

which, after substituting $p_{Li}$ by $1 - p_{Di}$ leads to

$$ \frac{dp_{Di}}{dt} = p_{LDi} - (p_{LDi} + p_{DLi}) p_{Di} \quad (24), $$

with a unique fixed point

$$ p_{Di} = \frac{p_{LDi}}{p_{LDi} + p_{DLi}} \quad (25), $$

so that replacing $p_{LDi}$ and $p_{DLi}$ by their softmax expression, leads to

$$ p_{Di} = p_{LDi} = p_{DDi} = \frac{e^{\beta q_{LDi}}}{e^{\beta q_{LDi}} + e^{\beta q_{DLi}}} \quad (26) $$
and

$$p_{Li} = p_{DLi} = p_{LLi} = \frac{e^{\beta q_{DLi}}}{e^{\beta q_{LDi}} + e^{\beta q_{DLi}}} \quad (27),$$

which, thanks to equation (15), can simply be written

$$p_{Di} = \frac{e^{\beta q_{Di}}}{e^{\beta q_{Li}} + e^{\beta q_{Di}}} \quad (28)$$

and

$$p_{Li} = \frac{e^{\beta q_{Li}}}{e^{\beta q_{Li}} + e^{\beta q_{Di}}} \quad (29).$$

Combining equations (21) and (28) allowed to express the residence probability as a function of the q value at dispenser:

$$p_{Di} = \left( 1 + e^{-\frac{\beta}{1+\gamma} q_{Di}} \right)^{-1} \quad (30).$$

Then, combining equations (22) and (30) yields

$$\left\{ \begin{aligned} \frac{dq_{D1}}{dt} &= \alpha \left( \frac{1+\gamma}{1+e^{-\frac{\beta}{1+\gamma} q_{D1}}} \left[ 2 - \frac{1}{1+e^{-\frac{\beta}{1+\gamma} q_{D1}}} - \frac{2}{1+e^{-\frac{\beta}{1+\gamma} q_{D2}}} \right]^{+} r - q_{D1} \right) \\ \frac{dq_{D2}}{dt} &= \alpha \left( \frac{1+\gamma}{1+e^{-\frac{\beta}{1+\gamma} q_{D2}}} \left[ 2 - \frac{1}{1+e^{-\frac{\beta}{1+\gamma} q_{D2}}} - \frac{2}{1+e^{-\frac{\beta}{1+\gamma} q_{D1}}} \right]^{+} r - q_{D2} \right) \end{aligned} \right. \quad (31).$$

Note that this system describes the social interaction between two e-mice having similar β values (i.e., two males or femelles). This two-dimensional EDO system describes the time evolution of learning and behavior variables (as residence, lever press and direct sequence probabilities can be derived from q values, see below). Unfortunately, this system (even studied by parts, depending on whether probabilities of food at dispenser ate positive or null) has no analytical solution. We therefore assessed a version of the EDOs with linearized approximations of residence probabilities. This linearized version could be solved analytically, which allowed to assess the number and stability of fixed points as a function of β. To do so, we used a linearized approximation of residence probabilities at position X in e-mouse i, i.e.

$$p_{Xi} = \frac{1}{2} + \frac{\beta}{4(1+\gamma)} q_{Xi} \quad (32).$$
Such approximation was typically of the order of one to a few percents. Combined with equation (22), it allowed to provide the tractable system we analyzed to generate Figs. 4 et Supp 5 and which, after a few lines of algebraic manipulation, expresses as:

$$ \left\{ \begin{aligned} \frac{dq_{D1}}{dt} &= \alpha \left( \frac{1+\gamma}{4} \left( 1 + \frac{\beta q_{D1}}{2(1+\gamma)} \right) \left[ 1 - \frac{\beta(q_{D1} + 2q_{D2})}{2(1+\gamma)} \right]^+ r - q_{D1} \right) \\ \frac{dq_{D2}}{dt} &= \alpha \left( \frac{1+\gamma}{4} \left( 1 + \frac{\beta q_{D2}}{2(1+\gamma)} \right) \left[ 1 - \frac{\beta(q_{D2} + 2q_{D1})}{2(1+\gamma)} \right]^+ r - q_{D2} \right) \end{aligned} \right. \quad (33). $$

In the reduced model, the propensity of e-mouse i to press the lever was expressed as a probability, which scaled both with the probability to be at the lever position and the probability to press, once at that position:

$$ P_{LPi} = P_{Press} P_{Li} = P_{Press} (1 - P_{Di}) \quad (34). $$

The propensity to achieve direct sequences was also expressed as a probability, scaling both with the be at the lever position, to press the lever and then move toward the dispenser position (which is equal to the residence probability at dispenser (equation (26)):

$$ P_{Si} = P_{Press} P_{Li} P_{Di} = P_{Press} (1 - P_{Di}) P_{Di} \quad (35). $$

## Simulation of the reduced model

The reduced model was simulated using a noisy version of system (19):

$$ \left\{ \begin{aligned} \frac{dq_{L1}}{dt} &= \alpha \left( \gamma p_{Eat} p_{D1} \left[ p_{Press} (p_{L1} + p_{L2}) - p_{Eat} p_{D2} \right]^+ r - q_{L1} \right) + \sigma_{Noise} \sqrt{dt} \epsilon(t) \\ \frac{dq_{D1}}{dt} &= \alpha \left( (1+\gamma) p_{Eat} p_{D1} \left[ p_{Press} (p_{L1} + p_{L2}) - p_{Eat} p_{D2} \right]^+ r - q_{D1} \right) + \sigma_{Noise} \sqrt{dt} \epsilon(t) \end{aligned} \right. \quad (36), $$

where $\sigma_{Noise}$ is the noise standard deviation and $\epsilon(t)$ is a Gaussian stochastic process with zero mean and unit standard deviation, similar equations for e-mouse #2 and residence probabilities computed according to equation (30). Simulations were used to check the validity of the qualitative analyse achieved to produce bifurcation diagrams and illustrate system's dynamics in phase portraits (see below).

## Qualitative analysis of the reduced model

The system (33) admitted one or several fixed points $(q_{D1}^*, q_{D2}^*)$, which number and nature depended on parameters. We focused on the dependence of $q_{Di}^*$ steady-states of individual e-mice on $\beta$, i.e. $q_{Di}^*(\beta)$, denoted branches in the following (see figures XXX). Indeed, as in the full model, the value of $\beta$ separated distinct regimes of social interaction between e-mice. The main results of the qualitative analysis of the system are recapitulated in Table 1.
We found that for $\beta$ values below a critical value $\beta_{Bifurc}$ (see below, for the analytical calculation of critical values), only one fixed point, $q_{Mono}^{FP}$ existed, with both e-mice admitting a similar $q$ steady-state:

$$q_{Mono}^{FP} = (q_{Mono}^{\star}, q_{Mono}^{\star}) \quad (37)$$

with

$$q_{Mono}^{\star} = \frac{2(1+\gamma)}{3\beta^2 r} \left( 2\sqrt{\beta^2 r^2 + 2\beta r + 4} - (\beta r + 4) \right) \quad (38).$$

Its stability was set by

$$\lambda_{Mono}^{1,2} = -\alpha \left[ 1 + \frac{\beta r}{8(1+\gamma)} \left( 2\beta q_{Mono}^{\star} \pm (2(1+\gamma) + \beta q_{Mono}^{\star}) \right) \right] \quad (39),$$

with both Eigen values negative, so that $q_{Mono}^{FP}$ was a stable node.

Above $\beta_{Bifurc}^{+}$, $\lambda_{Mono}^{2}$ was positive, so that $q_{Mono}^{FP}$ became a saddle node. Also, two other symmetric fixed points appeared, $q_{Bist}^{FP-A} = (q_{Bist}^{Up\star}, q_{Bist}^{Down\star})$ and $q_{Bist}^{FP-B} = (q_{Bist}^{Down\star}, q_{Bist}^{Up\star})$, with

$$q_{Bist}^{Up,Down\star} = \frac{2(1+\gamma)}{\beta^2 r} \left( \beta r - 4 \pm 2\sqrt{\beta^2 r^2 - 6\beta r + 4} \right) \quad (40).$$

The fixed point $q_{Bist}^{FP-A}$ corresponded to the case in which e-mouse #1 displayed a higher $q$ value at dispenser ($q_{Bist}^{Up\star}$), and, as a consequence, a higher residence probability ($p_{Bist}^{Up\star}$) at dispenser, and lower probabilities of lever presses ($p_{Bist}^{LP\star}$) and direct sequences ($p_{Bist}^{CS\star}$), compared to e-mouse #2 (see below). In other words, $q_{Bist}^{FP-A}$ was a fixed point at which e-mouse #1 behaved as a scrounger and e-mouse #2 as a worker, while behaviors were inverted at $q_{Bist}^{FP-B}$. In this context, the $q_{Mono}^{FP}$ saddle node branch could be interpreted in this reduced model as corresponding to independent e-mice found in social male boxes in experiments and the full model (starting from (0,0) with small levels of noise yields animals to this saddle branche, see below). The stability of both fixed points was set by

$$\lambda_{Bist}^{1,2} = \alpha \left[ 1 - \frac{\beta r}{2} \pm \frac{\beta r}{8(1+\gamma)} \sqrt{(2(1+\gamma) + \beta q_{Bist}^{Up\star})(2(1+\gamma) + \beta q_{Bist}^{Down\star})} \right] \quad (41),$$

with both Eigen values negative, so that both $q_{Bist}^{FP-A}$ and $q_{Bist}^{FP-B}$ were stable nodes.

For $\beta$ values above a second critical value $\beta_{NoFood}$ (with $\beta_{NoFood} > \beta_{Bifurc}$, see below), the probability of available food at the dispenser became null for the worker e-mouse (i.e., the argument within the positive part function, $[\cdot]^{+}$, becomes negative for e-mouse #1 or 2 in equation (33)). The system thus writes, if the worker is e-mouse #2,
$$ egin{cases} \frac{dq_{D1}}{dt} = \alpha \left( \frac{1 + \gamma}{4} \left( 1 + \frac{\beta q_{D1}}{2(1 + \gamma)} \right) \left[ 1 - \frac{\beta(q_{D1} + 2q_{D2})}{2(1 + \gamma)} \right]^+ r - q_{D1} \right) \\ \\ \frac{dq_{D2}}{dt} = -\alpha q_{D2} \end{cases} \quad (42), $$

A single fixed point $q_{NoFood}^{FP-A} = (q_{NoFood}^{Up\star}, q_{NoFood}^{Down\star})$ existed

$$ egin{cases} q_{NoFood}^{Up\star} = \frac{2(1 + \gamma)}{\beta^2 r} \left( \sqrt{\beta^2 r^2 + 16} - 4 \right) \\ q_{NoFood}^{Down\star} = 0 \end{cases} \quad (43), $$

which stability was determined by

$$ egin{cases} \lambda_{NoFood}^{Up} = -\frac{\alpha}{4} \sqrt{\beta^2 r^2 + 16} \\ \lambda_{NoFood}^{Down} = -\alpha \end{cases} \quad (44), $$

with both Eigen values negative, so that $q_{NoFood}^{FP-A}$ was a stable node. Obviously, a symmetrical stable node, $q_{NoFood}^{FP-B} = (q_{NoFood}^{Down\star}, q_{NoFood}^{Up\star})$, co-existed with $q_{NoFood}^{FP-A}$ when the worker was e-mouse #1.

Both $\beta_{Bifurc}$ and $\beta_{NoFood}$ critical values could be found analytically. At the bifurcation point where $q_{Mono}^{FP}$ becomes a saddle node, the two symmetric fixed points $q_{Bist}^{FP-A} = (q_{Bist}^{Up\star}, q_{Bist}^{Down\star})$ and $q_{Bist}^{FP-B} = (q_{Bist}^{Down\star}, q_{Bist}^{Up\star})$ switch from complex ($\beta < \beta_{Bifurc}$) to real ($\beta \ge \beta_{Bifurc}$) values. Thus, this point was set by $\beta$ by a simple criterion, i.e., canceling of the term under the square root function in equation (40):

$$ \beta_{Bifurc}^2 r^2 - 6\beta_{Bifurc} r + 4 = 0 \quad (45), $$

which yields

$$ \beta_{Bifurc}^{+/-} = r^{-1} \frac{6 \pm \sqrt{20}}{2} \quad (46), $$

with $\beta_{Bifurc} = r^{-1} \frac{6 - \sqrt{20}}{2}$ separating branches of real negative $q_{Bist}^{Up,Down\star}$ values from the branch of complex values between $\beta_{Bifurc}^-$ and $\beta_{Bifurc}^+$. The branch of real negative $q_{Bist}^{Up,Down\star}$ values below $\beta_{Bifurc}^-$ corresponded to impossible solutions, as $q$ values are positive (with a null initial condition (equation (13)) and positive reward expectations). Finally,

$$ \beta_{Bifurc} \equiv \beta_{Bifurc}^{+} = r^{-1} \frac{6 + \sqrt{20}}{2} \approx 5.2361 \quad (47), $$
was the critical point (assuming $r = 1$; see parameter values, below) at which the system switched from one to three fixed points (two stable ones separated by an unstable one), at a supercritical pitchfork bifurcation (« trifurcation »; Strogatz, 1994). At this point, where $q_{Bist}^{Up\star} = q_{Bist}^{Down\star}$, i.e., at which $q_{Bist}^{Up,Down\star} = \frac{2(1+\gamma)}{\beta^2 r}(\beta r - 4)$, one could easily check, using equation (39), that $\lambda_{Mono}^2 = 0$, i.e. $q_{Mono}^{FP}$ loses its stability along the second Eigen direction (i.e. switching from a stable to a saddle node).

The other critical value, $\beta_{NoFood}$, could be found by considering that $q_{Bist}^{Down\star} = 0$ at that point, i.e.,

$$ \beta_{NoFood}r - 4 \pm 2\sqrt{\beta_{NoFood}^2 r^2 - 6\beta_{NoFood}r + 4} = 0 \quad (48) $$

which yields

$$ \beta_{NoFood} = r^{-1}\frac{16}{3} \approx 5.3333 \quad (49). $$

Bifurcation diagrams were built (Figs. 4 et Supp 5) by representing $q_{Mono}^*$ stable and saddle branches (below and above $\beta_{Bifurc}$, respectively, $q_{Bist}^{Up\star}$ and $q_{Bist}^{Down\star}$ branches for $\beta_{Bifurc} \le \beta < \beta_{NoFood}$, and $q_{NoFood}^{Up\star}$ and $q_{NoFood}^{Down\star}$ branches for $\beta \ge \beta_{NoFood}$.

We verified that bifurcation diagrams obtained using linearized approximation of residence probabilities (equation (32)) were indeed qualitatively correct and quantitatively consistent with the original reduced model, using numerical simulations of the reduced model. To do so, the stable and bistable branches were computed by taking steady-state $q_{Bi}$ values of these simulations, and computing $p_{Bi}$, $p_{Li}$, and $p_{Si}$ from equations (30), (34) and (35), respectively. The $q_{Mono}^{FP}$ saddle node branch was acquired by doing the same, in the absence of noise ($\sigma_{Noise} = 0$), so that system's state converged to the saddle node from its null initial condition (i.e., along the first bisector in the $(q_{B1}, q_{B2})$ plane, due to system's symmetry) and stayed at it, without further diverging toward $q_{Bist}^{FP-A}$ or $q_{Bist}^{FP-B}$ (which occurred in presence of noise), as Eigen directions are stable along the first bisector and unstable along the second bisector (see below).

## Eigen Vectors of the reduced model

Eigen Vectors in the general case expressed as

$$ V_{1,2}^{Mono,Bist} = \left( \pm \frac{\sqrt{(2(1+\gamma) + \beta q_{B1}^*)(2(1+\gamma) + \beta q_{B2}^*)}}{2(1+\gamma) + \beta q_{B2}^*} \right) \quad (50), $$

which, in the case of the $q_{Mono}^{FP} = (q_{Mono}^*, q_{Mono}^*)$ branch, trivially simplified to

$$ V_{1,2}^{Mono} = \left( \pm \frac{1}{1} \right) \quad (51). $$

In the case with no food for one e-mouse, Eigen Vectors expressed as
$$V_1^{\text{NoFood}} = \left( -\frac{2(1 + \gamma) + \beta q_{B1}^*}{\beta(q_{B1}^* + q_{B2}^*)} \right) \quad V_2^{\text{NoFood}} = \begin{pmatrix} 1 \\ 0 \end{pmatrix} \quad (52)$$

which, in the case of the $q_{NoFood}^{FP-A}$ simplified to

$$V_1^{\text{NoFood}} = \left( -\frac{\beta r + 4 + \sqrt{\beta^2 r^2 + 16}}{\beta r} \right) \quad V_2^{\text{NoFood}} = \begin{pmatrix} 1 \\ 0 \end{pmatrix} \quad (53).$$

## Phase portraits of the reduced model

In phase portraits (Figs. 4 et Supp 5), the $q$ vector field was numerically computed from equation (19), with $q_{Di}$ values set within a fixed range and $q_{Li}$ computed from equation (21).

The $p$ vector field was computed with $p_{Di}$ values set within a fixed range, $q_{Di}$ values computed as

$$q_{Di} = \frac{1 + \gamma}{\beta} \log \left( \frac{p_{Di}}{1 - p_{Di}} \right) \quad (54),$$

i.e., the reciprocal equation of equation (30), $q_{Li}$ computed from $q_{Di}$ with equation (21) and $p$ derivatives computed thanks to the chain rule, i.e.,

$$\frac{dp_{Di}}{dt} = \frac{dp_{Di}}{dq_{Di}} \frac{dq_{Di}}{dt} \quad (55),$$

where

$$\frac{dp_{Di}}{dq_{Di}} = \frac{\beta e^{\frac{\gamma q_{D1}}{1+\gamma}} e^{\beta q_{D1}}}{(1+\gamma) \left( e^{\frac{\gamma\beta}{1+\gamma} q_{D1}} + e^{\beta q_{D1}} \right)^2} \quad (56)$$

was obtained by differentiating equation (30) and $\frac{dq_{Di}}{dt}$ taken from system (19).

The $p_{LP}$ vector field was computed with $p_{Li}$ values set within a fixed range, $p_{Di}$ values computed as

$$p_{Di} = 1 - \frac{p_{LPi}}{p_{Press}} \quad (57),$$

which was derived from equation (34), $q_{Di}$ values with equation (21), $\frac{dq_{Di}}{dt}$ from system (19), $\frac{dp_{Di}}{dq_{Di}}$ from equation (56) and by applying the chain rule, i.e.,

$$\frac{dp_{LPi}}{dt} = \frac{dp_{LPi}}{dp_{Di}} \frac{dp_{Di}}{dq_{Di}} \frac{dq_{Di}}{dt} \quad (58),$$
with

$$ \frac{dp_{LPi}}{dp_{Di}} = -p_{Press} (59), $$

which was obtained by differentiating equation (34).

The $p_S$ vector field was computed as the $p_{LP}$ vector field, using equations (34), (21), (19) and (56) and applying the chain rule

$$ \frac{dp_{Si}}{dt} = \frac{dp_{Si}}{dp_{Di}} \frac{dp_{Di}}{dq_{Di}} \frac{dq_{Di}}{dt} (60), $$

with

$$ \frac{dp_{Si}}{dp_{Di}} = -p_{Press} (1 - 2p_{Di}) (61), $$

which was obtained by differentiating equation (35).

In phase portraits, vector fields were represented both by a direction field plot and through a color coded map. Nullclines were estimated as zero level contours of the vector fields. Fixed points were taken with bistable Up and Down steady-state values obtained from the bifurcation diagram for $\beta = \beta_{Female}$ (social female condition), $\beta = \beta_{Male}$ (social male condition). In the mixed box condition, we considered one female and one male and the fixed point was taken as the intersection of both nullclines. Eigen vector were not represented for clarity. An example of simulation of reduced model's dynamics was also represented in each condition to show dynamics from initial conditions toward the – or one of the – fixed points.

## Numerical procedures

Numerical simulations and qualitative analysis were achieved thanks to Matlab custom scripts that are available online. The qualitative analysis of the reduced model was done using the Symbolic toolbox of Matlab. Numerical simulations of the full and reduced model were performed using Euler forward integration. In the full model, standard parameters, unless mentioned, were $\alpha_0 = 0.05$, $\gamma = 0.99$, $\alpha_{Sat} = 20$, $n_{Sat} = 150$, $\tau_{PressFat} = 50s$, $p_{PressFat}^{max} = 1$, $\Delta p_{PressFat} = 0.5$, $\tau_{EatFat} = 50s$, $p_{EatFat}^{max} = 1$, $\Delta p_{EatFat} = 0.5$, and the time step was $\Delta t = 2s$. Simulations lasted $t_{max} = 20000s$. $\beta$ distributions are discussed in the main text. Parameters of the reduced model were as following: $\alpha = 0.01$, $\beta = [1,10]$, $\beta_{Female} = 1$, $\beta_{Male} = 10$, $\gamma = 0.99$, $r = 1$, $p_{Press} = 1$, $p_{Eat} = 1$, $\sigma_{Noise} = 1e^{-4}$, $dt = 1s$. Simulations lasted $t_{max} = 1500s$ for example simulations shown in phase portraits and $t_{max} = 100000s$ in simulations with noise, i.e., a much longer time to estimate steady-states in the numerical estimation of phase portraits. The distribution of e-mice
behavioral profiles and box type distributions were obtained on 3000 repetitions of 16 boxes of either the lone condition, or single sex or mixed sex conditions, with 16 animals per box (Fig. 3c, g, and Supp. 4 h, i).

## References

*Nonlinear dynamics and chaos, Strogatz, CRC Press, Chapman & Hall Book. Second edition. (1994) Steven H. Strogatz, ISBN 0-201-54344-3*
823 Extended Data

824 Table 1: Model

825

<table><caption>826 Table 1</caption><thead><tr><th rowspan="2"></th><th rowspan="2">model</th><th colspan="3">general case<br>(food available for both e-mice)</th><th rowspan="2">no food available<br>for e-mouse #2 (e.g.)</th></tr><tr><th>β</th><th>β &lt; β<sub>Bifurc</sub></th><th>β<sub>Bifurc</sub> ≤ β &lt; β<sub>NoFood</sub></th><th>β ≥ β<sub>NoFood</sub></th></tr></thead><tbody><tr><td rowspan="4">q<sup>PF</sup><sub>Mono</sub></td><td>q<sup>*</sup><sub>B1</sub> = q<sup>*</sup><sub>Mono</sub></td><td rowspan="2">&gt; 0</td><td rowspan="2"></td><td rowspan="2">&gt; 0</td><td rowspan="2"></td></tr><tr><td>q<sup>*</sup><sub>B2</sub> = q<sup>*</sup><sub>Mono</sub></td></tr><tr><td>λ<sup>1</sup><sub>Mono</sub></td><td rowspan="2">&lt; 0</td><td rowspan="2"></td><td>&lt; 0</td><td rowspan="2"></td></tr><tr><td>λ<sup>2</sup><sub>Mono</sub></td><td>&gt; 0</td></tr><tr><td></td><td>type</td><td>stable node</td><td></td><td>saddle node</td><td></td></tr><tr><td rowspan="4">q<sup>PF-A</sup><sub>Bist</sub></td><td>q<sup>*</sup><sub>B1</sub> = q<sup>Up*</sup><sub>Bist</sub></td><td rowspan="4">no real positive<br>solution</td><td rowspan="2">&gt; 0</td><td rowspan="4"></td><td rowspan="4">NA</td></tr><tr><td>q<sup>*</sup><sub>B2</sub> = q<sup>Down*</sup><sub>Bist</sub></td></tr><tr><td>λ<sup>1</sup><sub>Bist</sub></td><td rowspan="2">&lt; 0</td></tr><tr><td>λ<sup>2</sup><sub>Bist</sub></td></tr><tr><td></td><td>type</td><td></td><td>stable node</td><td></td><td></td></tr><tr><td rowspan="4">q<sup>PF-A</sup><sub>BistNoF</sub></td><td>q<sup>*</sup><sub>B1</sub> = q<sup>Up*</sup><sub>BistNoF</sub></td><td rowspan="4">NA</td><td rowspan="4"></td><td>&gt; 0</td><td rowspan="4"></td></tr><tr><td>q<sup>*</sup><sub>B2</sub> = q<sup>Down*</sup><sub>BistNoF</sub></td><td>0</td></tr><tr><td>λ<sup>1</sup><sub>BistNoF</sub></td><td>&lt; 0</td></tr><tr><td>λ<sup>2</sup><sub>BistNoF</sub></td><td></td></tr><tr><td></td><td>type</td><td></td><td></td><td>stable node</td><td></td></tr></tbody></table>

827

828 The table recapitulates the name of the fixed-points (first column), steady-states for both e-mice and their eigenvalues (second column), as well as their value or type, depending of the ranges of β. NA: not applicable.

831

832 Table 2: Figures statistical tests

833
<table><tbody><tr><td rowspan="5">Figure 1</td><td>Figure 1b</td><td>One-way RM ANOVA</td><td>Time main effect LP : F<sub>(3,61)</sub> = 35.10; p &lt; 0.0001<br>Time main effect BB : F<sub>(3,61)</sub> = 86.59; p &lt; 0.0001</td></tr><tr><td>Figure 1c (right)</td><td>One-way RM ANOVA</td><td>Time main effect: F<sub>(3,61)</sub> = 30.37; p &lt; 0.0001</td></tr><tr><td>Figure 1f</td><td>Two-way RM ANOVA</td><td>Time main effect: F<sub>(3,60)</sub> = 49.97; p &lt; 0.0001<br>Group main effect : F<sub>(1,60)</sub> = 83.01; p &lt; 0.0001<br>Interaction Time x Group : F<sub>(3,60)</sub> = 40.35; p &lt; 0.0001</td></tr><tr><td>Figure 1g</td><td>χ2-Test</td><td>χ2 = 6.62; p = 0.0101</td></tr><tr><td>Figure 1h</td><td>Wilcoxon rank test</td><td>p = 0.0009</td></tr><tr><td rowspan="7">Figure 2</td><td>Figure 2c (up left)</td><td>Two-way RM ANOVA<br>followed by Bonferroni-Holm<br>correction</td><td>Time main effect: F<sub>(6,80)</sub> = 19.40; p &lt; 0.0001<br>Group main effect : F<sub>(3,80)</sub> = 44.41; p &lt; 0.0001<br>Interaction Time x Group : F<sub>(18,80)</sub> = 7.58; p &lt; 0.0001</td></tr><tr><td>Figure 2c (up right)</td><td>One-way ANOVA followed<br>by Bonferroni-Holm<br>correction</td><td>Group main effect: F<sub>(3,80)</sub> = 42.23; p &lt; 0.0001</td></tr><tr><td>Figure 2c (down left)</td><td>Two-way RM ANOVA<br>followed by Bonferroni-Holm<br>correction</td><td>Time main effect: F<sub>(6,80)</sub> = 44.80; p &lt; 0.0001<br>Group main effect : F<sub>(3,80)</sub> = 105.18; p &lt; 0.0001<br>Interaction Time x Group : F<sub>(18,80)</sub> = 15.61; p &lt; 0.0001</td></tr><tr><td>Figure 2c (down right)</td><td>One-way ANOVA followed<br>by Bonferroni-Holm<br>correction</td><td>Group main effect: F<sub>(3,80)</sub> = 101.29; p &lt; 0.0001</td></tr><tr><td>Figure 2d (left)</td><td>Unpaired t-test</td><td>t<sub>(82)</sub> = -1.50; p = 0.14</td></tr><tr><td>Figure 2d (right)</td><td>Unpaired t-test</td><td>t<sub>(82)</sub> = 12.08; p &lt; 0.0001</td></tr><tr><td>Figure 2e</td><td>χ2-Test</td><td>χ2 = 71.37; p &lt; 0.0001</td></tr><tr><td>Figure 5a (left)</td><td>Unpaired t-test</td><td>t<sub>(34)</sub> = -5.61; p &lt; 0.0001</td></tr><tr><td>Figure 5a (right)</td><td>Unpaired t-test</td><td>t<sub>(34)</sub> = -1.49; p = 0.15</td></tr><tr><td>Figure 5b (up)</td><td>One-way ANOVA followed<br>by Bonferroni-Holm<br>correction</td><td>Group main effect: F<sub>(2,33)</sub> = 43.6; p &lt; 0.0001</td></tr></tbody></table>
<table><tr><td rowspan="7">Figure 5</td><td>Figure 5b (down)</td><td>One-way ANOVA followed by Bonferroni-Holm correction</td><td>Group main effect: F(2,33) = 14; p &lt; 0.0001</td></tr><tr><td>Figure 5c</td><td>χ2-Test</td><td>χ2 = 19.42; p &lt; 0.0001</td></tr><tr><td>Figure 5i (up)</td><td>Two-way RM ANOVA</td><td>Time main effect: F(6,110) = 29.98; p &lt; 0.0001<br>Group main effect : F(3,110) = 58.80; p &lt; 0.0001<br>Interaction Time x Group : F(18,110) = 8.75; p &lt; 0.0001</td></tr><tr><td>Figure 5i (down)</td><td>Two-way RM ANOVA</td><td>Time main effect: F(6,110) = 80.57; p &lt; 0.0001<br>Group main effect : F(3,110) = 43.41; p &lt; 0.0001<br>Interaction Time x Group : F(18,110) = 3.17; p &lt; 0.0001</td></tr><tr><td>Figure 5j (up)</td><td>One-way ANOVA followed by Bonferroni-Holm correction</td><td>Group main effect: F(3,110) = 65.60; p &lt; 0.0001</td></tr><tr><td>Figure 5j (down)</td><td>One-way ANOVA followed by Bonferroni-Holm correction</td><td>Group main effect: F(3,110) = 52.88; p &lt; 0.0001</td></tr><tr><td>Figure 5k</td><td>χ2-Test</td><td>χ2 = 7.75; p = 0.0513<br>χ2 = 18.87; p = 2.91x10<sup>-4</sup></td></tr></table>

<table><tr><td rowspan="5">Supp. Figure 1</td><td>Supp Figure 1d</td><td>Wilcoxon rank test</td><td>p = 0.18</td></tr><tr><td rowspan="2">Supp Figure 1e</td><td>Kolmogorov-Smirnov 2 sample test</td><td>K = 0.23; p = 0.36<br>K = 1; p &lt; 0.0001</td></tr><tr><td>Two-way RM ANOVA followed by Bonferroni-Holm correction</td><td>Time main effect: F(23,77) = 35.49; p &lt; 0.0001<br>Group main effect : F(3,77) = 15.89; p &lt; 0.0001<br>Interaction Time x Group : F(69,77) = 4.76; p &lt; 0.0001</td></tr><tr><td>Supp Figure 2a (left)</td><td>One-way ANOVA followed by Bonferroni-Holm correction</td><td>Group main effect: F(3,80) = 15.89; p &lt; 0.0001</td></tr><tr><td>Supp Figure 2a (middle-left)</td><td>Two-way RM ANOVA</td><td>Time main effect: F(23,77) = 5.26; p &lt; 0.0001</td></tr></table>
<table><tr><td rowspan="4">Supp. Figure 2</td><td>Supp Figure 2a (middle-right)</td><td>followed by Bonferroni-Holm correction</td><td>Group main effect : F(3,77) = 3.77 ; p = 0.0139<br>Interaction Time x Group : F(69,77) = 1.85 ; p &lt; 0.0001</td></tr><tr><td>Supp Figure 2a (right)</td><td>One-way ANOVA followed by Bonferroni-Holm correction</td><td>Group main effect: F(3,80) = 3.77; p = 0.0139</td></tr><tr><td>Supp Figure 2b</td><td>Pearson correlation test</td><td>see Supp Figure 2b</td></tr><tr><td>Supp Figure 2c</td><td>Two-way RM ANOVA followed by Bonferroni-Holm correction</td><td>Time main effect: F(23,77) = 19.31; p &lt; 0.0001<br>Group main effect : F(3,77) = 7.86 ; p = 0.0001<br>Interaction Time x Group : F(69,77) = 1.87 ; p &lt; 0.0001</td></tr><tr><td rowspan="2">Supp. Figure 3</td><td>Supp Figure 2d</td><td>Two-way RM ANOVA followed by Bonferroni-Holm correction</td><td>Time main effect: F(23,77) = 22.94 ; p &lt; 0.0001<br>Group main effect : F(3,77) = 5.40 ; p = 0.0023<br>Interaction Time x Group : F(69,77) = 1.88 ; p &lt; 0.0001</td></tr><tr><td>Supp Figure 2e</td><td>Kolmogorov-Smirnov 2 sample test</td><td>K = 0.33; p = 0.0140<br>K = 0.83; p &lt; 0.0001</td></tr><tr><td rowspan="4">Supp. Figure 6</td><td>Supp Figure 3a</td><td>Wilcoxon rank test</td><td>p = 0.0009</td></tr><tr><td>Supp Figure 3b</td><td>Wilcoxon rank test</td><td>p = 0.0009</td></tr><tr><td>Supp Figure 3c</td><td>One-way ANOVA</td><td>p &gt; 0.05</td></tr><tr><td>Supp Figure 6c</td><td>Wilcoxon rank test</td><td>p = 0.0238</td></tr><tr><td></td><td>Supp Figure 6d</td><td>Wilcoxon rank test</td><td>p = 0.49</td></tr><tr><td></td><td>Supp Figure 6e</td><td>χ2-Test</td><td>χ2 = 1.97; p = 0.58</td></tr></table>