Oligodendrocyte-mediated myelin plasticity and its role in neural synchronization
Abstract
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.
https://doi.org/10.7554/eLife.81982.sa0Introduction
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, , the sign of which determines whether long-term potentiation or depression occurs, with marking the sharp transition between the two. For MP, such discontinuous learning curves would be unstable and hence must be smoothly ramped across the 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 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.
The first essential element of OMP is the release of a global intracellular signaling factor, , after each neural spike on any of the axons myelinated by a given OL. This response has a characteristic transient temporal profile, (Step 1 in Figure 1D). When the OL encounters a sequence of spike trains, the resulting global intracellular signal, , 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, , is also required. It is released after each neuronal spike (Step 2) on a given axon , but its release is catalyzed by and hence is proportional to the global signal, . The OL process that myelinates axon has a dynamically changing concentration of such local factor, , 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 allows for selective modification of the CV and the conduction delays, . One can envision more elaborate OMP models (e.g. OMP-n in Figure 1—figure supplement 1), in which the release of can depend on several different global factors, , which potentially are released via a cascade of events triggered by the original factor . Here, we use the simplest model of this kind, OMP-1, which has only one global signal, , that does both, responding to neuronal spikes and modulating the local release of at axon , 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 into myelin with some addition rate, , and the second being the steady removal of myelin with rate Dutta et al., 2019, resulting in time-varying conduction delays, . The dynamics of the OMP model and its main variables , , and , 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, , 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 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, , to be compensated via corresponding adaptive delays of the myelinated axonal segments, , 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.
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, , 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 . The expression for , which is a contribution to coming from spikes on axon , can then be estimated by recognizing that the sum over spikes in Equation 10 can equivalently be seen as sampling the response function, , in reverse time (see Figure 1G). Here, the important parameter is the time to the last spike, , as the remaining spikes are just the sequential samples drawn from , effectively making a function of , where indicates the spike times on axon for a given OL. If we label the cumulative sum of the subsequent interspike intervals as , the expression for can be written as
For uncorrelated signals, due to symmetry, any of the axons is equally likely to produce a spike, resulting in equidistributed concentrations of guided only by the average concentration, . 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 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 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 and . For example, for the case of two axons () and a given fixed delay, td, we need to consider separately the cases where and where (see Figure 1—figure supplement 3A), in order to derive the expressions for , , , as a function of , and td. For , 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 , regular spiking, and without jitter. The calculations become more cumbersome with increasing , 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 (), we can utilize its memoryless property to derive simple expressions for and . The expression for , at any of the locations along the axons, is simply
The expectation for in the case of Poisson spiking is also simple. For example, when , the average concentration of at the leading edge axon is while for the lagging axon , where . Hence, the ratio of their concentrations is always lower than one,
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 , which happens at low firing rates (large ) and a fast spike response time, (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 on axon as
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 ms (Figure 1—figure supplement 4).
Results
Our measure for quantifying the synchronization properties of the OMP model is the spread in the spike arrival times across all axons, , where stands for pre-specified fixed delays (normalized to some prescribed value, ), and are the adaptive delays (see Figure 2D). The spread starts with some large value, mostly due to the spread among fixed delays, , 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, . We study the performance of the OMP model by characterizing the synchronization profiles, , 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, , 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, which represents the long-time ability of the OMP model to synchronize signals for a given parameter setting.
We show in Figure 3A the estimated distribution of , as well as the model selection chart, when fitting a large number () of synchronization profiles, , 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, , was not monotonic and sometimes appeared oscillatory, which can be seen upon inspecting individual runs, particularly for small (), short , and small jitter, (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 and . 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 profiles grouped by are shown. The runs for small 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, . 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 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 . They all indicate that, with increasing , synchrony is improved and instabilities disappear. For 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 ms to below 3 ms; for ms this increases to 97%, with 74% synchronized below 1 ms).
For correlated spikes, the effectiveness of synchronization strongly depended on two temporal parameters: the characteristic response time of the OL, , and the mean inter-spike interval, , which was assumed to be the same for all axons. In Figure 3E, we explored the synchronization effect as a function of and ; the results shown are for , (see Appendix 1—table 1). For very short ms, performance was unstable, reflecting the fact that for the short-lasting spike responses, , the comparison window between spikes in different axons is too narrow. On the other hand, having spike responses last too long, i.e., ms, makes the resulting too smooth to differentially release , particularly for small (Figure 1—figure supplement 3B, C). Accordingly, we found that for intermediate values of ms, synchronization was predictably achieved and was highly efficient for firing rates Hz.
In Figure 3F, we study the dependence of and the synchronization time constant, , on the number of axons that OLs myelinate, . The results are grouped by the firing rate simulated. For a low firing rate of 5 Hz, optimal performance was achieved for , whereas synchronization became progressively harder with increase in firing rate and number of axons to synchronize. In general, increasing sped up synchronization for all firing rates explored. The fastest synchronization achieved at 5 Hz for an intermediate number of axons () also coincided with the highest fraction of stable, E1, synchronization profiles (peaking at , light-blue). We note that myelinating more than 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 , 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.
For OMP to be operational requires an overall balance between myelin removal (controlled by ) and myelin addition (controlled by and ). For example, if 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 as a variable and control its rate of change with the homeostatic rate, , 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 ( ms), increasing the value of has some detrimental effect on model performance, this effect is not seen for ms. The oscillatory dynamics of depends on and (see Figure 4—figure supplement 3) which can propagate into synchronization profiles. Nevertheless, the pie-chart in Figure 4E does not indicate that strongly influences the model selection, independent of what value of 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 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 , , and . We also explored scenarios in which we use the ‘true’ balancing homeostatic value of , obtained via a trial run, and set . 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 , , while greatly influencing the learning rate, , has negligible influence on (Figure 4—figure supplement 2A and B). Similarly, the conversion rate 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, . For correlated signals, synchronization always occurs but its efficiency decreases in terms of when becomes large (for 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.
Discussion
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 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 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 ) 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 , , and 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 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, , need to have their own homeostatic regulation, so that the axons with higher firing rate will down-regulate their . 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 are mainly implemented here through the rise time, , however, an increase in 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 and the global signal . Our usage of the terms somewhat implies that is a molecular factor, while 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 , , 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 too short can be detrimental in some situations. For firing rates below 10 Hz, should ideally be in the range 10 ms 40 ms for effective synchronization of correlated inputs even for large fixed delays (e.g., ms). These timing requirements make intracellular Ca2+ a good candidate for the role of , 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 ( and ) 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 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 in our model is outside the scope of the present work.
Materials and methods
OL equation and spike-response curves
Request a detailed protocolThe spike response curve, , 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, and ; for a spike occurring at time ts, it can be written as,
where , represents the single release amount/quantity, which is constant and independent of the current value of a global signal, . The peak of the response is happening at time , reaching the value . This form allows more general explorations (e.g., , yields the exponentially decaying ), and most importantly it is the impulse response curve of the second-order linear system,
where represents the signal that drives the system, and and are the first and second time derivatives of .
In order to simplify extensive explorations of all OMP parameters, here we chose a simplified form for that uses a single characteristic time of the OL spike response, , yielding,
In Figure 1E, we show a set of such curves for varying values of . With this simplification, the differential equation for , corresponding to Equation 6, becomes
which we use to simulate the global response of an OL to input .
OMP is a continuous time model and the input signal, , 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, , 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, . We use two main forms for distributions: (a) the exponential distribution with the rate with added constant refractory time, , yielding the refractory Poisson process and (b) regular spiking, spaced at constant intervals, . These s cover two extremes: for Poisson spiking, with , 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 we also add jitter, specified by , which spreads the spike times, such that these temporal shifts are normally distributed according to . 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 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,
where indicates the time of the spike on axon , and the sum goes over all spikes that have occurred prior to time . When Equation 9 is applied to the linear system in Equation 8 the analytical solution for becomes the sum of the responses to spikes on all axons that it myelinates, i.e.,
which can also be viewed as the sampling of 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, , obtained via Equation 8, serves as a catalyst for the local myelin-promoting factor, . We model this by making the increase in proportional to , as well as to the signal strength in the axon it myelinates. Since is also continuously converted to myelin with some rate, , the differential equation for its concentration on axon , , can be written as,
where specifies the production rate of (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 in a given axon and its effects will lead to the increase in myelin sheath thickness (with rate ) and will compete with another continuous process of myelin removal which, in the absence of any activity, decreases with some rate, .
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, , and removal, . The OMP equation for the latency on axon , , can be written as
where and , is the normalized ramp function, is the Heaviside (unit step) function, and 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 (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 protocolThe conversion rate, , appears not to play an important role (Figure 4—figure supplement 2D), particularly if OLs operate far from the saturation limits, as is then just a currency for conversion into myelin, or, changes in CV. We can eliminate Equation 11 by taking the limit , that is, presume that is instantly converted to myelin, resulting in an immediate change in the time delay, . The solution for is obtained by convolving the impulse response of the left side of Equation 11, , with the expression on the right side, i.e.,
Replacing Equation 13 in Equation 12 and using the fact that is a spike train we obtain
The expression , appearing in Equation 14, can be interpreted as a Dirac delta function when , since . The ‘instantaneous’ equivalent of Equation 12 can then be written as
where it is indicated that in the sum will only depend on its values at spike times , after Equation 15 is integrated. To have stable integration in this case, it is important to set sufficiently small, particularly when used with homeostatic regulation described below.
Homeostatic Regulation
Request a detailed protocolThe 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,
where 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, , 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, , such that any deviation from it will lead to a slow change in according to,
where is a homeostatic rate. We set it to a very small value () so that the time scale for changes in 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 to 0, after guessing or finding the value for 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, .
OL Chain
Request a detailed protocolThe 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 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 , 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 sets of OL equations, each having its own myelin-promoting factors, , and its own local delays, . As already emphasized, the OC depicted in Figure 2A does not imply literally that there are OL cells along the axons, but rather that there are segments, representing 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 "effective" OL cells, each modifying its local fraction of the total delay along axon, , so that a total conduction delay on the axon is just the sum of all local delays .
OL-axon Connectivity
Request a detailed protocolIn 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,
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 OLs that are fully myelinating a given bundle of axons, that is, is simply an 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 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 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 , that now becomes an OMP variable. Omitting the equation for homeostasis requires fine-tuning the balancing value for 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 be chosen sufficiently small to insure that the integration is stable. If is large, sudden jumps in the value of will make Equation 17 stochastic. This can lead to negative values for or even the delays themselves, which is not realistic.
The standard variables of the model, for a single OL, are , , and , , giving a total of variables per OL, excluding . In practice, solving the model will require variables per OL, since Equation 8 is a -order differential equation, requiring an auxiliary variable (see the next section). and are ‘fast’ variables that change on a milliseconds time scale (dictated by and impulse responses to the spikes, with mean ISI, ), while, and , are ‘slow’ variables whose rate of change is controlled by and , respectively (). 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 (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 , , and a single trial, except for (dashed lines) and (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, , as shown in Figure 3A (), but now for three different values of , indicating that the oscillations seen in are a result of homeostatic control. We note here, again, that 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 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 .
Implementation of the OMP model
Request a detailed protocolThe 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, ,
where constants, , , and qs, are defined using the parameters of the general spike response curve (Equation 5, i.e., the characteristic rise, , and decay, , time constants, and the amplitude of the release, ) 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., (). We also set = 1 for all spikes, in which case the constants in Equation 19 simplify to . 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, ). 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 protocolWe 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 and , 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 . Here is the standard deviation among fixed delays on all axons and is an important parameter in our simulations, while 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, , on myelinated axons, so we often normalize them to zero mean (, and we label such normalized fixed delay on axon as . The spread of arrival times of synchronized spikes at the target will then simply be the standard deviation of , that is, . We use extensively in this manuscript, as a measure of such spread in arrival times, making it an inverse measure of synchronization with 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 , i.e., .
We initialized a set of random fixed delays, parameterized by their spread, , for each trial (nr total). We initialized the two matrices, one carrying the information about 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 ) and the latter is initialized based on the specified mean, and the percent spread, , i.e., . For simplicity, we set , and we used %. We created spike trains with prescribed dynamics, based on a given , and parameterized by (examples shown in Figure 2B). For time-locked signals, we generated a single train with a given , and others were shifted versions of the same, based on the values of fixed delays (, or ), which are subsequently randomized using jitter. The jitter spreads the location of all spikes, such that these temporal deviations are normally distributed according to . For each trial, we ran ne epochs of learning, each with duration , and collected the values of the ‘slow’ variables, , , as well as the mean value during the epoch of and , 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 (), we allowed 1 or 2 extra (‘warm up’) epochs to run, during which modifications to were disabled, allowing to be closer to its equilibrium value when we start to track the spread, . Due to oscillatory behavior of , in most cases, this was not very important to do and did not affect .
Model fitting for synchronization profiles
Request a detailed protocolThe averages over a large number of runs can obscure the details in synchronization profiles, , 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’, , which is the inverse of the learning rate, . 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,
where the damped oscillatory behavior of is described via functions f0 and f1 (with long-time limit 0 and 1, respectively), used to quantify the observed instabilities in learning. For and , 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
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:
E1:
E2:
E2C:
E2C2: full, unrestricted model, shown in Equation 21.
Our consists of data points which we use to estimate parameters for all models, via independent fits. We constrained the parameters during the fitting as follows: , , , , , , , , We used these limits to have better stability and to avoid extreme outliers (since the total number of runs was close to , 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 (unless it is model C, for which ). 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 protocolAll 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, , given by
where is the residual sum of squares of model . This is -distributed, with (, ) 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, . This was not surprising, as the synchronization profiles were never truly statistically constant or exponential. For example, in Figure 3—figure supplement 2C, when 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 , we introduced a parameter, , with which we essentially control what level of noise or deviations we deem tolerable. The expresses this level of tolerance as the percentage of the initial, , or approximately, . Hence, we declare the minimal amount of RSS in any fit, and . 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 and . The numerator portion is not affected, as it is a difference between two RSS. This yields the modified F-statistic, ,
that we use for our tests. Here we use the modified F-statistic (Equation 22, in most cases with %). 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 (. The choice of the model was in many cases not strongly influenced by the choice of , but was strongly dependent on the choice . For %, even at extremely low , the unrestricted model would always be chosen (see Figure 3—figure supplement 2C). Using %, 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 %, 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 and , the derived values of the parameters and is not significantly changed when different values of (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 , 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: , , ms, ms, ms-1, ms-1, ms-2, sec, and , ms, ms, %, ms, ms, and ms, 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, , for which , , and ( runs). Runs with ms () are used in creating Figure 3F. The second explores the influence of parameter, i.e., fixed delays, shown in Figure 4F, for which ms and (). Below we show the tables for the other 8 sets, where the same units are used as above.
Appendix 2
Predictions for two axons with regular spiking
In the case of regular spiking, where the ISI times are distributed as , Equation 1 reduces to evaluating the sum,
The sum in Equation 23 can be easily evaluated if we note that for our choice of (Equation 7), the response function can be written as the difference of two exponential functions, (for clarity, abbreviating as just ).
which reduces to evaluating two geometric sums, yielding,
where is the time to the most recent spike, and is replaced with , indicating that the form of will be the same for different axons, but on separate axons will differ. In other words, the contributions to the global from different axons is only going to differ through the difference of their times to the most recent spikes, , i.e., . Assuming coincides with one of the spikes on a given axon we can write, . When combining signals from all axons, generally only one of them can be chosen as such reference, while 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, , for a given fixed delay, td, between two regular spiking trains (Figure 1—figure supplement 3A), we need to distinguish between the cases versus , where . 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 . When ,
and for ,
Since is set for "red" (lagging) spikes, they will occur at times, , while "green" spikes will occur at , and, according to the OMP model, the amount of will be proportional to 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 than the "green" (leading), which will synchronize them. This, however, switches when , i.e., for critical spiking frequency, , which in this case 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 contour occurs for . Evaluating for 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.
References
-
Motor learning promotes remyelination via new and surviving oligodendrocytesNature Neuroscience 23:819–831.https://doi.org/10.1038/s41593-020-0637-3
-
Dynamics of self-organized delay adaptationPhysical Review Letters 82:1594–1597.https://doi.org/10.1103/PhysRevLett.82.1594
-
Myelination: an overlooked mechanism of synaptic plasticity?The Neuroscientist 11:528–531.https://doi.org/10.1177/1073858405282304
-
White matter in learning, cognition and psychiatric disordersTrends in Neurosciences 31:361–370.https://doi.org/10.1016/j.tins.2008.04.001
-
BookRegulation of myelination by functional activityIn: Kettenmann H, Ransom BR, editors. Neuroglia (3rd edition). Oxford, UK: Oxford University Press. pp. 573–585.https://doi.org/10.1093/med/9780199794591.001.0001
-
A new mechanism of nervous system plasticity: activity-dependent myelinationNature Reviews. Neuroscience 16:756–767.https://doi.org/10.1038/nrn4023
-
Neuronal activity biases axon selection for myelination in vivoNature Neuroscience 18:683–689.https://doi.org/10.1038/nn.3992
-
Leukemia inhibitory factor regulates the timing of oligodendrocyte development and myelination in the postnatal optic nerveJournal of Neuroscience Research 87:3343–3355.https://doi.org/10.1002/jnr.22173
-
Vesicular glutamate release from axons in white matterNature Neuroscience 10:311–320.https://doi.org/10.1038/nn1850
-
Calcium signaling in the oligodendrocyte lineage: regulators and consequencesAnnual Review of Neuroscience 43:163–186.https://doi.org/10.1146/annurev-neuro-100719-093305
-
ConferenceModels of plasticity and learning employing adaptive temporal delaysAnnual Meeting of The Society for Neuroscience.
-
Preservation of a remote fear memory requires new myelin formationNature Neuroscience 23:487–499.https://doi.org/10.1038/s41593-019-0582-1
-
The axo-myelinic synapseTrends in Neurosciences 34:393–400.https://doi.org/10.1016/j.tins.2011.06.004
-
BookPython reference manualThe Netherlands: Centrum Voor Wiskunde En Informatica Amsterdam.
-
Vesicular release of glutamate from unmyelinated axons in white matterNature Neuroscience 10:321–330.https://doi.org/10.1038/nn1854
Article and author information
Author details
Funding
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.
Acknowledgements
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).
Copyright
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.
Metrics
-
- 1,812
- views
-
- 335
- downloads
-
- 19
- citations
Views, downloads and citations are aggregated across all versions of this paper published by eLife.
Download links
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)
Further reading
-
- Computational and Systems Biology
- Neuroscience
Diffusional kurtosis imaging (DKI) is a methodology for measuring the extent of non-Gaussian diffusion in biological tissue, which has shown great promise in clinical diagnosis, treatment planning, and monitoring of many neurological diseases and disorders. However, robust, fast, and accurate estimation of kurtosis from clinically feasible data acquisitions remains a challenge. In this study, we first outline a new accurate approach of estimating mean kurtosis via the sub-diffusion mathematical framework. Crucially, this extension of the conventional DKI overcomes the limitation on the maximum b-value of the latter. Kurtosis and diffusivity can now be simply computed as functions of the sub-diffusion model parameters. Second, we propose a new fast and robust fitting procedure to estimate the sub-diffusion model parameters using two diffusion times without increasing acquisition time as for the conventional DKI. Third, our sub-diffusion-based kurtosis mapping method is evaluated using both simulations and the Connectome 1.0 human brain data. Exquisite tissue contrast is achieved even when the diffusion encoded data is collected in only minutes. In summary, our findings suggest robust, fast, and accurate estimation of mean kurtosis can be realised within a clinically feasible diffusion-weighted magnetic resonance imaging data acquisition time.
-
- Neuroscience
The retina transforms patterns of light into visual feature representations supporting behaviour. These representations are distributed across various types of retinal ganglion cells (RGCs), whose spatial and temporal tuning properties have been studied extensively in many model organisms, including the mouse. However, it has been difficult to link the potentially nonlinear retinal transformations of natural visual inputs to specific ethological purposes. Here, we discover a nonlinear selectivity to chromatic contrast in an RGC type that allows the detection of changes in visual context. We trained a convolutional neural network (CNN) model on large-scale functional recordings of RGC responses to natural mouse movies, and then used this model to search in silico for stimuli that maximally excite distinct types of RGCs. This procedure predicted centre colour opponency in transient suppressed-by-contrast (tSbC) RGCs, a cell type whose function is being debated. We confirmed experimentally that these cells indeed responded very selectively to Green-OFF, UV-ON contrasts. This type of chromatic contrast was characteristic of transitions from ground to sky in the visual scene, as might be elicited by head or eye movements across the horizon. Because tSbC cells performed best among all RGC types at reliably detecting these transitions, we suggest a role for this RGC type in providing contextual information (i.e. sky or ground) necessary for the selection of appropriate behavioural responses to other stimuli, such as looming objects. Our work showcases how a combination of experiments with natural stimuli and computational modelling allows discovering novel types of stimulus selectivity and identifying their potential ethological relevance.