Oligodendrocyte-mediated myelin plasticity and its role in neural synchronization

  1. Sinisa Pajevic  Is a corresponding author
  2. Dietmar Plenz
  3. Peter J Basser
  4. R Douglas Fields
  1. Section on Critical Brain Dynamics, National Institute of Mental Health, NIH, United States
  2. Section on Quantitative Imaging and Tissue Sciences, Eunice Kennedy Shriver National Institute of Child Health and Human Development, NIH, United States
  3. Nervous System Development and Plasticity Section, Eunice Kennedy Shriver National Institute of Child Health and Human Development, NIH, United States


Temporal synchrony of signals arriving from different neurons or brain regions is essential for proper neural processing. Nevertheless, it is not well understood how such synchrony is achieved and maintained in a complex network of time-delayed neural interactions. Myelin plasticity, accomplished by oligodendrocytes (OLs), has been suggested as an efficient mechanism for controlling timing in brain communications through adaptive changes of axonal conduction velocity and consequently conduction time delays, or latencies; however, local rules and feedback mechanisms that OLs use to achieve synchronization are not known. We propose a mathematical model of oligodendrocyte-mediated myelin plasticity (OMP) in which OLs play an active role in providing such feedback. This is achieved without using arrival times at the synapse or modulatory signaling from astrocytes; instead, it relies on the presence of global and transient OL responses to local action potentials in the axons they myelinate. While inspired by OL morphology, we provide the theoretical underpinnings that motivated the model and explore its performance for a wide range of its parameters. Our results indicate that when the characteristic time of OL’s transient intracellular responses to neural spikes is between 10 and 40 ms and the firing rates in individual axons are relatively low (10 Hz), the OMP model efficiently synchronizes correlated and time-locked signals while latencies in axons carrying independent signals are unaffected. This suggests a novel form of selective synchronization in the CNS in which oligodendrocytes play an active role by modulating the conduction delays of correlated spike trains as they traverse to their targets.

Editor's evaluation

This paper presents a new mathematical model describing biologically plausible feedback that glial cells might use to properly modify the conduction velocity in axons and promote optimal timing of neural impulses through changes in myelination. This work provides an important step forward by providing the theory for myelin-mediated neuronal plasticity.



Temporal precision required in neural processing can range from sub-millisecond in sound localization and echolocation tasks to milliseconds and hundreds of milliseconds in perceptual and motor system signal processing. Often, this is a consequence of individual neural cells or brain regions requiring a narrow time window to integrate signals arriving from multiple sources. Signals traveling from distant regions commonly traverse complex conduction paths along which conduction velocity (CV) is not constant and undergoes dynamical changes and perturbations, particularly during the development. This will alter the arrival times of action potentials which may undermine the required temporal precision for information processing. It has been argued for more than a decade that a solution to this problem is the adaptive adjustment of the CV through a mechanism of myelin plasticity (MP) Fields, 2005; Fields, 2008, which postulates that myelination is an adaptive and neural activity-dependent process. Modifying myelin sheath thickness and node of Ranvier structure provides the most efficient means to alter conduction delays through changes in CV. There is growing evidence that myelin plasticity is important for fear conditioning Pan et al., 2020; Steadman et al., 2020, spatial learning Wang et al., 2020; Steadman et al., 2020, and is shown to be essential for motor skill learning Bacmeister et al., 2020; Kato et al., 2020; McKenzie et al., 2014; Xiao et al., 2016. Yet, very little progress has been made in understanding the local learning rules in this new form of plasticity and what feedback oligodendrocytes (OLs) use to properly adjust myelination in the CNS. OLs are mostly located far from the target neurons and lack direct feedback on what the desired CV is because the information about the arrival times of the action potentials, that is spikes, is not available at these intermediate locations. Moreover, in most studies of activity-dependent myelination (ADM) Pajevic et al., 2014; Fields, 2015; Dutta et al., 2019; Stevens et al., 2002; Talidou et al., 2021 the precise timing of individual spikes is ignored. It is well known that the introduction of time delays can change both stability and synchronizability at a system level, which then provides an indirect mechanism for MP to affect both the stability Pajevic et al., 2014 and synchronization, for example in a network of spiking neurons Talidou et al., 2021. However, these schemes are based on the activity rate in the connections and do not explicitly include spike timing information in their local rules.

In principle, the arrival times at the target neuron can be explicitly used as the feedback signal, via learning curves similar to that of spike-timing-dependent plasticity (STDP) Bi and Poo, 1998 but with some important differences. In STDP, the crucial parameter is the pre- and post-synaptic spike time difference, Δt, the sign of which determines whether long-term potentiation or depression occurs, with Δt=0 marking the sharp transition between the two. For MP, such discontinuous learning curves would be unstable and hence must be smoothly ramped across the Δt=0 line Eurich et al., 1999; Pajevic et al., 2015. More importantly, any such feedback information at the target will have to be passed in a retrograde fashion, which is problematic since OLs are mostly located very far from the post-synaptic targets of the axons they myelinate. The same problem applies to schemes in which a network of Kuramoto oscillators is studied and the feedback is based on the phase differences Noori et al., 2020.

To develop spike-timing-dependent myelination (STDM) rules, it becomes important to consider schemes in which the mediators of feedback have to act locally and adjust the delays only based on local signal timing information, where the final arrival time error is not available. In this work, we propose models in which OLs are not only the myelinating agents but also serve as the mediators, providing feedback through the interaction with spikes in different axons. We call this form of STDM oligodendrocyte-mediated myelin plasticity (OMP). Specifically, we focus on a particular type of OMP, which uses the transient temporal profile of the OL responses to neural spikes as a reference for adjusting CV and relative time delays. We use theoretical arguments, mathematical modeling, and simulations to show that even the simplest form of OMP models can lead to effective synchronization of correlated and time-locked spikes while leaving temporally uncorrelated spikes unaffected.

Oligodendrocyte-mediated myelin plasticity (OMP)

Our OMP model was inspired by the fact that OL morphology differs drastically from that of Schwann cell (SC), which myelinate axons in the PNS. The most important morphological difference between these two myelinating cell types is that the OL extends many of its processes to myelinate multiple axons, while the SC myelinates only a single axon (Figure 1A and B). We hypothesize that this difference comes from the distinct functional roles myelination plays in the PNS and CNS. In the PNS, the goal is to maximize the CV, while in the CNS the presumed goal is to optimize the synchrony of signals arriving from multiple sources. The OL-axon connectivity is schematically depicted in Figure 1C and quantified via the myelination matrix, . A single OL can have up to 50 such processes extended to axons within a 100–200 μm distance from the soma, but with a tendency to maximize the number of axons it can myelinate, making it unlikely that a given OL would myelinate the same axon twice Dumas et al., 2015; Walsh et al., 2016. We postulate that this morphology enables a single OL to integrate and compare signals from different axons and act as a mediator in providing the needed feedback for adjusting the relative timing of different signals via dynamic regulation of the CV along these axons. In Figure 1D, we outline a general form for such spike-response OMP models, consisting of three essential steps, described below. In these continuous-time models, the transient nature of the temporal profile of the OL response to an action potential plays a crucial role in creating the needed reference for adjusting the relative delays on different axons.

Figure 1 with 4 supplements see all
OMP motivation and model description.

(A) Cross-section of myelinated axon bundle in CNS illustrates that myelin thickness on adjacent axons differs. (B) Single OL (green) myelinates many different axons (purple), which is in stark contrast to Schwann cells in the PNS, which myelinate only one axon. (C) OL-axon connectivity: OLs tend to avoid placing multiple processes on a single axon, thus maximizing the number of axons it myelinates. (D) Schematic depiction of the OMP-1 model which contains only the three basic steps required for this type of OMP to work. (E) Spike-response curves R(t) for increasing values of τG. (F) The equation governing the release of G(t) is linear and the response to multiple spikes (vertical lines), is the linear sum of individual responses (dashed lines). The release of M after each spike at any given OL process/axon will be proportional to the amplitude of G(t) at the time of spike (red dot for red spike). (G) The sum over responses can also be viewed as sampling of R(t) in reverse time. Panels A and B have been adapted from Chapter 45 Figure 1 in Fields, 2013, used with permission.

The first essential element of OMP is the release of a global intracellular signaling factor, G, after each neural spike on any of the axons myelinated by a given OL. This response has a characteristic transient temporal profile, R(t) (Step 1 in Figure 1D). When the OL encounters a sequence of spike trains, the resulting global intracellular signal, G(t), will simply be the sum of the individual responses and will fluctuate in time, thus providing a common and time-dependent reference to all of its processes. To allow for differential myelination between different axons, a local myelin-promoting factor, M, is also required. It is released after each neuronal spike (Step 2) on a given axon a, but its release is catalyzed by G and hence is proportional to the global signal, G(t). The OL process that myelinates axon a has a dynamically changing concentration of such local factor, Ma(t), that depends on the temporal profile of the spike trains and will generally differ between the axons myelinated by a given OL. This difference in local concentration of M allows for selective modification of the CV and the conduction delays, τa. One can envision more elaborate OMP models (e.g. OMP-n in Figure 1—figure supplement 1), in which the release of M can depend on several different global factors, GM, which potentially are released via a cascade of events triggered by the original factor G. Here, we use the simplest model of this kind, OMP-1, which has only one global signal, G, that does both, responding to neuronal spikes and modulating the local release of Ma at axon a, and for the remainder of this work, we are going to refer to the OMP-1 model simply as the OMP model. The last step (Step 3) represents two continuous processes, one being the conversion of M into myelin with some addition rate, λA, and the second being the steady removal of myelin with rate λR Dutta et al., 2019, resulting in time-varying conduction delays, τa(t). The dynamics of the OMP model and its main variables G(t), Ma(t), and τa(t), is governed by a set of equations that implement steps 1–3 (Equation 8, Equation 11, and Equation 12), but also include presumed homeostatic regulation of the myelin conversion/removal rates (Equation 17), which is needed for the long-time stability of the model. An example of their time progression is shown in Figure 1—figure supplement 2. These equations govern the behavior of a population of OL at a particular segment located at distance xo along the axonal bundle (see Figure 1G). The full OMP model simulates the net delays across a discrete number of such segments, NO, each containing its own set of equations. We call this sequence of OL segments the oligo-chain (OC), which is graphically depicted in Figure 2A. Based on theoretical arguments elaborated in the next subsection, we expect OMP to synchronize correlated/time-locked signals on different axons while leaving the axons carrying independent signals unaltered (Figure 2B), even in situations where individual OL myelinates axons carrying both types of signals (Figure 2C). We study its ability to synchronize signals arriving from multiple sources that are temporally dispersed by fixed delays, representing persistent relative temporal delays among NA axons that arise either from developmental and other structural disturbances in the conduction pathways, or are due to fixed temporal sequences of activations between different sources. We expect those fixed delays, td(a),a=1,,NA, to be compensated via corresponding adaptive delays of the myelinated axonal segments, τa,a=1,,NA, to facilitate synchronous arrival, as depicted graphically in Figure 2D. We test this ability of the OMP model in our simulations that are, together with the details of the model, described in Materials and Methods.

Figure 2 with 2 supplements see all
OMP simulations with action potentials/spikes conducted along ‘oligodendrocyte’ chain (OC).

(A) OC with NO ‘effective oligodendrocytes’ (segments), myelinating NA axons ( is an NO×NA matrix of ones). (B) examples of ‘pure’ spikes used in our simulations: correlated versus independent spikes, and Poisson versus regular spiking. Gray lines are moving averages of all spikes using a sliding Gaussian window with 10 ms RMS width. (C) example of ‘mixed’ signals in which axons in Group 1 are conducting independent spikes and those in Group 2 carry correlated spikes. The two groups can potentially interfere with the expected synchronization behavior of each ‘pure’ group. Different groups could also contain spikes time-locked within each group, but independent between the groups. (D) schematic depiction of time delays between the signal source and the target. Fixed delays are not modifiable and essentially represent the spread in spike times as they enter the axonal bundle of myelinated axons, whose delays are adaptive. Note, that the horizontal bars represent the magnitude of the fixed and adaptive delays, not the axons.

OMP theory

While OL morphology inspired our OMP model, the motivation was also guided by basic theoretical considerations which provide hints about the model’s synchronization performance. We presume that the spike trains in individual axons have the inter-spike intervals (ISIs) that are independent and identically distributed (i.i.d.) with density, pISI(t(i)), that is, are generated by a renewal process. We denote the mean spiking/firing rate as fs and the mean inter-spike time interval with τs=1/fs. The expression for Ga(t), which is a contribution to G(t) coming from spikes on axon a, can then be estimated by recognizing that the sum over spikes in Equation 10 can equivalently be seen as sampling the response function, R(t), in reverse time (see Figure 1G). Here, the important parameter is the time to the last spike, δt, as the remaining spikes are just the sequential samples drawn from pISI(t(i)), effectively making Ga(t) a function of δta=min|ttka|,tka<t, where tka indicates the spike times on axon a for a given OL. If we label the cumulative sum of the subsequent interspike intervals as tkc=ikt(i), the expression for Ga(t) can be written as

(1) Ga(δta)=R(δta)+k=1R(δta+tkc).

For uncorrelated signals, due to symmetry, any of the axons is equally likely to produce a spike, resulting in equidistributed concentrations of M guided only by the average concentration, Gav=G(t)t. Predicting the synchronization effects for time-locked signals when the OL myelinates many axons is a more difficult task, particularly when Gaussian jitter with spread σj is added to the specified renewal process ISIs and also due to non-linear saturation effects in the learning equation (see Materials and methods).

When combining Ga from all axons, mutual ordering of spikes will generally need to be considered. In most cases, this derivation will depend on the exact form of R(t) and pISI(t(i)). For example, for the case of two axons (NA=2) and a given fixed delay, td, we need to consider separately the cases where δttb=τstd and where δt>tb (see Figure 1—figure supplement 3A), in order to derive the expressions for Gav, M1, M2, as a function of τs,τG, and td. For τs<2td, the "leading edge" axon becomes the follower, in which case the myelination pattern is reversed, making its CV faster instead of slower (Figure 1—figure supplement 3B, C). We demonstrate such calculation in the Appendix B for a simple case with NA=2, regular spiking, and without jitter. The calculations become more cumbersome with increasing NA, even when jitter is ignored, and this is true for most of the renewal processes governing the spiking dynamics on a given axon. However, in the case of a pure Poisson process (pISI(t)=et/τs), we can utilize its memoryless property to derive simple expressions for Gav and Ma. The expression for Gav, at any of the locations along the axons, is simply

(2) Gav=G(t)t=NA/τs.

The expectation for M in the case of Poisson spiking is also simple. For example, when NA=2, the average concentration of M at the leading edge axon is M1(t)=CMGav while for the lagging axon M2(t)=CM(Gav+R(td)), where CM=λM/(λAτs). Hence, the ratio of their concentrations is always lower than one,

(3) r2=M1M2=GavGav+R(td)=22+τsR(td)<1,

which will consistently lead to a greater increase in CV of the lagging axon, hence supporting synchronization. This also reveals the landscape of OMP synchronization, expected to be most efficient for r21, which happens at low firing rates (large τs) and a fast spike response time, τG (Figure 1—figure supplement 3C). In the case of multiple axons, the above argument retains its simplicity and we can write the expected time-averaged concentration of M on axon a as

(4) Mat=CM(Gav+td(i)>td(a)R(td(i)td(a))).

The Equation 4 is an important result that indicates that with Poisson spiking and fixed time delays along different axons, the differential expression of the myelination factor will always myelinate the lagging axons more than the preceding ones. The expression in Equation 4 assumes no jitter but still agrees well with the simulated values for σj<3 ms (Figure 1—figure supplement 4).


Our measure for quantifying the synchronization properties of the OMP model is the spread in the spike arrival times across all axons, στ=SDa(Da+τa), where Da stands for pre-specified fixed delays (normalized to some prescribed value, σD), and τa are the adaptive delays (see Figure 2D). The spread starts with some large value, mostly due to the spread among fixed delays, στ(0)σD=SDa(Da), and in the course of time, due to the OMP dynamics, is reduced to lower values, ideally to zero, indicating perfect synchronization. We collect the time course of στ during learning, and we call it a synchronization profile, στ(t). We study the performance of the OMP model by characterizing the synchronization profiles, στ(t), obtained for a wide range of the OMP model parameter values, utilizing a grid-search-like exploration. The sets of parameter values explored in these simulations are provided in Appendix 1, while the summary description of the OMP model parameters is given in Table 1. In most studies, we obtained a large number of synchronization profiles, στ(t), one for each parameter setting/set, which were then reported as averages over all runs, but grouped by a given parameter of interest, or characterized using a model fitting and selection procedure described in Materials and methods (see also Figure 3—figure supplement 1C). In particular, we focus on the long-time baseline parameter, στ=limtστ(t) which represents the long-time ability of the OMP model to synchronize signals for a given parameter setting.

Table 1
List of symbols used in the manuscript: OMP model parameters, OMP variables, spiking signal parameters, and data analysis and quantification parameters and other symbols.

For each we provide a short description and the range of the parameter values, or its dimensionality, in the case of the variables.

List of symbols used in this manuscript
SymbolDescriptionExplored Values/Dimensionality
OMP model parameters
R(t)OL transient response curvesee τG
τGcharacteristic time for R(t) (τG=τr=τd)[2-100] ms
λMproduction rate for factor M[0.01, 0.5] ms-1
λAmyelin conversion/addition rate[0.01, 0.5] ms-1
λRmyelin removal rateEquation 16; variable forλH>0
λHhomeostatic rate[0, 10-2] ms-2
NOnumber of OL in OC[1-20]
NAnumber of axons[2-150]
Mmyelination matrixNO×NA: {full connectivity}
τminminimal delay attainable on axons3 ms
τmaxmaximal delay attainable on axons100 ms
τnomnominal/homeostatic delay[10-90] ms
τ0initial adaptive delays parameterτ0=τnom
pτpercent spread of initial adaptive delays5%
OMP model variables
G(t)global OL signal (Equation 8)NO variables
Ma(t)concentration of M on axon a (Equation 11)NO×NA  variables
τa(t)conduction delay(s) for axon a (Equation 12)NA or NO×NA variables
τao(t)τa for OC segment o (Equation 12)NO×NA variables
λR(t)myelin removal rate (Equation 17)OMP parameter when λH=0
Spiking signal parameters
τsinter-spike time interval parameter[10-250] ms
fsmean firing rate[4-100] Hz
tRrefractory period in refractory Poisson process[0–100] ms
σjamount of jitter given to spikes[0–10] ms
σDSD for fixed delays[0–20 ] ms
σspercent variability in firing rates[0–20] %
Texptotal duration of simulation[10 min - 10 hrs ]
nenumber of recorded epochs during simulation[20–500 ]
Teduration of each recorded epochsTexp/ne [10 sec - 5 min]
nrnumber of replicated simulations/trials[3 - 10]
Quantification parameters/other symbols
στOC spread SD (Da+τa)evaluated after each epoch
στ(t)synchronization profile during OMP learningστ values for ne epochs
στ(0)initial spread before learning, στ(0)=στ(0)determined by σD, τ0 , and pτ
στlong-time baseline, στ=limtστ(t)estimated via model fitting
τLcharacteristic time for synchronizationestimated via model fitting
Lτlearning/synchronization rate = 1/τL
tdgeneric name for fixed delaysNA
Dapre-specified/normalized fixed delay on axon arandom values, N(0,σD)
t(i)ISI for the ith intervalNA
pISI(t(i))probability density function of ISINA

We show in Figure 3A the estimated distribution of στ, as well as the model selection chart, when fitting a large number (n=17280) of synchronization profiles, στ(t), obtained using a wide range of OMP parameters (Appendix 1—table 1, Figure 3—figure supplement 2, and Figure 3—figure supplement 1 for details). The results indicate that the OMP model behaves as desired. When the myelinated axons conducted correlated spikes, we observed a significant reduction in the conduction delay spread, στ, that is, a significant increase in synchronization. No change was observed for independent spikes resulting in synchronization profiles best fit to the constant model, C (purple) (see Model Fitting for synchronization profiles, Figure 3—figure supplement 2). For correlated spikes, a single exponential approach to synchrony was most commonly observed (E1). In several instances, στ(t), was not monotonic and sometimes appeared oscillatory, which can be seen upon inspecting individual runs, particularly for small NO (NO<3), short τG, and small jitter, σj (Figure 3B). As suspected, the OMP model with only a single OL has an inherent instability due to the presence of fixed delays, which could be alleviated by increasing τG and σj. This instability, however, rapidly disappeared for any level of jitter when longer OCs are used (Figure 3B, solid lines; see also Figure 1—figure supplement 1 and Figure 4—figure supplement 3). These trends can be seen in Figure 3C, where the averages of all στ(t) profiles grouped by NO are shown. The runs for small NO were longer, to match them in terms of the number of spikes processed. When plotted in actual time, it becomes evident that having more OLs in an OC greatly increases OC stability as well as the synchronization/learning rate, Lτ. Since the averages were obtained over a large number of runs, the standard error (SE) is small and the oscillatory or other forms of instabilities average out. To better quantify all στ(t) s we fit them to five different models, as described in Methods. In Figure 3D, we show the στ distributions and the model selection chart for different values of NO. They all indicate that, with increasing NO, synchrony is improved and instabilities disappear. For NO=10 most of the OMP simulations yielded a stable exponential decay to a synchronized state (73% of all runs reduced the arrival time spread from the initial στ(0)=10 ms to below 3 ms; for τG=10 ms this increases to 97%, with 74% synchronized below 1 ms).

Figure 3 with 2 supplements see all
OMP model behavior as a function of OMP and signal parameters.

(A) When spikes are independent, στ(t) consistently show no change in overall synchronization (purple; innermost circle in the model selection pie-chart); for correlated signals significant synchronization occurs, largely depending on τs and τG. (B) individual OMP synchronization profiles, στ(t), simulations for NA=10, λH=10-6, λM=0.02, with τG=30 ms (top panel) vs. τG=10 ms (bottom). Dashed lines are for NO=1 and solid lines for NO=10. Colors indicate the jitter level: σj=1 ms (red), σj=3 ms (blue), σj=5 ms (green). (C) dependence of στ(t) on the number of OLs, NO, in the OC. We show both the comparison based on the raw time, as well as when matched in terms of the total number of neural spikes encountered by the OLs (inset). (D) density of στ estimates dependence of στ on the number of OL in OC. (E) Percent reduction in στ vs τG and τs for σD=10 ms; (inset) ratio of the M concentrations for two axon case, NA=2, which matches well the overall pattern of synchronization. (F) dependence of στ (bar plot) and the synchronization time constant, τL=1/Lτ, (semi-log plot, above) on the number of axons, NA values are only for runs declared as E1 blue curve inset indicates percent of E1, for different NA (for all spike rates).

For correlated spikes, the effectiveness of synchronization strongly depended on two temporal parameters: the characteristic response time of the OL, τG, and the mean inter-spike interval, τs, which was assumed to be the same for all axons. In Figure 3E, we explored the synchronization effect as a function of τG and τs; the results shown are for NO=5, σj=1 (see Appendix 1—table 1). For very short τG<10 ms, performance was unstable, reflecting the fact that for the short-lasting spike responses, R(t), the comparison window between spikes in different axons is too narrow. On the other hand, having spike responses last too long, i.e., τG>40 ms, makes the resulting G(t) too smooth to differentially release M, particularly for small τs (Figure 1—figure supplement 3B, C). Accordingly, we found that for intermediate values of τG[10-40] ms, synchronization was predictably achieved and was highly efficient for firing rates fs<10 Hz.

In Figure 3F, we study the dependence of στ and the synchronization time constant, τL=1/Lτ, on the number of axons that OLs myelinate, NA. The results are grouped by the firing rate simulated. For a low firing rate of 5 Hz, optimal performance was achieved for NA=20, whereas synchronization became progressively harder with increase in firing rate and number of axons to synchronize. In general, increasing NA sped up synchronization for all firing rates explored. The fastest synchronization achieved at 5 Hz for an intermediate number of axons (30<NA<60) also coincided with the highest fraction of stable, E1, synchronization profiles (peaking at NA=50, light-blue). We note that myelinating more than NA=50 axons did not improve synchronization efficiency for any of the spiking rates, which, perhaps, hints to why OLs rarely extend more than 50 processes.

In Figure 4A-C, we tested the ability of the OMP model to selectively handle mixed sources of signals by having non-overlapping groups of axons carry spikes with different pISI, or different mutual correlations. Spikes from different groups are inducing responses in the same OL, and hence mutually interfere and can potentially corrupt the expected behavior for the equivalent ‘pure’ group, that is, disrupt OMP’s ability to synchronize the correlated group of signals, or erroneously synchronize the independent signals. In Figure 4A and B, we show the results when one group of axons carries correlated and another carries independent spikes. In Figure 4C, we explore the effect of different groups carrying signals that are correlated within each group, but not between groups. We evaluate στ within each group and compare it with the equivalent pure signal sources, the number of axons being matched in all comparisons. We found an increase in synchronization only for correlated groups, with only slight ‘jamming’ of the synchronization performance due to the presence of other groups with potentially corrupting signals. The independent signals remain unaffected, which is a desired behavior for the selective synchronization of spikes from different neuronal populations.

Figure 4 with 3 supplements see all
Selective synchronization with mixed signals; effects of λH and σD parameters.

(A) Average στ(t) (grouped by σj; shaded regions indicate SE over the parameter sets for NA=10 and σD=5; see Appendix 1—table 2). In the left panel is the average for the set of NA=10 axons conducting only ‘pure’ correlated signals, and on the right are the averages within two groups of 10 axons, out of NA=20 total that OL myelinates, carrying ‘mixed’ signals – one group conducting correlated signals and the other independent signals. The behavior of the correlated groups is clearly distinguishable from the independent ones (indicated by arrows). (B) Estimated density for στ comparing pure (left) vs mixed signals (right panel) for a wider range of parameters, including NA=20 and NA=50, but with the number of axons carrying correlated signals matched in all cases (NA=10, NA=25). (C) Top panel shows στ(t) for two groups of correlated signals (colored lines) with Poisson spiking (correlated within group, but mutually independent); bottom panel shows results for four separately correlated groups + 1 independent group (cyan). Black lines represent the spread over all axons. (D) Investigating the influence of homeostatic regulation and homeostatic rate, λH, grouped by different firing rates and different σD. (E) Model selection chart depicting the relative proportions of oscillatory (E2C2, E2C) vs stable, single exponential (E1) στ(t), for two different values of the F-test fudge-factor pMSE; the innermost circle is for λH=10-2 and the outermost for λH=10-8. (F) OMP synchronization performance for different magnitudes of fixed delays, σD.

For OMP to be operational requires an overall balance between myelin removal (controlled by λR) and myelin addition (controlled by λM and λA). For example, if λR is too large the long-term behavior of the OMP model would lead to complete myelin removal. Here, we achieve such an operational regime by treating λR as a variable and control its rate of change with the homeostatic rate, λH, as defined in the homeostatic equation (Equation 17). The results in Figure 4—figure supplement 2C indicate that this homeostatic regulation, while keeping the model operational, does not play an essential role in OMP synchronization. Although it appears that for low firing rate and for large fixed delays (σD=10 ms), increasing the value of λH has some detrimental effect on model performance, this effect is not seen for σD=5 ms. The oscillatory dynamics of λR depends on λH and λM (see Figure 4—figure supplement 3) which can propagate into synchronization profiles. Nevertheless, the pie-chart in Figure 4E does not indicate that λH strongly influences the model selection, independent of what value of pMSE was used. Figure 2—figure supplement 1, Figure 4—figure supplement 3 show that the oscillations roughly average to the same value and thus do not significantly affect synchronizability in the long term. We note that the average values of λR deviate from the expected value given by Equation 16, even for the case of independent spikes (Figure 2—figure supplement 2). This difference arises from ignoring saturation effects and depends on the values of τnom, τmin, and τmax. We also explored scenarios in which we use the ‘true’ balancing homeostatic value of λR, obtained via a trial run, and set λH=0. The corresponding synchronization profiles were not significantly affected, further indicating that the homeostatic process is not an essential element for achieving spike synchronization in the OMP model.

The production rate of M, λM, while greatly influencing the learning rate, Lτ, has negligible influence on στ (Figure 4—figure supplement 2A and B). Similarly, the conversion rate λA had very little influence on the outcome of OMP (Figure 4—figure supplement 2D), indicating that the simplified OMP model with instantaneous myelination could be a more efficient way of studying its behavior (Figure 3—figure supplement 1E). In Figure 4F, we explored the dependence of στ on the magnitude of fixed delays, σD. For correlated signals, synchronization always occurs but its efficiency decreases in terms of στ/στ(0) when σD becomes large (for σD>10 ms see Figure 4—figure supplement 1B). Such large delays might be commensurate with the delay corrections needed during the development but are presumably much larger than the timing corrections needed in the adult brain.


Here, we report the development of a simple, biologically plausible model of oligodendrocyte-mediated myelin plasticity that synchronizes temporally correlated neuronal spikes, as they travel along an axonal bundle, while leaving temporally independent spikes unaffected. This enables the OLs in the model to counteract the temporal dispersion arising from heterogeneous conduction delays and facilitate synchronous arrival times of spikes coming from distant neuronal populations with correlated activity. The general idea of myelin plasticity is not new Fields, 2005; Fields, 2008, however, a local STDM mechanism by which the brain could robustly and selectively adjust axonal latencies has been missing. Our OMP model introduces robust local learning rules and feedback mechanisms for adaptive changes that yield desired results under a wide range of biologically realistic parameters. The adaptive changes to axonal delays are selectively applied to groups of axons that carry correlated spikes so that they arrive at their targets simultaneously, while those carrying independent or uncorrelated spikes are not affected (Figure 4A–C). This selectivity conforms to known relationship between circuit anatomy and function, in which neighboring axons share similar temporal firing and functional properties; for example, the tonotopic organization of auditory cortex and the cortical homunculus in somatosensory cortex. Such correlated firing also drives refinement of connections between the retina and LGN during development Meister et al., 1991. Myelination usually begins after axons reach their target and become functional, starting from the cell body and proceeding toward the axon terminal. This is clearly evident in the optic nerve, where OL progenitor cells migrate out of the brain and into the optic nerve during development, yet axonal myelination proceeds in the reverse direction, beginning at the retina on retinal ganglion cells and proceeding toward the optic chiasm Ishibashi et al., 2009. Such a proximal-to-distal gradient agrees well with our OMP model as the OLs at the source end of the OC will experience less synchronized signals. Previous work that use phase- and time-dependent models of myelin plasticity Noori et al., 2020; Pajevic et al., 2015 presume a priori that the temporal difference feedback at the target is available to OLs, however, this local information at the target, that is, synapse, will have to be transported in a retrograde fashion. Besides being slow, such a process would suggest that more myelin will be found close to the target rather than far away from it. We note, that with the existence of targeted fast axonal transport, for example via mitochondria, it is still possible that any myelin-promoting factor is transported in a retrograde fashion far from the synapse.

The fact that OLs with more than 50 processes are uncommon is also in accord with our observation that having NA>50 is not advantageous, according to our OMP results. Another testable prediction of this model would be that structural deformation or lesion to a portion of axons carrying correlated activity would lead to re-myelination downstream or immediately after the site of the lesion, while myelin upstream will not be affected. Results in Figure 4 suggest that the efficiency of synchronization is slightly reduced when a fraction of axons carry uncorrelated signals. Hence, we postulate that in the brain, OLs might be removing their processes from such axons while keeping or growing new processes only on those axons that carry synchronized signals. OMP predicts that synchronization is most effective at lower firing rates, and thus we expect it to occur when the brain operates in low spiking rate regimes.

It was observed that OLs can undergo sudden depolarization which leads to significant changes (10%) in Yamazaki et al., 2007, which might suggest that such events are essential for efficient MP. While this can be incorporated into a general OMP scheme (Figure 1D), the OMP-1 shows that passive responses of OL cells are sufficient for synchronizing neural signals, as well that astrocytes, commonly considered the actors in providing the feedback in MP, are not needed. Our results also demonstrate that the presence of noise acts as a stabilizer, both in terms of jitter removing the instabilities in the OMP model as well as Poisson spike dynamics always biasing larger concentrations of the M to be on the axons with larger fixed delays, as opposed to regular spiking which can reverse this pattern and destabilize the system. Hence, Poisson spiking, even though slower and in many situations less efficient than regular spiking, is more reliable for the MP. Our model of synchronization provides further support for the narrative that the noisy brain is a healthy brain, and too much regularity/synchronization in the brain dynamics can lead to its failure.

The OMP model described here, while simple, is not the simplest form that can serve as a proof of concept. An instantaneous myelination model described in Methods, also shows robust synchronization performance (Figure 3—figure supplement 1E), with shortcomings coming only at high learning rates (large λM) and when using homeostatic regulation in Equation 17, making it then a stochastic equation. But, we kept a more general form, since we expect that in future developments a more sophisticated regulation of λM, λA, and λR will be needed, to account for elaborate time-locking patterns where the firing frequencies are different but matched via integer multiples. One of the weaknesses of the current model is that λM is constant and the same for all axons, making the model sensitive to rate differences, which can override the synchronization effects. To make it work, the production rates, λM, need to have their own homeostatic regulation, so that the axons with higher firing rate will down-regulate their λM. In more elaborate OMP models a mix of activity- and time-dependent MP might be needed. The net result will depend on the relative strengths/learning rates between the time-dependent and the activity-dependent learning, something that will require its own independent study. For the current model, we explored its sensitivity to firing rate variability (Figure 4—figure supplement 2E and F) which showed that variations in firing rate greater than 5% are highly detrimental to OMP’s ability to adjust fixed delays properly. When the firing rates across axons are very different, it might also not be desirable to synchronize those groups of axons, nor is easy to define the temporal synchrony in such a situation.

There are other aspects of OMP that are not explored here. For example, the effects of inhomogeneous OL-axon connectivity are not addressed, as we use the same myelination matrix along the OC. The stochasticity in our models comes mainly from the stochasticity of the spikes and their jitter. Future work will address complex patterns of connectivity, other independent sources of noise for both global and local factors, as well as more sophisticated homeostatic regulation discussed above. Delays in factor G are mainly implemented here through the rise time, τr=τG, however, an increase in G coming from a particular OL process will not affect all processes simultaneously, and the relative delays between different processes can have significant effects on the resulting synchronization, which will need to be explored via delay-differential equations, or using a discrete implementation. Some of the proposed mechanisms for myelin plasticity are discrete in nature, for example, the treadmilling model Dutta et al., 2019, but the discreteness is only in the state of myelination; the feedback mechanisms needed for temporal adjustments in any time-dependent MP might still need the continuous-time comparisons. Discretizing time in the simulations can introduce its own effects, which could dominate synchronization performance since MP mechanisms usually require long-time simulations (the MP time scale is many orders of magnitude longer than the neuronal spiking time scale).

In summary, we demonstrated that the simple and biologically plausible adaptive dynamics of the OMP model leads to efficient and selective synchronization of correlated and time-locked signals, without affecting mutually independent streams. This is a novel perspective in brain organization in which white matter conduction properties and myelin plasticity act as a temporal ‘lens’ to ‘focus’ multiple spike trains as they target a particular brain region. Our model also addresses lingering questions about the illusive local learning rules and feedback mechanisms in MP, as it circumvents the need for direct information about the actual arrival times at the target. With its precisely spelled-out dynamics and learning rules, it serves as a useful basis for designing future experimental tests aiming to elucidate the nature of MP in the CNS.

Ideas and speculation

The main goal of our OMP model is to conceptualize a general but biologically plausible myelin plasticity mechanism by which synchronization can be achieved. In the Discussion, we made suggestions about some testable predictions of our model and here we additionally speculate on the biological aspects of our mathematical model, in particular, we discuss the potential candidates for the factor M and the global signal G. Our usage of the terms somewhat implies that M is a molecular factor, while G is some fast propagating signal, for example, intracellular potential or ionic concentration, which is mainly based on the timing constraints implied by our model. It suggests that the release and clearance of a global intracellular signal cannot be too slow since, according to the OMP model, synchronization is not very efficient when the characteristic time for the release and clearance of G, τG, is larger than 80 ms, even for slow firing rates. But it also suggests that the clearance does not have to be too rapid; in fact, having τG too short can be detrimental in some situations. For firing rates below 10 Hz, τG should ideally be in the range 10 ms <τG< 40 ms for effective synchronization of correlated inputs even for large fixed delays (e.g., σD=10 ms). These timing requirements make intracellular Ca2+ a good candidate for the role of G, since it is also a catalyst for many bio-molecular reactions. While at present we only speculate that this is the case, we also emphasize here some of the established biological evidence regarding signaling mechanisms between axons and OL, as well as highlight some of the difficulties in conducting experimental tests of our model.

OLs express many of the same neurotransmitter receptors and ion channels that are expressed by neurons, enabling robust activity-dependent axon-OL communication through several signaling mechanisms that differ depending on the developmental stage of the cells. Several types of neurotransmitter receptors have been identified on the axon underlying compact myelin Stys, 2011; however, the detection of local signaling events between axons and mature OLs with compact myelin is difficult with current methods. Depolarization at the distal tips of the OL processes that are in contact with axons is not accessible for measurement by patch electrode recording at the cell body because of the electrotonic decay over the long slender cell process. Calcium imaging using genetically encoded reporters or fluorescent dyes is limited by the slow kinetics of the indicators. Signaling with axons beneath the compacted layers of the myelin sheath is inaccessible by electrophysiological and live-cell imaging methods, which are obscured by the thick layers of the compacted myelin membrane. Mathematical modeling of the kinetics of local and global signaling in activity-dependent myelin plasticity is hence an important tool for guiding the determination of the types of inter- and intracellular signaling molecules involved in the mechanisms of myelin plasticity.

Confocal imaging of local calcium responses in myelinating oligodendrocytes in cell culture, together with the imaging of local translation of myelin basic protein, show that action potentials in axons cause local calcium transients in OL via vesicular release of glutamate from axons acting on NMDA and glutamate receptors (mGluR) on the oligodendrocyte processes Wake et al., 2011. This promotes the formation of an axo-glial signaling complex, triggering the local translation of myelin basic protein to initiate myelination Wake et al., 2011, preferentially on the electrically active axon Wake et al., 2015. This local signaling can trigger myelin synthesis rapidly, within minutes Wake et al., 2011. The latency of local calcium signaling in OL was slower than the 80 ms required by our model but faster than 500 ms (Figure S7 in Wake et al., 2011). This reflects the slow kinetics of GCaMP2 calcium indicator and the actual signaling kinetics is likely much faster, not occurring via synaptic vesicles but rather by the action-potential-induced exocytosis on glutamatergic vesicles at axon varicosities Wake et al., 2015.

Similar responses are observed using in vivo imaging in zebrafish Hines et al., 2015; Mensch et al., 2015; Krasnow et al., 2018 which further indicate that myelin sheath elongation during development is regulated by the kinetics of calcium transients in oligodendrocytes that are evoked by neuronal activity. In Krasnow et al., 2018 they show that local calcium transients in oligodendrocyte cell processes are independent from one another and that the local calcium signals can be integrated within the cell to trigger global calcium responses in the cell body via temporal summation (Supp Figure 1 in Krasnow et al., 2018). The same study shows that myelin sheath elongation is promoted by high-frequency calcium transients, and sheath shortening is associated with low-frequency calcium transients. It also shows that the elongation occurs approximately one hour after Ca2+ while sheath shortening happens on a much longer time scale. Note the similar asymmetry in the OMP model between the myelin addition and removal rates (λA and λR) but with the important difference that the myelin removal in our model does not depend on activity directly and is controlled only homeostatically. This further emphasizes the need to develop an activity-dependent mechanism of myelin homeostasis in future OMP models.

In spite of accumulating evidence, our suggestion that the role of G might be played by calcium is only speculative and requires further investigation. There are many mechanisms of axon-OL interactions and a recent review of those can be found in Munyeshyaka and Fields, 2022. The OLs express a wide range of calcium channels that can regulate OL formation and function, and the diverse roles they play have already been investigated Paez and Lyons, 2020; Wake et al., 2011; Wake et al., 2015; Fields, 2015; Krasnow et al., 2018. For example, in Wake et al., 2011 it was shown that the local calcium transients, released via glutamatergic vesicles and in response to action potentials firing in OL processes, can be blocked by botulinum toxin, but a global somatic calcium response persists, due to purinergic receptors that are expressed throughout the OL cell membrane. It is also important to note that oligodendrocyte progenitor cells (OPCs) often couple synaptically to axons via both, the excitatory, glutamatergic Kukley et al., 2007; Ziskin et al., 2007 and inhibitory, GABAergic connections Maldonado and Angulo, 2015. This enables the cells to respond to different patterns of action potentials with different functional effects on cell differentiation and proliferation Nagy et al., 2017, for example, by increasing the number of OL available in the OC. However, while the OPCs can play an important role in modifying myelin content they do not have the needed feedback for adjusting the CV in a manner that eventually leads to a synchronized arrival of spikes at axonal targets. There are other activity-dependent signaling molecules released from axons firing action potentials, notably ATP and adenosine, that activate purinergic receptors, leading to global increases in cell calcium and activation of the myelination promoting genes Stevens et al., 2002; Ishibashi et al., 2006; Fields and Ni, 2010. The mechanisms by which gene expression can exert local changes in a given OL process are complex, as is the case for synaptic modifications, hence suggesting the identity of the reactions and reactants that play the role of the factor M in our model is outside the scope of the present work.

Materials and methods

OL equation and spike-response curves

Request a detailed protocol

The spike response curve, R(t), represents the global response of an OL to a single neuronal spike on any of its myelinated axons and gives OMP its timing-dependent character. In a more general OMP model, there can be several such curves responding to any triggering event in a cascade of responses, as depicted in Figure 1—figure supplement 1. Each of the responses represents the change in the concentration or amplitude of some global factor/signal released in the OL after each spike. A general spike response curve is parameterized by the separate rise and decay times, τr and τd; for a spike occurring at time ts, it can be written as,

(5) R(t,τr,τd,Q | ts=0)=Qτr+τdτd2et/τd(1et/τr),

where Q, represents the single release amount/quantity, which is constant and independent of the current value of a global signal, G(t). The peak of the response is happening at time tmax=τrlogτr+τdτr, reaching the value Gmax(t)=Qτd(τr+τdτr)τrτd. This form allows more general explorations (e.g., τr0, yields the exponentially decaying R(t)), and most importantly it is the impulse response curve of the second-order linear system,

(6) τdτrG¨(t)+(τd+2τr)G˙(t)+(τd+τr)τdG(t)=s(t),

where s(t) represents the signal that drives the system, and G˙(t) and G¨(t) are the first and second time derivatives of G(t).

In order to simplify extensive explorations of all OMP parameters, here we chose a simplified form for R(t) that uses a single characteristic time of the OL spike response, τG=τr=τd, yielding,

(7) R(t)={2τGet/τG(1et/τG),t>00,otherwise.

In Figure 1E, we show a set of such curves for varying values of τG. With this simplification, the differential equation for G(t), corresponding to Equation 6, becomes

(8) τG2G¨(t)+3τGG˙(t)+2G(t)=s(t),

which we use to simulate the global response of an OL to input s(t).

OMP is a continuous time model and the input signal, s(t), appearing in Equation 8, can technically be any integrable input. However, the time-dependent MP will have to rely on events that are sharply defined in time. We use trains of neuronal spikes for individual axons which are prescribed using a particular interspike interval (ISI) distribution, pISI(t(i)), that is, they are generated with a renewal process. An important parameter that characterizes these trains is their mean firing rate, fs, or equivalently the mean ISI, τs=t(i)t=1/fs. We use two main forms for pISI distributions: (a) the exponential distribution with the rate fs=1/τs with added constant refractory time, tR, yielding the refractory Poisson process and (b) regular spiking, spaced at constant intervals, τs. These pISI s cover two extremes: for Poisson spiking, with tR=0, the appearance of the next spike is completely independent of the previous spikes (memoryless process), and for regular spiking, the appearance of the next spike is precisely determined by the last spike. To both of the forms of pISI we also add jitter, specified by σj, which spreads the spike times, such that these temporal shifts are normally distributed according to N(0,σj). Such a mix of refractory Poisson and regular spiking with added jitter seems to cover the spike dynamics for communication between many areas of the brain Maimon and Assad, 2009.

When the spikes on different axons follow the same renewal process, that is, obey the same pISI distribution, we call these ‘pure’ signals, and among them distinguish two different cases: (1) the renewal processes for different axons are fully independent, and (2) they are time-locked via imposed relative time shifts between spikes in different axons, that is the fixed delays. Due to added jitter, which is always independent between different axons, the spikes will not be precisely time-locked, but will still be ‘correlated’. In Figure 2B, we show examples of correlated and independent spike trains, prior to the imposition of fixed delays, for both, Poisson and regular spiking. Mathematically, spike trains are generally formulated as a sequence of Dirac delta functions,

(9) sa(t)=kaδ(ttka),

where tka indicates the time of the kth spike on axon a, and the sum goes over all spikes that have occurred prior to time t. When Equation 9 is applied to the linear system in Equation 8 the analytical solution for G(t) becomes the sum of the responses to spikes on all axons that it myelinates, i.e.,

(10) G(t)=atka<tR(ttka),

which can also be viewed as the sampling of R(t) in reverse time, which we use in our theoretical derivations (see Figure 1F). We do not use Equation 10 directly, but rather solve Equation 8 numerically, as described in the OMP Implementation and Simulations section.

OMP model learning equations

The fluctuating global signal, G(t), obtained via Equation 8, serves as a catalyst for the local myelin-promoting factor, M. We model this by making the increase in M proportional to G(t), as well as to the signal strength in the axon it myelinates. Since M is also continuously converted to myelin with some rate, λA, the differential equation for its concentration on axon a, Ma, can be written as,

(11) Ma˙+λAMa(t)=λMG(t)sa(t),

where λM specifies the production rate of M (Figure 1D, Table 1) and it is the same for all processes. It is an OMP parameter that effectively controls the rate of adaptive changes in our model. The presence of the factor M in a given axon and its effects will lead to the increase in myelin sheath thickness (with rate λA) and will compete with another continuous process of myelin removal which, in the absence of any activity, decreases with some rate, λR.

For most neuronal plasticity models, saturation functions need to be introduced to stabilize the learning process. Similarly, in our implementation of OMP, we introduced two separate saturation functions for myelin addition, FsA, and removal, FsR. The OMP equation for the latency on axon a, τa, can be written as

(12) τa˙=λRFsR(τa(t))-λAFsA(τa(t))Ma(t),

where FsR(τ)=Hr(τmax-τ) and FsA(τ)=Hr(τ-τmin) , Hr(x)=xH(x)/(τmax-τmin) is the normalized ramp function, H(x) is the Heaviside (unit step) function, τmax and τmin are the parameters of the model which specify the maximal and minimal delays that are attainable on any axonal connection, respectively.

We make two modifications to this basic model, one being the case of λA (instantaneous myelination), for which Equation 11 is not needed, and another case includes a homeostatic equation (Equation 17), which presumes that overall myelination for each OL reflects a long-term homeostatic steady-state between adding and removing myelin.

‘Instantaneous’ Myelination

Request a detailed protocol

The conversion rate, λA, appears not to play an important role (Figure 4—figure supplement 2D), particularly if OLs operate far from the saturation limits, as M is then just a currency for conversion into myelin, or, changes in CV. We can eliminate Equation 11 by taking the limit λA, that is, presume that M is instantly converted to myelin, resulting in an immediate change in the time delay, τa. The solution for Ma(t) is obtained by convolving the impulse response of the left side of Equation 11, e-λAt, with the expression on the right side, i.e.,

(13) Ma(t)=λM0teλA(tt)G(t)sa(t)dt+Ma(0)eλAt.

Replacing Equation 13 in Equation 12 and using the fact that sa(t) is a spike train we obtain

(14) τa˙=λRFsR(τa(t))λMtka<tG(tka)FsA(τa(t))λAeλA(ttka).

The expression fλA(t)=λAeλAt, appearing in Equation 14, can be interpreted as a Dirac delta function when λA, since limλAg(t)fλA(ta)dt=g(a). The ‘instantaneous’ equivalent of Equation 12 can then be written as

(15) τa˙=λRFsR(τa(t))λMtka<tFsA(τa(tka))G(tka)δ(ttka),

where it is indicated that τa(t) in the sum will only depend on its values at spike times tka, after Equation 15 is integrated. To have stable integration in this case, it is important to set λM sufficiently small, particularly when used with homeostatic regulation described below.

Homeostatic Regulation

Request a detailed protocol

The form in Equation 15 illustrates that the essence of the adaptive process for timing adjustments is the balance between continuous myelin removal (longer delay) and the discrete increments induced by spikes. For independent Poisson spikes and ignoring saturation effects, the balance condition is,

(16) λR=λMNA/τs2,

where τs is the average inter-spike interval of the Poisson process on a single axon. However, with the saturation functions and when correlated signals are introduced, this homeostatic balance can be disturbed. In order to keep the system in balance, we make the removal rate of myelin, λR, another time dependent variable in the dynamics of the OMP model. To do so, we assume that each OL operates with some local and nominal homeostatic level of myelination, parameterized by some nominal delay, τnom, such that any deviation from it will lead to a slow change in λR according to,

(17) λR˙=λHλR(t)(τnom1NAa=1NAτa)

where λH is a homeostatic rate. We set it to a very small value (λH10-5) so that the time scale for changes in λR is much slower than the time scale of individual spikes or the changes in myelination. In some instances, in order to test the importance of having Equation 17, we simply set λH to 0, after guessing or finding the value for λR that balances the increase in myelin content for a given input signal. This homeostatic rule can be interpreted as a tendency of each OL to conserve its overall amount of myelin, while re-distributing it over different axons. We assume that such an activity-dependent MP process is in place and simulate only its ability to adjust the overall rate of myelin removal, that is the myelin removal rate, λR.

OL Chain

Request a detailed protocol

The consistency of Equation 4 can, for small jitter and fixed delays, cause instabilities as the ‘leading’ axon will keep losing myelin more than other axons. Generally, the spikes on less delayed axons will produce less Ma than spikes on more delayed ones, creating an imbalance. This trend will stop only after the rate of change is significantly slowed down by the saturation limits and the decreasing λR, due to homeostatic regulation. This problem is not surprising as most of the plasticity models are inherently unstable (e.g. Hebbian) and rely on saturation mechanisms. In the case of OMP models, this problem can also be resolved by using a natural assumption that the final arrival time will depend on the action of all OLs along a given axon bundle; the temporal differences of spikes across the bundle will become smaller for OLs close to the target compared to OLs close to the sources. For this reason, it is worth, if not necessary, to simulate the sequential action of multiple OLs, in which the preceding OL-axon bundle segment can pass its modified spike arrival times to the next segment. Hence, we simulate a sequence of OMP equations, each feeding its output to the next segment in the OC (see Figure 2A). The OC will have NO sets of OL equations, each having its own myelin-promoting factors, Ma, and its own local delays, τa(o). As already emphasized, the OC depicted in Figure 2A does not imply literally that there are NO OL cells along the axons, but rather that there are NO segments, representing NO different populations of oligodendrocyte cells myelinating different portions of the axonal bundle, which modulate the delays locally. Assuming that all cells within the same segment will receive the same pattern of spikes and respond to it in the same way, they all can be governed by a single OMP equation. Individual oligodendrocyte cells, in fact, would not be able to modify the delays effectively and independently from other oligodendrocyte cells in the same location, as they would not be able to form tight nodes of Ranvier, considering that OLs prefer not to myelinate the same axon multiple times. Neighboring OL cells are then needed to stack their processes, reducing the width of the nodes of Ranvier and in this way greatly increasing the CV, that is, reducing the conduction delays. Hence, we have a sequence/chain of NO"effective" OL cells, each modifying its local fraction of the total delay along ath axon, τa(o), so that a total conduction delay on the axon is just the sum of all local delays τa=oτa(o).

OL-axon Connectivity

Request a detailed protocol

In general, each OL myelinates many axons and can have multiple processes on a single axon (Figure 1C), or can have none on others. The OL-axon connectivity can be mathematically described with the myelination matrix, , with OLs as rows and axons as columns, and indicates the number of processes a given OL places on each of the axons, for example, for the ‘general’ connectivity in Figure 1C, the matrix is,

(18) M=[324221220302],

while the observed avoidance of OLs myelinating single axons with multiple processes Dumas et al., 2015; Walsh et al., 2016, makes likely to consist of ones and zeros. Here, we use an ‘effective’ OL, which represents a population of OLs behaving in an identical manner, and hence we also treat the OL-axon connectivity in an ‘effective’ way. We currently only address a simple situation in which the OC contains NO OLs that are fully myelinating a given bundle of NA axons, that is, is simply an NO×NA matrix of ones. We assume stable connectivity and that modifications of CV are achieved only through the remodeling of myelin sheaths and the nodes of Ranvier. The issue of precise, cellular-level OL-axon connectivity becomes important in more detailed OMP models in which individual OLs are simulated and in which combined effects on the CV of multiple OLs myelinating the same location along an axon can be addressed; however, such models will be computationally extremely demanding. The issues of random or partial connectivity (the OLs in the chain myelinating different subsets of the NA axons considered) can still be addressed and we expect those to have only less effective but not detrimental synchronization effects. The lower efficiency when using more realistic might not be an issue, considering that the number of OLs in the OC that we simulate is vastly lower than the number of actual cells myelinating axons along a given pathway.

OMP implementation and simulations

In its fullest form, the OMP model can have up to 13 scalar parameters plus the NO×NA myelination matrix . We provide an overview of all OMP parameters in Table 1, together with the OMP variables and other symbols used for specifying the spiking signals conducted along the axons and for the data analysis.

The nominal OMP model is governed by three main Equations 8; 11; 12, expanded with Equation 17, which adds a homeostatic control of λR, that now becomes an OMP variable. Omitting the equation for homeostasis requires fine-tuning the balancing value for λR for a given parameter setting. This model can be simplified further by assuming ‘instantaneous myelination’, in which case Equation 11 can be omitted, when using Equation 15 instead of Equation 12; however, using Equation 15 requires, particularly if used together with Equation 17, that λM be chosen sufficiently small to insure that the integration is stable. If λM is large, sudden jumps in the value of τa will make Equation 17 stochastic. This can lead to negative values for λR or even the delays themselves, which is not realistic.

The standard variables of the model, for a single OL, are G(t), Ma(t), and τa(t), a=1NA, giving a total of 2NA+1 variables per OL, excluding λR. In practice, solving the model will require 2NA+2 variables per OL, since Equation 8 is a 2nd-order differential equation, requiring an auxiliary variable (see the next section). G(t) and Ma are ‘fast’ variables that change on a milliseconds time scale (dictated by τG and impulse responses to the spikes, with mean ISI, τs), while, λR and τa,a=1NA, are ‘slow’ variables whose rate of change is controlled by λH and λA, respectively (λHλA). In Figure 1—figure supplement 2A, B, we show an example of time progression for both fast and slow variables, as well as the synchrony measure, στ (Figure 1—figure supplement 2B, top row) and the epoch-averages of Ma (bottom row), when a set of correlated signals is conducted along OC. Examples shown are for an OL at the beginning of the oligo-chain with NO=5, NA=10, and a single trial, except for στ (dashed lines) and λR (different colors), where the results from independent trials (3 total) are also shown. In Figure 4—figure supplement 3 we similarly show the time-progression for the unstable case, NO=1, as shown in Figure 3A (λH=106), but now for three different values of λH, indicating that the oscillations seen in στ(t) are a result of homeostatic control. We note here, again, that NO=1 does not mean literally that there is a single oligodendrocyte acting on a given axonal bundle but rather a population of oligodendrocytes acting at a particular location along the axon receiving the same pattern of activations and responding to it in the same way. When independent signals are conducted along OC, the time progression of the slow OMP variables displayed only stochastic variations and no clear trends, as shown in Figure 2—figure supplement 2 (apart from initial ‘refocusing’ of all τa to the same mean value, due to homeostatic constraints). This result was consistent across all simulations in which independent signals are used, with either Poisson or regular spiking pISI.

Implementation of the OMP model

Request a detailed protocol

The crucial element of our model implementation is solving a system of differential Equations 8; 11; 12. To solve Equation 8 we implemented a more general case, given by Equation 6, which can be written as a set of first-order equations using an auxiliary variable, v(t),

(19) v˙(t)=a b G(t)(a+b)v(t)+qs >Σtkδ(ttk),G˙(t)=v(t),withG(0)=0,

where constants, a, b, and qs, are defined using the parameters of the general spike response curve (Equation 5, i.e., the characteristic rise, τr, and decay, τd, time constants, and the amplitude of the release, Q) as follows


To use Equation 8 with single characteristic time, for simplicity, as is described in this manuscript, we set the characteristic rise and decay times to be equal, i.e., τr=τd=τG (tmax=τGln2). We also set Q = 1 for all spikes, in which case the constants in Equation 19 simplify to a=2/τG,b=1/τG,qs=2/τG2. Implementation of Equation 11 is straightforward, as it just adds one more first-order differential equation, while implementation of Equation 12 requires special handling of Dirac delta functions. The simplest way to do this is to integrate the OMP equations between subsequent spikes, as the spikes are specified externally, and handle the discontinuities at each spike separately. For more sophisticated/alternative OMP schemes, which might contain a cascade of triggered events, as depicted in Figure 1D, thus, requiring some kind of "events" functionality (functions evaluated when a set of conditions on the variables are satisfied), which is available in many integration packages. Van Rossum and Drake, 1995 solver in Virtanen et al., 2020, scipy.integrate.solve_ivp, and the full simulation code, including the specification of parameters as well as subsequent data analysis, was implemented in Python. Depending on the parameters used, simulations took anywhere between a few hours to more than 10 days (including the repeats, nr[3,10]). Implementation of the same code in C, or using JIT Ansmann, 2018, could make evaluations substantially faster, but might sacrifice some flexibility and ease in implementing the event handling.

OMP model evaluation

Request a detailed protocol

We conducted our simulations on a highly parallel National Institutes of Health Biowulf cluster (http://hpc.nih.gov) and divided them into more than a dozen studies. For each simulation study, we specify a set of values to be explored for each of the OMP parameters, which are listed in Appendix 1. In the early and exploratory phase, we coarsely identified the working range for the most important parameters and usually chose a few values, usually three, indicating the low, mid, and high values of its "working" range, and include sometimes explorations outside that range (as was done in our early, coarse exploration of the model). Due to a large number of parameters that could influence the performance, we could not afford to exhaustively explore all of them on a fine grid for a large range of values. Instead, in each study, we explored the influence of a particular parameter on a finer grid, or in some cases a group of parameters, for example, the exploration of τG and τs, shown in Figure 3E (see the right panel in Appendix 1—table 1).

In our simulations, we chose fixed delays to be randomly drawn from a normal distribution N(Dmean,σD). Here σD is the standard deviation among fixed delays td(a) on all axons and is an important parameter in our simulations, while Dmean is an arbitrary offset, insuring that the delays are positive. This positivity constraint is not very important since we were only interested in the spread of the synchronized spikes on different axons. These fixed delays are always combined with adaptive delays, that is, the conduction delays, τa, on myelinated axons, so we often normalize them to zero mean (Dmean=0), and we label such normalized fixed delay on axon a as Da. The spread of arrival times of synchronized spikes at the target will then simply be the standard deviation of Da+τa, that is, στ=SD(Da+τa). We use στ extensively in this manuscript, as a measure of such spread in arrival times, making it an inverse measure of synchronization with στ=0 indicating perfect synchronization. In order to allow for easier comparison between different sets of parameters, particularly when plotting the average στ over many runs (see Figure 3—figure supplement 1A and the left panel in C), we introduce a normalization step, so that the initial spread due to fixed delays is exactly σD, i.e., Danorm=σDDa/SD(Da).

We initialized a set of NA random fixed delays, parameterized by their spread, σD, for each trial (nr total). We initialized the two NO×NA matrices, one carrying the information about M concentration for each OL process in OC, and the other carries the local delays on each axon for each OL in the OC. The former is initialized to zero (all Ma=0) and the latter is initialized based on the specified mean, τ0 and the percent spread, pτ, i.e., τa˙N(τ,pττ/100). For simplicity, we set τ0=τnomo=τnom/NO, and we used pτ=5 %. We created NA spike trains with prescribed dynamics, based on a given pISI, and parameterized by τs (examples shown in Figure 2B). For time-locked signals, we generated a single train with a given pISI, and others were shifted versions of the same, based on the values of fixed delays (Da, or td(a)), which are subsequently randomized using jitter. The jitter spreads the location of all spikes, such that these temporal deviations are normally distributed according to N(0,σj). For each trial, we ran ne epochs of learning, each with duration Te, and collected the values of the ‘slow’ variables, στ, λR, as well as the mean value during the epoch of τa and Ma, for every OL in the OC. Due to the large number of runs we conducted, only στ information was saved for every run, while other variables were saved only if the analysis required it. When the homeostatic equation was used (λH>0), we allowed 1 or 2 extra (‘warm up’) epochs to run, during which modifications to τa were disabled, allowing λR to be closer to its equilibrium value when we start to track the spread, στ. Due to oscillatory behavior of λR, in most cases, this was not very important to do and did not affect στ(t).

Model fitting for synchronization profiles

Request a detailed protocol

The averages over a large number of runs can obscure the details in synchronization profiles, στ(t), obtained from single runs, such as oscillations and other instabilities. Hence, it is useful to summarize and condense thousands of obtained profiles in an automated fashion and obtain distributions of critical parameters, most importantly the long-time baseline parameter, στ, and the ‘learning time’, τL, which is the inverse of the learning rate, Lτ. We do this by fitting a sufficiently rich model, with additional parameters, pi, that is able to capture most of those profiles reasonably well. The general form of the full (unrestricted) model we use can be written as,

(20) στ(t)στ+(p3exp(t/τL)p4exp(t/p5))f1(t)+f0(t),

where the damped oscillatory behavior of στ(t) is described via functions f0 and f1 (with long-time limit 0 and 1, respectively), used to quantify the observed instabilities in learning. For f0(t) and f1(t), we have used cosine functions with additional parameters describing their amplitude, period, and phase, as well as their damping factor. Hence, the most general fitting model we use, containing 12 parameters, is a double exponential function containing both, the multiplicative and damped additive cosine oscillation. We label the full model as E2C2 and its formula is

(21) στ(t)στ+(p3exp(t/τL)p4exp(t/p5))(1+p6cos(2πt/p7+p8)+p9exp(t/p10)cos(2πt/p11+p12))

We then explore a set of restricted models, all nested within the full model above, starting with the simplest, constant model (C), having only one free parameter, στ, which, consequently is also contained in all other models, with progressively more parameters. These are: single exponential model (E1), double-exponential (E2), double-exponential with multiplicative oscillation (E2C). The full set of models that we fit is then,

  • C:στ(t)στ

  • E1:στ(t)στ+p3exp(-t/τL)

  • E2:στ(t)στ+(p3exp(-t/τL)-p4exp(-t/p5))

  • E2C:στ(t)στ+(p3exp(t/τL)p4exp(t/p5))(1+p6cos(2πt/p7+p8))

  • E2C2: full, unrestricted model, shown in Equation 21.

Our στ(t) consists of npts=ne data points which we use to estimate parameters for all models, via independent fits. We constrained the parameters during the fitting as follows: 0στ2στ(0), Texp/1000<τL<1000Texp, 0<p3<2στ(0), 0<p4<2στ(0), Texp/1000<p5<5Texp, 0<p6<στ(0)/2, Texp/25<p7<, 0<p8<2π, We used these limits to have better stability and to avoid extreme outliers (since the total number of runs was close to 100k, we did not inspect every single fit, but only a small fraction of all fitted parameters values were obtained at the boundary). An example of such fits is shown in Figure 3—figure supplement 2A. We conducted a large number of randomly initialized fits, in order to insure that the best possible fit is obtained (Figure 3—figure supplement 2A, B). Due to the complexity of the unrestricted model, E2C2 might not end up finding the true global minimum and having the lowest MSE, but that happened very rarely (<0.3% or runs, and in all those cases E2C was the one with the minimal MSE). In most cases, different fitting models give very similar estimates for the most important parameter, στ, but in many cases the estimates can differ substantially, depending on what model is selected (see Figure 3—figure supplement 2A, B). Our desire is to opt for a simpler model, when possible, as it often provides more reliable estimates of στ and also τL (unless it is model C, for which τL). This allowed us to categorize better the behavior of OMP and particularly to identify single runs in which the approach to synchrony was truly unstable, or oscillatory, or if στ was diverging away from synchrony. This categorization was implemented by our model selection procedure described below.

Model selection

Request a detailed protocol

All of the simpler/nested models we call restricted since they are essentially equivalent to the unrestricted model with the coefficients of all extra explanatory variables being restricted to zero. Of course, the unrestricted model, having more parameters will then always be able to fit the data better or at least as well as the restricted model, in terms of the mean-squared error (MSE). The question is, whether this improvement is sufficiently large to warrant sacrificing the level of parsimony of the restricted model. One approach to this problem is to use an F-test, which compares two models, the unrestricted (U) vs a particular restricted version (R). Given the obtained fits to npts points for both models and their corresponding residual sum of squares (RSS) one can calculate the F-statistic, F, given by


where RSSm is the residual sum of squares of model m. This F is F-distributed, with (nU-nR, nptsnU) degrees of freedom, which we can use for our statistical tests. The null hypothesis in these tests is that the unrestricted model does not provide a significantly better fit than the restricted model, hence rejecting it at a given significance level means that a more complicated model is needed.

Using the conventional F-test for model selection did not work well for our purposes, since the undulations we observe are nearly always statistically significant even if the most stringent tests with extreme significance levels were used, for example, p<1015. This was not surprising, as the synchronization profiles were never truly statistically constant or exponential. For example, in Figure 3—figure supplement 2C, when στ(t) is observed on a full scale, one would expect that the constant model is the best description of the behavior observed. However, the model selection with regular F-test chooses the less restricted models, in fact E2C2 in this case. After zooming in, we see that the results were indeed not constant, as those deviations from constancy were not just due to noise, but were significant. If the obtained simulation profiles were noisier then the restricted models would have a reasonable chance of being selected. Instead of adding artificial noise to our στ(t), we introduced a parameter, pMSE, with which we essentially control what level of noise or deviations we deem tolerable. The pMSE expresses this level of tolerance as the percentage of the initial, στ, or approximately, σD. Hence, we declare the minimal amount of RSS in any fit, RSSmin=npts*MSEmin and MSEmin=(pMSE*στ(0)/100)2. This sets the level of MSE that is presumed by default, that is, some minimal amount of noise present in residuals of any model. This is essentially specifying how much of RMS error can be tolerated in the restricted model, in order to reject the null hypothesis that the unrestricted model is better. Since the numerator will remain unchanged, this essentially only modifies the denominator,


where RSSmin=npts*MSEmin and MSEmin=(pMSE*στ(0)/100)2. The numerator portion is not affected, as it is a difference between two RSS. This yields the modified F-statistic, Fm,

(22) Fm=FnumFdenm

that we use for our tests. Here we use the modified F-statistic (Equation 22, in most cases with pMSE=2 %). Under the modified test, the curve shown in Figure 3—figure supplement 2C is now declared as model C ("constant") for the three smallest significance levels used (see below).

In our procedure, we start with the model with minimal MSE (usually E2C2), test it against all of the restricted models, and choose the most restricted model for which the null hypothesis is not rejected. We performed the tests separately at different but very low significance levels (α[0.01,0.00001,1010,1015]). The choice of the model was in many cases not strongly influenced by the choice of α, but was strongly dependent on the choice pMSE. For pMSE=0 %, even at extremely low α, the unrestricted model would always be chosen (see Figure 3—figure supplement 2C). Using pMSE=2 %, allowed the choice of model to depend on α in most cases, as is indicated in Figure 3—figure supplement 2A, where E2C2, E2C, or E1 would be chosen, depending on how stringent the test was. We chose pMSE=2 %, based on a set of 100 random examples in which the best model is chosen manually by visual inspection (e.g., the constant model in Figure 3—figure supplement 2C). Note that our model selection is largely ad-hoc and we emphasize that our modified F-test does not aim to provide a quantitative statistical analysis, as use of such absurdly small significance levels indicates, but only to provide a useful quantitative tool for summarizing tens of thousands of runs that we have performed. While the distribution of different models changes significantly for different choices of pMSE and α, the derived values of the parameters στ and τL is not significantly changed when different values of pMSE (but >1%) and α (but smaller than 0.01) were used. The same holds for the ad-hoc rule, of reverting to a less restricted model when MSER>500×MSEU, which happened very infrequently (see Figure 3—figure supplement 2D).

Appendix 1

Tables of parameter values explored in simulations

The values explored for some OMP parameters were similar in different simulations and, to save space, in the tables below we only report the parameter values that differed from these most commonly used default values: NA=10, NO=5, τG={10,20,30} ms, τs={50,100,200} ms, λM={0.01,0.02,0.05,0.1} ms-1, λA={0.1,0.01} ms-1, λH=10-6 ms-2, Te=10 sec, nr=4 and ne=100, σD={5,10} ms, σj={1,3,5} ms, σs=0 %, τmax=100 ms, τmin=3 ms, and τnom=50 ms, tR=0 ms. In two cases the differences with the defaults were small, so we report them only here in the text. The first is studying the influence of the number of axons, NA, for which λA=0.01, NO={5,10}, and NA={5,10,20,30,,70,80,100,150} (n=4752 runs). Runs with σD=5 ms (n=2376) are used in creating Figure 3F. The second explores the influence of σD parameter, i.e., fixed delays, shown in Figure 4F, for which σD={2,3,5,7,10,15,20} ms and nr=5 (n=1512). Below we show the tables for the other 8 sets, where the same units are used as above.

Appendix 1—table 1
Parameter sets used in our OMP model simulations.

Parameter sets in the table on the left were used for exploring the general OMP model behavior with both, correlated and independent spikes, yielding n=2×8640=17280 runs. Simulations with σD=10 ms were used to create Figure 3A, while other runs from this set are used in Figure 3—figure supplement 1 and Figure 4—figure supplement 2. The parameter sets on the right (n=3200) were used in Figure 3E, in which we comprehensively examined the two most important parameters, τG and τs.

λM0.01, 0.02, 0.05, 0.1, 0.2λM0.02, 0.05
λA0.01, 0.1, 1λA0.01
λH10-5, 10-6τG4, 6, 8, …, 82
NA10, 25σD10
NO5, 10τS5, 10, 15, … 195, 200
τs25, 50, 100, 200σj1
Appendix 1—table 2
Parameter sets used in our OMP model simulations.

On the left are parameter sets used to create Figure 3B, C and D exploring the influence of the number of OL, NO, in the OC (n=576 runs). On the right are parameter sets exploring mixed signals shown in Figure 4A, B and C, which were repeated in three scenarios: (1) two equally sized groups of axons, one carrying correlated and the other independent spikes, (2) two equal and separately correlated groups, and (3) one independent and four equal separately correlated groups. Scenario 1 was used in Figure 4A and B, and scenarios 2 and 3 are shown in Figure 4C. Two additional matched sets of runs are performed for "pure" correlated or independent spikes, matched in terms of NA within a group, yielding a total of n=6912+2×2304=11520 runs.

λA0.01, 0.1, 1NA20, 50
NO1, 2, 5, 10NO5, 10
σD5τs25, 50, 100, 200
tR0, 30nr5
τs25, 50, 100, 200
Appendix 1—table 3
Parameter sets used in our OMP model simulations.

On the left are parameter sets used to explore the influence of homeostatic rate (n=2520) shown in Figure 4D and E and Figure 4—figure supplement 2. Parameters used in Figure 4—figure supplement 2E, F, exploring the influence of the spiking rate variability on synchronization (n=1296).

λM0.01, 0.02, 0.05, 0.1, 0.2λM0.02, 0.05, 0.1
λH10−2, 10−3, … , 10−7, 10−8λA0.1
NA10, 20NA10, 20
σD10τG10, 20, 30, 50
σj1, 3σD5, 3, 7
σs0, 1, 2, 5, 10, 20
Appendix 1—table 4
Two parameter sets for simulations used in Figure 1—figure supplement 4A (left) and Figure 1—figure supplement 4B (right), comparing theoretical predictions for the pure Poisson synchronized spikes with σj=0 to the values obtained in simulations with non-zero σj.

For the table on the left, with tR=0 ms there were n=2916 runs and for the one on the right, with refractory Poisson spikes (tR=30 and 80 ms) there were n=288 runs.

λM0.01, 0.1λM0.01, 0.1
λA0.1, 0.01, 0.001λH10-7
λH10-7τG5, 10, 30
NA2, 5, 10σD5
NO1, 5, 10tR30, 80
τG5, 10, 30τs100, 200
σj0.1, 1, 3, 5σj0.1, 1, 3, 5

Appendix 2

Predictions for two axons with regular spiking

In the case of regular spiking, where the ISI times are distributed as pISI(t)=kδ(tτs), Equation 1 reduces to evaluating the sum,

(23) Ga(δta)=k=0R(kτs+δta).

The sum in Equation 23 can be easily evaluated if we note that for our choice of R(t) (Equation 7), the response function can be written as the difference of two exponential functions, (for clarity, abbreviating τG as just τ).

(24) R(t)=etτe2tτ,

which reduces to evaluating two geometric sums, yielding,

(25) GA(δt|τ,τs)=2eτs2δtτ(eδt/τ+eδt+τsτeτs/τ)τ(e2τsτ1),

where δt is the time to the most recent spike, and Ga is replaced with GA, indicating that the form of GA will be the same for different axons, but δt on separate axons will differ. In other words, the contributions to the global G(t) from different axons is only going to differ through the difference of their times to the most recent spikes, δta, i.e., Ga(δta)GA(δta). Assuming t=0 coincides with one of the spikes on a given axon we can write, δtt(modτs). When combining GA(δta) signals from all axons, generally only one of them can be chosen as such reference, while δt for others will be expressed in terms of fixed delays between them, and will require considering separately different orderings, depending on the fixed temporal delays between them. For example, in the simple case, NA=2, for a given fixed delay, td, between two regular spiking trains (Figure 1—figure supplement 3A), we need to distinguish between the cases δttb versus δt>tb, where tb=τstd. In Figure 1—figure supplement 3A and B, we color the "leading" axon (the one whose spikes arrive first) as "green", while the other one, referred to as the "lagging" axon, is colored red. We can see that the time to the last spike will depend on the magnitude of δt1. When δt1tb,


and for δt>tb,


Since t=0 is set for "red" (lagging) spikes, they will occur at times, t=kτs, while "green" spikes will occur at t=kτs+td, and, according to the OMP model, the amount of M will be proportional to G(t) at those times. This is plotted in Figure 1—figure supplement 3B, and we see that the lagging axon will higher concentration of myelin-promoting factor M than the "green" (leading), which will synchronize them. This, however, switches when τs<2td, i.e., for critical spiking frequency, fscrit=1/(2td), which in this case fscrit=50 Hz. Such "reversed myelination" can also be seen in Figure 1—figure supplement 3B, indicated by "orange" and "red" regions. The dashed lines indicate the contours where the ratios are equal to 0.25, 0.75, and 1, as labeled. The r2=1 contour occurs for fscrit. Evaluating Ma for NA> 2 and in the presence of jitter requires more elaborate calculations, that are left to be addressed outside of the current manuscript.

Data availability

All results shown in our figures are produced via computer simulations of our model. The code and the scripts that generated these results are provided on GitHub https://github.com/pajevic/OMPmodel (copy archived at Pajevic et al., 2023). The only data shown that are not the result of our simulations are the images in panels A and B of Figure 1, which are reused with permission, as they also appeared in another publication.


  1. Book
    1. Fields RD
    (2013) Regulation of myelination by functional activity
    In: Kettenmann H, Ransom BR, editors. Neuroglia (3rd edition). Oxford, UK: Oxford University Press. pp. 573–585.
  2. Conference
    1. Pajevic S
    2. Basser PJ
    3. Fields RD
    Models of plasticity and learning employing adaptive temporal delays
    Annual Meeting of The Society for Neuroscience.
  3. Book
    1. Van Rossum G
    2. Drake FL
    Python reference manual
    The Netherlands: Centrum Voor Wiskunde En Informatica Amsterdam.

Article and author information

Author details

  1. Sinisa Pajevic

    Section on Critical Brain Dynamics, National Institute of Mental Health, NIH, Bethesda, United States
    Conceptualization, Resources, Data curation, Software, Formal analysis, Supervision, Validation, Investigation, Visualization, Methodology, Writing – original draft, Project administration, Writing – review and editing
    For correspondence
    Competing interests
    No competing interests declared
    ORCID icon "This ORCID iD identifies the author of this article:" 0000-0002-8880-3320
  2. Dietmar Plenz

    Section on Critical Brain Dynamics, National Institute of Mental Health, NIH, Bethesda, United States
    Resources, Supervision, Funding acquisition, Visualization, Project administration, Writing – review and editing
    Competing interests
    No competing interests declared
    ORCID icon "This ORCID iD identifies the author of this article:" 0000-0002-0008-3657
  3. Peter J Basser

    Section on Quantitative Imaging and Tissue Sciences, Eunice Kennedy Shriver National Institute of Child Health and Human Development, NIH, Bethesda, United States
    Resources, Supervision, Project administration, Writing – review and editing
    Competing interests
    No competing interests declared
    ORCID icon "This ORCID iD identifies the author of this article:" 0000-0003-4795-6088
  4. R Douglas Fields

    Nervous System Development and Plasticity Section, Eunice Kennedy Shriver National Institute of Child Health and Human Development, NIH, Bethesda, United States
    Conceptualization, Supervision, Writing – review and editing
    Competing interests
    No competing interests declared


Intramural Research Program of NIMH/NIH (ZIAMH002797)

  • Sinisa Pajevic
  • Dietmar Plenz

Intramural Research Program of NICHD/NIH (ZIAHD000713)

  • R Douglas Fields

Intramural Research Program of NICHD/NIH (1ZIAHD008972-04)

  • Peter J Basser

The funders had no role in study design, data collection and interpretation, or the decision to submit the work for publication.


This research was supported by the Division of the Intramural Research Programs (DIRP) of the National Institute of Mental Health (NIMH), USA, (SP, DP, ZIAMH002797) and the Eunice Kennedy Shriver National Institute of Child Health and Human Development (NICHD), USA, (RDF, ZIAHD000713; PJB, 1ZIAHD008972-04). This work utilized the computational resources of the NIH HPC Biowulf cluster (http://hpc.nih.gov).

Version history

  1. Preprint posted: July 3, 2022 (view preprint)
  2. Received: July 19, 2022
  3. Accepted: March 15, 2023
  4. Accepted Manuscript published: March 28, 2023 (version 1)
  5. Version of Record published: April 26, 2023 (version 2)


This is an open-access article, free of all copyright, and may be freely reproduced, distributed, transmitted, modified, built upon, or otherwise used by anyone for any lawful purpose. The work is made available under the Creative Commons CC0 public domain dedication.


  • 1,676
  • 317
  • 15

Views, downloads and citations are aggregated across all versions of this paper published by eLife.

Download links

A two-part list of links to download the article, or parts of the article, in various formats.

Downloads (link to download the article as PDF)

Open citations (links to open the citations from this article in various online reference manager services)

Cite this article (links to download the citations from this article in formats compatible with various reference manager tools)

  1. Sinisa Pajevic
  2. Dietmar Plenz
  3. Peter J Basser
  4. R Douglas Fields
Oligodendrocyte-mediated myelin plasticity and its role in neural synchronization
eLife 12:e81982.

Share this article


Further reading

    1. Neuroscience
    Cristina Sáenz de Miera, Nicole Bellefontaine ... Carol F Elias
    Research Article

    The hypothalamic ventral premammillary nucleus (PMv) is a glutamatergic nucleus essential for the metabolic control of reproduction. However, conditional deletion of leptin receptor long form (LepRb) in vesicular glutamate transporter 2 (Vglut2) expressing neurons results in virtually no reproductive deficits. In this study, we determined the role of glutamatergic neurotransmission from leptin responsive PMv neurons on puberty and fertility. We first assessed if stimulation of PMv neurons induces luteinizing hormone (LH) release in fed adult females. We used the stimulatory form of designer receptor exclusively activated by designer drugs (DREADDs) in LeprCre (LepRb-Cre) mice. We collected blood sequentially before and for 1 hr after intravenous clozapine-N-oxide injection. LH level increased in animals correctly targeted to the PMv, and LH level was correlated to the number of Fos immunoreactive neurons in the PMv. Next, females with deletion of Slc17a6 (Vglut2) in LepRb neurons (LeprΔVGlut2) showed delayed age of puberty, disrupted estrous cycles, increased gonadotropin-releasing hormone (GnRH) concentration in the axon terminals, and disrupted LH secretion, suggesting impaired GnRH release. To assess if glutamate is required for PMv actions in pubertal development, we generated a Cre-induced reexpression of endogenous LepRb (LeprloxTB) with concomitant deletion of Slc17a6 (Vglut2flox) mice. Rescue of Lepr and deletion of Slc17a6 in the PMv was obtained by stereotaxic injection of an adeno-associated virus vector expressing Cre recombinase. Control LeprloxTB mice with PMv LepRb rescue showed vaginal opening, follicle maturation, and became pregnant, while LeprloxTB;Vglut2flox mice showed no pubertal development. Our results indicate that glutamatergic neurotransmission from leptin sensitive neurons regulates the reproductive axis, and that leptin action on pubertal development via PMv neurons requires Vglut2.

    1. Neuroscience
    Zahra Ghasemahmad, Aaron Mrvelj ... Jeffrey J Wenstrup
    Research Article

    The basolateral amygdala (BLA), a brain center of emotional expression, contributes to acoustic communication by first interpreting the meaning of social sounds in the context of the listener’s internal state, then organizing the appropriate behavioral responses. We propose that modulatory neurochemicals such as acetylcholine (ACh) and dopamine (DA) provide internal-state signals to the BLA while an animal listens to social vocalizations. We tested this in a vocal playback experiment utilizing highly affective vocal sequences associated with either mating or restraint, then sampled and analyzed fluids within the BLA for a broad range of neurochemicals and observed behavioral responses of adult male and female mice. In male mice, playback of restraint vocalizations increased ACh release and usually decreased DA release, while playback of mating sequences evoked the opposite neurochemical release patterns. In non-estrus female mice, patterns of ACh and DA release with mating playback were similar to males. Estrus females, however, showed increased ACh, associated with vigilance, as well as increased DA, associated with reward-seeking. Experimental groups that showed increased ACh release also showed the largest increases in an aversive behavior. These neurochemical release patterns and several behavioral responses depended on a single prior experience with the mating and restraint behaviors. Our results support a model in which ACh and DA provide contextual information to sound analyzing BLA neurons that modulate their output to downstream brain regions controlling behavioral responses to social vocalizations.