## Abstract

Understanding the effect of spike-timing-dependent plasticity (STDP) is key to elucidating how neural networks change over long timescales and to design interventions aimed at modulating such networks in neurological disorders. However, progress is restricted by the significant computational cost associated with simulating neural network models with STDP and by the lack of low-dimensional description that could provide analytical insights. Phase-difference-dependent plasticity (PDDP) rules approximate STDP in phase oscillator networks, which prescribe synaptic changes based on phase differences of neuron pairs rather than differences in spike timing. Here we construct mean-field approximations for phase oscillator networks with STDP to describe part of the phase space for this very high-dimensional system. We first show that single-harmonic PDDP rules can approximate a simple form of symmetric STDP, while multiharmonic rules are required to accurately approximate causal STDP. We then derive exact expressions for the evolution of the average PDDP coupling weight in terms of network synchrony. For adaptive networks of Kuramoto oscillators that form clusters, we formulate a family of low-dimensional descriptions based on the mean-field dynamics of each cluster and average coupling weights between and within clusters. Finally, we show that such a two-cluster mean-field model can be fitted to synthetic data to provide a low-dimensional approximation of a full adaptive network with symmetric STDP. Our framework represents a step toward a low-dimensional description of adaptive networks with STDP, and could for example inform the development of new therapies aimed at maximizing the long-lasting effects of brain stimulation.

## 1 Introduction

Synaptic plasticity is considered the primary mechanism for learning and memory consolidation. Neurons with similar activity patterns strengthen their synaptic connections, while others connections may weaken. Spike-timing-dependent plasticity (STDP) has been suggested as an unsupervised, local learning rule in neural networks (Gerstner et al., 1996; Song et al., 2000) motivated by experimental findings (Markram et al., 1997; Bi & Poo, 1998; Feldman, 2000; Froemke & Dan, 2002; Cassenaer & Laurent, 2007; Sgritta et al., 2017). These experimental studies reported synaptic strengthening or weakening depending on the order and timing of pre- and postsynaptic spikes. In causal STDP (see Figure 1A), long-term potentiation (LTP) occurs when the postsynaptic neuron fires shortly after the presynaptic neuron. Conversely, long-term depression (LTD) occurs when the postsynaptic neuron fires shortly before the presynaptic neuron. The closer the spike times are, the larger the effect. Noncausal, symmetric STDP has also been reported in the hippocampus (Abbott & Nelson, 2000; Mishra et al., 2016; see Figure 1B). Although Donald Hebb emphasized the importance of causality in his theory of adaptive synaptic connections in 1949 (Hebb, 2005), such noncausal rules, which neglect temporal precedence, are sometimes referred to as Hebbian plasticity rules.^{1}

As a major contributor to long-term plasticity (timescale of 1 s or longer), STDP is key to long-term neural processes in the healthy brain, such as memory (Litwin-Kumar & Doiron, 2014) and sensory encoding (Coulon et al., 2011), but also to modulate networks affected by neurological disorders (Madadi Asl et al., 2022). In particular, long-term plasticity is critical to the design of effective therapies for neurological disorders based on invasive and noninvasive brain stimulation. Paired associative stimulation using transcranial magnetic stimulation (TMS) has been shown to trigger STDP-like changes (Müller-Dahlhaus et al., 2010; Johnen et al., 2015; Wiratman et al., 2022), and STDP models have been used to design electrical stimulation for stroke rehabilitation (Kim et al., 2021). Coordinated reset deep brain stimulation (DBS) for Parkinson’s disease was designed to induce long-term plastic changes outlasting stimulation using STDP models (Tass & Majtanik, 2006) and was later validated in nonhuman primates (Tass et al., 2012; Wang et al., 2016) and patients (Adamchic et al., 2014).

However, progress in these areas is restricted by the significant computational cost associated with simulating neural network models with STDP and by the lack of low-dimensional description of such networks. Indeed, simulating networks with STDP of even moderate size for more than a couple of minutes of biological time can be impractical because the number of weight updates scales with the square of the network size. Low-dimensional mean-field approximations have been developed for networks with short-term plasticity (Tsodyks et al., 1998; Taher et al., 2020; Gast et al., 2020; Schmutz et al., 2020; Gast et al., 2021), but reductions for STDP have assumed that plasticity does not change firing rates and spike covariances (Ocker et al., 2015) or that the network is in a balanced state at every point in time (Akil et al., 2021). Even if node dynamics are given by a simple Kuramoto-type model, it is a challenge to understand the dynamics of large networks with adaptivity. While the continuum limit of Kuramoto-type networks with STDP-like adaptivity can be described by integro-differential equations from a theoretical perspective (Gkogkas et al., 2021), these do not necessarily elucidate the resulting network dynamics or yield a computational advantage. Approaches like the Ott-Antonsen reduction (Ott & Antonsen, 2008), which have been instrumental to deriving low-dimensional descriptions of phase oscillators (see Bick et al., 2020) are not directly applicable to adaptive networks where all connection weights evolve independently of one another. Indeed, one would not expect an exact mean-field description with only a few degrees of freedom to be possible (as for the Kuramoto model with static connectivity) without making further assumptions, as any exact low-dimensional description would have to reflect the high, adaptivity-induced multistability (Berner et al., 2019).

Here, we construct mean-field approximations for coupled Kuramoto phase oscillators subject to PDDP and compare these approximations to fully adaptive networks where every edge evolves according to STDP. Specifically, we consider two types of PDDP rules: rules that update connection weights continuously, as was done in previous studies, and event-based rules, where weights are updated according to phase differences at a particular phase corresponding to spiking. In section 2, we show that single-harmonic PDDP rules can indeed approximate symmetric STDP in adaptively coupled networks of Kuramoto oscillators, while a multiharmonic rule is required to accurately approximate causal STDP. For PDDP rules, we derive exact equations describing the evolution of the mean coupling strength in terms of the Kuramoto/Daido order parameters that encode synchrony of the oscillators’ phases (see section 3). We then focus on networks with symmetric adaptive coupling that naturally form clusters (see section 2). For such networks, we construct mean-field approximations (see section 4) for the emergent coupling topologies, where each cluster corresponds to a coupled population. If we assume that coupling between clusters is through the mean coupling strength—rather than by individual weights between oscillators—we obtain low-dimensional Ott-Antonsen equations for the mean-field limit. We explicitly analyze the dynamics of the reduced equations for adaptive networks for one and two clusters. Note that these mean-field descriptions are not valid globally (i.e., there is no single reduced equation that is valid on all of state space) but rather aim to capture the dynamics on part of overall phase space determined by the initial conditions. In other words, we have a family of low-dimensional dynamics that can describe part of the phase space for this very high-dimensional system. Finally, we show that the dynamics of the full network can be approximated by such a family of low-dimensional dynamics by extracting the mean field description from the emergent clustering (see section 5). Since brain activity is transient, we focus on transients and consider additive plasticity without bounds on individual weights for causal STDP. In line with previous studies, however, we include a weight decay term when considering symmetric STDP (Seliger et al., 2002; Berner et al., 2019).

## 2 PDDP Can Approximate STDP in Kuramoto Networks

Adaptation of network connections through STDP rely—as the name suggests—on the timing of action potentials of the coupled neurons. If the state of each neuron can be described by a single phase variable (e.g., if the coupling is weak; Ashwin et al., 2016), then it may be possible to approximate STDP by an adaptation rule that depends on the phase differences between oscillators such as equation 1.3. In this section, we now consider general PDDP rules, which can update weights continuously (as in equation 1.3) or update at discrete time points (spiking events). We show that for a network of phase oscillators, these PDDP rules can approximate both symmetric and causal STDP. For causal STPD, the accuracy increases substantially as the number of harmonics in the PDDP rule is increased.

### 2.1 Principles to Convert STDP to PDDP

PDDP rules prescribe synaptic changes based on differences in phase of neuron pairs rather than differences in spike timing and can be used to approximate both symmetric and causal STDP in phase oscillator networks. In particular, the approximation is expected to hold under the assumption that the evolution of phase differences is slower than the phase dynamics (Lücken et al., 2016). Under this assumption, spike time differences are approximated by dividing phase differences by the mean angular frequency of the network $\Omega =1N\u2211k=1N\omega k$.

As the phase difference is continuous in time, the discrete weight updates, based on spike-time differences in the case of STDP, can be converted to a continuous-time differential equation in terms of the phase differences. As a result, the coupling weight between each pair of neurons updates continuously based on the phase difference between the pre- and postsynaptic oscillators; we refer to this as continuously updating PDDP or simply PDDP when there is no ambiguity. STDP updates occur every time a neuron spikes, while PDDP updates occur continuously at every point in time. To ensure that STDP and continuously updating PDDP scale similarly, we scale the discrete STDP updates by the average number of spikes per unit time $\Omega /2\pi $.

Rather than updating weights continuously, we can restrict weight updates to occur only at spiking events. We say that oscillator $k$ spikes if its phase increases through $\theta k=0$, and let $tkq$ be the $qth$ firing time of neuron $k$. At each spiking event, we update the coupling weight between each pair of neurons based on their phase difference; we give explicit examples of the functional form of the updates below. We refer to this type of PDDP rule as event-based PDDP (ebPDDP). While there is an explicit phase dependence through the events, the actual change only depends on the phase difference. To ensure appropriate scaling, we again multiply by the average number of spikes per unit time and introduce an additional factor that is only nonzero when either the pre- or postsynaptic neuron spikes. This factor is defined as $C=\u2211q\delta (t-tkq)+\u2211q\delta (t-tlq)/2$, where $\delta $ denotes the Dirac delta function.^{2}

### 2.2 Symmetric STDP and Single-Harmonic PDDP

As shown in Figure 3, both single-harmonic PDDP rules (the Seliger rule and equation 2.2) can approximate symmetric STDP (equation 1.2) in networks of Kuramoto oscillators. Simulation details can be found in appendix A.1. The time evolution of the weight distribution (panel A), the coupling matrix at the last simulation time point (panel B), as well as the the time evolution of the average coupling (panel C1) and network synchrony (panels C2–C3) are comparable across learning rules. In particular, the Pearson’s correlation between the STDP coupling matrix and the PDDP coupling matrix is 0.88, while the correlation between the STDP coupling matrix and the ebPDDP coupling matrix is 0.90. For the parameter set shown in Figure 3, we see the emergence of two synchronized clusters. For this type of dynamics, the average coupling as a function of time may be better captured by ebPDDP than PDDP (panel C1). However, for the desynchronized state (shown in Figure 11 in appendix C), there is little difference between continuous PDDP and ebPDDP. In both states, the accuracy of the approximation could be improved by considering a Fourier expansion of the Mexican hat function rather than a single cosine term. We explore this for causal, nonsymmetric STDP in the next section, as we found that a single sine term is generally a poor approximation of the causal STDP kernel.

### 2.3 Causal STDP and Multiharmonic PDDP

In the previous section, we considered PDDP rules with a single harmonic in the phase difference. To get a better approximation of causal STDP, one can take more harmonics into account.

#### 2.3.1 Obtaining Multiharmonic PDDP Rules from Causal STDP

#### 2.3.2 Comparison of Learning Rules

Multiharmonic PDDP can approximate causal STDP in Kuramoto networks (simulation details can be found in appendix A.2). As the number of harmonics included in the PDDP rules is increased, the dynamics for the networks with PDDP begin to match those of the STDP network (see Figure 5). The time evolution of the weight distribution (panel A), the coupling matrix at the last simulation time point (panel B), as well as the the time evolution of the average coupling (panel C1) and network synchrony (panels C2 to C7) are closely matched for causal STDP and multiharmonic PDDP or ebPDDP when enough Fourier components are included ($Nf=25$). Simulations were also performed for a range of parameter values, and the findings were similar (see Figures 12 to 15 in appendix C). Single-harmonic PDDP or ebPDDP ($Nf=1$) can provide a first level of approximation of causal STDP in certain cases when the time evolution of the weight distribution is simple (see Figure 14). However, single-harmonic rules are unable to describe more complex cases, even qualitatively (as seen in Figure 5). In all cases studied for a network frequency of 5 Hz, 25 Fourier components are deemed sufficient to approximate the dynamics of the network with causal STDP. While the performance of ebPDDP is similar to PDDP, ebPDDP is slightly more accurate than PDDP. To quantitatively compare STDP to PDDP and ebPDDP, we construct error metrics for the time evolution of the weight distribution $ehist(\kappa kl)$, the average coupling $e\kappa ^$, and the network synchrony $e\rho $, and consider the Pearson’s correlation between coupling matrices at the last simulation time point $r\kappa kl\u221e$. These metrics are defined in appendix A.2. In general, these metrics improve with increasing $Nf$ (see Figure 6A), although stagnation or a slight worsening can be seen when the error is already low. Though it is of no consequence in Kuramoto networks with sine coupling, self-coupling weights are consistently different between causal STDP and multiharmonic PDDP or ebPDDP in our simulations. This is due to the truncated Fourier expansions of $F$ not being zero when the phase difference is zero (see Figure 4B3).

#### 2.3.3 Parameter Dependence

Model parameters influence the four metrics described in the previous section (see Figure 6B), although the impact on $r\kappa kl\u221e$ is minor (always stays $>0.96$). Increasing the standard deviation of the frequency distribution ($\Delta $) tends to improve the error metrics, except $r\kappa kl\u221e$, which gets slightly lower. This is expected since a larger $\Delta $ reduces synchrony and makes the weight distribution more unimodal. The impact of the standard deviation of the initial coupling distribution ($\sigma \kappa $) on the metrics is smaller, except for the time evolution of the weight distribution. The effect of the ratio of the scales of LTD to LTP ($\beta =A-/A+$) depends on the metric considered. The largest effect is a lowering of $e\rho $ when LTD dominates ($\beta =1$) compared to the balanced situation ($\beta =0.5$). This is due to the fact that the time evolution of $\rho $ is closely approximated with $Nf=1$ when LTD dominates, whereas in the balanced state, more Fourier components are required to obtain a good approximation. Since dominant LTD leads to lower synchrony, this matches the previous observation that lower synchrony is associated with lower error metrics. Time courses of $\rho $ for both states can be found in appendix C; Figure 12 shows the balanced state, and Figure 14 shows the LTD dominant regime.

To test the robustness of our findings, we studied a network with a mean frequency four times higher (20 Hz) and considered the least favorable part of parameter space (lowest $\Delta $, lowest $\sigma \kappa $, and $\beta =0.5$) (see Figure 15 in appendix C). The order of magnitude of the error metrics is the same as for 5 Hz, except for $r\kappa kl\u221e$ (see Figure 18C in appendix C). The lower value for $r\kappa kl\u221e$ may be explained by the greater complexity and finer structures in the coupling matrix. However, the coupling matrix obtained at the last stimulation point with multiharmonic ebPDDP for $Nf=75$ (see Figure 15 in appendix C, panel B7) is qualitatively very similar to the coupling matrix obtained with STDP (ses panel B1). At 20 Hz, the period of the oscillators is comparable to the STDP time constants; hence, the weight distribution patterns unfolding in time are more highly multimodal than at 5 Hz. These patterns are still well approximated by multiharmonic PDDP and ebPDDP, but a larger number of Fourier components than at 5 Hz is warranted for accurate results. Simulating this higher-frequency network required a smaller simulation time step ($\Delta t=0.1$ ms). However, the time step has overall little impact on the error metrics at lower frequencies (see Figure 18 in appendix C, panel B).

## 3 Evolution of the Average Coupling Strength in Networks with PDDP

### 3.1 Average Coupling for General PDDP

### 3.2 Average Coupling for General Event-Based PDDP

Simulations show that average weights obtained from equation 3.6 are similar to average weights obtained from equation 2.4 as shown in Figure 7. Although equation 3.6 is an exact description of the average weight in phase oscillator networks with adaptivity given by equation 3.5 (where the Dirac deltas are functions of phase), our simulations are based on equation 2.4 (where the Dirac deltas are functions of time). The approximation $\u2211q\delta (t-tkq)\u2248\Omega \delta (\theta k)$ used to derive equation 3.5 from equation 2.4 gives rise to the discrepancies visible in Figure 7.

## 4 Mean-Field Dynamics of Oscillator Populations with Adaptive Coupling

### 4.1 Low-Dimensional Dynamics for Homogeneous Coupling

which describes homogeneously coupled populations.

### 4.2 Single Harmonic PDDP, One Population

^{3}determined by

There is a trivial fixed point at $\kappa ^=0$, $\rho =0$. While $\kappa ^=\lambda \rho 2$, $\rho =1-2\Delta \kappa ^$ defines a pair of nontrivial fixed points, which exist for $\lambda >8\Delta $. Computing the Jacobian, we find that the trivial solution is stable for all physical parameter values ($\Delta >0$, $\epsilon >0$), and for the nontrivial fixed points, one is stable and the other unstable.

Using XPPAUT (Ermentrout, 2002), we performed a one-parameter continuation in the heterogeneity parameter $\Delta $ (see Figure 8). When the heterogeneity is low, there exists a nontrivial stable fixed point, where the mean coupling does not decay to zero, whereas after the saddle-node bifurcation at $\Delta =0.125$, the mean coupling will always decay to zero and the oscillators will be asynchronous ($\rho =0$). As the strength of the plasticity rule $\lambda $ is increased, the saddle node moves to the right, and the region of bistability (where the trivial and non-trivial fixed points coexist) increases. These results agree with analytical results outlined above.

### 4.3 Single Harmonic PDDP, Two Populations

Setting $q=0.5$ (two equally sized populations), we perform a one-parameter continuation in the intrinsic frequency difference $\Delta \Omega $ (see Figure 9A). We find that for a small range of $\Delta \Omega $ values ($\Delta \Omega \u2208[-0.23,0.23]$) the two populations synchronize in frequency. This frequency-locked solution lies on an invariant set, where the level of synchrony of each population is identical ($\rho 1=\rho 2$) and the coupling strengths are symmetric ($\kappa ^11=\kappa ^22$, $\kappa ^12=\kappa ^21$). The system itself, however, is generally not symmetric (unless $\Delta \Omega =0$) due to distinct intrinsic frequencies resulting in distinct mean phases of the two populations. The two nontrivial equilibrium branches in the interpopulation coupling strength (see Figure 9, panel A3), correspond to the in-phase/symmetric solution (positive $\kappa ^\mu \nu $) and the antiphase/asymmetric solution (negative $\kappa ^\mu \nu $).

For general population sizes $q$, an additional branch of solution can emerge, where one population is fully asynchronous and the other is partially synchronized with nontrivial dynamics (e.g., $\rho 1=0$, $\rho 2\u22600$); we call this the *decoupled solution*. More specifically, the structure of the equations implies that $\rho 1=\kappa ^12=\kappa ^12=0$ defines a dynamically invariant set. On this set, the coupling strength $\kappa ^11$ decays exponentially, and the second population evolves independent of the first one (as the network is decoupled) according to the equations of motion for a single population with adaptive coupling. The decoupled solution now corresponds to the nontrivial solution branch for a single population (see Figure 8), which exists for $\lambda >8\Delta $. The only difference here is that the coupling $\kappa ^$ is replaced by an effective coupling $q\mu \kappa ^\mu \mu $ scaled with the (relative) population size. Hence, decoupled solutions exist only for $q<2\Delta /\lambda $ (population 2 has nontrivial dynamics $\rho 2\u22600$) and $q>8\Delta /\lambda $ (population 1 has nontrivial dynamics $\rho 1\u22600$). Note that these are equilibrium points for the effective dynamics but periodic solutions in the full system.

Letting $q=0.1$, we find and continue the decoupled solution where population 1 is asynchronous and population 2 has nontrivial dynamics ($\rho 1=0$, $\rho 2\u22600$) (see Figure 9B). The dynamics of population 1 (population 2) is shown in red (black). The interpopulation coupling strengths are equal for all values of $\Delta \Omega $. Hence, the curve in Figure 9, panel B3 corresponds to $\kappa ^12$ and $\kappa ^21$. The decoupled solution goes unstable at a torus bifurcation (blue stars) $\Delta \Omega \u2248-0.3$ and restabilizes at a second torus bifurcation at $\Delta \Omega \u22480.3$. The frequency-locked solution branches off the decoupled solution at a pitchfork bifurcation at $\Delta \Omega \u2248-0.15$ and $\Delta \Omega \u22480.15$. As in the $q=0.5$ case, for the frequency-locked solution, the two populations have nonzero within-population synchrony and intra-/interpopulation coupling strengths. However, the within-population synchrony and the intrapopulation coupling strength are no longer equal ($\rho 1\u2260\rho 2$, $\kappa ^11\u2260\kappa ^22$). We note that the trivial solution ($\rho \mu =0$, $\kappa ^\mu \nu =0$), which is stable for all values of $\Delta \Omega $, is not shown in Figure 9B, as it would have obscured the unstable region of the decoupled solution.

Finally, we performed a two-parameter continuation in $\Delta \Omega $ and $q$ (see Figure 9C). The saddle-node curves, which demarcate the region of existence for the coupled solution, are shown in green. The coupled solution exists between the two green curves. The orange curve corresponds to the pitchfork bifurcation, where the frequency-locked solution branches off the decoupled solution. There is a global bifurcation at $q=0.2=2\Delta /\lambda $ and $q=0.8=8\Delta /\lambda $, when the decoupled solution ceases to exist, and so there is no longer a pitchfork bifurcation. The torus bifurcation curve, where the decoupled solution changes stability, is plotted in blue. In the region between the torus bifurcation curve and the saddle-node bifurcation curve, the two populations have nontrivial dynamics. The interpopulation coupling is weak enough that the populations do not entrain. Hence, the mean phases precess at different rates, and as such, the phase difference oscillates in time. As a result, the coupling strengths oscillate as the populations move in- and out-of-phase with each other. Given that the coupling strengths and the synchrony are interdependent, the synchrony variables also oscillate in time.

## 5 Describing the Full Adaptive Network Using Coupled Populations Evolving according to Mean-Field Dynamics

In this section, we aim to approximate a network of Kuramoto oscillators (see equation 2.1) with adaptivity given by symmetric PDDP (see equation 1.3 with $\varphi =0$) using two coupled populations evolving according to the mean-field dynamics. As shown in section 2.2, the full adaptive network with symmetric PDDP is itself a good approximation of the full adaptive network with symmetric STDP. Here, we optimize parameters of the coupled mean-field equations to best approximate the full network.

### 5.1 Optimizing the Parameters of the Two-Population Mean-Field to Approximate the Full Adaptive Network

To describe the full adaptive network using the two-population mean-field approximation, equation 5.1, we optimize $\epsilon \mu $ and $\lambda \mu $ with $\mu \u2208{1,2}$, as well as the proportion of oscillators allocated to the first population denoted by $q1$. Note that adaptivity parameters within and between populations ($\lambda 1$, $\lambda 2$, $\epsilon 1$, $\epsilon 2$) are not the same as adaptivity parameters between oscillators in the full network ($\lambda $ and $\epsilon $). The optimal proportion of oscillators to allocate to each population is also unclear; we therefore optimize the allocation of oscillators as follows. We sort oscillators by the average value of their mean outgoing coupling (i.e., $\u2211k=1N\kappa kl/N$) and allocate the first $100\xd7q1$% to the first population and the rest to the second population ($q2=1-q1$). This results in populations of size $N1$ and $N2$, and initial conditions for equation 5.1 can be obtained from the initial conditions used to simulate the full adaptive network. For each initial condition, the parameters $\Omega \mu $ and $\Delta \mu $ are obtained as the median frequency and half of the interquartile range of the frequency of oscillators in each population, respectively. The optimization is performed as a sweep over $q1={0.1,0.2,0.3,\u2026,0.9}$, and for each value of $q1$, 36 local optimizations over the remaining parameters are carried out. Each local optimization starts from a random set of parameters and consists of successive optimizations using patternsearch and fminsearch (MatlabR2021a) over all initial conditions in the training set.

### 5.2 Performance on Training and Test Sets

The optimized two-population mean-field approximation with the lowest cost was found for $q1=0.6$ (dark blue in Figure 10), and its parameters are given in Table 1. To test the model performance, we construct three control models; (1) a constant model defined from initial conditions by $\rho =\rho 0$ and $\kappa ^=\kappa ^0$, (2) a two-population mean-field approximation with $q1=0.6$ as in the fully optimized model but without adaptivity ($\epsilon 1=\epsilon 2=0$), and (3) a two-population mean-field approximation where $\lambda 1\lambda 2=25$ and $\epsilon 1\epsilon 2=0.5$ are chosen to match $\lambda $ and $\epsilon $, respectively. For the third control model, we performed a sweep over $q1$ and found the best fit for $q1=0.95$.

The fully optimized two-population mean-field approximation with the lowest cost on the training set is better at describing the dynamics of $\rho $ and $\kappa ^$ on both the training set and the test set than the controls shown in Figure 10 in appendix C. This shows that to best reproduce synthetic data from the full adaptive model, it is necessary to consider both the evolution of the order parameter and of the average weights and to optimize the within- and between-population adaptivity parameters. The constant model is shown in light orange, the two-population mean-field approximation with no adaptivity is shown in dark orange, and the two-population mean-field approximation where the plasticity parameters are not optimized is in light blue. In the test set, the greatest improvement of the optimized two-population mean-field approximation over the controls is for $\kappa ^$ (Figure 10D). With the exception of the constant model, the controls already approximate $\rho $ well in the test set, see Figure 10C. Representative trajectories of the two-population mean-field approximation obtained from initial conditions in the test set approximate the network synchrony, phase, and average coupling of synthetic data from the full adaptive network (see Figure 20 in appendix C). The approximation of $\kappa ^$ is overall better when $\epsilon \mu $ and $\lambda \mu $ are optimized. We note that while transients were not included in the cost function used to fit to synthetic data, they are relatively well described by the model when starting from random connectivity with different means in the test set (see Figure 20).

## 6 Discussion

Using simulations, we showed that PDDP and ebPDDP can provide useful approximations of STDP. In particular, single-harmonic rules can approximate simple forms of symmetric STDP, while multiharmonic rules are required to accurately approximate causal STDP. In the latter case, the accuracy of the approximation increases with the number of Fourier coefficients before reaching a plateau. The Fourier coefficients can be easily computed using analytical expressions only involving causal STDP parameters (see equation 2.6). We found ebPDDP to be a slightly better approximation than PDDP in the case of both symmetric STDP and causal STDP. This is expected since contrary to PDDP, ebPDDP restricts synaptic weight updates to spiking events, which is conceptually closer to STDP. One limitation of approximating STDP using plasticity rules based on phase difference is that the evolution of phase differences is required to be slower than the phase dynamics (Lücken et al., 2016). Under this assumption, the intricate evolution of coupling weights can be well approximated by PDDP and ebPDDP (see Figures 11 and 15 in appendix C). For best results, the STDP function should also decay to zero faster than the network mean frequency for both positive and negative spike-timing differences (in the case of causal STDP, to avoid both components in equation 2.3 strongly overlapping).

We derived exact expressions for the evolution of the average coupling weight in networks of phase oscillators with PDDP and ebPDDP. These expressions make no assumption about the underlying network (in particular, no assumption about the strength and type of coupling between oscillators), and are compatible with any plasticity rule based on phase difference as long as it can be expanded as a Fourier series of the phase difference. As a proof of principle, we focused on mean-field approximations based on the average coupling evolution for two-cluster states in a population of adaptive Kuramoto oscillators and performed a bifurcation analysis to highlight the different possible behaviors. Here we have focused predominantly on the Kuramoto model with adaptivity rules that mimic the adaptation between individual neural cells. However, as a prototypical model to study synchronization in neuroscience, the Kuramoto model has also found application on different scales including whole-brain modeling, epilepsy, and Parkinson’s disease (Cumin & Unsworth, 2007; Breakspear et al., 2010; Cabral et al., 2014; Schmidt et al., 2015; Ponce-Alvarez et al., 2015; Finger et al., 2016; Asllani et al., 2018; Weerasinghe et al., 2019, 2021; Bick et al., 2020; Duchet et al., 2022). While here each Kuramoto oscillator typically represents neural populations rather than individual cells, our results may still help understand the role of adaptivity in such networks—albeit with potentially different adaptivity rules on different timescales.

Our framework could easily be adapted to consider more biologically realistic oscillator models, such as the $\theta $-neuron model or the formally equivalent quadratic integrate-and-fire (QIF) neuron model. Like the Kuramoto model, the $\theta $-neuron model is amenable to the Ott-Antonsen ansatz and as such, is amenable to exact mean-field description (Luke et al., 2013). For the QIF model, the equivalent Lorentzian ansatz should be used (Montbrió et al., 2015). The model analysis would be identical, but given the explicit phase dependency in the coupling, we can expect a richer set of dynamics for a network of $\theta $-neurons than observed here for Kuramoto oscillators. We could also include explicit synaptic variables, as in Byrne et al. (2020), to more accurately model synaptic processing. More generally, under the assumption of weak coupling, even more detailed neuron models, such as the Hodgkin-Huxley model, can be approximated by networks of phase oscillators through phase reduction (Brown et al., 2004; Pietras & Daffertshofer, 2019). While the resulting phase equations will generically contain higher-order terms that affect the global dynamics, a truncation to first order (i.e., Kuramoto-type coupling) can still provide suitable approximation where the first harmonics are dominant and shape the dynamics. While generic higher-order terms may break the Ott-Antonsen ansatz, its extensions (see, e.g., Vlasov et al., 2016; Tyulkina et al., 2018) or more general moment closure approaches (Kuehn, 2016) can provide suitable approximations that do not rely on a single order parameter to capture clustering; the dynamics of the mean coupling strength may then depend on more than a single order parameter. Indeed, previous studies of adaptive oscillators (Berner et al., 2019) and our full network simulations point to the existence of three-, four-, and five-cluster states that motivate extending the results presented here. Although computationally expensive and, perhaps, numerically challenging, both the bifurcation analysis and model fitting could be extended to consider more than two clusters.

Combining theoretical insight and data-driven inference, we fitted a two-cluster mean-field approximation to a full adaptive network of Kuramoto oscillators in order to obtain a low-dimensional representation of the full system. The goal here is not to find a relationship between model parameters and (microscopic) parameters that could be measured in the brain (which is hardly ever possible with limited data), but to find a mean-field model with plasticity that can be used to describe population-level recordings, such as local field potentials. While the two-cluster mean-field model can approximate the full adaptive network on the test set, the accuracy of the approximation could be improved in several ways. First, a richer training set could be used. Second, the mean-field description of the order parameter is not exact, and a moment or cumulant approach may capture more of the full system dynamics (Tyulkina et al., 2018). Third, considering more than two populations may provide a more accurate approximation. Fourth, rather than allocating oscillators to clusters by thresholding their mean outgoing coupling, a finer partition could be learned (Snyder et al., 2020). Nevertheless, data-driven inference of low-dimensional representations of phase-oscillator networks (Thiem et al., 2020; Snyder et al., 2020; Fialkowski et al., 2022) is a promising approach to approximate the behavior of networks with STDP. Once clusters have been identified, a very recent mean-field technique based on the collective coordinate method can be used (Fialkowski et al., 2022).

More general types of STDP call for extensions of our framework. First, one would naturally expect that the strength of plastic connections is bounded. We included soft bounds in our investigation of symmetric STDP through a damping term (see equation 1.3), which is conserved in the average coupling (see equation 3.4), and is therefore present in the mean-field analyses carried out in sections 4 and 5. While neither soft nor hard bounds can easily be included in the average coupling derivation for causal STDP, we show that short transients of causal STDP with hard bounds are well captured by causal PPDP/ebPPD without bounds (see Figures 16 and 17 in appendix C). Second, we primarily focused on symmetric STDP rules, that is, interchanging the order of spikes will have little to no effect on the change of connection strength. For such adaptation, the network naturally forms clusters of strongly connected units (see Figure 3). By contrast, the causality of traditional asymmetric STDP rules will be reflected in the network structure as shown in Figure 5 and highlighted very recently in Thiele et al. (2023). To analyze the mean-field dynamics of such networks, a natural approach would be to computationally identify emerging feedforward structures in such networks instead of looking for clusters. In this case, finding a corresponding low-dimensional description as in section 4 is more challenging. Third, we consider STDP adaptation rules that depend on pairs of oscillator states. More elaborate, spike-based rules such as triplet interactions have recently attracted attention (Pfister & Gerstner, 2006; Montangie et al., 2020). It would be interesting to have such adaptation reflected in the STDP model as these could be interpreted as “higher-order” network effects (see Bick et al., 2021).

As a step toward low-dimensional description of adaptive networks with STDP, our framework has implications for the study of long-term neural processes. In particular, the effects of clinically available DBS quickly disappear when stimulation is turned off; thus, stimulation needs to be provided continuously. To spare physiological activity as much as possible, it would therefore be highly desirable to design stimuli aimed at eliciting long-lasting effects. Continuous stimulation can also lead to habituation, where stimulation benefits diminish considerably over the years in some patients with essential tremor (Fasano & Helmich, 2019). Optimizing brain stimulation to have long-lasting effects so far has relied on computational studies (Tass & Majtanik, 2006; Popovych & Tass, 2012; Ebert et al., 2014; Popovych et al., 2015; Manos et al., 2018) or on analytical insights under the simplifying assumption that spiking is triggered only by stimulation pulses (Kromer & Tass, 2020). Other analytical approaches based on mean-field models are focused on short-term changes due to stimulation and do not consider plastic changes (Duchet et al., 2020; Weerasinghe et al., 2019, 2021). Our framework offers an alternative, where exact evolution equations for the average coupling within and between neural populations could inform the development of new therapies aimed at maximizing the long-term effects of brain stimulation. Although microscopic connectivity is unknown, stimulation can be designed to change average connectivity and impact synchrony as beneficial.

## Appendix A: Simulation Methods Used to Compare STDP to PDDP

Using simulations, we compare networks of Kuramoto oscillators with adaptivity given by symmetric or causal STDP, to the same networks with adaptivity given by the corresponding single- or multiharmonic PDDP and ebPDDP rules.

### A.1 Symmetric Learning Rules

We simulate networks of $N=60$ Kuramoto oscillators evolving according to equation 2.1, where the natural frequencies $\omega k$ are sampled from a normal distribution of mean $\Omega =10\pi $ (5 Hz) unless otherwise stated and standard deviation $\Delta $. A normal distribution was chosen over a Lorentzian distribution to avoid extreme natural frequencies, which lead to increased variability in simulation repeats. Initial conditions are also sampled from normal distributions. Phases $\theta k(t=0s)$ are sampled from the circular equivalent of $N0,\pi 2/9$ and couplings $\kappa kl(t=0s)$ from $N5,\sigma \kappa 2$. The evolution of the coupling weights $\kappa kl$ is simulated according to equation 1.2 (symmetric STDP), together with the weight decay $d\kappa kldt=-\epsilon \kappa kl$, equation 1.3, with $\varphi =0$ (symmetric PDDP), or equation 2.2 (symmetric ebPDDP). For ebPDDP, we use $\delta (t-tkq)\u22481Ikq(t)/\Delta t$ where $1Ikq$ is the indicator function of $Ikq=[tkq,tkq+\Delta t)$ and $\Delta t$ is the simulation time step. The network is simulated for 150 s, using the Euler method with $\Delta t=1$ ms unless otherwise stated. We use $\epsilon =0.5$, and the values of other parameters are detailed in Figures 3 and 11.

### A.2 Causal Learning Rules

Error metrics (described below) are computed between multiharmonic PDDP or ebPDDP and causal STDP for selected regions of parameter space. To limit computational cost, we constrain causal STDP parameters based on data and biological motivations. We use $\tau +=16.8$ ms and $\tau -=33.7$ ms, which were obtained by fitting equation 1.1 to experimental data (data published in Bi & Poo, 1998, and fit in Bi & Poo, 2001). Several experimental studies have reported the LTD time constant $\tau -$ to be larger than the LTP time constant $\tau +$, for example, in the rat hippocampus (Bi & Poo, 1998), somatosensory cortex (Feldman, 2000), and visual cortex (Froemke & Dan, 2002). We take $A+=0.2$ and $A-/A+=\beta $. Since the time window for LTD is twice as long as the time window for LTP, we choose $\beta =0.5$ to study a balanced situation where LTP dominates for shorter spike-timing differences and LTD dominates for longer spike-timing differences, and $\beta =1$ to study a situation where LTD always dominates.

We simulate networks of Kuramoto oscillators as for symmetric STDP (see section A.2) with the following differences. Initial couplings $\kappa kl(t=0s)$ are sampled from $N12,\sigma \kappa 2$. The evolution of the coupling weights $\kappa kl$ is simulated according to equation 1.1 (causal STDP), equation 2.3 (causal PDDP), or equation 2.4 (causal ebPDDP). For PDDP and ebPDDP, the Fourier expansion of $F$ is truncated after $Nf$ coefficients. No weight decay term is included.

## Appendix B: Methodological Details Pertaining to the Optimization of the Two-Population Mean-Field Approximation

Synthetic data are generated by stimulating $N=100$ Kuramoto oscillators evolving according to equation 2.1, with adaptivity given by symmetric PDDP (equation 1.3 with $\varphi =0$, $\lambda =25$, and $\epsilon =0.5$). Oscillator natural frequencies $\omega k$ are sampled from a Lorentzian distribution of center $\Omega =10\pi $ (5 Hz) and width $\Delta =0.6\pi $. Initial phases are sampled from a Von Mises distribution of standard deviation $\pi /4$. Initial couplings $\kappa kl(t=0s)$ are sampled from $N\kappa ^0,0.52$. The network is simulated for 20 s, using the Euler method with $\Delta t=0.1$ ms.

We create a training set and a test set based on different initial conditions and in particular various $\kappa ^0$. The training set corresponds to $\kappa ^0={2$, 2, 2, 3, 4, 5, 5, 10, 10, 15, 15, 15, 20$}$ and the test set to $\kappa ^0={3.5$, 3.5, 3.5, 3.5, 7, 7, 7, 7, 12, 12, 12, 12, 17, 17, 17, 17, 22, 22, 22, 22$}$. For each trajectory in the training and test sets, natural frequencies, initial phases, and initial couplings are sampled from their respective distributions. Repeated $\kappa ^0$ values therefore correspond to different systems with different initial synchrony.

For optimization speed and accuracy, the two-population mean-field approximation is simulated using the variable order solver ode113 in Matlab (variable-step, variable-order Adams-Bashforth-Moulton solver of orders 1 to 13).

## Appendix C: Supplementary Figures and Tables

Parameter . | $\epsilon 1$ . | $\epsilon 2$ . | $\lambda 1$ . | $\lambda 2$ . |
---|---|---|---|---|

Value | 0.0263 | 3.0262 | 5.2985 | 5.4759 |

Parameter . | $\epsilon 1$ . | $\epsilon 2$ . | $\lambda 1$ . | $\lambda 2$ . |
---|---|---|---|---|

Value | 0.0263 | 3.0262 | 5.2985 | 5.4759 |

Notes: Parameters correspond to equation 5.1. Full network parameters used to generate synthetic data for the optimization are given in appendix B.

## Acknowledgments

We are grateful to Rafal Bogacz and Alessandro Torcini for helpful discussions. We acknowledge the use of the University of Oxford Advanced Research Computing facility in carrying out this work http://dx.doi.org/10.5281/zenodo.22558

## Funding Information

B.D. was supported by Medical Research Council grant MC_UU_00003/1. C.B. acknowledges support from the Engineering and Physical Sciences Research Council through grant EP/T013613/1.

## Data Availability Statement

No new data were created or analyzed in this study. We are sharing code at https://github.com/benoit-du/STDPvsPDDP to compare the evolution of a Kuramoto network with STDP to the evolution of a Kuramoto network with PDDP and ebPDDP.

## Notes

^{1}

In this work, we call plasticity rules *causal* when temporal precedence is enforced (as in “classical” STDP rules; Bi & Poo, 1998), and *symmetric* when inverting spike timings and phase differences leads to the same type of adaptation. We stay clear of the ambiguous term *Hebbian* STDP, which is also sometimes used to refer to causal STDP in the experimental literature (Cassenaer & Laurent, 2007; Sgritta et al., 2017).

^{2}

The division by 2 ensures that this rule scales similarly to the STDP rule and the PDDP rule with continuous updates.

^{3}

Note that equilibrium points of the effective dynamics correspond to periodic solutions of the full system.