# Myelin-induced gain control in nonlinear neural networks

Jérémie Lefebvre, Andrew Clappison, André Longtin, Axel Hutt

► To cite this version:

Jérémie Lefebvre, Andrew Clappison, André Longtin, Axel Hutt. Myelin-induced gain control in nonlinear neural networks. Communications Physics, 2025, 8 (1), pp.145. 10.1038/s42005-025-02055-8. hal-04835140

HAL Id: hal-04835140

https://inria.hal.science/hal-04835140v1

Submitted on 13 Dec 2024

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
# Myelin-induced gain control in nonlinear neural networks

Jérémie Lefebvre<sup>1,2,3,4*†</sup>, Andrew Clappison<sup>5</sup>, André Longtin<sup>2†</sup>, Axel Hutt<sup>6,7†</sup>

<sup>1*</sup>Department of Biology, University of Ottawa, Ottawa, ON K1N 6N5, Ontario, Canada.

<sup>2</sup>Department of Physics, University of Ottawa, Ottawa, ON K1N 6N5, Ontario, Canada.

<sup>3</sup>Krembil Brain Institute, University Health Network, Toronto, ON M5T 0S8, Ontario, Canada.

<sup>4</sup>Department of Mathematics, University of Toronto, Toronto, ON M5S 2E4, Ontario, Canada.

<sup>5</sup>School of Electrical Engineering and Computer Science, University of Ottawa, Ottawa, ON K1N 6N5, Ontario, Canada.

<sup>6</sup>Mimesis, Inria Centre at Université de Lorraine, Strasbourg, 67000, France.

<sup>7</sup>MLMS, Université de Strasbourg, CNRS, iCube, Strasbourg, 67000, France.

*Corresponding author(s). E-mail(s): jeremie.lefebvre@uottawa.ca;

Contributing authors: aclap051@uottawa.ca; alongtin@uottawa.ca;

axel.hutt@inria.fr;

†These authors contributed equally to this work.

## Abstract

Myelin surrounds axonal membranes to increase the conduction velocity of nerve impulses and thus reduce communication delays in neural signaling. Changes in myelination alter the distribution of delays in neural circuits, but the implications for their operation are poorly understood. We present a joint computational and non-linear dynamical method to explain how myelin-induced changes in axonal conduction velocity impact the firing rate statistics and spectral response properties of recurrent neural networks. Using a network of spiking neurons with distributed conduction delays driven by a spatially homogeneous noise, combined
probabilistic and mean field approaches reveal that myelin implements a gain control mechanism, optimizing neural coding while stabilizing neural dynamics away from oscillatory regimes. The effect of myelin-induced changes in conduction velocity on network dynamics was found to be more pronounced in presence of correlated stochastic stimuli. Further, computational and theoretical power spectral analyses reveals a paradoxical effect where the loss of myelin promotes oscillatory responses to broadband time-varying stimuli. Taken together, our findings show that myelination can play a fundamental role in neural computation and its impairment in myelin pathologies such as epilepsy and multiple sclerosis.

**Keywords:** Myelin, Delays, Recurrent networks, Computation, Nonlinear dynamics, Noise, Mean-field

Nonlinear systems with delayed feedback have widespread applications across the fields of physics[1], engineering [2] and biology [3]. In biological neural networks, time delays reflect the combined effect of slow synaptic transmission and the finite propagation velocity of nerve impulses along axons. Moreover, delayed neural signaling has been shown to support a wide range of nonlinear phenomena in the brain, ranging from collective synchronization [4, 5], stochastic resonance [6] and induction of brain oscillations [7] to various types of noise-induced bifurcations [8–13]. For this reason, time delays have been implicated with a rich repertoire of physiological dynamics supporting brain function [12, 14–16].

In the nervous system, axonal conduction delays are not free parameters but are strictly regulated by myelin, a fatty substance enwrapping neural membranes. Myelin enables saltatory conduction, and greatly accelerates the propagation of action potentials (APs) along axons. Produced by glia cells, myelin compensates for variability in length and caliber between different axons [17–19], adjusting the time required for APs to propagate towards their post-synaptic targets (known as conduction delay). Conduction delays are of prime importance for the faithful transmission of correlated and/or temporally structured fluctuations resulting from recurrent connectivity [20, 21], sensory stimuli and/or functionally relevant afferent inputs. Such correlated stimuli challenge both the limited coding capacity of neurons[22], as well as the temporal coordination of neural signaling [17, 19, 21, 23]. To face this challenge, glial cells can adaptively adjust the amount of myelin along axons, speeding up or slowing down the propagation of APs by tuning conduction velocity, thereby stabilizing neuronal trafficking [12, 13, 24–27]. Indeed, myelin has recently been shown to be plastic, supporting and imparting flexibility for various neural computations, notably by optimizing conduction delay distributions amongst neuronal assemblies [27–31] with important consequences on neural dynamics, synchrony and brain function[13, 21, 27, 32, 33]. Research on glial cells, especially myelin, has emerged as one of the most effervescent and rapidly evolving areas in neuroscience (see [34] and references therein).

In the light of the manifest richness of delayed feedback systems [3], it is no surprise that myelin-induced changes in axonal conduction velocity (and delays) influence neural signaling and transform the dynamics of neural circuits and how they respond to time-varying fluctuations. Yet, networks with time delays are amongst the most
challenging to characterize mathematically, which severely limits data interpretation. We here explore this problem using a spiking neural network of excitatory neurons, with distributed delay and correlated additive noise. In our analysis, we considered conduction velocities over a range aligned with experimental observations made in the central nervous system of across multiple mammalian species [19, 35, 36].

We develop a computational strategy for stochastic spiking networks, and associated non-Markovian theory, that lifts the restriction of single delays to reveal a host of new delay density-dependent dynamical phenomena of physiological and clinical relevance. Leveraging mean field theory, we develop a method to determine the statistics of the stationary probability density function associated with the collective neuron activity in this network. Building on previous results on stochastic delay equations [37–41], this approach permits the derivation of implicit expressions for statistical moments of the mean field dynamics. Most importantly, it reveals that delay distributions may regulate the activity of cells in these network by tuning their individual dynamic range and encoding capability for noise-like stimuli. Using this approach, we characterize changes in neural activity resulting from variation in axonal conduction velocity reflecting various levels of myelination.

While it is generally assumed that delays mainly influence synchrony and oscillations, our mathematical and numerical results show that myelin-induced changes in conduction velocity implement a gain control mechanism amplifying neural activity, independently of synaptic coupling, and/or irrespective of the presence of oscillations. This gain control mechanism was further found to be especially salient in the presence of dynamic stimuli by amplifying the collective response of the neurons to time-varying signals. Changes in conduction velocity were found to shape the information coding capacity of the network, suggesting that myelination may be tuned in a contextual way to optimize the faithful encoding of stimuli of a given amplitude and baseline noise into network firing properties [21, 42]. In addition, power spectral analysis reveals that the dispersion of time delays, resulting from either variability in conduction velocities and/or axonal lengths, regularizes neural activity, suppressing oscillatory activity and promoting asynchronous dynamics. These findings echo effects seen experimentally under the influence of activity-dependent myelination [17, 24, 26, 27, 29]. Taken together, these results allow us to draw conclusions about how myelination may solicit nonlinear neural interactions to implement gain control and optimize the encoding capability of neural circuits. Further, they enable us to make predictions about potential dynamical consequences of myelin disruption in pathological states, such as those seen in epilepsy [43] and/or multiple sclerosis[44].

## Results

Our analysis starts with the observation that the collective activity of neurons in our network model (cf. Methods section) is highly sensitive to myelination. Exposing individual neurons to a spatially correlated stimuli, a consistent increase in the mean network firing rate could be observed whenever conduction velocity - reflecting myelination - is increased. As shown in Fig. 1, faster conduction was found to amplify neural activity and firing rates, and the network response to stimuli became more salient,
<figure><img src="image_2.png" /><figcaption>**Figure 1** Gain control induced by axonal conduction velocity A. Network of Poisson spiking neurons with low ($c = 0.5$m/s) conduction velocity, reflecting poor myelination. The diagram illustrates neurons (black circles) connected through axons (edges) expressing different conduction delays due to various degrees of myelination (illustrated in shades of red). A raster plot of the network spiking activity shown alongside the network mean firing rate above shows the response of the network to a spatially correlated stochastic stimulus $S(t)$. The histogram corresponds to the firing rate distribution of the network (right insert). B For the same stimulus $S(t)$, if conduction velocity is increased ($c = 50$m/s) to reflect higher myelination, the response of the network is amplified and mean firing rates increase. This also manifests as a change in the network firing rate distribution, which is now heavily skewed (right insert). In A and B, axonal conduction delays $\tau_{ij} = l_{ij}/c$ are Gamma distributed with shape $k = 2$ and scale $\theta = 4/c$ ms and axonal length distribution $l_{ij}$. Except from conduction velocity (i.e., c), all other parameters are identical. The network is fully (i.e., all-to-all) connected and synaptic weights are both positive and identical. The other parameters in the model (cf. Methods section) used for these exemplar simulations are $N = 100$, $J_{ij} = \bar{J} = 0.9$, $I = 0$ and $D = 0.05$. The activation function is nonlinear with parameters $\beta = 2$, $h = 0.5$.</figcaption></figure>

suggestive of a gain control effect. We note that this effect is entirely due to a shift in delay distribution within the network: aside from conduction velocity, all other parameters - such as synaptic weights - were kept constant. The presence of this response amplification could also be confirmed by computing the network mean firing rate distribution, where faster conduction promotes increased probability of high firing rates (see below).

## Mean field statistics

To better understand this phenomenon, we examined how conduction velocity (i.e. c), reflecting the amount of myelin covering axons, shapes the dynamical properties of cells in the network. This was accomplished by computing the probability density function of the network-averaged instantaneous potential $\bar{u}$, a quantity that fluctuates in time as determined by the network collective dynamics. This was done for two different categories of delay distributions: Gamma, and uniformly distributed. This was
A $\tau_{ij} \rightarrow \Gamma(k, \theta/c)$

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

B $\tau_{ij} \rightarrow U[0, \ell_{max}/c]$

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

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

<figure><img src="image_6.png" /><figcaption>**Figure 2** Axonal conduction velocity shapes delay distributions and corresponding mean field statistics A Mean field probability density function (top) computed numerically for a network of spiking neurons with Gamma-distributed delays (bottom). The mean field probability density function may be well approximated by a Gaussian (solid lines), and its mean and variance changes as conduction velocity increases. The associated delay distributions are plotted alongside their analytical probability distributions. Here, the shape parameter is $k = 2$ and the scale $\theta = 4/c$ ms. B Corresponding mean field distribution for uniformly distributed delays (top) and associated delay distributions (bottom). Here, $l_{max} = 10$mm. Numerical simulations [45] of the network ensemble average $\bar{u}$ (shaded areas) are shown alongside normal distributions with associated mean and variance (curves) for $c = 5.0$ m/s (dark red), $c = 2.5$ m/s (red) and $c = 1.0$ m/s (black), respectively. Simulations were performed using the Euler-Maruyama scheme, with a population of $N = 100$ neurons over trials of duration $T = 500$ms with steps of 1ms and $D = 0.1$. Other parameters are $J_{ij} = \bar{J} = 0.9$, $I = 0$, and the activation function is nonlinear with parameters $\beta = 2$, $h = 0.5$.</figcaption></figure>

first computed numerically through simulations (see Subsection *Mean Field Moments* in the Methods section for details of the mathematical analysis), which then motivated our theoretical work. Within the range of parameters examined here, Fig. 2 indicates that the probability density function of the network mean may be well approximated by a Gaussian - a feature that unlocks a variety of properties that makes this problem amenable to a thorough analytical characterization. These simulations also reveal changes in mean field statistics: increasing the conduction velocity increases both the mean and variance of the network mean potential, or conversely, larger conduction delays tend to suppress both the mean and variance of the network mean activity.

As a consequence of this Gaussian property, high-order moments of the probability density function (i.e., $\rho_S(\bar{u})$) can be assumed to be small, so that the mean field $\bar{u}$ in
Eq. (2) may be used to develop equations for the dynamics of the mean (i.e., $\mu$) and variance (i.e., $\sigma^2$). This allows the derivation of the stationary distribution without relying on a Fokker-Planck formalism [38, 40, 41, 46]. We used this strategy to compute the stationary mean and variance of the mean field and characterized their respective dependence on myelin - that is, on conduction velocity and resulting delay distribution - and on the noise variance $D$. As shown in Fig. 3A-B for both Gamma- and uniformly distributed delays, conduction velocity increases both the mean and variance of the mean field $\bar{u}$, confirming the numerical observations in Figs. 1 and 2. While the changes in mean and variance induced by faster conduction velocity portrayed in Fig. 3 might appear small, they in fact support significant variations in mean firing rate (see Fig. 1 and below). We could further confirm the alignment between our analytical derivations and analytical results, which were found to be in good agreement over multiple orders of magnitude of the noise amplitude. We note that deviations between our analytical and numerical results were however observed for higher values of the noise variance, as the system approaches a regime characterized by stochastic/coherence resonance [21, 47]. However, it is interesting to note that the theoretical predictions for the variance were accurate over a broader range of noise variance compared to the mean.

The changes in statistics portrayed in Figs. 1-3 together with our theoretical analysis suggest that myelination impacts the recurrent interactions of the network by modulating neural activity and increasing the saliency of stimuli through augmented responses. This results directly from the expression of mean field moments. Notably, the stationary variance $\sigma_s^2$ is shown to scale with both the stimulus amplitude $D$, and the summed contribution of delay-dependent correlations $C_k$ in mean potential measured across the network (see Methods section; Eq. (9)). As conduction velocities are decreased, these delay-dependent correlations are suppressed: the network behaves as if it was effectively disconnected (i.e., $\sigma_S^2 = D$). Consequently, by increasing conduction velocity, myelination solicits recurrent interactions in the network, promoting correlated activity and amplifies the temporal coordination between cells. We note that this gain effect occurs here because the network is purely excitatory (see Discussion).

## Gain control and susceptibility depend on conduction velocity

Taken together, the above results suggest that neural networks may utilize changes in myelination and axonal conduction velocity, as well as resulting shifts in activity portrayed in Fig. 3, to optimize the representation of information about stimuli - such as noise variance and bias - which in turn could be favorably read out by a downstream network. Indeed, as we now show, the variance $D$ of the spatially correlated noise input and the input strength $I$ influence network firing statistics and sensitivity to these variables in a manner that is contingent on the distribution of delays between neurons, and specifically, on the mean conduction velocity. The aforementioned changes in stationary mean and variance of the network ensemble average $\bar{u}$ in Fig. 3 may hence be recast as the manifestation of a gain control effect, in which the saliency of neuronal responses can be amplified through adjustments in axonal conduction velocity [21].

To characterize the dynamic range and sensitivity of the network, we computed the susceptibility or gain $G$ of the network for different conduction velocities. This susceptibility quantifies how the mean firing rate changes with the noise variance (i.e. $D$) or
<figure><img src="image_7.png" /><figcaption>**Figure 3** Mean field stationary statistics. A Stationary mean of the network ensemble average $\bar{u}$ computed numerically (triangles) plotted alongside the stationary mean $\mu_S$ computed analytically (dotted lines, see Methods section Eqs. (18), (22)) for Gamma distributed delays, and across various conduction velocities considered ($c = 5.0$m/s, dark red; $c = 2.5$m/s, red; $c = 1.0$m/s, black) and various orders of noise variances $D$. The grey shaded area is bounded by the limit cases where there are no delays (top, $c \to \infty$) or alternatively if delays become infinite (bottom, $c \to 0$). Both of these limits were computed analytically. Significant deviations between numerical and analytical results emerge for strong noise variances $D$ (marked by stars). B Equivalent plot for uniformly distributed delays. C Stationary variance (i.e., $\sigma^2$) of the network ensemble average $\bar{u}$ for Gamma distributed delays, computed both numerically (triangles) and analytically (dotted lines), see Methods section for details. Here, deviations appear for large noise variances $D$ (marked as stars). D Same as panel C for uniformly distributed delays. In panels A and C, $k = 2$ and $\theta = 4$m/cm. In panels B and D, $l_{max} = 10$m. Simulations were performed using the Euler-Maruyama scheme, with a population of $N = 100$ neurons over trials of duration $T = 500$ secs with steps of 1ms. Other parameters are $J_{ij} = \bar{J} = 0.9$, $I = 0$, and the activation function is nonlinear with parameters $\beta = 2$, $h = 0.5$.</figcaption></figure>

the constant input strength (i.e. $I$) for different values of $c$. This gain corresponds to the slope of the network activation function $F$ averaged over the probability density of the mean potential (cf. Eq. (21) in the Methods section). We remind here that as $D$ is swept across the shown range, the values of the mean $\mu_S$ and variance $\sigma_S$ have to be continually and self-consistently recomputed in order to evaluate the gain function.

As can be seen in Fig. 4A, the mean firing rate scales monotonically with noise amplitude, as one intuitively expects. Moreover, as conduction velocity $c$ increases, the mean firing rates of the individual neurons - and thus of the whole network - increase as well, a direct consequence of the phenomenon portrayed in Fig. 3. We note that despite seemingly small changes in mean membrane potential (i.e., Fig. 3), firing rates
<figure><img src="image_8.png" /></figure>

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

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

<figure><img src="image_11.png" /><figcaption>Figure 4 Delay-induced firing and gain control as a function of noise and input strength.</figcaption></figure>

A Mean network firing rate as a function of noise variance *D* for various axonal conduction velocities (*c* = 5.0m/s, dark red; *c* = 2.5m/s, red; *c* = 1.0m/s, black). Firing rates are computed by multiplying the network firing probability *F* by a factor of 100. These rates expectedly scale with *D*, but are also amplified by higher conduction velocity. The grey shaded area is bounded by the limit cases where there are no delays (top, *c* → ∞) or alternatively if delays become infinite (bottom, *c* → 0). Both of these limits were computed analytically. The theory (dashed line) agrees well with the simulation of the full spiking network. B Network gain *G* or susceptibility as a function of *D* for different conduction velocities. These theoretical curves represent the slope of the curves in (A), and reveal an optimal noise strength that decreases with conduction velocity. Grey lines represent the same cases as in panel A. C Mean firing rate as a function of constant input strength *I* for different values of *c*. Theory (dashed line) is accurate up to approximately *I*=0.045, after which discrepancies grow due to the nonlinear behaviour of the network and the fact that the theory does not account for increased spiking noise variance in the network simulations. D Susceptibility as a function of input strength *I* for different *c* values. The theory is accurate again over the same range as in C. Network simulations for higher values of *I* reveal a susceptibility peak due to the inflection point present in panel C. The model assumes a distribution of axonal lengths and delays are Gamma-distributed with *k* = 2 and *θ* = (4mm/c). Other parameters are *J* = 0.9, and the activation function is nonlinear with parameters *β* = 2, *h* = 0.5.

experience important changes, increasing by about 20% as conduction velocity is in-
creased from 1m/s to 5m/s, and even more if conduction velocity is increased further.

Such gain control effect must however be regarded from the perspective of the neu-
rons' limited dynamic range and coding capacity: stimuli that are too weak may fail
to elicit a visible postsynaptic response in terms of network firing rate. If their ampli-
tude is too strong, such stimuli may in contrast lead to saturation (i.e., postsynaptic
neuron fires continuously). As direct manifestation of this, Figure 4B shows that the
gain peaks as a function of the noise strength, and this peak shifts to smaller values
of D as the velocity increases. This indicates that while myelin might control the re-
sponse gain, it is the stimulus amplitude D that dictates whether the susceptibility is
<figure><img src="image_12.png" /><figcaption>**Figure 5** Stabilizing effect of delay variability. A-C Gaussian delay distributions (left) determine the shape of the mean field power spectral density $S(\omega)$ (right). The Gaussian mean and variance are independent from each other. Here, the mean delay is kept constant and given by $m = k\theta/c$, while the variance is increased from $s^2 = 0.001$ to 0.5 (lines with increased shading). To compute the power spectral density, we assumed a stimulus $S(t)$ with noise variance $D = 0.05$ and zero bias i.e., $I = 0$. The dashed line corresponds to the Lorentzian power spectral density associated with $\tilde{J} = 0$. Increasing the conduction velocities shifts the mean delay, and are given by $c = 1.0\text{m/s}$ (black lines, panel A); $c = 2.5\text{m/s}$ (red lines, panel B); $c = 5.0\text{m/s}$ (dark red lines, panel C). D-F Eigenvalues associated with the characteristic polynomial (cf. Eq. (26) in the Methods section) whenever conduction velocity $c$ is fixed and the delay variance $s^2$ is varied. Increasing the delay dispersion $s^2$ shifts the oscillatory eigenvalues away from the imaginary axis and stabilizes the dynamics. Different conduction velocities are given by $c = 1.0\text{m/s}$ (black shading, panel D), $c = 2.5\text{m/s}$ (red shading, panel E) and $c = 5.0\text{m/s}$ (dark red shading, panel F). Other parameters are $\tilde{J} = 0.9$, $k = 2$, $\theta = 4$ mm, $I = 0$, and the activation function is nonlinear with parameters $\beta = 2$, $h = 0.5$.</figcaption></figure>

increased or decreased over the neurons dynamic range. This suggests that myelination may be tuned in a contextual way to optimize the faithful encoding of stimuli of a given amplitude and baseline noise into network firing properties [21, 42].

Figure 4C shows the corresponding results when the (constant) input bias (i.e., $I$) is increased. While the mean firing rate expectedly changes as $I$ is increased, the effect of conduction velocity and/or delay distribution does not appear to impact the response significantly compared to the noise strength $D$. The same trend can be seen in Fig. 4D, where the susceptibility monotonically increase with input strength but does not vary significantly with conduction velocity. These observations suggest that conduction velocity mainly influences the response of the network to dynamic, time-varying stimuli and not constant inputs.

These observations together suggest that conduction velocity may be tuned in a contextual way to optimize the faithful encoding of stimuli into network firing properties [21, 42]. We further note the good agreement of our theory with the numerical simulations of the spiking network, but only as long as $I$ does not exceed the range
shown. Beyond that, the susceptibility also peaks, albeit with a relative invariance to
conduction velocity. We hypothesize that this is due to the system being pushed into
a more nonlinear regime, in particular one in which neglecting the increase in Pois-
son noise leads to a substantial error. Coherence resonance may further amplify this
discrepancy (see below).

Conduction delay variability stabilizes network dynamics

We concluded our study by examining the effect of delay dispersion on the response
of our network to time-varying stimuli. To do so, we computed the power spectral
density $S(\omega)$ of the mean field $\bar{u}$ for various delay distributions. For that purpose,
we first relaxed our previous assumption by which eigenvalues associated with the
dynamics of our network are purely real. We examined purely stochastic inputs (i.e.
setting $I=0$) and deliberately considered a limit case in which the delay distribution
can be approximated by a Gaussian, so that the mean and variance may be modu-
lated independently from one another (see Methods section *Power spectral density*).

As can be seen in Fig. 5A-C, increasing the mean delay - i.e., decreasing the degree
of myelination and thus the conduction velocity - expectedly introduces resonant os-
cillatory modes in the network, whose frequencies scale with conduction velocity. Such
frequencies can easily be determined through linear stability analysis, corresponding
to imaginary eigenmodes. Increasing the variance of the delay distribution (while
keeping the mean delay constant) suppresses these resonances, as the power spectral
density converges towards a purely Lorentzian distribution $S(\omega) = D (\pi(1 + \omega^2))^{-1}$.
We highlight that the latter corresponds to the uncoupled case where $J_{ij} = 0$, in
which the mean field in Eq. (23) simplifies to a non-delayed Langevin equation. To
better understand the mechanism at play, we performed linear stability analysis of our
mean field system, computing explicitly the eigenvalues and how they depend on con-
duction velocity and delay variance (see Eq. (26)). As shown in Fig. 5D-F, increasing
the dispersion of the delays while keeping their mean constant shifts complex eigen-
value away from the imaginary axis. This indicates that delay dispersion quenches
the oscillatory modes of the system, making the system increasingly linear. These ob-
servations echo previous findings [48–50] showcasing the strongly stabilizing influence
of delay dispersion, by which the system becomes effectively decoupled through the
suppression of oscillatory modes, thereby promoting irregular, asynchronous activity.

Discussion

The influence of time delays on nonlinear stochastic systems has been thoroughly
studied [38, 40, 41], especially from the perspective of synchrony [48, 51, 52], yet their
consequence on the computational properties of neural networks has received compar-
atively less attention. Action potential conduction velocity is strongly influenced by
the presence (and amount) of myelin covering the axon, which impacts the temporal
coordination of neural signaling with important consequences on network dynamics
[13, 18, 19, 21, 27, 29]. Glial-mediated changes in conduction velocity along axons have
been shown to be titrated autonomously to promote AP coincidence [17], reinforcing
the probability of post-synaptic response and hence to enhance information transfer [21, 53] as well as the resilience of brain function [13, 26]. We leveraged mean field theory for stochastic delay equations [47, 54–56] to investigate how conduction velocity - through its influence on conduction delays - influences the activity of recurrent networks and response to stimuli, in order to unveil the potential functional role of myelin in neural computations.

Using mean field theory, we analyzed an excitatory network of spiking neurons with distributed delay and additive, spatially correlated stimuli. The shape of the delay distribution in this network was allowed to change with conduction velocity, meant to reflect the net effect of myelin on AP propagation. Building on a set of assumptions supported numerically, we derived the stationary moments for the probability density function of the network ensemble average. Considering these delay distributions, we quantified the mean response of the network and its variability to changes in stimulus bias and noise strength. This analysis revealed that conduction velocity - through its effect on delay distribution - operates as a gain control mechanism, amplifying network firing rates to shared (correlated) input. Indeed, faster conduction velocity decreases delay dispersion (i.e., rendering conduction delays more similar to each other), increasing correlations in neural activity and, consequently, network firing rates. The network gain (i.e. susceptibility) was found to exhibit a peak as a function of noise strength as well as input strength, revealing that conduction velocity may be tuned to optimize the response of the network to these changes. Further, short delays, i.e. slow conduction, was found to promote oscillatory activity in the recurrent network. Our theory reproduces the gain properties seen in our numerical simulations of the spiking network, provided that the system is not pushed too far into its nonlinear regime.

Taken together, our computational and theoretical results suggest that myelinating glia in even simplified recurrent neural networks can play an important role in controlling firing rate statistics, along with their sensitivity to input changes beyond the induction of oscillations. Factoring in delay distributions also expands our knowledge about the determinants of firing rate distributions measured in real neural systems [57, 58]. Over the past decades, there has been a surge of interest in understanding the role played by myelin in brain function, health and disease [26, 27, 30, 59–61]. Oligodendrocytes have since been implicated with a plethora of biological functions across the CNS, ranging from learning, memory consolidation, stimulus discrimination to functional resilience [13, 28, 31–33, 59, 61, 62]. Yet, understanding how myelin - and by extension glia - shape the computational properties of neural networks remains a formidable challenge, notably because of the complexity of non-linear delayed feedback systems.

In light of these results, both insufficient and excessive conduction velocity (i.e., myelination) represent potential disruptors of neural function. Recent experiments have causally linked excessive myelination to seizures in mice models of generalized epilepsy[43]. The gain control effect described here may thus be involved in supporting pathologically high firing rates and correlated neural activity, both of which accompany seizures [63]. Conversely, our model also predicts that decrease in conduction velocity predisposes the network to lower firing rates and gain, which may be associated with impairments observed in demyelinating diseases such as multiple sclerosis
(MS) [28, 44, 64–67].

Our results are an extension to the distributed delay case of previous theoretical work on oscillations in single delay rate-based models with finite-size, spatially uncorrelated noise [68], as well as in leaky integrate-and-fire stochastic networks without [69] and with [70] spatial noise correlations. Although our work has considered that the delay distribution resulted from a distribution in axonal lengths, the theory can not distinguish this case from that where the delay distribution arises due to changes in axonal diameters but fixed lengths, since it is known experimentally [19] that the axonal conduction velocity is proportional to axonal diameter for myelinated fibers. The more realistic situation may involve both significant diameter and length variations, which could be taken into account in our theory.

Our approach suffers from limitations which should be developed in future studies. First, the network model considered here is purely excitatory. Future work addressing the impact of inhibition on our results is of great interest. How conduction velocity along both excitatory and inhibitory axons influences the dynamics of such networks would provide important information about the functional role of myelin, especially given excitatory and inhibitory neurons are known to display different myelination profiles [71, 72].

Moreover, we have introduced variability in conduction delays by sampling axonal lengths from a given probability (i.e., either Gamma or uniform) distribution, while keeping conduction velocity the same for all network axons. However, myelin distribution, caliber and/or g-ratio may vary along and between axons. Including this additional source of variability would lead to axon-specific conduction velocities. We note that, however, that resulting conduction delays would remain distributed - and hence the current results would remain qualitatively unaffected.

The network under study and associated mean field (see Methods section) are known to exhibit oscillations and/or stochastic/coherence resonance [21, 47, 51, 52, 54–56, 73]. Such dynamics yield non-Gaussian probability density profiles, and are further characterized by complex correlation rates $\gamma$. As such, our approach holds for positive values of synaptic weights ($\bar{J} > 0$) below the multistability threshold, as well as for small values of the response gain $\beta$ and noise amplitude $D$. We have also used a small delay expansion [41] to approximate the correlation rate $\gamma$, which would otherwise be both complex and satisfy a transcendental equation, greatly complicating the analysis. Despite this approximation, we highlight that our analysis remains in good agreement with the numerical simulations.

In addition, our network is fully connected (i.e., all-to-all) which was done for simplicity. We orient the reader to [21], where a similar network was examined with sparse connectivity. In our mean field equation (Eq. (2) in the Methods section), we have neglected fluctuations in membrane potential resulting from Poisson spiking. Such random fluctuations, whose amplitude scale with firing rate, have a multiplicative effect on the dynamics of Eq. (1) and may have profound influence on network stability [37]; it likely explains the limitation of our theory at high input strength. We hypothesize that this additional source of variability might be involved in the deviations observed for high noise intensity in Fig. 3A-B. Lastly, our model is abstract and neglects multiple relevant physiological details in order to preserve mathematical
tractability. Further work is required to augment the biological accuracy of our model and its ability to mimic real brain circuits and the influence of glia on neural function.

Methods

Model

We considered a generic and widely used class of spiking recurrent neural networks in which the potential $u_i$, $i = 1...N$, of neurons is governed by the set of stochastic delayed equations with additive correlated noise

$$ \tau_s \frac{d}{dt} u_i = L[u_i] + \frac{1}{N} \sum_{j=1}^{N} J_{ij} X_j(t - \tau_{ij}) + S(t), \quad (1) $$

where $\tau_s$ corresponds to the synaptic time constant, $L[u_i] = -u_i$ is a linear relaxation operator, and $J_{ij}$ are synaptic weights. Spike trains $X_j = \sum_{t_{ij}} \delta(t - t_j)$ are assumed to obey non-homogeneous Poisson processes with rate $r[u_j]$, where $r$ is a given sigmoidal response function [54] and $\delta(\cdot)$ is the Dirac delta function. An additive, spatially correlated input $S(t) = I + \sqrt{2D}\xi(t)$ is interpreted as being a functionally meaningful stimulus applied equally to all neurons, for instance resulting from coherent afferent input. This stimulus is composed of an input $I$, which is constant over the network and over time, and a spatially homogeneous Gaussian white noise term $\xi(t)$ of variance $2D$, satisfying $\langle \xi(t)\xi(t') \rangle = \delta(t - t')$. Here, $\langle \cdot \rangle$ denotes the ensemble average. Since it is additive, the stimulus $S(t)$ does not depend on any of the system's state variables, so that $\langle S(t) \rangle = I$.

Equation (1) describes the dynamics of a spatially distributed network wherein neurons can either be clustered spatially (e.g., cortical microcircuit) or distant from each other (e.g., located in different brain areas). As such, the interactions within this network occur via axonal pathways which may extend globally through white matter, or alternatively locally between cortical layers, for instance. Due to this spatial separation, action potentials transit along axons of various lengths $\ell$ with axonal conduction velocity $c$ and thus experience a conduction delay $\tau = \ell/c$. The axonal conduction velocity depends on the amount of myelin covering the corresponding axon. Specifically, resulting conduction delays $\tau_{ij} > 0$ between two neurons $i$ and $j$ in Eq. (1) are assumed to be random and collectively sampled from a probability distribution $q(\tau)$ (see section below).

Under specific assumptions [21, 47, 54, 55] and whenever $\mathbb{E}_{N,N}[J_{ij}] \equiv \bar{J} \neq 0$, Eq. (1) admits a well defined mean field for the network ensemble average potential $\bar{u} \equiv \mathbb{E}_{N}[u_i]$. Here $\mathbb{E}_{N,N}[J_{ij}]$ and $\mathbb{E}_{N}$ are the network average across the $N \times N$ connections and $N$ neurons of the network, respectively. Then

$$ \tau_s \frac{d}{dt} \bar{u} = L[\bar{u}] + \frac{\bar{J}}{K} \sum_{k=1}^{K} f[\bar{u}_{\tau_k}] + S(t) \quad (2) $$
for $\bar{u}_{\tau_k} = \bar{u}(t - \tau_k)$ and conduction delay times $\{\tau_k\}$, $\tau_1 \le \tau_2 \le \dots \le \tau_K$. We orient the reader to the Appendix and [21, 54, 55] for a thorough derivation. In this mean field equation (2), we have introduced the network average response function $f[\bar{u}] \equiv \mathbb{E}_N[r]$ [47, 54, 55] and expressed the sum over the $K = N^2$ conduction delays for notational convenience. The delay times are linked to the delay distribution by $q(\tau) = \sum_{k=1}^K \delta(\tau - \tau_k)/K$.

Equation (2) is our starting point. As defined, this mean field corresponds to a nonlinear stochastic delay differential equation (SDDE), a well-studied [38, 40, 41, 46, 51, 52, 73], yet challenging family of dynamical systems. We are here interested in determining the first two moments of the stationary probability density function $\rho_S(\bar{u})$ for the mean activity $\bar{u}$, and how such moments depend on the delay distribution $q(\tau)$. By doing so, our goal is to complement our understanding of the generic effect of delay distribution on the dynamics of non-linear networks, and formulate conclusions about the underlying functional role of myelin towards gain control in neural circuits.

## Mean Field Moments

In what follows, we set the membrane time constant $\tau_s = 1$ for simplicity and consider a purely excitatory network, i.e. $\bar{J} > 0$, for which the dynamics of Eq. (2) possess a single equilibrium. This regime holds whenever the response gain $df/d\bar{u}$ is small, and whenever synaptic weights $\bar{J}$ remain below the multistability threshold. Previous studies have successfully used approximating approaches by which the Fokker-Plank equation for the probability density function $\rho(\bar{u}, t)$ may be formulated in presence of time delay, despite non-Markovian dynamics [38, 40, 41, 46]. We here opt for an alternative (yet still approximate) strategy. We will assume that the a priori unknown stationary probability density function $\rho_S(\bar{u})$ is normally distributed with mean $\mu_S$ and variance $\sigma_S^2$, that is

$$ \rho_S(\bar{u}) \approx \mathcal{N}(\mu_S, \sigma_S^2) \qquad (3) $$

This assumption implies that the stationary state of Eq. (2) in absence of noise is non-oscillatory, and that associated eigenvalues of the corresponding linearized system are purely real. While Eq. (2) may admit imaginary eigenvalues, this approximation remains valid over a wide range of parameter values, and is well supported by numerical observations (see Fig. 2C-D). By virtue of this Gaussian approximation, all high order moments of $\rho_S(\bar{u})$ vanish, and we may thus characterize the mean field $\bar{u}$ in Eq. (2) by developing explicit differential equations for the first and second moment, whose fixed points correspond to the stationary statistics of $\rho_S(\bar{u})$.

The first moment, defined by

$$ \mu(t) \equiv \langle \bar{u}(t) \rangle = \int_{-\infty}^{\infty} \bar{u} \rho(\bar{u}, t) d\bar{u}, $$

is the mean potential at time $t$. Taking this expectation over the density of the mean potential at time $t$ on both sides of Eq. (2), the time-dependent mean potential $\mu(t)$ is found to obey the following delayed (noise-free) differential equation:
$$ \frac{d}{dt}\mu = L[\mu] + \frac{\bar{J}}{K} \sum_{k=1}^{K} \langle f[\bar{u}_{\tau_k}] \rangle + I, \quad (4) $$

since $\langle S \rangle = I$. We here introduce the expectation of the activation function $f$.

$$ F[\mu, \sigma^2] \equiv \langle f[\bar{u}] \rangle = \int_{-\infty}^{\infty} f[\bar{u}] \rho(\bar{u}, t) d\bar{u}, \quad (5) $$

which also corresponds to the ensemble average firing rate of neurons in the network.

One may hence determine the stationary mean $\mu_S$ through the equilibria of Eq. (4)

above, which are implicitly defined by

$$ \mu_S = \bar{J}F[\mu_S, \sigma_S^2] + I. \quad (6) $$

Similarly, the second moment defined by

$$ \langle \bar{u}^2(t) \rangle = \int_{-\infty}^{\infty} \bar{u}^2 \rho(\bar{u}, t) d\bar{u} $$

can be shown to obey the following differential equation [39, 41]

$$ \frac{d}{dt}\langle\bar{u}^2\rangle = 2L[\langle\bar{u}^2\rangle] + \frac{2\bar{J}}{K}\sum_{k=1}^{K}\langle\bar{u}f[\bar{u}_{\tau_k}]\rangle + 2\langle\bar{u}I\rangle + 2\sqrt{2D}\langle\bar{u}\xi\rangle, \quad (7) $$

applying the chain rule $\frac{d}{dt}\langle\bar{u}^2\rangle = 2\langle\bar{u}\frac{d\bar{u}}{dt}\rangle$. Introducing the covariance terms

$$ C_k[\mu, \sigma^2] = \langle \bar{u} f[\bar{u}_{\tau_k}] \rangle, $$

one may hence consider the equilibrium of Eq. (7) to obtain an expression for the stationary second moment

$$ \langle \bar{u}^2 \rangle_S = \frac{\bar{J}}{K} \sum_{k=1}^{K} C_k[\mu_S, \sigma_S^2] + \langle \bar{u}I \rangle_S + \sqrt{2D} \langle \bar{u}\xi \rangle_S \quad (8) $$

where the stationary average $\langle U \rangle_S = \int U \rho_S(\bar{u}) d\bar{u}$ is performed over the stationary

probability density function $\rho_S(\bar{u}) = \lim_{t\to\infty} \rho(\bar{u}, t)$. The stationary value of the third

term on the right hand side of Eq. (8) may be approximated through the property of

the Ornstein-Uhlenbeck process [39]

$$ \langle \bar{u}\xi \rangle_S \approx \sqrt{\frac{D}{2}}. $$

Since $\sigma_S^2 = \langle \bar{u}^2 \rangle_S - \mu_S^2$ and $\langle \bar{u}I \rangle_S = \mu_S I$, replacing the above in Eq. (8) and

substracting $\mu_S^2$ on both sides gives us an implicit expression for the stationary
439 variance,

$$
\sigma_S^2 = \frac{\bar{J}}{K} \sum_{k=1}^{K} C_k[\mu_S, \sigma_S^2] + \mu_S I + D - \mu_S^2 . \qquad (9)
$$

440 The covariance terms $C_k[\mu_S, \sigma_S^2]$ correspond to the average over the stationary joint

441 probability density function $\rho_S(\bar{u}, \bar{u}_{\tau_k})$ [40]

$$
C_k[\mu_S, \sigma_S^2] = \int_{-\infty}^{+\infty} \int_{-\infty}^{+\infty} \bar{u} f[\bar{u}_{\tau_k}] \rho_S(\bar{u}, \bar{u}_{\tau_k}) d\bar{u} d\bar{u}_{\tau_k} . \quad (10)
$$

442 The above integral is to be performed over two random variables $\bar{u}$ and $\bar{u}_{\tau_k}$, which are

443 assumed to be Gaussian with mean $\mu_S$, variance $\sigma_S^2$ and covariance $\langle \bar{u}\bar{u}_{\tau_k} \rangle = \sigma_S^2 p_k$,

444 where $p_k = e^{-\gamma\tau_k}$ is the autocorrelation and $\gamma$ is the correlation rate (see Appendix B

445 for more details). We discuss how to determine the correlation rate in further detail

446 in the subsequent sections below. Hence, as per Eq. (3), the joint probability density

447 function may be written explicitly as

$$
\begin{align}
\rho_S(\bar{u}, \bar{u}_{\tau_k}) &= \rho_S(\bar{u}|\bar{u}_{\tau_k})\rho_S(\bar{u}_{\tau_k}), \tag{11} \\
&= \frac{e^{-\left[\frac{Q}{2(1-p_k^2)} - \frac{(\bar{u}_{\tau_k}-\mu_S)^2}{2\sigma_S^2}\right]}}{2\pi\sigma_S^2} \sqrt{1-p_k^2}, \nonumber
\end{align}
$$

448 where $Q = ((\bar{u} - \mu_S)/\sigma_S - p_k(\bar{u}_{\tau_k} - \mu_S)/\sigma_S)^2$. Let us now examine specific examples.

449

Linear Networks

450

To illustrate the approach above, let us first consider the simplistic case in which the
response function is given by $f[\bar{u}] = \bar{u}$ with the constant input $I = 0$. The mean field
dynamics in Eq. (2) may be rewritten as

451
452

453

$$
\frac{d}{dt} \bar{u} = L[\bar{u}] + \frac{\bar{J}}{K} \sum_{k=1}^{K} \bar{u}_{\tau_k} + \sqrt{2D\xi}(t) \quad (12)
$$

454 The stationary mean $\mu_S$, which corresponds to the equilibrium of Eq. (12), is here

455

given by

$$
\mu_S = \bar{J}F[\mu_S, \sigma_S^2] = \bar{J}\mu_S = 0,
$$

456 since $\bar{J} > 0$. Determining the variance $\sigma_S^2$ of the stationary probability density function

457 requires to compute the K covariance terms $C_k[\mu_S, \sigma_S^2]$ which as per Eq. (10) simplify

458 in the linear case to

$$
C_k[\mu_S, \sigma_S^2] = \langle \bar{u} \bar{u}_{\tau_k} \rangle = \sigma_S^2 p_k = \sigma_S^2 e^{-\gamma \tau_k} .
$$
As such, since $\mu_S = 0$, Eq. (9) becomes

$$ \sigma_S^2 = \frac{\bar{J}\sigma_S^2}{K} \sum_{k=1}^{K} e^{-\gamma\tau_k} + D. \qquad (13) $$

We note that whenever $K \gg 1$, the above sum can be approximated by an integral to obtain our final expression for the linear case

$$ \sigma_S^2 = \bar{J}\sigma_S^2\mathcal{L}[-\gamma] + D, \qquad (14) $$

where $\mathcal{L}[-\gamma] = \mathbb{E}[e^{-\gamma\tau}] = \int_0^\infty e^{-\gamma\tau}q(\tau)d\tau$ is the Laplace transform of the delay probability density function $q(\tau)$. To determine the correlation rate $\gamma$, we first observe that the dynamics of Eq. (12) may be well described for $K \gg 1$ by the following Ornstein-Uhlenbeck process [41]

$$ \frac{d}{dt}\bar{u} = -\alpha\bar{u} + \sqrt{2D}\xi(t), \qquad (15) $$

where $\alpha = 1 - \bar{J}\mathcal{L}[-\gamma]$. This result follows from the analysis of the characteristic equation of Eq. (12), associated with the trivial fixed point $\bar{u}_o = 0$ in absence of noise, here given by

$$ \lambda = -1 + \bar{J}\mathcal{L}[-\lambda], \qquad (16) $$

where $\lambda \in \mathbb{C}$ is the eigenvalue. The eigenvalue $\lambda$ in the equation above generically depends on both the stationary mean $\mu_S$ and variance $\sigma_S^2$. However, for $\tau_k$ small, one may consider the first order expansion [41] $\lambda = -1 + \bar{J} + \epsilon$, $\epsilon \ll 1$ in Eq. (16), for which the correlation rate becomes $\gamma = -1 + \bar{J}$, which we note is purely real. Substituting into Eq. (14) yields an explicit expression for the variance in the linear case i.e.,

$$ \sigma_S^2 \approx \frac{D}{1 - \bar{J}\mathcal{L}[1 - \bar{J}]} = \frac{D}{\alpha}. $$

We note that this result corresponds to the approximation provided in previous work [41].

## 477 Nonlinear Networks

We now examine the case where the activation function in Eq. (2) is instead given by the sigmoid

$$ f[\bar{u}] = \frac{1}{2} (1 + \operatorname{erf}[\beta(\bar{u} - h)]) . \qquad (17) $$

Here, $\beta$ is the nonlinear response gain and $h$ the excitability threshold of the neurons.

Note that Eq. (17) corresponds to the individual mean firing rate of the neurons in
the network, i.e. $r_j = f[u_j]$. In this more involved case, the stationary mean is given by the following implicit equilibrium relationship

$$ \mu_S = \bar{J}F[\mu_S, \sigma_S^2] + I, \qquad (18) $$

where the functional $F$ at equilibrium reflecting the ensemble average firing rate in the network is given as per Eq. (5),

$$ F[\mu_S, \sigma_S^2] = \frac{1}{2} \left( 1 + \operatorname{erf} \left[ \frac{\beta(\mu_S - h)}{\sqrt{1 + 2\beta^2 \sigma_S^2}} \right] \right). \qquad (19) $$

We note that in contrast to the linear case, now the stationary mean of $\bar{u}$ implicitly depends on its variance. To compute $\sigma_S^2$, one needs to examine the covariance terms $C_k[\mu, \sigma^2]$ at equilibrium. Following a rather laborious technical derivation, one can obtain

$$ C_k[\mu_S, \sigma_S^2] = G[\mu_S, \sigma_S^2] \sigma_S^2 p_k + \mu_S F[\mu_S, \sigma_S^2] $$

with the average gain $G$ of the network

$$ G[\mu_S, \sigma_S^2] \equiv \langle f'(\bar{u}) \rangle = \int_{-\infty}^{\infty} f'(\bar{u}) \rho_S(\bar{u}) d\bar{u}. \qquad (20) $$

This mean network gain is expressed as an average of the single neuron gain $f'(\bar{u}) \equiv df/d\bar{u}$ over the density of the average potential. This gain of the network quantifies how the mean activity (i.e. firing rate) varies in response to changes in the mean potential, which itself is a function of all the parameters of the system. The network gain will be used later when gain control properties are investigated. It can be written in terms of network parameters and $(\mu_S, \sigma_S)$ as

$$ G[\mu_S, \sigma_S^2] = \frac{\beta \exp \left[ - \frac{\beta^2 (\mu_S - h)^2}{1 + 2\sigma_S^2 \beta^2} \right]}{\sqrt{\pi (1 + 2\sigma_S^2 \beta^2)}}. \qquad (21) $$

Replacing into Eq. (9) yields the variance of the mean potential

$$ \sigma_S^2 = \frac{\bar{J}\sigma_S^2 G[\mu_S, \sigma_S^2]}{K} \sum_{k=1}^{K} e^{-\gamma\tau_k} + \bar{J}\mu_S F[\mu_S, \sigma_S^2] + \mu_S I + D - \mu_S^2. $$

Using the fact that $\mu_S \bar{J} F[\mu_S, \sigma_S^2] + \mu_S I = \mu_S^2$ as per Eq. (18), we may consider $K \gg 1$ to obtain our final expression for the stationary variance in the nonlinear case

$$ \sigma_S^2 = \bar{J}\sigma_S^2 G[\mu_S, \sigma_S^2] \mathcal{L}[-\gamma] + D. \qquad (22) $$

We note that this equation together with Eq. (18) implicitly define the stationary probability density function in Eq. (3). In other words, $(\mu_S, \sigma_S)$ must be computed (numerically) as the solution to these two equations. If one was to scan a parameter in
the model, as we do below for the noise and input strength parameters D and I, the solution set would define a curve of such solutions in the $(\mu_S, \sigma_S)$ space, along which the functions F and G can be evaluated.

Here, the correlation rate $\gamma$ is still undefined. As in the section above, the correlation rate is also associated with the eigenvalue $\lambda \in \mathbb{C}$ of Eq. (2), which in the non-linear case corresponds to

$$ \lambda = -1 + \bar{J}f'[\bar{u}_o]\mathcal{L}[-\lambda], $$

where the fixed point $\bar{u}_o$ satisfies the implicit equilibrium relationship $\bar{u}_o = \bar{J}f[\bar{u}_o] + I$ for $D = 0$. We here adopt the same approach as in the linear case and consider $\tau_k$ small, for which the first order expansion yields the real-valued eigenvalue $\lambda \approx -1 + \bar{J}f'[\bar{u}_o] + \epsilon$ for $\epsilon \ll 1$, resulting in a correlation rate given by $\gamma = -1 + \bar{J}$ to the first order. Substituting into Eq. (22) leads to an implicit expression for the variance

$$ \sigma_S^2 = \bar{J}\sigma_S^2 G[\mu_S, \sigma_S^2]\mathcal{L}[1 - \bar{J}] + D. $$

The above expressions for the statistics of the mean field agree well with numerical simulations (see Fig. 3A-D).

## Myelin and conduction delay distributions

To model the effect of myelin, we examined the effect of different conduction velocities $c$ on the propagation of action potentials along axons of given lengths $\ell$. We considered two specific cases, in which resulting conduction delays $\tau = \ell/c$ adopt different distributions:

**Gamma distributed delays.** Let us assume that the $K = N^2$ axons separating neurons in Eq. (1) possess individual lengths $\ell_{ij}$ randomly sampled from a Gamma distribution i.e. $\ell_{ij} \rightarrow \Gamma[k, \theta]$, where $k$ and $\theta$ refer to the shape and scale parameters, respectively. Further, the APs propagate along those axons with constant conduction velocity $c$, a value that scales with level of myelination along axons. Resulting conduction delays $\tau_{ij} = \ell_{ij}/c$, for a given conduction velocity $c$ are hence also Gamma distributed with density $q(\tau) = \Gamma[k, \theta/c]$. We point out that the mean connection strength $\bar{J}$ is independent of the connection topology given by the distribution of connection lengths $\ell_{ij}$. In the limit of large $N$, the associated Laplace transform $\mathcal{L}[-\gamma]$ is given by

$$ \mathcal{L}[-\gamma] = \mathbb{E}[e^{-\gamma\tau}] = (1 + \theta\gamma/c)^{-k}. $$

**Uniformly distributed delays.** Alternatively, axons separating neurons may display individual lengths $\ell_{ij}$ randomly sampled from a uniform distribution i.e. $\ell_{ij} \sim U[0, \ell_{\max}]$. Then resulting conduction delays $\tau_{ij} = \ell_{ij}/c$ are also uniformly distributed
$$ \mathcal{L}[-\gamma] = \mathbb{E}[e^{-\gamma\tau}] = \frac{1 - e^{-\gamma\ell_{\max}/c}}{\gamma\ell_{\max}/c} . $$

## Power spectral density

$$ \frac{d}{dt} \bar{u} = L[\bar{u}] + \bar{J} f'[u_o] \int_0^\infty \bar{u}(t-\tau) q(\tau) d\tau + \sqrt{2D} \xi(t), \quad (23) $$

$$ i\omega\hat{u} = L[\hat{u}] + \bar{J}f'[u_o]\hat{u}\hat{q}(\omega) + \sqrt{2D}\hat{\xi}(\omega) \quad (24) $$

$$ \langle \hat{u}(\omega) \hat{u}^*(\omega') \rangle = \delta(\omega - \omega') S(\omega) $$

$$ S(\omega) = \frac{D}{\pi (\mathrm{Re}(P[0, \omega])^2 + \mathrm{Im}(P[0, \omega])^2)} . \quad (25) $$

$$ P[a, \omega] = a + i\omega + 1 - Jf'[\bar{u}_o]\mathcal{L}[-\lambda], \quad (26) $$

$$ P[0, \omega] = i\omega + 1 - Jf'[\bar{u}_o]\hat{q}(\omega) . \quad (27) $$

<footnote n="539">539 with density $q(\tau) = U[0, \ell_{\max}/c]$. In this case, the associated Laplace transform $\mathcal{L}[-\gamma]$ reads</footnote>

<footnote>540</footnote>

<footnote n="542">542 To gain better insight as to the effect of distributed delays on the response of the network to time-varying stimuli, we computed the power spectral density $S(\omega)$ of the mean field $\bar{u}$. Specifically, we linearized Eq. (2) around the deterministic equilibrium $\bar{u}_o$ to obtain the scalar SDDE,</footnote>

<footnote>543</footnote>

<footnote>544</footnote>

<footnote>545</footnote>

<footnote n="546">546 where we assumed $K \gg 1$. Taking the Fourier transform yields</footnote>

<footnote n="547">547 where $\hat{u} = \hat{u}(\omega)$ and $\hat{\xi}(\omega)$ is the Fourier transform of the mean field and the external driving noise, respectively. Given that $\langle \hat{\xi}(\omega) \hat{\xi}^*(\omega') \rangle = \delta(w - w')/2\pi$, where * denotes the complex conjugate,</footnote>

<footnote>548</footnote>

<footnote>549</footnote>

<footnote n="550">550 with the power spectral density</footnote>

<footnote n="551">551 Here, $\mathrm{Re}(P[a, \omega])$ and $\mathrm{Im}(P[a, \omega])$ correspond to the real and imaginary parts of the characteristic polynomial $P[a, \omega]$ associated with the eigenvalue $\lambda = a + i\omega \in \mathbb{C}$, respectively. The polynomial $P[a, \omega]$ associated to Eq. (23) for $D = 0$ is given by</footnote>

<footnote>552</footnote>

<footnote>553</footnote>

<footnote n="554">554 which simplifies for $a = 0$ to</footnote>

<footnote n="555">555 To separate the (potentially distinct) contributions of the mean delay versus delay variance on the power spectrum in Eq. (25), we considered Gamma distributed delays i.e., $q(\tau) = \Gamma(k, \theta/c)$ in the limiting case for large mean and variance. Indeed, whenever $k \gg 1$, the delay distribution $q$ approaches a normal distribution i.e. $q(\tau) \approx N(m, s^2)$ with mean delay $m = k\theta$ and variance $s^2 = k\theta^2$. This allows to express the Fourier transform $\hat{q}$ in Eq. (27) in a separable form i.e., $\hat{q}(\omega) = e^{-\frac{1}{2}\omega^2 s^2} e^{i\omega m}$. We relaxed</footnote>

<footnote>557</footnote>

<footnote>558</footnote>

<footnote>559</footnote>

<footnote>560</footnote>
the premises above, and assumed that the mean delay $m$ and variance $s^2$ may be modulated independently in a regime where $m > s$. While this clearly does not hold for the Gamma distributed case (since the mean and variance are not independent), this approach remains insightful for delay distributions, that adopt a Gaussian shaped profile. We note, that one may alternatively consider uniformly distributed delays within the interval $[m - \delta\tau, m + \delta\tau]$, in which case the delay distribution adopts a mean $m$ and variance $s^2 = 3\delta\tau^2$ instead. This yields qualitatively similar results as those portrayed in Fig 5.

## Declarations

### Funding

We thank the National Sciences and Engineering Research Council of Canada (NSERC Grants RGPIN-2017-06662) for support of this research.

### Data availability

The data, algorithms and numerical methods that support the findings of this study are available from the corresponding author upon reasonable request.

### Author contribution

J.L. designed the study. A.C. provided support for the simulations. J.L., A.H. and A.L. performed the analysis and wrote the manuscript.

## Appendix A Mean Field Derivation

The evolution equation (1) describes how the activity of single neurons evolve on short time scales. This short time scale is observed by the neurons' spike-coupling. The equation can be written as

$$ \tau_s \frac{d}{dt} u_i = L[u_i] + I_i + \sqrt{2D\xi}(t) \quad (A1) $$

with

$$ \begin{aligned} I_i &= \frac{1}{N} \sum_{j=1}^{N} J_{ij} X_j(t - \tau_{ij}) \\ &= \frac{1}{N} \int_{0}^{\infty} dT \sum_{j=1}^{N} J_{ij} \delta(T - \tau_{ij}) X_j(t - T) . \end{aligned} $$

Let us consider the eigenvalue spectrum of the connectivity matrix **J** [47] with the basis vectors {$\psi_j$}, {$\phi_i$}, $\psi_j, \phi_j \in C^N$ and the eigenvalues $\lambda_j \in C$, which obey

$$ \psi_j^\dagger \phi_i = \delta_{ji} $$
$$ \psi_j^\dagger \mathbf{J} = \lambda_j \psi_j^\dagger . $$

Here, $\dagger$ denotes the complex conjugate transposition and $\delta_{ji}$ is the Kronecker symbol with $\delta_{ii} = 1$ and $\delta_{ij} = 0$, $i \neq j$. Then

$$ \psi_k^\dagger \mathbf{I} = \frac{1}{N} \int_0^\infty dT \sum_{i=1}^N \sum_{j=1}^N (\psi_k^\dagger)_i J_{ij} \delta(T - \tau_{ij}) X_j(t-T) \quad (A2) $$

and

$$ \begin{aligned} & \sum_{i=1}^{N} (\psi_k^\dagger)_i J_{ij} \delta(T - \tau_{ij}) \\ &= N \frac{1}{N} \sum_{i=1}^{N} (\psi_k^\dagger)_i J_{ij} \delta(T - \tau_{ij}) \\ &= N \left( \frac{1}{N} \sum_{i=1}^{N} (\psi_k^\dagger)_i J_{ij} \right) \left( \frac{1}{N} \sum_{i=1}^{N} \delta(T - \tau_{ij}) \right) + \varepsilon_j \end{aligned} \quad (A3) $$

$$ = \lambda_k(\psi_k^\dagger)_j \left( \int_0^\infty \delta(T - x_j) p_j^\tau(x_j) dx_j + \varepsilon_j \right) + \varepsilon_j \quad (A4) $$

$$ \begin{aligned} &= \lambda_k(\psi_k^\dagger)_j (p_j^\tau(T) + \varepsilon_j^\tau) + \varepsilon_j \\ &\approx \lambda_k(\psi_k^\dagger)_j p_j^\tau(T) . \end{aligned} \quad (A5) $$

The terms $\varepsilon_j$ denote the error of the approximation $E[XY] \approx E[X]E[Y]$ applied in Eq. (A3) and $\varepsilon_j^\tau$ represents the approximation error made by replacing the finite sum by an integral in Eq. (A4) [47]. Both error terms vanish for infinite networks, i.e. $N \to \infty$, yielding Eq. (A5). The term $p_j^\tau(T)$ denotes the probability density function of the conduction delays $\tau_{ij}$ in connections to the node j.

Inserting Eq. (A5) into Eq. (A2) yields

$$ \begin{aligned} \psi_k^\dagger \mathbf{I} &= \frac{\lambda_k}{N} \int_0^\infty dT \sum_{j=1}^N (\psi_k^\dagger)_j p_j^\tau(T) X_j(t-T) \\ &\approx \frac{\lambda_k}{N} \int_0^\infty dT \left( \frac{1}{N} \sum_{j=1}^N (\psi_k^\dagger)_j \right) \times \\ & \qquad \left( \frac{1}{N} \sum_{j=1}^N X_j(t-T) \right) q(T) \\ &= \frac{\lambda_k \bar{\psi}_k}{N} \int_0^\infty dT \bar{X}(t-T) q(T), \end{aligned} \quad (A6) $$
$$ \tau_s \frac{d}{dt} \psi_k^\dagger \mathbf{u}(t) = L[\psi_k^\dagger \mathbf{u}] + \frac{\lambda_k \bar{\psi}_k}{N} \int_0^\infty dT \bar{X}(t-T) q(T) + \sqrt{2D} \psi_k^\dagger \xi(t) $$

$$ \tau_s \frac{d}{dt} U(t) = L[U] + \frac{\lambda_1}{N^2} \int_0^\infty dT \bar{X}(t-T) q(T) + \sqrt{2D} \bar{\xi}(t) $$

$$ \bar{u}(t) = \frac{1}{\Delta t} \int_{t-\Delta t}^{t} U(t') dt' $$

$$ \tau_s \frac{d}{dt} \bar{u}(t) = L[\bar{u}] + \frac{\lambda_1}{N^2} \int_0^\infty dT f[\bar{u}(t-T)] q(T) + \sqrt{2D} \bar{\xi}(t). $$

$$ \begin{aligned} f[x(t)] &= \frac{1}{N} \sum_{j=1}^{N} \int_{t-\Delta t}^{t} \frac{X_j(t')}{\Delta t} dt' \\ &= \frac{1}{N} \sum_{j=1}^{N} M_j[\bar{u}(t)], \end{aligned} $$

<footnote n="595">595 where we assumed homogeneous conduction delays with $p_j^\tau(x) = q(x)$ and we approximated the network mean over $\psi_k$ and $X_j(t-T)$ by the product of their respective mean values $\bar{\psi}_k$ and $\bar{X}(t-T)$. This approximation holds well for large networks.</footnote>

<footnote n="597">597 Multiplying Eq. (A1) by $(\psi_k^\dagger)_i$, taking the corresponding sum and inserting Eq. (A6) leads to</footnote>

<footnote n="600">600 It can be shown [47, 74], that imbalanced infinite Erdoes-Renyi networks can be chosen to have an edge spectrum with a single eigenvalue $\lambda_1$ and the eigenvector $\psi_1 = (1/N, \dots, 1/N)^t$ and a vanishing bulk spectrum. Then for $k=1$</footnote>

<footnote n="603">603 with $U(t) = \psi_1^\dagger \mathbf{u}(t)$. If the synaptic time constant $\tau_s$ is large, i.e. tens of milliseconds, then the equation considers slow synaptic response activity $U(t)$ and fast incoming $\delta$-function-like spike trains (included in $\bar{X}(t)$). Now choosing the slow time scale as the primary scale of description of the network's dynamics, we average the system activity over a short time window $\Delta t$ with</footnote>

<footnote n="608">608 yielding</footnote>

<footnote n="609">609 Here,</footnote>

<footnote n="610">610 where the last equation assumes asynchronous neural activity in the population in line with previous studies [69, 75] and $M_j$ is the number of spikes at neuron j in the time interval $[t - \Delta t; t]$. 611 612</footnote>
Since the delay times $\tau_{ij}$ of number $K = N^2$ are discrete values, we can write

$$q(T) = \frac{1}{K} \sum_{k=1}^{K} \delta(T - \tau_k),$$

which leads to the final equation (2) with $\bar{J} = \lambda_1/K$.

## Appendix B Stationary conditional probability density

We assume two random normal distributed variables $X$ and $Y$, which are correlated

with the correlation function $\sigma_{XY}^2$. Then the conditional probability density function

$p(y|x)$ is also normal distributed with

$$p(y|x) = \frac{1}{\sqrt{2\pi}\sigma_{y|x}^2} e^{(y-\mu_{y|x})^2/2\sigma_{y|x}^2} \quad (B7)$$

with mean $\mu_{y|x}$ and variance $\sigma_{y|x}^2$.

Moreover, we assume that there is some kind of linear relationship between $X$ and $Y$.

This leads to the linear ansatz

$$\mu_{y|x} = \mu_y + \rho \frac{\sigma_y}{\sigma_x} (x - \mu_x) \quad (B8)$$

with the mean value $\mu_x, \mu_y$ and variance $\sigma_x^2, \sigma_y^2$ of $X, Y$ accompanied by the variance

$$\sigma_{y|x}^2 = \sigma_y^2 (1 - \sigma_{XY}^2) \quad (B9)$$

Inserting Eq. (B8) and (B9) into Eq. (B7) yields the conditional probability density function.

In the context of Eqs. (10) and (11), $\mu_x = \mu_y = \mu_s$, $\sigma_x = \sigma_y = \sigma_s$, $Y = \bar{u}$, $X = \bar{u}_\tau$ we gain $\rho_s(\bar{u}|\bar{u}_\tau) = p(y|x)$. The correlation function is

$$\sigma_{XY}^2 = \frac{E[\bar{u}(t)\bar{u}(t - \tau_k)]}{\sigma_s^2}.$$

## References

[1] Buldú, J.M., Garcia-Ojalvo, J., Torrent, M.C.: Delay-induced resonances in an optical system with feedback. Phys. Rev. E 69(4), 046207 (2004) https://doi.org/10.1103/PhysRevE.69.046207

[2] Eurich, C.W., Milton, J.G.: Noise-induced transitions in human postural sway. Phys. Rev. E 54(6), 6681 (1996) https://doi.org/10.1103/PhysRevE.54.6681

[3] Erneux, T.: Applied Delay Differential Equations. Springer, Heidelberg (2009)
635

636

637

[4] Hutt, A., Bestehorn, M., Wennekers, T.: Pattern formation in intracortical neu- ronal fields. Network: Comput. Neural Syst. **14**, 351–368 (2003) https://doi.org/ 10.1088/0954-898X_14_2_310

638

639

640

[5] Hutt, A., Rougier, N.: Activity spread and breathers induced by finite trans- mission speeds in two-dimensional neural fields. Phys. Rev. E **82**, 055701 (2010) https://doi.org/10.1103/PhysRevE.82.055701

641

642

643

644

[6] Liu, C., Wang, J., Yu, H., Deng, B., Tsang, K.M., Chan, W.L., Wong, Y.K.: The effects of time delay on the stochastic resonance in feed-forward-loop neuronal network motifs. Comm. Nonlin. Sci. Num. Sim. **19**(4), 1088–1096 (2014) https://doi.org/10.1016/j.cnsns.2013.08.021

645

646

647

648

[7] Hashemi, M., Hutt, A., Hight, D., Sleigh, J.: Anesthetic action on the transmission delay between cortex and thalamus explains the beta buzz observed under propo- fol anesthesia. PLoS **12**(6), 0179286 (2017) https://doi.org/10.1371/journal.pone. 0179286

649

650

651

[8] Choi, M.Y., Kim, H.J., Kim, D., Hong, H.: Synchronization in a system of globally coupled oscillators with time delay. Phys. Rev. E **61**(4), 371 (2000) https://doi. org/10.1103/PhysRevE.61.371

652

653

654

[9] Dhamala, M., Jirsa, V.K., Ding, M.: Enhancement of neural synchrony by time de- lay. Phys. Rev. Lett. **92**(7), 074104 (2004) https://doi.org/10.1103/PhysRevLett. 92.074104

655

656

657

[10] Sen, M.K., Bag, B.C., Petrosyan, K.G., Hu, C.-K.: Effect of time delay on the onset of synchronization of the stochastic kuramoto model. J. Stat. Mech., 08018 (2010) https://doi.org/10.1088/1742-5468/2010/08/P08018

658

659

660

[11] Ray, S., Sen, M.K., Baura, A., Bag, B.C.: Lifetime of the incoherent state of coupled phase oscillators. Eur. Phys. J. B **85**(11), 306 (2012) https://doi.org/10. 1140/epjb/e2012-30412-8

661

662

663

[12] Deco, G., Jirsa, V., McIntosh, A.R., Sporns, O., Kötter, R.: Key role of coupling, delay, and noise in resting brain fluctuations. Proc. Natl. Acad. Sci. USA **106**(25), 10302–10307 (2009) https://doi.org/10.1073/pnas.0901831106

664

665

666

667

[13] Noori, R., Park, D., Griffiths, J.D., Bells, S., Frankland, P.W., Mabbott, D., Lefebvre, J.: Activity-dependent myelination: A glial mechanism of oscillatory self-organization in large-scale brain networks. Proc. Natl. Acad. Sci. USA **117**(24), 13227–13237 (2020) https://doi.org/10.1073/pnas.1916646117

668

669

670

[14] Roxin, A., Brunel, N., Hansel, D.: Role of delays in shaping spatiotemporal dy- namics of neuronal activity in large networks. Phys. Rev. Lett. **94**, 238103 (2005) https://doi.org/10.1103/PhysRevLett.94.238103
671

672

673

[15] Deco, G., Jirsa, V.K., Robinson, P.A., Breakspear, M., Friston, K.: The dynamic brain: from spiking neurons to neural masses and cortical fields. PLoS Comput. Biol. 4(8), 1000092 (2008) https://doi.org/10.1371/journal.pcbi.1000092

674

675

676

677

[16] Cabral, J., Castaldo, F., Vohryzek, J., Litvak, V., Bick, C., Lambiotte, R., Friston, K., Kringelbach, M.L., Deco, G.: Metastable oscillatory modes emerge from synchronization in the brain spacetime connectome. Commun. Phys. 5(1), 184 (2022) https://doi.org/10.1038/s42005-022-00950-y

678

679

680

681

[17] Salami, M., Itami, C., Tsumoto, T., Kimura, F.: Change of conduction velocity by regional myelination yields constant latency irrespective of distance between thalamus and cortex. Proceed. Natl. Acad. Sci. 100(10), 6174–6179 (2003) https://doi.org/10.1073/pnas.09373801

682

683

684

[18] Swadlow, H.A., Rosene, D.L., Waxman, S.G.: Characteristics of interhemispheric impulse conduction between prelunate gyri of the rhesus monkey. Exp. Brain Res. 33, 455–467 (1978) https://doi.org/10.1007/BF00235567

685

686

[19] Swadlow, H.A., Waxman, S.G.: Axonal conduction delays. Scholarpedia 7(6), 1451 (2012) https://doi.org/10.4249/scholarpedia.1451

687

688

689

[20] Doiron, B., Litwin-Kumar, A., Rosenbaum, R., Ocker, G.K., Josić, K.: The mechanics of state-dependent neural correlations. Nature Neurosci. 19(3), 383–393 (2016) https://doi.org/10.1038/nn.4242

690

691

692

693

[21] Talidou, A., Frankland, P.W., Mabbott, D., Lefebvre, J.: Homeostatic coordination and up-regulation of neural activity by activity-dependent myelination. Nature Comput. Sci. 2(10), 665–676 (2022) https://doi.org/10.1038/s43588-022-00315-z

694

695

696

697

[22] Shadlen, M.N., Newsome, W.T.: The variable discharge of cortical neurons: implications for connectivity, computation, and information coding. J. Neurosci. 18(10), 3870–3896 (1998) https://doi.org/10.1523/JNEUROSCI.18-10-03870.1998

698

699

[23] Izhikevich, E.M.: Polychronization: computation with spikes. Neural Comput. 18(2), 245–282 (2006) https://doi.org/10.1162/089976606775093882

700

701

702

703

704

[24] Gibson, E.M., Purger, D., Mount, C.W., Goldstein, A.K., Lin, G.L., Wood, L.S., Inema, I., Miller, S.E., Bieri, G., Zuchero, J.B., Barres, B.A., Woo, P.J., Vogel, H., Monje, M.: Neuronal activity promotes oligodendrogenesis and adaptive myelination in the mammalian brain. Science 344(6183), 1252304 (2014) https://doi.org/10.1126/science.1252304

705

706

707

[25] Auer, F., Vagionitis, S., Czopka, T.: Evidence for myelin sheath remodeling in the cns revealed by in vivo imaging. Curr. Biol. 28(4), 549–559 (2018) https://doi.org/10.1038/s41593-018-0121-5
708

709

[26] Mount, C.W., Monje, M.: Wrapped to adapt: experience-dependent myelination. Neuron **95**(4), 743–756 (2017) https://doi.org/10.1016/j.neuron.2017.07.009

710

711

712

[27] Fields, R.D.: A new mechanism of nervous system plasticity: activity-dependent myelination. Nature Rev. Neurosci. **16**(12), 756–767 (2015) https://doi.org/10.1038/nrn4023

713

714

715

716

[28] Steadman, P.E., Xia, F., Ahmed, M., Mogle, A.J., Penning, A.R.A., Geraghty, A.C., Steenland, H.W., Monje, M., Josselyn, S.A., Frankland, P.W.: Disruption of oligodendrogenesis impairs memory consolidation in adult mice. Neuron **105**(1), 150–164 (2020) https://doi.org/10.1016/j.neuron.2019.10.013

717

718

[29] Fields, R.D., Bukalo, O.: Myelin makes memories. Nature Neurosci. **23**(4), 469–470 (2020) https://doi.org/10.1038/s41593-020-0606-x

719

720

721

[30] Xin, W., Chan, J.R.: Myelin plasticity: sculpting circuits in learning and memory. Nature Rev. Neurosci. **21**(12), 682–694 (2020) https://doi.org/10.1038/s41583-020-00379-8

722

723

724

[31] Cheng, S.-M., Carr, C.E.: Functional delay of myelination of auditory delay lines in the nucleus laminaris of the barn owl. Develop. Neurobiol. **67**(14), 1957–1974 (2007) https://doi.org/10.1002/dneu.20541

725

726

727

728

[32] Bells, S., Lefebvre, J., Prescott, S.A., Dockstader, C., Bouffet, E., Skocic, J., Laughlin, S., Mabbott, D.J.: Changes in white matter microstructure impact cognition by disrupting the ability of neural assemblies to synchronize. J. Neurosci. **37**(34), 8227–8238 (2017) https://doi.org/10.1523/JNEUROSCI.0560-17.2017

729

730

731

[33] Bells, S., Lefebvre, J., Longoni, G., Narayanan, S., Arnold, D.L., Yeh, E.A., Mabbott, D.J.: White matter plasticity and maturation in human cognition. Glia **67**(11), 2020–2037 (2019) https://doi.org/10.1002/glia.23661

732

733

[34] Nature Neuroscience: Glia move to the foreground. Nature Neuroscience **27**(8), 1427–1427 (2024) https://doi.org/10.1038/s41593-024-01735-y

734

735

736

[35] Ferraina, S., Paré, M., Wurtz, R.H.: Comparison of cortico-cortical and cortico-collicular signals for the generation of saccadic eye movements. Journal of neurophysiology **87**(2), 845–858 (2002) https://doi.org/10.1152/jn.00317.2001

737

738

739

740

[36] Swadlow, H.A.: Efferent neurons and suspected interneurons in s-1 forelimb representation of the awake rabbit: receptive fields and axonal properties. Journal of Neurophysiology **63**(6), 1477–1498 (1990) https://doi.org/10.1152/jn.1990.63.6.1477

741

742

743

[37] Horsthemke, W., Lefever, R.: Noise-Induced Transitions: Theory and Applications in Physics, Chemistry, and Biology. Springer Series in Synergetics, vol. 15. Springer, Heidelberg (1984). https://doi.org/10.1007/978-3-642-69689-3
744

745

746

[38] Frank, T.D., Beek, P.J.: Stationary solutions of linear stochastic delay differential equations: Applications to biological systems. Phys. Rev. E **64**(2), 021917 (2001) https://doi.org/10.1103/PhysRevE.64.021917

747

748

749

[39] Frank, T.D., Beek, P.J., Friedrich, R.: Fokker-Planck perspective on stochastic delay systems: Exact solutions and data analysis of biological systems. Phys. Rev. E **68**(2), 021912 (2003) https://doi.org/10.1103/PhysRevE.68.021912

750

751

752

[40] Frank, T.D.: Delay Fokker-Planck equations, perturbation theory, and data analysis for nonlinear stochastic systems with time delays. Phys. Rev. E **71**(3), 031106 (2005) https://doi.org/10.1103/PhysRevE.71.031106

753

754

755

756

[41] Frank, T.D.: Kramers–Moyal expansion for stochastic differential equations with single and multiple delays: Applications to financial physics and neurophysics. Phys. Lett. A **360**(4-5), 552–562 (2007) https://doi.org/10.1016/j.physleta.2006.08.062

757

758

759

[42] Pajevic, S., Plenz, D., Basser, P.J., Fields, R.D.: Oligodendrocyte-mediated myelin plasticity and its role in neural synchronization. eLife **12**, 81982 (2023) https://doi.org/10.7554/eLife.81982

760

761

762

763

[43] Knowles, J.K., Xu, H., Soane, C., Batra, A., Saucedo, T., Frost, E., Tam, L.T., Fraga, D., Ni, L., Villar, K., Talmi, S., Huguenard, J.R., Monje, M.: Maladaptive myelination promotes generalized epilepsy progression. Nature Neurosci. **25**(5), 596–606 (2022) https://doi.org/10.1038/s41593-022-01052-2

764

765

766

767

768

769

[44] Rahmanzadeh, R., Lu, P.-J., Barakovic, M., Weigel, M., Maggi, P., Nguyen, T.D., Schiavi, S., Daducci, A., La Rosa, F., Schaedelin, S., Absinta, M., Reich, D.S., Sati, P., Wang, Y., Bach Cuadra, M., Radue, E.-W., Kuhle, J., Kappos, L., Granziera, C.: Myelin and axon pathology in multiple sclerosis assessed by myelin water and multi-shell diffusion imaging. Brain **144**(6), 1684–1696 (2021) https://doi.org/10.1093/brain/awab088

770

771

772

773

[45] Lam, S.K., Pitrou, A., Seibert, S.: Numba: a llvm-based python jit compiler. In: Proceedings of the Second Workshop on the LLVM Compiler Infrastructure in HPC. LLVM '15. Association for Computing Machinery, New York, NY, USA (2015). https://doi.org/10.1145/2833157.2833162

774

775

776

[46] Guillouzic, S., L'Heureux, I., Longtin, A.: Small delay approximation of stochastic delay differential equations. Phys. Rev. E **59**(4), 3970 (1999) https://doi.org/10.1103/PhysRevE.59.3970

777

778

779

[47] Lefebvre, J., Hutt, A.: Induced synchronization by endogenous noise modulation in finite-size random neural networks: A stochastic mean-field study. Chaos **33**(12), 123110 (2023) https://doi.org/10.1063/5.0167771

780

[48] Atay, F.M.: Distributed delays facilitate amplitude death of coupled oscillators.
781

782

Physical Rev. Lett. 91(9), 094101 (2003) https://doi.org/10.1103/PhysRevLett.91.094101

783

784

785

[49] Meyer, U., Shao, J., Chakrabarty, S., Brandt, S.F., Luksch, H., Wessel, R.: Distributed delays stabilize neural feedback systems. Biol. Cybern. 99, 79-87 (2008) https://doi.org/10.1007/s00422-008-0239-8

786

787

[50] Tavakoli, S.K., Longtin, A.: Multi-delay complexity collapse. Phys. Rev. Res. 2(3), 033485 (2020) https://doi.org/10.1103/PhysRevResearch.2.033485

788

789

790

[51] Hutt, A., Lefebvre, J., Longtin, A.: Delay stabilizes stochastic systems near a non-oscillatory instability. Europhys. Lett. 98(2), 20004 (2012) https://doi.org/10.1209/0295-5075/98/20004

791

792

[52] Lefebvre, J., Hutt, A.: Additive noise quenches delay-induced oscillations. Europhys. Lett. 102(6), 60003 (2013) https://doi.org/10.1209/0295-5075/102/60003

793

794

795

[53] Naud, R., Longtin, A.: Spike-diffuse-spike model of axonal propagation links demyelination to compound action potential dispersion. J. Math. Neurosci. 9, 3 (2019) https://doi.org/10.1186/s13408-019-0071-6

796

797

798

[54] Hutt, A., Mierau, A., Lefebvre, J.: Dynamic control of synchronous activity in networks of spiking neurons. PloS One 11(9), 0161488 (2016) https://doi.org/10.1371/journal.pone.0161488

799

800

801

[55] Lefebvre, J., Hutt, A., Knebel, J.-F., Whittingstall, K., Murray, M.M.: Stimulus statistics shape oscillations in nonlinear recurrent neural networks. J. Neurosci. 35(7), 2895-2903 (2015) https://doi.org/10.1523/JNEUROSCI.3609-14.2015

802

803

804

805

[56] Hutt, A., Lefebvre, J., Hight, D., Kaiser, H.A.: Phase coherence induced by additive Gaussian and non-Gaussian noise in excitable networks with application to burst suppression-like brain signals. Front. Appl. Math. Stat. 5, 69 (2020) https://doi.org/10.3389/fams.2019.00069

806

807

808

809

[57] Lim, S., McKee, J.L., Woloszyn, L., Amit, Y., Freedman, D.J., Sheinberg, D.L., Brunel, N.: Inferring learning rules from distributions of firing rates in cortical neurons. Nature Neurosci. 18, 1804-1810 (2015) https://doi.org/10.1038/nn.4158

810

811

812

[58] Ly, C., Marsat, G.: Variable synaptic strengths controls the firing rate distribution in feedforward neural networks. J. Comput. Neurosci. 44, 75-95 (2018) https://doi.org/10.1007/s10827-017-0670-8

813

814

815

[59] Williamson, J.M., Lyons, D.A.: Myelin dynamics throughout life: an ever-changing landscape ? Front. Cell. Neurosci. 12, 424 (2018) https://doi.org/10.3389/fncel.2018.00424
[60] Allen, N.J., Lyons, D.A.: Glia as architects of central nervous system formation and function. Science **362**(6411), 181–185 (2018) https://doi.org/10.1126/science.aat0473

[61] Suminaite, D., Lyons, D.A., Livesey, M.R.: Myelinated axon physiology and regulation of neural circuit function. Glia **67**(11), 2050–2062 (2019) https://doi.org/10.1002/glia.23665

[62] Seidl, A.H.: Regulation of conduction time along axons. Neuroscience **276**, 126–134 (2014) https://doi.org/10.1016/j.neuroscience.2013.06.047

[63] Jasper, H.H.: Jasper's Basic Mechanisms of the Epilepsies vol. 80. OUP USA, Bethesda (MD) (2012)

[64] Bacmeister, C.M., Huang, R., Osso, L.A., Thornton, M.A., Conant, L., Chavez, A.R., Poleg-Polsky, A., Hughes, E.G.: Motor learning drives dynamic patterns of intermittent myelination on learning-activated axons. Nature Neurosci. **25**(10), 1300–1313 (2022) https://doi.org/10.1038/s41593-022-01169-4

[65] Pan, S., Mayoral, S.R., Choi, H.S., Chan, J.R., Kheirbek, M.A.: Preservation of a remote fear memory requires new myelin formation. Nature Neurosci. **23**(4), 487–499 (2020) https://doi.org/10.1038/s41593-019-0582-1

[66] Kohama, S.G., Rosene, D.L., Sherman, L.S.: Age-related changes in human and non-human primate white matter: from myelination disturbances to cognitive decline. Age **34**, 1093–1110 (2012) https://doi.org/10.1007/s11357-011-9357-7

[67] Wang, F., Ren, S.-Y., Chen, J.-F., Liu, K., Li, R.-X., Li, Z.-F., Hu, B., Niu, J.-Q., Xiao, L., Chan, J.R., Mei, F.: Myelin degeneration and diminished myelin renewal contribute to age-related deficits in memory. Nature Neurosci. **23**(4), 481–486 (2020) https://doi.org/10.1038/s41593-020-0588-8

[68] Dumont, G., Northoff, G., Longtin, A.: Linear noise approximation for oscillations in a stochastic inhibitory network with delay. Phys. Rev. E **90**, 012702 (2014) https://doi.org/10.1103/PhysRevE.90.012702

[69] Brunel, N., Hakim, V.: Fast global oscillations in networks of integrate-and-fire neurons with low firing rates. Neural Comput. **11**, 1621–1671 (1999) https://doi.org/10.1162/089976699300016179

[70] Doiron, B., Lindner, B., Longtin, A., Maler, L., Bastian, J.: Oscillatory activity in electrosensory neurons increases with the spatial correlation of the stochastic input stimulus. Phys. Rev. Lett. **93**, 048101 (2004) https://doi.org/10.1103/PhysRevLett.93.048101

[71] Yang, S.M., Michel, K., Jokhi, V., Nedivi, E., Arlotta, P.: Neuron class-specific responses govern adaptive myelin remodeling in the neocortex. Science **370**(6523),
[72] Call, C.L., Bergles, D.E.: Cortical neurons exhibit diverse myelination patterns that scale between mouse brain regions and regenerate after demyelination. Nat.

[73] René, A., Longtin, A.: Mean, covariance, and effective dimension of stochastic distributed delay dynamics. Chaos **27**(11) (2017) https://doi.org/10.1063/1.5007866

[74] Hutt, A., Wahl, T., Voges, N., Hausmann, J., Lefebvre, J.: Coherence resonance in random Erdos-Renyi neural networks : mean-field theory. Front. Appl. Math.

[75] Brunel, N.: Dynamics of sparsely connected networks of excitatory and inhibitory spiking neurons. J. Comput. Neurosci. **8**, 183–208 (2000) https://doi.org/10.