Consider identical Kuramoto oscillators on a random graph. Specifically, consider Erdős–Rényi random graphs in which any two oscillators are bidirectionally coupled with unit strength, independently and at random, with probability . We say that a network is globally synchronizing if the oscillators converge to the all-in-phase synchronous state for almost all initial conditions. Is there a critical threshold for above which global synchrony is extremely likely but below which it is extremely rare? It is suspected that a critical threshold exists and is close to the so-called connectivity threshold, namely, for . Ling, Xu, and Bandeira made the first progress toward proving a result in this direction: they showed that if , then Erdős–Rényi networks of Kuramoto oscillators are globally synchronizing with high probability as . Here, we improve that result by showing that suffices. Our estimates are explicit: for example, we can say that there is more than a chance that a random network with and is globally synchronizing.
Random graphs are fascinating topologies on which to study the dynamics of coupled oscillators. Despite the statistical nature of random graphs, coherent synchronization seems to be ubiquitous above a critical threshold. To investigate this, we consider the homogeneous version of the Kuramoto model in which all oscillators have the same intrinsic frequency. For simplicity, the oscillators are arranged on an Erdős–Rényi random graph in which any two oscillators are coupled with unit strength by an undirected edge, independently at random, with some probability ; however, our proof strategy can be applied to other random graph models too. We say that a network is globally synchronizing if the oscillators converge to the all-in-phase synchronous state for almost all initial conditions. For what values of is an Erdős–Rényi random network very likely to be globally synchronizing? Here, we prove that as is a sufficient condition. Our proof uses trigonometric inequalities and an amplification argument involving the first two moments of the oscillator phase distribution that must hold for any stable phase-locked state. Specifically, we show that the spectral norms of the mean-centered adjacency and graph Laplacian matrix can be used to guarantee that a network is globally synchronizing. Our analysis is explicit, and we can reason about random networks of finite, practical size.
Networks of coupled oscillators have long been studied in biology, physics, engineering, and nonlinear dynamics.1–13 Recently, they have begun to attract the attention of other communities as well. For example, oscillator networks have been recognized as having the potential to yield “beyond-Moore’s law” computational devices14 for graph coloring,15 image segmentation,16 and approximate maximum graph cuts.17 They have also become a model problem for understanding the global convergence of gradient descent in nonlinear optimization.18 In such settings, global issues come to the fore. When performing gradient descent, for instance, one typically wants to avoid getting stuck in local minima. Conditions to enforce convergence to the global minimum then become desirable. Likewise, conditions to enforce the global synchrony of oscillators play an analogous role in dynamical systems theory.
We say that a network of oscillators globally synchronizes if it converges to a state for which all the oscillators are in phase, starting from all initial conditions except a set of measure zero. Until recently, only a few global synchronization results were known for networks of oscillators.4,5 These results were restricted to complete graphs, in which each oscillator is coupled to all the others. In the past decade, however, several advances have been made for a wider class of network structures, starting with work by Taylor,19 who proved that if each oscillator in a network of identical Kuramoto oscillators is connected to at least 94% of the others, the network will fall into perfect synchrony for almost all initial conditions, no matter what the topology of the network is like in other respects. Taylor’s result was strengthened by Ling et al.18 in 2018, and further progress has been made since then.20,21
Ling et al.18 also made a seminal advance in the study of random networks. They considered identical Kuramoto oscillators on an Erdős–Rényi random graph,22 in which any two oscillators are coupled with unit strength, independently and at random, with probability ; otherwise, the oscillators are uncoupled. They showed that if , then with high probability, the network is globally synchronizing as .
The open question is to find and prove the sharpest result along these lines. Intuitively, as one increases the value of from to , one expects to find a critical threshold above which global synchrony is extremely likely and below which it is extremely rare.
At the very least, for global synchrony to be ubiquitous, must be large enough to ensure that the random network is connected, and from the theory of random graphs,22 we know that connectedness occurs with high probability once for any . Therefore, the critical threshold for global synchronization cannot be any smaller than this connectivity threshold and is apt to be a little above. On the basis of numerical evidence, Ling et al.18 conjectured that if , then Erdős–Rényi graphs are globally synchronizing as , but nobody has proven that yet. The challenge is to see how close one can get. In this paper, we come within a factor of and prove that is sufficient to give global synchrony with high probability.
We study the homogeneous Kuramoto model19,23,24 in which each oscillator has the same frequency . By going into a rotating frame at this frequency, we can set without loss of generality. Then, phase-locked states in the original frame correspond to equilibrium states in the rotating frame. Therefore, to explore the question that concerns us, it suffices to study the following simplified system of identical Kuramoto oscillators:
where is the phase of oscillator (in the rotating frame) and the adjacency matrix is randomly generated. In particular, with probability , with otherwise, independently for . Thus, all interactions are assumed to be symmetric, equally attractive, and of unit strength. (The assumption that all interactions are equally attractive is made for convenience. We expect our arguments and results below can be generalized to attractive interactions of non-uniform strength.) We take the unusual convention that the network can have self-loops so that oscillator is connected to itself with probability ; i.e., . Since , this convention does not alter the dynamics of the oscillators, but it does make proof details easier to write down.
Finally, since the adjacency matrix is symmetric, we know that (1) is a gradient system.23,24 Thus, all the attractors of (1) are equilibrium points. Therefore, we do not need to concern ourselves with the possibility of more complicated attractors, such as limit cycles, tori, chimera states, or strange attractors. However, one slight complication is that the equilibria are not isolated; each equilibrium point lies on a subspace of equilibria because of the rotational symmetry of the Kuramoto model. Specifically, if is an equilibrium, then so is for any real number . Despite this non-generic feature, the long-term dynamics are simple. Since the dynamics of the oscillators is governed by a gradient system and the gradient is zero along these subspaces of equilibria, each solution of (1) will still converge to an equilibrium point (that happens to lie on a subspace of equilibria). In fact, since the time derivative of is zero, the average of the ’s remains constant over time with the constant being determined by the initial conditions. This means that one can assume that the dynamics of the oscillators is restricted to a hyperplane on which there is a unique all-in-phase state. When we say “the all-in-phase state,” we mean the unique one in this invariant hyperplane.
We find it helpful to visualize the oscillators on the unit circle, where oscillator is positioned at the coordinate . From this perspective, a network is globally synchronizing if starting from any initial positions (except a set of measure zero), all the oscillators eventually end up at the same point on the circle. Due to periodicity, we assume that the phases, , take values in the interval .
One cannot determine whether a network is globally synchronizing by numerical simulation of (1), as it is impossible to try all initial conditions. Of course, one can try millions of random initial conditions of the oscillators’ phases and then watch the dynamics of (1). However, even if all observed initial states eventually fall into the all-in-phase state, one cannot conclude that the system is globally synchronizing because other stable equilibria could still exist; their basins might be minuscule but could nevertheless have a positive measure. These small basins of attraction are very difficult to detect with random initial conditions.
With that caveat in mind, we note that such numerical experiments have been conducted, and they tentatively suggest that is sufficient for global synchronization.18 In this paper, we investigate global synchrony via a theoretical study. We show that is good enough to ensure global synchronization with high probability, improving on proved by Ling et al.,18 and bringing us closer to the connectivity threshold of Erdős–Rényi graphs.
Although we focus on Erdős–Rényi graphs, many of the inequalities we derive hold for any random or deterministic network. To highlight this, we state many of our findings for a general graph and a general parameter . In the end, we restrict ourselves back to Erdős–Rényi graphs and take to be the probability of a connection between any two oscillators.
Our results depend on both the adjacency matrix and the graph Laplacian matrix , where is a diagonal matrix and is the degree of vertex (counting self-loops). For any , denote the shifted adjacency and the graph Laplacian matrix by
where is the matrix of all ones and is the identity matrix. It is worth noting that for Erdős–Rényi graphs, the shifts are precisely the expectation of the matrices as and . Remarkably, we show that the global synchrony of a network can be guaranteed by ensuring that the spectral norms and satisfy particular inequalities. The spectral norm of a symmetric matrix is the maximum eigenvalue in absolute magnitude, and and are extensively studied in the random matrix literature.25,26 We also find it appealing that the spectral norm of the graph Laplacian matrix appears naturally in our analysis, as it has been used previously to study the dynamics of networks of oscillators27 as well as diffusion on graphs. From a high-level viewpoint, we prove that global synchronization of a graph is a spectral property, like many other graph properties.
While we focus in this paper on global synchronization of identical oscillators, Medvedev and his collaborators have considered other aspects of synchronization for non-identical oscillators on Erdős–Rényi graphs, as well as on other graphs, such as Cayley and Watts–Strogatz graphs, often in the continuum limit.28–31 Their findings have a similar overarching message that synchronization occurs spontaneously above a critical threshold.
II. OVERVIEW OF OUR PROOF STRATEGY
To prove that a random graph is globally synchronizing with high probability, we bridge the gap between spectral graph theory and coupled oscillator theory. The literature contains many good probabilistic estimates for the spectral norm of a random graph’s adjacency and graph Laplacian matrices, which we use to control the long-time dynamics of (1). The key to our proof is to establish conditions on these two spectral norms that force any stable equilibrium to have phases that lie within a half-circle. Confining the phases in this way then guarantees global synchronization because it is known that the only stable equilibrium of (1) with phases confined to a half-circle is the all-in-phase synchronized state.8,23
Our first attempt at controlling the phases is the inequality we state below as (9), which can be used to guarantee that very dense Erdős–Rényi graphs are globally synchronizing with high probability. As an aside, when inequality (9) together with (4) provides a new proof that a complete graph of identical Kuramoto oscillators is globally synchronizing.5 A similar inequality to (9) was derived by Ling et al.18 to show that is sufficient for global synchrony with high probability.
To improve on (9), our argument becomes more intricate. We carefully examine the possible distribution of edges between oscillators whose phases lie on different arcs of the circle and show that any equilibrium is destabilized if there are too many edges between oscillators that have disparate phases (see Lemma 8 and Fig. 2).
III. BOUNDS ON THE ORDER PARAMETER AND ITS HIGHER-ORDER MOMENTS
An important quantity in the study of Kuramoto oscillators is the so-called complex order parameter, . The magnitude of is between and and measures the synchrony of the oscillators in the network. We find it useful to also look at second-order moments of the oscillator distribution for analyzing the synchrony of random networks. Higher-order moments are also called Daido order parameters and can be used to analyze oscillators with all-to-all coupling, corresponding to a complete graph.9,32–34
For an equilibrium , we define the first- and second-order moments as
(For convenience, we use the notation to mean .) Without loss of generality, we may assume that the complex order parameter is real-valued and nonnegative. To see this, write for some . Then, is also an equilibrium of (1) with the same stability properties as since (1) is invariant under a global shift of all phases by . Therefore, for the rest of this paper, we assume that for any equilibrium of interest with . This allows us to select a representative equilibrium point from its associated subspace of equilibria. When , the oscillators are in pure synchrony, and when , the phases are scattered around the unit circle without a dominant phase.
To avoid working with complex numbers, it is convenient to consider the quantity . For , we have
By analyzing and , one hopes to witness the rough statistics of an equilibrium to understand its potential for synchrony without concern for its precise pattern of phases.
Let and note that
where is the adjacency matrix, is the complex conjugate of , and denotes the transpose of a vector. Since , we have
where is the vector of all ones and .
We would like to know when a stable equilibrium is close to the all-in-phase state, and we know this when is close to . Therefore, we begin by deriving a lower bound on for any stable equilibrium of (1). Similar inequalities involving and have been used to demonstrate that sufficiently dense Kuramoto networks are globally synchronizing.18,21
To maximize the lower bound on from Lemma 1, one can optimize over . For a random graph where each edge has a fixed probability of being present, independently of the other edges, we usually just select to be that probability. Regardless, to make the lower bound on in Lemma 1 useful, we need to find a nontrivial lower bound on since this quantity appears in the right hand side of (4). To obtain such a bound, we use the following technical lemma.
The inequality in the lemma follows by squaring the above inequalities, summing over , and noting that .
To see how Lemma 2 can be used to derive a lower bound on for any stable equilibrium state, we start by dropping the second sum in Lemma 2. Using the upper bound , we find that . Since , we have the following lower bound on :
From (4) and (7), we can deduce that when , then and must both be close to . Intuitively, this should mean that the corresponding stable equilibrium, , must be close to the all-in-phase state with the possible exception of a small number of stray oscillators. However, our goal is to prove global synchrony, which is a more stringent condition, and we must completely rule out the existence of any stray oscillators.
To precisely control the number of stray oscillators, we define a set of indices for oscillators whose phases lie outside of a sector of half-angle (centered about the all-in-phase state),
for any angle (see Fig. 1). If we can prove that is empty, then we know that all the phases in the equilibrium state lie strictly inside a half-circle. That would give us what we want because by a basic theorem for the homogeneous Kuramoto model, the only stable equilibrium with all its phases confined to a half-circle is the all-in-phase state.8,23 In terms of bounds, since is integer-valued, if we can show that , then we know that is the empty set.
From Lemma 2 and , we find that
where is if and otherwise. Therefore, by plugging in , we see that . Thus, if , then must be the empty set and the network is globally synchronizing.
Unfortunately, this kind of reasoning is not sufficient to prove global synchrony for graphs of interest to us here. For example, for an Erdős–Rényi random graph, we know that with high probability for large (see Sec. VI). Therefore, the upper bound on is approximately , which for is certainly not good enough to conclude that is empty. Instead, we must further improve our bounds on by using a recursive refinement strategy that we refer to as an “amplification” argument (see Sec. V).
IV. BOUNDS ON THE NUMBER OF EDGES AND SIZES OF SETS IN GRAPHS
The precise amplification argument that we use requires bounds on the number of edges and sizes of vertex sets of a graph expressed in terms of the spectral norms of and . It is worth noting that these bounds hold for both deterministic and random graphs. For a vertex set of a graph , we denote the characteristic vector by ; i.e., if and if . We use to denote the vector transpose of . We write to be the number of vertices in and denote the number of (directed) edges in starting at vertex set and ending at as . Therefore, is twice the number of (undirected) edges between vertices in . For Erdős–Rényi graphs, one expects to have . However, for our argument to work, expectations are not adequate; instead, we need to have bounds on the difference between and . The results in this section are proved by classical techniques, and closely related bounds are well known. We give the proofs of the precise inequalities that we need so that the paper is self-contained.
We first show that for any vertex set , the number of edges connecting a vertex in to another vertex in deviates from by at most for any .
Let be the characteristic vector for . By the min-max theorem,35 we have . Finally, we note that , , , and .
Lemma 3 controls the number of edges connecting a set of vertices. The following result controls the number of edges leaving a vertex set. We denote the vertices of that are not in by .
When a bound on is available, Lemma 4 can be used to ensure that is connected. In particular, Lemma 4 tells us that a graph of size is connected if for some .
Since , we know that . Therefore, Lemma 4 is only a useful bound when . We can take Lemma 4 a step further and bound the number of edges between any two disjoint sets of vertices of using . We denote the union of the vertex sets and by .
For any partitioning of the vertices of into three sets , , and , we have . By Lemma 4, is bounded between , is bounded between , and is bounded between . Hence, deviates from by less than .
Lemmas 3 and 5 can be combined to obtain the following result. We do not regard Theorem 6 as having significant independent interest, but later, it is an important statement in the proof of Corollary 9. It only deserves the status of a theorem in the sense that it is a key technical step in our overall argument.
V. AMPLIFICATION ARGUMENT
We are finally ready for our amplification argument, which is a way to improve the bounds on from (9). We first write down a new inequality that holds for any stable equilibrium. We write it down using a kernel function that later allows us to improve our argument with the aid of a computer (see Sec. IV D). For any stable equilibrium of (1), we know that it satisfies equilibrium and stability conditions. The equilibrium conditions are given by for all . In addition, stability requires that the Hessian matrix [see (2.3) of Ling et al.18] associated with is a nonnegative definite matrix; i.e., it has nonnegative eigenvalues.
The fact that the Hessian matrix has nonnegative eigenvalues can be used to form many useful inequalities. For example, if is the Hessian matrix, then a stable equilibrium must satisfy the Rayleigh quotient for any vector . In particular, selecting (the th unit vector) implies that the th diagonal entry of is nonnegative. This yields the inequalities for . Next, one can combine and in tricky ways and then use trigonometric identities to simplify the inequalities. After manipulating these inequalities for a while, we find that it is beneficial to massage them differently depending on the values of and that appear in them. However, then, having combined the inequalities in diverse ways, we notice that they all share a common structure. To highlight this commonality, we express the various inequalities as a unified family of the form , where is defined piecewise to reflect its dependence on its arguments.
Let be an integer between and . Due to periodicity, we may assume that the phases, , take values in the interval . We split the proof into three cases depending on the value of .
Case 1: . We first show that for all by checking the three possible subcases: (i) If , then . (ii) If , then . (iii) If , then , where the inequality holds because and . The inequality follows from the equilibrium condition .
Case 2: . For this case, we begin by showing that for all by checking the three possible subcases: (i) If , then , where the inequality holds because and . (ii) If , then . (iii) If , then . The inequality follows from the equilibrium condition .
In preliminary work, we have found indications that Lemma 7 can be used to prove stronger results than those we report below; see Sec. IV D for further discussion. However, we are not sure yet how to write down an argument that uses the full strength of Lemma 7 in a readily digestible fashion. So for now, we use the following simplified lemma instead. It is a key step in our amplification argument.
Lemma 7 shows us that if is a stable equilibrium, then any oscillator with phase in absolute value cannot be coupled to too many oscillators with phases . This is because when and . If there are too many negative contributions in the sum , then it will end up negative, which is not allowed for a stable equilibrium. To make this precise, we can use our sets in (8).
We now show that if is a stable equilibrium and is small, then must be even smaller for ; otherwise, the oscillators in the set would destabilize the equilibrium (see Fig. 2). Since we know that , the next bound is useful when .
Note that it is only possible to have an and such that in Corollary 9 when . Corollary 9 can be used in a recursive fashion to improve the bound on . Below, we start at and incrementally increase to conclude that is empty.
Finally, we summarize our findings. In particular, we can now provide a list of technical criteria that ensure that the network is globally synchronizing.
Let be a graph with vertices and . If
(ii) , and
Theorem 11 shows that a graph’s global synchrony can be ensured by the size of and alone. This is particularly beneficial for random networks as and are quantities that are studied in the random matrix literature. For random networks, hypotheses (i)–(iii) become inequalities involving random variables, namely, and . It is crucial to appreciate that these inequalities are not independent of each other, and we must be careful about this dependence when applying Theorem 11.
VI. THE GLOBAL SYNCHRONY OF ERDŐS–RÉNYI GRAPHS
The random matrix literature involves a large array of bounds on and , and navigating the literature can be daunting to outsiders. The main thing to keep in mind is that some of these bounds are asymptotic in nature and only hold for enormously large , while others are explicit (but weaker) in . We find that both kinds of bounds are useful. If we want to say something about a random graph for a given , then we need explicit bounds; whereas if we are interested in the asymptotic properties, then we may use either explicit or asymptotic bounds.
A. Explicit bounds to derive explicit statements
where . By applying the union bound to the pair of inequalities in (12), we conclude that both and hold with probability .
Next, we use this conclusion to prove that an Erdős–Rényi graph for a particular and globally synchronizes with high probability. For example, let . Then, we try a range of from and test for which all three hypotheses in Theorem 11 hold. If the three hypotheses all hold, then the graph is globally synchronizing. By following this procedure, we find that for , the corresponding Erdős–Rényi graph is globally synchronizing with probability . More specifically, when and , one can calculate that and . Hence, with probability , hypotheses (i) and (ii) in Theorem 11 hold. [In fact, hypotheses (i) and (ii) hold for much lower values of , and hypothesis (iii) is the reason that we need .]
It is at this point that we need to be careful about the lack of independence of the three hypotheses in Theorem 11. In particular, we do not check (iii) directly as this would involve messy conditional probability calculations. Instead, we find it easier to check hypothesis (iii) by checking the following deterministic inequality:
This deterministic inequality implies (iii) provided that and , which we already know hold simultaneously with probability . Next, using Theorem 11 with replaced by and on the left and right of the inequality, respectively, we conclude that an Erdős–Rényi graph with and is globally synchronizing with probability . By the same argument, we find that for , suffices.
B. Explicit bounds to derive asymptotic statements
Next, we leverage the explicit probability bounds in (12) to derive an asymptotic statement about the global synchrony of Erdős–Rényi graphs.
If we can establish that all three hypotheses in Theorem 11 hold, then the Erdős–Rényi graph is globally synchronizing. As before, by applying the union bound to the pair of inequalities in (12), we conclude that both and hold with probability . Then, by taking of the form for some and , we see that as . Hence, we find that , which shrinks to as . We conclude that hypotheses (i) and (ii) hold for any provided (12) holds. For hypothesis (iii) in Theorem 11, we again check the deterministic inequality in (13). Since for small , we see that the left hand side of (13) grows asymptotically like , while the right hand side grows like . We find that is sufficient to ensure that , and the statement of the theorem then follows by applying Theorem 11.
C. Asymptotic bounds to derive asymptotic statements
Next, we turn to the stronger asymptotic probability bounds. There are many potentially useful asymptotic bounds in the work of Füredi and Komlós,25 Vu,26 Feige and Ofek,37 and Lei and Rinaldo.38 Taken together, this literature shows that with high probability, we have (see Theorem 5.2 in Ref. 38) and (see Ref. 37). By the same argument as in the proof of Theorem 12, but with these stronger probability bounds, we arrive at our main result.
D. Using explicit bounds and optimizing with a computer
We now return to using the explicit bounds in (12) and aim to optimize our amplification argument using a computer. For a given , we find that one can significantly improve the range of for which the corresponding Erdős–Rényi graph is globally synchronizing (see Table I). Our computer program can be turned into a proof, and thus, for above the thresholds in Table I, the Erdős–Rényi networks are globally synchronizing with probability . However, writing the proofs down is unwieldy since the program works with bounds on for 1000 different values of and iteratively refines those bounds 100 000 times over.
|n .||104 .||105 .||106 .||107 .|
|n .||104 .||105 .||106 .||107 .|
By starting with and , one can alternate between (4) and (7)—in an iterative fashion—to obtain a lower bound on . The lower bound on can be substituted in (9) to give initial upper bounds on for . One can then use Corollary 9 to progressively improve the bounds on for . (Technically, we are using a slight strengthening of Corollary 9, where we redo the proof keeping all the inequalities as tight as possible. However, using Corollary 9 gives extremely similar results.) To do so, one selects and attempts to apply Corollary 9. If the application of Corollary 9 is successful, then one also has for . Since is integer-valued, any upper bound that is implies that is empty. We repeat this procedure 100 000 times to refine the upper bounds on for . If at any point we have , then we conclude that the Erdős–Rényi graph is globally synchronizing with probability . In Table I, we used the explicit value of to bound the spectral norms of and .
There are several further improvements to our computer program that we tried: (1) Using stronger probability bounds25,26,37,38 for ; (2) using Lemma 7 instead of Lemma 8, and (3) doing additional optimizations to improve the upper bounds for . For example, by selecting triples and proving a generalization of Corollary 9, we get bounds for in terms of and and bounds for in terms of and . There are similar generalizations for more than three angles. For example, when , using Lemma 8, we can only show that guarantees that an Erdős–Rényi graph is globally synchronizing with high probability, but with these extra improvements, we find that suffices. These improvements also provide good evidence—but not a proof—that Erdős–Rényi networks with globally synchronize with high probability as .
We have demonstrated how spectral properties of a graph’s adjacency and graph Laplacian matrix can be used to understand the global synchrony of a Kuramoto model with identical oscillators coupled according to a network. For Erdős–Rényi graphs, we prove that is sufficient to ensure global synchrony with high probability as . As conjectured by Ling et al.,18 we also believe that the global synchrony threshold is close to the connectivity threshold of . With the aid of a computer and Lemma 7, we have convincing evidence that Erdős–Rényi networks with are globally synchronizing with high probability as , and it is a future challenge to write down a formal proof. While Sec. VI focuses on Erdős–Rényi graphs, most of our analysis applies to any network, and we hope that it can deliver intriguing results for other random graph models.
This research was supported by the Simons Foundation grant (No. 713557) (M.K.), the NSF grants (Nos. DMS-1513179 and CCF-1522054) (S.H.S.), and the NSF grants (Nos. DMS-1818757, DMS-1952757, and DMS-2045646) (A.T.). We thank Lionel Levine and Mikael de la Salle on MathOverFlow for references on bounding and for Erdős–Rényi graphs.
Conflict of Interest
The authors have no conflicts to disclose.
Martin Kassabov: Conceptualization (equal); Formal analysis (equal); Methodology (equal); Software (equal); Writing – original draft (equal); Writing – review & editing (equal). Steven H. Strogatz: Conceptualization (equal); Formal analysis (equal); Methodology (equal); Software (equal); Writing – original draft (equal); Writing – review & editing (equal). Alex Townsend: Conceptualization (equal); Formal analysis (equal); Methodology (equal); Software (equal); Writing – original draft (equal); Writing – review & editing (equal).
The data that support the findings of this study are available within the article.