Phase dependence of response curves to stimulation and their relationship: from a Wilson-Cowan model to essential tremor patient data

Essential tremor manifests predominantly as a tremor of the upper limbs. One therapy option is high-frequency deep brain stimulation, which continuously delivers electrical stimulation to the ventral intermediate nucleus of the thalamus at about 130 Hz. Investigators have been looking at stimulating less, chiefly to reduce side effects. One strategy, phase-locked deep brain stimulation, consists of stimulating according to the phase of the tremor, once per period. In this study, we aim to reproduce the phase dependent effects of stimulation seen in patient data with a biologically inspired Wilson-Cowan model. To this end, we first analyse patient data, and conclude that half of the datasets have response curves that are better described by sinusoidal curves than by straight lines, while an effect of phase cannot be consistently identified in the remaining half. Using the Hilbert phase we derive analytical expressions for phase and amplitude responses to phase-dependent stimulation and study their relationship in the linearisation of a stable focus model, a simplification of the Wilson-Cowan model in the stable focus regime. Analytical results provide a good approximation for response curves observed in patients with consistent significance. Additionally, we fitted the full non-linear Wilson-Cowan model to these patients, and we show that the model can fit in each case to the dynamics of patient tremor as well as the phase response curve, and the best fits are found to be stable foci for each patients (tied best fit in one instance). The model provides satisfactory prediction of how patient tremor will react to phase-locked stimulation by predicting patient amplitude response curves although they were not explicitly fitted. This can be partially explained by the relationship between the response curves in the model being compatible with what is found in the data. We also note that the non-linear Wilson-Cowan model is able to describe response to stimulation more precisely than the linearisation.


Introduction
Essential tremor (ET) is the most common movement disorder, affecting 0.9% of the population [1]. It predominantly manifests as a tremor of the upper limbs, and can severely affect daily-life. When medications are ineffective or not tolerated, thalamic deep brain stimulation (DBS) is a well-established therapy option. Clinically available DBS continuously delivers high-frequency (≈ 130 Hz) electrical stimulation to deep structures within the brain via an electrode connected to a pulse generator implanted in the chest. There is no agreement in the research community on the mechanisms of action of high-frequency DBS [2], but it is believed there is room for improvement in terms of efficacy, decrease in power usage, and most importantly reduction of side effects [3]. Reported side effects of high-frequency thalamic DBS include speech impairment (dysarthria), gait disorders, depression, and abnormal dermal sensations (paresthesia). Improving high-frequency DBS generally means stimulating less by closing the loop on a signal related to motor symptoms, while maintaining clinical efficacy. One example of closed-loop DBS is adaptive DBS, whereby stimulation is triggered in Parkinson's disease (PD) patients when pathological neural oscillation amplitude in the beta band is higher than a threshold. Compared to high-frequency DBS, it has been shown to improve motor performance, and reduce speech side-effects in humans [4,5,6]. Another example is phase dependent stimulation, which has been investigated in a computational model of PD [7], and in PD patients [8,9]. Phase-locked DBS has recently been studied as a new therapy for ET [10]. Hand tremor is recorded, and the reduction in stimulation comes from stimulating only once per period of the tremor, rather than continuously, according to the phase of tremor. In some patients, the strategy only requires half the energy delivered by high-frequency DBS for the same effect. Optimising phase-locked DBS requires a detailed understanding of the phase dependence of the response across patients, but so far data collection from phase-locked stimulation experiments has been restricted to small datasets because patients fatigue quickly. While direct analysis of the data has proven insightful [10], modelling phase-locked stimulation would allow predictions to be made from analytic and computational studies regarding the phase dependence of the response to stimulation, and would open the door to supplement scarcely available patient data with synthetic data. The ability to easily generate large amounts of synthetic data could come in handy to help devise and test control algorithms, or when trying to predict an effect that, because of noise in recordings, can only be deciphered when a large number of trials is available.
Tremulous hand movements are believed to be closely related to thalamic activity [11,12], and it is believed that ET originates in the cerebellar-thalamic-cortical pathway [13]. However detailed knowledge of how ET comes about is missing, which makes simple, canonical models natural candidates to study ET. Recently, phaselocked DBS was studied using a Kuramoto oscillator model which conveniently captures phases, but does not model interacting neural populations with distinct properties [14]. In the present work, we focus on a neural mass model, the Wilson-Cowan (WC) model, whose architecture can be mapped onto the neural populations thought to be involved in the generation of ET, and allows for strong coupling between the populations. Additionally, stimulation can target the most common stimulation site for ET, the ventral intermediate nucleus (VIM). The model describes the firing rates of an excitatory and an inhibitory population, and only has a few parameters, which makes it less prone to overfitting and significantly easier to constrain than more detailed models. The WC model has been shown to be adept at describing beta oscillations in PD [15,16]. Moreover, the work presented in [17] provides evidence that the effects of high-frequency DBS for ET in a WC model are similar to the description given by conductance-based models. While the WC model has been used to design closed-loop strategies for PD [18,19], whether a firing-rate model such as the WC model can model the effects of phase-locked DBS has not been approached in the literature. Based on strong assumptions, Polina et al. reduced a WC model to a one-dimensional ordinary differential equation and looked at periodic forcing, but not in the context of DBS, and without attending to the dependence on the phase of stimulation [20]. The present work will focus on reproducing the phase-dependent effects of phase-locked DBS measured in human data with a WC model.
Stimulation changes the phase and the amplitude of tremor and the dependence of these changes on the phase of stimulation can be quantified by the phase response curve (PRC, in this study change in tremor phase as a function of tremor phase) and the amplitude response curve (ARC, in this study change in tremor amplitude as a function of tremor phase). The ARC directly measures the change in tremor, hence the change in patient handicap, but both the ARC and PRC are important to understand the effects of phase-locked DBS and potentially optimise the stimulation pattern. Theoretically, PRCs and ARCs have been defined differently, mostly in the context of limit cycle models concerned with asymptotic response to infinitesimal perturbations, see for example [21,22,23,24,25]. In patients, DBS stimulation is not infinitesimal, and tremor data is very variable so stimulation happens in transient states. Therefore rather than considering an asymptotic description of the changes in phase and amplitude, we will be focusing on the experimental response curve measurement methodology applied to blocks of stimulation in [10], which we will hereafter refer to as the "block method". It provides a finite time response and relies on the changes in the Hilbert phase and amplitude following blocks of phase-locked stimulation (more details in 2.1). The only exception to this will be in analytical derivations (section 4), where a first order measurement of the response curves (i.e. measurement at the end of the stimulation period) will be used for tractability, as a simplified first approach to the model. For coherence with the experimental response curve measurement methodology, the notion of phase and amplitude used throughout will be the Hilbert phase and amplitude. It should also be noted that we are considering population response curves and not single neuron response curves. The best performing WC models in reproducing patient data are found in this work to be stable foci (tied best fit in one instance), where tremor dynamics are being reproduced by adding noise to the system, so we restrict our analytical considerations to stable foci.
The main contributions of this work are the following. We first analyse patient response curves, identify a subset of datasets passing appropriate statistical tests, and characterise the relationship between PRC and ARC in these patients (section 2). Following the introduction of our biologically motivated WC model (section 3), we derive approximate analytical expressions that delineate the response to stimulation of a 2D dynamical system described by a linearised focus, with the goal of better understanding the constraints built in the model (section 4). The derived response curves are close to sinusoidal, and a relationship between them is found, revealing similarities in shape and phase shift with patients who have significant PRCs and ARCs. We then show that for these patients, the WC model can be fitted to the data and can reproduce the dependence of the effects of stimulation on the phase of stimulation. The model is fitted to the PRC and can reasonably predict the ARC, and notably what is approximately the best phase to stimulate (section 5). We then proceed to compare the relationship between response curves in the linearised and the full model and conclude that non-linearity is important to better reproduce the relationship found in patients (section 6). Finally a discussion is provided (section 7).

Patient response curves and their relationship
We begin by computing response curves from patient data, providing a statistical analysis of the curves, and describing the PRC-ARC relationship when applicable.

Analysis method
In [10], ET patients are fitted with an accelerometer to record their tremor and DBS locked to the phase of tremor acceleration is provided in blocks of 5 s to the VIM of the thalamus, with 1 s without stimulation between blocks. Each block targets a stimulation phase randomly selected out of 12 tremor phases (e.g. 7 degrees for block number 6 in Figure 1). Stimulation is delivered once per period at the target phase, in the form of a burst of 4 to 6 pulses at high frequency (130 Hz or higher). There are about 10 trials available per phase (about 120 blocks per patient). Patient's PRC and ARC are obtained via the "block method" as in [10], which was specifically developed for this type of data, and is implemented as follows. Tremor frequency is around 5 Hz, and the dominant axis tremor recordings are bandpass-filtered (4 Hz band encompassing the patient tremor frequency content) by means of a backwards and forwards Butterworth second order filter (zero-phase filtering) and z-scored.
Change in phase For each block, a straight line is fitted to the evolution of the Hilbert phase during the 1 s period without stimulation before the block. The change in phase ∆ϕ due to the block is given by the difference between the phase of the fitted reference line evaluated at the end of the block and the actual Hilbert phase at the end of the block (see Figure 1). This phase response is divided by the number of pulses in the block (on the basis of 4 pulses per burst for patient 4R and 4L, and 6 pulses per burst for the rest). The target phase at which stimulation is supposed to occur is known for each block, but phase tracking not being perfect, the actual Hilbert phase at which stimulation occurred is determined for each burst of stimulation as the circular mean of the Hilbert phase during the burst. We take the circular mean of these burst angles for a given block as the actual mean phase of stimulation for the block. These values are then binned into 12 phases bins, and the change in phase is averaged within bins.
Change in amplitude For each block, the change in amplitude is given by the difference between the mean of the Hilbert amplitude during the last second of the block and the mean of the Hilbert amplitude during the 1 second without stimulation before the block (see Figure 1). This amplitude response is divided by the number of pulses in the block, and averaged across blocks in the same actual mean phase of stimulation bin as for the change in phase.
Measuring response curves significance and PRC-ARC shift We performed two statistical analysis on patients' response curves. First, PRCs and ARCs were tested for a main effect of phase by means of a Kruskal-Wallis ANOVA (12 phase bins) to differentiate patients' response curves that may be dominated by noise. Second, since we are expecting response curves to have a dominant first harmonic, the cosine model y = c 1 + |c 2 | cos(x + c 3 ) was fitted to patients' phase response and amplitude response curves. We assessed via F-tests whether the cosine model was better at describing the data than a straight line at the mean. Including the less specific ANOVA test allows for more generality, as we do not wish to exclude patients with significant, but non-sinusoidal response curves. On the other hand, the cosine test is a more powerful test for sinusoidal curves. We therefore define the following criterion for selection of a patient for further study in the rest of the manuscript. Significance criterion: having both PRC and ARC deemend significant under FDR control (see below) by at least one of the two tests -ANOVA test for a main effect of phase or cosine model F-test.
In both cases, the adaptive linear step-up procedure modified by Storey et al. in [26] and reviewed in [27] was used to control the false discovery rate (FDR) so that F DR ≤ 0.05. It improves on the original Benjamini and Hochberg procedure [28] by using an estimatorm 0 for the number of true null hypothesis m 0 (total number of tests m = 12 for each analyses). Controlling the FDR at 5% guarantees that the expectation of the number of false positives over the number of positives is less than 5%. Additionally, in datasets where both PRC and ARC are significant according to the cosine F-test, the shift in phase between the cosine model fits to the PRC and the ARC is calculated as

Results of the analysis
We analysed 6 datasets from 5 patients (patient 4R and 4L are for the right and left upper limbs of the same patient). PRCs and ARCs obtained are shown in Figure 14 in Appendix, and results of the statistical tests are presented in Table  1. Based on the significance criterion defined in the previous section, patients 1, 5 and 6 are selected for further study, as both their PRCs and ARCs are found significant by the cosine F-test under FDR control. We note that patient 5 also has both his response curves deemed significant by the ANOVA test under FDR control. Datasets 3, 4R and 4L do not satisfy our selection criterion. In Figure 2, the shift φ P RC − φ ARC is plotted for patients for whom the cosine model was deemed significant in describing both their PRC and ARC (which happens to be the same subset as patients satisfying our significance criterion). They have a shift in π 2 , π , patients 5 and 6 being quite close to π 2 .

Implementation of the Wilson Cowan model
We introduce in this section the WC model as implemented in the present work. We map a 2 population WC model without delays as described in [29] onto the anatomy of the Thalamus (Figure 3). The circuit we are about to describe is a good candidate, but not the only biologically plausible mapping of an E/I loop in the context of tremor. In our candidate mapping, the VIM is modelled as an excitatory population E, connected to an inhibitory population of the thalamus I, the reticular nucleus (nRT). The high coherence between ventral thalamic activity and electromyographic recording of the contralateral wrist flexors [11,12] is our justification for modelling tremor by the activity of the E population. VIM and nRT are reciprocally connected (the excitatory projections from VIM to nRT are via Cortex). The VIM receive a constant input from the deep cerebellar nuclei (DCN) and is part of a self-excitatory loop via Cortex. nRT receives a constant cortical input. We add Gaussian white noise to this two-population WC, and the activity of the VIM, E, and the activity of the nRT, I, are described by the following stochastic differential equations (SDE): where dW E and dW I are Wiener processes, σ the noise standard deviation, and with w P R the weight of the projection from population "P" to population "R", θ P the constant input to population "P", τ a time constant (assumed to be the same for both populations). f is the sigmoid function used in [29] (parametrized by a steepness parameter β): The VIM is the most common target of DBS for ET, which why we model stimulation as a direct increase in E. Analytical expressions for response curves are out of reach for the full non-linear model, and we will study next a linearisation of a deterministic stable focus model. Such a linearisation can be applied to the deterministic WC model (with σ = 0 in equation (1)) in the focus regime.

Response curves and their relationship in a focus model
The goal here is to derive approximate analytic expressions for the first order phase and amplitude responses to one pulse of stimulation in a 2D dynamical system that is described by a (stable) focus. We follow the previous section in modelling the tremor signal as the first coordinate of the dynamical system, and in providing stimulation pulses along the first dimension. The results will provide a basis for understanding how the effects of stimulation on phase and amplitude are coupled in the WC model, and for comparison with experimental data.

Linearisation for a focus
To distinguish scalars and vectors more easily, vectors will be denoted in bold. Leṫ Z = F (Z) be a dynamical system, with Z ∈ R 2 . Let J be the Jacobian of the system (assuming F is differentiable): We will assume that F creates fixed point dynamics, let Z * be the fixed point. The dynamics of X = Z − Z * are approximated in the vicinity of the equilibrium X = 0 by the linear equatioṅ where J(Z * ) is the Jacobian evaluated at the fixed point. We will treat the case of Jocabians having complex conjugate eigenvalues λ ± = σ ± iω. In particular, we are interested in stable foci, which imply σ < 0. The WC model can operate in that regime [29]. The case of centers (σ = 0, purely imaginary complex conjugate eigenvalues) will also be described, although it is of little interest for patient fits. If k = a + ib is the right eigenvector associated with λ + , the real valued solution of (3) reads as follows: We will be using the following notations for the coordinates of the eigenvector k: Equation (4) and what follows are not valid in the case of multiple or repeated real eigenvalues, which are of no interest for our purposes (no rotation).

Phase definition
The notion of phase is central to phase-locked stimulation, and phase can be defined in various ways: asymptotic phase (based on isochrons), geometric phase in the phase plane for 2D systems, or Hilbert phase to cite a few. A typical signal only has one component, and the Hilbert transform provides a convenient way of reconstructing a phase from a single component. The Hilbert phase is therefore widely used to analyse experimental data, and this is the reason why we choose in our linearised system a phase definition approximately equivalent to it. We define a phase variable as φ = ωt with a zero phase point defined as the maximum of X 1 (t) (similarly to the Hilbert phase). This phase definition is different from the trajectory polar angle and from usual asymptotic phase definitions, and we will demonstrate next that it is very close to the Hilbert phase of X 1 for slow decay. It should be noted that this is generally only true for the linearisation. As the Hilbert phase is also the phase definition used in the other sections of this manuscript, the following proof will ensure consistency.
To establish equivalence with the Hilbert phase, a first step is to calculate the Hilbert transform of our signal X 1 (t). The Hilbert transform is a linear operator, and X 1 (t) is a linear combination of s(t)s c (t) and s(t)s s (t). So let us consider s(t) = e σ|t| , s c (t) = cos ωt, and s s (t) = sin ωt. We will start by showing that H(s(t)s j (t)) ≈ s(t)H(s j (t)) for j = c, s. The Bedrosian identity [30] states that the Hilbert transform of the product of a low-pass and a high-pass signal with nonoverlapping spectra is the product of the low-pass signal and the Hilbert transform of the high-pass signal. The spectrum support of s is R, but for low decay compared to the rotation, the spectrum of s is very small where it overlaps with the spectra of s c or s s . The equality given by the Bedrosian identity turns into an approximation, and inspired by the proof in [30], we can calculate an error term. Let S and S c be the Fourier transforms of s and s c respectively: The Fourier transform of s c is given by where Γ(u) = π i sgn(u + ω)e iωt + sgn(u − ω)e −iωt , which can be simplified as follows: The Fourier transform S(u) = 2σ σ 2 +u 2 is even, therefore with H(s c (t)) = sin ωt with H(s s (t)) = − cos ωt Numerical integration proves that for ω |σ|, and in particular in the case of the patients we are interested in, I sc and I ss are under a couple of percent of the signal scale for enough periods for our purposes (see Figure 15 in Appendix) -only one period is needed to derive response curves. It is therefore reasonable to ignore I sc and I ss , and the Hilbert phase of X 1 is given by In the case of ω |σ|, it can be shown that α and β are such that φ Hilbert is referenced to the maximum of X 1 .

Reference trajectory and stimulated trajectory
We will consider a reference trajectory without stimulation, and a trajectory that underwent an instantaneous stimulation pulse δX 1 at a stimulation phase φ 0 (See Figure 4). The effects of stimulation on phase and amplitude will be measured at the next maximum of X 1 for both trajectories to obtain the PRC and the ARC (first order response). Expressions for the coefficients K ref and K ref of the reference trajectory are derived in Appendix A. We want to study the effects of stimulating at the phase φ 0 . The point of stimulation X 1 − at phase φ 0 is expressed as An instantaneous stimulation δX 1 is applied at X 1 − as follows: The trajectory after stimulation is still constrained by the dynamics given by equation (4), which allows for expressions for the coefficients on this new trajectory K stim and K stim to be found (see Appendix B). To measure the change in phase and amplitude between the next peaks of the stimulated trajectory and the reference trajectory, the phase φ max of the next maximum of the first coordinate on the stimulated trajectory X stim 1 is needed (the phase of the next maximum of X 1 on the reference trajectory is 2π). A derivation for φ max is provided in Appendix C.

Phase response
The first order phase response curve (PRC) can be calculated based on the reference trajectory period T 0 and the stimulated trajectory period T stim , which is given by the sum of the time spent on the reference trajectory before stimulation and the time spent on the new trajectory after stimulation.
For a PRC in radian, we obtain φ max depends on δX 1 through K stim and K stim , and a software assisted Taylor expansion around 0 yields, to lowest order in δX 1 : The PRC is found to be proportional to the inverse of the peak amplitude of the oscillations at the beginning of the stimulation period X 0 1 and to the stimulation amplitude δX 1 . It directly depends on phase via sinusoidal functions and a factor Figure 4: Illustration of the approach taken to derive expressions for the phase and amplitude response in the linearisation of a 2D focus model. Top: phase plane, bottom: time-series of X 1 . The tremor is modelled by X 1 , and the stimulation δX 1 is applied to X 1 . The system shown corresponds to the linearised fit of patient 1 as described in section 6.1. related to the decay. A, B, and C only depend on the real and imaginary parts of the eigenvalue λ + and the associated eigenvector k.

Amplitude response
For our purposes we are interested in the amplitude of the first coordinate, and the first order amplitude response curve (ARC) is obtained as the difference in first coordinates between the stimulated and the reference trajectories evaluated at their respective next peak after stimulation. It should be noted this is equivalent to a first order change in Hilbert amplitude, at least for ω |σ|.
A software assisted Taylor expansion around 0 yields, to lowest order in δX 1 : Interestingly, the ARC close to the fixed point does not depends on the amplitude of the oscillations X 0 1 . As expected, the ARC is directly proportional to the stimulation amplitude δX 1 . Similarly to the PRC, it directly depends on phase via sinusoidal functions and a factor related to the decay, and D only depends on k. The obvious similarities between the PRC and the ARC suggest there may be a relationship between the two.

Relationship between PRC and ARC
From equation (17), the derivative of the PRC is given by Let us study the case where ω |σ| Therefore in that case the ARC is approximately the opposite of the derivative of the PRC scaled by the peak amplitude at the beginning of the stimulation period (in general, the scaling factor is F X 0 1 ):

Applications to simple systems
We turn to simple systems to illustrate the results of the previous sections. In what follows, response curves are given for δX 1 = 2 × 10 −4 and X 0 1 = 10 −3 (these only act as scaling factors of the response curves and will not change their shape).
Circular flow without decay As an introductory example, let us consider a simple circular flow for which the matrix J is as The eigenvalues of J circ are ±i so the results of the previous sections can be applied.
Equations (17) and (19) are simulated for this system with our choice of δX 1 and X 0 1 . The result for the PRC is shown in Figure 5, panel A2, and for the ARC in panel A3. For this system, the PRC is simply the opposite of a sine, the ARC simply a cosine. σ = 0, hence G = D (see section 4.6) and equation (21) is exact, as exemplified in Figure 5, panel A3. The ARC is obtained by only scaling the derivative of the PRC by −X 0 1 as a 2 = b 1 = 0 and F = 1. Note that WC parameters for which the system's Jacobian at the fixed point is J circ cannot be found as the second diagonal term cannot be 0, at least in the version of the WC model used in this work (see (39) in Appendix).
Circular flow with decay We can introduce a slow decay ( Figure 5, panel B) and then a fast decay ( Figure 5, panel C) in the circular flow.
The slow decay leads to a scaling factor F ≈ 1, and the approximation of equation (21) is very good, as ω |σ| (see Figure 5, panel B3 -ω = 200|σ|). The case of the fast decay corresponds to ω = 5|σ|. The PRC and ARC no longer look like pure sinusoids and the approximation relating PRC and ARC is less accurate (Figure 5, panel C3). It is possible to find WC parameters for which the system's Jacobian at the fixed point is J slow circ or J f ast circ . How such parameters are found is explained in Appendix D, and the results are presented in Table (3) in Appendix. In both cases, w IE = w IE , and w EE = 0.  (17). Third column: ARC as per equation (19) and opposite of the derivative of the PRC scaled by F X 0 1 . Panel A corresponds to J circ (circular flow, no decay), panel B to J slow circ (circular flow, slow decay), panel C to J f ast circ (circular flow, fast decay), and panel D to J ellip (tilted elliptic flow, no decay) Tilted elliptic flow without decay The tilted elliptic flow without decay of Figure  5, panel D, corresponds to the J matrix The PRC and ARC are sums of a sine and a cosine, which brings a horizontal shift in phase compared to circular flow without decay. The eigenvalues are still purely imaginary, but F is no longer one. Because σ = 0, the relationship of equation (21) is still exact (see Figure 5, panel D3). It is possible to find WC parameters for which the system's Jacobian at the fixed point is J ellip (see Table (3) in Appendix).
Patient fits fall in the category of (potentially tilted) elliptic flows with decay, and will be dealt with in section 6.1. The linearised stable focus model exhibits close to sinusoidal response curves and a PRC-ARC shift close to π 2 , which when contrasted with patient data (response curves passing the cosine model F test and PRC-ARC shifts in π 2 , π as shown in Figure 2), provides a strong motivation to fit the WC model to data.
5 Fitting the full Wilson Cowan model to patient data and response to phase-locked stimulation

Fitting procedure
We now turn to fitting our stochastic neural mass model (equation (1)) to patient data. The model is fitted to features extracted from patient tremor recordings. The parameters we fit are shown in Table 2, and include model parameters, stimulation magnitude and stimulation delay (time between the stimulation trigger is recorded and stimulation is actually provided to the E population). Stimulation is implemented directly in the Euler update of our integration scheme. We aim at reproducing tremor dynamics and fit to 3 dynamical features: the power spectrum density of the data (PSD), its Hilbert envelope probability density function (PDF), and its Hilbert envelope PSD. While the envelope PDF captures the range of amplitudes present in the tremor, the envelope PSD describes how quickly tremor amplitude changes. But we also aim at reproducing response to stimulation, and fit to the patient PRC. The data dynamical features are obtained after filtering and z-scoring the data as described in section 2.1. The data PRC is obtained as in section 2.1. The fitting procedure is as follows (summarized in Figure 6). For each patient, random sets of parameters are picked from uniform distributions (bounds in Table  5 in Appendix). To improve the efficiency of the optimisation, we accept parameters only if the PSD peak of the corresponding model (without stimulation) is within 1 Hz and 25% in magnitude of the data PSD peak. Once 2500 parameters have been accepted, we put them through local optimisations. Local optimisations are carried out using gradient free optimisation, specifically a direct search algorithm called the generalized pattern search algorithm (implementation details in Appendix E). Parameters are put on a similar scale to improve search robustness, and hard limits are given to the optimiser (see Table 5 in Appendix). Optimisations were performed in parallel on a supercomputer. A time step of 1 ms was used for the fits (a period is about 200 ms). At the end of this process, the 20 best performing sets of parameters were put through more local optimisations with a finer time step of 0.1 ms and stop criteria leaving room for more steps. Finally, the best performing set of parameters is selected from the results of this second optimisation as the highest R 2 (defined later in this section). The finer time step is also used to produce the results shown in section 5.2).
In order to measure response to stimulation as in the data, each local optimisation step needs to simulate the model with phase-locked blocks of stimulation. This requires integrating the differential equations of the model while tracking the phase and providing stimulation at the right time, which is done by monitoring the zerocrossing phase alongside a Euler integration scheme. One simulation consists of 600 trials with 12 blocks of phase-locked stimulation each. As in the experimental paradigm, blocks last 25s, and inter-block intervals are 1s. Inter-trial intervals are 5s, and the first trial starts after about 200 periods. During this initial period, the mean of E and the standard deviation of E, σ sim , are obtained from about 20 periods after a ramp-up time of about 40 periods. Phase-tracking subsequently starts: E is centered and a threshold T = 0.2σ sim is used to track positive zero-crossings. The use of hysteresis via a threshold was found critical to handle the noise included in the model. We define a positive zero-crossing as happening when These conditions are constantly monitored, and if found true, a positive zerocrossing is declared to have happened at time step χ = n+p 2 . We evolve the zerocrossing phase according to a frequency based on the previous period, and if χ k is the last positive zero-crossing to have occurred, the current value of the zero-crossing phase is given by If the value of 2π is reached, the phase value is set to 0 until the next positive zerocrossing is detected. Stimulation is provided after ϕ reaches the target phase for the block, and the stimulation trigger is recorded ∆t stim before stimulation occurs. If the zero-crossing phase hasn't reached the target stimulation phase yet when the next positive zero-crossing is detected, stimulation is provided right then. As in [10], a pulse of stimulation consists of 6 quick bursts at 130 Hz.
The 4 features are computed on the model output at each optimisation step. The same method is used as for the data features, with three differences. The first is that for increased stability of the optimisation, the model PRC is averaged over a much greater number of trials (600 trials), while the more robust dynamical features are obtained from 9 trials only to reduce computation cost. The second is that the model output is not filtered to compute the dynamical features (only z-scored), as we want the model output to primarily generate the filtered tremor signal (a model generating mostly 1Hz oscillations but reproducing patient tremor when filtered at 5Hz wouldn't be desirable). Computing the PRC still requires filtering, as it relies on the Hilbert transform. The third difference is that the filtering window for the PRC can't be adjusted manually in optimisation steps, so a 4Hz band centered on the model PSD peak is used.
The zero-crossing phase is only needed to enable phase-locked stimulation in the model, and the actual Hilbert phase at which stimulation occurred is used to compute response curves via the re-binning process described in section 2.1. Phasetracking performance is illustrated in Figure 16 in Appendix. The simulation of the model at each optimisation step requires to track the zerocrossing phase in order to provide stimulation at the right phase. The phase-tracking ability of the scheme is satisfactory when compared to the actual Hilbert phase (left, detailed in Figure 16 in Appendix). The optimiser minimises a cost function that includes the comparison of 3 tremor dynamics features (tremor PSD, tremor envelope PSD, tremor envelope PDF) plus the PRC against the data (middle). Following a second optimisation of the 20 best results with a finer time step, a best set of parameter comes out of the procedure, and the model ARC can be compared against the data ARC.
At each step, once the 4 features have been computed on the model output, the following cost is returned to the optimiser: With y i , i ∈ {1, 2, 3, 4} being the 4 features considered. The fit with the lowest R 2 = 1 − c for each patient is deemed the best fit. In case of a tie (difference in cost lower than sem error bars), foci are preferred over limit cycles.

Results of the fits
Patients passing our significance criterion (section 2) are fitted to, namely patient 1, 5, and 6. For these patients, we find that the model successfully reproduces tremor dynamics, including tremors with sudden bursts, and can fit to patient phase response to stimulation. The best fits obtained upon completion of the optimisation procedure are shown in Figures 7, 8, and 9. In addition to reproducing tremor dynamics and being able to fit to patient PRCs, the model seems to be able to reasonably predict patient ARCs (obtained as in section 2.1, but not fitted to), and in particular which phases are approximately the best phase to stimulate, i.e. the phases at which the maximum decrease in tremor happens. Because of averaging   Checking fitted stimulation magnitude Cagnan et al. [10] report what the device settings are, and in particular the total electrical energy delivered (TEED) per unit time. We can build a quantity based on the fitted stimulation that should scale with the TEED per unit time. Because of z-scoring along the E dimension, we have to divide the fitted stimulation which is measured along the E dimension by the standard deviation of E (before z-scoring), and turn the fitted stimulation into an effective stimulation. As bursts are delivered once per period, this effective stimulation should be multiplied by the mean frequency of E to obtain a quantity proportional to energy per unit time (the number of pulses per burst is the same for the 3 patients). Figure 10 shows the effective stimulation times the mean frequency  for the 15 best performing fits against the TEED per unit time for each patient (correlation coefficient for fit averages r = 0.98). Under the assumption that patient intrinsic sensitivities to stimulation are somewhat similar, we can conclude from the correlation that the fitting procedure successfully captures the scale of stimulation across patients. PRC-ARC shift in WC synthetic data The PRC-ARC shift is computed on WC synthetic data with phased-locked blocks of stimulation generated by the full model fitted to each patient. This time we can take full advantage of the model and compute response curves from more trials than for patient data or model data in optimisation steps, and perform 10 repeats of 600 trials for the top 15 fits for each patient. The PRC-ARC shift is then measured as in section 2.1 for each of the 10 repeats, and shown in Figure 11. The large filled circles represent the mean of the 10 repeats for each patient fit. It appears that the PRC-ARC shift obtained for synthetic data of top patient fits mostly lie in the upper-left quadrant of the unit circle for all 3 patients π 2 , π , similarly to patient data. One fit to patient 6 is an outlier in terms of its shift, due to high effective stimulation. While the model can allow for a larger shift than π 2 , this is not the case for the linearised model, and the difference is the focus of the next section.

PRC-ARC shift in the model
The linearised model makes different predictions for patient fits than the full model, in particular in terms PRC-ARC shift. The present section will look at the linearisation of patient fits, and then contrast it with the full model.

Relationship between analytic response curves in the linearised fitted WC models
The PRC and ARC expressions derived in section 4 can be applied to the linearisation of the best WC models fitted to data from the 3 selected patients (Jacobians at the fixed point shown below). In the fits b 1 = 0 or b 2 = 0, which marginally  (17) and (19). The curves obtained are shown in Figure 12. The same values as in section 4.7 are used for X 0 1 and δX 1 . Note that the stimulation delay ∆t stim is not shownit affects both the PRC and the ARC and does not play a role in the PRC-ARC shift. More interestingly, we observe that ω |σ| in the 3 fits (see Table 4 in Appendix), suggesting that the PRC-ARC relationship described by equation (21) should approximately hold. This is indeed the case as shown in the third column of Figure 12. |σ| is higher for patient 5 and as expected, the approximation is slightly worse for this patient (panel B3 in Figure 12). For small stimulation, the deterministic picture with patient parameters close to the fixed-point is that the PRC-ARC shift should be very close to π 2 . In what follows, we investigate the difference between this idealised picture and what is observed in synthetic data.

Accounting for the difference in shift between focus model analytic expressions
and WC synthetic data Four factors could account for the difference in PRC-ARC shift between the idealised picture given by analytic response curves with patient parameters (previous section) and what is observed in WC synthetic data (section 5.2). First, stimulation may be large enough that the Taylor expansions used to derive the analytic PRC and ARC expressions are no longer approximately valid. Second, tremor in patient fits may correspond to a regime not so close to the fixed point, compromising the linearisation validity. Third, the introduction of noise in the model may result in effects on the PRC-ARC shift that do not average out to zero. Fourth, in synthetic data, the response to stimulation is measured by the block method, which differs Second column: PRC as per equation (17). Third column: ARC as per equation (19) and opposite of the derivative of the PRC scaled by F X 0 1 . Panel A, B, and C respectively correspond to patient 1, 5, and 6. from the first order approach taken in our derivations. We next show that for the 3 best fits considered non-linearity is the main driver.
Ten repeats of 600 trials of synthetic data are generated for the linearisation of the best fits to each patients. The integration scheme with live phase tracking and stimulation is the same as described in section 5.1, only the SDEs are now Where dW E and dW I are Wiener processes, σ the noise standard deviation (same values as in the non-linear case), E * and I * are the coordinates of the fixed point, and J is the Jacobian at the fixed point of the patient fit. The same values as in the non-linear case are used for the stimulation magnitude and stimulation delay, with the exception of patient 5, for which the stimulation magnitude is set to a fifth of its value in the non-linear case, as higher values were seen to cause a breakdown of phase tracking, and result in unreliable response curves. For each patient and for each of the 10 repeats, response curves are obtained via the block method and the PRC-ARC shift is then measured as in section 2.1. The results are shown in Figure 13 (middle), alongside the shifts measured from the response curves presented in section 6.1 (left), and the shifts measured in the full WC model (right). It can be seen that going from the analytic response curves to the linearised model (i.e. adding noise, measuring the response to stimulation via the block method and not a first order method, and using a finite stimulation magnitude rather than a infinitesimal stimulation), doesn't affect the shift much. However, a substantial increase in the shift is obtained by introducing the non-linearity, which brings the shift in the upper-left quadrant, where patient data lie. Figure 13: Non-linearity accounts for most of the difference in PRC-ARC shift seen in synthetic data (middle and right), when compared to the PRC-ARC shift derived in the focus model (left). When computed from synthetic data, the PRC-ARC shift of all 10 repeats is shown (smaller circles), as well as the repeat mean (larger circles). One repeat corresponds to 600 trials, only showing best fits for each patient.

Discussion
We showed that in a 2D linearised stable focus model, PRC and ARC are close to sinusoidal, in particular for small decay. Moreover, the PRC-ARC shift is close to π 2 . Half of the patients in our dataset had significant sinusoidal PRCs and ARCs (an effect of phase could not be found in other patients in at least one of their response curves), and the significant patients have a PRC-ARC shift in the interval π 2 , π . A full WC model can be fitted to tremor dynamics features and to the PRC for these patients, and as hinted at by the similarities seen in the linearised focus model and the data, the best fits -found to be stable foci -can reproduce the dependence of the effects of stimulation on the phase of stimulation. The best fits also reasonably predict the ARC, and notably what is approximately the best phase to stimulate. Compared to the 2D linearised focus, the non-linearities of the full WC model allow for a better reproduction of the phase dependence found in patient data, in particular as far as the PRC-ARC shift is concerned. Our full model describes the neural populations involved in ET, which, together with its success in reproducing phase response and predicting amplitude response in patients, makes it a strong candidate for further study of phase-locked DBS.
Phase definition While an asymptotic phase definition is common in theoretical studies, experimental studies tend to favour finite time phase definitions such as Hilbert phase. To reproduce the data, an instantaneous phase seems more appropriate than an asymptotic phase, as there is no indication of stimulation happening on or close to an attractor. It has been shown recently in [31] how an operational definition of the phase can describe transient spiking, when an asymptotic phase does not capture the phase dependence of transients. In this study, our phase definition is the Hilbert phase of the tremor data. It is therefore referenced to the maximum of the tremor oscillations (represented by the first coordinate of the dynamical system), and does not require a limit cycle. The Hilbert phase is an angle in the analytic signal space, and is a protophase [32] that does not generally grow linearly with time. This is not a concern from the perspective of describing patient data, as this is the observable choice we are making for both the data and the model. Commonly used with data, the Hilbert transform has also been proposed as a robust method to measure steady state PRCs in single neuron models [33]. Moreover, stimulation is assumed to be small in our analytical expressions (section 4), but not in the full model, contrary to standard asymptotic phase reduction strategies.
Linearisation The response curves derived for the linearisation of a 2D focus in section 4 can be related to previously published expressions. In particular, the infinitesimal PRC for radial isochron clocks has been derived in [34], and has been recently included in [35] under the larger umbrella of general radial isochron clocks. The radial clock case (K(φ) = ω in [35]) perturbed along the first dimension agrees with our equation (17) for the case of a circular flow (see section 4.7). For this simple system, the asymptotic phase response is the same as the first order Hilbert phase response.
Moreover, for small decay, the best phase to stimulate corresponds to the maximum positive slope of the PRC in the response curves derived. In fact, the ARC is simply a scaled version of the opposite of the PRC derivative. A similar relationship has been first reported in a theoretical study in the context of an individual oscillator [36], and more recently in [14] in the context of population response curves of a Kuramoto model. It is noteworthy that a similar result is found for 2D focus models with slow decay, and in particular for the linearisation of the WC model, another popular neuroscience model very different in essence from coupled oscillator models.
Our derivations do not assume proximity to a limit cycle, and this allows the study of the dependence of the response to stimulation on the amplitude of the oscillations for a given model (limit cycles don't have an amplitude variable in the case of infinitesimal perturbations). In the linearisation, the PRC is found to be inversely proportional to the amplitude of the oscillations before stimulation (see equation (17)), while the ARC does not dependent on it.
Because the block method phase and amplitude response used in the rest of paper are normalised by the number of pulses and blocks are only about 25 period long, it seems legitimate to think that, although they are different objects, the first order response to a single pulse and the block method response could be related, and in particular that they might share similar PRC-ARC relationships. Part of the connection hinges on our proof that the phase definition in the linearisation of the focus model overlaps with the Hilbert phase when the decay is small compared to the rotation (section 4.2). And indeed, the PRC-ARC shift predicted by our expressions derived for the first order response to one pulse of stimulation in a linearised focus is very close to the shift obtained by the block method on linearised WC synthetic data (see Figure 13). Our analytical derivations provide a rationale to fit the full WC model to data and an intuition for why the model can predict patient ARC, but do not offer an exact analytic treatment of the block method. Specifically, individual pulses in a block may have different effects depending on where they are located in the block and depending on stimulation history within the block [10].
Fit and shift The full WC model is fitted to data with Gaussian white noise (equation (1)). The best performing fits are stable foci for all 3 patients, and very few limit cycles are found in the top 15 fits for all 3 patients -1 for patient 1 (shares the 1st place with a stable focus -distance between costs only 30% of sem error bars), 1 for patient 5 and none for patient 6. In the stable focus regime, noise brings the system away from the stable fixed point, and the interaction of the noise with the dynamics of the system makes the reproduction of patient tremor possible. In the absence of noise, the system would converge to the stable fixed point and no tremor would be generated, so symptoms are related to the noise level in this model. Instead of noise, tremor-like activity may be obtained by exploiting chaotic dynamics arising from coupling several WC models together [37], but this would significantly increase the complexity of the model (more on increasing complexity in the last part of this section).
In fitting our thalamic model to tremor acceleration, we are assuming thalamic activity and tremor are directly related as mentioned before (see section 5.1). Tremor activity is however expected to lag thalamic activity due to conduction delays. In addition to that, thalamic activity itself may be affected by stimulation with a delay. In the model, we allow for a stimulation delay ∆t stim between the stimulation trigger and the time when stimulation is actually delivered to the E population. This parameter is fitted to the data, and can account for both delays, which are virtually indistinguishable with the data available. It's interesting to note that the stimulation delay of the best performing model for patient 5 is longer than one period (see Table 2). This is found consistently in the top 3 best fits, and reducing the delay to its value modulo the average period substantially reduces the quality of the PRC fit. Besides this short term delay, our model does not include medium or long term plasticity effects, which are not expected to be strongly present in the recordings as stimulation is only delivered for periods of 5 seconds in a row. In our model, stimulation is provided to the E population via a direct increase in the population activity. While stimulation is provided via the sigmoid function of the E population in other studies [17], we found this approach too restrictive due to sigmoid saturation, and inadequate to reproduce the full extend of the response to phase-locked DBS in some patients. As a reminder, the choice of stimulating the E population rather than I is made for biological consistency, as the VIM is the most common stimulation target in ET DBS.
Fits were performed using the generalized pattern search algorithm on many sets of random initial parameters. This approach was chosen for its robustness and computational efficiency in a non-smooth landscape with 4 non-linear features and 10 parameters, despite requiring the use of a supercomputer. In particular it has been deemed superior in finding better fits to the simplex algorithm. The implementation used also has the additional benefit of being able to handle failed simulations (which occasionally happen as response curves with 12 phase bins can't be obtained for some parameters). However the fitting procedure results in many "good" local optima. What these "good" sets of parameters have in common and what they can tell us about the patients we are fitting to is not easily addressed with our current fitting strategy. Even real biological networks may have redundancies, and may exhibit the same behavior under different network configurations. Recent developments in Bayesian inference methods such as [38,39] allow to recover the posterior distribution over parameters, hence to answer the question what is the space of parameters consistent with the data. Whether Bayesian inference could successfully tackle a complicated landscape and provide more meaningful insight on fitted model parameters in the setting of the present work is a promising avenue for further research. Additionally, the trade off between features could be studied from the perspective of a Pareto front. A limitation of our fitting method is related to the integration scheme: to reduce computation cost, the Euler step used in the first optimisation process is 10 −3 . The top 20 best fits are then re-optimised based a Euler step of 10 −4 , and results are produced with this finer time step, as dynamics can be qualitatively different (further reduction in the Euler step does not change the dynamics). While the necessities of phase-locked stimulation precludes the use of built-in, powerful integration schemes, a more advanced integration scheme could remove the need for a second optimisation while keeping the computation cost down. The performance of our simple phase-tracking strategy is good for patient 1 and 6 and satisfactory for patient 5 (see Figure 16 in Appendix). Response curves are obtained based on the actual Hilbert phase of stimulation in a post-hoc manner, which makes up for the reduced performance observed for patient 5. Still, more accurate algorithms could be explored. Our zero-crossing strategy would benefit from a better frequency estimate for the current period (currently simply the period of the previous period) and more robustness to noise. The proposition of [40] involving autoregressive forward prediction and the Hilbert transform is attractive, although its computational cost may be high, and some parameters need to be adjusted for each time series.
The success of the WC model in predicting patient ARCs when fitted to their PRCs is partially explained by its ability to modulate the PRC-ARC shift. The PRC-ARC shift in the full model can reach the range found in patients while the linearised version of the WC is limited to the close vicinity of π 2 ). The response curves of the full WC model are also better at reproducing the data and can vary from pure sinusoids. However there is still some room for improvement in reproducing the shift, in particular as far as patient 1 is concerned (patient shift quite a bit larger than the model). The model can allow for a larger shift as shown by a fit selected in the top 15 shown in Figure 17 in Appendix. The PRC-ARC shift could be selected as an additional feature to fit to to improve ARC reproduction.
In terms of describing the neural circuits involved in ET, the model could be further improved by including the cortex and the DCN as populations of their own, which would make the model 4 dimensional. As suggested in [17], the inferior olive which provides input to the DCN could also be modelled, and the spatial extent of the VIM could be accounted for by splitting it in 2 populations or more. Increasing the number of populations would however increase the number of parameters of the model, and make the optimisation process more computationally intensive, and the model more prone to over-fitting. The model seems to be able to reproduce the data in its current state, which suggests an increase in complexity is not warranted.

Conclusion
The focus WC model with noise can be fitted to ET patients with significant sinusoidal response curves (no effect of phase was found in at least one of the response curves of other patients). The model reproduces the phase dependence of the response to stimulation as well as predicts the amplitude response to stimulation, which directly relates to tremor reduction. Phase-locked stimulation promises less stimulation hence less side effects, for the same clinical benefits, which would be highly desirable for patients. Our study positions the WC as a strong candidate to model the effects of phase-locked DBS. Its ability to describe all patients with both response curves significant in at least one of our tests should be re-assessed as more data becomes available, both in terms of number of patients and recording length. Phase dependent activity is thought to play a central role in physiological information processing [41,42], and in our analytical derivations, the phase of the linearised model was defined in a way that does not depend on modelling oscillations by a limit cycle, and that for small decay overlaps with the Hilbert phase, which is widely used in experiments. While our analytical derivations have been pursued for a first order response to a single pulse, the experimental measurement paradigm, the block method, was used throughout the rest of the paper. Our linearised model results have delivered an informative approximation of the block method response, but it remains to be seen whether different approaches could provide additional insights on the block method response.

Appendices
We include in appendix details of the derivations leading to response curves analytical expressions in the linearised system (Appendices A to C), the procedure used to obtain WC parameters from a given Jacobian (Appendix D), implementation details of the optimiser used for fitting to patient data (Appendix E), as well as supplementary figures (Appendix F) and tables (Appendix G).

A Reference trajectory without stimulation
Let us find the coefficients K ref and K ref of the trajectory starting at t = 0 at a maximum of the first coordinate X 1 = X 0 1 > 0. With the choice φ = ωt, this will ensure we are referencing the phase to the maximum of X 1 . It should be noted at this point that we are not using the nullcline equations in what follows as we are interested in the dependence of the response on the rotation ω and the decay σ.
From the initial condition at t = 0 : Additionally, X 0 1 being a maximum requires that dX1 dt = 0 at t = 0: Using the condition at t = 0: We are excluding the case where the denominator in (29) is equal to zero, which corresponds to both a 1 and b 1 being zero, which would imply X 1 (t) = 0. Also note that by picking a positive X 0 1 , we are ensuring that the null derivative corresponds to a maximum of X 1 rather than a minimum.

B Trajectory with stimulation
Let us determine what the coefficients K stim and K stim are for the stimulated trajectory (still constrained by the dynamics of equation (4)).

D Finding WC parameters corresponding to a given J
The Jacobian of (1) evaluated at (E * , I * ) can be simplified by making use of f (x) = βf (x)(1 − f (x)). We also have f (Θ 2 ) = I * with We are interested in finding WC parameters so that the linearisation of the WC model at the fixed point will be characterised by a given J matrix If we pick values for β, E * and I * , the remaining parameters can obtained by equating (39) and (40), and by re-arranging equations (37) and (38). This is how parameters were obtained in Table 3.

E Implementation details of optimisation algorithm
The implementation used is Matlab's patternsearch optimiser with the following stop criteria: • main optimisation: mesh size of 10 −4 , function call budget of 800 • second optimisation: mesh size of 10 −5 , function call budget of 1000 Figure 15: Relative error made across patients in estimating H(s(t)s c (t)) by s(t)H(s c (t)) (solid lines) and H(s(t)s s (t)) by s(t)H(s s (t)) (dashed lines). The error is calculated as the ratio of I sc (respectively I ss ) over the modulus of the numerical Hilbert transform of the signal, which is the envelope of the signal. The relative error is under 5% in all cases for at least 12 periods.  Averages are obtained using circular means. The effect of the stimulation delay was removed, and phases are reference to positive zero-crossings. Phase tracking is satisfactory for all patients, although tracking is less precise for later phases in patient 5.