# Optimal responsiveness and information flow in networks of heterogeneous neurons

Matteo Di Volo, Alain Destexhe

► To cite this version:

Matteo Di Volo, Alain Destexhe. Optimal responsiveness and information flow in networks of heterogeneous neurons. Scientific Reports, 2021, 11 (1), pp.17611. 10.1038/s41598-021-96745-2. hal-03352455

HAL Id: hal-03352455

https://hal.science/hal-03352455

Submitted on 17 Nov 2021

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>

Check for updates

OPEN

# Optimal responsiveness and information flow in networks of heterogeneous neurons

Matteo Di Volo<fnref n="1" />✉ & Alain Destexhe²

Cerebral cortex is characterized by a strong neuron-to-neuron heterogeneity, but it is unclear what consequences this may have for cortical computations, while most computational models consider networks of identical units. Here, we study network models of spiking neurons endowed with heterogeneity, that we treat independently for excitatory and inhibitory neurons. We find that heterogeneous networks are generally more responsive, with an optimal responsiveness occurring for levels of heterogeneity found experimentally in different published datasets, for both excitatory and inhibitory neurons. To investigate the underlying mechanisms, we introduce a mean-field model of heterogeneous networks. This mean-field model captures optimal responsiveness and suggests that it is related to the stability of the spontaneous asynchronous state. The mean-field model also predicts that new dynamical states can emerge from heterogeneity, a prediction which is confirmed by network simulations. Finally we show that heterogeneous networks maximise the information flow in large-scale networks, through recurrent connections. We conclude that neuronal heterogeneity confers different responsiveness to neural networks, which should be taken into account to investigate their information processing capabilities.

Studying the collective behavior of large numbers of units interacting non-linearly is a classical theme in physical and computational sciences. In biology, such studies are complicated by the fact that the units are usually non identical, but rather display considerable heterogeneity. This particularly apparent in cerebral cortex, where neuronal size and properties are highly heterogeneous<fnref n="1" />⁻¹². Neuronal heterogeneity is particularly high for inhibitory neurons, for which many cell classes were observed¹³⁻¹⁶. Recent works have shown also that heterogeneity in neural response, together with synaptic plasticity, can have a strong impact for response of single neurons as part of a network and can be functional to maintain low firing rates¹⁷⁻¹⁹. In these cases it is shown that networks of neurons can synchronize even in presence of heterogeneity²⁰. In general, in networks of oscillators, as in neuronal networks, such heterogeneity across units (bare frequencies or neurons' excitability) induces typically desynchronization at a population scale²¹⁻²⁹. As one would expect, the more the neurons are different the less they are able to synchronize and to correlate their reciprocal activity. Nevertheless, despite this heterogeneity, cortical populations are able to respond coherently to external stimuli and also to generate synchronous collective oscillations³⁰,³¹.

In the present paper, we consider sparse networks of spiking neurons endowed with heterogeneity in neurons' intrinsic excitability, that we treat independently for excitatory and inhibitory neurons. We investigate the responsiveness of such networks to external inputs, and also compare the heterogeneity to experimental estimates. To yield insights on the causes of the responsiveness properties observed numerically, we employ theoretical methods to derive a mean-field model of heterogeneous networks. We then use the mean-field model to predict the emergence of new activity states due to heterogeneity, a prediction which is tested by numerical simulations of network dynamics. Finally, to determine if the enhanced responsiveness also applies to recurrent connectivity, we considered the propagation of activity in large-scale networks of mean-field models.

## Results

We first investigate the responsiveness of heterogeneous networks of excitatory and inhibitory neurons numerically, and next we use a theoretical approach to understand the mechanisms underlying responsiveness.

<footnote n="1">¹Laboratoire de Physique Théorique et Modélisation, Université de Cergy-Pontoise, CNRS, UMR 8089, 95302 Cergy-Pontoise cedex, France. ²Paris-Saclay University, Institute of Neuroscience, CNRS, Gif sur Yvette, France. ✉email: matteo.di-volo@u-cergy.fr</footnote>
<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" /><figcaption>f</figcaption></figure>

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

**Figure 1.** Inhibitory neuron heterogeneity optimizes network responsiveness. (a) An afferent excitatory input (Poissonian spike train at a time-dependent frequency $v_{ext}(t)$, *Input*) was submitted to a network of excitatory and inhibitory neurons, whose population spiking activity (amount of spikes per time unit, *Output*) is measured. The excess of activity in response to the input (blue area) measures the Responsiveness *R* of the network to the external stimulation. (b) Histograms of the resting potential $e_L = E_L/E_L$ of excitatoy (inhibitory) neurons in blue (red), top (low) row. Cells originate from the adult human brain, cortical layers 5/6 and 2/3, and mouse cortical layers 6 and 2/3 (Allen Brain Atlas³²). The continuous line is a Gaussian distribution with the same standard deviation as measured from the data. (c) Network activity in response to a time varying input $v_{ext}(t)$ of amplitude *A* = 1Hz and baseline value $v_0$ = 1.5Hz (see “Methods” section). Excitatory neurons population rate, normalised to pre-stimulus baseline activity, (top row, blue line) and the corresponding raster plot (bottom row), i.e. spiking times of excitatory (inhibitory) neurons marked with blue (red) dots. Panel (c) corresponds to an homogeneous network $\sigma_I = 0$, panel (d) to $\sigma_I = 0.1$ and panel (e) to $\sigma_I = 0.15$ (in these panels (c)–(e) $\sigma_E = 0$). (f) The evoked response *R* is reported in function of the heterogeneity of inhibitory (excitatory) neurons $\sigma = \sigma_I$ ($\sigma = \sigma_E$), red (blue) dots (squares). Error bars are estimated as the standard deviation over 20 different realisations. Continuous lines report the prediction based on the mean field model (see main text). In the inset, *R* is reported in function of the spontaneous excitatory firing rate $r_E$ without the stimulus (average over 10 s). Red dots are obtained by varying $\sigma_I$ (same data as the main panel). Green dots show *R* in function of $r_E$ for an homogeneous network ($\sigma_I = \sigma_E = 0$) and different values of the average resting potential of inhibitory neurons $E_L^I$.

**Optimal responsiveness for heterogeneous networks.** We studied networks of sparsely connected excitatory and inhibitory spiking neurons (see “Methods” for details). Each neuron receives afferent excitatory spike trains at a frequency $v_{ext}(t)$. To characterise the responsiveness of the network to external stimulation, we considered $v_{ext}(t)$ composed by a time constant baseline value $v_0$ and a time variation of amplitude *A* (see “Methods”, *Input* in Fig. 1a). The network responds to the stimulation by increasing its population spiking activity (see *Output* in Fig. 1a). The network responsiveness *R* can be estimated as the total amount of evoked spikes by the whole network (see the blue area in the *Output* of Fig. 1a). We compared homogeneous networks to networks constructed from measurements in the adult human and mouse brain from the Allen Brain Atlas³². In Fig. 1b, we report the histogram of the resting membrane potential $E_L$ measured experimentally from excitatory (blue, top panels) and inhibitory neurons (red, lower panels) in human and mouse cortical layers. By calling $\bar{E}_L$ the average value of $E_L$, one can see a heterogeneous distribution of the re-scaled resting potential $e_L = E_L/E_L^I$. The distribu-
<figure><img src="image_9.png" /><figcaption>**Figure 2.** Excitatory and inhibitory heterogeneity determines network responsiveness. (a) Responsiveness R as a function of heterogeneity of excitatory ($\sigma_E$) and inhibitory ($\sigma_I$) neurons for networks displaying relatively high levels of spontaneous activity (input baseline $v_0$ = 1.5 Hz and input amplitude $A$ = 1 Hz). (b) Responsiveness R for networks with a lower level of spontaneous activity (input baseline $v_0$ = 0.4 Hz and input amplitude $A$ = 0.6 Hz). Different markers correspond to heterogeneity values estimated from Human cortex layer 2, 3 (white circle), Human cortex layer 5, 6 (grey circle), mouse cortex layer 2/3 (grey star) and mouse cortex layer 6 (white diamond) estimated from Allen Brain database³², see Fig. 1.</figcaption></figure>

tions have a comparable level of variability, quantified by the standard deviations $\sigma_I$ and $\sigma_E$ (e.g. $\sigma_I \sim 0.07$ and $\sigma_E \sim 0.055$ for Human cortical layers 2/3). Note that heterogeneity in inhibitory cells is generally higher with respect to that in excitatory cells.

To investigate the effect of neuronal heterogeneity, we first studied the effect of excitatory neurons' heterogeneity $\sigma_E$ and inhibitory neurons' heterogeneity $\sigma_I$ separately. Panels c,d and e of Fig. 1 show that $\sigma_I$ strongly affects network responsiveness. In these simulations, the network displays spontaneous asynchronous activity where neurons fire irregularly, a dynamical regime due to the balance between excitation and inhibition typically observed in the cortex of awake animals, called asynchronous irregular¹⁷,³³,³⁴. The presence of external excitatory stimuli of short duration (see caption for details) produces a network response (an increase in population spiking activity) corresponding to a transient synchronous event, as it can be observed in the raster plots of Fig. 1c-e. Panel c corresponds to an homogeneous network, while panel d and e correspond to heterogeneous networks with increasing values of $\sigma_I$ ($\sigma_I$ = 0.1 in panel d and $\sigma_I$ = 0.15 in panel e, the resting potential $e_L$ follows a Gaussian distribution with fixed average, see “Methods” for details). By increasing $\sigma_I$ we observe a clear increase in the intensity of the burst of response (see panel d), that decreases for larger values of $\sigma_I$ (panel e in Fig 1). This can be quantified by estimating the responsiveness R as the amount of evoked spikes by the stimuli (see Fig. 1a and “Methods” section). By looking at Fig. 1e we observe that the same input induces a bell-shaped response R in function of the heterogeneity $\sigma_I$, indicating that there is an optimal heterogeneity level ($\sigma_I \sim 0.1$), where the stimuli provokes a strong population response. Eventually, when $\sigma_I$ is too large the response is very weak. On the other side, as it can be noticed from the raster plots in Fig. 1c-e, increasing heterogeneity in inhibitory neurons decreases excitatory neurons spontaneous activity. This is due to the presence of a fraction of inhibitory neurons with high values of the resting potential (closer to the firing threshold) that are more excitable and thus inhibit the excitatory population. In the inset of Fig. 1f we report R, as in the main panel, as a function of the excitatory neurons spontaneous activity pre-stimulus $r_E$. We can observe that the responsiveness is maximum, in correspondence of $\sigma_I \sim 0.1$, for a relatively low value of excitatory spontaneous activity (around $r_E \sim 0.3$ Hz). Interestingly, if we consider a homogeneous network ($\sigma_I = 0$) and we increase the average resting potential of inhibitory neurons, we observe a decrease of excitatory neurons spontaneous activity pre-stimulus $r_E$ but not an increase of responsiveness R (see green dots in the inset of Fig. 1f). This result shows that the increase in the size of the synchronous response is due to the presence of heterogeneity and cannot be replaced by a modification of the average excitability in the corresponding homogeneous network. If we now consider an homogeneous inhibitory population ($\sigma_I = 0$), we observe that increasing only the heterogeneity of excitatory neurons $\sigma_E$ has the effect to decrease network responsiveness (blue squares in Fig. 1f). To further investigate the role of heterogeneity in excitatory neurons we study, in the next section, the combined effect of heterogeneity in both populations.

**Responsiveness depends on heterogeneity in excitatory and inhibitory neurons.** To explore the combined effect of heterogeneity in both excitatory and inhibitory populations we calculated the responsiveness R as a function of both $\sigma_E$ and $\sigma_I$ (see Fig. 2a). Apart from fluctuations due to the endogenous noise in network dynamics, we observe a region of optimal responsiveness for heterogeneous networks (warm colors in Fig. 2).

Importantly, the experimental values of heterogeneity (Fig. 1b) fall close to the predicted optimal region of network responsiveness (see symbols in Fig. 2a). More specifically, we observe that, if heterogeneity in inhibitory neurons is very high (and thus responsiveness very low), increasing heterogeneity in excitatory neurons permits to increase responsiveness R. In another way, both heterogeneity cooperate to increase network responsiveness
<figure><img src="image_10.png" /><figcaption>**Figure 3.** Enhanced responsiveness corresponds to regions closer to instability. (a) The largest stability Lyapunov exponent $\lambda$ (real part) of asynchronous dynamics as a function of the heterogeneity of inhibitory neurons ($\sigma_I$). Different colors indicate different parameters of the baseline external drive and the strength of excitatory-excitatory quantal conductance ($v_0$, $Q_{EE}$), i.e. black (1.5 Hz, 1.5 nS), red (3 Hz, 1.5 nS), blue (2 Hz, 1.5 nS) and orange (3 Hz, 1.65 nS). Symbols are located at the value of $\sigma_I$ for which the responsiveness $R$ is maximum (same color code as continuous line). Different symbols indicate different amplitudes $A$ of the input, diamond ($A = 0.1$ Hz), star ($A = 0.5$ Hz) and dot ($A = 1$ Hz). (b) The larger stability Lyapunov exponent $\lambda$ of asynchronous dynamics as a function of the heterogeneity of inhibitory neurons ($\sigma_I$) and the strength of excitatory-excitatory quantal conductance $Q_{EE}$ for a baseline external drive $v_0 = 1.5$ Hz. The dotted (diamond) line is the responsiveness $R$ for an input amplitude $A = 0.1$ Hz ($A = 1$ Hz) and $Q_{EE} = 1.5$ nS (as in Figs. 1 and 2). Such responsiveness has been properly rescaled on the y-axes (i.e. multiplied by an ad-hoc factor) in order to fit in the image.</figcaption></figure>

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

$R$, that is low for homogeneous or too heterogeneous networks. Thus, the model predicts that the optimal heterogeneity matches that found in real neural networks, suggesting that this aspect is crucial to understand their responsiveness. In order to study the generality of this scenario, we consider a very different setup, with a lower background excitatory afferent input to the network (see caption in Fig. 2). In this parameter setup the homogeneous network has a lower spontaneous activity. In Fig. 2b we report the responsiveness $R$ for this setup. We observe indeed that also in this case a heterogeneous network yields a higher responsiveness $R$ and that the experimental amount of heterogeneity falls in the region of high responsiveness $R$. This result shows that, also across different parameter setups (other parameters are also studied in next sections), a heterogeneous network corresponds to an optimal responsiveness. Moreover, the optimal responsiveness appears for values of heterogeneity close to experimental amounts. In order to test the robustness of this observation, we have verified that changing network realization and initial conditions the results are the same (see Supplementary Information).

**Optimal responsiveness comes from pushing the network at the edge of a dynamical transition.** To determine the mechanism at the origin of the enhanced responsiveness due to neuronal heterogeneity, we developed a mean-field approach explicitly taking into account diversity. We started from a mean-field model previously introduced for homogeneous neural populations³⁵⁻³⁷. We extend this approach to heterogeneous systems by employing a technique, called heterogeneous mean field (HMF), successfully applied previously to model networks with heterogeneous connectivity³⁸,³⁹ and extended here to networks with heterogeneous cell properties (see “Methods” section). This procedure describes the collective dynamics of large sparsely connected networks of heterogeneous neurons with a relatively simple three dimensional model. This HMF model predicts the time evolution of excitatory neurons population firing rate $r_E(t)$ and of inhibitory neurons population firing rate $r_I(t)$. By comparing the prediction of the HMF model to estimations from direct numerical simulations of a large network of neurons, we observed a very good agreement, that can be appreciated from the prediction of the response $R$ in function of $\sigma_I(\sigma_E)$ (see continuous lines in Fig. 1d).

The excellent match of the continuous curves in Fig. 1 shows that the HMF correctly predicts the numerical observations of an optimal amount of heterogeneity in inhibitory neurons for population response. We can go further, and exploit the simplicity of the HMF model to understand what are the differences between more or less responsive asynchronous states. By computing the stationary solution from the HMF we estimate the relative stability eigenvalues $\{\lambda_i\}$. The maximum exponent $\lambda$ is reported in Fig. 3a for different parameter setups. The closer $\lambda$ to the critical value $\lambda_c = 0$ the less the asynchronous state is stable to perturbations. We observe that, independently of model parameters it exists an intermediate value of heterogeneity $\sigma_I$ for which $\lambda$ is maximum and closer to the transition point $\lambda = 0$. We have then estimated the responsiveness $R$ for these setups and we observe that the maximum responsiveness appears at values of the heterogeneity $\sigma_I$ at which the asynchronous state is less stable (see symbols in Fig. 3a). In Fig. 3b then we report $\lambda$ in function of $\sigma_I$ and the intensity of excitatory-excitatory neurons synaptic strength $Q_{EE}$. We observe that, increasing $Q_{EE}$ there exists a
<figure><img src="image_12.png" /><figcaption>**Figure 4.** Heterogeneous networks admit more dynamical regimes compared to homogeneous networks. (a) Inhibitory neurons population stationary firing rate $r_I^*$ in function of $\sigma_I$ for $E_L^I = -70$ mV. Black line indicates stable ($\lambda < 0$) asynchronous state, dashed blue line indicates unstable ($\lambda > 0$) asynchronous state where a limit cycle appears (red lines indicates maximum and minimum value of $r_I$ in time). (b) The largest stability Lyapunov exponent $\lambda^-$ (real part) of the asynchronous state as a function of the average resting potential of inhibitory neurons $E_L^I$ and their standard deviation $\sigma_I$. Whenever $\lambda > 0$ the asynchronous state is unstable and a stable limit cycle appears. In direct network simulations we observe sparsely synchronous oscillations (see raster plots in panel (c) where we use $E_L^I = -70$ mV and $\sigma_I = 0, 0.12, 0.2$ from top to bottom). In these simulations $Q_{EE} = 1.53$ nS.</figcaption></figure>

region where the asynchronous state becomes marginally stable (red region, $\lambda$ close to zero). By superimposing the responsiveness $R$ from Fig. 1d (dotted line, in this case for fixed $Q_{EE} = 1.5$ nS), we can see that the value of $\sigma_I$ for which the responsiveness is maximum corresponds to a network that is closer to the transition point where the asynchronous regime becomes unstable.

**Heterogeneity can induce new dynamical regimes.** Finally, to determine what the unstable regions correspond to, we investigated heterogeneous networks both numerically and theoretically using the mean-field model. Figure 4a shows that, for intermediate amount of heterogeneity the asynchronous state can be unstable (blue dashed line in panel a). This dynamical phase corresponds to an oscillatory regime at the level of the whole network, arising from a super-critical Hopf bifurcation where the amplitude of the limit cycle (red lines in panel a) smoothly increases in function $\sigma_I$. In panel b of Fig. 4 we observe that this oscillating synchronous regime (present whenever the asynchronous state is unstable, $\lambda > 0$) appears only for heterogeneous networks ($\sigma_I>0$) but the location of the synchronous region depends on the average value of neurons' resting potential $E_L^I$ (notice that in panel a of Fig. 4 we used $E_L^I = -70$ mV and in Figs. 1, 2 and 3 we used $E_L^I = -65$ mV). In order to verify these theoretical predictions, in panel c we finally perform numerical simulation of the neuronal network where we observe that the limit cycle corresponds indeed to collective network oscillations (middle panel) that disappear in homogeneous or too heterogeneous setups (higher and lower panels). This shows that not only heterogeneity can induce higher responsiveness, but it can also induce new dynamical regimes. The dynamics of heterogeneous networks is therefore richer.

**Heterogeneity boosts information flow in large-scale networks.** To determine if the enhanced responsiveness applies not only to external inputs but also to the flow of information through recurrent connectivity, we considered large-scale networks where each unit is a population of heterogeneous neurons. We considered a two dimensional lattice of mean-field models, that are interconnected to each other via Gaussian connectivity profiles (see Fig. 5a and “Methods” section). According to anatomical connectivity estimates, inhibitory neurons have a short-range connectivity (red curves in Fig. 5a) at variance with excitatory connectivity (blue curves in Fig. 5a)$^{40}$. Importantly, we integrate distance-dependent propagation delays due to the finite velocity of axonal conduction of action potentials (we considered here an axonal conduction velocity of 0.3 m/s, see “Methods” section).

We stimulated this large-scale model by an external input (see “Methods” section) and we compared a model with locally heterogeneous neurons (each node is modeled with an heterogeneous mean field with $\sigma_I = 0.15$) versus the same model with locally homogeneous neurons (each node is modeled with an homogeneous mean field with $\sigma_I = 0$). As it can be appreciated from Fig. 5b, both the heterogeneous and the homogeneous model respond with an activation that propagates in space as a wave. We observed that the intensity of the wave is much larger in the model with locally heterogeneous networks. This is a direct consequence of the enhanced responsiveness of each node. But does local heterogeneity also impacts the spatial extent of the wave through recurrent connections? To answer this question we compared the homogeneous and the heterogeneous system
<figure><img src="image_13.png" /><figcaption>**Figure 5.** Enhanced activity propagation in large-scale heterogeneous networks. (a) Two dimensional lattice of connected mean-field models. The connectivity between excitatory (inhibitory) cells is drawn from a Gaussian distribution with standard deviation of length $l_{exc} = 2$ mm ($l_{inh} = 1$ mm), see blue (red) curves. (b) Response of the system (firing rate of excitatory neurons) following a Gaussian afferent stimulus (see "Methods"). The upper row stands for locally homogeneous networks ($\sigma_I = 0$) while lower row for locally heterogeneous networks ($\sigma_I = 0.15$). (c) Spatio-temporal profile of excitatory neurons firing rate normalised by its maximum in space and time ($\sigma_I = 0$). (d) Same as (c) but for $\sigma_I = 0.15$. (e) Normalised firing rate in function of the distance $d$ from stimulus onset at a specific time [see dashed white line in panel (d)]. Different colors stand for different levels of heterogeneity (see the legend inside the panel). The dashed red line shows a fit of the exponential decay of activity found in the heterogeneous system (red line, $\sigma_I = 0.15$), with $\lambda_d \sim 3.5$ mm.</figcaption></figure>

by normalising the firing rate activity with its maximum in space and time, in order to highlight the spatio-temporal profile of the response. As it can be observed from panels c and d in Fig. 5, the homogeneous and the heterogeneous systems have different spatio-temporal profiles, indicating that local heterogeneity influences not only the intensity of the response to the stimulus, but also its spatio-temporal pattern through recurrent interactions. Remarkably the heterogeneous system has much longer-range spatial propagation of activity. More specifically, for an optimal amount of heterogeneity we observed a much enhanced propagation of the activity through the network (red curve in Fig. 5e), compared with homogeneous and highly heterogeneous networks.

These results show that an optimal amount of local heterogeneity has a double effect at large scale. First, it amplifies the intensity of the response to the external input (Fig. 5b, a consequence of local enhanced responsiveness) and, second, enhanced responsiveness also applies to recurrent inputs, which results in a long range propagation of activity across the network (Fig. 5e).

## Discussion

In conclusion, we report here four findings. First, we have found that the heterogeneity of inhibitory neurons, which has been well documented experimentally$^{13-16}$, optimises the responsiveness of spontaneously active networks to external stimuli. There appears a resonance peak as a function of the level of heterogeneity. We have here studied the heterogeneity in neurons' resting potential, that is experimentally quantified (see Fig. 1b and Supplementary Information) and determines neurons' proximity to firing threshold (and thus neurons' excitability). On the other side our results do not limit to only this cellular parameter. We have indeed found, both theoretically and by performing network simulations, that an optimal responsiveness appears also for intermediate heterogeneity in neurons' leakage conductance, while almost no effect is observed for heterogeneity in neurons' membrane capacitance (see Supplementary Information).

A similar effect of diversity-induced resonance was previously observed in excitable or bistable systems$^{41-45}$, where heterogeneity creates active excitation clusters which were absent in the quiescent homogeneous system. We showed here a different type of heterogeneity, that of inhibitory neurons, can induce optimal population responses. More generally, we found that the heterogeneity of both excitatory and inhibitory populations induces optimal responsiveness in spontaneously active sparse networks with irregular firing activity. Importantly, we found that the level of heterogeneity measured experimentally across different cortical layers and species (Human
and Mouse) corresponds to the resonance peak, which suggests that cortical networks may have naturally evolved towards optimal responsiveness by adjusting their heterogeneity. Moreover, while several studies reported that heterogeneity can enhance coding in uncoupled networks$^{46,47}$ and decrease neuronal correlations$^{48–50}$, we report here that a higher input-output population response is linked to an increased tendency to synchronization in heterogeneous networks. The coding capabilities of neural networks will therefore be largely affected by neuronal heterogeneity, which opens interesting perspectives for future studies.

Second, we found that the enhanced responsiveness of heterogeneous networks is paralleled with a decreased stability of the spontaneous activity regime. To obtain this result, we designed a mean-field model that explicitly includes heterogeneity, and which can capture this diversity-induced resonance. This new mean-field formulation keeps track of microscopic complexity, compared to traditional mean-field approaches which implicitly assume homogeneous systems and would not predict the correct responsiveness. This also shows that responsiveness must be understood using the knowledge of the spontaneous activity of the network—and in particular its level of stability. The relation between instability and responsiveness also constitutes a promising subject to further explore in the future.

In this work we have studied heterogeneity in the intrinsic proprieties of neurons, but it is known that heterogeneity can be found also in the structure of connections (say for example in neurons’ in-degrees). Structural Heterogeneity (i.e. in the topology of connections) can have a strong impact in network dynamics$^{11,51–53}$. Our mean field approach is general and can be employed in the future to study the impact of heterogeneity in neurons’ connectivity.

Third, we have shown that neuronal heterogeneity is not only important for responsiveness, but also can induce new dynamical regimes. Using the mean-field models, we could predict a transition to a sparsely synchronous collective oscillation regime, which was confirmed by network simulations. This type of diversity-induced oscillations reminds some aspects found in noise-induced transitions in dynamical systems$^{54,55}$. Whether the effects of heterogeneity could be considered as analogous to the effect of noise (“quenched noise”) in neural networks is also an interesting direction for future studies.

Finally, we have found that the enhanced responsiveness enables activity to propagate much easier in large-scale heterogeneous networks. The enhanced response not only applies to the stimulus, but also it applies to information flow through recurrent excitatory inputs. As a result, comparing heterogeneous to homogeneous networks, a given stimulus produces a larger local response, and in addition, this response also propagates to a larger spatial extent because the effect of recurrent excitatory inputs is also amplified by heterogeneity. Thus, we conclude that heterogeneous networks can provide activity propagation at a level much superior compared to a homogeneous system, and thus will necessarily better propagate information and make it available to larger brain areas.

## Methods

**Network model.** We examined networks of excitatory and inhibitory neurons connected through conductance based synapses. We used networks of $N = 10,000$ neurons, 80% of excitatory ($N_E = 0.8 N$) and 20% of inhibitory ($N_I = 0.27 N$) neurons. The membrane potential $V_i$ of each neuron evolves according to the Adaptive Exponential integrate and fire model (AdExp)$^{56}$:

$$C_m \dot{V}_i = g_L(E_L^i - V_i) + g_L e^{\frac{V_i - v_{th}}{\Delta}} + I_s^i - w_i \quad (1)$$

$$\tau_w \dot{w}_i = -w_i + b \sum_{\{t_j^{sp}\}} \delta(t - t_i^{sp}), \quad (2)$$

where $C_m = 200$ pF is the membrane capacitance, $g_L = 15$ nS the leakage conductance, $v_{th} = 50$ mV the effective threshold and $\Delta$ defines the action potential rise ($\Delta = 0.5$ mV for inhibitory neurons and $\Delta = 2$ mV for excitatory neurons). The adaptation current $w_i$ increases of an amount $b = 60$ nS at each spike emitted by neuron $i$ at times $\{t_i^{sp}\}$ and has an exponential decay with time scale $\tau_w = 500$ ms. Only excitatory neurons have spike frequency adaptation, while for inhibitory neurons $b = 0$. The current $I_s^i$ is the current received by neuron $i$ from other neurons in the network. We consider a random graph where each couple of neurons is connected with probability $p = 0.05$. By calling $\{t_j^{sp}\}$ the ensemble of spiking times of neuron $j$ we have, for an excitatory post-synaptic neuron $i$:

$$I_s^i = g_i^{EE}(V_i - E_E) + g_i^{EI}(V_i - E_I) \quad (3)$$

$$\tau_s \dot{g}_i^{EE} = -g_i^{EE} + Q_{EE} \sum_{\{t_j^{sp}\} \in (E)} \delta(t - t_j^{sp}) \quad (4)$$

$$\tau_s \dot{g}_i^{EI} = -g_i^{EI} + Q_{EI} \sum_{\{t_j^{sp}\} \in (I)} \delta(t - t_j^{sp}), \quad (5)$$

where $E_{E,I}$ is the reversal for excitatory ($E_E = 0$ mV) and inhibitory synapses ($E_I = -80$ mV), $\tau_s = 5$ ms the synaptic decay time and $Q_{EE}$ ($Q_{EI}$) is the interaction strength of excitatory (inhibitory) synapses to excitatory neurons. The same equations (with $g_i^{EE} \rightarrow g_i^{IE}$, $g_i^{EI} \rightarrow g_i^{II}$, $Q_{EE} \rightarrow Q_{IE}$ and $Q_{EI} \rightarrow Q_{II}$) apply for inhibitory
post-synaptic neurons. We fix $Q_{EI} = Q_{II} = 5$ nS and $Q_{EE} = Q_{IE} = 1.5$ nS (we employed different values of $Q_{EE}$ in Fig. 3 and in Fig. 4, see the relative caption).

Each neuron receives an external Poissonian train of excitatory spikes at a rate $v_{ext}$. The value of $v_{ext}$ determines the amount of ongoing spontaneous activity in the network. In order to study the response to external stimuli we considered a time varying $v_{ext}(t)$ of the form $v_{ext}(t) = v_0 + Ae^{-\frac{(t-t_0)^2}{2\tau^2}}$, where A is the input amplitude, $t_0$ the time when input is maximum and $T = 50$ ms measures the duration of the input. We have employed $v_0 = 1.5$ Hz in all the numerical simulations, apart from Figs. 2b and 3 (see Figure caption for details). In Fig. 1 we employed $t_0 = 6$ s, $v_0 = 1.5$ Hz and A = 1 Hz. The responsiveness R is estimated by computing the amount of spikes of excitatory neurons while the input is on (i.e. between $t = t_0 - 3T$ and $t = t_0 - 3T$) minus the baseline activity (average of excitatory spike rate for A = 0). Responsiveness R is estimated by averaging over 20 different repetitions of this procedure.

To model heterogeneity, we considered a Gaussian distribution of the resting potential $E_L^I$ of inhibitory (excitatory) population $N(E_L^{E,I}, \sigma_{E,I}^2)$ with average $\overline{E_L^{E,I}}$ and standard deviation $\tilde{\sigma}_{E,I}$. We considered $\overline{E_L^E} = \overline{E_L^I} = -65$ mV if not stated otherwise (e.g. we used different values of $E_L^I$ in the inset of Fig. 1d and in Fig. 4). The re-scaled standard deviation $\sigma_{E,I} = \tilde{\sigma}_{E,I}/\overline{E_L^{E,I}}$ is the main parameter to quantify heterogeneity. The same definitions apply for heterogeneity in other parameters (see Supplementary Information for $g_L$ and $C_m$).

**Mean field model.** In the homogeneous case a mean field model for this network has been recently developed³⁷. By employing a Markovian approximation over a time scale $\tau = 15$ ms and by considering a sufficiently slow time scale for the dynamics of adaptation $\tau_w$, mean field equations read:

$$\tau \dot{r}_I = F^I(r_E + v_{ext}, r_I) - r_I \qquad (6)$$

$$\tau \dot{r}_E = F^E(r_E + v_{ext}, r_I, W) - r_E \qquad (7)$$

$$\tau_w \dot{W} = -W + br_E, \qquad (8)$$

where $r_E(r_I)$ is excitatory (inhibitory) neurons population firing rate, W is excitatory neurons average spike frequency adaptation, $F^I(v_E, v_I)$ and $F^E(v_E, v_I, W)$ are the transfer functions of inhibitory and excitatory neurons, respectively. They measure the stationary firing rate of one neuron when receiving an excitatory (inhibitory) Poissonian spike train at a rate $v_E(v_I)$.

In the case of heterogeneous inhibitory neurons, a parameter x, say neurons reversal potential $E_L$, is distributed according to a probability density function $P(x)$. In the limit of large networks we can decompose inhibitory neurons in classes, each one characterized by a parameter x and by its own transfer function $F_x^I(r_E + v_{ext}, r_I)$. We indicate with $r_x^I$ the firing rate of the class of neurons with parameter x. The whole population rate is $r_I = \int dxP(x)r_x^I$. The equations are closed by a self consistency equation for the mean input received by one neuron and we need to replace Eq. (6) with

$$\tau \dot{r}_I = \int dxP(x)F_x^I(r_E + v_{ext}, r_I) - r_I. \qquad (9)$$

Notice that the model still stays three dimensional but keeps track of the distribution $P(x)$ of heterogeneity. In the case $P(x) = \delta(x - x_0)$ we recover the homogeneous model.

The same procedure applies for excitatory neurons. Nevertheless, in this case each class x is characterised by its own adaptation variable $w_x$. For each class with parameter x we get:

$$\tau \dot{r}_x^E = F_x^E(r_E + v_{ext}, r_I, w_x) - r_E \qquad (10)$$

$$\tau_w \dot{w}_x = -w_x + br_x^E. \qquad (11)$$

In this case the population quantities can be written as $W = \int P(x)w_x$ and $r_E = \int P(x)r_x^E$. These equations can be solved by sampling $P(x)$ and we found that a sampling of around 50 points gives an accurate precision. In our work we employ the mean field for heterogeneous excitatory neurons only in Fig. 1d, for which we employed a sampling of 50 points.

Nevertheless, it is possible to reduce the dimensionality of this model by making the hypothesis that $w_x$ is slow enough. The stationary solution is $w_x = br_x^E$. In order to follow the dynamics of adaptation we evolve the population variable W as $\tau_w \dot{W} = -W + br^E$ and estimate the adaptation of each class from the equation $w_x = W r_x^E / r_E = (W/r_E)F_x^E(r_E + v_{ext}, r_I, w_x)$. In this way the model stays three dimensional. Preliminary results indicate that this is a good approximation for the population dynamics in the heterogeneous case.

Finally, the estimation of neuron transfer functions for the AdExp model $F^{E/I}(r_E, r_I, W)$ is done through a semi-analytical fitting procedure³⁷ (see Supplementary Information).

**Spatially extended model.** To model large-scale networks, we considered a square lattice of length $L = 20$ mm composed of $M \times M$ nodes, we employed $M = 100$ in the simulations of Fig. 5. Each node is modeled as a mean field model. In the limit of large M the dynamics of each node at a location (x, y) follows:
$$ \tau \frac{\partial r_E(x, y)}{\partial t} = -r_E(x, y) + F^E(v_E(x, y) + v_{\text{inp}}(x, y, t), v_I(x, y), W(x, y)) \quad (12) $$

$$ \tau \frac{\partial r_I(x, y)}{\partial t} = -r_I(x, y) + \int dz P(z) F_z^I(v_E(x, y) + v_{\text{inp}}(x, y, t), v_I(x, y)) \quad (13) $$

$$ \tau_w \dot{W}(x, y) = -W(x, y) + br_E(x, y), \quad (14) $$

where $P(z)$ is the distribution of heterogeneity, i.e. a Gaussian distribution with average $\overline{E}_I^I$ and rescaled standard deviation $\sigma_I$, and $v_E(x, y)$ ($v_I(x, y)$) is the excitatory (inhibitory) input incoming in $(x, y)$ from the other lattice locations, i.e.

$$ v_E(x, y) = \int dx_1 \int dy_1 G_E(x - x_1, y - y_1) r_E \left( x_1, y_1, t - \frac{d}{v_c} \right) \quad (15) $$

$$ v_I(x, y) = \int dx_1 \int dy_1 G_I(x - x_1, y - y_1) r_I \left( x_1, y_1, t - \frac{d}{v_c} \right), \quad (16) $$

where $d = \sqrt{(x - x_1)^2 + (y - y_1)^2}$ is the distance between two points in the lattice and $G_E(G_I)$ is the excitatory (inhibitory) connectivity in space. We consider a Gaussian connectivity in both directions with standard deviation $\Delta_E = 2$ mm ($\Delta_I = 1$ mm) for excitatory (inhibitory) connections. We considered an axonal conduction velocity $v_c = 0.3$ m/s. We considered an external stimulation as a stationary input $v_{\text{inp}}(x, y, t)$:

$$ v_{\text{inp}}(x, y, t) = \frac{A}{2\pi \Delta_{\text{inp}}^2} e^{-\frac{(t-t_0)^2}{2T^2}} e^{-\frac{(x-x_0)^2}{2\Delta_{\text{inp}}^2}} e^{-\frac{(y-y_0)^2}{2\Delta_{\text{inp}}^2}}, \quad (17) $$

with $T = 50$ ms, $\Delta_{\text{inp}} = 1.5$ mm, $A = 1$ Hz, $t_0 = 100$ ms and $x_0 = y_0 = 10$ mm.

Received: 30 April 2021; Accepted: 11 August 2021

Published online: 02 September 2021

## References

1. Tsodyks, M., Mitkov, I. & Sompolinsky, H. Pattern of synchrony in inhomogeneous networks of oscillators with pulse interactions. *Phys. Rev. Lett.* **71**, 1280 (1993).

2. Roxin, A., Brunel, N., Hansel, D., Mongillo, G. & van Vreeswijk, C. On the distribution of firing rates in networks of cortical neurons. *J. Neurosci.* **31**, 16217–16226 (2011).

3. Roxin, A. The role of degree distribution in shaping the dynamics in networks of sparsely connected spiking neurons. *Front. Comput. Neurosci.* **5**, 8 (2011).

4. Denker, M., Timme, M., Diesmann, M., Wolf, F. & Geisel, T. Breaking synchrony by heterogeneity in complex networks. *Phys. Rev. Lett.* **92**, 074103 (2004).

5. Neltner, L., Hansel, D., Mato, G. & Meunier, C. Synchrony in heterogeneous networks of spiking neurons. *Neural Comput.* **12**, 1607–1641 (2000).

6. Golomb, D. & Rinzel, J. Dynamics of globally coupled inhibitory neurons with heterogeneity. *Phys. Rev. E* **48**, 4810 (1993).

7. Zerlaut, Y. et al. Heterogeneous firing rate response of mouse layer v pyramidal neurons in the fluctuation-driven regime. *J. Physiol.* **594**, 3791–3808 (2016).

8. Tseng, G.-F. & Prince, D. A. Heterogeneity of rat corticospinal neurons. *J. Comp. Neurol.* **335**, 92–108 (1993).

9. Pospischil, M. et al. Minimal Hodgkin–Huxley type models for different classes of cortical and thalamic neurons. *Biol. Cybern.* **99**, 427–441 (2008).

10. Sharpee, T. O. Toward functional classification of neuronal types. *Neuron* **83**, 1329–1334 (2014).

11. Landau, I. D., Egger, R., Dercksen, V. J., Oberlaender, M. & Sompolinsky, H. The impact of structural heterogeneity on excitation-inhibition balance in cortical networks. *Neuron* **92**, 1106–1121 (2016).

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

13. Gupta, A., Wang, Y. & Markram, H. Organizing principles for a diversity of GABAergic interneurons and synapses in the neocortex. *Science* **287**, 273–278 (2000).

14. DeFelipe, J. et al. New insights into the classification and nomenclature of cortical GABAergic interneurons. *Nat. Rev. Neurosci.* **14**, 202–216 (2013).

15. Jiang, X. et al. Principles of connectivity among morphologically defined cell types in adult neocortex. *Science* **350**, aac9462 (2015).

16. Mihaljevic, B., Benavides-Piccione, R., Bielza, C., Larraaga, P. & DeFelipe, J. Classification of GABAergic interneurons by leading neuroscientists. *Sci. Data* **6**, 221 (2019).

17. Van Vreeswijk, C. & Sompolinsky, H. Chaos in neuronal networks with balanced excitatory and inhibitory activity. *Science* **274**, 1724–1726 (1996).

18. Vardi, R. et al. Sudden synchrony leaps accompanied by frequency multiplications in neuronal activity. *Front. Neural Circuits* **7**, 176 (2013).

19. Vardi, R. et al. Neuronal response impedance mechanism implementing cooperative networks with low firing rates and µs precision. *Front. Neural Circuits* **9**, 29 (2015).

20. Goldental, A., Uzan, H., Sardi, S. & Kanter, I. Oscillations in networks of networks stem from adaptive nodes with memory. *Sci. Rep.* **7**, 1–8 (2017).

21. Chow, C. C. Phase-locking in weakly heterogeneous neuronal networks. *Phys. D* **118**, 343–370 (1998).

22. White, J. A., Chow, C. C., Rit, J., Soto-Treviño, C. & Kopell, N. Synchronization and oscillatory dynamics in heterogeneous, mutually inhibited neurons. *J. Comput. Neurosci.* **5**, 5–16 (1998).
23. Traub, R. D. et al. Single-column thalamocortical network model exhibiting gamma oscillations, sleep spindles, and epileptogenic bursts. *J. Neurophysiol.* **93**, 2194–2232 (2005).

24. Wang, X.-J. & Buzsáki, G. Gamma oscillation by synaptic inhibition in a hippocampal interneuronal network model. *J. Neurosci.* **16**, 6402–6413 (1996).

25. Tiesinga, P. & José, J. V. Robust gamma oscillations in networks of inhibitory hippocampal interneurons. *Netw. Comput. Neural Syst.* **11**, 1–23 (2000).

26. Devalle, F., Montbrió, E. & Pazó, D. Dynamics of a large system of spiking neurons with synaptic delay. *Phys. Rev. E* **98**, 042214 (2018).

27. Luccioli, S., Angulo-Garcia, D. & Torcini, A. Neural activity of heterogeneous inhibitory spiking networks with delay. *Phys. Rev. E* **99**, 052412 (2019).

28. Ostojic, S., Brunel, N. & Hakim, V. Synchronization properties of networks of electrically coupled neurons in the presence of noise and heterogeneities. *J. Comput. Neurosci.* **26**, 369 (2009).

29. Palmigiano, A., Geisel, T., Wolf, F. & Battaglia, D. Flexible information routing by transient synchrony. *Nat. Neurosci.* **20**, 1014 (2017).

30. Gray, C. M. Synchronous oscillations in neuronal systems: Mechanisms and functions. *J. Comput. Neurosci.* **1**, 11–38 (1994).

31. Llinas, R. & Ribary, U. Coherent 40-hz oscillation characterizes dream state in humans. *Proc. Natl. Acad. Sci.* **90**, 2078–2081 (1993).

32. Allen Brain Institute. Allen cell types database. https://celltypes.brain-map.org/data?donor__species=Homo%20Sapiens (Accessed Dec 2020) (2015).

33. Monteforte, M. & Wolf, F. Dynamical entropy production in spiking neuron networks in the balanced state. *Phys. Rev. Lett.* **105**, 268104 (2010).

34. Destexhe, A., Rudolph, M. & Paré, D. The high-conductance state of neocortical neurons in vivo. *Nat. Rev. Neurosci.* **4**, 739–751 (2003).

35. El Boustani, S. & Destexhe, A. A master equation formalism for macroscopic modeling of asynchronous irregular activity states. *Neural Comput.* **21**, 46–100 (2009).

36. Zerlaut, Y., Chemla, S., Chavane, F. & Destexhe, A. Modeling mesoscopic cortical dynamics using a mean-field model of conductance-based networks of adaptive exponential integrate-and-fire neurons. *J. Comput. Neurosci.* **44**, 45–61 (2018).

37. di Volo, M., Romagnoni, A., Capone, C. & Destexhe, A. Biologically realistic mean-field models of conductance-based networks of spiking neurons with adaptation. *Neural Comput.* **31**, 653–680 (2019).

38. di Volo, M., Burioni, R., Casartelli, M., Livi, R. & Vezzani, A. Heterogeneous mean field for neural networks with short-term plasticity. *Phys. Rev. E* **90**, 022811 (2014).

39. Burioni, R., Casartelli, M., Di Volo, M., Livi, R. & Vezzani, A. Average synaptic activity and neural networks topology: A global inverse problem. *Sci. Rep.* **4**, 1–7 (2014).

40. Buzás, P. et al. Model-based analysis of excitatory lateral connections in the visual cortex. *J. Comp. Neurol.* **499**, 861–881 (2006).

41. Tessone, C. J., Mirasso, C. R., Toral, R. & Gunton, J. D. Diversity-induced resonance. *Phys. Rev. Lett.* **97**, 194101 (2006).

42. Assisi, C. G., Jursa, V. K. & Kelso, J. S. Synchrony and clustering in heterogeneous networks with global coupling and parameter dispersion. *Phys. Rev. Lett.* **94**, 018106 (2005).

43. Lafuerza, L. F., Colet, P. & Toral, R. Nonuniversal results induced by diversity distribution in coupled excitable systems. *Phys. Rev. Lett.* **105**, 084101 (2010).

44. Mejias, J. & Longtin, A. Optimal heterogeneity for coding in spiking neural networks. *Phys. Rev. Lett.* **108**, 228102 (2012).

45. Mejias, J. F. & Longtin, A. Differential effects of excitatory and inhibitory heterogeneity on the gain and asynchronous state of sparse cortical networks. *Front. Comput. Neurosci.* **8**, 107 (2014).

46. Tripathy, S. J., Padmanabhan, K., Gerkin, R. C. & Urban, N. N. Intermediate intrinsic diversity enhances neural population coding. *Proc. Natl. Acad. Sci.* **110**, 8248–8253 (2013).

47. Beiran, M., Kruscha, A., Benda, J. & Lindner, B. Coding of time-dependent stimuli in homogeneous and heterogeneous neural populations. *J. Comput. Neurosci.* **44**, 189–202 (2018).

48. Padmanabhan, K. & Urban, N. N. Intrinsic biophysical diversity decorrelates neuronal firing while increasing information content. *Nat. Neurosci.* **13**, 1276 (2010).

49. Shamir, M. & Sompolinsky, H. Implications of neuronal diversity on population coding. *Neural Comput.* **18**, 1951–1986 (2006).

50. Pfeil, T. et al. Effect of heterogeneity on decorrelation mechanisms in spiking neural networks: A neuromorphic-hardware study. *Phys. Rev. X* **6**, 021023 (2016).

51. Brama, H., Guberman, S., Abeles, M., Stern, E. & Kanter, I. Synchronization among neuronal pools without common inputs: In vivo study. *Brain Struct. Funct.* **220**, 3721–3731 (2015).

52. Kanter, I. et al. Nonlocal mechanism for cluster synchronization in neural circuits. *EPL (Europhys. Lett.)* **93**, 66001 (2011).

53. di Volo, M., Burioni, R., Casartelli, M., Livi, R. & Vezzani, A. Neural networks with excitatory and inhibitory components: Direct and inverse problems by a mean-field approach. *Phys. Rev. E* **93**, 012305 (2016).

54. Horsthemke, W. & Lefever, R. *Noise-Induced Transitions* 2nd edn. (Springer, 2006).

55. Lindner, B., Garcia-Ojalvo, J., Neiman, A. & Schimansky-Geier, L. Effects of noise in excitable systems. *Phys. Rep.* **392**, 321–424 (2004).

56. Brette, R. & Gerstner, W. Adaptive exponential integrate-and-fire model as an effective description of neuronal activity. *J. Neurophysiol.* **94**, 3637–3642 (2005).

## Acknowledgements

Research funded by the CNRS, the European Community (Human Brain Project H2020-785907 and H2020-945539) and the French National Research Agency (ANR Projects ERMUNDY, PARADOX and SHOOTINGSTAR).

## Author contributions

M.V. performed numerical simulations, analysed data and prepared figures. M.V. and A.D. wrote the manuscript and conceived the work.

## Competing interests

The authors declare no competing interests.

## Additional information

**Supplementary Information** The online version contains supplementary material available at https://doi.org/10.1038/s41598-021-96745-2.

Correspondence and requests for materials should be addressed to M.V.
**Reprints and permissions information** is available at www.nature.com/reprints.

**Publisher's note** Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.

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

**Open Access** This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons.org/licenses/by/4.0/.

© The Author(s) 2021