# Higher-order interactions in neuronal function: From genes to ionic currents in biophysical models

Maria Reva, Alexis Arnaudon, Mickael Zbili, Abdullah Makkeh, Henry

Markram, Jean-Marc Goaillard, Werner van Geit

► To cite this version:

Maria Reva, Alexis Arnaudon, Mickael Zbili, Abdullah Makkeh, Henry Markram, et al. Higher-order interactions in neuronal function: From genes to ionic currents in biophysical models. Proceedings of the National Academy of Sciences of the United States of America, 2025, 122 (40), 10.1073/pnas.2500048122. hal-05289784

HAL Id: hal-05289784

https://amu.hal.science/hal-05289784v1

Submitted on 30 Sep 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
<figure><img src="image_2.png" /></figure>

# Higher-order interactions in neuronal function: From genes to ionic currents in biophysical models

Maria Reva<sup>a,1</sup> <sup>1</sup> ID, Alexis Arnaudon<sup>a</sup>, Mickael Zbili<sup>a,2</sup> <sup>1</sup> ID, Abdullah Makkeh<sup>b,c</sup> <sup>1</sup> ID, Henry Markram<sup>a,d</sup>, Jean-Marc Goaillard<sup>e,3</sup>, and Werner Van Geit<sup>a,3,4</sup> <sup>1</sup> ID

Affiliations are included on p. 11.

Edited by Eve Marder, Brandeis University Department of Biology, Waltham, MA; received January 15, 2025; accepted August 25, 2025

Neuronal firing patterns are the consequence of precise variations in neuronal membrane potential, which are themselves shaped by multiple ionic currents. In this study, we use biophysical models, statistical methods, and information theory to explore the interaction between these ionic currents and neuron electrophysiological phenotype. We created numerous electrical models with diverse firing patterns. By analyzing these models, we identified intricate relationships between model parameters and electrical features. Our findings show that neuronal activity is often influenced by multiple biophysical model parameters, in a nonadditive (i.e. synergistic) fashion. When comparing this with single-cell RNAseq data, we found a contrasting structure: gene expression profiles were dominated by redundancy, reflecting differences in regulatory constraints and sampling diversity. This research sheds light on the complex links between biophysical parameters and neuronal phenotypes.

cellular variability | high-order interactions | biophysical detailed neurons models | gene expression

While neurons are often categorized under broad classes, a number of studies have demonstrated that, even in seemingly homogeneous neuronal populations, neurons exhibit profound heterogeneity across multiple dimensions, ranging from genetic and molecular differences to variations in electrophysiological properties (1–5). Thus neuronal diversity appears as an intricate and multifaceted phenomenon, involving not only variations between neuronal types but also within neuronal types.

Genetic heterogeneity, for instance, is reflected in the differential expression of genes, particularly those involved in the regulation of ion channels. Interestingly, the landscape of ion channel genetic expression does not seem to be random. Several studies found correlations between the levels of expression of various ion channels (6–9), the patterns of these correlations giving rise to coexpression modules. Moreover, these coexpression modules appear to be specific of neuronal types. This suggests that, despite the variability in ion channel expression, neuronal phenotype might be defined by a specific coexpression module rather than by individual levels of expression of ion channels. A study utilizing single-cell multiplex RT-PCR on neocortical neurons identified 26 ion channel genes clustering with three calcium-binding proteins—calbindin, parvalbumin, and calretinin—suggesting a potential link between gene expression profiles and calcium dynamics (10). In fact, single-cell transcriptomics has consistently uncovered significant variations in mRNA and noncoding RNA levels between individual neurons of the same neuronal type, emphasizing the significance of molecular heterogeneity. Electrophysiological diversity adds yet another layer of complexity to neuronal phenotype. The combination of ion channel gene expression, along with the factors governing mRNA transcription, protein translation, and posttranslational modifications, plays a crucial role in shaping the electrical phenotype of a neuron (11), such that variations in electrophysiological behaviors across neurons are often associated with genetic and molecular variations (12–15).

To decipher this complex link, tools like information theory offer valuable insights by helping researchers analyze large-scale data and distinguish meaningful heterogeneity from noise. For example, the Partial Information Decomposition for Cell regulatory networks enables the identification of gene regulatory networks from single-cell data, moving beyond pairwise correlations to uncover higher-order relationships (16). Higher-order interactions in complex systems such as gene regulatory networks can exhibit both synergistic and redundant characteristics. Synergistic interactions occur when the combined effect of multiple components is greater than the sum of their individual effects. This means that certain information is only accessible when considering the joint state of a group of components, rather than any single component alone. In contrast,

## Significance

How neurons acquire their electrical identities is a central question in neuroscience. This study shows that combinations of ion channels interact in complex, high-dimensional ways to shape neuronal activity. By combining neuronal models with single-cell gene expression data, we reveal that the balance between synergy and redundancy depends on the statistical structure of underlying variability. Synergy emerges in diverse, decorrelated ensembles, while redundancy dominates when gene expression or model parameters are tightly coregulated or form subclusters. Notably, a small subset of ion channel genes captures essential variability across neuronal subtypes, highlighting how compact molecular programs can give rise to diverse electrical identities. These findings link statistical organization to neuronal function, offering insight into what makes each neuron type unique.

The authors declare no competing interest.

This article is a PNAS Direct Submission.

Copyright © 2025 the Author(s). Published by PNAS. This open access article is distributed under Creative Commons Attribution License 4.0 (CC BY).

¹To whom correspondence may be addressed. Email: mreva.compneuro@gmail.com.

²Present address: Foundation for Research on Information Technologies in Society, Zurich CH-8004, Switzerland.

³J.-M.G. and W.V.G. contributed equally to this work.

⁴Present address: Université Clermont-Auvergne, Centre hospitalier universitaire Clermont-Ferrand, Inserm, Neuro-Dol, Clermont-Ferrand F-63000, France.

This article contains supporting information online at https://www.pnas.org/lookup/suppl/doi:10.1073/pnas.2500048122/-/DCSupplemental.

Published September 29, 2025.
redundant interactions involve multiple components providing overlapping or duplicate information (such as strongly covarying variables). Synergistic and redundant interactions between genes play crucial roles in biological systems, influencing phenotype and drug response. Studies have demonstrated synergistic interactions can help elucidate cooperative dynamics within gene networks and predict nonadditive effects observed in drug combinations, while redundancy in gene copies can provide advantages under conditions promoting rapid growth and translation, highlighting the complex interplay between these types of interactions in genetic networks (17, 18). One example of synergy is cooperativity, where system components (e.g., channel conductances, gene expression) influence excitability through interdependent, often multiplicative effects. While cooperativity can enhance sensitivity, it may also introduce fragility. In contrast, degeneracy describes the system's capacity to achieve a similar pattern of activity through multiple distinct combinations of underlying components. This implies a high-dimensional solution manifold, where both synergistic and redundant interactions can coexist. Complementing this, degeneracy can also arise from partial functional redundancy among system components. Each of these modes of interaction (degeneracy, synergy, including cooperativity) supports robustness through diverse, distributed, and compensatory mechanisms.

The focus of our study is to investigate how high-dimensional interactions among biophysical model parameters and gene expression patterns relate to neuronal electrical phenotypes. We apply a unified information-theoretic framework to both conductance-based neuron models and single-cell transcriptomic data to quantify higher-order dependencies that shape excitability and firing behavior. Building on our previous work (19, 20), we generate populations of excitatory and inhibitory cortical

neuron models that replicate similar electrophysiological features despite varying in their underlying parameters. This enables us to assess how combinations of ionic conductances contribute to specific electrical properties. To quantify these dependencies, we apply Organization information (O-info) (21) to assess global dependencies among parameters, and the Redundancy-Synergy Index (RSI) (22) to relate these interactions to specific electro-physiological features. Additionally, we analyze publicly available Patch-seq and single-cell mRNA sequencing data from cortical interneuron and pyramidal neuron populations (12, 23) to investigate high-order structure in gene expression and compare it with the interaction patterns observed in model parameters. We hypothesize that the population of conductance-based models generated through stochastic sampling can serve as a "null model" for biological variability. In this framework, deviations from the "null model" may reflect biological constraints, such as transcriptional coregulation of ion channels, or methodological factors, such as sampling strategy. This perspective underscores the complexity and variability of the relationships between gene expression and neuronal function, and highlights the importance of comparing observed structure against appropriately sampled generative baselines.

## Results

**Synergistic High-Order Interactions of Ion Currents in Biophysically Detailed Models.** To investigate the presence and impact of high-order interactions on the electrical phenotype (e-phenotype) of neurons, we used a detailed biophysical model of Layer 5 Pyramidal Cells (L5PC) of the somatosensory cortex (SSCx) (Fig. 1A) (19). This model, incorporating region-specific distributions of ion channels and calcium dynamics across

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

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

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

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

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

**Fig. 1.** Analysis of a biophysically detailed L5PC neuron model and parameter interaction. (A) Morphology of the L5PC neuron model, showing the distribution of ion channels across soma, dendrites, AIS, and axon. (B) Parameter distributions from models generated using MCMC. (C) MI map between parameter pairs. The colors of the nodes correspond to the color of the compartments of the neural morphology from (A and B). The color of the edges corresponds to the MI value. (D) O-info plot showing the relationship between tuple number (number of parameters) and O-info values. Each dot represents the O-info value for a specific tuple. The bars and whiskers illustrate the mean and SD of O-info values at each tuple order. The color encodes positive (green) and negative (purple) O-info values. Shaded regions indicate the 95% CIs computed via bootstrapping (50 repetitions, sample size = 1,000). Dashed lines represent the maximum observed positive and negative O-info values obtained from 100 shuffled surrogates of the original data. (E) Normalized frequency map of parameters frequently appearing in top 5% of high-synergy tuples.
the soma, dendrites, and axon (present in the model as AIS followed by a myelinated compartment), underwent parameter exploration using a Markov Chain Monte Carlo (MCMC) approach (20). Through this process, we generated 110128 models capturing 60 electrical features of L5PC, presenting a wide range of ion channel conductances, calcium dynamics, and other biophysical properties. This process created a population of parameter sets representing a single e-phenotype of L5PC neuron. Importantly, this population of models display a variability in electrical features similar to the biological diversity observed in L5PC biological population (20).

The resulting parameter distributions revealed substantial variability in conductances such as gIh and gK_Pst, as well as in the parameters reflecting calcium decay and free calcium levels across neuronal compartments (Fig. 1B). This parameter variability underscores both variability in e-phenotype but also the degeneracy of the parameter space. Moreover, the parameter space was homogeneous and independent of the model fitness (SI Appendix, Fig. S1A). Next, we examined the relationships between these parameters by calculating pairwise mutual information (MI), which revealed significant interdependencies (Fig. 1C). Notable interactions included axonal sodium conductance (gNaT in axon) and its decay constant along the apical trunk, correlations between axonal and somatic sodium conductances, and strong interactions between calcium conductances (gCa_HVA, gCa_LVA) and intracellular calcium levels. Overall, 17 out of 30 parameters demonstrated statistically significant MI (Materials and Methods), indicating that the e-phenotype is supported by a degenerate solution space in which extensive parameter interdependencies emerge as signatures of the system's flexibility.

To move beyond pairwise interactions, we applied high-order information measures to characterize the presence and nature of multiparameter interactions (21). We first used O-info, a measure able to reveal redundant (O-info > 0) or synergistic (O-info < 0) interactions between variables. We computed O-info for combinations of 3 parameters, 4 parameters, and so forth up to 30 parameters, thereby spanning the full range of possible multiparameter interactions (Fig. 1D). The minimum of O-info values decreased consistently as the tuple size increased, reaching a minimum for 25 parameters. Beyond this point, adding more parameters did not further change the O-info measure. The negative value of O-info indicated a strong prevalence of synergistic interactions between model parameters. Synergistic interactions, in this context, are akin to collective computations within the model, where multiple parameters interact in a nonadditive manner to create outcomes that would not be predictable by examining each parameter individually. These synergistic effects suggest that certain groups of ion channels and calcium dynamics act together to influence the models' behavior in a way that goes beyond the sum of individual influences. Notably, the O-info associated with redundancy (O-info > 0) was minimal, approximately 10 times lower than that for synergy, highlighting the dominance of synergistic interactions (Fig. 1D). These patterns contrasted with those observed in the shuffled control ensembles, where column-wise randomization resulted in O-info values fluctuating around zero. Interestingly, when restricting the analysis to higher-cost L5PC models (optimization score > 7), we observed a marked increase in redundancy and a decrease in synergy values (SI Appendix, Fig. S1C), indicating that model quality influences the balance of high-order interactions.

To identify the specific parameters that contribute most significantly to the minimum O-info values, we analyzed the composition of parameter combinations associated with the

highest levels of synergy. Specifically, we focused on the parameters that appeared most frequently in the top 5% of tuples with the lowest O-info values for each tuple order. For each tuple size (e.g., combinations of 3 parameters, 4 parameters, and so on), we calculated how often each parameter appeared in the most synergistic tuples, providing a frequency score that reflects each parameter's contribution to high-synergy interactions across different tuple orders (Fig. 1E) (for the redundancy related map see SI Appendix, Fig. S1B). Key contributors to these synergistic interactions included parameters related to calcium dynamics, such as Free_Ca2+ (in both soma and dendrites) and Ca2+ decay, which frequently appeared in the most synergistic tuples. Additionally, potassium channel conductances, specifically gK_Pst and gKv3.1, were also frequently present in top synergistic tuples. This indicates that potassium conductance is another major factor shaping the collective dynamics, potentially influencing the model's capacity to generate or shape its electrical features. Other channels, including gNaT, gCa_LVA, and gCa_HVA, showed consistent but moderate contributions to synergy. Interestingly, the parameter that appeared least frequently in the most synergistic tuples was axonal gK_Tst. This aligns with our parameter sensitivity analysis (19), which showed that removing gK_Tst from the model had minimal impact on the model's output. This finding suggests that gK_Tst plays a limited role in driving the simulated neuron model's collective, synergy-based dynamics and may not be critical for the model's overall behavior.

**Dependencies Between Electrical Features and Model Parameters in Biophysically Detailed Models.** To map out the relationships between the electrical phenotype and the underlying ionic currents, we examined the relationship between the neuronal model outputs (electrical features, or “e-features”) and the model parameters. First, we extracted e-features from voltage traces produced by each generated model (Fig. 2 A, Left). We then computed MI between each pair of e-feature and biophysical parameters to quantify the dependencies. All e-features showed statistically significant MI values, ranging from 0.01 to 1. For further analysis, we focused on e-feature/parameter pairs with MI greater than 0.1 (Fig. 2 A, Right).

High MI values revealed distinct relationships between e-features and specific parameters. Calcium concentrations in the axon and backpropagating action potential (bAP) features in the basal apical dendrites displayed strong associations with calcium-related parameters in the respective compartments. Somatic and axonal sodium conductances each exhibited dependencies with rheobase (I_Threshold), while action potential (AP) amplitude and afterhyperpolarization depth were influenced by axonal sodium conductances. Mean spiking frequency showed a strong relationship with axonal calcium decay, whereas passive conductances, including g pas, shared high MI with the time to first spike (the feature used is the inverse of time to first spike 1/TFS). Input resistance (Rin) and resting membrane potential (RMP) both displayed high MI with g pas (Fig. 2 A, Right).

In order to determine whether specific sets of parameters collectively determine specific e-features, we moved beyond pairwise interactions and explored higher-order interactions between parameters and specific e-features. We focused on two e-features: RMP and the first interspike interval (ISI). Intuitively, these e-features are expected to involve distinct parameter sets due to their different nature: RMP should mainly rely on passive conductances while ISI is influenced by active conductances. Additionally, while RMP shared high pairwise MI with two distinct parameters, the first ISI did not have
<figure><img src="image_8.png" /><figcaption>**Fig. 3.** High-order interaction in the models of bNAC and cAC interneurons. (A) Morphology of the interneuron model used for optimization, illustrating the distribution of ion channels across the soma, dendrites, AIS, and axon used in the model. (B) Top: Example voltage traces from the continuous accommodating (cAC) and bursting nonaccommodating (bNAC) models (in response to IDrest protocol, 150%/140% from the rheobase). Bottom: Distribution of optimized parameters, with values for cAC shown in blue and bNAC in orange. (C) O-info values as a function of tuple order, showing the progression of high-order interactions for bNAC (Right) and cAC (Left) models. Shaded regions indicate the 95% CIs computed via bootstrapping. Dashed lines represent extrema values of O-info for shuffled data. (D) Left: Mean firing rate distribution for bNAC (orange) and cAC (blue). Middle: RSI plot for bNAC, showing the minimum parameter set driving high synergy (with colors corresponding to their anatomical location in the neuron model). Right: RSI plot and minimum parameter set for cAC, similar to the Middle panel. (E) Similar analysis as in (D) for the inverse First ISI feature.</figcaption></figure>

synergy. For the bNAC model, the minimum tuple order that achieved the most negative RSI (indicating maximum synergy) was 5 (4 parameters), whereas for the cAC model, it was 8 (7 parameters). In the bNAC model, axonal LVA and HVA calcium channels, SK channels, and calcium dynamics formed the minimum parameter set that produced maximum synergy. For the cAC model, the minimum synergistic module included the same parameters as for bNAC, along with somatodendritic LVA calcium channels, axonal Kv3.1 channels, and free calcium. Interestingly, at a tuple order of 6, both models shared the same core parameters, indicating close similarities between their underlying interaction structures.

To examine whether this similarity extends to other features, we tested whether the covariation modules for bNAC would

be a subset of those for cAC for the inverse first interspike interval (ISI). Consistent with their electrical phenotypes, the inverse first ISI was larger for bNAC (29.9 ± 10.3) than for cAC (9.6 ± 6.7) (Fig. 3 E, Right). For the bNAC model, we observed a steep and prolonged decay in RSI values, requiring 12 parameters to reach the minimum RSI, indicating significant synergy (Fig. 3 E, Middle). In contrast, only 7 parameters were needed to achieve minimal RSI values in the cAC model (Fig. 3 E, Left). Notably, only a quarter of the parameter sets overlapped between bNAC and cAC for this feature, indicating that the biophysical modules are highly specific to the electrical phenotype.

In summary, cAC and bNAC models, built from the same parameter space and morphology, exhibit nearly identical parameter distributions. Both display a predominance of synergy, highlighting the importance of cooperative parameter interactions in shaping distinct firing phenotypes. While cAC and bNAC models share similarities in their synergistic parameter modules for certain e-features, these modules are distinct and phenotype-specific for other e-features. The specificity of biophysical modules defining specific e-features highlights that the complexity of electrical phenotypes emerges from unique sets of synergistic interactions, tailored to the firing behavior characteristic of each neuronal type.

**Ion Channel mRNAs As Key Elements Identifying Neuronal Types.** To investigate high-order interactions beyond in silico model populations, we leveraged patch-seq data (12) to explore the relationship between mRNA expression levels and electrical features in biological cortical GABAergic interneurons. Since our in silico results focused on interactions between ion currents, our first step was to determine the role of ion channel expression levels in defining principal interneuron subtypes.

We began by isolating ion channel genes (a set of 150 among 45 768 mRNA) from the gene expression data (12) and performed clustering, aiming to capture the main interneuron types. When we compared our clustering results with the ground truth classifications provided by the dataset (Fig. 4A), we found that most primary clusters were well replicated using only these 150 ion channel genes. Notably, we observed distinct clusters for PV, SST, and LAMP5, indicating that ion channel gene expression alone effectively identifies these major interneuron subtypes. However, the clustering of ion channel genes merged SNCG and VIP into a single cluster, and an additional two clusters emerged—one corresponding to a mixed population of the primary cell types and the other to a distinct SST subtype.

To test whether this clustering was specific to ion channel genes, we repeated the clustering process with randomly selected mRNA subsets. We randomly sampled 150 genes 500 times and used these samples for clustering (Fig. 4B), comparing the Adjusted Rand Index (ARI) of the resulting clusters with the ground truth (Fig. 4C). The ARI of the ion channel gene clustering was 0.69, at least seven times higher than the average ARI obtained from random gene sets (mean ARI of 0.05; gray violin plot) (Fig. 4D). This difference suggests that ion channel genes have substantial predictive power for interneuron subtypes compared to random mRNA subsets.

An additional validation was performed by examining differentially expressed genes (DEGs) (Materials and Methods) to control if a more specific subset of mRNAs could match the clustering power of ion channel genes. Using the ground truth, we identified DEGs across cell types and noticed that only a subset of ion channel mRNAs was comprised in these DEGs
<figure><img src="image_9.png" /><figcaption>**Fig. 4.** Transcriptomic clustering of interneurons using Ion channel genes vs. random gene sets and differential expression analysis. (A) UMAP representation of IN transcriptomic clustering based on ion channel-encoding genes. Each point represents a single cell. **Left:** clusters labeled according to ground truth cell types for comparison. **Right:** same UMAP, with clusters identified by ion channel gene expression, each cell is labeled with a unique identifier. (B) UMAP representation of interneuron transcriptomic clustering based on a randomly selected gene subset. **Left:** cells labeled according to ground truth cell types. **Right:** cells labeled according to clustering results. (C) Venn diagram showing the overlap of DEGs identified through ground truth clustering and ion channel encoding genes. (D) Distribution of the ARI values for clustering based on randomly selected genes. Gray distribution represents ARI values for clustering with random genes sampled from the entire gene pool, while blue distribution represents ARI values for clustering with random genes from the DEGs. Each dot corresponds to a single random sampling of the corresponding condition. The ARI value for clustering based on ion channel-encoding genes is marked by a red dot.</figcaption></figure>

(Fig. 4C). We then sampled 150 DEGs 500 times and performed clustering. The average ARI for DEGs (mean ARI of 0.4; blue violin plot) was higher than the ARI for randomly sampled genes from the full gene pool (mean ARI of 0.05), showing that DEGs generally have better clustering performance. Importantly, the ARI for ion channel genes was within the top 1% of the ARI distribution generated by DEG clusters, confirming that ion channel genes have a predictive power greater than 99% of DEG subsets (Fig. 4D).

In summary, these findings demonstrate that ion channel gene expression alone provides a robust basis for distinguishing interneuron subtypes, outperforming both random mRNA samples and most DEG subsets. This result emphasizes the unique and predictive role of ion channels in defining neuronal identity and capturing the main structural and functional variations among interneurons.

**High-Order Interactions Between Ion Channel mRNAs Specific to the Electrical Phenotype.** To investigate whether ion channel genes share statistical interactions specific to neuronal subtypes, we focused on two electrically different interneuron types: Lamp5 and Pvalb. Pvalb interneurons are characterized by a fast-spiking pattern of activity, while Lamp5 neurons exhibit slower, delayed spiking (12). We selected 50 ion channel genes with nonzero median expression in both subtypes, including 2 HCN channels, 6 calcium channels and subunits, 6 chloride channels and subunits, 4 sodium channels, and 32 potassium channels.

First, we analyzed pairwise interactions between these ion channel genes using mutual information (Fig. 5A). The Lamp5 neurons exhibited higher MI values, reaching up to 0.26, compared to Pvalb neurons, which peaked at 0.19. In total, 58% of the ion channels in Lamp5 displayed statistically significant pairwise interactions, whereas only 34% of ion channels in Pvalb did so. In Lamp5 neurons, these interactions involved a broad range of potassium channels (Kcnq, Kcnk, Kcnj, Kcnt, Kcnc, Kcnn, and Kncd families), as well as calcium channels

(Cacna1 family), chloride channels (Clcn family), and sodium channels (Scn family). In contrast, the Pvalb neurons exhibited a more focused pattern, with strong correlations primarily among potassium channels from the Kcn family, along with select sodium channels from the Scn family.

We then examined high-order interactions among ion channel genes for both cell types. Both cell types showed a strong prevalence of redundant interactions (positive O-info), as indicated by negative O-info values fluctuating around zero. The O-info increased monotonically with tuple size, with potassium-encoding genes dominating the tuples that exhibited the highest redundancy (SI Appendix, Fig. S3).

To explore the link between electrical features and ion channel gene expression, we analyzed electrical features from patch-seq data alongside gene expression profiles from the same individual cells for both Lamp5 and Pvalb neurons. We focused on features that distinguish these cell types, particularly action potential half-width (AP half-width), known to be shorter in Pvalb neurons (thus allowing fast-spiking). The AP width in Pvalb was threefold shorter (0.4 ± 0.1 ms) than in Lamp5 (1.1 ± 0.5 ms). Redundancy-Synergy Index (RSI) analysis again revealed predominantly redundant interactions among ion channel genes in both types. The tuple order yielding the highest RSI was similar between cell types (15 for Lamp5 and 14 for Pvalb). However, the genes driving maximal RSI for AP width differed: in Lamp5, potassium channel genes predominated, along with sodium channels and one chloride channel, while in Pvalb, we observed contributions from Hcn1, multiple calcium channels, and sodium channels. We further extended our analysis on the rheobase feature and similarly found implication of distinct coexpression modules (SI Appendix, Figs. S4 and S5).

These findings suggest that, although both cell types exhibit redundant high-order gene interactions, the specific ion channel genes contributing to these interactions differ between Lamp5 and Pvalb, reflecting their distinct electrophysiological properties. Moreover, while in silico populations displayed
<figure><img src="image_10.png" /><figcaption>**Fig. 5.** High-order interactions in the Ion Gene expression data. (A) A map of 50 ion channel-encoding genes. Node color represents the general current type: potassium (coral), sodium (blue), chloride (green), calcium (yellow), and HCN (purple) families. The color bar indicates the MI values. *Left*: MI map for Lamp5 ion channel gene expression. *Right*: same for Pvalb. (B) O-info values plotted as a function of tuple order for Lamp5 (*Left*) and Pvalb (*Right*) neurons. Shaded regions indicate the 95% CIs computed via bootstrapping. Dashed lines represent extrema values of O-info for shuffled data. (C) Distribution of action potential (AP) widths (in milliseconds) for cells from Lamp5 and Pvalb neuron types. (D) RSI values as a function of tuple number (*Right*) and corresponding set of the genes that result in highest value of RSI (*Left*) for Lamp5 (*Top*) and Pvalb (*Bottom*) neurons.</figcaption></figure>

mostly synergistic interactions between ion conductances, biological populations displayed strong redundant interactions between ion channel genes.

**Bridging the Gap Between the Models and Transcriptomics Results.** While our initial analyses revealed a dominance of synergistic interactions in uniformly sampled conductance-based models, and a contrasting dominance of redundancy in gene expression data, the origin of this discrepancy remained unclear. Thus we investigated whether this observed discrepancy can be explained by the covariance structure or by the sampling strategy of the parameter space. To address this, we employed two complementary strategies: 1) reanalyzing a dataset of biophysical models fitted to pyramidal neurons using a different sampling strategy (26), and 2) subsampling high-covariance models from a uniformly sampled interneuron population.

In the first approach, we asked whether a shift toward redundancy could be observed in models generated through an entirely different sampling strategy, specifically, one that

reflects targeted optimization around biological data rather than stochastic exploration. To this end, we analyzed a previously published dataset of conductance-based models fitted to L4/L5 pyramidal neurons from the mouse visual cortex using a multiobjective evolutionary algorithm (IBEA) (26, 27). Although these models were not generated in the present study, they offer a valuable point of comparison, as they exemplify a data-constrained, nonrandom modeling approach in which multiple solutions are clustered around fixed electrophysiological targets.

To ensure comparability, we selected only those neurons whose electrophysiological properties were tightly clustered into a single electrical type (as classified in ref. 28), yielding a subset of 72 neurons. For each of these neurons, a population of 40 models had been generated, each representing a distinct solution from the final generation of the evolutionary optimization, known as the Pareto-optimal “Hall of Fame” (HOF) set. These 40 solutions per neuron are all consistent with the target electrophysiological features, but differ in their underlying ion channel parameters.

We first analyzed the full set of 2,880 HOF models (72 neurons × 40 models) and computed O-info across increasing tuple orders. Despite all models matching their electrophysiological targets, the interaction profiles exhibited a clear emergence of both redundancy and synergy, with redundancy dominating at higher tuple orders (Fig. 6A).

We next examined how sampling strategy affects the prevalence of synergy and redundancy. Specifically, we compared two equally sized subsets of 72 models: 1) a “Unique group” comprising one randomly selected HOF model from each of the 72 neurons, and 2) an “Oversampled group” comprising 36 models from each of two neurons. The Unique group exhibited synergy-dominated interaction profiles, whereas the Oversampled group showed pronounced redundancy (Fig. 6B). This effect was not due to sample size alone: both groups contained the same number of models. Rather, it reflected the structure of sampling: the Unique group combined models from diverse clusters, forming a broad distribution, whereas the Oversampled group isolated two tight clusters, where constrained parameter relationships drove redundancy. Dimensionality reduction of the full model ensemble further illustrates this point: PCA reveals that models from the same optimization run (i.e., 40 HOF models per cell) form tight, cell-specific clusters, while models from different cells are more broadly dispersed (Fig. 6G).

To quantify these differences, we computed the absolute area under the curve (AUC) for O-info trajectories across tuple orders. The Unique group showed significantly higher synergy AUC than redundancy (*P* < 0.001), while the Oversampled group showed the opposite pattern (*P* < 0.001) (Fig. 6C). Both groups significantly differed from their respective shuffled controls (SI Appendix, Fig. S6C), confirming that the observed effects were not due to sampling noise. These findings suggest that oversampling parameter sets optimized to individual neurons creates local parameter-space clusters, leading to redundancy-dominated structures. In contrast, broader sampling across diverse neurons avoids these clusters, allowing the emergence of synergy.

We then asked whether similar principles apply to biological gene expression. Using single-cell transcriptomic data from the mouse visual cortex (23), we analyzed 466 excitatory neurons from layers 4 and 5 and focused on ion channel-encoding genes. O-info analysis of expression profiles revealed a clear dominance of redundancy across all interaction orders (Fig. 6D), closely paralleling the redundancy observed in oversampled biophysical model populations. Although these transcriptomic data originate
<figure><img src="image_11.png" /><figcaption>**Fig. 6.** Model parameter and gene expression space shapes high-order interactions. (A) O-info across tuple orders for all 2,880 models. Green (positive) and purple (negative) values indicate redundancy- or synergy-dominated interactions, respectively; lines show mean ± SD. Shaded areas denote 95% CIs from bootstrapping. (B) O-info for two model groups: “Unique” (72 models from distinct HOF sets) and “Oversampled” (72 models from two cells, 36 each). Gray and blue dashed lines represent control groups with parameter-wise shuffled models. (C) Distribution of absolute O-info AUC values (minima for synergy, maxima for redundancy) in Unique vs. Oversampled groups. Violin plots show variability; Mann-Whitney *U* test, *P* < 0.001. (D) O-info across tuple orders for all L4/5 pyramidal neuron models from VISp. (E) O-info across tuple orders restricted to models from the L4 Sparcl1+; curves as in (D). (F) O-info across tuple order for a combined ensemble of four excitatory subtypes (L4 Sparcl1+, L5 Chrna6+, L5 Samd3+, and L5 Cdh13+). (G) PCA projection of the 2,880 models’ parameter space [from panel (A)], with five randomly selected cells highlighted in color (one dot = one HOF model). The gray background shows the distribution of other models. (H) PCA projection of the ion gene expression of the cells shown in panel (D). Each dot represents one cell, with color indicating its excitatory subtype cluster.</figcaption></figure>

from a separate study and are not directly linked to the modeled neurons, they sample from the same cortical layers (L4/L5 VISp) and offer a meaningful point of comparison.

To probe subtype-specific differences, we examined four well-defined excitatory neuron subtypes (23): L4 Sparcl1+ ($n = 75$), L5 Chrna6+ ($n = 7$), L5 Samd3+ ($n = 25$), and L5 Cdh13+ ($n = 28$). When analyzing the L4 Sparcl1+ subtype alone, we found a surprisingly strong synergy signature (Fig. 6E), suggesting that within-subtype variability may encode complementary, high-order information. Importantly, this analysis was conducted at the single-cell level and the observed synergy reflects intrinsic high-order relationships present in the full single-cell ensemble of the Sparcl1+ population. However, when we combined the L4 and L5 subtypes into a pooled ensemble, redundancy became dominant (Fig. 6F). This shift echoes the pattern observed in oversampled model sets and suggests that mixing neuronal subtypes can dilute synergistic interactions and promote redundancy. It is further supported by the low dimensional projections of the combined transcriptomic dataset showing partial clustering by subtype (Fig. 6H), confirming that the underlying gene expression distributions retain subtype-specific structure even after pooling. Notably, this effect was not observed when we combined two closely related L4 subtypes [Sparcl1+ and Scnn1a+ ($n = 95$)], which appear intermixed in low-dimensional space; in this case, synergy remained slightly dominant rather than giving way to redundancy (SI Appendix, Fig. S6 A and B). Together, these results demonstrate that both synergy and redundancy can emerge in biophysical and biological systems, depending on how sampling captures the structure of underlying variability. Redundancy prevails when sampling isolates tightly clustered solutions—whether through repeated optimization around a fixed target or by focusing on biological subtypes. In contrast, synergy emerges when sampling

spans across distinct clusters, revealing broader, less constrained parameter relationships.

To further examine how structural properties of a dataset influence high-order interactions, we turned to our uniformly sampled cAC interneuron models, which previously exhibited strong synergy (Fig. 3 *B* and *C*). Unlike the Nandi et al. dataset or transcriptomic subtype comparisons, these models were not clustered around specific targets and served as an ideal background to isolate statistical effects. We asked whether it would be possible to reveal redundancy within this synergy-rich population purely by selecting models with high interparameter covariance. To do so, we selected a subset of models with high interparameter covariance, using PCA and a covariance-based scoring metric (*Materials and Methods*). Importantly, this selection preserved both the electrophysiological targets and the model structure—only the statistical dependencies among parameters were altered.

Strikingly, despite being drawn from a synergy-rich ensemble, the high-covariance subset showed a clear shift toward redundancy (Fig. 7 *A*, *Middle* row), in stark contrast to equivalently sized random subsets, which retained synergistic interactions (Fig. 7 *A*, *Top* row). This divergence was evident both in the O-info trajectories and in the MI matrices (Fig. 7 *B*): While the random subset showed minimal pairwise dependencies (Fig. 7 *B*, *Top* row), the high-covariance subset revealed a pronounced off-diagonal structure (Fig. 7 *B*, *Middle* row), indicative of widespread parameter covariances. A similar pattern of prevalence of the off-diagonal structure can be observed in the MI matrix of the ionic gene expression of Lamp5 (Fig. 7 *B*, *Bottom* row).

To visualize these dependencies, we represented the 4D relationship among four example parameters (gNaT, gSk, gKv3.1, glh) within each group (Fig. 7 *C*). In the randomly sampled subset, the joint distribution appeared diffuse and unstructured,
A

Randomly sampled subset (models)

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

B

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

C

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

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

Highly covariant subset (models)

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

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

Gene expression (Lamp5)

<figure><img src="image_18.png" /><figcaption>**Fig. 7.** Covariance structure and the balance between synergy and redundancy in biophysical models and gene expression. (A) O-info trajectories across tuple orders for three distinct datasets: a randomly sampled (n = 1,000) cAC model subset (Top), a highly covariant model cAC subset (Middle, n = 1,000), and Lamp5 gene expression data (Bottom). Positive values (green) indicate redundancy-dominated interactions; negative values (purple) indicate synergy. Shaded regions show the full bootstrap range (n = 50 repetitions of 100 samples each for cAC models, repetitions of 240 for Lamp5), and lines with error bars denote the mean ± SD (B) Pairwise mutual information (MI) matrices for the same parameter or gene sets shown in (A). Only the lower triangle is shown, with parameter names rotated and placed along the diagonal. Color of the parameter names correspond to their anatomical location in the neuron model as in (Fig. 3A). (C) Example 4D scatterplots highlighting relationships among four parameters or genes from each subset. In the Top two panels (model subsets), the same set of four parameters is shown (e.g., gSK, gNaT, gKv3.1, glh), with color representing one dimension. In the Bottom panel (gene expression), four ion channel genes (e.g., Hcn1, Kcnn2, Scn1a, Kcnh7) are shown for a Lamp5 interneuron population. All values were normalized and vary from 0 to 1.</figcaption></figure>

consistent with high-dimensional independence. In contrast, the high-covariance subset revealed a clear geometric organization, suggestive of tightly coupled parameter relationships. Interestingly, this structured pattern was very similar to the one observed in gene expression data from Lamp5 interneurons (Fig. 7 C, Bottom), where a selected set of ion channel genes (Scn1a, Kcnh7, Kcnn2, Hcn1) showed clear interdependence in both the MI matrix and the 4D projection.

In summary, these results demonstrate that redundancy arises under two distinct structural conditions: when sampling is clustered around narrow regions of parameter or expression space and when variables exhibit strong covariation across the population. Conversely, synergy is most prominent in diverse or decorrelated ensembles. This parallel between the biophysical and transcriptomic domains suggests that the observed interaction structure reflects population-level statistical organization, rather than intrinsic properties of the individual components.

Discussion

This study explores the connection between biophysical parameters, gene expression, and neuronal electrophysiological phenotypes. Using biophysical models, statistical tools, and information theory, we examined how ionic currents shape

neuronal behavior and linked these insights to transcriptomic profiles. We generated diverse electrical models via MCMC sampling, enabling systematic analysis of how combinations of ion channel parameters influence neuronal features. Our findings show that these features rarely depend on a single current, but instead emerge from complex interactions among multiple channels. This highlights the need for multivariate tools, such as O-info and RSI, to reveal the integrated effects of conductance interactions. Furthermore, our analysis of transcriptomic data identified gene expression modules associated with specific interneuron types, reinforcing the link between molecular and electrical identity. The prevalence of high-order redundancy among ion channel genes suggests a mechanism for robustness, degeneracy, and functional compensation. Potassium channel genes were especially dominant in redundant tuples, highlighting their role in stabilizing neuronal excitability. Notably, the specific genes driving maximal redundancy differed between Lamp5 and Pvalb interneurons, reflecting their distinct electrophysiological properties. This difference was particularly evident in the action potential width, a key distinguishing feature between these two interneuron types. Pvalb neurons are characterized by significantly shorter AP widths and lower rheobase values compared to Lamp5 neurons, aligning with their fast-spiking capabilities and higher excitability (29, 30). The AP width influences the
timing precision and frequency of neuronal firing, affecting synaptic transmission and the synchronization of neural networks (31). A narrower AP width in Pvalb neurons allows for rapid and high-frequency firing, which is essential for their role in controlling the timing of cortical circuitry (29). The shorter action potential width in Pvalb neurons, indicative of their fast-spiking nature, can be attributed to the specific combination of ion channels contributing to the maximal redundancy in this cell type. It is crucial to note that mRNA expression does not directly translate into protein expression (32). Future studies incorporating proteomic data would provide a more direct link between gene expression and functional ion channels.

We would like to note that our modeling strategy uses conductance-based ensembles as a null framework to investigate high-order interactions in populations sharing a single functional phenotype. Specifically, we fit models to the mean electrophysiological features of neurons representing one canonical electrical type. This averaging approach allows us to isolate parameter interactions that are intrinsic to the phenotype, without confounding variability from individual cells or subtypes. While this strategy omits potentially informative within-cell parameter structure, explicitly fitting individual neurons introduces statistical and conceptual issues, notably artificial clustering in parameter space. As a result, our null model provides a principled baseline for interpreting synergy and redundancy. Future work may explore alternative sampling strategies (transcriptomic constraints) to better balance biological realism with statistical interpretability.

The synergistic nature of the models and the redundancy observed in gene expression patterns raise interesting questions about the evolutionary advantages of these arrangements. Redundant elements in gene expression may provide robustness to the system, allowing for functional compensation if one component fails. In E.coli it was shown that the presence of multiple rRNA/tRNA gene copies can be advantageous under conditions that promote rapid growth and translation (33). This redundancy allows organisms to adapt more quickly to favorable environmental conditions (33). At the same time, redundancy offers significant advantages for both resilience and evolution. It has been shown that biological redundancy in protein-protein interaction networks enables robustness and increased tolerance for perturbations, enhancing their ability to maintain functionality despite genetic mutations or environmental stresses (34, 35). Conversely, synergetic interactions in the model might offer more efficient information processing but at the cost of increased vulnerability if a critical component is disrupted.

In this study, we directly addressed the possible causes of redundancy and synergy in both models and biological neurons, and provide answers about the discrepancies observed in these two datasets: While synergy between biophysical parameters dominated in most models, redundancy in ion channel expression dominated in most biological neuronal types. We show that this balance is shaped by interparameter covariance: even in broadly sampled MCMC models, selecting high-covariance subsets was sufficient to shift interactions toward redundancy. This result demonstrates that parameter covariation alone, independent of biological constraints or fitting procedures, can suppress synergy by effectively reducing the dimensionality of the solution space. One possible interpretation is that biological solutions occupy a covariation-constrained subspace of a larger computationally adequate parameter space. The existence of redundancy-dominated subsets of models (Fig. 7) suggests that biological neurons may use a subspace of the theoretically possible and much larger space of gene expression fitting their electrical phenotype. A putative explanation for this fact is that covarying solutions might be

much easier to generate (for instance through transcription regulatory mechanisms targeting ion channels belonging to same coexpression modules) and might significantly decrease the expense of energy necessary to achieve a target phenotype.

Another potential source of redundancy dominance comes from our observation that models derived through repeated optimization around fixed biological targets (26) exhibit increased redundancy despite matching the same electrophysiological features. A similar redundancy-dominated structure was evident in our analysis of ion channel mRNA from excitatory neurons (23). In both cases, oversampling within local basins of the parameter or gene expression space leads to clustered solutions, further shaping the redundancy-synergy balance.

Our findings also indicate that excitability can emerge from multiple distinct configurations of underlying parameters, consistent with the principle of degeneracy. In this light, synergy observed among ion channels does not necessarily imply fragility, but may reflect this degeneracy: high-dimensional, flexible solution spaces that support stable excitability through diverse compensatory strategies. Such parameter interdependencies have been extensively described in modeling studies (36–41) and are thought to underlie homeostatic regulation in real neurons, allowing stable function despite genetic variability, developmental changes, or perturbations (42). Thus, degeneracy, as revealed by synergistic parameter combinations, may be a key design principle of robust neural excitability.

While our study demonstrates the significant predictive power of ion channel genes in distinguishing interneuron types, it is important to acknowledge that a comprehensive understanding of neuronal diversity requires consideration of synaptic genes and other groups of genes. Recent research has highlighted the strong predictive power of synaptic and GPCR genes in defining neuronal subtypes and their electrical properties (3).

This research contributes to the growing body of evidence linking mRNA expression to observable neuronal phenotypes (1, 7, 15). By integrating computational modeling, statistical analysis, and transcriptomic data, we provide a framework for understanding the complex relationships between gene expression and the electrophysiological properties that define neuronal function. This approach holds significant potential for future studies investigating the molecular basis of neuronal diversity and the pathogenesis of neurological disorders arising from disrupted neuronal signaling.

## Materials and Methods

**Detailed Biophysical Models.** We constructed conductance-based models of mouse cortical pyramidal neurons and interneurons using morphologically reconstructed cells and a set of predefined ion channel mechanisms in NEURON environment (43), as previously described and available online (19). Electrophysiological features were extracted from Patch Clamp recordings across various stimulation protocols using the eFEL library (19, 44). Each model was optimized to match the canonical electrophysiological profile-defined as the mean feature vector of a given cell type-using a z-score-based cost function and a MCMC procedure (20). Simulations were performed in NEURON, with ion channel distributions assigned by compartment and e-type. The resulting model populations included over 63,408 cAC, 86,361 bNAC, and 110,128 L5PC models that satisfied optimization thresholds (3 or 5 SDs). Full methodological details, including protocols, features, and channel distributions, are provided in SI Appendix.

**Information theory measures and computations.** We quantified high-order dependencies using the O-info (O-info) and the RSI. O-info captures net redundancy or synergy among a set of variables $X = \{X_1, ..., X_N\}$, based on their joint and marginal entropies (21). RSI measures the net interaction
among predictors S relative to a target variable Y, defined as the difference between joint and individual mutual information terms (22). Both metrics were computed under a multivariate Gaussian assumption, enabling closed-form solutions based on the covariance matrix. Full definitions and derivations are provided in SI Appendix.

**Gaussian Formula for O-Info.** For a multivariate Gaussian system with covariance matrix Σ, the O-info is approximated as

$$O(\mathbf{X}) = \frac{1}{2} \ln \left( \frac{|\Sigma|^{|N-2|} \prod_{i=1}^{N} \sigma_i^2}{\prod_{i=1}^{N} |\Sigma_{-i}|} \right),$$

where:

* $\sigma_i^2$ is the variance of the i-th variable (diagonal of Σ),

* $|\Sigma|$ is the determinant of the covariance matrix.

**Gaussian Formula for RSI.** For RSI under Gaussian assumptions, the mutual information terms $I(S; Y)$ and $I(X_i; Y)$ are calculated directly from the covariance matrices:

$$I(S; Y) = \frac{1}{2} \log \left( \frac{|\Sigma_S|}{|\Sigma_S| |Y|} \right),$$

$$I(X_i; Y) = \frac{1}{2} \log \left( \frac{\sigma_{X_i}^2}{\sigma_{X_i|Y}^2} \right).$$

Substituting into the RSI formula:

$$\text{RSI}_{\text{Gaussian}}(S; Y) = \frac{1}{2} \log \left( \frac{|\Sigma_S| \prod_{X_i \in S} \sigma_{X_i|Y}^2}{|\Sigma_S| |Y| \prod_{X_i \in S} \sigma_{X_i}^2} \right).$$

**Detection of Stabilization Points.** To identify stabilization points in O-info or RSI curves, we applied a derivative-based algorithm. For each tuple order, we computed first and second discrete derivatives of the O-info/RSI series and defined the stabilization point as the index with the steepest concavity (minimum second derivative). This marks the transition from rapid decay to plateau. Full derivative expressions and rationale are provided in SI Appendix.

**Clustering and Analysis of the Single Cell mRNA from Patch-Seq Data.** We analyzed single-cell RNA-seq data using Scanpy (45), applying standard preprocessing, log-transformation, and gene-wise standardization. Gene subsets included ion channel genes, DEGs, and 1,000 random sets of matched size. PCA was used for dimensionality reduction, followed by kNN graph construction and Leiden clustering. Cluster quality was assessed using ARI against known cell types. Full details on gene selection, normalization, and evaluation appear in SI Appendix.

**Selection of the High-Covariance Model Subset.** To assess the role of structured parameter covariation, we selected a high-covariance subset from the cAC model ensemble (cost < 3) using a PCA-based method. Models

were projected onto the top five principal components (explaining 70% of variance), reconstructed, and ranked by residual L2 error. The 1,000 models with lowest reconstruction error were retained. The complete selection procedure is described in SI Appendix.

**Mutual Information Computation and Significance Testing.** We computed pairwise MI between parameters and electrophysiological features using histogram-based discretization. For each pair, significance was assessed by shuffling one variable (1,000 permutations) and comparing to a null distribution. P-values were corrected for multiple comparisons using Benjamini-Hochberg FDR (alpha = 0.05). Only MI values passing FDR and a minimum effect size threshold (MI > 0.1) were considered significant. Details are provided in SI Appendix.

**Bootstrap Estimation of High-Order Measurements.** To estimate Cls for O-info and RSI metrics, we performed 50 bootstrap repetitions per condition. Samples of 1,000 models (for MCMC data) or 60% of cells (for transcriptomic data) were drawn without replacement. Metrics were computed across tuple orders to capture variability. Datasets with fewer than 100 samples or strongly imbalanced clusters were excluded due to resampling artifacts. Full rationale and bootstrapping parameters are given in SI Appendix.

**Data, Materials, and Software Availability.** Single-cell RNA-seq data from Tasic et al. (23) and Patch-seq data from Gouwens et al. (12) are available via the Allen Brain Map portal: https://portal.brain-map.org. Conductance-based models from Nandi et al. (26): https://github.com/AllenInstitute/All-active-Manuscript. Code for computing O-info and RSI: https://github.com/BlueBrain/emodel-generalisation. The e-models (cAC, bNAC, L5PC): https://doi.org/10.5281/zenodo.16813044. Models generated have been deposited in MCMC-Generated Neuronal Models for cAC, bNAC, and L5PC Cell Types with Features, Parameters, and Scores (DOI: 10.5281/zenodo.16813045). Previously published data were used for this work (https://knowledge.brain-map.org/data/1HEYEW7GMUKWIQW37BO/summary).

**ACKNOWLEDGMENTS.** This study was supported by funding to the Blue Brain Project, a research center of the École Polytechnique Fédérale de Lausanne, from the Swiss government's Eidgenössische Technische Hochschule Board of the Swiss Federal Institutes of Technology. We would like to thank Dr. Fernando Rosas from the University of Sussex for the initial discussions and suggestions regarding Information theory tools.

Author affiliations: aBlue Brain Project, École polytechnique fédérale de Lausanne (EPFL), Geneva 1202, Switzerland; bDepartment of Data-driven Analysis of Biological Networks, Göttingen Campus Institute for Dynamics of Biological Networks, University of Göttingen, Göttingen 37077, Germany; cComplex Systems Theory, Max Planck Institute for Dynamics and Self-Organization, Göttingen 37018, Germany; dLaboratory of Neural Microcircuitry, Brain Mind Institute, École Polytechnique Fédérale de Lausanne (EPFL), Lausanne 1015, Switzerland; and eSystems Approaches to Neural Excitability, Inst. de Neurosciences de la Timone, UMR7289 CNRS & Aix Marseille Université, Marsellie 13385, France

Author contributions: M.R., A.A., M.Z., and J.-M.G. designed research; M.R. performed research; M.R., A.A., and A.M. contributed new reagents/analytic tools; M.R. analyzed data; H.M., J.-M.G., and W.V.G. supervised; M.R. wrote original draft; and M.R., A.A., M.Z., A.M., J.-M.G., and W.V.G. revised and edited manuscript.

1. Y. Bernaerts et al., Combined statistical-mechanistic modeling links ion channel genes to physiology of cortical neuron types. bioRxiv [Preprint] (2023). https://doi.org/10.1101/2023.03.02.530774 (Accessed 20 March 2024).

2. J. Eberwine, J. Kim, Cellular deconstruction: Finding meaning in individual cell variation. Trends Cell Biol. 25, 569-578 (2015).

3. R. M. Kendrick, S. Linderman, S. F. Owen, Transcriptomically-measured gene expression predicts physiological variation across single neurons in humans and mice. bioRxiv [Preprint] (2024). https://doi.org/10.1101/2024.08.26.609746 (Accessed 5 October 2024).

4. A. M. Swensen, B. P. Bean, Robustness of burst firing in dissociated purkinje neurons with acute or long-term reductions in sodium conductance. J. Neurosci. 25, 3509-3520 (2005).

5. D. J. Schulz, J. M. Goaillard, E. Marder, Variable channel expression in identified single and electrically coupled neurons in different animals. Nat. Neurosci. 9, 356-362 (2006).

6. T. Kodama et al., Graded coexpression of ion channel, neurofilament, and synaptic genes in fast-spiking vestibular nucleus neurons. J. Neurosci. 40, 496-508 (2020).

7. M. Tapia et al., Neurotransmitter identity and electrophysiological phenotype are genetically coupled in midbrain dopaminergic neurons. Sci. Rep. 8, 13637 (2018).

8. D. J. Schulz, J. M. Goaillard, E. E. Marder, Quantitative expression profiling of identified neurons reveals cell-specific constraints on highly variable levels of gene expression. Proc. Natl. Acad. Sci. U.S.A. 104, 13187-13191 (2007).

9. J. M. Goaillard, A. L. Taylor, D. J. Schulz, E. Marder, Functional consequences of animal-to-animal variation in circuit parameters. Nat. Neurosci. 12, 1424-1430 (2009).

10. M. Toledo-Rodriguez et al., Correlation maps allow neuronal electrical properties to be predicted from single-cell gene expression profiles in rat neocortex. Cereb. Cortex 14, 1310-1327 (2004).

11. F. Buettner et al., Computational analysis of cell-to-cell heterogeneity in single-cell RNA-sequencing data reveals hidden subpopulations of cells. Nat. Biotechnol. 33, 155-160 (2015).

12. N. W. Gouwens et al., Integrated morphoelectric and transcriptomic classification of cortical gabaergic cells. Cell 183, 935-953 (2020).

13. F. Scala et al., Phenotypic variation of transcriptomic cell types in mouse motor cortex. Nature 598, 144-150 (2021).
14. A. L. Taylor, J. M. Goaillard, E. Marder, How multiple conductances determine electrophysiological properties in a multicompartment model. *J. Neurosci.* **29**, 5573-5586 (2009).

15. C. Bomkamp et al., Transcriptomic correlates of electrophysiological and morphological diversity within and across excitatory and inhibitory neuron classes. *PLoS Comput. Biol.* **15**, e1007113 (2019).

16. T. E. Chan, M. P. Stumpf, A. C. Babtie, Gene regulatory network inference from single-cell data using multivariate information measures. *Cell Syst.* **5**, 251–267 (2017).

17. J. Wildenhain et al., Prediction of synergism from chemical-genetic interactions by machine learning. *Cell Syst.* **1**, 383–395 (2015).

18. D. Anastassiou, Computational analysis of the synergy among multiple interacting genes. *Mol. Syst. Biol.* **3**, 83 (2007).

19. M. Reva et al., A universal workflow for creation, validation, and generalization of detailed neuronal models. *Patterns* **4**, 100855 (2023).

20. A. Arnaudon et al., Controlling morpho-electrophysiological variability of neurons with detailed biophysical models. *IsScience* **26**, 108222 (2023).

21. F. E. Rosas, P. A. Mediano, M. Gastpar, H. J. Jensen, Quantifying high-order interdependencies via multivariate extensions of the mutual information. *Phys. Rev. E* **100**, 032305 (2019).

22. G. Chechik et al., Group redundancy measures reveal redundancy reduction in the auditory pathway. *Adv. Neural Inf. Process. Syst.* **14**, 173–180 (2001).

23. B. Tasic et al., Adult mouse cortical cell taxonomy revealed by single cell transcriptomics. *Nat. Neurosci.* **19**, 335–346 (2016).

24. W. Hu, B. P. Bean, Differential control of axonal and somatic resting potential by voltage-dependent conductances in cortical layer 5 pyramidal neurons. *Neuron* **97**, 1315–1326 (2018).

25. H. Markram et al., Interneurons of the neocortical inhibitory system. *Nat. Rev. Neurosci.* **5**, 793–807 (2004).

26. A. Nandi et al., Single-neuron models linking electrophysiology, morphology, and transcriptomics across cortical cell types. *Cell Rep.* **40**, 111176 (2022).

27. W. Van Geit et al., Bluepyopt: Leveraging open source software and cloud infrastructure to optimise model parameters in neuroscience. *Front. Neuroinf.* **10**, 17 (2016).

28. N. W. Gouwens et al., Classification of electrophysiological and morphological neuron types in the mouse visual cortex. *Nat. Neurosci.* **22**, 1182–1195 (2019).

29. H. Hu, J. Gan, P. Jonas, Fast-spiking, parvalbumin+ gabaergic interneurons: From cellular design to microcircuit function. *Science* **345**, 1255263 (2014).

30. R. Tremblay, S. Lee, B. Rudy, Gabaergic interneurons in the neocortex: From cellular properties to circuits. *Neuron* **91**, 260–292 (2016).

31. B. P. Bean, The action potential in mammalian central neurons. *Nat. Rev. Neurosci.* **8**, 451–465 (2007).

32. J. Reimegård et al., A combined approach for single-cell mRNA and intracellular protein expression analysis. *Commun. Biol.* **4**, 624 (2021).

33. P. K. Raval, W. Y. Ngan, J. Gallie, D. Agashe, The layered costs and benefits of translational redundancy. *eLife* **12**, e81005 (2023).

34. B. Klein et al., A computational exploration of resilience and evolvability of protein-protein interaction networks. *Commun. Biol.* **4**, 1352 (2021).

35. P. Hunter, Understanding redundancy and resilience: Redundancy in life is provided by distributing functions across networks rather than back-up systems. *EMBO Rep.* **23**, e54742 (2022).

36. A. A. Prinz, D. Bucher, E. Marder, Similar network activity from disparate circuit parameters. *Nat. Neurosci.* **7**, 1345–1352 (2004).

37. J. Golowasch, M. S. Goldman, L. Abbott, E. Marder, Failure of averaging in the construction of a conductance-based neuron model. *J. Neurophysiol.* **87**, 1129–1131 (2002).

38. G. Dion, T. O'Leary, E. Marder, Ion channel degeneracy enables robust and tunable neuronal firing rates. *Proc. Natl. Acad. Sci. U.S.A.* **112**, E5361–E5370 (2015).

39. A. Fyon, A. Franci, P. Sacré, G. Dion, Dimensionality reduction of neuronal degeneracy reveals two interfering physiological mechanisms. *PNAS Nexus* **3**, pgae415 (2024).

40. A. Jain, R. Narayanan, Degeneracy in the emergence of spike-triggered average of hippocampal pyramidal neurons. *Sci. Rep.* **10**, 374 (2020).

41. R. K. Rathour, R. Narayanan, Homeostasis of functional maps in active dendrites emerges in the absence of individual channelostasis. *Proc. Natl. Acad. Sci. U.S.A.* **111**, E1787–E1796 (2014).

42. T. O'Leary, A. H. Williams, A. Franci, E. Marder, Cell types, network homeostasis, and pathological compensation from a biologically plausible ion channel expression model. *Neuron* **82**, 809–821 (2014).

43. N. T. Carnevale, M. L. Hines, *The NEURON Book* (Cambridge University Press, 2006).

44. H. Markram et al., Reconstruction and simulation of neocortical microcircuitry. *Cell* **163**, 456–492 (2015).

45. F. A. Wolf, P. Angerer, F. J. Theis, Scanpy: Large-scale single-cell gene expression data analysis. *Genome Biol.* **19**, 1–5 (2018).