A minimally invasive neurostimulation method for controlling abnormal synchronisation in the neuronal activity

Many collective phenomena in Nature emerge from the -partial- synchronisation of the units comprising a system. In the case of the brain, this self-organised process allows groups of neurons to fire in highly intricate partially synchronised patterns and eventually lead to high level cognitive outputs and control over the human body. However, when the synchronisation patterns are altered and hypersynchronisation occurs, undesirable effects can occur. This is particularly striking and well documented in the case of epileptic seizures and tremors in neurodegenerative diseases such as Parkinson’s disease. In this paper, we propose an innovative, minimally invasive, control method that can effectively desynchronise misfiring brain regions and thus mitigate and even eliminate the symptoms of the diseases. The control strategy, grounded in the Hamiltonian control theory, is applied to ensembles of neurons modelled via the Kuramoto or the Stuart-Landau models and allows for heterogeneous coupling among the interacting unities. The theory has been complemented with dedicated numerical simulations performed using the small-world Newman-Watts network and the random Erdős-Rényi network. Finally the method has been compared with the gold-standard Proportional-Differential Feedback control technique. Our method is shown to achieve equivalent levels of desynchronisation using lesser control strength and/or fewer controllers, being thus minimally invasive.

where W ij = W ji > 0 is the weighted adjacency matrix. The original Kuramoto model is hence recovered for the choice W ij = K/N for all i and j, where K is the strength of the non-linearity and N -the number of oscillators -a normalising factor to make the coupling intensive.
In the case under scrutiny there exist several ways to normalise the interaction term [4], they mainly differ when we are interested in comparing networks with different sizes or considering the behaviour of the system in the thermodynamics limit, N → ∞. For the scope of this work it is enough to normalise the intensity of the interaction term by using the node degree, k i = j A ij , where A ij is the binary adjacency matrix, i.e. A ij = 1 if and only if W ij = 0. The model can thus be rewritten asφ To go one step forward let us define an order parameter, as we did for the Kuramoto model, allowing to rewrite the system in a self-contained way using a mean-field approach. We observe however that in the present case we must use a local order parameter, depending of the node index i: from which one can straightforwardly obtain and thus rewrites Eq. (2) asφ The next step consists in observing that the equation above can still be obtained from a Hamiltonian formalism (see Eqs. (5) and (6) in the main text where A ij → W ij /k i ) and thus one can build a control term in a very similar way. A long but straightforward computation eventually provides the following control term: which can be further simplified after the introduction of We thus eventually get where B (w) i has been defined in the last equality. The control term (8) can be further simplified to make it operational along the same ideas developed in the main text, namely by replacing the phase Ψ i with the angle φ i , and neglect the term B (w) i . We also take into account the presence of an electromagnetic field potential [5,6] created by the controlled microelectrodes acting on the neighbouring neurons. In conclusion, we introduce a new stimulation signal, S stim,w i , generated on the position of the i-th neuron by the potential produced by the microelectrodes located in all the controlled neuronal patches: where r il and c s are respectively the distance of node i from the origin of the electromagnetic field l, and the strength of the potential which in our case is taken to be c s = 1, finally M N is the number of directly controlled nodes. The proposed control strategy taking into account the weights is thus: (10)

SOME MODELS OF WEIGHTED NETWORKS
The proposed control strategy is independent from the network topology considered and the weights distribution. However for a sake of completeness and for the relevance with the proposed application, we hereby consider three classes of weighted networks: an all-to-all connectivity topology, an Erdős-Rényi network [7] and a Newman-Watts [8] coupling. In all the cases the weights are uniformly randomly distributed.
The results reported in the following Fig. 1SM show that the improved control method is able to achieve a strong degree of desynchronisation for non-linear oscillators evolving according to the Kuramoto model. We then define the (de)synchronisation level of the system using the Kuramoto order parameter R(t) = | j e ιφj (t) |/N . All the considered networks are made by N = 100 nodes, M = 20 microelectrodes are used to measure and to control the whole system and the control parameter has been fixed to γ/4 = 20. The natural frequencies are drawn from the normal distribution N (1, 0.01), and the initial angles are uniformly randomly drawn in [0, 2π]. In the left panel we report the case of a complete topology where all nodes are connected with all other nodes and the weights are drawn from an uniform distribution U[0.7, 0.9], and then normalised by the node degree W ij → W ij /k i . One can observe that the original Kuramoto model exhibits a strong synchronisation with the order parameter R very close to one (blue curve), on the other hand the controlled model (red curve) shows a very small R, thus denoting a high level of desynchronisation. In the middle panel, we consider as underlying topology an Erdős-Rényi network with parameter p = 0.8, the probability that a link exists between two nodes; existing links have their weights drawn from the same uniform distribution used above. One can again observe that after a very short transient phase during which the controlled system exhibits an order parameter of the same order of the original Kuramoto model, the controlled system desynchronises (red curve) while the original Kuramoto model does the opposite (blue curve). Finally, in the right panel, we present the results concerning the Newman-Watts network with parameter p = 0.85 and the conclusions are similar to the previous ones, the controlled system (red curve) is strongly desynchronised while the original model doesn't.

THE ROLE OF THE NETWORK TOPOLOGY IN THE DESYNCHRONISATION ISSUE
In the main text, we presented the result concerning the impact of the control parameters M and γ to achieve the desynchronised state, under the assumption of an all-to-all unweighted topology ( Fig. 2 main text, and here reported in panel (d) of Fig. 2SM to help the comparison). The outcome is that a high level of desynchronisation (dark blue) can be achieved using few microelectrodes, i.e. take small values for M , provided a strong enough control is used, namely assume a large enough value for γ. The opposite is also true, one can use a small coupling γ but then more mictroelectrodes are needed. In the remaining panels of Fig. 2SM, we analyse the same question assuming a Newman-Watts topology with binary weights and we study the impact of the parameter p responsible for the creation of shortcuts. Once p is small, the network is close to a 1D-regular lattice (with coordination number k = 6 in the present case), i.e. very few new links are added; for larger p the network becomes more dense and in the limit p = 1 we recover a complete network. For any fixed parameters values (M, γ/4), we report the asymptotic value of the Kuramoto order parameter R = | j e ιφj |/N averaged over 25 independent realisation. We repeated this construction for three Newman-Watts networks, with a high link density (p = 0.85, panel c), with a medium link density (p = 0.5, panel b) and a diluted link density (p = 0.25 panel a). The following colour coding has been used: the dark (blue) corresponds to small values of R meaning thus a desynchronised state, while the light (yellow) is associated with large value for the average order parameter, i.e. synchronisation. Taking into account that the all-to-all network of Fig. 2 (main text) and panel d of Fig 2SM is the limit of a Newman-Watts network with parameter p = 1, our results confirm the intuition that the smaller the M the larger should be the strength to achieve a desynchronised state; however we can also emphasise the fact that diluted networks are more easily desynchronised than denser ones, the dark (blue) zone in panel a is much larger than in the remaining panels. Finally, we observe that for small p the system almost never synchronises (no light yellow spots) for all the used values of (M, γ/4); indeed in the case of a regular lattice corresponding to p = 0 the system never synchronises (data not shown).

COMPARISON WITH THE PROPORTIONAL-DIFFERENTIAL FEEDBACK (PDF) CONTROL TECHNIQUE
The aim of this section is to compare the control method proposed here (for unweighted and complete networks) with the proportional-differential feedback (PDF) control method [9], whose importance has been already proved in the literature.
Let us start by briefly recalling the PDF method and fix, for the sake of completeness, the Kuramoto model on binary complete networks as working problem. In this way we can easily compare the method with the one presented in the main text. The method relies on the hypothesis of dividing the population of oscillators into two groups, the first one made by N 1 oscillators whose mean-field signal will be measured during the time evolution and then used to directly control the second group of N 2 oscillators. This translates into the following equations for the oscillators phases φ j :φ where Θ(·) is the Heaviside function, Θ(k) = 0 if k ≤ 0 and Θ(k) = 1 if k ≥ 1, its role being thus to limit the feedback action only to nodes whose index satisfies j ≥ N 1 + 1, i.e. nodes of the second group. The control signal is composed by two terms, the leftmost one being proportional to the mean-field signal of the first group of oscillators, being P the proportionality constant. The rightmost term contains the time derivative of the mean-field signal of the first N 1 oscillators, namely it is the average of the derivative of the phases of the first group modulated with a signal retarder by a fourth of the period, i.e. involving the cosine function instead of the sine one. D is a parameter defining the strength of this second term. Let us observe that this empirically based functional dependence, (cos), is consistent with the control strategy we derived based on the Hamiltonian control theory, which thus can support the goodness of this choice.
Having fixed the strength of the non-linear term in the Kuramoto model, K, (and of course the number N of oscillators) the PD controlled systems (11) depends on three free parameters: N 2 the number of controlled oscillators, P the strength of the proportional term and D the strength of the differential term. Let us observe a main difference with respect to our scheme; if one uses few controllers in the PDF strategy, i.e. N 2 is small, then automatically one assumes to be able to acquire the signal from a large number of neurons, namely N 1 = N − N 2 is large, this reflects the distinction between measuring nodes and controlled nodes and can be an issue in real applications. This division does not exist in our setting, where the same microelectrodes are used to acquire the signal and then to re-inject the controlled signal; in this way we can use a small number of controllers without the need to measure the signal from a large portion of the system. In other words, our model has more local flavour.
One can show [9] that the PD control algorithm involving both the proportional and differential feedback components is more efficient than the case D = 0, moreover one can compute [9] the desynchronisation thresholds,P andD, to achieve a desynchronised regime and prove that the latter increase if the number of oscillators in the second group, N 2 , decreases. One needs larger parameters values for P and D to destabilise the synchronous state once N 2 is small. This behaviour is similar to the one of our model.

COMPARE THE TWO CONTROL STRATEGIES.
The goal of this section is to compare the PDF control method and the one we proposed in the main text based on the Hamiltonian control theory. Let us recall that the PD-control is based on a division of the oscillators ensemble into two groups, the first one made by N 1 oscillators whose signal will be measured and the second one composed by N 2 = N − N − 1 oscillators upon which the control will directly act. As already observed our model select M oscillators among N whose signal is measured and the same M oscillators are used to control the system. So to compare the two methods we can assume that N 2 = M , this value being a proxy for the invasiveness of the method. A second proxy for the invasiveness, is the strength of the injected signal; in our scheme this is given by γ/4 (see Eq. (12) in the main text) while in the PD control one can assume this strength given by P + D, being the control signal the sum of bounded functions whose maximal amplitude are P and D.
The results of the comparisons are reported in Fig. 3SM where we show the action of the two control strategies for several values of the strength parameters, P , D and γ/4, and number of controlled nodes, N 2 = M , for N = 100 oscillators coupled using an all-to-all setting. The Kuramoto parameter has been set to K = 0.5 and the natural frequencies drawn from the distribution N (1, 0.01), in such a way the uncontrolled system does synchronise (blue line in the panels). In the left panel we report the case of a mild control, P = D = 10, hence γ/4 = 20, using few controlled oscillators, N 2 = M = 20; we can observe that the PD-control reduces the synchronisation by a factor 5 while our method is able to completely set the system into a strongly desynchronised state. The same result is achieved by both control schemes once the strength is increased (middle panel), P = D = 20 hence γ/4 = 40. Using again the mild control strength but increasing the number of controlled nodes (N 2 = M = 50, right panel) one can achieve again a strongly desynchronised state using both strategies. Such results support our claim that the control strategy hereby proposed is minimally invasive both in terms of number of microelectrodes needed and strength of the signal used to achieve a desynchronisation state. To gain a more global view of the PDF control scheme an compare it with our method, we perform a study of the impact of the parameters N 2 , P and D on the desynchronisation issue for the Kuramoto model (on a complete unweighted network). Using several combinations of the involved parameters, we realised independent simulations of the PDF control system and we measure the averaged (over the realisations) order parameter R . The results are reported in Fig. 4SM, where we varied P and N 2 , while D is fixed and equal to D = 0 (left panel) and to D = 10 (right panel). One observe that the latter control strategy gives the best results in terms of desynchronisation, the zone associated with a small R (green-blue) being larger. Such results should be compared with the ones reported in Fig. 2 (main text) or panel (d) of Fig. 2, where a similar analysis has been done using our control scheme. In particular the parameters used in the left panel of Fig. 4 corresponds to panel (d) of Fig. 2, γ/4 = P + D = P ranges in [2,10] (being D = 0) and N 2 = M range in [5,50]. While the behaviour is globally the same, the use of the same colour code easily allows to appreciate the fact that in panel (d) of Fig. 2 the desynchronised state is stronger, namely R smaller (darker blue). Moreover (see right panel), at fixed N 2 = M , the PD-control requires a larger control strength (observing the vertical scale and the choice D = 10, we conclude that the strength varies in [12, 50]) with respect to our method, to achieve an equivalent level of desynchronisation.