Analytical prediction of specific spatiotemporal patterns in nonlinear oscillator networks with distance-dependent time delays

We introduce an analytical approach that allows predictions and mechanistic insights into the dynamics of nonlinear oscillator networks with heterogeneous time delays. We demonstrate that time delays shape the spectrum of a matrix associated to the system, leading to the emergence of waves with a preferred direction. We then create analytical predictions for the specific spatiotemporal patterns observed in individual simulations of time-delayed Kuramoto networks. This approach generalizes to systems with heterogeneous time delays at finite scales, which permits the study of spatiotemporal dynamics in a broad range of applications.


I. INTRODUCTION
What is the effect of heterogeneous time delays in networked systems?This question is difficult to treat analytically in the context of multiple distributed time delays.In recent work [1], we studied intracranial electrophysiological recordings from human clinical patients during sleep.We found that the 11-15 Hz sleep "spindle" oscillation, a brain rhythm important for learning and memory [2], was not perfectly synchronized with zero phase difference across the cortex; rather, sleep spindles are organized into rotating waves that travel in a preferred direction (see Movie 1 in [1]).Importantly, the propagation speed of the observed waves is consistent with the axonal conduction speed of the long-range fiber network in the cortex (3-5 ms −1 [3]).This set of observations raises an important question: how do these fibers, with no major anisotropy, create a specific spatiotemporal structure with a preferred chirality?
In this work, we analyze a time-delay Kuramoto model to address this question.Utilizing a recently reported analytical approach to the Kuramoto dynamics [4], we introduce a complex-valued delay operator.This operator shapes the dynamics of the Kuramoto system into waves traveling across the network.The combination of this delay operator and the adjacency matrix determines these dynamics through their effect on eigenvalues in the complex plane, thus providing mechanistic insights into the effect of heterogeneous time delays.The approach introduced here offers a mathematical description for the dynamics of time-delayed networks, an important open problem in physics [5] with many applications in neuroscience [6], engineering [7], and technology [8].In general, approaches to systems with heterogeneous time delays center on numerical simulations, and no coherent analytical approach currently exists [9,10].Importantly, while this question first arose from observations of neural dynamics in the human cortex during sleep, the delay operator we introduce here is general to studying the effect of distributed time delays in networks at finite scales, potentially allowing insight into these dynamics in a broad range of systems [11][12][13].

II. DELAY OPERATOR
We start with the standard Kuramoto model (KM) [14][15][16] and then consider the model with distancedependent time delays [17][18][19].The original KM on a general network of N nodes is defined by: where θ i ∈ [−π, π) represents the state variable (phase) of oscillator i at time t, ω i is the intrinsic angular frequency, scales the coupling strength, and A ij ∈ {0, 1} represents the elements of the adjacency matrix.The coupling of two connected oscillators i and j causes their phases to attract [15,16,20,21].
Time delays have been observed to be an important mechanism underlying the generation of traveling waves in the brain [13,[22][23][24].With this in mind, we consider a time-delay Kuramoto model (dKM) with delays τ ij that depend on the distance between two oscillators i and j: The delay operator approach we introduce here generalizes to arbitrary adjacency matrices.In order to demonstrate this approach, we start by considering an undirected ring graph G RG , where N = 100 nodes are arranged on a one-dimensional ring with periodic boundary conditions.Each node in G RG is connected to For the parameters chosen in this work, the time delays range from approximately 2 to 62 ms, a timescale relevant to neural dynamics [10,25,26].We consider the case where all oscillators have the same frequency of 10 Hz (ω = 20π); however, our approach can be applied to the case of non-identical natural frequencies [27].
The time-delay term θ j (t − τ ij ) can be approximated by θ j (t) − ωτ ij [17,18,24].Using this approximation, in combination with the algebraic approach to the Kuramoto dynamics [4,28], we introduce a delay operator, which provides analytical insight into how heterogeneous time delays can create specific, sophisticated spatiotemporal structures in the resulting nonlinear dynamics.Applying this approximation to Eq. ( 2), we arrive at an equation that captures the time-delay dynamics in the dKM in heterogeneous phase lags [17,18,24].We can then use our algebraic approach to the Kuramoto dynamics and arrive at (see Appendix -Sec.V A -for details) where x ∈ C N , and the matrix W is given by where • represents the Hadamard (elementwise) product.This matrix has information about the coupling strength , the time delays η = ωτ present in the original dKM, and the connection scheme of the system on A. In previous work, we have shown this complex-valued equation, when evaluated through the procedure described below, precisely captures the trajectories of the original, nonlinear Kuramoto model [28].We now show that this approach generalizes to the case of heterogeneous time delays.
With this approach, we have two dynamical systems: the original, nonlinear KM and a complex-valued system with the explicit solution in Eq. (3) (details on the derivation can be found in [28] and in the Appendix -Sec.V A).In the complex-valued system, x ∈ C N has elements x i (t) ∈ C whose argument we compare with the numerical solution of the original Kuramoto model with heterogeneous time delays (dKM) θ i (t) ∈ R (obtained by Euler integration of Eq. ( 2) with high temporal precision).That is, Arg[x i (t)] is compared with θ i (t).When initialized with unit-modulus initial conditions |x i (0)| = 1 for all i, with arguments Arg[x i (0)] that match the initial phases θ i (0) in the original dKM, the trajectories in the original and complex-valued KM correspond for a non-trivial window of time [4].As mentioned above, in [28] we found that iterating the explicit expression Eq. (3) in a specific manner produces trajectories in the complex-valued system that precisely match those in the original, nonlinear Kuramoto model.Specifically, we can evaluate: where ς is small but finite, t ∈ [0, ς, 2ς, • • • , nς], and Λ represents an elementwise operator mapping the modulus of each state vector element x i (t) to unity.This approach represents an iterative analytical procedure, defined by the application of the linear matrix exponential and Λ.Note that Eq. (5) propagates the solution at discrete time intervals defined by ς, Eq. ( 3) can be applied within intervals defined by ς, and ς > dt.Critically, while this iterative procedure does not represent a closed-form, all-time solution for the dynamics of the original nonlinear Kuramoto system, all evolution of the arguments Arg[x i ] (which, again, correspond with θ i (t) ∀ i in the original KM) is governed under the linear matrix exponential operator, and it is clear that the elementwise Λ operator only changes the moduli.In this work, we show that this approach applies also in the case of heterogeneous time delays and provides analytical insight into how distance-dependent time delays create specific spatiotemporal patterns.

III. RESULTS
We first study phase synchronization in networks with (dKM) and without (KM) time-delays on G RG , as a function of the coupling strength (Fig. 1).We use the Kuramoto order parameter: and its time average R for 10-second simulations to measure the level of phase synchronization.As the coupling strength increases in the non-delayed case (original KM and complex-valued KM), R begins at a low value and increases until approaching unity (representing phase synchronization).In the case with heterogeneous time delay (original dKM and complex-valued dKM), the order parameter remains low (red squares and green triangles), reflecting the fact that time delays induce a range of spatiotemporal patterns, as observed previously [17,18,[29][30][31][32][33].Here, we observe that the complex-valued model is able to capture the average dynamics that the original Kuramoto model depicts, for both the non-delayed and delayed cases, for different coupling strength across different initial conditions (Fig. 1).
We next study dynamics in the KM and dKM considering an individual realization, for a fixed coupling strength ( = 0.5), and compare the dynamics of the original dKM to the evaluation of the complex-valued approach.Without time delays, the original KM exhibits a quick transi- for the non-delayed case (blue dots: original KM, orange diamonds: complex-valued KM) and for the delayed case (red squares: original dKM, green triangles: complex-valued dKM).Each dot represents one 10-second simulation with random initial conditions (U(−π, π)), which are the same for the complex-valued case and for the numerical simulation at each point.
tion from random initial conditions to a phase synchronized state (horizontal lines, Fig. 2a).With time delays, however, phase synchronization is not reached, and the original dKM exhibits a transition from random initial conditions to a traveling wave state (diagonal structures, Fig. 2b).The evaluation of the complex-valued dKM captures both the transient dynamics and traveling wave state exhibited in the original dKM (Fig. 2c), as well as the dynamics of the Kuramoto order parameter R(t) (Fig. 2d).
Our approach to systems with heterogeneous time delay provides insight into the mechanism for these dynamics in terms of the spectrum of W -Eq. ( 4).If A and τ are circulant, W is also circulant (see Appendix -Sec.V B), hence W and A share the same eigenvectors (which form an orthonormal basis).We can then write Eq. ( 3) using the eigenspectrum of W , which results in , where α i can be written in terms of initial conditions.Importantly, we can also write Eq. ( 5) in a similar fashion, which results in where α i can again be written in terms of the state of the system at time t ∈ [0, ς, 2ς, • • • , nς].Thus, while it is in general a very difficult problem to understand the dynamics of nonlinear networks in terms of eigenspectra, this approach provides a unique insight into the connection between the spectrum of W -Eq. ( 4) -and the spatiotemporal dynamics of the nonlinear oscillator network -Eq.( 2).Critically, our approach uses familiar mathematical techniques from linear algebra matrix theory in a distinct way: while previous approaches in nonlinear dynamics have sought to describe the dynamics using the spectrum of the Laplacian matrix [34][35][36], the focus on the complex-valued system in our approach enables the insight that the argument of the eigenvectors of the matrix W provides analytical predictions about the resulting nonlinear dynamics.
Following this idea, Fig. 2e shows the eigenmode contributions, here represented by log |µ i |, as a function of time, for the dynamics in Fig. 2c.Here, the eigenmode contributions are given by the projection of the complex-valued approach solution x(t) onto the eigenvectors of W .The eigenmode contributions are obtained as µ k (t) = x(t), v k , where .denotes the standard complex inner product.Figure 2e shows that, when the network exhibits incoherent dynamics, the eigenmode contributions remain uniform across µ i .When the traveling wave pattern is reached, on the other hand, the 3 rd eigenmode becomes dominant (note the log scale).These results demonstrate that the change from incoherent dynamics to a traveling wave can be understood quite directly through the geometry of the eigenmodes.Further, in the case of circulant networks, we can evaluate eigenvalues and eigenvectors analytically using the circulant diagonalization theorem (CDT) [37], and in this case, the 1 st eigenvector represents the solution where all oscillators have the same phase (phase synchronization), and higher modes represent wave patterns, given by Fourier modes (see Appendix -Sec.V C).
The effect of heterogeneous time delays on the dynamics of the dKM can be understood through the geometry of eigenvalues in the complex plane.Figure 2f illustrates the eigenvalues of A (non-delayed) and W (delayed).While the non-delayed case (blue line and dots) has purely real eigenvalues, the effect of the heterogeneous time delays (red line and squares) can be understood in our framework in terms of the Hadamard (elementwise) product of the delay operator τ and A (see Eq. ( 4), and Appendix -Sec.V B).The effect of this operation is to provide a specific rotation of the eigenvalues in the complex plane.This rotation allows the system to access higher modes and, therefore, to exhibit different traveling wave patterns.Further, the rotation is not the same for all eigenvalues because the delays are heterogeneous.In this particular case, the rotation leads to eigenvalues associated to the 3 rd and 99 th modes to have the largest real part, allowing the system to reach traveling wave states associated with the 3 rd and 99 th modes.In the particular example of Fig. 2, the network evolves to a wave given by the 3 rd mode, but different (random) initial conditions can either evolve to the dynamics described by the 3 rd or 99 th mode [38].Moreover, when different time delays are considered, different modes can be dominant, and therefore the system evolves to a different wave pattern [39].
We can now uncover how the combination of network structure, time delays, and node state can create specific spatiotemporal patterns.By using our delay operator approach, we can analytically predict the specific pattern to which the original dKM evolves.Figure 3a shows the wave pattern given by θ obtained from the The eigenvalues of W (delayed) and A (non-delayed) provide further analytical insights into the effect of the delay in the system: it rotates the eigenvalues in the complex plane, which allows the system to access different modes.In the non-delayed case, the leading eigenvalue (in real part) is associated with an eigenvector with a zero phase difference configuration (1 st mode).In the delayed case, otherwise, there are two leading eigenvalues that are associated to the eigenvectors v3 and v99, which have phase configurations representing traveling waves.
original dKM (blue line) and the argument (elementwise) of the 3 rd eigenvector (red line), which predicts the observed dynamics [40].In this case, phases increase in the clockwise direction around the ring, which we define to be the positive direction (+1).It is important to note that, in our approach, the argument of each eigenvec- ) directly relates with the phase offset in the resulting network dynamics.Because of the correspondence between trajectories in the complex-valued model and the original dKM, this approach creates a direct link between eigenvectors of the adjacency matrix and the specific spatiotemporal dynamics that result.For the dynamics in Fig. 3a, the eigenmode contribution is given by µ 3 (see Fig. 2e), and the phase configuration matches the argument of v 3 .In the example considered here, two eigenvalues are dominant (i.e.having the largest real part) -λ 3 and λ 99 (Fig. 2f).Different initial conditions can thus evolve to the phase pattern given by the 99 th mode, which is predicted by v 99 (Fig. 3b).In this case, the spatial frequency is the same as observed in the previous case, but the direction of the wave pattern is the opposite [41].These results show a clear connection between the spectrum of the network (described by W ) and the dynamics on the original dKM, where the wave pattern (solution) can be described by the phase configuration of the eigenvector associated to the dominant mode.We take counterclockwise increases in phase to be in the negative direction, and clockwise increases to be positive.Because the network considered here has two dominant eigenvalues equal in their real parts, random initial conditions evolve equally either to the phase pattern of v 3 or v 99 in individual simulations (Fig. 3c).To quantify the spatiotemporal dynamics, the spatial frequency, and the direction of propagation, we compare the phases obtained from the original dKM and the argument of the eigenvectors of W . Specifically, we evaluate: where θ i (t) is the phase of the oscillators i at time t obtained from the original dKM, N is the number of oscillators in the network, i is the imaginary unit, and v k is the k th eigenvector of W . Here, we use v 3 , and ρ (k) = 1 means the phase configuration of the network given by the θ(t) is the same as the one given by the argument of the eigenvector v k .In the case shown in Fig. 3c, approximately half of the simulations evolve to the positive direction, indicating the dynamics matches the argument of v 3 , and approximately half evolve to the negative, indicating the dynamics is given by the argument of v 99 .
A small fraction of initial conditions exhibit inner products approximately ±0.5, corresponding to a wave with a different spatial frequency.
Using the insights from this approach, we can now design initial conditions that generate waves in a preferred direction.To do this, we started from the phase pattern specified by v 3 and randomized the phases by nearly a full cycle (0.8U[−π, π], then wrapped in [−π, π]).While this initial condition is nearly random (Fig. 3d, bottom right, where the red line represents Arg[v 3 ]; compare with Fig. 3c, bottom right), nearly all simulations evolve to the positive direction.These results demonstrate that the combination of connectivity, time delays, and network state can generate specific spatiotemporal patterns in oscillator networks -here, traveling waves with a chirality in a preferred direction.
The framework for systems with heterogeneous time delays introduced in the work generalizes to many types of networks.This approach can be applied to very sophisticated networks obtained from experimental data.In particular, this approach can successfully predict traveling wave patterns arising in an oscillator network based on connectivity in the human brain.Figure 4 illustrates simulations and the analytical prediction resulting from our approach for networks where the connectivity data is based on the Human Connectome Project (HCP) [42].In this case, N = 998 cortical regions are given at a point in 3-space, with connections between areas derived from neuroimaging data.Connection weights between regions are determined by the number of fibers [42,43], which we use to build the adjacency matrix A. Here, the coupling strength is scaled with = 200, and the initial conditions for each analysis are given by random phases [−π, π].Further, time delays are obtained by τ ij = d ij /ν, where the distances d ij are determined by the average length of these fibers, and the known axonal conduction speed is given by ν = 5 m/s [44].The dynamics of each node is represented by the Kuramoto model, given either by Eq. ( 1) in the nondelayed case and by Eq. ( 2) in the delayed case.The natural frequency of each oscillator is given by 10 Hz (simulating, for example, a specific drive from the thalamus).Using the delay operator, we construct the matrix W for these systems -Eq.( 4) -which allows us to obtain analytical predictions of the spatiotemporal patterns that emerge.
First, we consider the case without time delays, where τ ij = 0. We then obtain the eigenspectrum of the matrix W and plot the argument (elementwise) of the eigenvector associated with the leading eigenvalue (Fig. 4a).In this case, this eigenvector shows a zero phase difference across nodes, predicting phase synchronization.We then perform the numerical simulation of the Kuramoto model (without delay), given by Eq. ( 1) and plot the phase of each node in color-code (Fig. 4b), where we observe a phase synchronized behavior [45].
On the other hand, when we consider time delays in the interaction between cortical areas, the scenario is different.In this case, the argument (elementwise) of the eigenvector associated with the leading eigenvalue depicts a phase offset increasing from the bottom left to top right (in this projection), predicting a wave propagating along that direction (Fig. 4c).We then perform numerical simulations of the Kuramoto model with heterogeneous time delays -Eq.( 2) and we observe the wave pattern that is predicted by our approach, as shown in Fig. 4d [46]. .Analytical predictions of spatiotemporal patterns in brain networks.We use our approach to investigate networks based on the Human Connectome Project (HCP) [42].We consider the case without delay in the coupling between nodes and also the case with distance dependent delays (heterogeneous delay).We use our delay operator to create the matrix W , which allows analytical predictions of the dynamics.We show the phase of each node, given by the Kuramoto model, using a color-code (dark colors are values close to −π, and light colors are values close to π).(a) In the case without delay, the argument of the leading eigenvector depicts zero phase difference, which predicts phase synchronization.(b) We then study the numerical simulation for the network without delay, given by Eq. ( 1), which shows a phase synchronized state.(c) In the case with heterogeneous time delays, the argument of the leading eigenvector shows a phase offset from bottom left to the top right (in the projection), which predicts a wave pattern.(d) We perform the numerical simulation with the delayed Kuramoto model, given by Eq. ( 2), and the network depicts the wave pattern that is predicted by our approach.This shows that we are able to predict the dynamics observed in the simulations using our delay operator.
This example now clearly demonstrates the advantage of this analytical approach: when we numerically evaluate the eigenspectrum of W in this case, the leading eigenvector for the case without delays predicts phase synchrony, while the leading eigenvector for the case with delays predicts the precise wave pattern observed in the simulation.This result shows that our approach is able to predict the spatiotemporal pattern that results from connectivity and time delays in a highly relevant, realworld case.

IV. CONCLUSION
In this work, we have introduced an analytical approach to the dynamics of nonlinear oscillator networks with heterogeneous time delays, an important open problem in physics with many potential applications.The advance in this work is based on an algebraic approach to the Kuramoto model introduced in [28].Importantly, the flexibility of this framework allowed us to introduce a delay operator, which provides rigorous analytical predictions for the specific traveling wave patterns induced by distance-dependent time delays.Using this approach, we can explain the effect of time delays in terms of a rotation of the eigenvalues of the matrix describing the system, which provides a clear and precise way to understand heterogeneous time delays in terms of the geometry of eigenmodes.Our approach therefore allows analytical predictions for the specific spatiotemporal patterns exhibited by the original dKM.This framework allows us to understand how the combination of isotropic connectivity and time delays can produce traveling waves propagating in a preferred direction, as observed in experimental data [1].Importantly, while this question first arose in our study of neural dynamics in human cortex during sleep, the approach we have introduced here is general to networks of oscillators at finite scales.The results shown in this work, together with the results in [4,28], represent a coherent and general framework for nonlinear oscillator networks.
The central advance of this framework is to consider the dynamics in an individual simulation, taking into account both the initial conditions and the specific connectivity pattern in the network.This framework thus provides an opportunity to connect an individual adjacency matrix, for example a single network taken from experimental data or a single realization of a random graph model, to the specific spatiotemporal pattern that results in a simulation.This approach has important potential applications, for example in linking an experimentally reconstructed brain network to dynamics and computation in a neural system [47] or in linking the connections in a large-scale power grid to potential large and transient disruptions [48,49].In this work, we have generalized this framework to systems with heterogeneous time delays, which demonstrates the utility of this algebraic, operator-based approach to nonlinear dynamical systems at finite scales.

CODE AVAILABILITY
An open-source code repository for this work is available on GitHub: http://mullerlab.github.io.
V. APPENDIX
Based on [4,28], we introduce the complex-valued approach to the Kuramoto model described by Eq. ( 8).To do that, we introduce a new dynamical system, described by the variable ψ ∈ C: Next, multiplying both sides by i and applying Euler's formula yields i ψi (t) = iω + e −iψi(t) We define W as: where • represents the Hadamard product (or elementwise product), and η ij = ωτ ij .This results in the following matrix form of Eq. ( 10): where we note explicitly that T .Furthermore, we can write the previous equation as: Lastly, letting x(t) = e iψ(t) , we have whose general solution is x(t) = e iωt e tW x(0) .
In this work, the dynamics of the complex-valued approach is studied by considering the elementwise argument of x(t), i.e.Arg[x i (t)] ∀ i ∈ [1, N ].As shown in [28], when |xj | |xi| ≈ 1, the dynamics of Arg[x(t)] precisely matches the trajectories of the Kuramoto model given by Eq. ( 8).This allows us to use the eigenspectrum of W to understand and predict the dynamics of the Kuramoto model with heterogeneous time delays.

B. Circulant networks and Hadamard product
The definition of the Hadamard product can be described as follows: Definition 1.Let A, B be two n × n matrices.The Hadamard product A • B is a matrix of dimension n × n with elements given by For a complex number λ, we also define e •(λA) to be the matrix of dimension n × n with elements given by (e •λA ) ij = e λAij .
We have the following observation.Therefore, we conclude that both A • B and e λ•A are circulant.

C. The Circulant Diagonalization Theorem
In the case of circulant networks, we can use the circulant diagonalization theorem (CDT) to obtain the eigenspectrum of the adjacency matrix analytically [37].In our work, both the non-delayed network A and the delayed one W are circulant (see Prop. 1).The CDT states that all circulant matrices, say H = circ(h), where circ(h) is the circulant matrix constructed from the generating vector h = (h 1 , • • • , h N ), are diagonalized by the same unitary matrix U with components where k, s ∈ [1, N ] and that the N eigenvalues are given by We let Eq. (17) determine the ordering of the eigenvalues throughout this work.The argument of the eigenvectors associated with these eigenvalues correspond to the columns of the discrete Fourier transform (DFT) matrix, which range from low to high spatial frequencies.Figure 5 shows the argument of the eigenvectors in color-code.Here, Arg[(v 1 ) i ] = 0 ∀ i ∈ [1, N ] (as shown in Fig. 2), which represents zero phase difference across oscillators, or phase synchronization.The other eigenvectors represent Fourier modes (waves) with different spatial frequencies.Figure 2 shows the cases of the eigenvectors v 3 and v 99 .

Figure 1 .
Figure 1.Synchronization level for non-delayed and delayed networks.The time-average Kuramoto order parameter R is plotted as a function of the coupling strengthfor the non-delayed case (blue dots: original KM, orange diamonds: complex-valued KM) and for the delayed case (red squares: original dKM, green triangles: complex-valued dKM).Each dot represents one 10-second simulation with random initial conditions (U(−π, π)), which are the same for the complex-valued case and for the numerical simulation at each point.

Figure 2 .
Figure 2. Analytical and geometric view on the effect of time delays.The spatiotemporal dynamics of the system is represented in color-code, where the phase of each oscillator is plotted as a function of time for (a) the original KM, (b) the original dKM, (c) and the complex-valued dKM.Dark colors represent phases close to −π and light colors phases close to π.(a) Without delay, the network transitions to phase synchronization, which is represented by the horizontal lines.The effect of the delay, however, induces wave patterns in the system, whose dynamics are represented (b) in the original dKM and also captured (c) by the complex-valued model.(d) These dynamical characteristics are corroborated by the Kuramoto order parameter R(t).(e) The eigenmodes offer a geometric perspective to such dynamics, where the waves are represented by a single eigenmode contribution (3 rd mode in this case).(f) The eigenvalues of W (delayed) and A (non-delayed) provide further analytical insights into the effect of the delay in the system: it rotates the eigenvalues in the complex plane, which allows the system to access different modes.In the non-delayed case, the leading eigenvalue (in real part) is associated with an eigenvector with a zero phase difference configuration (1 st mode).In the delayed case, otherwise, there are two leading eigenvalues that are associated to the eigenvectors v3 and v99, which have phase configurations representing traveling waves.

Figure 3 .
Figure 3. Analytical predictions of specific wave patterns.(a) The phase configuration for the original dKM (blue line) matches the argument (elementwise) of the 3 rd eigenvector Arg[v3] -analytical prediction -(red dotted line).A representation on the circle using a color-code reveals the wave pattern (right).(b) Different initial conditions lead to the wave pattern that matches the argument of the 99 th eigenvector Arg[v99].These waves can propagate either counterclockwise (negative) or clockwise (positive).(c) With random initial conditions, due to the dominance of two eigenvalues (3 rd and 99 th ), the system exhibits waves propagating in both directions -with approximately half of the initial conditions evolving to each direction (top right).(d) With biased initial conditions, starting from Arg[v3] (red line, bottom right) and adding uniform random phases 0.8(U(−π, π)), we obtain a preferred direction of propagation.

Figure 4
Figure 4. Analytical predictions of spatiotemporal patterns in brain networks.We use our approach to investigate networks based on the Human Connectome Project (HCP)[42].We consider the case without delay in the coupling between nodes and also the case with distance dependent delays (heterogeneous delay).We use our delay operator to create the matrix W , which allows analytical predictions of the dynamics.We show the phase of each node, given by the Kuramoto model, using a color-code (dark colors are values close to −π, and light colors are values close to π).(a) In the case without delay, the argument of the leading eigenvector depicts zero phase difference, which predicts phase synchronization.(b) We then study the numerical simulation for the network without delay, given by Eq. (1), which shows a phase synchronized state.(c) In the case with heterogeneous time delays, the argument of the leading eigenvector shows a phase offset from bottom left to the top right (in the projection), which predicts a wave pattern.(d) We perform the numerical simulation with the delayed Kuramoto model, given by Eq. (2), and the network depicts the wave pattern that is predicted by our approach.This shows that we are able to predict the dynamics observed in the simulations using our delay operator.

Figure 5 .
Figure 5.A graphical representation of the matrix with the phase configuration of the eigenvectors of W .The k th column is the color-coded argument (elementwise) of the k th eigenvector.
k = 25 nearest neighbors in each direction, and A ij ∈ {0, 1} is 1 if oscillators i and j are connected, and 0 otherwise.The time delay τ ij = d ij /ν between two nodes i and j grows linearly with distance (d ij ) with respect to the periodic boundary conditions on the ring ( arXiv:2207.13785v2 [math.DS] 17 Jan 2023 the