STDP and the distribution of preferred phases in the whisker system

Rats and mice use their whiskers to probe the environment. By rhythmically swiping their whiskers back and forth they can detect the existence of an object, locate it, and identify its texture. Localization can be accomplished by inferring the whisker’s position. Rhythmic neurons that track the phase of the whisking cycle encode information about the azimuthal location of the whisker. These neurons are characterized by preferred phases of firing that are narrowly distributed. Consequently, pooling the rhythmic signal from several upstream neurons is expected to result in a much narrower distribution of preferred phases in the downstream population, which however has not been observed empirically. Here, we show how spike timing dependent plasticity (STDP) can provide a solution to this conundrum. We investigated the effect of STDP on the utility of a neural population to transmit rhythmic information downstream using the framework of a modeling study. We found that under a wide range of parameters, STDP facilitated the transfer of rhythmic information despite the fact that all the synaptic weights remained dynamic. As a result, the preferred phase of the downstream neuron was not fixed, but rather drifted in time at a drift velocity that depended on the preferred phase, thus inducing a distribution of preferred phases. We further analyzed how the STDP rule governs the distribution of preferred phases in the downstream population. This link between the STDP rule and the distribution of preferred phases constitutes a natural test for our theory.

During whisking, the animal moves its vibrissae back and forth in a rhythmic manner, Fig 1a-1c. Neurons that track the azimuthal position of the whisker by firing in a preferential manner to the phase of the whisking cycle are termed whisking neurons. Whisking neurons in the ventral posteromedial nucleus (VPM) of the thalamus as well as inhibitory whisking neurons in layer 4 of the barrel cortex are characterized by a preferred phase at which they fire with the highest rate during the whisking cycle [1,6,[12][13][14][15][16] , Fig 1d and 1e. The distribution of preferred phases is non-uniform and can be approximated by the circular normal (Von-Mises) distribution Prð�Þ ¼ e k cosð�À cÞ 2pI 0 ðkÞ ð1Þ where ψ is the mean phase and I 0 (κ) is the modified Bessel function of order 0, Fig 2a. The parameter κ quantifies the width of the distribution; κ = 0 yields a uniform (flat) distribution, whereas in the limit of κ ! 1 the distribution converges to a delta function. Typical values for κ in the thalamus and for layer 4 inhibitory (L4I) whisking neurons are κ VPM � κ L4I � 1 where ψ VPM � 5π/6 [rad] and κ L4I � 0.5 [rad] [2,13]. Throughout this work we ignore the possible contribution of recurrent connections to the shaping of phase selectivity in layer 4 (but see Discussion). Assuming the rhythmic input to L4I neurons originates solely from the VPM, the distribution of preferred phases of L4I neurons is determined by the distribution in the thalamus and the profile of thalamo-cortical synaptic weights. However, synaptic weights have been reported to be highly volatile [17][18][19][20]. If synaptic plasticity is activity independent, then one would expect the synaptic weight profile to become unstructured. In the naive pooling model, which lacks activity-dependent plasticity, we consider two types of unstructured connectivity; namely, uniform pooling and random pooling. Fig 2b shows the distribution of preferred phases in a naive model that pools the responses of N = 100 VPM neurons. The source of variability in the preferred phases of the downstream neuron is the random choice of N VPM neurons, each characterized by a preferred phase that is distributed in an i.i.d. manner following Eq (1). As a result, one would expect that the number, N, of pooled thalamic neurons would govern the width of the distribution. Fig 2c shows the expected distribution width, κ L4I , as a function of N for uniform (squares) and random (circles) pooling. As can be seen from the figure, even a random pooling of N = 10 results in a distribution of preferred phases that is considerably narrower than empirically observed. This result is particularly surprising since the number of thalamic neurons synapsing onto a single L4 neuron was estimated to be on the order of 100, see e.g. [16,21]. Moreover, in the naive pooling model, the mean preferred phase of the downstream population is given by the mean phase of the upstream population. The delay in the response of downstream neurons to their input should also be added. It is estimated to be around 1-10ms; hence, for whisking at 10Hz translates into ψ VPM − ψ L4I = 0.06-0. 6 [rad]. This result also runs counter empirical findings reporting ψ VPM − ψ L4I � 2.2 [rad]. Thus, the synaptic weight profile in a model that lacks activity-dependent plasticity is expected to converge to a random pooling (due to constant synaptic remodelling), which is inconsistent with the observed distribution of preferred phases.
Recently, the effects of rhythmic activity on the spike timing dependent plasticity (STDP) dynamics of feed-forward synaptic connections have been examined [22,23]. It was shown that in this case the synaptic weights remain dynamic. As a result, the phase of the downstream neuron is arbitrary and drifts in time; thus, effectively, inducing a distribution of preferred phases in the downstream population. However, if the phases of the downstream population are arbitrary and drift in time, how can information about the whisking phase be transmitted?
Here we investigated the hypothesis that the distribution of phases in the downstream layer is governed by the interplay of the distribution in the upstream layer and the STDP rule. The The preferred phases in the upstream population were drawn in an i.i.d. manner following Eq (1), with κ UpStream = 1 and ψ UpStream = 0. The dashed red line depicts the mean phase. (c) Distribution width in the naive uniform/random pooling model. The width of the distribution of preferred phases, κ L4I , in the downstream layer (L4I) is shown as a function of the number of pooled VPM neurons, N, for the uniform and random pooling in squares and circles, respectively. The width of the distribution, κ L4I , was estimated from 10000 repetitions of drawing N preferred phases of the upstream population with κ = 1 (blue) and κ = 2 (red).
https://doi.org/10.1371/journal.pcbi.1009353.g002 remainder of this article is organized as follows. First, we define the network architecture and the STDP learning rule. We then derive a mean-field approximation for the STDP dynamics in the limit of a slow learning rate for a threshold-linear Poisson downstream neuron model. Next, we show that STDP dynamics can generate non-trivial distributions of preferred phases in the downstream population and analyze how the parameters characterizing the STDP govern this distribution. Finally, we summarize the results, discuss the limitations of this study and suggest essential empirical predictions deriving from our theory.

The upstream thalamic population
We model a system of N thalamic excitatory whisking neurons, synapsing in a feed-forward manner onto a single inhibitory L4 barrel cortical neuron. Unless stated otherwise, following [16], in our numerical simulations we used N = 150. The spiking activity of the thalamic neurons is modelled by independent inhomogeneous Poisson processes with an instantaneous firing rate that follows the whisking cycle: . .N}, is the spike train of the kth thalamic neuron, with ft k;i g 1 i¼1 denoting its spike times. The parameter D is the mean firing rate during whisking (averaged over one cycle), γ is the modulation depth, ν is the angular frequency of the whisking, and ϕ k is the preferred phase of firing of the kth thalamic neuron. We further assume that the preferred phases in the thalamic population, f� k g N k¼1 , are distributed i.i.d. according to the von-Mises distribution, Eq (1).

The downstream layer 4 inhibitory neuron model
To facilitate the analysis, we model the response of the downstream layer 4 inhibitory (L4I) neuron to its thalamic inputs by the linear Poisson model, which has been used frequently in the past [22,[24][25][26][27][28][29][30][31]. Given the thalamic responses, the firing of the L4I neuron follows inhomogeneous Poisson process statistics with an instantaneous firing rate where d > 0 represents the characteristic delay, and w k is the synaptic weight of the kth thalamic neuron. Due to the rhythmic nature of the thalamic inputs, Eq (2), and the linearity of the L4I neuron, Eq (3), the L4I neuron will exhibit rhythmic activity: with a mean, D L4I , a modulation depth, γ L4I , and a preferred phase ψ L4I , that depend on global order parameters that characterize the thalamocortical synaptic weights population. For large N these order parameters are given by: andw where � w is the mean synaptic weight andwe ic is its first Fourier component. The phase ψ is determined by the condition thatw is real and non-negative. Consequently, the L4I neurons in our model respond to whisking with a mean D L4I ¼ D � w, a modulation depth of g L4I ¼ gw= � w, and a preferred phase ψ L4I = ψ + νd.

The STDP rule
We model the modification of the synaptic weight, Δw, following either a pre-or post-synaptic spike as the sum of two processes: potentiation (+) and depression (-) [31][32][33], as The parameter λ denotes the learning rate. We further assume a separation of variables and write each term (potentiation and depression) as the product of the function of the synaptic weight, f ± (w), and the temporal kernel of the STDP rule, K ± (Δt), where Δt = t post − t pre is the time difference between pre-and post-synaptic spike times. Following Gütig et al. [33] the weight dependence functions, f ± (w), were chosen to be: where α > 0 is the relative strength of depression and μ 2 [0, 1] controls the non-linearity of the learning rule. The temporal kernels of the STDP rule are normalized: i.e., R K ± (Δt)dΔt = 1. Here, for simplicity, we assume that all pairs of pre and post spike times contribute additively to the learning process via Eq (7).

STDP dynamics in the limit of slow learning
In the limit of a slow learning rate, λ ! 0, we obtain deterministic dynamics for the mean synaptic weights (see [32] for a detailed derivation) where Γ j, L4I (Δ) is the cross-correlation function between the jth thalamic pre-synaptic neuron and the L4I post-synaptic neuron; see detailed derivation in Temporal correlations.

STDP dynamics of thalamocortical connectivity
We simulated the STDP dynamics of 150 upstream thalamic neurons synapsing onto a single L4I neuron in the barrel cortex, see Details of the numerical simulations & statistical analysis. Fig 3a shows the temporal evolution of the synaptic weights (color-coded by their preferred phases). As can be seen from the figure, the synaptic weights do not relax to a fixed point; instead there is a continuous remodelling of the entire synaptic population. Examining the order parameters, Fig 3b and 3c, reveals that the STDP dynamics converge to a limit cycle.
The continuous remodelling of the synaptic weights causes the preferred phase of the downstream neuron, κ L4I (see Eq (4)) to drift in time, Fig 3b. As can be seen from the figure, The distribution of preferred phases in the thalamic population that served as input to the downstream L4I neuron is presented as a polar histogram. The distribution followed Eq (1) with κ VPM = 1 and ψ VPM = 5π/6 rad. (e) The temporal distribution of preferred phases of the downstream L4I neuron is shown in a polar plot. Fitting the von-Mises distribution yielded κ L4I = 1.1 and κ L4I = 0.8 rad. In this simulation we used the following parameters: N = 150, � n ¼ n=ð2pÞ ¼ 7hz, γ VPM = 1, D VPM = 10hz, and d = 3ms. The temporally asymmetric STDP rule, Eq (9), was used with τ − = 50ms, τ + = 22ms, μ = 0.01, α = 1.1, and λ = 0.01.
https://doi.org/10.1371/journal.pcbi.1009353.g003 the drift velocity is not constant. Consequently, the downstream (L4I) neuron 'spends' more time in certain phases than others. Thus, the STDP dynamics induce a distribution over time for the preferred phases of the downstream neuron. One can estimate the distribution of the preferred phases of L4I neurons by tracking the phase of a single neuron over time. Alternatively, since our model is purely feed-forward, the preferred phases of different L4I neurons are independent; hence, this distribution can also be estimated by sampling the preferred phases of different L4I neurons at the same time. Fig 3d and 3e show the distribution of preferred phases of thalamic and L4I whisking neurons, respectively. Thus, STDP induces a distribution of preferred phases in the L4I population, which is linked to the temporal distribution of single L4I neurons. Fig 4a shows the (non-normalized) distribution of preferred phases induced by STDP as a function of the number of pooled thalamic neurons, N. For large N, the STDP dynamics converge to the continuum limit of Eq (11), and the distribution converges to a limit that is independent of N. This contrasts with the naive pooling model lacking plasticity (cf. Fig 4b). In that model, the distribution of preferred phases in the cortical population results from a random process of pooling phases from the thalamic population. We shall refer to this type of variability as quenched disorder, since this randomness is frozen and does not fluctuate over time. This process is characterized by a distribution of preferred phases with vanishing width, κ L4I ! 1, in the large N limit,   (1). (c) The distribution due to both quenched disorder and STDP dynamics is shown. Unless stated otherwise, the parameters used in these graphs are as follows: � n ¼ n=ð2pÞ ¼ 7Hz, κ VPM = 1, γ = 0.9, D = 10hz, d = 3ms, ϕ 0 = 5π/6. For the STDP we used the temporally exponential asymmetric kernel, Eq (9), with τ − = 50ms, τ + = 22ms, μ = 0.01, α = 1.1, and λ = 0.01. https://doi.org/10.1371/journal.pcbi.1009353.g004 As N increases, the distribution widens and is centered around a preferred phase that is determined by the STDP. Thus, activity dependent plasticity helps to shape the distribution of preferred phases in the downstream population. Below, we examine how these different parameters affect the ways in which STDP shapes the distribution; hence, we will not average over the quenched disorder.
Parameters characterizing the upstream input. Fig 5a depicts the distribution of preferred phases as a function of the distribution width of their thalamic input, κ VPM . For a uniform input distribution, κ VPM = 0, the downstream distribution is also uniform, κ L4I = 0, see [22]. As the distribution in the VPM becomes narrower, so does the distribution in the L4I population. If the distribution of the thalamic population is narrower than a certain critical value, STDP will converge to a fixed point and κ L4I will diverge. Typically, we find that the width of the cortical distribution, κ L4I , is similar to or larger than that of the upstream distribution, κ VPM , Fig 5b. This sharpening is obtained via STDP by selectively amplifying certain phases while attenuating others. Consequently the rhythmic component, in terms of the modulation depth, γ L4I , is also typically amplified relative to the uniform pooling model, Fig 5c. The effect of the whisking frequency is shown in Fig 6a. For moderate rhythms that are on a similar timescale to that of the STDP rule, STDP dynamics can generate a wide distribution of preferred phases. However, in the high frequency limit, ν ! 1, the synaptic weights converge to a uniform solution with w(ϕ) = (1 + α 1/μ ) −1 , 8ϕ (see [22]). In this limit, due to the uniform pooling, there is no selective amplification of phases and the whisking signal is transmitted downstream due to the selectivity of the thalamic population, κ VPM > 0. Consequently, at high frequencies, the distribution of preferred phases will be extremely narrow, Oð1= ffi ffi ffi ffi N p Þ because of the quenched disorder, and the rhythmic signal will be attenuated, γ L4I = γ VPM × I 1 (κ)/I 0 (κ) (where I j (κ) is the modified Bessel function of order j) , Fig 6b. The rate of convergence to the high frequency limit is governed by the smoothness of the STDP rule. Discontinuity in the STPD rule, such as in our choice of a temporally asymmetric rule, will induce algebraic convergence in ν to the high frequency limit, whereas a smooth rule, such as our choice of a temporally symmetric rule, will manifest in exponential convergence; compare Fig 6 with Fig 7, see also [22]. In our choice of parameters, τ + � 20ms and τ − � 50ms, STDP dynamics induced a The effects of synaptic weight dependence, μ. Previous studies have shown that increasing μ weakens the positive feedback of the STDP dynamics, which generates multi-stability, and stabilizes more homogeneous solutions [22,23,33,46]. This transition is illustrated in Fig  8: at low values of μ, the STDP dynamics converge to a limit cycle in which both the synaptic weights and the phase of the the L4I neuron cover their entire dynamic range , Fig 8a and 8b. As μ is increased, the synaptic weights become restricted in the limit cycle and no longer span their entire dynamic range , Fig 8c and 8d. A further increase of μ also restricts the phase of the L4I neuron along the limit cycle, Fig 8e and 8f. Finally, when μ is sufficiently large, STDP dynamics converge to a fixed point , Fig 8g and 8h. This fixed point selectively amplifies certain phases, yielding a higher value of γ than in the uniform solution. These results are summarized in Fig 9 that shows the (non-normalized) distribution of preferred phases and the relative strength of the rhythmic component, γ L4I , as a function of μ. Note that except for a small value  range of μ, around which the distribution of preferred phases is bi-modal (see, e.g . Fig 8e and  8f and μ � 0.07 in Fig 9), STDP yields higher γ L4I values than in the uniform pooling model.
The relative strength of depression, α. As one might expect, decreasing α beyond a certain value will result in a potentiation dominated STDP dynamics, saturating all synapses close to their maximal value, thus, approaching the uniform solution, which is characterized by a narrow preferred phase distribution centered around the mean preferred phase of the VPM neurons shifted by the delay, and low values of γ L4I , Fig 10. Increasing α strengthens the competitive winner-take-all like nature of the STDP dynamics. Initially, this competition will generate a fixed point with non-uniform synaptic weights, thus increasing γ L4I . A further increase of α results in a limit cycle solution to the STDP dynamics that widens the distribution of the preferred phases. Increasing α beyond a certain critical value will result in depression dominated STDP dynamics, thus driving the synaptic weights to zero, Fig 10. Parameters characterizing the temporal structure of STDP. Fig 11 shows the effect of the temporal structure of the STDP on the distribution of preferred phases in the downstream population. As can be seen from the figure, varying the characteristic timescales of potentiation and depression, τ + and τ − , induces quantitative effects on the distribution of the preferred phases. However, qualitative effects result from changing the nature of the STDP rule. Comparing the temporally asymmetric rule , Fig 11a and 11b, to the temporally symmetric rule ,  Fig 11c and 11d, reveals a dramatic difference in the mean preferred phase of L4I neurons, dashed black lines.

Discussion
We examined the possible contribution of STDP to the diversity of preferred phases of whisking neurons. Whisking neurons can be found along different stations throughout the somatosensory information processing pathway [1, 4-6, 9, 10, 13, 16, 47]. Here we focused on L4I whisking neurons that receive their whisking input mainly via excitatory inputs from the VPM, and suggested that the non-trivial distribution of preferred phases of L4I neurons results from a continuous process of STDP.  STDP has been reported in thalamocortical connections in the barrel system. However, STDP has only been observed during the early stages of development [43][44][45]. Is it possible that the contribution of STDP to shaping thalamocortical connectivity is restricted to the early developmental stages? This is an empirical question that can only be answered experimentally. Nevertheless, several comments can be made.
First, Inglebert and colleagues recently showed that pairing single pre-and post-synaptic spikes, under normal physiological conditions, does not induce plastic changes [48]. On the other hand, activity dependent plasticity was observed when stronger activity was induced. Thus, whisking activity, which is a strong collective and synchronized activity, may induce STDP of thalamocortical connectivity.
Second, in light of the considerable volatility of synaptic connections observed in the central nervous system [17][18][19][20], it is hard to imagine that thalamocortical connectivity remains unchanged throughout the lifetime of an animal. Third, if activity independent plasticity alone underlies synaptic volatility, thalamocortical synaptic weights would be expected to be random. In this case, thalamic whisking input to layer 4 should be characterized by an extremely narrow distribution centered around the delayed mean thalamic preferred phase, cf. Fig 4b. Consequently, neglecting the possible contributions due to recurrent interactions with layer 4, a non-trivial distribution of L4I phases, with κ L4I * 1, can be obtained either via STDP or by pooling the whisking signal from an extremely small VPM population of N < 10 neurons. However, in the latter scenario (without STDP) the mean preferred phase of L4I neurons is expected to be determined by the (delayed) mean phase of VPM neurons (see e.g . Fig 4b), thus raising serious doubts as to the viability of the latter solution.
An alternative hypothesis to our proposed mechanism is that the recurrent dynamics within layer 4 generate the empirically observed wide distribution of preferred phases centered on κ L4I � ψ VPM − 2.2 [rad]. As layer 4 excitatory neurons have been reported to exhibit very weak rhythmic activity [13,16], such a mechanism would mainly be based on recurrent interaction within the L4I population. The putative ways in which inhibitory interactions widen the distribution remain unclear. Nevertheless, this hypothesis should be explored further.
Our hypothesis views STDP as a continuous process. As a result of STDP, the preferred phase of a L4I neuron, in our model, is a stochastic process with a non-uniform drift, Fig 3b, which, in turn, induces a distribution over time for the phase , Fig 3e. As our model lacks recurrent interactions, phases of different L4I neurons will be i.i.d. random processes, albeit with different random initial conditions. Consequently, simultaneous sampling of the distribution of preferred phases of a population of L4I neurons will yield the same distribution as obtained by sampling the preferred phase of a single neuron over time. Thus, functionality, in terms of transmission of the whisking signal and retaining the stationary distribution of the preferred phases in the downstream population is obtained as a result of continuous remodelling of the entire population of synaptic weights.
Several key features of our hypothesis provide clear empirical predictions. Our theory links two empirical observations. One is the distribution of preferred phases in up-and downstream populations. The other is the STDP rule. The literature can provide a good assessment of the distribution of preferred phases of thalamic and of cortical whisking neurons. Further effort is required to estimate the STDP rule of the thalamocortical connections to L4I neurons. Thus, we can reject our hypothesis if the predicted distribution, based on our theory and the estimated learning rule, is not consistent with the empirically observed distribution.
Our theory also posits that the distribution of preferred phases results from their stochastic dynamics. Our theory also incorporates a direct relationship between the drift velocity and the distribution of preferred phases. Essentially, the system will spend more time around phases characterized by slower velocity; hence, a higher probability of preferred phases should correlate with slower drift velocity. To the best of our knowledge, the dynamics of the preferred phases of whisking neurons has not been studied. Monitoring the preferred phases of single L4I neurons over long periods of time (days and weeks) could serve as additional tests of our theory. Specifically, it would address the question of whether the preferred phases drift in time, and if so, whether the drift velocity and the distribution of preferred phases are consistent with our theory.
Further empirical exploration can test the degree in which our approximation of a purely feed-forward architecture (i.e., neglecting recurrent interactions) holds. Simultaneous recordings of a large population of fast spiking layer 4 whisking neurons can be harnessed to compare the distribution of preferred phases over space and time, which are expected to be identical in the absence of recurrent interactions. In addition, in the visual system, Lien and Scanziani combined intracellular recordings with optogenetic methods to assess the contribution of thalamic inputs to the selectivity of cortical neurons [49]. In the context of the current work, we explored the thalamic inputs to whisking neurons in the barrel cortex. Thus, applying ideas and approaches from the work of Lien and Scanziani could (i) provide a measurement of the thalamic input, and (ii) determine the relative contribution of feed-forward and recurrent connections to the selectivity to the phase of whisking in L4I neurons.
More broadly, our findings raise the question of whether the constant remodelling of synaptic efficacies is an artefact of the 'biological hardware' the brain must overcome, or whether it reflects an underlying principle of the central nervous system. Though we cannot provide an answer to this question at this point in time, we may speculate on the possible advantages of a dynamic solution. One interesting hypothesis to consider is that the constant remodelling makes the neuronal network more flexible and consequently enables a faster response to a changing environment.
In the current study, we made several simplifying assumptions. The spiking activity of the thalamic population was modeled as a rhythmic activity with a well-defined frequency. However, whisking activity spans a frequency band of several Hertz. Moreover, the thalamic input relays additional signals, such as touch and texture. These signals modify the cross-correlation structure and add 'noise' to the dynamics of the preferred phase of the downstream neuron. As a result, the distribution of preferred phases in the downstream population is likely to be wider. In addition, our analysis used a purely feed-forward architecture ignoring recurrent connections in layer 4, which may also affect the preferred phases in layer 4. A quantitative application of our theory to the whisker system should consider all these effects. Nevertheless, the theory presented provides the basic foundation to address these effects.

Temporal correlations
The cross-correlation between pre-synaptic neurons j and k at time difference Δt is given by: In the linear Poisson model, Eq (3), the cross-correlation between a pre-synaptic neuron and the post-synaptic neuron can be written as a linear combination of the cross-correlations in the upstream population; hence, the cross-correlation between the jth VPM neuron and the post-synaptic neuron is Where � w andwe ic are order parameters characterizing the synaptic weights profile, as defined in Eqs (5) and (6).

The mean field Fokker-Planck dynamics
For large N we obtain the continuum limit from Eq (11) and � K � ,K � e iO Z � are the Fourier transforms of the STDP kernels Note that in our specific choice of kernels, � K � ¼ 1, by construction.

Fixed points of the mean field dynamics
The fixed point solution of Eq (15) is given by where X � �w � w g 2 2K � cosð� À nd À O � À cÞ: ð20Þ Note that, from Eqs (19) and (20), the fixed point solution, w(ϕ) � , will depend on ϕ, for κ VPM > 0. As μ grows to 1 the fixed point solution will become more uniform, see [22].

Details of the numerical simulations & statistical analysis
Scripts of the numerical simulations were written in Matlab. The numerical results presented in this paper were obtained by solving Eq (15) with the Euler method with a time step Δt = 0.1 and λ = 0.01. The cftool with the Nonlinear Least Squares method and a confidence level of 95% for all parameters was used for fitting the data presented in Figs 3, 5b, 6b and 7b.

Modeling pre-synaptic phase distributions
Unless stated otherwise, STDP dynamics in the mean field limit was simulated without quenched disorder. To this end, the preferred phase, ϕ k , of the kth neuron in a population of N pre-synaptic VPM neurons was set by the condition R φ k À p PrðφÞd� ¼ k=N. In Fig 4b and 4c we used the accept-reject method [50,51] to sample the phases.