# Activity-dependent myelination: A glial mechanism of oscillatory self-organization in large-scale brain networks

Edited by Marcus E. Raichle, Washington University in St. Louis, St. Louis, MO, and approved April 30, 2020 (received for review September 24, 2019)

## Significance

The structure and conductive properties of white matter play a key role in brain dynamics and cognitive processes such as memory. The coordination between different brain areas, essential for brain function, requires an adaptive traffic control system. This role is mediated by oligodendrocytes—glial cells which form myelin sheaths around axons to speed up saltatory conduction of nerve impulses across the brain. Whilst previously thought to be largely static after development, recent results show that white matter rewires itself in an activity-dependent way with experience and learning. Using a network model based on primate white matter data, our interdisciplinary approach reveals how activity-dependent myelination promotes neural phase synchronization, endowing white matter with self-organizing properties, where conduction delay statistics are autonomously adjusted to ensure efficient neural communication.

## Abstract

Communication and oscillatory synchrony between distributed neural populations are believed to play a key role in multiple cognitive and neural functions. These interactions are mediated by long-range myelinated axonal fiber bundles, collectively termed as white matter. While traditionally considered to be static after development, white matter properties have been shown to change in an activity-dependent way through learning and behavior—a phenomenon known as white matter plasticity. In the central nervous system, this plasticity stems from oligodendroglia, which form myelin sheaths to regulate the conduction of nerve impulses across the brain, hence critically impacting neural communication. We here shift the focus from neural to glial contribution to brain synchronization and examine the impact of adaptive, activity-dependent changes in conduction velocity on the large-scale phase synchronization of neural oscillators. Using a network model based on primate large-scale white matter neuroanatomy, our computational and mathematical results show that such plasticity endows white matter with self-organizing properties, where conduction delay statistics are autonomously adjusted to ensure efficient neural communication. Our analysis shows that this mechanism stabilizes oscillatory neural activity across a wide range of connectivity gain and frequency bands, making phase-locked states more resilient to damage as reflected by diffuse decreases in connectivity. Critically, our work suggests that adaptive myelination may be a mechanism that enables brain networks with a means of temporal self-organization, resilience, and homeostasis.

### Sign up for PNAS alerts.

Get alerts for new articles, or get an alert when an article is cited.

The rich repertoire of oscillatory activity observed in brain dynamics is thought to play a central role in neural communication and brain functions such as attention, large-scale integration, and memory (1–3). To implement these critical functions, distributed populations of neurons—behaving like nonlinear oscillators with dynamic frequencies and phases—need to coordinate in a flexible yet reliable fashion. Critical to this complex coordination task are long-range axons interconnecting spatially distant brain regions. In the brains of humans (and indeed, all higher organisms), these long-range, myelinated axons aggregate into large fiber bundles, collectively forming a dense and extensive set of tissues termed white matter (WM). Data about WM micro- and macrostructure have been used to understand recurrent patterns in the resting state (4–6), the emergence of neurological disorders (7), and even to inform personalized computational models (8, 9). Consistent across species and imaging modalities, studies have confirmed the fundamental role played by WM in the orchestration of oscillations involved in brain function (7, 8).

Like many other neurophysiological and neuroanatomical structures, WM properties are tuned to promote functionally relevant dynamics. This wiring is also flexible and changes through developmental stages and learning (10–12). Brain plasticity is driven by a combination of neurophysiological processes that unfold at different spatial and temporal scales. Synaptic plasticity, through which the effective strength of synaptic connections changes in time, plays a central role in adaptive brain wiring, but nonsynaptic processes are also involved. Notably, conductive properties of axons are also evolving in time—a feature that likely becomes even more significant as spatial distances increase and conduction delays become longer and more variable. While traditionally seen as a fixed template through which neural populations interact, recent experiments show that WM is, to the contrary, a dynamic and adaptive structure whose properties evolve across multiple temporal and spatial scales (11, 13–15). WM is not a static network of cables: through the continuous and dynamic action of glia, WM conductive properties change not only after development but well into adulthood.

Oligodendrocyte precursor cells are the largest population of dividing cells in the adult brain and generate myelin-forming oligodendrocytes that maintain and regulate the myelin sheath surrounding axonal membranes (16–20). Oligodendrogenesis is regulated by experience and by neural activity (10–15, 19, 21) and therefore, can adjust the local axonal conduction velocity (and thereby, conduction delays) of propagating action potentials (22–24). Myelin remodeling is not exclusive to WM and can impact gray matter as well (11, 13, 14). Through this process, axonal conduction velocities can span multiple orders of magnitude in the nervous system of mammals (25), allowing them to adaptively promote coincident signaling and efficient neural coding, despite highly variable axonal lengths (16, 18, 23, 26). Interestingly, this temporal regulatory mechanism fails across numerous conditions such as injury and/or demyelinating diseases (7, 19, 27–29)—leading to a loss of synchrony and/or oscillatory coherence between neural populations, resulting in hindered communication (29, 30).

Taken together, these results suggest that neuronal activity plays an important role in oligodendrogenesis and myelination (11, 13–15). Previous research has developed computational models of local myelin formation, notably by examining how node of Ranvier distributions impact axonal conduction velocities to influence action potential timing (31–33). However, the large-scale contributions of glia on brain dynamics are seldom considered in modeling studies. To address this gap, we combined whole-brain WM connectivity data with computational modeling to study the impact of adaptive myelination—specifically adaptive and phase-dependent changes in axonal conduction velocity—on self-organization of brain oscillations. Using a simple model of coupled neural oscillators (34, 35), we constructed a whole-brain network simulation that is augmented with plastic conduction velocities over a range of physiologically realistic values to reflect adaptive myelination (Fig. 1). Our simulations and mathematical analyses reveal that this addition of WM plasticity in the model confers enhanced synchronization properties. Under these conditions, our brain network model autonomously synchronized after a prolonged period of adaptation. Conduction velocity was found to compensate for variability in tract length, yielding highly similar conduction delays throughout the network. We further investigated the relative effect of neural and nonneural influences on network synchronization. While restricted phase locking did occur due to an increase in synaptic gain alone, adaptive myelination resulted in more robust synchronization across a broader range of frequency bands. Finally, our results show that adaptive myelination not only stabilizes coherent phase synchrony but makes it more resilient to sudden changes in connectivity resulting from diffuse injury or damage resulting from trauma. We further examine oscillatory dynamics resulting from bidirectional myelin remodeling. Our work provides insight about the regulatory role played by glia on brain oscillations, self-organization, and homeostasis.

Fig. 1.

## Results

### A Whole-Brain Network Model Built Using the Human WM Connectome.

As a first step to explore the impact of glia and adaptive myelination on brain synchrony, we built a brain-scale model of interacting oscillatory neural populations. Neural population dynamics were modeled by the well-known coupled Kuramoto phase oscillator equations (34, 36). This canonical model has been used as a stepping-stone to understand the mechanisms of synchronization across a wide range of complex systems in biology (35). Oscillatory units were placed at each node of a 96-node network (Fig. 2), composed of 80 cortical and 16 subcortical regions covering the entire human brain. Interregional connectivity weights within this network are based on macaque chemical tracer injections, and tract lengths are approximated by interregional Euclidean distances (

*Methods*has full details).Fig. 2.

In its simplest form, a given connection between two phase oscillators can be characterized by 1) how strong the connection is and 2) how fast signals travel along that connection. Thus, for each connection, two experimentally determined parameters were used to inform the model: 1) the synaptic gain (either zero or one representing the existence of a connection between two networks) and 2) the axonal tract length (representing the length of that connection and, for a given conduction speed, its associated conduction delay). To perform the simulations, we randomized the network phases within the model, and conduction velocities along each WM connection were set to an initial value of 3 m/s throughout the network. This velocity lies in the low to middle range of myelinated axons (37, 38). We initially focused on low-frequency neural oscillations (i.e., ${\omega}_{o}\approx 10\hspace{0.17em}\text{Hz}$) whose spectral features are known to display high sensitivity to WM structural properties (7, 8). However, we also examined other frequency bands (see below). Given the set of initial parameters and random phases, the network was first found in the incoherent state (36) (Fig. 3

*A*): despite continuous interactions between the different brain networks, no consistent phase relationship characterized the dynamics. The Kuramoto order parameter $\left(r\right)$, which is a standard metric of synchronization for this type of system (see*Methods*) and scales with the degree of synchronization across neural populations, was also found to be low across the oscillators in the network.Fig. 3.

### Inclusion of Adaptive Myelination Drives Spontaneous Synchronization in the Network.

Starting with the human connectome network model in the incoherent state, we then investigated the influence of gradual, activity-dependent changes in conduction velocities on self-organized brain synchronization. To do this, we enabled the model with a phenomenological myelination rule that increases conduction speeds based on the relative phase differences between connected neural populations. While simplistic, this glia plasticity rule mimics activity-dependent influences on myelin formation (see

*Methods*). We first considered the case of pure myelin accumulation over time where the conduction velocities can only increase (bidirectional myelin remodeling is discussed later). Through this phenomenological rule, changes in conduction velocities reflect adaptive myelination over the length of the axon, leading to conduction delays of various durations. Augmented with this adaptive myelination rule, and after a transient period of stabilization where the dynamics were also incoherent (i.e., not phase locked), we let the system evolve autonomously for a period of 10 h (3.6 × 10^{4}s) (Fig. 3*B*). During this time, conduction velocities increased slowly $\left({\alpha}_{c}\ll 1\right)$ as a function of the local phase offset between pairs of oscillatory neural populations. After this plasticity period of 10 h, we examined the distribution of conduction velocities over the network, as well as the Kuramoto order parameter $\left(r\right)$. The network was found to spontaneously synchronize: the value of the order parameter gradually increased until phase synchrony was achieved (Fig. 3*C*and*D*). To disambiguate the effects of the modeled WM plasticity mechanism from the generic tendency of a dynamic network to move toward a synchronous regime, we repeated the simulations without adaptive myelination by setting the myelination rate to zero (see*Methods*). Without the ability to adapt its conduction velocities, the network failed to synchronize: the order parameter remained close to zero throughout (i.e., in the incoherent state) (Fig. 3*A*). This indicates that plastic changes in axonal properties in the network cause phase locking through the self-regulation of its conduction times (i.e., delays), which became collectively shorter. However, this leads us to wonder: what is the statistical structure of the resulting connections after this period of adaptation?### Variability in Conduction Velocity Compensates for Variability in Tract Length.

If conduction velocity is constant and uniform across all network connections, conduction delays scale linearly with tract length (Fig. 4). Consequently, before plasticity occurs and where conduction velocity is the same throughout the network, conduction delays are highly variable. This can be seen in Fig. 3

*A*. However, this is not the case after the adaptation period (i.e., Fig. 3*B*), so how are conduction delays distributed? We investigated this question by considering the statistics of conduction velocities and delays across the network after an adaptation period of 10 h, to see how conduction delays (after the WM plasticity period) correlate with tract lengths. Our simulations indicate that there is a high degree of variability in conduction velocity—and thus, degrees of myelination—but also, a negative correlation between conduction velocity and tract length. As can be seen in Fig. 4*A*, high variability in the conduction velocity can be seen, as well as a statistical trend: longer connections tend to possess slower conduction velocities. This negative correlation is a direct consequence of the length-dependent metabolic drag coefficient (*Methods*) that tends to favor the myelinated connections with smaller tract lengths. Despite the variability in conduction velocities, unfolding the distribution of conduction delays according to tract length revealed an interesting bimodal configuration. Indeed, as seen in Fig. 4*B*, the distribution of conduction delays was dominated by two main modes: first, a linear distribution of delays corresponding to the connections that possess minimal myelination with the baseline conduction velocity (i.e., $\tau =l/{c}^{o}$). Another contribution could be observed, where statistically similar and small delays were found to be uniformly distributed, independently from tract length. Within this region, the model results suggest a degeneracy in tract lengths: multiple combinations of tract lengths and/or conduction velocities lead to conduction delays that are statistically similar. This means that within this region of parameter space, tract lengths are poor predictors of conduction delays and outline the multiplexing role played by glia. The statistics of these delays were not only aligned with our theoretical predictions (*Methods*) but also shed some light on the mechanism involved in the autonomous transition to phase-locked synchrony we observed in our simulations. Results presented in Fig. 4 show that synchrony is achieved in our model because of continuous changes in conduction delay statistics through a compensation of axonal tract length variability. Through this degeneracy, our results show that glia may, through adaptive myelination, lead to an effective decrease in conduction delay variability. Despite variable axonal distances, conduction speeds have been adjusted so that timing (i.e., delays) between two populations remains statistically the same, promoting phase locking. This dependence of synchrony on delay statistics has been observed in other biological networks, where delay distribution variance has been found to hinder phase locking of neural oscillators (39, 40). This consistency in conduction delays is furthermore in line with previous experimental observations where similar timings were observed across connections, irrespective of the variability of traveling distances (16–18). While this did not apply in our model to all network connections, their statistical weight indicates that this degeneracy in conduction delays represents an emergent feature of our model and was confirmed by theoretical calculations. However, what is this mechanism underlying spontaneous transition to synchrony, and how does adaptive myelination relate to other plasticity mechanisms?Fig. 4.

### Myelination Stabilizes Phase-Locked Solutions.

Adaptive synchrony in the brain results from a combination of physiological mechanisms, both neural and nonneural. However, how do they relate to each other? To disambiguate the contribution of local synaptic coupling and long-range adaptive myelination, we computed the Kuramoto order parameter $\left(r\right)$ as a function of varying synaptic coupling ($g$; scaling how strong local synaptic connections are) and oscillatory frequency band $\left({\omega}_{o}\right)$ for various rates of myelination. We found that the adaptive nature of our model stabilized phase-locked states across a wide range of synaptic strengths and frequency bands. As expected, the order parameter was found to increase as a function of synaptic coupling. In the absence of plasticity $\left(\mathit{\u03f5}=0\right)$, the synaptic strength required to achieve phase locking scales with the frequency band: stronger connections are required to stabilize higher-frequency phase synchrony (Fig. 5

*A*). We found that plastic conduction velocities increased significantly the repertoire of possible phase-locked states. Increasing the myelination rate (from $\u03f5=0.1\hspace{0.17em}\text{to 0.}2$) was found to enable autonomous synchronization of the system over much larger portions of parameter space. As seen in Fig. 5*B*and*C*, synchronization was spontaneously achieved over an extended range of frequency bands as well as weaker synaptic coupling gains. Interestingly, adaptive myelination was found to permit synchrony in regions of parameter space that were otherwise unable to do so. These regions, delimited in Fig. 4 for visibility, support synchronization not only for weaker synaptic gains $\left(g\right)$ but also, over a wide range of frequency bands.Fig. 5.

These results indicate that adaptive myelination, through the plastic modulation of conduction velocities, increases significantly the range of frequencies at which phase-locked synchrony can be observed and decreases the dependence of those states on the strength of synaptic coupling. Our results also indicate that in absence of this plasticity, and for weaker connection gains, the network loses its ability to self-organize and phase lock. This is in line with multiple experimental studies in both humans and animals (11, 12, 23, 29, 30), which have demonstrated the role of glia (specifically myelin formation and regeneration) in maintaining synchrony and/or phase coherence among neural populations.

### Resilience.

Our previous results indicate that our network model, enhanced with plastic conduction velocities to represent the impact of glial myelination, promotes synchrony over a larger portion of parameter space by compensating for weaker synaptic coupling and axonal length variability. We next asked, as a corollary, whether adaptive myelination could make the network more resilient to perturbations in connectivity—thus implementing a homeostatic mechanism preserving oscillatory coherence between neural populations subjected to structural changes (i.e., damage). To do this, we first considered strong synaptic gain, in which coherent synchrony can be achieved even in the absence of plasticity. As shown in Fig. 5, this occurs when both the synaptic gain is strong (frequency band is low) and when the conduction velocity is high (i.e., the conduction delays are less variable). We first fixed the conduction velocity for all connections to $3\hspace{0.17em}\mathrm{m}/\mathrm{s}$ and let the system evolve for 5 h (1.8 × 10

^{4}s). Then, 80% of existing connections were destroyed (i.e., connection probability of those connections was set to zero) (see*Methods*) after phase-locked synchronization was reached in the network, and the system was evolved for another 5 h. We then compared the dynamics exhibited by the network with and without plasticity (i.e., the dynamics where conduction velocities remain stable and those where conduction velocity can take any value within the defined range).Despite strong coupling, the damage inflicted by destroying network connections destabilized the phase-locked state (Fig. 6

*B*). The dynamics became incoherent, and the Kuramoto order parameter $\left(r\right)$ decreased dramatically. In the absence of plasticity, the network does not possess the ability to adapt (i.e., the statistics of the conduction velocities [and delays] remain unchanged), and the sudden decrease in network synaptic coupling due to damage destabilizes phase-locked synchrony. In contrast, plasticity allowed the system to maintain synchrony, albeit at lower levels than before damage was applied. After an adaptation period of 5 h, the order parameter recovered most of its amplitude (Fig. 6*B*). When exposed to this perturbation, the statistics of conduction velocities and delays undergo radical changes: they adapt to reserve synchronous phase locking over the remaining connections. As can be seen from Fig. 6*C*, many connections adjust their conduction velocity to compensate and decrease the delay to align the phases. We then explored how large the myelination rate had to be to preserve phase-locked synchrony for various degrees of insults. We gradually altered the network connectivity by changing the insult index and observed the robustness of the phase-locked state by measuring the order parameter $r$ in each case. An increasingly large myelination rate was required to preserve synchrony whenever damage severity increased (Fig. 6*D*). Taken together, these results show that glia, through the activity-dependent adjustment of conduction velocity, act like a homeostatic mechanism preserving phase-locked states in the presence of changes in network synaptic coupling and hence, damage.Fig. 6.

### Bidirectional Myelin Remodeling.

While the majority of the experimental evidence supports experience-dependent increases in myelination, there is emerging evidence for experience-dependent myelin retraction/thinning (e.g., refs. 41 and 42). Indeed, mature myelin sheath thickness and nodal gap length can be reversibly modulated via perinodal astrocytes, and this may represent one mechanism for these experience-dependent changes (43). Whether this process is regulated in an activity-dependent manner and on timescales similar to activity-regulated de novo myelination remains unclear. Nonetheless, we explored the potential contribution of bidirectional activity-regulated changes in oscillatory coherence. To do this, we modified the activity-dependent myelination rule to allow both myelin formation and retraction (see

*Methods*). Through this modification, illustrated in Fig. 7, the conduction velocity over a specific axonal tract is allowed to fluctuate bidirectionally as a function of phase, and can do so at different rates. With this modification, the ratio of the rate of retraction over the rate of stabilization (i.e., $R/S$) can be tuned to reflect a specific balance in myelin remodeling. Using this adjusted myelination rule, we repeated the simulations performed in Fig. 5 and computed the Kuramoto order parameter $\left(r\right)$ as a function of varying synaptic coupling $\left(g\right)$ and oscillatory frequency band $\left({\omega}_{o}\right)$ but this time for various retraction vs. stabilization ratios. We found that, even if myelin retraction is introduced, phase-locked states are still stabilized over a wide range of synaptic strengths and frequency bands. However, the scope of this synchronization was found to depend on the retraction vs. stabilization ratio: as retraction rates increased compared with stabilization, phase synchronization occurs over restricted regions of parameter space. This can be seen by comparing Fig. 7*B*and*C*with Fig. 5*C*: if retraction occurs at half the rate of stabilization (i.e., $R/S$ = 0.5), the adaptive synchronization property of the network is diminished. This is even more salient if retraction and stabilization would hypothetically occur at the same rate (i.e., $R/S=1.0$). These results can be understood through the lens of dynamical systems theory: the stronger myelin stabilization is, the higher the conduction velocity is and—as a corollary—the shorter the conduction delays. Favoring shorter delays results in enhanced and widespread synchronization, which is what is observed in Fig. 5*C*where there is no retraction (i.e., $R/S=0$).Fig. 7.

## Discussion

The active influence of glia represents one important aspect of brain plasticity that spans multiple spatial and temporal scales. However, the impact of adaptive myelination on brain synchronization remains poorly understood and seldom considered. To gauge the influence of myelination on large-scale brain dynamics, we built and analyzed a brain-scale network of neural oscillators, informed by whole-brain anatomical connectivity data, and enabled WM plasticity via a mechanism of phase-dependent modification of myelination level. Using this model, we characterized the effect of adaptive and time-dependent changes in conduction velocity on network synchronization. Our results reveal that phase-dependent myelination enables our brain-scale network model to establish a configuration of conduction delays that stabilizes synchronous phase-locked states across wider regions of parameter space. Variability in conduction velocities was found to compensate for tract length, yielding conduction delays that were statistically equivalent. As a corollary, activity-dependent myelination was found to make the network more resilient to spontaneous changes in network connectivity, preserving phase-locked synchrony in the presence of increasing, spatially diffuse damage. Even when bidirectional changes in conduction velocity were considered, oscillatory coherence was found to be maintained, indicating that myelin net accumulation is dominant. Taken together, our results indicate that the temporal structure of neural interactions, mediated by adaptive myelination (and the mechanisms that regulate it), is as important as synaptic strength in brain oscillatory self-organization and dynamics.

### Homeostasis.

Our results are revealing when considered from the perspective of brain homeostasis. Indeed, the simulations of Figs. 5 and 7 show that adaptive conduction velocities enable synchronization over an extended region of parameter space. Specifically, phase synchrony is facilitated for weaker network synaptic coupling and over a wider range of frequencies. This self-organizing feature of the model, mediated by glia, provides the opportunity to reorganize the timing of neural interactions to preserve specific network target dynamics (i.e., an attractor; here, phase-locked synchrony). When considered from the perspective of WM plasticity and larger-scale brain synchronization, this feature of our model represents an example of a dynamic implementation of homeostatic optimization—specifically, the self-organization and optimization of conduction velocities toward the preservation of a functionally relevant state. The need for such homeostatic optimization is further supported by experimental findings in mice. Indeed, experiments show that the absence of adaptive myelination contributes to methotrexate chemotherapy-related cognitive impairments in mice (20), which demonstrates the importance of adaptive myelination after brain injury in preserving cognitive function. In our model, this was achieved by the introduction of a glial-mediated phase-dependent myelination that acts as a trade-off between metabolic demands and local promotion of phase locking. Our results highlight that through this simple rule, it is possible to implement a form of oscillatory homeostatic control that corresponds, even in simplified systems like the one we studied here, to a form of parametric resilience (44, 45) implemented by glia.

Our results in Fig. 6 highlight the importance of adaptive myelination, and WM plasticity generally, in preserving functionally relevant dynamics in the presence of brain damage or injury. Experimental findings have shown that pathological states preventing adaptive myelination, for instance by injury-induced impaired glial reactivity, cause cognitive and behavioral impairment in animal models (e.g., ref. 20). From this perspective, our results thus suggest that adaptive myelination represents a critical therapeutic target to maintain brain function in the context of brain injury.

In any network with finite conduction velocity, delays are inevitable. The question thus arises as to what the benefit is of adaptive, activity-dependent myelination—as opposed to uniform and strong myelination (18, 22, 23, 26, 46). It would prove interesting to revisit our results from the perspective of large-scale metabolic and/or energetic constraints to see whether adaptive self-regulation of conduction velocities implements some form of metabolic optimization (e.g., ref. 22). Furthermore, the results in Fig. 7 also raise the question as to what the functional implications of bidirectional myelin remodeling would be. Interarea phase synchrony is thought to play a central role in neural communication and information routing (47). To be properly directed toward the right neural targets, oscillatory coherence must be supported by sets of axonal tracts and suppressed by others. From this perspective, synchrony is connectivity specific, and some form of architectural pruning should take place to implement selective functional relationships (i.e., synchrony) between related brain areas. Could adaptive myelination be used to prevent synchronization between functionally decoupled or unrelated regions? While there is insufficient experimental evidence to support these claims at the moment, they nonetheless suggest that WM modeling can in some cases promote, and some other cases prevent, phase synchronization.

Our model does not capture all features of glia and/or adaptive myelination processes and thus, remains an approximation meant to address large-scale consequences mediated by changes in conduction velocity configurations. It thus suffers from numerous limitations. First, the phenomenological myelination rule we consider here, despite time rescaling, occurs at a very fast pace compared with experimental observations. Indeed, the changes in conduction velocity in our model occur much faster than the myelin formation reported in animals, which stabilizes after a few hours (48). Nonetheless, our work attempts to bridge the timescale separation that characterizes the feedback loop between neural activity (milliseconds) and glia/myelination (hours) first by using a simplified neural oscillator model and also, by applying asymptotic (i.e., time-rescaling) analysis toward resulting conduction delay statistics. We further note that an asymptotic and/or adiabatic analysis of our model could provide key insight as to the net effect of activity on myelination. Second, the myelination model we implemented remains highly simplified. The different physiological constraints that dictate activity-dependent myelination in the brain are diverse, spatially distributed, and correlated (11, 18, 23, 24, 28, 49). To circumvent this, we have developed a myelination rule that possesses a local phase offset correcting term. This implies that a given connection “knows” the state of the neural oscillators it connects and tries to minimize phase misalignment between them. However, the physiological mechanisms involved in activity-dependent myelination in the brain remain poorly understood (11, 15, 20, 23, 49, 50). Third, the use of phase oscillators to model local neural dynamics, while well supported in the literature (ref. 35 and references therein), remains limited. The local synchronization features, and how they are impacted by afferent inputs, are likely to impact the myelination process and hence, greatly enhance the richness of the dynamics, as well as the complexity of the analysis.

Through a combination of experimental, computational, and mathematical analysis, our results highlight the importance of glia and adaptive myelination in large-scale brain synchronization. Our results show that through the plastic modulation of conduction velocities, adaptive myelination significantly increases the range of phase-locking frequencies, decreasing the dependence of those states on synaptic coupling. Our results further suggest that this plasticity reinforces the stability—and hence, resilience—of network synchrony to sudden and diffuse changes in connectivity. Taken together, our results emphasize the important role played by glia in implementing a homeostatic mechanism promoting phase-locked synchrony across large spatial distances and variable connectivity.

## Methods

### Human Whole-Brain Network Model.

To examine the influence of glia on brain synchronization, we revisited a well-studied network of coupled oscillators (34, 36) and built a brain-scale simulation of activity-dependent myelination. The Kuramoto model has been used as a stepping stone to understand the mechanisms of synchronization across a wide range of complex systems in biology (35, 36). The Kuramoto model has been used before to mimic localized oscillatory responses by focusing on phase-relevant dynamics. In neuroscience, it has proven to be a powerful tool to model macroscopic brain phenomena by quantifying the dynamic interactions between the relative phases of neural oscillators (51–54). Furthermore, its dynamics have been shown to capture synchronization-mediated interactions of a wide range of nonlinear coupled oscillators (55). Through this choice, our goal was to preserve a minimalistic neural model in order to focus on glial contributions, especially at larger spatial scales (Fig. 2).

According to the model below, oscillatory neural activity is generated locally, and interactions between neural populations occur via phase-dependent coupling. While our model remains agnostic as to the mechanism responsible for generating these local oscillations, we assume they are mediated by the interplay of local excitatory and inhibitory populations. We instead focus on how these oscillations organize through glia-supported connections at the macroscopic scale. Specifically, we considered a network of $N$ coupled nonlinear and delayed Kuramoto oscillators whose phases evolve according to the following set of equations (35, 39, 56):

$$\frac{d}{dt}{\theta}_{i}\left(t\right)=\hspace{0.17em}{\omega}_{o}+\hspace{0.17em}{N}^{-1}{\displaystyle \sum _{j=1}^{N}{w}_{ij}\hspace{0.17em}\mathrm{sin}\left[{\theta}_{j}\left(t-{\tau}_{ij}\left(t\right)\right)-{\theta}_{i}\left(t\right)\right]},$$

where ${\omega}_{o}$ represents a shared local frequency. Note that the delays are time dependent. The synaptic connectivity matrix ${w}_{ij}={\left[W\right]}_{ij}=g\hspace{0.17em}{p}_{ij}$ defines the connectivity of the network. The synaptic gain $g$ scales the strength of the connections, while the connection probability ${p}_{ij}=0$ or 1 refers to the existence or absence of a specific WM connection between populations $i$ and $j$. Note that through the combination of connection probability and synaptic gain, the effective weight of a connection (i.e., the net effect of a upstream network on a downstream network) is ${w}_{ij}=g{p}_{ij}$ and can be varied.

We have considered in the analysis the time evolution of the mean phase across the network: that is,

$$\overline{\theta}\left(t\right)={N}^{-1}{\displaystyle \sum _{i=1}^{N}{\theta}_{i}\left(t\right)}.$$

### Influence of Glia: Phase-Dependent Myelination.

Axonal conduction delays correspond to the time taken by an afferent signal to travel via saltatory conduction along a myelinated fiber, where the propagation velocity ${c}_{ij}$ changes along the length of the axon across successive nodes of Ranvier, and also changes in time. The net effective conduction delay over a given axonal connection may be approximated as the ratio of the axonal length over the mean conduction velocity: that is,

$${\tau}_{ij}\left(t\right)={\displaystyle \underset{A}{\int}{c}_{ij}^{-1}\left(\ell ,t\right)d\ell \approx \frac{{l}_{ij}}{{c}_{ij}\left(t\right)}},$$

where ${l}_{ij}={\left[L\right]}_{ij}$ is the total axonal tract length of the axonal path $A$ between populations $i$ and $j$.

We are interested in describing network effects of glia during a phase called “smart-wiring” (23, 26) in which WM is remodeled in an activity-dependent way. As such, we have enabled our human brain network model with a slow feedback loop linking oscillatory phase coupling and conduction velocity, which mediates a continuous and state-dependent modification of conduction delay statistics, thought to be implemented by glia. We consider here the following phenomenological phase-dependent myelination rule that governs the evolution of the conduction velocity ${c}_{ij}$,

$${\alpha}_{c}^{-1}\frac{d}{dt}{c}_{ij}=M\cdot \left(-A+{B}_{+}\right),$$

where ${c}_{ij}$ represents the conduction velocity along the axonal connection between neural populations $i$ and $j$. In this model, changes in conduction velocity mediated by glia are governed by the interplay between two general mechanisms: 1) an inertial drag toward a baseline conduction velocity (i.e., minimal myelination) to model the negative influence of metabolic demands on myelination ($A$) and 2) an activity-dependent myelination term that shapes conduction velocity in a phase-dependent way $\left({B}_{+}\right)$. Both mechanisms (

*A*and ${B}_{+}$) compete to regulate the local conduction velocity along the axonal connection between the nodes $i$ and $j$. An additional term, $M$, ensures that conduction velocity remains positive and bounded between 3 and 100 m/s [within realistic ranges for myelinated axons (e.g., refs. 25, 37, and 38)]. Myelin-induced changes in conduction velocities occur on a much slower timescale compared with neural oscillatory dynamics, and thus, the rate constant is chosen such that ${\alpha}_{c}\ll 1$. While this rule neglects a wide range of local neurophysiological processes underlying myelin formation by oligodendrocytes, its focus is on the macroscopic effect of myelination on axonal conduction (e.g., refs. 22, 49) and its influence on neural phase locking (23, 46).The inertial drift term is mathematically defined as

$$A={k}_{ij}\left({c}_{ij}\left(t\right)-{c}^{o}\right).$$

The baseline conduction velocity ${c}^{o}=3\hspace{0.17em}\text{m}/\text{s}$ has been set to represent unmyelinated and/or minimally myelinated axons. This baseline value has been chosen to match the lower bound for unmyelinated axons, to represent “minimal” myelination (25, 37, 38) (Table 1). The coefficient ${k}_{ij}$ represents the net effect of a combination of metabolic and/or biophysical forces preventing and/or slowing myelin formation. We have made this coefficient dependent on axonal tract length to represent the metabolic burden of myelinating longer axons: that is,

$${k}_{ij}={k}_{o}\frac{{l}_{ij}}{\mathrm{arg}\hspace{0.17em}\mathrm{max}{\left[L\right]}_{i,j}},$$

Table 1.

Symbol | Definition | Value |
---|---|---|

N | Network size | 96 |

g | Synaptic gain | Variable |

${\omega}_{o}$ | Local oscillation frequency | Variable |

${w}_{ij}$ | Synaptic weight | Variable |

${p}_{ij}$ | Connection probability | [0,1] |

${l}_{ij}$ | Connection length | Variable |

${\tau}_{ij}$ | Conduction delay | Variable |

${c}_{ij}$ | Conduction velocity | Variable |

${c}^{o}$ | Baseline conduction velocity | 3 m/s |

${c}_{max}$ | Maximal conduction velocity | 100 m/s |

dt | Integration time step | 1 ms |

$\gamma $ | Insult index | Variable |

${k}_{o}$ | Metabolic inertia constant | 0.01 arbitrary units |

${\alpha}_{c}$ | Plasticity rate constant | 0.001 s^{−1} |

$r$ | Order parameter | Variable |

$\mathit{\u03f5}$ | Myelination rate | Variable |

where ${k}_{o}$ is a scaling coefficient and $\mathrm{arg}\hspace{0.17em}\mathrm{max}{\left[L\right]}_{i,j}$ represents a normalization factor. Through this form, connections exhibiting longer spatial lengths will demyelinate faster, compared with shorter ones.

The activity-dependent myelination term has the form

$${B}_{+}=\{\begin{array}{cc}-\mathit{\u03f5}{p}_{ij}\mathrm{sin}{\mathrm{\Delta}}_{ij}& \text{if}\hspace{0.17em}2n\pi \le {\mathrm{\Delta}}_{ij}\le \left(2n+1\right)\pi \\ 0& \text{otherwise}\end{array}.$$

The myelination rate $\mathit{\u03f5}$ sets the gain or magnitude of the conduction velocity fluctuations and $n=\mathrm{0,1,2}\dots $. The connection probability ${w}_{ij}$, which we recall is either 0 or 1, is also present to ensure that only existing connections change (i.e., whenever the connection probability ${p}_{ij}>0$). The phase difference defined as ${\mathrm{\Delta}}_{ij}={\theta}_{j}-{\theta}_{i}$ represents a myelination-promoting forcing that depends on the phase difference between two oscillators $i$ and $j$.

The activity-dependent myelination rule ${B}_{+}$ only allows positive changes in conduction velocity. Given that myelin has been shown to accumulate over time (14), the rule above has been calculated so that that phase differences have a net positive influence on conduction velocity, or otherwise no effect (to reflect accumulation of myelin). A schematic of the phase-dependent myelination rule is plotted in Fig. 1, where changes to conduction velocity lead to phase locking of network oscillators. One can see that whenever the phase offset is negative, ${\theta}_{j}<{\theta}_{i}$, the phase lag between the two oscillators causes the conduction velocity to increase; myelin undergoes formation or stabilization. If the phase is the same between two oscillators, or the network has achieved phase-locked synchrony, then no changes in conduction velocity occur. Any further modification of the conduction velocity requires a nonzero phase difference, which can occur spontaneously or due to some input. This implies that fully phase-locked solutions (in which all phases are locked around one frequency and where no phase difference exists) would not emerge in our system unless the baseline conduction velocity ${c}^{o}$ is high enough. In the simulations, we set this baseline conduction velocity to a high value (3 m/s), in order not to bias our analysis and compare cases in which resulting delays become too large and unrealistic.

While experimental studies clearly support a net increase in activity-induced myelination, evidence of myelin retraction and/or thinning also exists (10, 43), although at a lower rate. Our computational model can provide predictions about the consequences of activity-dependent myelin thinning/demyelination on large-scale brain oscillatory coherence, notably bidirectional remodeling that would occur during learning (10), for instance. As such, in Fig. 7, we modified the activity-dependent myelination rule to reflect bidirectional (i.e., increase and decrease) phase-dependent changes in conduction velocity that occur at different rates. This modified rule is given by

$${B}_{\pm}=\{\begin{array}{cc}-\mathit{\u03f5}\hspace{0.17em}{p}_{ij}\mathrm{sin}{\mathrm{\Delta}}_{ij}& \text{if}\hspace{0.17em}2n\pi \le {\mathrm{\Delta}}_{ij}\le \left(2n+1\right)\pi \\ -\frac{R}{S}\mathit{\u03f5}\hspace{0.17em}{p}_{ij}\mathrm{sin}{\mathrm{\Delta}}_{ij}& \text{otherwise}\end{array},$$

where $0\le \frac{R}{S}\le 1$ reflects the retraction over stabilization ratio. If $\frac{R}{S}=0$, no retraction occurs, and we recover the case of unidirectional myelin accumulation (i.e., ${B}_{\pm}={B}_{+}$). In the hypothetical case where $R/S=1$, myelin retraction would occur at the same rate as stabilization. Other parameters in ${B}_{\pm}$ are identical to those used previously. In this bidirectional framework, whenever the phase offset ${\mathrm{\Delta}}_{ij}={\theta}_{j}-{\theta}_{i}$ is positive, the phase precession between two interconnected neural oscillators causes conduction velocity to decrease: myelin undergoes retraction and/or thinning. While myelin retraction and thinning occur at much smaller rates compared with formation and stabilization (57), they would nonetheless influence emerging oscillatory dynamics. We thus examined in Fig. 4 different stabilization/retraction $\left(R/S\right)$ ratios and their consequences on phase synchronization. We repeated the numerical experiments in Fig. 5 and examined synchronization properties as a function of stabilization/retraction ratios.

### Time Rescaling and Asymptotic Analysis.

The phenomenological network model above is a multiscale system: WM plasticity occurs on a much slower timescale compared with the dynamics of its oscillators. To circumvent this limitation for the purpose of numerical simulations, by combining those two scales in a common numerical framework, we have considered a time rescaling in which the time constant ${\alpha}_{c}$ is increased to match the dynamics of the neural oscillators. To validate this approach, and ensure that it does not alter the equilibrium of the system or the resulting delay and/or conduction velocity statistics, we analyzed the asymptotic dynamics of the myelination rule. Introducing the adiabatic approximation (i.e., setting ${{\alpha}_{c}}^{-1}\to 0$), if we consider ${c}^{o}<{c}_{ij}<{c}_{max}$,

$$0\approx -{k}_{ij}\left({c}_{ij}\left(t\right)-{c}^{o}\right)+{B}_{+}.$$

Solving for the conduction velocity yields, given that all $\mathit{\u03f5}$, ${w}_{ij}$, and ${k}_{ij}$ are positive,

$${c}_{ij}=\{\begin{array}{cc}{c}^{o}+\frac{\mathit{\u03f5}\hspace{0.17em}{p}_{ij}}{{k}_{ij}}\mathrm{sin}{\mathrm{\Delta}}_{ij}& \text{if}\hspace{0.17em}2n\pi \le {\mathrm{\Delta}}_{ij}\le \left(2n+1\right)\pi \\ {c}^{o}& \text{otherwise}\end{array}.$$

Given that we do not have access to the individual phase differences ${\mathrm{\Delta}}_{ij}$, and their short-time temporal fluctuations (which would be equivalent to solving analytically the full system), we can nonetheless assume that they are exponentially distributed (an approximation inspired by numerical observations): that is,

$$\rho \left(\mathrm{\Delta}\right)\approx \mathrm{exp}\left[-\left|\mathrm{\Delta}\right|\right].$$

Assuming that $\mathrm{\Delta}$ is small, we can approximate

$${c}_{ij}\approx \{\begin{array}{cc}{c}^{o}+\frac{\mathit{\u03f5}{p}_{ij}}{{k}_{ij}}{\displaystyle \underset{0}{\overset{\text{\pi}}{\int}}\mathrm{sin}\mathrm{\Delta}\rho \left(\mathrm{\Delta}\right)d\mathrm{\Delta}={c}^{o}+\left(\frac{1}{2}+\frac{{e}^{-\pi}}{2}\right)\left(\frac{\mathit{\u03f5}{p}_{ij}}{{k}_{ij}}\right)}& \text{if}\hspace{0.17em}0\le {\mathrm{\Delta}}_{ij}\le \pi \\ {c}^{o}& \text{otherwise}\end{array}.$$

According to this adiabatic and small phase difference approximation, by virtue of the symmetry of the $\mathrm{\Delta}$ distribution, half of the connections possess a baseline conduction velocity ${c}_{ij}={c}^{o}$ (i.e., ${\tau}_{ij}={l}_{ij}{c}_{o}^{-1}$), while the other half possess velocities that are distributed along the curve ${c}_{ij}={c}^{o}+\left(\frac{1}{2}+\frac{{e}^{-\pi}}{2}\right)\frac{\mathit{\u03f5}{p}_{ij}}{{k}_{ij}}$ [i.e., ${\tau}_{ij}={l}_{ij}{\left({c}^{o}+\left(\frac{1}{2}+\frac{{e}^{-\pi}}{2}\right)\frac{\mathit{\u03f5}{p}_{ij}}{{k}_{ij}}\right)}^{-1}$]. Delays resulting from these approximations are plotted in Fig. 4

*A*, dashed line alongside the simulated data generated with the time rescaling. One can see that this analysis provides a good approximation of the two main statistical trends of conduction speed distributions, and hence that the asymptotic properties of the system are preserved.### Order Parameter.

To measure the degree of phase synchronization, we used the well-established Kuramoto order parameter $r$ (34, 36) that quantifies the degree of phase locking across a population of $N$ oscillators. The order parameter is equivalent to the phase-locking value or phase-locking index used to quantify phase coherence between neural time series (29, 30):

$$r\left(t\right)={N}^{-1}{\displaystyle \sum _{j=1}^{N}{e}^{i{\theta}_{j}\left(t\right)}}.$$

This parameter, whose values are between zero (fully incoherent phases) and one (full synchronization), fluctuates in time and scales with the degree of phase locking in the network.

### Damage.

To model damage and test the resilience of our network, we introduced the insult index $\gamma $, which takes values between zero (no damage) and one (fully disconnected). We considered $\gamma $ as a disconnection probability and modified elements of the matrix $W$ according to the following rule:

$${\stackrel{\sim}{w}}_{ij}=\{\begin{array}{c}{w}_{ij}|p>\gamma \\ \hspace{0.17em}0\hspace{0.17em}|p\le \gamma \end{array},$$

where $p=\left[\mathrm{0,1}\right]$ is a uniform random deviate. Collectively, the application of this rule sets a total of $\gamma N$ connection weights to zero.

### Connectivity and Tract Lengths.

Our model represents long-range WM pathways interconnecting distal gray matter regions using two adjacency matrices, one describing anatomical connectivity and the other average axonal conduction delays between pairs of gray matter regions. The delays matrix is, in turn, derived from a fixed matrix of tract lengths and an adjustable per-connection conduction velocity.

The weights and length matrices were taken from the 96-node connectivity dataset (“connectivity96”) distributed with The Virtual Brain neuroinformatics and brain simulation platform (https://github.com/the-virtual-brain). This connectivity is a subset of the Collation of Connectivity data on the Macaque brain (CoCoMac), a database containing chemical tract tracer measurements of the macaque brain (58). The weights are mapped from the macaque brain onto the human neuroanatomy using a regional map (RM) parcellation scheme (59). The RM scheme defines brain areas topographically and is suitable for mapping data between brains of related species with similar parcellation maps. The weights in the network are denoted as connection probabilities ${p}_{ij}$ and span 80 cortical and 16 subcortical regions.

Because this connectivity is adapted from macaque tracing data, the physical length of the corresponding fiber tracts in the human brain is not available from the CoCoMac database. Thus, the matrix of tract lengths is approximated using Euclidean distances between coordinate centroids of each gray matter parcel of the human brain. Refs. 60 and 61 have a full description of this connectivity and recent improvements to it using diffusion-weighted MRI.

### Data Availability

The anatomical connectivity data used in this study is distributed freely with the open source neuroinformatics and brain simulation software library The Virtual Brain (https://github.com/the-virtual-brain). Information on reproducing the analysis in this study, including the code used for generating simulations, figures, and data, is fully accessible at https://github.com/Lefebvrelab/ActivityDependentMyelinationModel.

## Data Availability

Data deposition: C++ and Python source code and documentation for the numerical simulations presented here are freely available in GitHub at https://github.com/lefebvrelab/ActivityDependentMyelinationModel.

## Acknowledgments

We thank M. Chantigny for computing support. We also thank National Research Council of Canada Grant RGPIN-2017-06662 and Canadian Institute for Health Research Grant NO PJT-156164 for funding.

## References

1

F. Varela, J.-P. Lachaux, E. Rodriguez, J. Martinerie, The brainweb: Phase synchronization and large-scale integration.

*Nat. Rev. Neurosci.***2**, 229–239 (2001).2

X.-J. Wang, Neurophysiological and computational principles of cortical rhythms in cognition.

*Physiol. Rev.***90**, 1195–1268 (2010).3

M. Siegel, T. H. Donner, A. K. Engel, Spectral fingerprints of large-scale neuronal interactions.

*Nat. Rev. Neurosci.***13**, 121–134 (2012).4

G. Deco, V. Jirsa, A. R. McIntosh, O. Sporns, R. Kötter, Key role of coupling, delay, and noise in resting brain fluctuations.

*Proc. Natl. Acad. Sci. U.S.A.***106**, 10302–10307 (2009).5

R. Hindriks et al., Role of white-matter pathways in coordinating alpha oscillations in resting visual cortex.

*Neuroimage***106**, 328–339 (2015).6

S. Atasoy, I. Donnelly, J. Pearson, Human brain networks function in connectome-specific harmonic waves.

*Nat. Commun.***7**, 10340 (2016).7

P. L. Nunez, R. Srinivasan, Neocortical dynamics due to axon propagation delays in cortico-cortical fibers: EEG traveling and standing waves with implications for top-down influences on local networks and white matter disease.

*Brain Res.***1542**, 138–166 (2014).8

J. Cabral et al., Exploring mechanisms of spontaneous functional connectivity in MEG: How delayed network interactions lead to structured amplitude envelopes of band-pass filtered oscillations.

*Neuroimage***90**, 423–435 (2014).9

G. Deco et al., Perturbation of whole-brain dynamics in silico reveals mechanistic differences between brain states.

*Neuroimage***169**, 46–56 (2018).10

J. Scholz, M. C. Klein, T. E. J. Behrens, H. Johansen-Berg, Training induces changes in white-matter architecture.

*Nat. Neurosci.***12**, 1370–1371 (2009).11

E. M. Gibson et al., Neuronal activity promotes oligodendrogenesis and adaptive myelination in the mammalian brain.

*Science***344**, 1252304 (2014).12

P. E. Steadman et al., Disruption of oligodendrogenesis impairs memory consolidation in adult mice.

*Neuron***105**, 150–164.e6 (2020).13

S. Mitew et al., Pharmacogenetic stimulation of neuronal activity increases myelination in an axon-specific manner.

*Nat. Commun.***9**, 306 (2018).14

E. G. Hughes, J. L. Orthmann-Murphy, A. J. Langseth, D. E. Bergles, Myelin remodeling through experience-dependent oligodendrogenesis in the adult somatosensory cortex.

*Nat. Neurosci.***21**, 696–706 (2018).15

M. Swire, Y. Kotelevtsev, D. J. Webb, D. A. Lyons, C. Ffrench-Constant, Endothelin signalling mediates experience-dependent myelination in the CNS.

*eLife***8**, e49493 (2019).16

M. Salami, C. Itami, T. Tsumoto, F. Kimura, Change of conduction velocity by regional myelination yields constant latency irrespective of distance between thalamus and cortex.

*Proc. Natl. Acad. Sci. U.S.A.***100**, 6174–6179 (2003).17

R. Vicente, L. L. Gollo, C. R. Mirasso, I. Fischer, G. Pipa, Dynamical relaying can yield zero time lag neuronal synchrony despite long conduction delays.

*Proc. Natl. Acad. Sci. U.S.A.***105**, 17157–17162 (2008).18

A. H. Seidl, Regulation of conduction time along axons.

*Neuroscience***276**, 126–134 (2014).19

M. S. Kaller, A. Lazari, C. Blanco-Duque, C. Sampaio-Baptista, H. Johansen-Berg, Myelin plasticity and behaviour-connecting the dots.

*Curr. Opin. Neurobiol.***47**, 86–92 (2017).20

A. C. Geraghty et al., Loss of adaptive myelination contributes to methotrexate chemotherapy-related cognitive impairment.

*Neuron***103**, 250–265.e8 (2019).21

I. A. McKenzie et al., Motor skill learning requires active central myelination.

*Science***346**, 318–322 (2014).22

T. Chomiak, B. Hu, What is the optimal value of the g-ratio for myelinated fibers in the rat CNS? A theoretical approach.

*PLoS One***4**, e7754 (2009).23

R. D. Fields, A new mechanism of nervous system plasticity: Activity-dependent myelination.

*Nat. Rev. Neurosci.***16**, 756–767 (2015).24

A. Klingseisen, D. A. Lyons, Axonal regulation of central nervous system myelination: Structure and function.

*Neuroscientist***24**, 7–21 (2018).25

H. A. Swadlow, J. D. Kocsis, S. G. Waxman, Modulation of impulse conduction along the axonal tree.

*Annu. Rev. Biophys. Bioeng.***9**, 143–179 (1980).26

M. E. Bechler, M. Swire, C. Ffrench-Constant, Intrinsic and adaptive myelination-A sequential mechanism for smart wiring in the brain.

*Dev. Neurobiol.***78**, 68–79 (2018).27

P. J. Uhlhaas, W. Singer, The development of neural synchrony and large-scale cortical networks during adolescence: Relevance for the pathophysiology of schizophrenia and neurodevelopmental hypothesis.

*Schizophr. Bull.***37**, 514–523 (2011).28

K.-A. Nave, H. Ehrenreich, Myelination and oligodendrocyte functions in psychiatric diseases.

*JAMA Psychiatry***71**, 582–584 (2014).29

S. Bells et al., Changes in white matter microstructure impact cognition by disrupting the ability of neural assemblies to synchronize.

*J. Neurosci.***37**, 8227–8238 (2017).30

S. Bells et al., White matter plasticity and maturation in human cognition.

*Glia***67**, 2020–2037 (2019).31

I. L. Arancibia-Cárcamo et al., Node of Ranvier length as a potential regulator of myelinated axon conduction speed.

*eLife***6**, e23329 (2017).32

M. C. Ford et al., Tuning of Ranvier node and internode properties in myelinated axons to adjust action potential timing.

*Nat. Commun.***6**, 8073 (2015).33

Y. Bakiri, R. Káradóttir, L. Cossell, D. Attwell, Morphological and electrical properties of oligodendrocytes in the white matter of the corpus callosum and cerebellum.

*J. Physiol.***589**, 559–573 (2011).34

Y. Kuramoto,

*Chemical Oscillations, Waves, and Turbulence*, (Springer Series in Synergetics, Springer-Verlag, Berlin, Germany, 1984).35

M. Breakspear, S. Heitmann, A. Daffertshofer, Generative models of cortical oscillations: Neurobiological implications of the Kuramoto model.

*Front. Hum. Neurosci.***4**, 190 (2010).36

S. H. Strogatz, From Kuramoto to Crawford: Exploring the onset of synchronization in populations of coupled oscillators.

*Physica D***143**, 1–20 (2000).37

F. Aboitiz, A. B. Scheibel, R. S. Fisher, E. Zaidel, Fiber composition of the human corpus callosum.

*Brain Res.***598**, 143–153 (1992).38

F. Aboitiz, D. Morales, J. Montiel, The evolutionary origin of the mammalian isocortex: Towards an integrated developmental and functional approach.

*Behav. Brain Sci.***26**, 535–552 (2003).39

G. B. Ermentrout, Oscillator death in populations of “all to all” coupled nonlinear oscillators.

*Physica D***41**, 219–231 (1990).40

F. M. Atay, Distributed delays facilitate amplitude death of coupled oscillators.

*Phys. Rev. Lett.***91**, 94101 (2003).41

J. Liu et al., Impaired adult myelination in the prefrontal cortex of socially isolated mice.

*Nat. Neurosci.***15**, 1621–1623 (2012).42

M. Makinodan, K. M. Rosen, S. Ito, G. Corfas, A critical period for social experience-dependent oligodendrocyte maturation and myelination.

*Science***337**, 1357–1360 (2012).43

D. J. Dutta et al., Regulation of myelin structure and conduction velocity by perinodal astrocytes.

*Proc. Natl. Acad. Sci. U.S.A.***115**, 11832–11837 (2018).44

E. Marder, J.-M. Goaillard, Variability, compensation and homeostasis in neuron and network function.

*Nat. Rev. Neurosci.***7**, 563–574 (2006).45

E. Marder, Variability, compensation, and modulation in neurons and circuits.

*Proc. Natl. Acad. Sci. U.S.A.***108**(suppl. 3), 15542–15548 (2011).46

S. Pajevic, P. J. Basser, R. D. Fields, Role of myelin plasticity in oscillations and synchrony of neuronal activity.

*Neuroscience***276**, 135–147 (2014).47

P. Fries, A mechanism for cognitive dynamics: Neuronal communication through neuronal coherence.

*Trends Cogn. Sci.***9**, 474–480 (2005).48

L. Xiao et al., Rapid production of new oligodendrocytes is required in the earliest stages of motor-skill learning.

*Nat. Neurosci.***19**, 1210–1217 (2016).49

R. G. Almeida, D. A. Lyons, On myelinated axon plasticity and neuronal circuit formation and function.

*J. Neurosci.***37**, 10023–10034 (2017).50

C. Demerens et al., Induction of myelination in the central nervous system by electrical activity.

*Proc. Natl. Acad. Sci. U.S.A.***93**, 9887–9892 (1996).51

J. Cabral, E. Hugues, O. Sporns, G. Deco, Role of local network oscillations in resting-state functional connectivity.

*Neuroimage***57**, 130–139 (2011).52

J. Cabral, E. Hugues, M. L. Kringelbach, G. Deco, Modeling the outcome of structural disconnection on resting-state functional connectivity.

*Neuroimage***62**, 1342–1353 (2012).53

A. Ponce-Alvarez et al., Resting-state temporal synchronization networks emerge from connectivity topology and heterogeneity.

*PLoS Comput. Biol.***11**, e1004100 (2015).54

S. Petkoski, J. M. Palva, V. K. Jirsa, Phase-lags in large scale brain synchronization: Methodological considerations and in-silico analysis.

*PLoS Comput. Biol.***14**, e1006160 (2018).55

S. Petkoski, V. K. Jirsa, Transmission time delays organize the brain network synchronization.

*Philos. Trans. Royal Soc. Math. Phys. Eng. Sci.***377**, 20180132 (2019).56

M. K. S. Yeung, S. Strogatz, Time delay in the kuramoto model of coupled oscillators.

*Phys. Rev. Lett.***82**, 648–651 (1999).57

F. Auer, S. Vagionitis, T. Czopka, Evidence for myelin sheath remodeling in the CNS revealed by

*in vivo*imaging.*Curr. Biol.***28**, 549–559.e3 (2018).58

K. E. Stephan et al., Advanced database methodology for the collation of connectivity data on the macaque brain (CoCoMac).

*Philos. Trans. R. Soc. Lond. B Biol. Sci.***356**, 1159–1186 (2001).59

R. Kötter, E. Wanke, Mapping brains without coordinates.

*Philos. Trans. R. Soc. Lond. B Biol. Sci.***360**, 751–766 (2005).60

G. Bezgin, A. Solodkin, R. Bakker, P. Ritter, A. R. McIntosh, Mapping complementary features of cross-species structural connectivity to construct realistic “Virtual Brains.”.

*Hum. Brain Mapp.***38**, 2080–2093 (2017).61

K. Shen et al., A macaque connectome for large-scale network simulations in TheVirtualBrain.

*Sci. Data***6**, 123 (2019).## Information & Authors

### Information

#### Published in

#### Classifications

#### Copyright

© 2020. Published under the PNAS license.

#### Data Availability

Data deposition: C++ and Python source code and documentation for the numerical simulations presented here are freely available in GitHub at https://github.com/lefebvrelab/ActivityDependentMyelinationModel.

#### Submission history

**Published online**: June 1, 2020

**Published in issue**: June 16, 2020

#### Keywords

#### Acknowledgments

We thank M. Chantigny for computing support. We also thank National Research Council of Canada Grant RGPIN-2017-06662 and Canadian Institute for Health Research Grant NO PJT-156164 for funding.

#### Notes

This article is a PNAS Direct Submission.

### Authors

#### Competing Interests

The authors declare no competing interest.

## Metrics & Citations

### Metrics

#### Citation statements

#### Altmetrics

### Citations

If you have the appropriate software installed, you can download article citation data to the citation manager of your choice. Simply select your manager software from the list below and click Download.

#### Cited by

Loading...

## View Options

### View options

#### PDF format

Download this article as a PDF file

DOWNLOAD PDF### Get Access

#### Login options

Check if you have access through your login credentials or your institution to get full access on this article.

Personal login Institutional Login#### Recommend to a librarian

Recommend PNAS to a Librarian#### Purchase options

Purchase this article to get full access to it.