Synchronization and resilience in the Kuramoto white matter network model with adaptive state-dependent delays.

White matter pathways form a complex network of myelinated axons that regulate signal transmission in the nervous system and play a key role in behaviour and cognition. Recent evidence reveals that white matter networks are adaptive and that myelin remodels itself in an activity-dependent way, during both developmental stages and later on through behaviour and learning. As a result, axonal conduction delays continuously adjust in order to regulate the timing of neural signals propagating between different brain areas. This delay plasticity mechanism has yet to be integrated in computational neural models, where conduction delays are oftentimes constant or simply ignored. As a first approach to adaptive white matter remodeling, we modified the canonical Kuramoto model by enabling all connections with adaptive, phase-dependent delays. We analyzed the equilibria and stability of this system, and applied our results to two-oscillator and large-dimensional networks. Our joint mathematical and numerical analysis demonstrates that plastic delays act as a stabilizing mechanism promoting the network's ability to maintain synchronous activity. Our work also shows that global synchronization is more resilient to perturbations and injury towards network architecture. Our results provide key insights about the analysis and potential significance of activity-dependent myelination in large-scale brain synchrony.

Introduction

Synchronization, the mechanism by which oscillatory processes collectively organize to align their phase, has been the focus of intense research across the field of non-linear dynamics [1], especially in the biological sciences due to its involvement in myriads of physiological processes. In the brain, such synchronized oscillatory patterns, observed in the rhythmic discharge of neuronal electrical impulses, play a central role in neural communication, information processing, and the implementation of higher cognitive function [2]. However, the mechanisms by which these oscillations emerge and interact within and across brain microcircuits and systems remain poorly understood.

In the mammalian brain, local synchronized oscillatory activity are coordinated through a vast network of myelinated axonal fibers called white matter. The intricate white matter structure is formed under a population of glial cells called oligodendrocytes. The oligodendrocytes produce an insulating myelin sheath coiling around axonal membranes to greatly increase the conduction velocities of propagating signals between neurons. As white matter develops to adopt a genetically programmed structure, oligodendrocytes determine to which extent specific axons are myelinated. The resulting myelin structure of the network is responsible for the emergence and evolution of a rich repertoire of spatiotemporal activity patterns, notably oscillations with various spectral properties [3].

While white matter is highly relevant in shaping brain network dynamics, the mechanisms governing myelin development and its relationship with neural activity remain mostly unknown [4]. Nevertheless, the general hypothesis surrounding this topic has been shifting away from traditional viewpoints. Recent studies suggest that white matter structure undergoes continuous changes past the stages of developmental myelination and well into adulthood, rather than remaining static as initially presumed [5]. Indeed, substantial myelin formation continues to occur within the fully mature central nervous system in an activity-dependent manner [6, 7]. Furthermore, newly discovered evidence implies that white matter structure is responsive to experiences such as learning and social interactions [8]. These findings unveil potential ways white matter can restructure itself to facilitate neurological function over time. In particular, white matter characterized as a plastic, adaptive structure, can be critical in maintaining the essential oscillatory and synchronous processes found in many neural systems [9].

Despite the complexity of the neurophysiological processes involved, it has been hypothesized that myelin remodeling collectively reinforces synchronous dynamics and oscillatory coordination in large-scale brain networks [8]. From this perspective, the temporal structure of these plastic changes can provide a higher level of corrective adjustments, improving the network’s ability to converge towards phase synchrony. As a first approach to this intricate problem, we wish to use a simplified model to determine whether activity-dependent delays impact the phase coordination of coupled oscillators when compared to static delays. We have used the Kuramoto model, a phenomenological model which has a long history of applications in neuroscience, notably to study the collective organization of oscillatory neural systems [10]. The Kuramoto model, with or without time delay, has been extensively studied [1, 11, 12], and is increasingly used to model local oscillatory neural activity in brain-scale simulations informed by anatomical white matter connectivity [13-15]. Appealing for its relative simplicity, we have enabled the Kuramoto model with phase-dependent delays to examine the collective impact of adaptive conduction on coupled phase oscillators. Specifically, we performed the stability analysis of this modified system and examined its response to structural perturbations. While preliminary, these results can provide new insights into the large-scale impact of white matter remodeling and its potential role in neural synchrony.

This paper is structured as follows. In Sect. 2, we first introduce our network model composed as a system of coupled phase-oscillators with state-dependent delays. In Sect. 3, we derive the equations for an N-oscillator network’s synchronous state and its respective stability criteria. Section 4 applies the derived equations from Sect. 3 to a reduced two-oscillator network. Section 5 applies the derived equations from Sect. 3 to a large-dimensional oscillator network through taking an N-limit approximation by handling the asymptotic phases statistically. In Sect. 6, we conduct numerical experiments and provide evidence supporting the analysis in Sects. 4 and 5 through simulated results. Section 7 explores the stabilization property in the context of spontaneous changes in network connectivity and compares the resilience of the synchronized state with and without plastic delays.

Model

We consider a prototypical model involving a network of N weakly-coupled non-linear and delayed Kuramoto oscillators [9] whose phases evolve according to the following system of differential equations: where is the natural frequency of oscillator i, and is a constant global coupling gain. The coefficients represent synaptic weights: if there is a white matter connection between oscillator j to i; otherwise, . The axonal conduction delays correspond to the conduction time taken by an action potential along a myelinated fiber. The propagation speed fluctuates with respect to a propagating signal’s position ℓ along the length of an axon as saltatory conduction occurs at successive nodes of Ranvier. The conduction velocity also changes in time. As such, each conduction delay may be written as That is, a line integral along the axonal pathway connecting oscillator j to i. It is uncertain how neural activity and/or oscillatory brain dynamics influences myelin formation. However, it is known that myelination correlates with information transmission within and across brain areas [5], which are known to involve long range synchronization [16]. To represent this in our model, we echo previous work on oscillatory neural communication [8, 9, 17] and make axonal conduction delays phase-dependent. Specifically, we assume that each conduction delay evolves under the following sublinear plasticity equation: where , . Here, is the homeostatic rate constant that sets the time scale of the evolving delays. The plasticity coefficient sets the gain of the conduction delay changes. The initial condition for delays at is provided by baseline lags . That is, for all i, j. The Heaviside function is defined as a smooth function such that if and if for some small . The Heaviside function is present to preserve non-negativity of the delays at all times . For details on the construction of the smooth Heaviside function, we refer the reader to Appendix A. According to this rule, fluctuations in conduction delays are governed by an interplay between a local drift component that represents metabolic inertia towards an initial baseline lag representing minimal myelination, and a forcing term that depends on the phase difference between oscillators i and j.

A schematic of the network structure as well as the activity-dependent plasticity rule are plotted in Fig. 1. There, given , one can see that whenever , the time delay increases due to an effective reduction in conduction velocity caused by myelin retraction. The opposite occurs if and the time delay decreases: the connection speeds up due to myelin stabilization [8]. If the phase is the same, no changes in conduction velocity are required, and the delay remains stable at its initial lag .

Figure 1

Schematic of a network of coupled oscillators with plastic conduction delays implementing temporal control of global synchronization. The network is built of recurrently connected phase oscillators, and connections are subjected to conduction delays. Those delays are regulated by a phase-dependent plasticity rule. Depending on whether the phase difference between two-oscillators is positive, zero, or negative, the delays are either increased (myelin retraction), unaltered, or decreased (stabilization)

Synchronized state and stability with plastic delays

We are interested in characterizing the influence of the delay plasticity rule Eq. (3) towards the stability of global phase-locked solutions and, more generally, in stabilizing synchrony between neural populations (i.e. oscillators). Enabled with adaptive delays, our model corresponds to an -dimensional functional differential equation with state-dependent delays. The analysis of such systems is technically challenging, especially in terms of stability where a modified approach must be used [18-20]. Mathematically, we analyze the network’s ability to asymptotically achieve the following synchronized state: for all as , where Ω is the global fixed frequency of the oscillators and is the asymptotic phase-offset of oscillator i as time . As we will see, adaptive delays following the plasticity rule equation (3) require non-zero phase differences in order to maintain its equilibrium values, from which the network can no longer become in-phase. Nevertheless, the oscillators are able to become entrained to a common frequency Ω. Hence, we say that our network becomes synchronized when the nodes globally oscillate at some stable frequency Ω and become phase-locked under individual offsets . As we are primarily concerned with the effects of delays changes under the plastic Eq. (3), we assume that all oscillators share the same natural frequency for all i in Eq. (1). To simplify the convergent behaviour of delays following Eq. (3), we set .

Applying the ansatz Eq. (4) onto both the phases in Eq. (1) and the delays in Eq. (3), we obtain the following expressions for the global frequency and individual offsets : for all i, j, where . It must be true that if Ω satisfies Eq. (5), then . To assess the network’s stability at the synchronized state satisfying Eqs. (5) and (6), we introduce the perturbation terms , around the equilibrium states Eq. (4) and Eq. (6) for and , respectively. That is, we write and examine the stability at the equilibrium point . The delay perturbations abide by the linearized form of Eq. (3) around all positive delays . For the equilibrium delays , we proceed by assuming the corresponding perturbation terms will act in an asymptotically stable manner. Specifically, for all phase-offset differences such that , by the nature of the Heaviside cutoff function in Eq. (3) there exists some such that for all , and consequently . That is, if , the perturbation term is asymptotically . This condition will be satisfied given the synchronous state is indeed stable. Since the purpose of linearization is to determine stability, the above assumption is valid to make before proceeding. Taken together, the linearized equations for the terms as are given by where . We can use the linearization approach in the context of equations with state-dependent delays [18-20] as follows. Inserting Eq. (7) and Eq. (8) into the Kuramoto equation (1), where we have made the approximation when is near . With satisfying Eq. (5), and by taking a first-order expansion around the term , we obtain the linear system for the phase perturbation term given by where we denote and . Hence, the network synchronizes at given that all perturbation terms , following the linear system of fixed delay differential Eqs. (9) and (11) converges to 0. From here, we can analyze the stability at the synchronous point by setting the ansatz , with respect to eigenvalue . By Eq. (9), the coefficients satisfy for all i, j, . Applying Eq. (12) to the coefficients , That is, λ is an eigenvalue if there exists an eigenvector satisfying Eq. (13). In matrix form, Eq. (13) can be expressed as , where is the N-dimensional square matrix with entries In summary, we acquire the following stability criterion: The system is stable around the synchronized state if for all eigenvalues satisfying , . Appendix A discusses the existence and uniqueness of solutions , , as well as the justification of the above linearization and respective stability analysis.

Of note is that our approach is an extension of the non-plastic case. Indeed, without plasticity , the delays remain fixed at the baseline lag , and as a result there is no need for the delays to establish an equilibrium with positive phase-offset differences . Hence, the oscillators become perfectly in-phase with during synchronization. Consequently, Eq. (5) for the global frequency Ω reduces to the one-dimensional fixed-point expression By employing similar steps as above, the phase perturbation terms follow the linear system resulting in the corresponding eigenvalue matrix defined with entries With non-plastic delays, stability analysis at has been accomplished through other means. Lyapunov functionals [21] have shown that a sufficient criterion for synchronization around is for all active connections i, j such that , and Ω satisfies Eq. (15). This criterion becomes necessary given a unique baseline lag [22]. That is, it can be shown that the oscillators synchronize if and only if .

Synchronization of two coupled oscillators with plastic delays

As a first illustrative approach to the general case when the number of coupled oscillators is large, let us first consider a reduced two-oscillator network. The following setup is similar to the system analyzed in [23]. We first let with and remove all self-integration terms by setting . We also set the baseline delays to be equal with . Then the phases , follow the Kuramoto system, with respective plasticity equations By symmetry we may let . We proceed by setting a sufficiently large plasticity gain with . Then it follows that the respective equilibrium delays are and , resulting in a single positive equilibrium delay. Hence, by Eq. (5) the frequency Ω and offset satisfies Solving for above, we obtain where positivity holds only when . This leads to the oscillators synchronizing at a lower frequency . Substituting in Eq. (20) with Eq. (22), if we define the root function then , synchronizes at a common frequency Ω implicitly satisfying . Since it was assumed that , for self-consistency the synchronous frequency Ω must also satisfy , which can be written as since we set . Without plasticity , the oscillators become in-phase with and synchronizes at a frequency Ω that is a root of the function as stated by Eq. (15). Figure 2(A) plots functions on the interval with respect to representative values of the plasticity gain. We observe that higher plasticity gains generally leads to a greater number of potential synchronization frequencies Ω for our system. In Fig. 2(A), one can see that the functions and have single roots within the interval , while has five.

Figure 2

Theoretical stability plots for two-oscillator system. A. Plot of error functions with varying fixed gain (magenta), (yellow), (blue). All roots of are potential synchronization frequencies for the two-oscillator system. The number of roots Ω for increase with larger κ. B. The plasticity gain is set to . Plot of the real part of the non-zero branches , (orange, cyan) of the polynomial root equation over . Ticks on the Ω-axis (blue) indicate the frequencies solving where the system can synchronize. The plotted branches imply that the oscillators will synchronize at , and avoid the unstable frequency with . C, D. Error heatmaps with , respectively, approximate the distribution of eigenvalues solving near , scaled and normalized for visibility. Spots near zero error (white) suggest potential eigenvalue locations. Markers plot the eigenvalues (blue, orange, cyan) for . The heatmap in D indicates an eigenvalue λ near , which implies instability at . All other eigenvalues λ appear to be distributed either at or on the left-side of the imaginary axis. Here, and . For all plots, , , , , and

The stability in the two-dimensional case can be easily determined. Indeed, as derived in Sect. 3, the stability of the oscillators at the synchronization state is determined by the distribution of eigenvalues that satisfy where , , . This results in the root equation , where is the single positive delay and , are polynomials given by where . Note that , which means we have neutral stability. If we make the approximation , the eigenvalues λ correspond to exact cubic roots of . It is possible for the stability of the system under to align with the stability under , particularly for small τ. If we ignore the neutral stability of our zero transcendental equation, Theorem 1 of [24] states that under certain conditions, the stability of the system at and at any does not change. For rigorous purposes however, we are still interested in finding the distribution of eigenvalues solving Eq. (25), particularly if it is accessible through numerical approximation.

Setting the plasticity gain , Fig. 2(B) shows the real parts of the two non-zero branches , of the roots , plotted with respect to the global frequency Ω. If any of the branches are positive, it implies by the above discussion that the oscillators will not synchronize at the state , where is given by Eq. (22). We denote to be the three largest synchronization frequencies corresponding to equilibria solutions of the root equation as presented in Fig. 2(A). The eigenvalue branches imply that frequencies , are stable with at , while is an unstable frequency with at . Figure 2(C), (D) shows heatmap approximations of complex roots λ of Eq. (24) at , , respectively. We can see that at , both non-zero cubic roots are located on the left-side of the imaginary axis, and the error heatmap shows that the eigenvalues are located on the left-hand side of the imaginary axis. At , all cubic roots are real with a single positive root . Consistent with the stability of our system at , the error heatmap implies there exists an eigenvalue satisfying Eq. (25) near the positive root . This highlights that under sufficiently large plasticity gain, adaptive delays introduce multiple stable states in a two-oscillator system, whose stability can be assessed under the zero delay approximation . Our stability analysis reveals that the oscillators can synchronize at multiple possible frequencies, which suggests a greater degree of adaptability in our system.

Synchronization in large-scale oscillator networks with plastic delays

Using inspiration by the two-dimensional case in Sect. 4, we consider here a large-dimensional system using N-limit approximations. Again for simplicity, we set the baseline lags to be constant with and the connection topology to be all-to-all with . We approach the synchronization state of the network with adaptive delays, given by Eqs. (5) and (6), in the following statistical sense. Suppose that the phase-offsets are i.i.d. under some distribution. We can set the offsets to be centered at 0 by defining the re-centered offsets , where ϕ̅ is the mean of ϕ. Then . We take to be i.i.d. under some density function . Setting a random in the global frequency equation (5) for each , and then taking the limit , we obtain the following N-limit approximation for frequency Ω and density : for all fixed , where . In order to apply Eq. (28) to obtain the global frequencies Ω, we parametrize the unknown density by assuming it is Gaussian under some small phase-offset variance . That is, we have and we set . If the variance is small, we can approximate the fixed offset in Eq. (28) as , as well as make the approximation since Δ is small. Hence, our large-scale network synchronizes near a global frequency Ω with Gaussian distributed phase-offsets with variance given that is a root of the function Without plasticity, we recall that we have in-phase synchronization at global frequency Ω which is a solution to . That is, Figure 3(A) plots the curve for various fixed-values of on . We find that there is a unique but different root Ω to at each fixed variance . Hence, we can obtain an implicit curve by parametrizing the level curve with respect to . The curve is plotted in Fig. 3(B), and shows a continuous range of potential synchronization states along to a large N-dimensional system of oscillators. The graph of the level curve also shows that for all , implying that the oscillators are drawn to synchronize at a lower frequency from their natural frequency.

Figure 3

Theoretical stability plots for large N-dim oscillator system. A. Plots of error function with varying fixed over . The function is truncated between interval for visibility. There is a unique root for each fixed . B. Colour map of over states , along with the implicit solution curve (purple) parametrizing level set . Stable regions correspond to (blue) and unstable regions correspond to (red). The network synchronizes near a state overlapping the stable region. C. Plot of stability term along the solution curve over . There is a small interval for which is in the stable region (blue). Other states are in the unstable region (red). D. Complex plot of non-zero eigenvalues of on solution states across varying , scaled by s.log for visibility. The eigenvalues in plot D were computed at respective states in plot C indicated by the same colour. Power terms for polynomial were computed up to degree . The parameters used for all plots are , , , , and

For small offset differences, the equilibrium delays are approximately if , and if . Before proceeding, we set κ to be sufficiently larger than the baseline lag such that . Then most negative offset differences fall below the Heaviside cutoff . From this, the Heaviside term becomes approximately dependent to the sign of with for all i, j. For large gain κ, two categories of equilibrium delays emerge that contribute to the global frequency Ω in Eq. (30): For half the connections with offset difference , the corresponding delays decay to as we have with . Otherwise, with the delay establishes a positive equilibrium at . The positive delays become widely distributed under large standard deviation κδ.

The coupled network can synchronize towards any stable point along the curve . To assess the stability at each state , consider our N-dimensional eigenvalue stability criterion for eigenvector and matrix whose entries are given by Eq. (13). That is, for all , where , , , and . Once again, if our offset variance is small for Eq. (31) we can approximate the terms by assuming at each i. We derive an N-limit version of the eigenvalue Eq. (32) as follows. Likewise to the frequency equation Ω, we can obtain a similar N-limit approximation to Eq. (32) as follows. We define a continuous eigenfunction such that . Then, taking the limit to Eq. (32) with , we obtain the continuous eigenvalue criterion for λ with respect to eigenfunction given by for every . Justification regarding the N-limit step to derive Eq. (33) is provided in Appendix B. Here, the only continuous eigenfunction solution to the above equations is the constant function .1 Hence, Eq. (33) simplifies to the following N-limit eigenvalue equation for λ: We claim that the synchronization state is (neutrally) stable given that for all satisfying Eq. (34), . Note that for all satisfying Eq. (34), where . By Eq. (35), if , , which is a contradiction. It follows that, as , for the distribution of eigenvalues λ satisfying Eq. (34). If , then the frequency Ω solving Eq. (31) is (neutrally) stable if all eigenvalues given by have non-positive real parts . As shown in [22], the non-plastic delay network synchronizes at Ω if and only if .

At this point, the N-limit approximate eigenvalue Eq. (34) has done little to improve the N-dimensional criterion , due to exponential blow-up that persists within the integrand term. We proceed to reduce the eigenvalue Eq. (34) into an exponential polynomial root equation as follows. Rescale by some radius , with rescaled small delay term . Expressing the power series of the exponential λ term up to degree M in Eq. (37), we obtain the approximate exponential polynomial root equation with polynomials and power term coefficients Choosing a large radius , the rescaled delay term and coefficients for larger degrees m become arbitrarily small. From this, we claim that the stability at is predominantly determined by the first few terms of the exponential power expansion. That is, the stability at synchronization state is determined by the finitely many polynomial roots satisfying Eq. (38) with and some degree M. Denoting as the roots of M-degree polynomial , the synchronization state is (neutrally) stable given that We note that these analytic results are reminiscent of what we obtained for two-dimensional systems in Sect. 4 with the cubic exponential polynomial equation as defined by Eqs. (26) and (27).

As we experienced large scale fluctuations of stability term , we processed its value with either the sign function or the sign logorithm defined as In Fig. 3(B), we plot for all and small . Stable regions are indicated where . In computing , eigenvalues satisfying were considered to be non-zero. We notice that a section of the level curve overlaps part of the stable region for . Indeed, Fig. 3(C) plots the stability term along all points on the implicit solution curve . As the plot shows, there is a small interval such that is stable for all such that . Figure 3(D) shows the transition of the eigenvalues of , corresponding to state as δ leaves the interval . We see that the eigenvalues shift towards the right-side of the imaginary axis, which shows the transition from stability to instability along the solution curve .

The results above were obtained by setting the polynomial degree to for , while the approximation failed for powers . In addition, it remains unresolved whether the distribution of eigenvalues satisfying the N-dimensional Eq. (32) with some eigenvector generally converge to the N-limit eigenvalues λ satisfying Eq. (33) with corresponding continuous eigenfunction . That is, whether for each there is some large N such that for all satisfying Eq. (32) each at N dimensions, there exists some eigenvalue satisfying Eq. (33) such that . Integral equation theory has proven some relevant theorems. For instance, it can be proven [25] that there exists a non-zero eigenvalue to Eq. (33) such that for some subsequence of eigenvalues satisfying Eq. (32) at dimension . Nevertheless, our N-limit analysis highlights an approximate region of points where a large-dimensional system of oscillators will synchronize. We demonstrate in the following section that this region aligns with synchronous behavior in numerical simulations.

Comparison to numerical simulations

Here, we validate the theoretical analysis committed in the previous sections through comparisons with numerical simulations. We obtain the global frequency Ω and asymptotic phase-offsets numerically, and systematically compare them with their analytical counterparts. Given numerical solution , to the Kuramoto equation (1) with corresponding derivative solution , we can obtain the asymptotic frequencies for each oscillator by summing over the time interval and taking the limit . That is, and estimate the global frequency Ω as the sample mean of individual frequencies Likewise, we can numerically estimate the asymptotic phase-offsets for each oscillator by taking the difference and defining the limit After modding so that , we can estimate the offset variance by taking the sample variance and the average (asymptotic) phase ϕ̅ by taking the sample mean, If our solutions synchronize towards some synchronous frequency Ω, offsets , then Ω̂, are estimators for Ω, , respectively. Numerically, we evaluate the estimators by taking the average over interval starting at some large time .

We set all baseline delays to be a unique value , so that . Since the plasticity rule Eq. (3) is an ordinary differential equation, it suffices for all delays to have an initial value at . Before , we set the phases to be positioned in accordance to some initial frequency and initial phase-offsets . That is, we define the linear initial function and set on .2 In order to have reasonably behaving solutions , we modify the function so that it satisfies the necessary condition for all i. This adjustment is needed in order to avoid numerical discontinuities for at . For details, refer to Appendix C.

Figure 4 plots the results of a series of numerical simulations in the reduced two-oscillator network as set up in Sect. 4. All trials were run from 0–200 s, with estimated values Ω̂, obtained by averaging the arrays over the last 20 seconds. The same parameter values used in Fig. 2 were applied here when running the numerical trials. We demonstrate the existence of two stable synchronization states with different frequencies. Figure 4(A), (B), (C) shows the asymptotic behaviour of two trials (purple, orange) that started with different initial functions, plotting over the first 50 seconds. Figure 4(A) plots the derivative arrays , of each trial. We observe in Fig. 4(A) that each pair of oscillator frequencies entrain to a common frequency over time. The two trials converge to different common frequencies, estimated to be (purple lines) and (orange lines), respectively. Figure 4(B) plots the offset arrays , of each trial. For each trial, the two-oscillators become phase-locked as converges asymptotically to estimated constant differences (purple lines) and (orange lines), where . Figure 4(C) plots the connection delays , over time. By choosing our plasticity gain , it follows that for both trials, converges to some positive equilibrium delay and decays to 0.

Figure 4

Numerical plots for two-oscillator system. For plots A, B, C, two trials with different initial conditions are graphed. Trial 1, 2 (purple, orange) starts with initial frequency and phase difference , respectively. A. Plots of derivatives , over time. For each trial, (dashed) and (dotted) converge to a common value Ω̂ asymptotically. Each trial of oscillators entrain to a different frequency, implying the existence of multiple synchronization frequencies. B. Plots of sine phases , where . For both trials, the oscillators asymptotically phase-lock with , (dotted), (dashed). The phase-lock difference is also different for the two trials. C. Plots of adaptive delays (dashed), (dotted) over time. For each trial, delay converges to some positive equilibrium , and delay decays to 0. D. Plots showing where two-oscillators with randomized initial conditions (orange) synchronize towards in terms of asymptotic frequency and phase-offset (magenta) across 80 trials. Each of the two trials in A, B, C have their initial condition plotted as a diamond marker of matching colour. Theoretical synchronization states given by Eqs. (20) and (22) are also plotted (blue). Trials converge to two states , which align with the two theoretically stable states shown in Fig. 2. The parameters used for all plots are , , , , ,

The above results imply that there exist at least two stable synchronization states, and that the frequency the system converges towards is dependent on the initial functions , . To provide further confirmation towards this proposition, the trials shown in Fig. 4 were repeated 80 times randomized across sampled initial frequency and initial phase difference . The point (circle markers) defines the initial functions , . The theoretically stable states corresponding to frequencies and are plotted as orange and purple stars, respectively, which are both close to the asymptotic frequency of the matching coloured trial discussed above. The rest of the theoretical synchronous states given by the roots of Eq. (23) are plotted as blue stars. Each trial’s solution arrays synchronized near one of the two stable states, as shown by the connecting coloured lines. The matching colours indicate which of the two stable frequencies the trial’s solution arrays synchronized towards, such that the . Each of the two trials graphed in Fig. 4(A), (B), (C) are also plotted in Fig. 4(D) with a diamond marker of matching colour. Convergence of the points is suggestive of a separatrix curve between the basins of attraction of both stable fixed points. It was also observed that the system synchronized towards the frequency faster than , which suggests that the state has a greater force of attraction. Hence, the experimental results align with the analysis outlined in Sect. 4. We draw the conclusion that in our reduced two-oscillator system, plastic delays are able to generate multiple synchronization states in comparison to non-plastic delays.

Figure 5 provides the numerical results of a similar experiment performed as in Fig. 4 with a large-dimensional network and all-to-all network. All trials were run from 0-100 s, with estimated values Ω̂, obtained by averaging the arrays over the last 10 seconds. The same parameter values are used as in Fig. 3. The initial function for each trial was set up by choosing some frequency and deviation . The initial phases were i.i.d. sampled uniformly from the interval . The ith oscillator was equipped with the initial linear function . Figure 5(A), (B), (C), (D) graphs a single numerical trial with , . Figure 5(A) plots the derivative arrays , which we note converges to some constant frequency estimated as . Figure 5(B) plots the offset arrays , which shows that each converges to some constant phase-offset estimated by . Hence, the oscillators become asymptotically phase-locked under distributed offsets with estimated variance . Figure 5(C) plots a sample of 50 adaptive delays , which become part of the positive equilibrium distribution or decay to 0. Figure 5(D) plots the density of centered phases . As we assumed that the centralized phase-offsets follow a Gaussian distribution, we perform a normality test on the numerical asymptotic offsets . The Shapiro–Wilk test for non-normality returned a p-value of 0.005 for , which suggests another distribution would be more accurate. Visually, a Gaussian curve (black line) is fit over the density in Fig. 5(D). Nevertheless, the relevance of the Gaussian approximation, which greatly simplifies the analysis, becomes apparent as the numerical and analytical results nearly coincide. Other approximations could be used to facilitate the analysis further, and are left for future work.

Figure 5

Numerical plots of N-oscillator system. A. Plots of derivatives over time. We see that all converge to a common frequency, estimated to be . B. Plots of sine phases , where . All oscillators appear to asymptotically phase-lock to one another. C. Plots of a sample of 50 adaptive delays over time. Some delays converge to some positive equilibrium , while others decay to 0. D. Density of centralized phase-offsets , which was assumed to be Gaussian. The Gaussian curve (black line) is fit over the density. E. Plot showing oscillators with randomized initial frequency and phase deviation (yellow) synchronizing towards respective estimated asymptotic frequency and phase deviation (magenta star) across 10 trials. The yellow diamond refers to the trial plotted in A, B, C. The numerical values are plotted directly over Fig. 3(B). As shown, the network synchronizes near the theoretical stable region. The parameters used were , , , , , ,

The numerical simulation, as presented above, was repeated 10 times with randomized initial conditions . For each trial, , was sampled uniformly from intervals and , respectively. Figure 5(E) plots the following convergence results. Each trial with initial condition (yellow markers) synchronized near the respective estimated point (magenta markers). We can see that every trial synchronized at approximately the same state . To determine whether the numerical results align with the analysis discussed in Sect. 5, the numerical values were plotted on top of Fig. 3(B). We observe that for each trial, the network synchronizes near the portion of the level curve within the stable region where . Hence, the numerical experiment for generated results that validates the theoretical N-limit stability analysis in Sect. 5.

Neuroscience application: resilience to injury with sparse and uniform connectivity

One of the most salient examples of white matter plasticity comes from neuroimaging in the presence of a lesion. In these cases, white matter remodeling takes place in order to restore and maintain function, a process that notably impacts neural synchronization [26]. Let us now investigate whether the plasticity mechanism can be used to stabilize phase-locked states in the presence of network damage. That is, we would like to know whether changes in time delays can be used to compensate for a reduction in effective connectivity and make the global synchronous state more resilient. To investigate this problem, we model injury as a loss in connections . Defining as the insult index, we introduce here the sparse synaptic connectivity weights given by where is a uniformly distributed i.i.d. sampling on . γ represents the connectivity damage through an increase in the sparseness of network connections. Note that if , then we obtain the all-to-all connection topology corresponding to no injury in the system. If , we have the trivial system which means no signals between the oscillators occur. For any γ we have the probability to retain the connection . Without plasticity, the mean phase dynamics for is given by see Eq. (31). By observation, one can see that we are in the presence of the same network dynamics, but with an effective coupling coefficient given by . Thus, damage simply reduces the net coupling. To demonstrate this point, we start with a strong coupling parameter and decrease it until stability is either lost or preserved by plasticity. According to Eq. (51), as . Hence, if the baseline lag is chosen such that and , the network without plastic delays is susceptible to injury destabilizing the synchronous state.

We ran numerical simulations by applying similar parameter values as Sect. 6 while introducing injury to the connections at time , set at . The initial condition of the network was fixed at for all trials. Figure 6(A) shows the destruction of existing connections following injury, comparing connection grids for before and after . Figure 6(B) shows the distribution of the delays with existing connections at timestamps , (pre-injury), (post-injury). With fewer connections available, the ability for the surviving delays to adjust themselves are crucial in re-stabilizing the system’s synchrony. Figure 6(C) (no gain) and Fig. 6(D) (with gain) show that both networks entrain to a global frequency successfully before and after inflicted injury towards the network. The entrainment frequencies pre-injury and post-injury were estimated by taking the average of at times 148–160 s and 304–320 s, respectively, for each . Figure 6(E) (no gain) and Fig. 6(F) (with gain) plots the sine phase-offsets over time, given by Following injury, from Fig. 6(E) the network without adaptive delays collectively falls out of phase. In contrast, Fig. 6(F) shows the network with adaptive delays demonstrating resilience against the injury as most oscillators are able to collectively phase-lock within close proximity to each other.

Figure 6

Comparing resilience against injury between plastic and non-plastic delays. Injury towards the connection topology is introduced at (red line). A. Plots of the connection matrix before injury () and after injury (), with indicated in white, black, respectively. B. The log histogram plots of delays at initial time (purple), midtime before injury (green), and at the end time following injury (red). The delays become distributed away from to either some largely varying equilibrium delays or decay to 0. Following injury, there are fewer delays available to stabilize the synchronous network. C, D. Plots of derivatives over time, without and with plasticity, respectively. Both networks entrain in frequency pre-injury. Following injury, both networks entrain to a new frequency closer to . E, F. Plots of over time, without and with plasticity, respectively. Following injury, the oscillators with plastic delays are able to coherently phase-lock within close proximity to each other, while the network without plastic delays remain in a non-convergent state. The parameters used were , , , , , , and

Figure 7 examines the effect of gradually increasing the severity of injury towards the system’s global frequency and its phase-offset variance . For each trial at injury γ, the same initial condition and parameters were used as in Fig. 6. Figure 7(A) shows that connection loss generally leads to the system’s synchronization frequency Ω becoming closer to the natural frequency . Figure 7(B) plots the estimated phase standard deviation δ̂ with respect to increasing injury. The network without plastic delays exhibits a significant loss in coherent synchrony with increasing as more connections are lost. In contrast, the network equipped with adaptive delays persistently displays phase coherence with until higher injury levels γ.

Figure 7

Comparing of post-injury network asymptotic behaviour with increasing injury between plastic and non-plastic delays. Each numerical trial was run over 320 s, and all asymptotic values were evaluated by averaging over the final 16 s. Injury was introduced at . Both plots show trials with adaptive delays (red) and fixed delays (purple). A. Plot of post-injury asymptotic frequency Ω̂ for trials with injury . As γ increases, . B. Plot of post-injury asymptotic standard deviation δ̂ for trials with injury . As γ increases, δ̂ has significant increase for trials without plasticity, while δ̂ remains relatively small for trials with plasticity until . The parameters used were , , , , , , and

Discussion

Our goal was to provide a mathematical framework that captures the synchronizing properties of networks with adaptive delays. We sought to implement the activity-dependent property of myelinated connection delays by modifying the Kuramoto model as proposed in [9]. The focus was to determine whether such adaptive delays significantly improve the oscillatory system’s ability to become in-phase and to entrain to a global frequency. Given that this is the case, the results of the model’s study reinforces the proposition that myelin plasticity is essential in maintaining the synchrony in the developing or injured brain.

White matter plays a critical role in maintaining brain function through the coordination of neural dynamics across multiple temporal and spatial scales. Recent evidence has shown that through the action of glia, white matter properties evolve continuously in time. Specifically, conduction velocity within and across brain areas is adjusted to promote efficient neural signaling. While the mechanism remains poorly understood, the consequences of such plastic processes on brain dynamics and synchronization can be readily examined and characterized using simplified mathematical models.

To accomplish this, we here examined the influence of adaptive conduction delays on the synchronization of neural oscillators. We developed a repertoire of mathematical tools to better examine the stability of phase-locked solutions. In theory, we derived implicit equations for the global frequency Ω and eigenvalues that provide a rigorous criterion for the stability around the synchronous state in two-dimensional and large-dimensional settings. Based on our model, flexibility in the white matter structure introduces an additional corrective dynamic next to the phase interactions that can further drive the network’s phase alignment. Higher stability with adaptive delays was demonstrated as the Kuramoto model had higher resilience against injury perturbations. However, adaptive delays improve the system’s synchronous features only when the delays adjust with a sufficiently high degree of plasticity, as represented by the plasticity gain κ.

There are many limitations in the prototypical model we used and its corresponding results. Myelination is bound by many physiological constraints, some of which remain uncertain [4]. It is established that white matter restructures itself in response to ongoing neural activity [8]. We primarily incorporated this fact in our plasticity rule in a manner that promotes local synchrony. Indeed, each connection delay changes at a rate proportional to the sine of the oscillator’s phase difference. This rule remains a tentative construction, as more research is needed to develop more biological relevant models in activity-dependent myelination. In addition, the use of phase oscillators to model local neural dynamics remain limited and is relevant mostly in the context of large-scale neural systems. In our analysis, we relied heavily on i.i.d. parametrical frameworks in order to establish our N-limit approach, which may not be feasible as network elements are correlative in nature.

Moving forward, we hope to build upon our analysis alongside newly found experimental results pertaining to myelin. Despite its shortcomings, the mathematical approaches used and its results can potentially be applied to more complex, biological relevant models. The conduction delays can be alternatively modelled with respect to a system of adaptive conduction velocities . In the realm of temporal equations, other parametric avenues have yet to be explored. For instance, the delays can exhibit slow convergence by setting the rate constant . The aforementioned concepts are some proposed examples that may further lead to uncharted dynamics in the scope of neurocomputational models.