Systems and methods for generating a stimulation signal for brain or nerve stimulation and for performing brain or nerve stimulation

Abstract

Systems and methods for generating stimulation signals for brain or nerve stimulation are disclosed. In one arrangement, a signal generation unit generates a stimulation signal and transmits the stimulation signal to a stimulation unit. The stimulation unit applies the stimulation to a biological system. The stimulation signal is defined by at least one stimulation parameter. A data receiving unit receives time series biomarker data comprising measurement samples representing one or more biomarkers indicative of a condition of the biological system that is affected by applied stimulation based on the stimulation signal. A controller controls generation of the stimulation signal based on the received biomarker data. The controller uses a model of the biological system to estimate an optimal value of the at least one stimulation parameter using the biomarker data. The model defines strengths of contribution of measurement samples to the estimate based on the times at which the measurement samples were obtained.

Description

[0001] SYSTEMS AND METHODS FOR GENERATING A STIMULATION SIGNAL FOR BRAIN OR NERVE STIMULATION AND FOR PERFORMING BRAIN OR

[0002] NERVE STIMULATION

[0003] The present disclosure relates to generating stimulation signals for brain or nerve stimulation.

[0004] Over the last four decades neuromodulation has been adopted as an adjunct therapy for a variety of treatment resistant neurological and psychiatric disorders, whereby electrical or magnetic stimulation of certain brain regions is used to reduce disease symptoms. Invasive and non-invasive stimulation-based therapies (e.g., Deep Brain Stimulation (DBS), Transcranial Magnetic Stimulation (TMS), etc.) have so far relied on clinical expertise to empirically select patient specific stimulation parameters during clinical appointments. Due to the potentially very high dimensionality of stimulation parameterization, manual identification of effective therapy parameters is a difficult and time-consuming process. Electrode contacts and stimulation parameters (e.g., amplitude and frequency) must be tuned to provide optimal therapeutic outcomes. Therapy performance is guided by both symptom severity and side-effects, where an ideal parameter set is chosen as one that minimizes symptoms whilst also avoiding stimulation induced side-effects. In recent years, there also has been an increasing acknowledgement of the limitations of titrating therapies based on the limited snapshot of disease characteristics available from daytime clinical assessments.

[0005] Optimization and machine-learning techniques have received considerable interest for more rapid and efficient identification of patient-specific stimulation parameters. These techniques often focus on the sequential evaluation of samples from the stimulation parameter space to characterize and guide future exploration of the parameter space or to identify an optimal parameter set for therapy. Implementation of these approaches in practice requires objective quantification of each sample selected from the parameter space in terms of its effect on patient’s symptom severity. Biomarkers (i.e., signals that correlate with specific states of health or disease) can be derived from data recorded invasively or non-invasively and then utilized to assess the relationship between therapy settings and performance. A variety of such signals have been explored in the literature for a range of neurological and psychiatric disorders, but selection of optimal therapy settings remains challenging.

[0006] It is an object of the present disclosure to improve generation of stimulation signals, by supporting selection of stimulation parameters that provide improved reduction of patient symptoms and / or side-effects induced due to stimulation.

[0007] In neurodegenerative disorders such as Parkinson’s disease, therapy efficacy is lost over time due to disease progression. To counter act this, either the medication dose is increased or for those patients treated with deep brain stimulation the stimulation amplitude is increased to capture a larger portion of the target nucleus. Other variations which influence stimulation efficacy include concurrent medication intake, which in a subset of Parkinson’s disease patients may induce dyskinesias when used on its own or in conjunction with deep brain stimulation. Beyond disease progression and medication intake, symptom variations over the course of the day, which may reflect circadian dependencies (e.g., Parkinson’s disease, Essential Tremor, etc), could also influence therapy efficacy.

[0008] Brain stimulation-based therapies currently take a time-invariant approach to stimulation parametrisation. Optimal stimulation settings are chosen during day-time clinical visits through a process of trial and error and these settings remain the same until the next clinic visit. Adaptive brain stimulation approaches which switch stimulation ON when symptom severity or a representative neural biomarker exceeds a threshold, also assume that the subject specific thresholds do not change over time. Other adaptive stimulation approaches, such as phase-specific stimulation, also do not consider dynamic variations in effective stimulation settings and use the same stimulation phase. However, taking a time invariant approach to stimulation parametrisation may either over stimulate or under stimulate patients, resulting in stimulation induced side-effects and sub-optimal symptom management, respectively.

[0009] Although a variety of optimization approaches have been explored for neuromodulation therapy in both clinical and computational settings, a limitation of the approaches suggested to date is that the therapy parameter spaces are assumed to be timeinvariant (i.e., they are fixed). In practice, patient’s symptom severity can vary across multiple timescales due to factors such as changes over several hours arising due to medication intake or the sleep / wake circadian rhythm, or fluctuations at a much longer timescale due to disease progression. For this reason, optimization approaches which focus on the identification of a time-invariant, or stationary, optimal parameter set may provide suboptimal therapy as patient symptoms, or the mapping between stimulation settings and performance, drifts in response to these temporal variations over time.

[0010] Arrangements of the present disclosure provide improved performance by addressing these limitations. A time varying controller is provided for the delivery of brain stimulation that can maintain optimal stimulation performance in the face of shifting physiological demands that might be expected from the factors introduced above.

[0011] According to an aspect of the invention, there is provided a system for generating a stimulation signal for brain or nerve stimulation, comprising: a signal generation unit configured to generate a stimulation signal for performing brain or nerve stimulation on a biological system and to transmit the stimulation signal to a stimulation unit configured to apply the stimulation to the biological system, wherein the stimulation signal is defined by at least one stimulation parameter; a data receiving unit configured to receive time series biomarker data comprising measurement samples obtained by measurements performed on the biological system at different respective times, the biomarker data representing one or more biomarkers indicative of a condition of the biological system that is affected by brain or nerve stimulation applied by the stimulation unit based on the stimulation signal; and a controller configured to control generation of the stimulation signal by the signal generation unit based on the received biomarker data, wherein: the controller is configured to use a model of the biological system to estimate an optimal value of the at least one stimulation parameter using the biomarker data, wherein the model is configured to define strengths of contribution of measurement samples to the estimate based on the times at which the measurement samples were obtained.

[0012] Configuring the model to define the strengths of contribution of samples to the estimate based on the times when the samples were obtained has been shown to improve how stimulation is parametrised at a later time, thereby improving the process of generating stimulation signals. Stimulation signals can be generated that take account of variations in the state of the biological system over time, for example allowing a medical treatment performed by the stimulation signals to remain effective over time. An approach to therapy optimisation is thus enabled that can consider variations in therapy efficacy over time, due to factors such as medication intake, circadian rhythms and disease progression. By learning these variations, the time-varying controller can sustain therapy efficacy in contrast to time-invariant approaches which show degraded performance.

[0013] In an embodiment, the model comprises a forgetting component configured to progressively decrease strengths of contribution of measurement samples as a function of increasing age of the measurement samples. Use of such a forgetting component is demonstrated to provide effective tracking of changes in the biological system over time, particularly where a strength of the forgetting component (represented for example by a forgetting factor e) is selected carefully to provide an appropriate balance between considering sufficient samples to gather information and lowering the strength of contribution of older samples.

[0014] In an embodiment, the model further comprises a structured temporal component comprising one or more periodic components. The inventors have found that combining a forgetting component and one or more periodic components provides particularly robust tracking of changes in the biological system over time.

[0015] As described above, while use of a forgetting component alone can provide good performance when the strength of the forgetting component is carefully selected, combining the forgetting component with a periodic component increases the range of strengths of the forgetting component for which good performance is achieved, thereby increasing ease of setup, robustness, and performance stability. Requirements for accurate selection of the strength of the forgetting component are relaxed. Good performance can be achieved, for example, even when the strength of the forgetting component is too low for the forgetting component to provide good performance when used on its own (see Figures 13 and 14 and accompanying discussion below).

[0016] Furthermore, as demonstrated and described below with reference to Figure 14, the inventors have additionally found that combining the forgetting component and the one or more periodic components provides performance that is more robust to variations in the period (or periods) of the one or more periodic components.

[0017] In an embodiment, the structured temporal component further comprises an aperiodic component. The aperiodic component may comprise a progressive drift. The progressive drift may correspond to progression of a disease state of the biological system. Providing such an aperiodic component may improve the model’s performance in tracking aperiodic changes in the state of the biological system (which may be represented by aperiodic changes in the shape of the objective function in cases where the model comprises a Gaussian process).

[0018] Embodiments of the disclosure will now be further described, merely by way of example, with reference to the accompanying drawings.

[0019] Figure 1 depicts a system for generating a stimulation signal and a stimulation unit for applying stimulation.

[0020] Figure 2 depicts an example framework for methods of generating a stimulation signal and performing stimulation using the signal.

[0021] Figures 3 and 4 depicts parameter space exploration (Figure 3) in comparison to exploitation (Figure 4) for a sequential selection of samples during Bayesian Optimization. The top panels illustrate the true (dotted black line) and estimated (solid line) shape of the underlying objective function. The confidence bounds for the estimated shape of the objective function are highlighted in grey while the associated samples used for the estimation are represented as circles. Selection of the next parameter value to be tested is subsequently determined as the minimum of an acquisition function calculated from the mean and confidence bounds of the current estimated shape of the objective function. The bottom panels illustrate a corresponding lower confidence bound acquisition function, where the value k determines whether the algorithm prioritizes exploring regions of the parameter space where there is greater uncertainty (exploration) or prioritizes selecting regions of the parameter space where there is less uncertainty (exploitation). The location of the minimum value in the acquisition function (circle) is subsequently selected as the next parameter to be tested by the Bayesian Optimization.

[0022] Figures 5 to 7 depict time-varying objective functions and their respective covariance functions. Various explored drifts are illustrated for a periodic objective function with respect to a parameter of interest (phase trigger). The objective function is allowed to drift according to three cases: gradual (Figure 5(a), periodic (Figure 5(b)), and the superposition of the two (Figure 5(c)). At each time-step in the optimisation, samples from the objective function are weighed by their respective temporal covariance function (depicted in Figures 6(a)-(c)). This resultant weighing is illustrated in Figures 7(a)-(c) in the “colour” (represented in greyscale) gradient of the drifting objective function, with brighter colours (lighter greys) representing a stronger sample weight.

[0023] Figure 8 depicts a Kuramoto deep brain stimulation (DBS) model response to phase-locked stimulation. Individual oscillators in the population are represented as small solid circles distributed around the unit circle where the position of each oscillator represents its phase value at the associated timestep. The population order parameter at different simulation timesteps is illustrated in Figure 8(b). In the top row of Figure 8(b) the radial line within the unit circle illustrates the magnitude of the order parameter, p, and the mean phase coherence of the population, xp. The radial line points in the direction of xp and the length of the line represents the synchronization level of the population p (a short line represents low synchrony, and a long line represents high synchrony). Figure 8(a) represents the distribution of oscillators prior to stimulation when the model is at its steady state behaviour. In Figure 8(b) the rows illustrate the change in population synchrony in response to stimulation at a specific phase value, xptarget. Figure 8(c) summarizes the change in population synchrony, d , as a function of phase-locked stimulation applied between xptarget= (—n, n) radians, d < 0 corresponded to phase-locked stimulation at the specified xptargetvalues desynchronizing the population, while dp > 0 corresponded to the phase-locked stimulation increasing population synchrony.

[0024] Figure 9 is an illustration of dynamic variation of a population ARC and PRC. The shape of the initial population ARC corresponded to a PRC o at the first optimization step (t = 1). By incrementally adding an offset value, d0, to the PRF, the shape of the ARC and location of the optimum phase value for phase-locked stimulation, xp*, were gradually varied over the optimization process.

[0025] Figures 10(a)-(e) depict TV-BayesOpt algorithm performance for tracking a gradual drift in the optimal stimulation phase for phase-locked stimulation, xp*. Figure 10(a) illustrates the tracking performance of the TV-BayesOpt algorithm (upper series of dots) at locating the true optimal phase value (lower series of dots) for population desynchronization. Figure 10(b) illustrates the associated average regret for the TV- BayesOpt algorithm (upper line) at tracking the true optimum phase value in comparison to when static BayesOpt was implemented alone (lower line). Figures 10(c)-(e) illustrate the true location of the optimal stimulation phase (black dot), the minimum predicted by the TV-BayesOpt algorithm (black circle and cross) and the shape of the GPR predicted by the TV-BayesOpt algorithm at 0.01-, 17.62- and 20.00-day time points respectively.

[0026] Figures 1 l(a)-(e) depict TV-BayesOpt algorithm performance for tracking a periodic drift in the optimal stimulation phase for phase-locked stimulation, xp*. Figure 11(a) illustrates the tracking performance of the TV-BayesOpt algorithm (dark grey dots) at locating the true optimal phase value (black dots) for population desynchronization. The dark grey dots and the black dots overlie each other closely. Figure 11(b) illustrates the associated average regret for the TV-BayesOpt algorithm (upper line) at tracking the true optimum phase value in comparison to when static BayesOpt was implemented alone (lower line). Figures 1 l(c)-(e) illustrate the true location of the optimal stimulation phase (black dot), the minimum predicted by the TV-BayesOpt algorithm (black circle and cross) and the shape of the GPR predicted by the TV-BayesOpt algorithm at 0.01-, 17.62- and 20.00-day time points respectively.

[0027] Figures 12(a)-(e) depict TV-BayesOpt algorithm performance for tracking a gradual and periodic drift in the optimal stimulation phase for phase-locked stimulation, xp*. Figure 12(a) illustrates the tracking performance of the TV-BayesOpt algorithm (dark grey dots) at locating the true optimal phase value (black dots) for population desynchronization. The dark grey dots and the black dots overlie each other closely. Figure 12(b) illustrates the associated average regret for the TV-BayesOpt algorithm (upper line) at tracking the true optimum phase value in comparison to when static BayesOpt was implemented alone (lower line). Figures 12(c)-(e) illustrate the true location of the optimal stimulation phase (black dot), the minimum predicted by the TV-BayesOpt algorithm (black circle and cross) and the shape of the GPR predicted by the TV-BayesOpt algorithm at 0.01-, 17.62- and 20.00-day time points respectively.

[0028] Figure 13 depicts a sensitivity analysis of the TV-BayesOpt algorithm with a forgetting (upper line) or a forgetting-periodic (lower line) covariance function for a range of e values. Incorporation of prior knowledge of the temporal variation in the objective function optimum value (lower line) resulted in improved TV-BayesOpt algorithm performance than implementing a forgetting covariance function (upper line) alone. Best algorithm performance was observed for an e value of 0.22; above this value the algorithm begins to forget previous samples too quickly to accurately track the optimal value in the objective function.

[0029] Figure 14 depicts a sensitivity analysis of the TV-BayesOpt algorithm performance at tracking a periodic temporal drift when the period of the drift is offset from the covariance function period anticipated by the TV-BayesOpt algorithm. The performance of a scheduler (uppermost line), the TV-BayesOpt with a periodic temporal covariance function (middle line line) and the TV-BayesOpt with a periodic and smooth forgetting covariance algorithm (lowermost line) were estimated by calculating the AUC value for the associated regret plot for each implementation at different offset in the true temporal period.

[0030] Methods of the present disclosure are computer-implemented. Each step of the disclosed methods may therefore be performed by a computer. The computer may comprise various combinations of computer hardware, including for example CPUs, RAM, SSDs, motherboards, network connections, firmware, software, and / or other elements known in the art that allow the computer hardware to perform the required computing operations. The required computing operations may be defined by one or more computer programs. The one or more computer programs may be provided in the form of media, optionally non- transitory media, storing computer readable instructions. When the computer readable instructions are read by the computer, the computer performs the required method steps. The computer may consist of a self-contained unit, such as a general-purpose desktop computer, laptop, tablet, mobile telephone, smart device, etc. Alternatively, the computer may consist of a distributed computing system having plural different computers connected to each other via a network such as the internet or an intranet.

[0031] As exemplified schematically in Figure 1, the present disclosure relates to a system 2 for generating a stimulation signal for brain or nerve stimulation. The system 2 may perform all or part of a method according to the framework depicted in Figure 2. The method comprises receiving biomarker data from a biomarker data source in step S 1 , generating a stimulation signal in step S2, and applying the stimulation signal in step S3. Measurements made during the applying of the stimulation signal in step S3 may lead to new data being added to the biomarker data in the biomarker data source.

[0032] The system 2 comprises a signal generation unit 4. The signal generation unit 4 generates a stimulation signal 5 (step S2) for performing brain or nerve stimulation on a biological system 12 (step S3). The signal generation unit 4 transmits the stimulation signal 5 to a stimulation unit 10. The stimulation unit 10 applies the stimulation to the biological system 12. In some arrangements, the simulation unit 10 comprises one or more electrodes 14 and applies the stimulation to the biological system 12 via the one or more electrodes 14. The stimulation unit 10 may comprise a power source for applying the stimulation based on the stimulation signal. The stimulation may be achieved, for example, by driving the one or more electrodes 14 with an electrical signal provided by the power source. The biological system 12 may comprise the brain or a portion of the nervous system of a human or animal other than the brain.

[0033] The stimulation signal is defined by at least one stimulation parameter. In some arrangements, as exemplified in the detailed example described below, the stimulation signal is configured to be phase-locked to a biomarker of the biological system and the at least one stimulation parameter comprises a stimulation phase relative to the biomarker. The biomarker may, for example, comprise a disease biomarker having rhythmicity, which in the case of Parkinson’s disease may range from 14-30 Hz, and for tremor may range from 3-8 Hz.

[0034] Alternatively, or additionally, the stimulation parameter may comprise one or more of the following: a stimulation amplitude; a stimulation pulse width; and a stimulation pulse frequency.

[0035] The system 2 further comprises a data receiving unit 6. The data receiving unit 6 receives time series biomarker data 7 (step SI). The biomarker data 7 comprises measurement samples (which may be referred to herein as samples) obtained by measurements performed on the biological system 12 at different respective times. The biomarker data 7 is thus derived from measurements performed on the biological system 12. The measurements may be performed using the one or more electrodes 14 and / or by other means. The biomarker data represents one or more biomarkers indicative of a condition of the biological system that is affected by brain or nerve stimulation applied by the stimulation unit 10 based on the stimulation signal 5. The one or more biomarkers may comprise a disease biomarker having rhythmicity, such as a biomarker affected by Parkinson’s disease and / or tremor. The biomarker data may indicate a degree to which the stimulation applied by the stimulation unit 10 is effective in suppressing a disease biomarker. The condition of the biological system 12 that is affected by the stimulation may thus be a disease condition such as Parkinson’s and / or tremor.

[0036] The system further comprises a controller 8. The controller 8 controls generation of the stimulation signal by the signal generation unit 4 based on the received biomarker data 5. The controller 8 uses a model of the biological system 12 to estimate an optimal value of the at least one stimulation parameter using the biomarker data 5. The model is configured to define strengths of contribution (which may be referred to herein as weightings) of measurement samples to the estimate based on the times at which the measurement samples were obtained.

[0037] In some arrangements, as exemplified in the detailed example described below, the model comprises a Gaussian process. The Gaussian process may be configured to estimate an objective function representing an expected variation of the biomarker data as a function of the at least one stimulation parameter of the stimulation signal. The objective function may for example be considered as a mapping between stimulation settings and symptom suppression. By estimating the objective function, the Gaussian process can thus enable estimation of a value of the at least one stimulation parameter that optimally suppresses a disease biomarker. In the detailed example described below, an amplitude response curve (ARC) representing a change in a population synchrony (p) (of a population of neurons) due to stimulation delivered at a particular phase (\| / ) is used as an objective function. By configuring the model to define strengths of contribution of measurement samples to the estimate based on the times at which the measurement samples were obtained, the model can take account of changes in the objective function over time.

[0038] The Gaussian process may be defined by a covariance function. As exemplified in the detailed example described below, the covariance function may comprise a temporal part. The covariance function may additionally comprise a time-invariant spatial part representing prior knowledge of a shape of the objective function. In some arrangements, the temporal part and spatial part are combined into a spatio-temporal covariance function using the Hadamard product, as described below with reference to equation (9).

[0039] In some arrangements, the spatial part is implemented as described below with reference to equation (10).

[0040] The model (e.g., Gaussian process) may comprise a forgetting component. In cases where the model comprises a Gaussian process, the forgetting component may be represented by the temporal part of the covariance function. The forgetting component is configured to progressively decrease strengths of contribution of measurement samples as a function of increasing age of the measurement samples. A strength of the forgetting component (e.g., how quickly the forgetting component causes contributions from older measurement samples to be suppressed) may be represented by a forgetting factor (referred to as e in the detailed example below). The forgetting factor may, for example, have a value between 0 and 1, where 0 represents no samples being forgotten and values near to 1 correspond to samples being forgotten so quickly that only a current sample is used to estimate the objective function.

[0041] In some arrangements, the forgetting component is implemented as described below with reference to equations (11) and (12). In some arrangements, the forgetting factor may be learned using past data (e.g., past samples).

[0042] The model (e.g., Gaussian process) may comprise a structured temporal component. In cases where the model comprises a Gaussian process, the structured temporal component may be represented by the temporal part of the covariance function. As exemplified in the detailed example described below, the structured temporal component may comprise one or more periodic components. Each periodic component is configured to contribute a variation in the strength of contribution of measurement samples to the estimate that is periodic as a function of time. Each periodic component may thus promote a greater degree of similarity between strengths of contributions of measurement samples taken at corresponding time points in different cycles. For example, if a periodic component has a period of one day, the periodic component will contribute to samples taken at the same time on different days having more similar weightings than samples taken at different times of day. Furthermore, when the model is estimating an optimal value of the at last one stimulation parameter at a particular time point in the cycle of a periodic component, the periodic component will promote increased strengths of contribution to the estimate from samples taken at the same time point in previous cycles. This is illustrated for example in Figures 6(b) and 7(b), discussed below. In some arrangements, the periodic component is implemented using an exponential sine squared covariance function as exemplified in the detailed example below.

[0043] In some arrangements, the periodic component is implemented as described below with reference to equation (13).

[0044] In some arrangements, the combination of the forgetting component and the one or more periodic component is implemented as described below with reference to equation (14).

[0045] The one or more periodic components may arise for example due to factors such as biological rhythms, such as the circadian sleep-wake cycle or medication cycles. The one or more periodic components may have a predetermined period (e.g., one day, or a period between administration of treatments to a patient. Alternatively, or additionally, the controller 8 may be configured to leam a period of one or more of the periodic components. The structured temporal component may additionally or alternatively comprise an aperiodic component. As exemplified in the detailed example described below, the aperiodic component may comprise a progressive (e.g., monotonic) drift. The progressive drift may be caused, for example, by progression of a disease state of the biological system 12.

[0046] Embodiments of the disclosure may be implemented as methods. The methods may include computer-implemented methods of generating a stimulation signal for brain or nerve stimulation. The methods may comprise a method of performing brain or nerve stimulation. Such methods may include generating a simulation signal and applying stimulation to the biological system based on the generated stimulation signal. The methods may include steps within the example framework method depicted in Figure 2. A computer program and / or a computer program product (e.g., a non-transitory computer program product) may be provided for executing computer-implemented method steps.

[0047] For example, the method may comprise receiving time series biomarker data (e.g., step SI in Figure 2). The biomarker data comprises measurement samples obtained by measurements performed on a biological system at different respective times. The biomarker data represents one or more biomarkers indicative of a condition of the biological system that is affected by brain or nerve stimulation applied by a stimulation unit based on a stimulation signal. The stimulation signal is defined by at least one stimulation parameter. The method may comprise generating a stimulation signal (e.g., step S2 in Figure 2) for the stimulation unit based on the received biomarker data by using a model of the biological system to estimate an optimal value of the at least one stimulation parameter using the biomarker data. The model is configured to define strengths of contribution of measurement samples to the estimate based on the times at which the measurement samples were obtained. The generating of the stimulation signal may comprise generating a stimulation signal for which the at least one stimulation parameter has the estimated optimal value.

[0048] Detailed Example

[0049] Overview

[0050] Reference is made below to an algorithm based upon Bayesian optimization (referred to as “BayesOpt”) with modifications to enable time-adaptive therapy optimization (Time-Varying Bayesian Optimisation, referred to as “TV- BayesOpt”). The TV-BayesOpt algorithm is an example of how the controller 8 may implement the estimation of the optimal value of the at least one stimulation parameter using biomarker data, as described above.

[0051] Further details are given below about implementation choices for TV- BayesOpt, namely the surrogate model (i.e., a model used to estimate the shape of the true objective function) and the acquisition function (i.e., how to select the next sample from the parameter space). The surrogate model is an example of the model of the biological system discussed above. The acquisition function is an example of how the controller may use the biomarker data, for example to update the model and / or exploit the model to provide an estimation of an optimal value of the at least one stimulation parameter.

[0052] Spatiotemporal covariance functions used for a time-varying implementation of the BayesOpt algorithm (TV-BayesOpt) are described that incorporate prior knowledge regarding the form of the objective function along with its anticipated drift in time. These covariance functions are examples of the covariance functions defining the Gaussian process implementation of the model of the biological system discussed above. A form is established of a dynamical model of neural synchrony and is used to evaluate the performance of the system 2 in delivering optimal stimulation to suppress a rhythmic biomarker as a proxy for symptom suppression.

[0053] Kuramoto Model

[0054] The Kuramoto model was used to investigate the performance of the TV-BayesOpt algorithm. The Kuramoto Model provides a simple mathematical model to describe synchronization in biological systems. In this specific case the model represents a population of neurons as N coupled oscillators that oscillate with an intrinsic frequency ^biomarker- This intrinsic frequency reflects the rhythmicity of the disease biomarker, which in the case of Parkinson’s disease ranges from 14-30 Hz, and for tremor ranges from 3-8 Hz. The strength of interaction between the oscillators is captured by a coupling term Kbiomarker where the influence of the jthoscillator on the ithoscillator is described in terms of phase progression (i.e. speeding up or slowing down) of the ithoscillator depending upon their difference in phase.

[0055] The model has been extended to account for the effects of exogenous perturbations which is given by: where I is the intensity of stimulation, X(t) is a function whose value equals 1 at time t if stimulation is applied and is 0 otherwise, and Z(0j) is the phase response curve (PRC), which describes the response of the ithoscillator to stimulation (i.e. speeding up or slowing down) depending on the precise timing of the stimulation with respect to the oscillator’s phase, 0,. Neurons have been classified into two distinct categories based on their PRCs that characterise their response to perturbation: type I PRCs either exclusively delay or advance spike firing, whereas type II PRCs can both advance or delay dependent upon the specific phase at which stimulation is delivered. It has been demonstrated that neural oscillators in tremor exhibit type II PRCs. Therefore, in the present example, the PRC is described as Z Ot) = — sin(0j).

[0056] Under certain assumptions, the Kuramoto model can be solved to express the behaviour of each oscillator in terms of the order parameters, p and xp, describing the population’s mean phase-coherence and mean phase respectively.

[0057] The population’s mean phase-coherence, p, represents synchrony on the interval [0, 1] where p = 0 and p = 1 correspond to either complete desynchrony or synchrony, respectively. In this manner, p can be used as a control signal used a proxy for population level synchrony similar to that currently used in adaptive deep brain stimulation (DBS) to capture instantaneous biomarker power, or level of synchronization, in specific frequency bands of a neural oscillation. Following from the introduction of oscillopathies, it is assumed that lower levels of p corresponds to reduced disease symptoms whereas higher levels of p reflects increased symptom severity. The performance of the TV -BayesOpt algorithm can then be assessed by quantifying how the precise stimulation timing, i.e., the optimal xp at which to apply stimulation, can minimize population synchrony (in terms of P)-

[0058] The objective function for the TV-BayesOpt algorithm is represented as the effect of phase-locked stimulation on synchrony (p). This can be quantified with an amplitude response curve (ARC) where the ARC represents the change in the population synchrony (p) due to stimulation delivered at a particular phase (xp). The ARC therefore represents the underlying objective function (mapping between stimulation settings and symptom suppression) where certain phases of stimulation may increase or decrease population synchrony. The shape of population ARC can be varied over time by updating the PRC (Z(0j)) for the oscillators to emulate gradual or periodic changes.

[0059] Bayesian Optimization

[0060] To estimate the form of the ARC that maps from stimulation parameters to population synchrony (i.e., the objective function), a surrogate model based on a Gaussian Process (GP) prior was implemented. The GP prior model can be thought of as a probability distribution over a set of non-linear regression functions (as opposed to a Gaussian distribution, which is a set over variables). This permits efficient estimation of a mean function, as well as its relative uncertainty across the parameter space, indicating the confidence in the evaluated objective function for each parameter value.

[0061] The objective function / (x) is fully defined by its mean ju(x) and covariance K (x, x ' ) , that form sufficient statistics such that f (x) ~ GP p X) , K (x, x ' ) ) , where x in this context are the stimulation parameters to be explored. The covariance function describes the relationship between points in the parameter space, and is vital in determining the form (e.g., smoothness, or periodicity) of the estimated objective function. Thus K(x, x') describes the spatial covariance (i.e., the space between parameter x and its neighbour x'). An example covariance function for use with arrangements of the present disclosure is described below.

[0062] The mean and covariance functions allow us to derive a joint distribution of the noisy outputs y (previously sampled simulation parameters) and the estimated value of the objective function (mapping between stimulation settings and symptom suppression) / * evaluated at a new sample point x* in the parameter space. This may be defined as:

[0063] (4) additional samples from the parameter space are evaluated, the prior is updated to form a posterior distribution to improve the model’s approximation of the shape of the objective function, (x). Conjugacy between the prior and likelihood allows the posterior to be computed analytically by omitting the evidence term. The predictive distribution can therefore be written as p(f^Dn, x ~N(ji*,aY) (5)

[0064] Where Dn= {x, f represents the previously sampled parameter values and their respective outcome values and the posterior distribution is defined as The acquisition function is used to identify regions of the parameter space which should be tested during the optimization process. To determine the next sample to test during the optimization process, the acquisition function cr(x) uses the current estimate of the underlying surrogate model and calculates the expected utility of all samples in the parameter space. The next sample to be tested is subsequently selected as the sample point x* which minimizes the acquisition function: x4= argmin a (x; Dn)

[0065] A Gaussian Process - Lower Confidence Bound (GP-LCB) acquisition function was chosen for the BayesOpt algorithm as it allows flexible tuning of the amount of parameter space exploration during the optimization process. The GP-LCB is defined as where the parameter Knrescales tft- the estimated uncertainty of the GP. To balance between exploration of under sampled regions of the parameter space, or exploit already sampled regions, the value of Knin the acquisition function is chosen according to the specific needs of the user where high Knvalues favour exploration, and low Knvalues favour exploitation, of the parameter space. The role of Knfor guiding this explorationexploitation trade-off when sampling from the parameter space between optimization steps is illustrated in Figure 3 (exploration) and Figure 4 (exploitation) in which examples are shown of either maximally explorative or exploitative acquisition functions.

[0066] To track a dynamically varying set of optimal therapy parameters, arrangements of the present disclosure eschewed the common assumption that the objective function is time invariant. When applied to time-varying systems the standard BayesOpt routine will exhibit diminished performance as it assumes all samples have equal contribution to the current estimate of the objective function. A TV-BayesOpt algorithm in contrast weighs the contribution of previously tested samples to the current estimate based on the time at which samples were taken in relation to the current estimation step. Sources of temporal variation which lead to changes in the optimal parameter(s) may be gradual (e.g., linear), such as during disease progression, or due to other factors such as biological rhythms, such as the circadian sleep-wake cycle or medication cycles (e.g., periodic). The goal of the TV- BayesOpt algorithm is thus to find the optimal parameter set for a given time, i.e., the parameters which provide the best suppression of symptoms, and to subsequently track the development of this parameter set over time.

[0067] To meet these requirements, a spatio-temporal covariance function was constructed. The spatio-temporal covariance function comprises two parts which are referred to as the spatial covariance function (or spatial part) and the temporal covariance function (or temporal part), respectively. The spatial covariance function in this context is equivalent to the covariance function used for time- invariant BayesOpt and captures the prior knowledge of the shape of the objective function. On the other hand, the temporal covariance function captures the anticipated variation in the shape of the objective function over time, i.e., rescaling the relevance of previously sampled data in relation to the most recent sample. The resulting spatio-temporal covariance function, Kt(x, x' , is formulated as the Hadamard product of the spatial, K(x, x'), and the temporal covariance functions,

[0068] A multitude of covariance functions exist across the field of machine learning that can describe complex topologies. The performance of the BayesOpt algorithm is thus dependent on the selection of an appropriately shaped covariance function for the problem that is being investigated.

[0069] In the present example, the problem investigated concerns the optimization of a target phase for phase-locked stimulation to suppress a rhythmic biomarker of disease. Therefore, the parameter space (target phase values) is periodic, where —n and n radians represent the same point in the parameter space. Therefore, the spatial covariance function (i.e., the covariance between target stimulation phases in the parameter space) was therefore selected to also be periodic using the exponential sine squared covariance function: where hyperparameters lxand Txare the length scale and period (incorporating our prior knowledge of the periodicity of the parameter space). For biological systems there are several sources which may lead to variations in the shape of the objective function over time, including disease progression (leading to monotonic changes) and biological rhythms such as the circadian sleep-wake cycle (leading to periodic changes), or a combination of the two. Examples of temporal covariance functions that overcome these temporal variations are outlined in the following sections. Topologies of these covariance functions are summarized below in Figures 5-7.

[0070] To accommodate gradual variations in the shape of the objective function a ‘forgetting’ covariance function that reduces the influence of older samples on the current estimate can be implemented, as depicted in Figures 5(a), 6(a), and 7(a). Bogunovic et al. (2016) introduced a covariance function that was defined to meet these requirements (Bogunovic I, Scarlett J, Cevher V. Time-varying Gaussian process bandit optimization. Artif Intell Stat. 2016;314-23), where the covariance function is defined as

[0071] It-t'l K(t, t') = (1 - e) 2 (H) where t is the time of the most current sample, t' is the time at which previous samples were taken and e is referred to as the forgetting factor. The forgetting factor, e, has a value between [0, 1) which determines how quickly the contribution of previously taken samples to the current estimate are reduced. A forgetting factor 6 = 0 corresponds to no data being forgotten and all previous samples contributing equally to the estimate of the objective function, while 6 values close to 1 correspond to previous data samples being forgotten so fast that only the current sample is used to estimate the objective function. Moreover, as |t — t'| -> oo the contribution of previous samples on the current estimate reduces to 0. This behaviour is graphically summarized in Figure 6(a).

[0072] An appropriate value of 6 can be learned using past data through maximum likelihood estimation or by empirically selecting 6 to reduce the contribution of previous samples based on a specified half-life:

[0073] In 2 e = A = - (12)

[0074] G / 2

[0075] Where A corresponds to the half-life decay factor and tq / 2 defines half-life as a specified number of samples, i.e., a G / 2value of 100 corresponds to the objective function value associated with the 100th last sample taken contributing only 50 % of its original value to the current estimate of the surrogate model.

[0076] To accommodate periodic variations in the shape of the objective function a periodic covariance function (i.e., a temporal covariance function part representing a structured temporal component that is periodic) can be implemented to weight the contribution of previously sampled data based on the specific time points in the cycle at which the samples were taken, as depicted in Figures 5(b), 6(b), and 7(b). In this manner, the contribution of samples taken at time points near the same parts of the cycle are weighted more similarly to each other (i.e., such that the strengths of contribution of these measurement samples to the estimate are more similar to each other) than those taken at different time points in the cycle. This behaviour can be captured in a similar manner to the spatial covariance function by using the exponential sine squared covariance function, where Ttis the temporal period of the underlying biomarker oscillation (defined with the same units as t) and lt p, the length scale, controls the amplitude of the temporal oscillation. In practice, the units of Ttwill depend on the sampling protocol and may be defined in either sample counts or physical timescales such as hours or minutes. Furthermore, Ttcan be estimated empirically by tracking the variation of symptom suppression under fixed stimulation parameters. Use of this covariance function captures temporal variations which repeat themselves exactly, as depicted in Figure 6(b).

[0077] For implementations in real-life scenarios where both gradual and periodic temporal variations may be superimposed due to both disease progression and biological rhythms, the forgetting component and the structured temporal component (e.g., periodic component) of the temporal covariance function can be combined, as depicted in Figures 5(c), 6(c), and 7(c). Doing so produces the temporal covariance function below which is the product of the forgetting covariance function (forgetting component) and the periodic covariance function (structured temporal component) where e, lt pand Ttare the covariance function hyperparameters as described above in the previous sections. This temporal covariance function is used to subsequently reduce the contribution of older samples on the current estimate while also flexibly tracking periodic variations in the shape of the objective function over time, as depicted in Figure 6(c).

[0078] Simulation Details

[0079] The behaviour of the TV-BayesOpt algorithm was investigated in the Kuramoto DBS model. The Kuramoto DBS model was initially simulated for a 200s period to allow the model to reach its steady state behaviour. After 200s, the model was subsequently simulated for 3000 optimization steps where each optimization step corresponded to a 58s simulation period. During each optimization step, phase-locked stimulation at a target phase value was applied only during the last 8s of the 58s simulation period. The initial 50s of each optimization step was included to allow the model to return to its steady state post application of phase-locked stimulation in the last optimization step. The output of each optimization step, y, was quantified as the normalized change in the mean order parameter during the stimulation period from its baseline value:

[0080] . Pstim ~ Pbaseline y = P = - - (15)

[0081] Pbaseline where pstimis the mean order parameter averaged over the 8 s stimulation period and Pbaseline is the baseline mean order parameter value and was calculated as the mean order parameter of the 25s period prior to stimulation application.

[0082] To initiate the TV-BayesOpt algorithm, data from the parameter space was needed. To this end, 12 equally spaced samples between [-n, n) were tested to coarsely characterize the population response to phase-locked stimulation at different stimulation phase values. Following this, selection of subsequent phase values to be tested at each optimization were selected by the TV-BayesOpt algorithm as described above.

[0083] Drifts in the objective function (mapping between stimulation settings and symptom suppression) shape were simulated by incrementally adding a phase offset value to the Kuramoto DBS model PRC over the steps of the optimization process. Three drifts in the optimum were simulated to explore the TV-BayesOpt algorithm performance at locating and tracking an optimal phase value for desynchronizing the Kuramoto DBS model during phase- locked stimulation. The performance of the algorithm at each timestep was quantified using the cumulative regret, defined as is the location of the true optimum value phase value for phase-locked stimulation at optimization step k, xptargetkis the stimulation phase value tested at the kthoptimization step, (*) is the measured objective function value, i.e. A , at the specified input value and n is the total number of optimization steps.

[0084] A gradual drift was simulated as a phase offset advancing from 0 to -TT over 3000 optimization steps, a periodic drift was simulated as a sinusoidal phase offset advancing from 0 to -7T and back to 0 over 100 optimization steps, and lastly the superimposed gradual and periodic drifts were simulated as the superposition of the gradual and periodic drifts.

[0085] Finally, sensitivity analyses were conducted to investigate the influence of the temporal covariance forgetting factor, e, and period, Tt, on the TV-BayesOpt algorithm’s performance at tracking a gradual and periodic temporal drift and purely periodic temporal drift, respectively.

[0086] Results

[0087] The performance of the TV-BayesOpt algorithm was investigated in the Kuramoto DBS model for locating and tracking an optimum phase value to reduce population synchrony during phase-locked stimulation. The behaviour of the Kuramoto DBS model in response to phase-locked stimulation at different population phase values was first investigated prior to simulation of gradual and periodic drifts in the location of the optimal stimulation phase. The performance of the algorithm in each scenario was investigated in comparison to the traditional time-invariant implementation of the BayesOpt algorithm.

[0088] After an initial 200s transient period the Kuramoto DBS model reached its steady state behaviour with the magnitude of the order parameter p (i.e., population synchrony which is a direct correlate of symptom severity) settling to a value around 0.8, reflecting a highly synchronized state. Once the model reached its steady state behaviour, the population’s response to phase-locked stimulation was investigated using a static PRC for the oscillator population that was set to — sin(0). Stimulation was applied to the population at 12 equally spaced mean population phase values, xptarget, between (— n, n) (Figure 8). As shown in Figure 8(c), phase-locked stimulation of the population between xptarget=(—71, — -) radians and ptarget=radians resulted in desynchronization of the population and a reduction in the magnitude of the order parameter (Ap < 0). In contrast, phase-locked stimulation applied between xptarget=radians resulted in increased population synchrony, which was reflected in an increase in the magnitude of the population order parameter (Ap > 0). The optimal phase value for phase- locked stimulation, xp*, was identified as the value that resulted in the greatest reduction in population synchrony. This value was identified as xp* = n radians. Following this, the location of xp* was dynamically varied as described in the subsequent sections.

[0089] The location of the optimal stimulation phase, xp*, achieving maximum population desynchrony, was varied by incrementally adding an offset value, A0, to the initial population PRF: + A0). This is represented in Figure 9 where the shape of the population ARC and PRC are dynamically varied by inclusion of A0. The location of xp* could thus be gradually varied over the course of the optimization process (Figure 9).

[0090] A gradual drift in the optimal stimulation phase value, xp*, needed to desynchronize the population of oscillators was tested using the TV-BayesOpt algorithm plus a temporal forgetting covariance function and a simulation in which the optimal phase drifts from TT to 0 radians (Figure 10). Such a gradual drift may be anticipated in a context such as disease progression. The TV-BayesOpt algorithm resulted in lower cumulative regret over the optimization process, calculated as the area under the curve (AUC) of the cumulative regret plot, when compared to the time-invariant implementation of the BayesOpt algorithm (Figure 10).

[0091] A periodic drift in the optimal phase value, xp*, for phase-locked stimulation was emulated by periodically varying xp* from TT to 0 radians over 100 optimization steps (representing a 1-day period). Implementation of the TV-BayesOpt algorithm with a temporal periodic covariance function (where Tt= 100) led to accurate tracking of the location of the drifting xp*, as depicted in Figure 11. Such a periodic change could be expected from regular medication intake, circadian cycles, and other repeating patterns which may influence therapy efficacy. The temporal periodic covariance function with a 1- day period led to better tracking of xp* in comparison to static BayesOpt, as evidenced by a lower cumulative regret AUC (Figure 11).

[0092] Superposition of gradual and periodic drifts in the optimal stimulation setting for phase-locked stimulation, xp*, of the Kuramoto oscillators required implementation of the temporal forgetting-periodic covariance function. Implementation of the TV-BayesOpt algorithm with this covariance function resulted in accurate tracking of the location of xp* and produced a lower cumulative regret AUC in comparison to the time-invariant BayesOpt algorithm (Figure 12).

[0093] A sensitivity analysis was conducted to investigate the influence of the temporal kernel’s forgetting factor value, e, on the performance of the TV-BayesOpt algorithm (Figure 13). For each e value investigated, the AUC of the cumulative regret plot was calculated for tracking the gradual and periodic drift illustrated in Figure 12, with high values indicating inaccurate optimization. Variation of the TV-BayesOpt algorithm performance due to prior knowledge of the temporal drift was investigated by comparing the algorithm performance when solely implementing a forgetting covariance function in comparison to a forgetting-periodic covariance function; Figure 13 (upper and lower curves respectively). For e values less than 0.04 the forgetting covariance function resulted in greater cumulative regret AUC values than the forgetting-periodic covariance function. For e values above 0.04, both temporal covariance functions resulted in equivalent average regret AUC measures, where the algorithm provided the best performance with an e value of 0.22 and corresponded to a data half-life (t^) of 3.15 samples. Incorporating prior knowledge of periodicity in temporal variation by using the temporal periodic covariance function led to more stable performance of the TV-BayesOpt algorithm than when the forgetting factor alone was implemented. Consistent cumulative regret AUC values lower than the forgetting covariance function performance alone were observed for the algorithm when incorporating the temporal periodic covariance function. e values above 0.22 resulted in increased cumulative regret AUC values as e values in this range led to the algorithm forgetting data too quickly to accurately estimate the location of the optimal value.

[0094] Thus, in summary, good algorithm performance is found to be achieved by implementation of the TV-BayesOpt algorithm with smooth forgetting alone, but this requires careful selection of an appropriate forgetting factor well-matched to the temporal oscillation being tracked, as shown in Figure 13. When too low a value of e is selected, the algorithm weighs the contribution of older samples too heavily which leads to poor tracking of the optimum stimulation phase over time, Figure 13. When the TV-BayesOpt algorithm additionally incorporates a periodic temporal covariance function whose period is aligned to the temporal drift the algorithm performance is greatly improved for low e values. Robust algorithm performance is subsequently observed for all e values below 0.4. Above this value the algorithm performance begins to worsen due to previously taken data samples being forgotten too quickly.

[0095] A sensitivity analysis was conducted to investigate the TV-BayesOpt algorithm performance when the temporal drift anticipated by the algorithm, Tt, differed from the true period of the example periodic temporal variation highlighted in Figure 11. The performance of the algorithm was also compared to that of a scheduler which made scheduled adjustments to maintain stimulation at the optimal phase value over a 24-hour period, Figure 14.

[0096] The scheduler and TV-BayesOpt with a periodic temporal covariance function were unable to track the location of the optimum stimulation phase when the temporal drift in the optimal stimulation phase was different from the anticipated period. This resulted in large AUC values for both the scheduler and algorithm when the offset to the true temporal period was not 0.

[0097] When the temporal covariance function for the TV-BayesOpt algorithm was implemented with smooth forgetting in addition to a periodic kernel the algorithm was able to accommodate the offset in the difference between the true and anticipated temporal periods, as evidenced by the lower regret AUC values.

[0098] Thus, in summary, the observations discussed above demonstrate that when there is no difference between the actual and expected temporal periods, the algorithm performs well with a periodic temporal covariance function alone and maintains tracking of the optimal stimulation phase. In this case, once the algorithm has locked to the location of the optimum stimulation phase its performance is equivalent to a scheduler that updates the location of the optimum stimulation phase value based on the anticipated temporal period (Figure 14). Worsened performance is observed by both the algorithm and scheduler however when there is a difference between the anticipated and actual temporal period (Figure 14). This worsened algorithm performance can be compensated for by use of the smooth forgetting covariance function which maintains algorithm performance at tracking the location of the optimum phase when there is a mismatch between the actual and anticipated temporal period by the TV-BayesOpt algorithm.

[0099] The above results highlight the relationship between smooth forgetting and temporal rhythmicity in the temporal covariance function. Incorporation of both smooth forgetting and knowledge of the rhythmic variations in therapy efficacy leads to better performance of the TV- BayesOpt algorithm than when either of these are utilized by the algorithm alone.

Claims

CLAIMS1. A system for generating a stimulation signal for brain or nerve stimulation, comprising: a signal generation unit configured to generate a stimulation signal for performing brain or nerve stimulation on a biological system and to transmit the stimulation signal to a stimulation unit configured to apply the stimulation to the biological system, wherein the stimulation signal is defined by at least one stimulation parameter; a data receiving unit configured to receive time series biomarker data comprising measurement samples obtained by measurements performed on the biological system at different respective times, the biomarker data representing one or more biomarkers indicative of a condition of the biological system that is affected by brain or nerve stimulation applied by the stimulation unit based on the stimulation signal; and a controller configured to control generation of the stimulation signal by the signal generation unit based on the received biomarker data, wherein: the controller is configured to use a model of the biological system to estimate an optimal value of the at least one stimulation parameter using the biomarker data, wherein the model is configured to define strengths of contribution of measurement samples to the estimate based on the times at which the measurement samples were obtained.

2. The system of claim 1, wherein the model comprises a forgetting component configured to progressively decrease strengths of contribution of measurement samples as a function of increasing age of the measurement samples.

3. The system of claim 2, wherein the model further comprises a structured temporal component.

4. The system of claim 3, wherein the structured temporal component comprises one or more periodic components.

5. The system of claim 4, wherein each periodic component is configured tocontribute a variation in the strength of contribution of measurement samples to the estimate that is periodic as a function of time.

6. The system of claim 4 or 5, wherein one or more of the periodic components has a predetermined period, optionally corresponding to a periodicity of an underlying biomarker oscillation of the biological system.

7. The system of any of claims 4 to 6, wherein the controller is configured to learn a period of one or more of the periodic components.

8. The system of any of claims 3 to 7, wherein the structured temporal component further comprises an aperiodic component.

9. The system of claim 8, wherein the aperiodic component comprises a progressive drift, the progressive drift optionally corresponding to progression of a disease state of the biological system.

10. The system of any of claims 3 to 9, wherein the model comprises a Gaussian process defined by a covariance function.

11. The system of claim 10, wherein the Gaussian process is configured to estimate an objective function representing an expected variation of the biomarker data as a function of the at least one stimulation parameter of the stimulation signal.

12. The system of claim 10 or 11, wherein the covariance function comprises a temporal part, the temporal part representing the forgetting component and the structured temporal component.

13. The system of claim 12, wherein the structured temporal component comprises a periodic component and the temporal part of the covariance function, K t, t' , comprises the following:where t is the time of a most current measurement sample, t' is a time at which a previous measurement sample was taken, e is a forgetting factor defining a strength of the forgetting component, lt pdefines a length scale of the periodic component, and Ttis a temporal period of the periodic component.

14. The system of any of claims 10 to 13, wherein the covariance function further comprises a time-invariant spatial part representing prior knowledge of a shape of the objective function, the temporal part and the spatial part optionally being combined using the Hadamard product.

15. The system of any preceding claim, wherein: the stimulation signal is configured to be phase-locked to a biomarker of the biological system and the at least one stimulation parameter comprises a stimulation phase relative to the biomarker; and / or the at least one stimulation parameter comprises a stimulation amplitude; and / or the at least one stimulation parameter comprises a stimulation pulse width; and / or the at least one stimulation parameter comprises a stimulation pulse frequency.

16. The system of any preceding claim, further comprising the stimulation unit.

17. The system of claim 16, wherein the stimulation unit comprises one or more electrodes for applying the stimulation signal to the biological system.

18. A computer- implemented method of generating a stimulation signal for brain or nerve stimulation, the method comprising: receiving time series biomarker data comprising measurement samples obtained by measurements performed on a biological system at different respective times, the biomarker data representing one or more biomarkers indicative of a condition of thebiological system that is affected by brain or nerve stimulation applied by a stimulation unit based on a stimulation signal, wherein the stimulation signal is defined by at least one stimulation parameter; and generating a stimulation signal for the stimulation unit based on the received biomarker data by using a model of the biological system to estimate an optimal value of the at least one stimulation parameter using the biomarker data, wherein the model is configured to define strengths of contribution of measurement samples to the estimate based on the times at which the measurement samples were obtained.

19. The method of claim 18, wherein the generating of the stimulation signal comprises generating a stimulation signal for which the at least one stimulation parameter has the estimated optimal value.

20. A method of performing brain or nerve stimulation, comprising: generating a stimulation signal using the method of claim 20; and applying stimulation to the biological system based on the generated stimulation signal.

21. A computer program comprising instructions which, when the program is executed by a computer, cause the computer to carry out the method of claim 20.

Images