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
4 figures, 5 tables and 1 additional file

Figures

Figure 1 with 4 supplements
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.

Figure 1—figure supplement 1
OMP-1 model vs OMP-n model.

A more complex OMP models can be constructed; however their structure differs from the simplest OMP-1 model only in step 2, where multiple global factors can participate, together with a cascade of events, which can be triggered by given threshold events, each potential having their own response functions, Rk(t). Ultimately, the whole cascade will lead to the production of the myelination-promoting factor M, as in the case of OMP-1.

Figure 1—figure supplement 2
Time progression of OMP fast and slow variables.

Time progression of OMP variables when correlated Poisson spikes are conducted through OC with a single OL, NO=1, and for NA=10 axons. The ordering of the spikes carries statistically the signature of the fixed delays. The OMP parameters for this simulation are τG=30 ms, τs=100 ms, λA=0.1 ms-1, σj=5 ms, λM=0.01 ms-1, and λH=10-6 ms-1. (A) A millisecond time-scale dynamics of G(t) and Ma(t), in response to incoming spikes (yellow vertical lines). For clarity, only five of the total NA=10 axons are shown. Time reported in milliseconds is just an offset from the beginning of the last Te=10 second epoch. (B) Time progression of στ (top row) and OMP model variables, as indicated, over Texp=5000 seconds, plotted as the averages over Te=10 second epochs, and different columns showing different values of λM, as indicated. The second row from the top shows the 10 delays for each axon, for a single trial, combined with the fixed delays (showing the convergence of the arrival times). The third row shows the time evolution of the myelin removal rate, λR, guided by the homeostatic equation (Equation 17) with black and red dashed lines indicating the simple homeostatic balance estimate (λMNA/τs2) and twice that value, respectively. Differently colored trajectories are the results from three independent trials, with the same parameters but different random initializations of fixed and adaptive delays (σD=5 ms). The bottom row shows the concentrations of M for each axon in a single trial. For the relatively large, τG= 30 ms, differences in concentrations are not as drastic as they are for τG= 10 ms (see Figure 2—figure supplement 1).

Figure 1—figure supplement 3
Theoretical predictions of OMP model.

Theoretical predictions for Ma for σj=0,NA=2. (A) Schematic depiction of the temporal relationship between spikes in a pair of axons (A1, A2; NA=2) in the case of regular spiking, illustrating different scenarios when calculating G(t). The ‘leading’ edge spikes (green) are offset by a fixed delay, td, over the ‘lagging’ spikes (red). (B) Examples of G(t) for td=10 ms and for five different spiking frequencies fs[10,30,40,50,60] Hz, as indicated, showing also the corresponding ratios, r2=M1/M2=Mgreen/Mred. Occurrences of ‘red’ and ‘green’ spikes are plotted together with G(t), indicating also the level of M released, as it is proportional to G(t). We see that up to the critical frequency, fscrit=50 Hz, the lagging axon will have more myelination factor, M, as desired, but this reverses above the critical frequency. (C) Theoretical predictions for the ratio, r2=M1M2., as a function of τG and τs for three different fixed time-delays, td= 5, 10, and 20 ms (rows, as indicated), in the case of regular (left column) and Poisson (right column) spiking. The regions with red/orange hues are the regions of instability, where, r2>1 which will lead to further desynchronization, and is only happening for regular spiking. For Poisson spiking, the ratio is always lt1, as is evident from Equation 3, hence, the ‘lagging’ axon will always have enhanced myelination relative to the ‘leading’ axon, due to a higher concentration of, M which is desired for stable synchronization in all regimes.

Figure 1—figure supplement 4
Theoretical predictions vs simulations.

Comparing theoretical predictions with simulations (σj>0,NA=10). The mean Ma values obtained in simulations (y-axis) are compared to the theoretical prediction (x-axis) given by Equation 4, which ignores jitter. (A) The mean Ma values are from a simulation sets specified in Appendix 1—table 4, for which correlated Poisson spiking was used. The results are shown as a grid plot for different values of jitter and characteristic spike response time, τG. Individual points from the same parameter set are colored the same, with different parameter sets/runs being colored by cycling through 15 different colors with no particular meaning attached to the coloration. Each set had 2, 5, or 10 points (one per axon). We see that the predictions made for σj=0 case hold very well for small jitter (σj1), but fail for larger jitters, particularly for σj=5 ms, making the differences in Ma across axons smaller. Nevertheless, even for large jitter, when τG is large enough, the prediction given by Equation 4 is in reasonable agreement with the observed/simulated values. (B) Same as in (A), but this time spikes were simulated using a refractory Poisson process with tR=30 and tR=80 ms. Because now it violates the assumption of a memoryless process, under which the derivation of Equation 4 is made, the predictions fail for any amount of jitter or τG. We note that theoretical predictions for a general renewal process with any amount of jitter are outside the scope of the current manuscript.

Figure 2 with 2 supplements
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.

Figure 2—figure supplement 1
Time progression of OMP variables for correlated spikes.

Time progression of the slow variables in the OMP model, similar to one shown in Figure 1—figure supplement 2B, but now with faster spike responses in the OL, τG=10 ms, and for different λH (values shown in columns), using only one myelin production rate, λM=0.02, and σj=1 ms. The independent trials, each have with their own set of randomly chosen fixed delays (σD=6 ms). The top row shows the time progression of στ for each trial (dotted red line is the mean over all trials). The second row are the 10 delays, τa, for each axon, and the third row are the same delays when combined with the fixed delays (showing the convergence of the arrival times). The fourth row shows the time evolution of λR, with black and red dashed lines indicating values λMNA/τs2 (the balancing estimate for independent spikes) and twice that value, respectively. The bottom row shows the concentrations of M for each axon.

Figure 2—figure supplement 2
Time progression of OMP variables for independent spikes.

Time progression of the slow variables (averages over Te=10 s epochs) in the OMP model when the spikes are independent with Poisson ISI. Here, the columns show three independent trials with the same set of OMP parameters but each with its own set of randomly chosen fixed delays (σD=5 ms). The OMP parameters for this simulation are τG=10 ms, τs=100 ms, λM=0.01, λH=10-7 σj=1 ms, NO=10 and NA=10 axons. The top row shows the time progression of στ for each trial (dotted red line is the mean over all trials). The second row are the 10 axonal delays, τa and the third row are the same delays when combined with the fixed delays (showing the delay profile that will be fed into the next OL segment in the chain). The fourth row shows the time evolution of λR, with black and red dashed lines indicating values λMNA/τs2 (the balancing estimate for independent spikes) and twice that value, respectively. The bottom row displays the concentrations of M for each axon.

Figure 3 with 2 supplements
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).

Figure 3—figure supplement 1
Exploring OMP properties and parameters.

Exploring OMP with a large number of parameter settings. (A) For any given setting of OMP parameters, simulations are repeated nr times each with a random choice of fixed delays (Da, parameterized with σD) and initial adaptive delays (parameterized with τ0 and pτ). These individual temporal profiles/trials (colored lines) are averaged and the mean is reported as the synchronization profile, στ(t) (that we refer to as a set or run). (B) A grid plot of στ(t) for correlated spikes using 12 parameter runs (3 trials and the mean), revealing its dependence on NO and σj (4 NO and 3 σj values, for NA = 10, τG=10 ms, λH=0.000001, λM=0.02). Some of these profiles are shown in Figure 3B to illustrate OMP instabilities with small NO and fixed delays, which are overcome with increasing σj and NO. (C) The center panel shows all runs specified in the left panel of Appendix 1—table 1. We summarize such large number of runs in two ways: (1) we average all στ(t) within a subset of runs grouped by some parameter, τG, as shown in the left panel; (2) we also fit all στ(t) to 5 models, as described in Materials and Methods, and plot the distribution of στ as shown in the right panel. In some cases we also plot the distribution of τL parameters. (D) Figure 3A results with σD=10 ms are shown here for σD=5 ms and also (E) compared to the equivalent runs using the instantaneous myelination OMP model. We used a smaller bandwidth for kernel density estimation to highlight the differences, which are not apparent at smoother KDE values. (F) Summary of average profiles grouped by NA, shown in Figure 3F for σD=5 ms, and additionally for σD=10 ms.

Figure 3—figure supplement 2
Model fitting for synchronization profiles.

Model fits for quantifying synchronization performance across different parameter runs.

(A) An example of fits for all five nested models. Dashed lines represent estimates of the long-time baseline, στ. (B) In most cases, comparing the dashed lines, the model selection does not affect the estimate of the baseline (top panel), but can sometimes have a significant impact. For example, in the bottom panel, the baselines differ by nearly 50%. For such a complex synchronization profile it is actually hard to determine the true long-time baseline and the only way to resolve those situations is to run really long simulations, which we could not afford in this early and extensive exploratory phase. (C) In the majority of runs, the less restricted models were deemed the best with the standard F-test. When magnified, it is clear that E2C2 offers a statistically better description than the constant model, which we label (C). On the other hand, when considered at the full scale (small graph on the right), we desire to classify such a profile simply as the constant model, (C). This is achieved by introducing a fudge factor, pMSE=2%, which under our modified F-test selects the model C. (D) We also introduced a ‘swapping’ rule. This rule reverts the fit back to the less restricted model when the true MSE in the restricted model is 500 times larger than the one for the restricted model. This did not seem to be consequential, as in most simulations very few swaps occurred, and when they occurred seemed to be justified, as shown in the large panel on the left, where the E1 model suggestion, given by the modified F-test, was reverted to a double exponential fit. The four panels on the right show all swaps obtained with runs specified in Oligodendrocyte-mediated Myelin Plasticity and its Role in Neural Synchronization Appendix 1—table 3, which was the case that had most of the swaps, and most of those were swaps from C to E1 (lower right panel) which had a very small range. Hence, the swap did not change στ significantly but had a large effect on the value of τL, which is hence less reliable in our summary statistics.

Figure 4 with 3 supplements
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.

Figure 4—figure supplement 1
Supplemental results for Figure 4.

(A) Individual runs for the averages reported in Figure 4A for mixed signals, one conducting correlated spikes and the other independent spikes. Individual runs for independent spikes here are manually offset for clarity, not to overlap with profiles obtained with the correlated spikes, but are mainly shown to demonstrate that the actual spread in the performance is much larger than what the shaded standard error (SE) indicates for the independent group, appearing only as lines when averaging over many runs (and drawn as a ‘thick’ purple line for pure independent signals in Figure 4A, where the SE was even smaller). (B) Influence of σD on synchronization efficiency. (C–D) Averages over all runs with σD=10 ms, grouped by jitter, comparing (C) pure signals with (D) mixed signals with one independent and one correlated group, as was similarly shown in Figure 4A for σD=5 ms.

Figure 4—figure supplement 2
Influence of additional OMP and signal parameters.

Exploring additional OMP and signal parameters. (A) Average στ(t), grouped by λM, which highly influences the learning rate, as it sets the time scale for changes in myelination. (B) Bottom panel shows this explicitly via ‘learning time’, τL=1/Lτ, as a function of λM for σD= 5 and 10 ms. The dependence of στ on λM (top panel) was not strong and it appears to slightly favor faster rates. This ‘slight’ trend does not appear consistently and depends on the range of parameters chosen for averaging. (C) The average στ(t) for σD=10 and σD=5 ms, respectively, show that λH does not have a strong influence on the behavior of the OMP model (also emphasized in Figure 4E). (D) Average στ(t) for σD=10 and σD=5 ms now grouped by the myelin conversion rate, λA, including the instantaneous myelination model, Equation 15, and showing that λA has a very small influence on synchronization. This might change when the saturation limits τmin, τmax become narrower, as well when different τnom values in Equation 17 are used. We use saturation limits that allow highly heterogeneous values for the delays (3 ms to 100 ms), making it harder to achieve synchronization. (E) The sensitivity of OMP to firing rate variability, specified via σs in the case of regular spiking, having ‘beating’ synchronized periods, but also even longer non-synchronized periods, potentially detrimental to the synchrony. Averages, grouped by σs, indicate that when the spiking-rate variability is smaller than 5%, synchronization still occurs but is very inefficient. This might not be a problem since myelin plasticity operates on a very slow time scale, measured in days and months, but only if στ does not get worse during non-synchronized epochs. (F) The sensitivity to σs quantified via στ, showing that variations greater than 5% exhibit de-synchronization effects but these might not arise for Poisson spiking.

Figure 4—figure supplement 3
Oscillations due to homeostatic regulation.

Oscillations in λR due to homeostatic regulation. The theoretical guess for the balance between addition and removal of myelin, that we initialize λR to at the start of each simulation, should work well for independent spikes, if saturation effects can be ignored. However, in most cases, this prediction does not match the true balancing value for λR, and hence for fast homeostatic rates, the approach to equilibrium will be too fast and will induce oscillations. (A) Time progression of λR in the OMP model for three different values of λH (columns) and three different values of λM when a set of correlated spikes with exponential ISI (Poisson) are conducted through = NA 10 axons. For this set of runs NO=1, τG=30, and σj=5. (B) Time progression of λR in the OMP model for three different values of λH (columns) and three different values of τG when a set of correlated spikes with exponential ISI (Poisson) are conducted through = NA 10 axons. For this set of runs, NO=1,σj=5 and λM=0.02. The spike-response time, τG, does not influence the oscillation frequency of λR as strongly as λM.

Tables

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
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.

ParameterValuesParameterValues
λ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
nr5nr10
Te60
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.

ParameterValuesParameterValues
λA0.01, 0.1, 1NA20, 50
NO1, 2, 5, 10NO5, 10
σD5τs25, 50, 100, 200
tR0, 30nr5
τs25, 50, 100, 200
ne500/NO
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).

ParameterValuesParameterValues
λ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
nr5σj1
σs0, 1, 2, 5, 10, 20
Te20
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.

ParameterValuesParameterValues
λ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
nr3nr2
ne1000ne1000
Te20Te5

Additional files

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
(2023)
Oligodendrocyte-mediated myelin plasticity and its role in neural synchronization
eLife 12:e81982.
https://doi.org/10.7554/eLife.81982