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 activitydependent 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][14][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 largedimensional 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 θ i (t) evolve according to the following system of differential equations: where ω i ∈ R is the natural frequency of oscillator i, and g > 0 is a constant global coupling gain. The coefficients a ij represent synaptic weights: a ij = 1 if there is a white matter connection between oscillator j to i; otherwise, a ij = 0. The axonal conduction delays τ ij correspond to the conduction time taken by an action potential along a myelinated fiber. The propagation speed c ij = c ij (t, ) 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 c ij (t, ) also changes in time. As such, each conduction delay may be written as That is, a line integral along the axonal pathway P ij 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 τ ij (t) evolves under the following sublinear plasticity equation: where 1 ≤ i, j ≤ N . Here, α τ is the homeostatic rate constant that sets the time scale of the evolving delays. The plasticity coefficient κ > 0 sets the gain of the conduction delay changes. The initial condition for delays τ ij (t) at t = 0 is provided by baseline lags τ 0 ij ≥ 0. That is, τ ij (0) = τ 0 ij for all i, j. The Heaviside function H(τ ) is defined as a smooth function such that H(τ ) = 0 if τ < 0 and H(τ ) = 1 if τ > ε for some small ε > 0. The Heaviside function H(τ ) is present to preserve non-negativity of the delays τ ij (t) at all times t ≥ 0. 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 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) initial baseline lag τ 0 ij representing minimal myelination, and a forcing term that depends on the phase difference ij := θ jθ i 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 |θ jθ i | < π , one can see that whenever θ j > θ i , the time delay increases due to an effective reduction in conduction velocity caused by myelin retraction. The opposite occurs if θ j < θ i 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 τ 0 ij .

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 N + N 2 -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][19][20]. Mathematically, we analyze the network's ability to asymptotically achieve the following synchronized state: for all i ≤ N as t → ∞, where Ω is the global fixed frequency of the oscillators and φ i is the asymptotic phase-offset of oscillator i as time t → ∞. As we will see, adaptive delays following the plasticity rule equation (3) require non-zero phase differences ij (t) := θ j (t)θ i (t) 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 θ i (t) globally oscillate at some stable frequency Ω and become phase-locked under individual offsets φ i . 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 ω i = ω 0 for all i in Eq. (1). To simplify the convergent behaviour of delays τ ij (t) following Eq. (3), we set α τ = 1.
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 Ω ∈ R and individual offsets φ = (φ 1 , . . . , φ N ) ∈ R N : To assess the network's stability at the synchronized state (Ω, φ) satisfying Eqs. (5) and (6), we introduce the perturbation terms i (t), η ij (t) around the equilibrium states Eq. (4) and Eq. (6) for θ i (t) and τ ij (t), respectively. That is, we write and examine the stability at the equilibrium point i (t) = η ij (t) = 0. The delay perturbations η ij (t) abide by the linearized form of Eq. (3) around all positive delays τ E ij = τ 0 ij + κ sin( ij ) > 0. For the equilibrium delays τ E ij = 0, we proceed by assuming the corresponding perturbation terms will act in an asymptotically stable manner. Specifically, for all phase-offset differences such that τ E ij = 0, by the nature of the Heaviside cutoff function H(τ ) in Eq. (3) there exists some t asy ∈ R such that for all t > t asy , τ ij (t) = 0 and consequently η ij (t) = 0. That is, if τ E ij = 0, the perturbation term is asymptotically η ij (t) = 0. 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 η ij (t) as τ ij (t) → τ E ij are given by where C 0 ij := cos( ij ). We can use the linearization approach in the context of equations with state-dependent delays [18][19][20] as follows. Inserting Eq. (7) and Eq. (8) into the Kuramoto equation (1), where we have made the approximation j (tτ ij (t)) ≈ j (tτ E ij ) when τ ij (t) is near τ E ij . With (Ω, φ) satisfying Eq. (5), and by taking a first-order expansion around the term -Ωτ E ij , we obtain the linear system for the phase perturbation term given by where we denote C ij = cos(-Ωτ E ij + ij ) and H ij = H(τ E ij ). Hence, the network synchronizes at (Ω, φ) given that all perturbation terms i (t), η ij (t) 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 i (t) = v i e λt , η ij (t) = w ij e λt with respect to eigenvalue λ ∈ C. By Eq. (9), the coefficients w ij satisfy for all i, j, λ = -1. Applying Eq. (12) to the coefficients w ij , That is, λ is an eigenvalue if there exists an eigenvector v = (v 1 , . . . , v N ) ∈ C N satisfying Eq. (13). In matrix form, Eq. (13) can be expressed as In summary, we acquire the following stability criterion: The system is stable around the synchronized state (Ω, φ) if for all eigenvalues λ ∈ C satisfying det M λ = 0, Re(λ) < 0. Appendix A discusses the existence and uniqueness of solutions θ i (t), τ ij (t), 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 κ = 0, the delays remain fixed at the baseline lag τ 0 ij , and as a result there is no need for the delays to establish an equilibrium with positive phase-offset differences ij . Hence, the oscillators become perfectly in-phase with φ i = 0 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 i (t) follow the linear system resulting in the corresponding eigenvalue matrix M λ defined with entries a ij cos -Ωτ 0 ij , With non-plastic delays, stability analysis at θ i (t) = Ωt has been accomplished through other means. Lyapunov functionals [21] have shown that a sufficient criterion for synchronization around θ i (t) = Ωt is cos(Ωτ 0 ij ) > 0 for all active connections i, j such that a ij = 1, and Ω satisfies Eq. (15). This criterion becomes necessary given a unique baseline lag τ 0 ij = τ 0 [22]. That is, it can be shown that the oscillators synchronize if and only if cos(Ωτ 0 ) > 0.

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 N = 2 with a 12 = a 21 = 1 and remove all selfintegration terms by setting a 11 = a 22 = 0. We also set the baseline delays to be equal with τ 0 12 = τ 0 21 = τ 0 . Then the phases θ 1 (t), θ 2 (t) follow the Kuramoto system, with respective plasticity equations By symmetry we may let 12 = φ 2φ 1 > 0. We proceed by setting a sufficiently large plasticity gain with κ τ 0 . Then it follows that the respective equilibrium delays are τ E 12 := τ 0 + κ sin( 12 ) and τ E 21 = 0, resulting in a single positive equilibrium delay. Hence, by Eq. (5) the frequency Ω and offset 12 satisfies Solving for 12 above, we obtain 12 = arcsin where positivity holds only when Ω < ω 0 . This leads to the oscillators synchronizing at a lower frequency Ω ∈ [ω 0g, ω 0 ). Substituting 12 in Eq. (20) with Eq. (22), if we define the root function then θ 1 (t), θ 2 (t) synchronizes at a common frequency Ω implicitly satisfying R κ (Ω) = 0.

Synchronization in large-scale oscillator networks with plastic delays
Using inspiration by the two-dimensional case in Sect. 4, we consider here a largedimensional system using N -limit approximations. Again for simplicity, we set the baseline lags to be constant with τ 0 ij ≡ τ 0 and the connection topology to be all-to-all with a ij ≡ 1. 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 φ i are i.i.d. under some distribution. We can set the offsets to be centered at 0 by defining the re-centered offsets i : We take i to be i.i.d. under some density function ρ( ). Setting a random i in the global frequency equation (5) for each i ≤ N , and then taking the limit N → ∞, we obtain the following N -limit approximation for frequency Ω and density ρ( ): for all fixed ∈ R, where τ E ( ) = H(τ 0 + κ sin ) · (τ 0 + κ sin ). 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 δ 2 . That is, we have and we set ρ( ) = ρ δ ( ). If the variance δ 2 is small, we can approximate the fixed offset in Eq. (28) as ≈ 0, 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 δ 2 given that (Ω, δ 2 ) is a root of the function Without plasticity, we recall that we have in-phase synchronization at global frequency Ω which is a solution to R 0 (Ω, 0) = 0. That is, We find that there is a unique but different root Ω to R(Ω, δ 2 ) = 0 at each fixed variance δ 2 > 0. Hence, we can obtain an implicit curve Ω = Ω(δ) by parametrizing the level curve R(Ω, δ 2 ) = 0 with respect to δ > 0. The curve is plotted in Fig. 3(B), and shows a continuous range of potential synchronization states (Ω(δ), δ 2 ) along δ > 0 to a large N -dimensional system of oscillators. The graph of the level curve R(Ω, δ 2 ) = 0 also shows that Ω(δ) < ω 0 for all δ > 0, implying that the oscillators are drawn to synchronize at a lower frequency from their natural frequency. For small offset differences, the equilibrium delays are approximately τ E ij ≈ τ 0 + κ ij if ij > 0, and τ E ij = 0 if ij < 0. Before proceeding, we set κ to be sufficiently larger than the baseline lag τ 0 such that τ 0 κ ≈ 0. Then most negative offset differences ij < 0 fall below the Heaviside cutoff τ 0 + κ ij < 0. From this, the Heaviside term becomes approximately dependent to the sign of ij with H ij ≈ H( ij ) 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 ij < 0, the corresponding delays τ ij (t) decay to τ E ij = 0 as we have τ 0 < -κ ij < 0 with τ 0 κ. Otherwise, with ij > 0 the delay τ ij (t) establishes a positive equilibrium at τ E ij = τ 0 + κ ij . The positive delays τ E ij become widely distributed under large standard deviation κδ. The coupled network can synchronize towards any stable point (Ω, δ 2 ) along the curve R(Ω, δ 2 ) = 0. To assess the stability at each state (Ω, δ 2 ), consider our N -dimensional eigenvalue stability criterion M λ v = 0 for eigenvector v ∈ C N and matrix M λ whose entries are given by Eq. (13). That is, for all i ≤ N , whereκ = Ωκ, C 0 ( ) = cos( ), C Ω ( ) = cos(-Ωτ E ( ) + ), and τ E ( ) = H( )(τ 0 + κ ). Once again, if our offset variance δ 2 is small for Eq. (31) we can approximate the terms ij = ji ≈ j by assuming i ≈ 0 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 v : [0, 1] → C such that v i = v(i/N). Then, taking the limit N → ∞ to Eq. (32) with ij ≈ j ∼ N(0, δ 2 ), we obtain the continuous eigenvalue criterion for λ with respect to eigenfunction v(x) given by for every x ∈ [0, 1]. 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 v(x) = 1. a Hence, Eq. (33) simplifies to the following N -limit eigenvalue equation for λ: We claim that the synchronization state (Ω, δ 2 ) is (neutrally) stable given that for all λ ∈ C satisfying Eq. (34), Re(λ) ≤ 0. Note that for all λ = u + iv satisfying Eq. (34), where which is a contradiction. It follows that, as |λ| → ∞, Re(λ) → -∞ for the distribution of eigenvalues λ satisfying Eq. (34). If κ = 0, then the frequency Ω solving Eq. (31) is (neutrally) stable if all eigenvalues λ ∈ C given by λ = g cos Ωτ 0 e -λτ 0 -1 have non-positive real parts Re(λ) ≤ 0. As shown in [22], the non-plastic delay network synchronizes at Ω if and only if cos(Ωτ 0 ) > 0. At this point, the N -limit approximate eigenvalue Eq. (34) has done little to improve the N -dimensional criterion det M λ = 0, 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 λ → Rλ by some radius R > κ, with rescaled small delay term τ R = τ 0 /R. Expressing the power series of the exponential λ term up to degree M in Eq. (37), we obtain the approximate exponential polynomial root equation and power term coefficients Choosing a large radius R > κ, the rescaled delay term τ R and coefficients I m for larger degrees m become arbitrarily small. From this, we claim that the stability at (Ω, δ 2 ) is predominantly determined by the first few terms of the exponential power expansion. That is, the stability at synchronization state (Ω, δ 2 ) is determined by the finitely many polynomial roots λ ∈ C satisfying Eq. (38) with τ R ≈ 0 and some degree M. Denoting Λ(Ω, δ 2 ) = {λ 1 , . . . , λ M } as the roots of M-degree polynomial P(λ | Ω, δ 2 ) + Q(λ | Ω, δ 2 ), the synchronization state (Ω, δ 2 ) is (neutrally) stable given that We note that these analytic results are reminiscent of what we obtained for twodimensional systems in Sect. 4 with the cubic exponential polynomial equation P Ω (λ) + Q Ω (λ)e -λτ = 0 as defined by Eqs. (26) and (27).
As we experienced large scale fluctuations of stability term E(Ω, δ 2 ), we processed its value with either the sign function sgn(x) or the sign logorithm s. log(x) defined as s. log(x) = sgn(x) log 1 + |x| .

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 φ i numerically, and systematically compare them with their analytical counterparts. Given numerical solution θ i (t), i ≤ N to the Kuramoto equation (1) with corresponding derivative solution θ i (t), i ≤ N , we can obtain the asymptotic frequencies for each oscillator by summing over the time interval [t, t + h] and taking the limit t → ∞. That is, and estimate the global frequency Ω as the sample mean of individual frequencies Likewise, we can numerically estimate the asymptotic phase-offsets φ i for each oscillator by taking the difference φ i (t) := θ i (t) -Ωt and defining the limit After modding φ i so that φ i ∈ [-π, π), we can estimate the offset variance δ 2 by taking the sample variance 2 (47) and the average (asymptotic) phase φ by taking the sample mean, If our solutions θ i (t) synchronize towards some synchronous frequency Ω, offsets φ i , then Ω, φ i are estimators for Ω, φ i , respectively. Numerically, we evaluate the estimators by taking the average over interval [t, t + h] starting at some large time t > 0. We set all baseline delays to be a unique value τ 0 > 0, so that τ ij (0) = τ 0 ij = τ 0 . Since the plasticity rule Eq. (3) is an ordinary differential equation, it suffices for all delays to have an initial value at t = 0. Before t = 0, we set the phases θ i (t) to be positioned in accordance to some initial frequency Ω 0 and initial phase-offsets φ 0 i . That is, we define the linear initial function ϕ i (t) = Ω 0 t + φ 0 i and set θ i (t) = ϕ i (t) on t ≤ 0. b In order to have reasonably behaving solutions θ (t), we modify the function ϕ i (t) so that it satisfies the necessary condition for all i. This adjustment is needed in order to avoid numerical discontinuities for θ i (t) at t = 0. For details, refer to Appendix C. Figure 4 plots the results of a series of numerical simulations in the reduced twooscillator network as set up in Sect. 4. All trials were run from 0-200 s, with estimated values Ω, φ i 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 θ 1 (t), θ 2 (t) 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 Ω = 0.916 (purple lines) and Ω = 0.625 (orange lines), respectively. Figure 4(B) plots the offset arrays sin φ 1 (t), sin φ 2 (t) of each trial. For  which align with the two theoretically stable states shown in Fig. 2. The parameters used for all plots are α τ = 0.5, ε = 0.01, g = 1.5/2, κ = 30, ω 0 = 1.0, τ 0 = 0.1 s each trial, the two-oscillators become phase-locked as φ 2 (t)φ 1 (t) converges asymptotically to estimated constant differences 12 = 0.111 (purple lines) and 12 = 0.523 (orange lines), where 12 = φ 2φ 1 . Figure 4(C) plots the connection delays τ 12 (t), τ 21 (t) over time. By choosing our plasticity gain κ = 30 0.1 = τ 0 , it follows that for both trials, τ 12 (t) converges to some positive equilibrium delay τ E and τ 21 (t) decays to 0.
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 ϕ 1 (t), ϕ 2 (t). To provide further confirmation towards this proposition, the trials shown in Fig. 4 were repeated 80 times randomized across sampled initial frequency Ω 0 ∈ [ω 0g, ω 0 + g] and initial phase difference 0 ∈ (0, 1). The point (Ω 0 , 0 ) (circle markers) defines the initial functions ϕ 1 (t) = Ω 0 t, ϕ 2 (t) = Ω 0 t + 0 . The theoretically stable states corresponding to frequencies Ω = Ω 1 and Ω = Ω 3 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 Ω i the trial's solution arrays synchronized towards, such that the | Ω -Ω i | < 5× 10 -3 . 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 Ω = Ω 3 faster than Ω = Ω 1 , which suggests that the state Ω = Ω 3 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 N = 50 and all-to-all network. All trials were run from 0-100 s, with estimated values Ω, φ i obtained by averaging the arrays over the last 10 seconds. The same parameter values are used as in Fig. 3. The initial function ϕ(t) for each trial was set up by choosing some frequency Ω 0 ∈ [ω 0 -g 4 , ω 0 + g 4 ] and deviation δ 0 ∈ (0, 0.5). The initial phases φ 0 i were i.i.d. sampled uniformly from the interval [-√ 3δ 0 , √ 3δ 0 ]. The ith oscillator was equipped with the initial linear function Figure 5(A), (B), (C), (D) graphs a single numerical trial with Ω 0 = 0.913, δ 0 = 0.295. Figure 5(A) plots the derivative arrays θ i (t), which we note converges to some constant frequency estimated as Ω = 0.839. Figure 5(B) plots the offset arrays sin φ i (t), which shows that each φ i (t) converges to some constant phase-offset estimated by φ i . Hence, the oscillators become asymptotically phase-locked under distributed offsets with estimated variance δ 2 = 0.050 2 . Figure 5(C) plots a sample of 50 adaptive delays τ ij (t), which become part of the positive equilibrium distribution τ E ij > 0 or decay to 0. Figure 5(D) plots the density of centered phases i = φ iφ. As we assumed that the centralized phase-offsets follow a Gaussian distribution, we perform a normality test on the numerical asymptotic offsets i . The Shapiro-Wilk test for non-normality returned a pvalue of 0.005 for i , which suggests another distribution would be more accurate. Visually, a Gaussian curve N(0, δ 2 ) (black line) is fit over the density in Fig. 5(D). Nevertheless, the relevance of the Gaussian approximation, which greatly simplifies the analysis, Figure 5 Numerical plots of N-oscillator system. A. Plots of derivatives θ i (t) over time. We see that all θ i (t) converge to a common frequency, estimated to be Ω = 0.839. B. Plots of sine phases sin( φ i (t)), where φ i (t) = θ i (t) -Ωt. All oscillators appear to asymptotically phase-lock to one another. C. Plots of a sample of 50 adaptive delays τ ij (t) over time. Some delays τ ij (t) converge to some positive equilibrium τ E ij , while others decay to 0. D. Density of centralized phase-offsets i = φ i -φ, which was assumed to be Gaussian. The Gaussian curve N(0, δ 2 ) (black line) is fit over the density. E. Plot showing oscillators with randomized initial frequency and phase deviation (Ω 0 , δ 0 ) (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 N = 50, ε = 0.01, α τ = 0.1, g = 1.5, κ = 80, ω 0 = 1.0, τ 0 = 0.1 s 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.
The numerical simulation, as presented above, was repeated 10 times with randomized initial conditions (Ω 0 , δ 2 0 ). For each trial, Ω 0 , δ 0 was sampled uniformly from intervals [ω 0 -g 4 , ω 0 + g 4 ] and (0, 1), respectively. Figure 5(E) plots the following convergence results. Each trial with initial condition (Ω 0 , δ 2 0 ) (yellow markers) synchronized near the respective estimated point ( Ω, δ 2 ) (magenta markers). We can see that every trial synchronized at approximately the same state (Ω, δ 2 ). 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 R(Ω, δ 2 ) = 0 within the stable region where E(Ω, δ 2 ) < 0. Hence, the numerical experiment for N = 50 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 phaselocked 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 a ij . Defining γ ∈ [0, 1] as the insult index, we introduce here the sparse synaptic connectivity weights given by where p ij is a uniformly distributed i.i.d. sampling on [0, 1]. γ represents the connectivity damage through an increase in the sparseness of network connections. Note that if γ = 0, then we obtain the all-to-all connection topology a ij ≡ 1 corresponding to no injury in the system. If γ = 1, we have the trivial system a ij = 0 which means no signals between the oscillators occur. For any γ we have the probability to retain the connection P(a ij = 1) = 1γ . Without plasticity, the mean phase dynamics for N → ∞ 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 g eff := (1γ )g. 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), Ω → ω 0 as γ → 1. Hence, if the baseline lag τ 0 is chosen such that cos(Ωτ 0 ) ≥ 0 and cos(ω 0 τ 0 ) < 0, 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 γ = 0.8 to the connections at time t = t inj , set at t inj = 80 s. The initial condition of the network was fixed at (Ω 0 , δ 0 ) = (ω 0 , 0.25) for all trials. Figure 6(A) shows the destruction of existing connections following injury, comparing connection grids for a ij before and after t inj . Figure 6(B) shows the distribution of the delays τ ij (t) with existing connections a ij = 1 at timestamps t = 0 s, t = 79 s (pre-injury), t = 160 s (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 Ω pre ≈ Ω pre and post-injury Ω post ≈ Ω post were estimated by taking the average of θ i (t) at times 148-160 s and 304-320 s, respectively, for each i ≤ N . Figure 6(E) (no gain) and Fig. 6(F) (with gain) plots the sine phase-offsets sin φ i (t) 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 a ij is introduced at t inj = 160 s (red line). A. Plots of the connection matrix A = (a ij ) before injury (γ = 0) and after injury (γ = 0.8), with a ij = 1, 0 indicated in white, black, respectively. B. The log histogram plots of delays at initial time t = 0 s (purple), midtime before injury t = 160 s (green), and at the end time following injury t = 320 s (red). The delays become distributed away from τ ij = τ 0 to either some largely varying equilibrium delays τ E ij > 0 or decay to 0. Following injury, there are fewer delays available to stabilize the synchronous network. C, D. Plots of derivatives θ i (t) 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 ω 0 . E, F. Plots of sin φ i (t) 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 N = 50, ε = 0.01, α τ = 1.0, g = 1.5, κ = 80, ω 0 = 1.0, and τ 0 = 2.0 s Figure 7 examines the effect of gradually increasing the severity of injury towards the system's global frequency Ω ≈ Ω and its phase-offset variance δ 2 ≈ δ 2 . 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 ω 0 . 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 δ > 0 as more connections are lost. In contrast, the network equipped with adaptive delays persistently displays phase coherence with δ 1 until higher injury levels γ .

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 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 t = 160 s. Both plots show trials with adaptive delays (red) and fixed delays (purple). A. Plot of post-injury asymptotic frequency Ω for trials with injury γ > 0. As γ increases, Ω → ω 0 . B. Plot of post-injury asymptotic standard deviation δ for trials with injury γ > 0. As γ increases, δ has significant increase for trials without plasticity, while δ remains relatively small for trials with plasticity until γ = 0.9. The parameters used were N = 50, ε = 0.01, α τ = 1.0, g = 1.5, κ = 80, ω 0 = 1.0, and τ 0 = 2.0 s 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 λ ∈ C 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 activitydependent 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.

Appendix B: Derivation of N-limit equations
Here, we justify the limit steps taken in Sect. 5 to derive the limit-N equations for global frequency Ω and stability eigenvalues λ around the synchronous state. First, denote Eg(X) as the expectation of random variable g(X). That is, given X has distribution ρ(x), The main proposition we applied is as follows.
Theorem 1 Let (v k ) k≥1 be a bounded sequence of real numbers, and X be a random variable such that EX = 0, E|X| < ∞. Suppose (X k ) k≥1 is an i.i.d. sequence of random variables such that each X i has the same distribution as X. Then What directly follows from Theorem 1 is the following limit result, which we utilize to derive the limit-N equations for the global frequency Ω and the stability eigenvalues λ. Corollary 1 Let (v k ) k≥1 be a bounded sequence in C, and X be a random variable such that E|X| < ∞. Suppose (X k ) k≥1 is an i.i.d. sequence of random variables such that each X i has the same distribution as X. Then where v is the well-defined asymptotic sample mean of the coefficients: v = lim We applied Corollary 1 to random sequences of the form X k = f ( k , τ 0 ) where k is Gaussian distributed, and τ 0 is the constant baseline lag. The eigenfunction coefficients v k are from ansatz k (t) = v k e λt as N → ∞.

Appendix C: Numerical setup
All numerical simulations in Sects. 6, 7 were done using MATLAB's ddesd function. We have a system of delay differential equations for the phases θ i (t), i ≤ N given by Eq. (1) and an N 2 system of ordinary differential equations for the state-dependent delays τ ij (t), 1 ≤ i, j ≤ N given by Eq. (3). The initial data before starting time t = 0 are as follows. We have a history function ϕ : [-r, 0] → R N for the phases, such that θ i (t) = ϕ i (t) for all t ≤ 0. We considered only the simple case where delays τ ij (t) have a unique baseline lag τ ij (0) = τ 0 ∈ R. As shown in Appendix A, it is sufficient to consider C 1 initial phase functions ϕ(t) to ensure unique solutions. However, we are not guaranteed to have reasonably behaving solutions unless we have continuity of θ (t) at t = 0 [19]. This is also a necessary condition for providing accurate numerical simulations, as interpolation steps in ddesd rely on θ (t) being continuous everywhere.
For our numerical experiments in Sects. 6 and 7, we consider the linear initial functions ϕ i (t) = Ω 0 t + φ 0 i , i ≤ N with respect to initial frequency Ω 0 ∈ R and phase-offsets (φ 0 1 , . . . , φ 0 N ) ∈ R N . As discussed above, we require that our linear initial function must satisfy the necessary condition Eq. (49). Given any initial C 1 function ϕ(t) for θ (t), we define the modified C 1 functionφ(t) of ϕ(t) that satisfies the necessary condition as follows. We can obtain the modified slope v ∈ R N by imposing the condition on ϕ(t). That is, sin ϕ j -τ 0ϕ i (0) .
Then for each i we define the cubic polynomial p i : R → R that interpolates ϕ(t) between t = 0, -τ 0 using the modified slope. That is, p i (t) = ϕ i (t) at t = 0, -τ 0 and p i (0) = v i , p i (-τ 0 ) = ϕ i (-τ 0 ). Then the modified initial functionφ(t) is given bỹ for i ≤ N . All numerical trials were conducted using the corresponding modified functions ϕ i (t) of ϕ i (t) = Ωt + φ 0 i , i ≤ N .

a(x, y)v(y) dy
with respect to eigenvalue σ . b Trials were done using other types of initial functions ϕ i (t), and yielded similar convergence results. However, the synchronization occurred at a slower rate.

Publisher's Note
Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.