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 $\mu m$ distance from the soma, but with a tendency to maximize the number of axons it can myelinate, making it unlikely that a given OL would myelinate the same axon twice Dumas et al., 2015; Walsh et al., 2016. We postulate that this morphology enables a single OL to integrate and compare signals from different axons and act as a mediator in providing the needed feedback for adjusting the relative timing of different signals via dynamic regulation of the CV along these axons. In Figure 1D, we outline a general form for such spike-response OMP models, consisting of three essential steps, described below. In these continuous-time models, the transient nature of the temporal profile of the OL response to an action potential plays a crucial role in creating the needed reference for adjusting the relative delays on different axons.

Figure 1 with 4 supplements see all

Download asset Open asset

The first essential element of OMP is the release of a global intracellular signaling factor, $G$, after each neural spike on any of the axons myelinated by a given OL. This response has a characteristic transient temporal profile, $R(t)$ (Step 1 in Figure 1D). When the OL encounters a sequence of spike trains, the resulting global intracellular signal, $G(t)$, will simply be the sum of the individual responses and will fluctuate in time, thus providing a common and time-dependent reference to all of its processes. To allow for differential myelination between different axons, a local myelin-promoting factor, $M$, is also required. It is released after each neuronal spike (Step 2) on a given axon $a$, but its release is catalyzed by $G$ and hence is proportional to the global signal, $G(t)$. The OL process that myelinates axon $a$ has a dynamically changing concentration of such local factor, $M_a(t)$, that depends on the temporal profile of the spike trains and will generally differ between the axons myelinated by a given OL. This difference in local concentration of $M$ allows for selective modification of the CV and the conduction delays, ${\tau }_a$. One can envision more elaborate OMP models (e.g. OMP-n in Figure 1—figure supplement 1), in which the release of $M$ can depend on several different global factors, $G_M$, which potentially are released via a cascade of events triggered by the original factor $G$. Here, we use the simplest model of this kind, OMP-1, which has only one global signal, $G$, that does both, responding to neuronal spikes and modulating the local release of $M_a$ at axon $a$, and for the remainder of this work, we are going to refer to the OMP-1 model simply as the OMP model. The last step (Step 3) represents two continuous processes, one being the conversion of $M$ into myelin with some addition rate, ${\lambda }_A$, and the second being the steady removal of myelin with rate ${\lambda }_R$ Dutta et al., 2019, resulting in time-varying conduction delays, ${\tau }_a(t)$. The dynamics of the OMP model and its main variables $G(t)$, $M_a(t)$, and ${\tau }_a(t)$, is governed by a set of equations that implement steps 1–3 (Equation 8, Equation 11, and Equation 12), but also include presumed homeostatic regulation of the myelin conversion/removal rates (Equation 17), which is needed for the long-time stability of the model. An example of their time progression is shown in Figure 1—figure supplement 2. These equations govern the behavior of a population of OL at a particular segment located at distance x_o along the axonal bundle (see Figure 1G). The full OMP model simulates the net delays across a discrete number of such segments, $N_O$, 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 $N_A$ 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, $t_d^{(a)},a=1,\mathrm{\dots },N_A$, to be compensated via corresponding adaptive delays of the myelinated axonal segments, ${\displaystyle {\tau }_a,a=1,\dots ,N_A}$, to facilitate synchronous arrival, as depicted graphically in Figure 2D. We test this ability of the OMP model in our simulations that are, together with the details of the model, described in Materials and Methods.

Figure 2 with 2 supplements see all

Download asset Open asset

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, $p_{\text{ISI}}(t_{(i)})$, that is, are generated by a renewal process. We denote the mean spiking/firing rate as f_s and the mean inter-spike time interval with ${\tau }_s=1/f_s$. The expression for $G_a(t)$, which is a contribution to $G(t)$ coming from spikes on axon $a$, can then be estimated by recognizing that the sum over spikes in Equation 10 can equivalently be seen as sampling the response function, $R(t)$, in reverse time (see Figure 1G). Here, the important parameter is the time to the last spike, $\delta t$, as the remaining spikes are just the sequential samples drawn from $p_{\text{ISI}}(t_{(i)})$, effectively making $G_a(t)$ a function of ${\displaystyle \delta t_a=\text{min}\left|t-t_{k_a}\right|,\mathrm{\forall }{t}_{k_a}<t}$, where $t_{k_a}$ indicates the spike times on axon $a$ for a given OL. If we label the cumulative sum of the subsequent interspike intervals as ${\displaystyle t_k^c=\sum _i^k{t}_{(i)}}$, the expression for $G_a(t)$ 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 $M$ guided only by the average concentration, ${\displaystyle G_{\text{av}}={⟨G(t)⟩}_t}$. Predicting the synchronization effects for time-locked signals when the OL myelinates many axons is a more difficult task, particularly when Gaussian jitter with spread ${\sigma }_j$ is added to the specified renewal process ISIs and also due to non-linear saturation effects in the learning equation (see Materials and methods).

When combining $G_a$ from all axons, mutual ordering of spikes will generally need to be considered. In most cases, this derivation will depend on the exact form of $R(t)$ and $p_{\text{ISI}}(t_{(i)})$. For example, for the case of two axons ($N_A=2$) and a given fixed delay, t_d, we need to consider separately the cases where ${\displaystyle \delta t\le t_b={\tau }_s-t_d}$ and where ${\displaystyle \delta t>t_b}$ (see Figure 1—figure supplement 3A), in order to derive the expressions for $G_{\text{av}}$, ${\displaystyle ⟨M_1⟩}$, ${\displaystyle ⟨M_2⟩}$, as a function of ${\tau }_s,{\tau }_G$, and t_d. For ${\displaystyle {\tau }_s<2t_d}$, 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 $N_A=2$, regular spiking, and without jitter. The calculations become more cumbersome with increasing $N_A$, 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 (${\displaystyle p_{\text{ISI}}(t)=e^{-t/{\tau }_s}}$), we can utilize its memoryless property to derive simple expressions for $G_{\text{av}}$ and ${\displaystyle ⟨M_a⟩}$. The expression for $G_{\text{av}}$, at any of the locations along the axons, is simply

The expectation for $M$ in the case of Poisson spiking is also simple. For example, when $N_A=2$, the average concentration of $M$ at the leading edge axon is ${\displaystyle ⟨M_1(t)⟩=C_M{G}_{av}}$ while for the lagging axon ${\displaystyle ⟨M_2(t)⟩=C_M(G_{av}+R(t_d))}$, where $C_M={\lambda }_M/({\lambda }_A{\tau }_s)$. 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 ${\displaystyle r_2\ll 1}$, which happens at low firing rates (large ${\tau }_s$) and a fast spike response time, ${\tau }_G$ (Figure 1—figure supplement 3C). In the case of multiple axons, the above argument retains its simplicity and we can write the expected time-averaged concentration of $M$ on axon $a$ as

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 ${\displaystyle {\sigma }_j<3}$ ms (Figure 1—figure supplement 4).

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

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

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

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

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

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

where $Q$, represents the single release amount/quantity, which is constant and independent of the current value of a global signal, $G(t)$. The peak of the response is happening at time $t_{\max }={\tau }_r\log \frac{{\tau }_r+{\tau }_d}{{\tau }_r}$, reaching the value ${\displaystyle G_{\mathrm{m}\mathrm{a}\mathrm{x}}(t)=\frac{Q}{{\tau }_d}{\left(\frac{{\tau }_r+{\tau }_d}{{\tau }_r}\right)}^{-\frac{{\tau }_r}{{\tau }_d}}}$. This form allows more general explorations (e.g., ${\tau }_r\approx 0$, yields the exponentially decaying $R(t)$), and most importantly it is the impulse response curve of the second-order linear system,

where $s(t)$ represents the signal that drives the system, and $\dot{G}(t)$ and $\ddot{G}(t)$ are the first and second time derivatives of $G(t)$.

In order to simplify extensive explorations of all OMP parameters, here we chose a simplified form for $R(t)$ that uses a single characteristic time of the OL spike response, ${\tau }_G={\tau }_r={\tau }_d$, yielding,

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

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

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

When the spikes on different axons follow the same renewal process, that is, obey the same $p_{\text{ISI}}$ 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 $t_{k_a}$ indicates the time of the $k^{th}$ spike on axon $a$, and the sum goes over all spikes that have occurred prior to time $t$. When Equation 9 is applied to the linear system in Equation 8 the analytical solution for $G(t)$ becomes the sum of the responses to spikes on all axons that it myelinates, i.e.,

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

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

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

For most neuronal plasticity models, saturation functions need to be introduced to stabilize the learning process. Similarly, in our implementation of OMP, we introduced two separate saturation functions for myelin addition, $F_s^A$, and removal, $F_s^R$. The OMP equation for the latency on axon $a$, ${\tau }_a$, can be written as

where $F_s^R(\tau )=H_r({\tau }_{\max }-\tau )$ and $F_s^A(\tau )=H_r(\tau -{\tau }_{\min })$ , $H_r(x)=xH(x)/({\tau }_{\max }-{\tau }_{\min })$ is the normalized ramp function, $H(x)$ is the Heaviside (unit step) function, ${\tau }_{\max }$ and ${\tau }_{\min }$ are the parameters of the model which specify the maximal and minimal delays that are attainable on any axonal connection, respectively.

We make two modifications to this basic model, one being the case of ${\lambda }_A\to \mathrm{\infty }$ (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.

In its fullest form, the OMP model can have up to 13 scalar parameters plus the $N_O\times N_A$ 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 ${\lambda }_R$, that now becomes an OMP variable. Omitting the equation for homeostasis requires fine-tuning the balancing value for ${\lambda }_R$ for a given parameter setting. This model can be simplified further by assuming ‘instantaneous myelination’, in which case Equation 11 can be omitted, when using Equation 15 instead of Equation 12; however, using Equation 15 requires, particularly if used together with Equation 17, that ${\lambda }_M$ be chosen sufficiently small to insure that the integration is stable. If ${\lambda }_M$ is large, sudden jumps in the value of ${\tau }_a$ will make Equation 17 stochastic. This can lead to negative values for ${\lambda }_R$ or even the delays themselves, which is not realistic.

The standard variables of the model, for a single OL, are $G(t)$, $M_a(t)$, and ${\tau }_a(t)$, $a=1\mathrm{\dots }{N}_A$, giving a total of $2N_A+1$ variables per OL, excluding ${\lambda }_R$. In practice, solving the model will require $2N_A+2$ variables per OL, since Equation 8 is a $2^{nd}$-order differential equation, requiring an auxiliary variable (see the next section). $G(t)$ and $M_a$ are ‘fast’ variables that change on a milliseconds time scale (dictated by ${\tau }_G$ and impulse responses to the spikes, with mean ISI, ${\tau }_s$), while, ${\lambda }_R$ and ${\tau }_a,a=1\mathrm{\dots }{N}_A$, are ‘slow’ variables whose rate of change is controlled by ${\lambda }_H$ and ${\lambda }_A$, respectively (${\lambda }_H\ll {\lambda }_A$). In Figure 1—figure supplement 2A, B, we show an example of time progression for both fast and slow variables, as well as the synchrony measure, ${\sigma }_{\tau }$ (Figure 1—figure supplement 2B, top row) and the epoch-averages of $M_a$ (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 $N_O=5$, $N_A=10$, and a single trial, except for ${\sigma }_{\tau }$ (dashed lines) and ${\lambda }_R$ (different colors), where the results from independent trials (3 total) are also shown. In Figure 4—figure supplement 3 we similarly show the time-progression for the unstable case, $N_O=1$, as shown in Figure 3A (${\displaystyle {\lambda }_H={10}^{-6}}$), but now for three different values of ${\lambda }_H$, indicating that the oscillations seen in ${\sigma }_{\tau }(t)$ are a result of homeostatic control. We note here, again, that $N_O=1$ does not mean literally that there is a single oligodendrocyte acting on a given axonal bundle but rather a population of oligodendrocytes acting at a particular location along the axon receiving the same pattern of activations and responding to it in the same way. When independent signals are conducted along OC, the time progression of the slow OMP variables displayed only stochastic variations and no clear trends, as shown in Figure 2—figure supplement 2 (apart from initial ‘refocusing’ of all ${\tau }_a$ to the same mean value, due to homeostatic constraints). This result was consistent across all simulations in which independent signals are used, with either Poisson or regular spiking $p_{\text{ISI}}$.

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

where ${\text{RSS}}_m$ is the residual sum of squares of model $m$. This $F$ is $F$-distributed, with ($n_U-n_R$, ${\displaystyle n_{\mathrm{p}\mathrm{t}\mathrm{s}}-n_U}$) 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, ${\displaystyle p<{10}^{-15}}$. This was not surprising, as the synchronization profiles were never truly statistically constant or exponential. For example, in Figure 3—figure supplement 2C, when ${\sigma }_{\tau }(t)$ is observed on a full scale, one would expect that the constant model is the best description of the behavior observed. However, the model selection with regular F-test chooses the less restricted models, in fact E2C2 in this case. After zooming in, we see that the results were indeed not constant, as those deviations from constancy were not just due to noise, but were significant. If the obtained simulation profiles were noisier then the restricted models would have a reasonable chance of being selected. Instead of adding artificial noise to our ${\sigma }_{\tau }(t)$, we introduced a parameter, $p_{\text{MSE}}$, with which we essentially control what level of noise or deviations we deem tolerable. The $p_{\text{MSE}}$ expresses this level of tolerance as the percentage of the initial, ${\sigma }_{\tau }$, or approximately, ${\sigma }_D$. Hence, we declare the minimal amount of RSS in any fit, ${\mathrm{RSS}}_{\min }=n_{\mathrm{pts}}*{\text{MSE}}_{\min }$ and ${\text{MSE}}_{\min }={(p_{\text{MSE}}*{\sigma }_{\tau }^{(0)}/100)}^2$. This sets the level of MSE that is presumed by default, that is, some minimal amount of noise present in residuals of any model. This is essentially specifying how much of RMS error can be tolerated in the restricted model, in order to reject the null hypothesis that the unrestricted model is better. Since the numerator will remain unchanged, this essentially only modifies the denominator,

where ${\mathrm{RSS}}_{\min }=n_{\mathrm{pts}}*{\text{MSE}}_{\min }$ and ${\text{MSE}}_{\min }={(p_{\text{MSE}}*{\sigma }_{\tau }^{(0)}/100)}^2$. The numerator portion is not affected, as it is a difference between two RSS. This yields the modified F-statistic, $F^m$,

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

In our procedure, we start with the model with minimal MSE (usually E2C2), test it against all of the restricted models, and choose the most restricted model for which the null hypothesis is not rejected. We performed the tests separately at different but very low significance levels (${\displaystyle \alpha \in [0.01,0.00001,{10}^{-10},{10}^{-15}])}$. The choice of the model was in many cases not strongly influenced by the choice of $\alpha$, but was strongly dependent on the choice $p_{\text{MSE}}$. For $p_{\text{MSE}}=0$ %, even at extremely low $\alpha$, the unrestricted model would always be chosen (see Figure 3—figure supplement 2C). Using $p_{\text{MSE}}=2$ %, allowed the choice of model to depend on $\alpha$ 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 $p_{\text{MSE}}=2$ %, based on a set of 100 random examples in which the best model is chosen manually by visual inspection (e.g., the constant model in Figure 3—figure supplement 2C). Note that our model selection is largely ad-hoc and we emphasize that our modified F-test does not aim to provide a quantitative statistical analysis, as use of such absurdly small significance levels indicates, but only to provide a useful quantitative tool for summarizing tens of thousands of runs that we have performed. While the distribution of different models changes significantly for different choices of $p_{\text{MSE}}$ and $\alpha$, the derived values of the parameters ${\sigma }_{\tau }^{\mathrm{\infty }}$ and ${\tau }_L$ is not significantly changed when different values of $p_{\text{MSE}}$ (but >1%) and $\alpha$ (but smaller than 0.01) were used. The same holds for the ad-hoc rule, of reverting to a less restricted model when ${\displaystyle {\text{MSE}}_R>500\times {\text{MSE}}_U}$, which happened very infrequently (see Figure 3—figure supplement 2D).

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: $N_A=10$, $N_O=5$, ${\tau }_G=\{10,20,30\}$ ms, ${\tau }_s=\{50,100,200\}$ ms, ${\lambda }_M=\{0.01,0.02,0.05,0.1\}$ ms^-1, ${\lambda }_A=\{0.1,0.01\}$ ms^-1, ${\lambda }_H={10}^{-6}$ ms^-2, $T_e=10$ sec, ${\displaystyle n_r=4}$ and $n_e=100$, ${\sigma }_D=\{5,10\}$ ms, ${\sigma }_j=\{1,3,5\}$ ms, ${\sigma }_s=0$ %, ${\tau }_{\max }=100$ ms, ${\tau }_{\min }=3$ ms, and ${\tau }_{\mathrm{nom}}=50$ ms, $t_R=0$ ms. In two cases the differences with the defaults were small, so we report them only here in the text. The first is studying the influence of the number of axons, $N_A$, for which ${\lambda }_A=0.01$, $N_O=\{5,10\}$, and $N_A=\{5,10,20,30,\mathrm{\dots },70,80,100,150\}$ ($n=4752$ runs). Runs with ${\sigma }_D=5$ ms ($n=2376$) are used in creating Figure 3F. The second explores the influence of ${\sigma }_D$ parameter, i.e., fixed delays, shown in Figure 4F, for which ${\sigma }_D=\{2,3,5,7,10,15,20\}$ ms and $n_r=5$ ($n=1512$). Below we show the tables for the other 8 sets, where the same units are used as above.

In the case of regular spiking, where the ISI times are distributed as ${\displaystyle p_{\text{ISI}}(t)=\sum _k\delta (t-{\tau }_s)}$, Equation 1 reduces to evaluating the sum,

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

which reduces to evaluating two geometric sums, yielding,

where $\delta t$ is the time to the most recent spike, and $G_a$ is replaced with $G_A$, indicating that the form of $G_A$ will be the same for different axons, but $\delta t$ on separate axons will differ. In other words, the contributions to the global $G(t)$ from different axons is only going to differ through the difference of their times to the most recent spikes, $\delta t_a$, i.e., ${\displaystyle G_a(\delta t_a)\equiv G_A(\delta t_a)}$. Assuming $t=0$ coincides with one of the spikes on a given axon we can write, $\delta t\equiv t(mod{\tau }_s)$. When combining $G_A(\delta t_a)$ signals from all axons, generally only one of them can be chosen as such reference, while $\delta t$ for others will be expressed in terms of fixed delays between them, and will require considering separately different orderings, depending on the fixed temporal delays between them. For example, in the simple case, $N_A=2$, for a given fixed delay, t_d, between two regular spiking trains (Figure 1—figure supplement 3A), we need to distinguish between the cases ${\displaystyle \delta t\le t_b}$ versus ${\displaystyle \delta t>t_b}$, where ${\displaystyle t_b={\tau }_s-t_d}$. 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 $\delta t_1$. When ${\displaystyle \delta t_1\le t_b}$,

and for ${\displaystyle \delta t>t_b}$,

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