Moving from phenomenological to predictive modelling: Progress and pitfalls of modelling brain stimulation in-silico

Brain stimulation is an increasingly popular neuromodulatory tool used in both clinical and research settings; however, the effects of brain stimulation, particularly those of non-invasive stimulation, are variable. This variability can be partially explained by an incomplete mechanistic understanding, coupled with a combinatorial explosion of possible stimulation parameters. Computational models constitute a useful tool to explore the vast sea of stimulation parameters and characterise their effects on brain activity. Yet the utility of modelling stimulation in-silico relies on its biophysical relevance, which needs to account for the dynamics of large and diverse neural populations and how underlying networks shape those collective dynamics. The large number of parameters to consider when constructing a model is no less than those needed to consider when planning empirical studies. This piece is centred on the application of phenomenological and biophysical models in non-invasive brain stimulation. We first introduce common forms of brain stimulation and computational models, and provide typical construction choices made when building phenomenological and biophysical models. Through the lens of four case studies, we provide an account of

Introduction "All models are wrong, but some are useful" -George E P Box. A fundamental question in neuroscience is how the brain's structure and function may be manipulated to counteract injury or disease. One effective way to establish the relationship between brain and behaviour is to employ targeted electrical brain stimulation and study the subsequent 2/38 effects. While direct manipulations applied at different neural scales (from cells to circuits) enrich our knowledge, the translation to clinical applications is by no means trivial. This is because many of the mechanisms by which electrical brain stimulation influence brain function and behaviour remain elusive. A better understanding of the relationship between stimulation parameters and their behavioural and physiological effects is necessary for the effective deployment of this promising clinical and research tool. Computational modelling can help us understand the mechanisms underlying the effects of stimulation and how to optimise stimulation parameters to increase the efficacy of treatments.
Here we discuss the current state-of-the art in computational models of the effects of electrical brain stimulation.
We first introduce the most common methods of electrical brain stimulation and computational model design. We then consider how recent work uses phenomenological and biophysical computational models to gain insight on the effects of stimulation using a "case study approach" Finally, we critically assess the current state of the field and suggest possible future directions.

Electrical brain stimulation
The purpose of electrical brain stimulation is to extraneously manipulate neural activity through the generation of local electric fields. The generated electric fields can result in both local neurochemical changes and network-level changes [1] [2] [3] [4] [5]. Electrical brain stimulation can be broadly separated into invasive and non-invasive brain stimulation.
Invasive brain stimulation is commonly applied in the form of deep brain stimulation (DBS). DBS uses surgically implanted electrodes, which typically target subcortical structures [6]. The stereotactically guided implantation of the electrodes to the stimulation target results in highly precise spatial localisation.
DBS has a long history of trial-and-error underpinning its current use. Animal studies conducted as early as the late 18th century investigated responses to cortical stimulation, and in the 1950's, thalamic ablations were shown to reduce tremor [7]. In the late twentieth century DBS had become established as a treatment for Parkinsonian tremors [8], and now is an FDA-approved treatment for movement and psychiatric disorders, such as dystonia [9] and obsessive-compulsive disorder [10], respectively. Yet, over a century of trial and error has not resulted in a complete mechanistic understanding of how DBS operates. Recent work has posited several theories [11], but none have been conclusively embraced. Extending the application of DBS to processes that rely on distributed functions, such as cognition, requires a stronger mechanistic understanding of an expanding suite of stimulation parameters and their effects on brain function.
Non-invasive brain stimulation (NIBS) has been lauded as a well-tolerated, safe, and cheap neuromodulatory 3/38 tool [12] [13]. NIBS encompasses an umbrella of techniques developed to manipulate brain function by applying current through the scalp. The most widely used methods include transcranial magnetic stimulation (TMS) and transcranial electrical stimulation (tES). TMS employs transcranial electromagnetic induction to generate an electric field [14] [15], whereas tES produces neuromodulatory effects via weak currents applied to the scalp [16]. Two popular methods of tES include transcranial alternating current stimulation (tACS) and transcranial direct current stimulation (tDCS) [17] [18] [19]. Temporal Interference (TI) is a novel NIBS approach for sculpting electrical fields to enable non-invasive neural stimulation at depth [20] [21]. While NIBS eliminates the risk associated with surgical implantation of electrodes, it is inherently associated with lower spatial resolution. One way to address this limitation is via the usage of multiple electrodes to enable more focal stimulation (commonly referred to as high-definition (HD)-tES) [22]. Another development in the delivery of NIBS is the more widespread availability of devices with several independent channels for current delivery to allow for multi-site stimulation. This move is supported by evidence that the effects of NIBS can be enhanced when stimulation is delivered to multiple regions in the same functional network [23] [24] [25]. Despite brain stimulation's utility as a clinical and research tool, an undeveloped mechanistic understanding of the parameters underpinning its function limit its effective deployment [26] [27]. Some parameters are outside of experimental or clinical control, including individual anatomy, the organization of local inhibitory and excitatory circuits, neurochemistry, and genetics [28] [29]. However, many parameters can be controlled, including the type of stimulation applied, the duration/intensity of stimulation, and the targeted region(s) and/or brain state [28] [29].
One way to both cultivate a mechanistic understanding and optimise delivery of brain stimulation is to utilise computational models of stimulation. Before discussing some specific examples, we first set the broader context of the types of models that have been used in neuroscience.
Introducing models of the effects of stimulation on brain function Brain stimulation manipulates neural activity at different topological levels (micro, meso and macroscale) and over different timescales (seconds to days or weeks). Across these levels, the effects of stimulation can be thought to follow a casual pathway, as framed by Bergman and Hartwigsen [26]: Stimulation is delivered → inducing changes in neural function → that can manifest in physiological changes that impact behavioural outcomes [Fig 1]. Computational models are able to simulate some or all stages of the causal path of stimulation.
Models that consider the first stage in the causal path, the delivery of stimulation, often concern the effects of applying current directly (via DBS) or indirectly (via NIBS) to the brain. Models of this nature often focus on how current spreads through skin, skull, and brain tissues, and estimate electric field distributions in the brain. Illustration of the causal path of stimulation. Brain stimulation can be administered invasively or non-invasively, both of which effect neurophysiology at the micro, meso, and macroscale. Then, the stimulation-induced changes in the brain exert downstream effects on behaviour. The brain images were made using BrainNet [30] Comprehensive reviews on these computational forward modeling studies using finite elements can be found in [31] [32] [33].
Models of the second and third steps dealing with changes in neural function and behaviour can be built at the micro, meso, and macrocale. Models focusing on the microscale level investigate neurochemical and cellular-level mechanisms [34]. At the mesoscale, models focus on brain structures in between the cellular and whole-brain resolution, including white matter tracts, small circuits, and cortical tissues [35]. Models that explore macroscale levels investigate the effects of brain stimulation on multiple brain regions, often in the form of neural networks, and/or their consequences on behaviour [36] [37] [38] [39]. These models will be explored in greater detail in here.
Macroscale models often employ phenomenological or biophysical models, which group populations of neurons into functional units and consider how they interact. Macroscale models often capture oscillatory behaviour, and come in two popular versions: the neural mass model and the neural field model. Neural mass models characterise activity over time only, whereas neural field models describe how properties of neural units (such as average depolarization of a neural population) evolve over both space and time. Three models that have played a seminal role in mathematical neuroscience and have been used to model the effects of electrical brain stimulation on (different) aspects of neural

5/38
The original Wilson Cowan Model (WCM) was defined in Wilson and Cowan's seminal 1972 paper [40]. It is a mean-field biophysical model that represents the properties of a spatially localised population of neurons, rather than individual neurons. The original WCM is composed of two coupled nonlinear equations that model the firing rate of excitatory neurons, E(t), and the firing rate of inhibitory neurons, I(t), within a brain region as a function of time t respectively. A sigmoidal transfer function, also known as an "input-output" function, δ, determines how the input current to the subpopulations is transformed into neural firing rates. Connectivity or coupling coefficients, c i , where i =1 to 4, represent the average number of excitatory (i =1,3) and inhibitory (i =2,4) synapses per population. There are different external inputs to the excitatory and inhibitory subpopulations -P and Q, respectively -to account for cell type-specific afferents to the population. The equations describing the time evolution of E(t) and I(t) are then where r is a constant that defines the refractory period, and α(t − t ) models the decay of cell inputs with time. The parameters of a WCM are often selected so that the system has three potential states: a low fixed point, a high fixed point, or a limit cycle. Fixed points, also known as steady-states, are largely uneventful [40]. At the low fixed point the firing rates of the subpopulations are stable and at a low level, creating a baseline state. At the high fixed point, the firing rates of the subpopulations are stable at a high level. Limit cycles represent oscillatory activity emerging from a balanced excitation and inhibition [40], and have been used as models for the gamma-frequency oscillations (≈40 Hz) detected from neuronal populations at the level of local field potentials. Additionally, these limit cycles have been extended to address pathological oscillations, reflecting excess synchrony that arises during epilepsy and essential tremor [43] [44]. In network models of the brain, the dynamics of each node may be represented as a WCM operating in the limit-cycle regime, resulting in a system of coupled oscillators [45].

6/38
The dynamics of a brain region exhibiting oscillations in local field potentials can be further reduced from the WCM to a simple phase oscillator [46]. Each phase oscillator in isolation has very simple dynamics, being described by its phase rotating at its natural frequency. However, when phase oscillators are coupled together, their interactions are described by the Kuramoto model of coupled oscillators: In the Kuramoto model of coupled oscillators, θ i (t) represents the phase at time t for node i, N is the total number oscillators, ω i is the intrinsic frequency of oscillator i, and K is the global coupling coefficient. As K increases the solutions of the model tend towards phase synchronisation, where "synchronisation" refers to the stability of phase relationships. If all nodes have the same phase so that θ i (t) − θ j (t) = 0 for all i and j then the oscillators are synchronised in-phase. If the phases of different nodes differ, but the phase differences do not vary with time, then the oscillators are synchronised out-of-phase, or "phase locked". In general, the phase coherence of a population of oscillators can be captured by the Kuramoto order parameter. The order parameter r defines the relationship between synchrony (ρ) and the average phase (ψ) for N nodes The order parameter ranges from 0 to 1. When r =1, nodes are considered in-phase, and when r =0, nodes are perfectly out of phase. For further reading, see the review written by Bick et al [47]. The review covers Ott-Antonson and Wantanabe-Strogaz reductions, which allow researchers to analytically solve solutions for simple systems of Kuramoto oscillators (rather than rely on computer simulations) [47]. While these reductions are helpful to bridge the microscale to the macroscale, they do not extend to systems with added complexity, such as in the presence of heterogeneous couplings and time delays, which have been shown to modulate network activity not only in space and time, but also in the frequency domain, which is a crucial feature of brain activity [48].

7/38
The Stuart-Landau oscillator is a reductionist approach to model population-level oscillations in neural activity considering both the phase and the amplitude dynamics. The original equation was proposed in 1944 by Landau [42] to describe a system exhibiting a Hopf bifurcation, i.e., a system that transitions from a stable fixed point to a limit-cycle upon a change in parameters [49]. The canonical form defined by Stuart and colleagues in 1960 [50] [51] [52] is defined as: Here A is a complex number expressing the evolution of the amplitude |A| and the phase φ over time as A = |A|e iφ . The parameter σ expresses the growth rate (or decay rate if σ is negative) and the real part of l is known as the Landau constant [42]. Similar to the WCM, for some values of σ and l the SLO exhibits either a fixed point equilibrium or a stable limit cycle. Importantly, if the system is in a fixed-point regime but sufficiently close to the bifurcation point (sub-critical), it exhibits a damped oscillation with decaying amplitude in response to input. As such, sub-critical SLOs perturbed with noise can be used to approximate the transient and short-lived oscillations detected in neural activity.

Four common construction choices for computational models
So far, we have introduced different types of brain stimulation and the causal path of stimulation's effects on the brain and behaviour. We described the different topological scales at which the brain is commonly modelled, and introduced three common macroscale models of brain function -WCM, KO, and SLO. Direct comparison of different modelling studies is difficult, as each step in constructing a model, from defining nodes and edges to parameter choices for the oscillatory equations, requires choosing from a number of options (see Fig. 2 for definition of nodes and edges). Nevertheless, there are commonalities, and here we list and briefly discuss four common construction choices used when building phenomenological or biophysical models of the effects of brain stimulation [Fig 2].

Network nodes describe brain regions
Most studies to date opted for atlas-based parcellations of the brain to relate brain regions to nodes in a network.
Studies have proposed that including cerebellar and subcortical structures may increase the biophysical relevance of models [53] [38], as well as accounting for differences in tissue excitability. For example, in Giannakakis et al [38] differences between healthy and epileptic brains were distinguished by modifying the internal weight values in the WC oscillators. Another consideration when defining model nodes is that atlas parcellations based on 8/38 Illustration of the four common construction choices made when developing computational oscillatory brain models. DWI = diffusion weighted imaging. Brain images were made using BrainNet [30] and MRIcroGL 9/38 cytoarchitectural characteristics may define regions with non-uniform functional characteristics. Models that group independent functional subunits into a region do not capture the individual contributions of each subunit, but rather a net misrepresentation of both. An alternative approach involves parcellating the brain into functional regions [54], which can be done on a group or individual level [55]. Functionally derived parcellations can reduce confounds of overlapping representations of regional activity that arise when independent functional subunits are grouped within the same structural/anatomical regions [55].
Another consideration when researchers decide whether to use a structurally or functionally defined atlas for dynamical models is provided by the work of Domhof et al, who used Kuramoto and Wilson Cowan models to systematically investigate how several parameters, including atlas parcellation, influenced the specificity and reliability of simulated data to participants' empirical data [56]. The authors showed there was a substantial, significant effect of atlas type (structurally or functionally defined), particularly on how well the simulated structural-functional relationships captured empirical ones. The simulated data generated using a structurally defined atlas more reliably related model outputs to the empirical data, whereas the functionally defined atlas showed improved specificity of simulated data to participants' empirical data [56]. Thus, the systematic influences of atlas parcellation technique observed by Domhof et al should be considered when designing and implementing modelling studies.

Network edges and their weights
Connections between nodes are often derived from diffusion weighted imaging (DWI) to create models that capture connectivity among functional units and how stimulation propagates through the brain. DWI measures white matter density as a function of water anisotropy across imaging voxels, and is well-suited for detecting long, continuous projections of axonal bundles, such as the corpus callosum [57], corticospinal, and frontal aslant tracts [58]. Work by Nowak and Bullier [59] showed that postsynaptic responses to grey matter stimulation are due to axonal, but not cell body, activation. Recent work confirms and expands these results, showing that DBS near white matter (WM) tracts are more likely to deliver excitatory effects [60] and the efficacy of non-invasive brain stimulation depends on WM structure [28] [5] [61]. Moreover, modelling studies found that stimulation of WM tracts with coherent fibre orientation result in more focused electric fields [62] [63]. Therefore, modelling weighted connections between regions should better simulate how stimulation will propagate through the brain. In some studies, the signal propagation time between brain areas is included by considering a time delay associated with each edge. Time delays can be approximated by considering the distance between brain areas [45] or the average length of the fibres [64] and assuming a conduction velocity. While accounting for the effect of time delays in conduction velocity and interregional communication increases the biological relevance of the work, there are a wide range of conduction speeds that generate simulated brain function that reproduce features of empirical brain function [65].
Models that use participant's structural or functional connectivity to inform edges and their weights are referred 10/38 to as personalised Brain Network Models (BNM) [66]. Simulations of brain stimulation on BNM thus reveal individual responses to stimulation, and provide insights as to how individual connectivity influences the effects of brain stimulation [66]. This enables the identification of key nodes that, when targeted by stimulation, innervate whole-brain networks. One such study was conducted by Muldoon et al and is explored later in this work [37].

Local brain activity is described by dynamical models
The dynamics of each node are then described by a mathematical model such as WCM, KO, or SLO, all of which are capable of reflecting the stochastic, complex variability in brain activity [53]. Although the Kuramoto model focuses solely on the interaction between the phases of oscillatory regions, it has the advantage of being analytically tractable. We note that spiking (integrate and fire) models, traditionally used for micro or mesoscale investigation, can also be used in macroscale models to generate dynamic asynchronous and/or oscillatory behaviour comparable to the phenomenological models covered in this work, but the additional degrees of complexity make the numerical and analytic solutions less tractable [67].

Parameter values are selected by fitting model outputs to empirical data
Parameter values are selected that align model outputs with empirical neural phenomena. Parameter values are often informed by literature (as in Muldoon et al [37]) or searching through a range of parameter values (as in Deco et al [39] or Weerasinghe et al [36]). The choice of parameter values greatly impacts model outputs [64] [65], and thus the validity of conclusions drawn from modelling studies. We end this work by introducing how machine learning approaches may improve parameter selection.
In the next two sections we provide case studies of work using computational models to answer either mechanistic or predictive questions of the effects of brain stimulation on brain function and/or behaviour. We have focused on a few particular exemplars which are representative of the use of WC, Kuramoto, or Stuart-Landau models, and highlight how the four common construction choices are implemented in current research.

Mechanistic modelling of stimulation effects
In this section, we consider two models that aim to develop a mechanistic understanding of the effects of stimulation. Specifically, the system-wide impact of regional stimulation [37] and the long-term effects of stimulation [38].

11/38
"How to maximise the system-wide impacts of regional stimulation?" Different psychiatric [68] and neurological [69] conditions are characterised by pathological brain function, the mitigation of which is the goal of neuromodulatory interventions. Indeed, a large body of empirical evidence shows that focal brain stimulation influences whole-brain network dynamics [70] [71] which attenuate the effects of neuropathology [8] [9] [10]. Therefore, the extent to which stimulation affects local and global brain dynamics is of important clinical utility. Further characterisation of how underlying structural connectivity influences the effects of regional brain stimulation on local and whole-brain function is an important step in identifying targets suitable for stimulation to attenuate pathological brain function.
Recent, state-of-the-art modelling [72] [73] and empirical [74] work shows focal deep brain stimulation influences cortical and subcortical network communication. For example, work by Bansal et al constructed personalised BNM to investigate the role of interregional connectivity on the whole-brain effects of stimulation, and showed that interindividual variability in structural connectivity explain local and global effects of brain stimulation [75].
Moreover, measures of the local and global effects of simulated stimulation explained differences in participant's performance of less and more complex tasks, respectively [75]. Collectively, the work by Bansal et al provides a compelling case that individual structural connectivity provides a scaffold upon which stimulation exerts its influence.
An example of research addressing the relationship between the functional and structural properties of stimulation targets and their local to global impacts is found in the first case study presented in this work, which features work by Muldoon et al. The authors modelled the effects of regional stimulation on whole-brain networks [37]. Their framework is based on Network Control Theory (NCT) [76] which identifies regions that are influential in shifting neural network dynamics; that is, a region's controllability [76] [37]. Controllability originated from engineering and physics research, and evaluates the potential that an external input to a part of a system can drive the whole system's dynamics to a target state [77]. In a neurological context, a target state refers to the desired pattern of activity or connectivity which the stimulation is intended to induce. For example, if stimulation is delivered in hopes of shifting aberrant dynamics to emulate those of healthy controls, then the pattern of activity or connectivity in healthy controls is considered to be the target state.
The authors first focused on the system-wide impacts of regional stimulation. They subsequently applied their results to the neural networks that underpin cognitive function. Here, we will focus on the first part of their study to provide an example of how computational models were constructed and employed to run a systematic assessment of how individual connectivity influences the controllability of stimulation targets.
The authors constructed their model by parcellating the brain into discrete regions using anatomical scans obtained with MRI [Fig 2Ai]. These regions formed the nodes of the network. The network's edges were defined 12/38 through the number of streamlines between regions estimated using DWI [Fig 2Aii]. This resulted in a weighted, undirected, structural connectivity matrix for each study participant, in line with previous literature [78]. This was used to inform a personalised BNM, where individual connectivity is used as a scaffold for the dynamical components of the model.
The dynamics of each brain region was described by a WCM [Box 1 ] [Fig 2Aii]. We note that the authors only included excitatory inputs to the excitatory subpopulation, assuming that most long-range connections are excitatory.
The time delay in communication between regions was set to be proportional to their physical distance to ensure activity propagated through the system in a biologically relevant manner.
Having defined nodes and edges, the authors then selected their model's parameters. They kept all region-specific parameters and the sigmoidal transfer function the same as in the canonical WCM [37] [40]. The authors selected the global coupling coefficient so that the dynamics of the system were close to the transition from steady, baseline activity, to a state where the firing rate at most nodes oscillated at around 20Hz. The transition point was dependent on structural connectivity and Muldoon et al found that the transition value of global coupling coefficient had low within-subject variability and high between-subject variability, which highlights the importance of personalised BNM [66].
Though Muldoon et al do not specify which modality of brain stimulation they are modelling, we infer that they are investigating the general effects of brain stimulation. Muldoon et al modelled the effects of brain stimulation by changing the current input to the system via the excitatory subpopulation, P, for a target region. The "stimulation" was delivered as a single pulse for 1 second, which caused the target region to transition from a baseline non-oscillatory state to an oscillatory state operating near 20 Hz [ Box 1 ].
In Wilson and Cowan's paper introducing the model, they related increased population firing due to an increase in the input current P to primate studies that show increased firing due to constant motor stimulation [40].
The local changes in the target region then resulted in changes in the dynamics of other regions. A functional connectivity matrix one second before and during the stimulation was constructed by considering the cross-correlation in firing rates between all pairs of regions.
The authors measured the effects of stimulation by defining the "functional effect" and the "structural effect".
The functional effect measured the difference between the functional connectivity matrix before and during stimulation. The structural effect measured the change in the difference between the structural and functional connectivity matrices before and during stimulation. The "functional effect" of a region is a measure of its ability to shift whole-brain activity and connectivity. The "structural effect" is a measure of the extent to which a region's structural connectivity accounts for its influence on neural network dynamics.
Key results from their study were that first, individual variation in structural connectivity lead to variations in 13/38 the neural/physiological response to stimulation, further supporting the use of BNM as a method to understand how structural constraints influence the functional effects of brain stimulation [66]. Second, NCT identified nodes with high controllability. Regions of low structural effect (regions with high structural connectivity) induced global changes in the brain state, whereas regions of high structural effect (which were more constrained due to sparse structural connectivity) induced more focal effects [37]. "What are the long-term effects of stimulation?" As evidence accumulates that tDCS [80] [81] and other brain stimulation modalities [82] [83] have clinically relevant benefits when delivered over more than a single session, insights characterising and understanding the source of the differential effects of brain stimulation in patients and healthy controls are increasingly important. For example, it is thought that the long-term effects of tDCS increase neural plasticity by shifting the balance of excitatory and inhibitory neurotransmitters [84] [85], though it is unclear how stimulation in clinical populations (who may already exhibit dysregulated excitatory/inhibitory dynamics) affects this balance. Yet, there are many practical limitations to characterising the effects of brain stimulation over 24 hours, let alone over week-or month-long interventions. Models of the long-term effects of brain stimulation in clinical populations are increasingly needed, as they mitigate many of the difficulties inherent to long-term studies and can inform how stimulation differentially influences clinical populations.
In our second case study answering a mechanistic question, Giannakakis et al modelled the effects of stimulation on structural connectivity and neural network dynamics over different time periods, throughout 30 minutes of stimulation ("short-term"), and 24 hours after ("long-term") [38]. Using structural data from healthy controls and people with epilepsy, they considered how differences in connectivity and inherent regional excitability influence the effects of stimulation [38].
The authors parcellated the brain into discrete cortical and subcortical regions as defined by the Desikan atlas [86] [ Fig 2Ai]. As in Muldoon et al, personalised BNM were created by using streamline tractography data to inform the coupling between regions [Fig 2Aii]. However, a key difference from [37] is that the structural connectivity was updated during and after the simulation via a Hebbian learning rule where the connection strength is increased 14/38 between two nodes if increases in firing in one node precede increased firing in the post-synaptic node; otherwise, pairwise connectivity decreased.
The dynamics of the nodes of the network were modelled as modified WCM based on their previous study [87] [ Fig 2Aiii]. The author's version of the WCM deviates from the canonical by including two inhibitory subpopulations as opposed to one, to model the effects of both subtractive and divisive inhibition. They cite two main reasons for including two inhibitory subpopulations. First, by including two inhibitory subpopulations, the model exhibits complex dynamics without abrupt transitions to chaos [88]. Papasavvas et al argue that this increased complexity better reflects the brain's ability to balance long-range synchronous activity [89] [87]. Secondly, this increased complexity also better models selective attention [90], sensory [91], and motor processing [92]. The study by Giannakakis et al was selected as a case study in this work because the inclusion of the second inhibitory population exemplifies how canonical models can be improved and adapted.
Next, the internal weight values between the excitatory and inhibitory subpopulations were chosen to model the differences between healthy and epileptic regions of the brain. The internal weights for the healthy regions were selected to give high amplitude oscillations for both the excitatory and inhibitory subpopulations, with the oscillatory amplitude tapering off after a few hours [38]. Epileptic regions were modelled as having a greater excitatory and lesser inhibitory weights than healthy regions, in line with the nature of the pathology [38].
Other differences from the canonical WCM were the use of a logistic instead of a sigmoid input-output function.
Within the input-output function, the inputs to both inhibitory populations are the same, but the inputs from the inhibitory populations to the excitatory reflects their differing physiological targets. Divisive inhibition occurs through targeting the soma, while subtractive inhibition occurs at dendrites [93] [94].
Similar to the study conducted by Muldoon et al [37], the authors modelled stimulation by changing the value of P, which represents excitatory afferents to the region. During stimulation, the authors decreased the value of P (i.e., excitatory input) by 50% for 3 target regions, simulating the possible inhibitory effects of cathodal tDCS [38] [ Box 1].
The target regions were the amygdala, hippocampus, and parahippocampal gyrus and were selected because of their known role in inciting seizures. While cathodal tDCS would not be a feasible method of stimulating subcortical regions, TI stimulation can non-invasively target stimulation to subcortical regions [20]. Recent work has shown TI stimulation of the anterior hippocampus modulates task-related hippocampal activity and connectivity [95]. Thus, inferences gained from Giannakakis et al's effects on subcortical regions may inform empirical work employing TI.
Three scenarios were considered. In this first scenario, the network structure was based on healthy controls and all regions were modelled as healthy regions. In the second scenario, the network structure was based on epilepsy patients and the target regions for stimulation were modelled as epileptic, more excitable, regions. For the third scenario, the network structure was based on epilepsy patients and all regions were modelled as healthy regions.

15/38
As a result of their simulations, Giannakakis et al predicted that the more excitable, epileptic group had a larger change in whole-brain and local connectivity than healthy controls at the end of the 30 minute stimulation session.
Over the next 24 hours, although the trajectory of effects of stimulation on whole-brain and local connectivity differed, ultimately, the changes in connectivity were similar between the healthy control and epileptic groups. They conclude that the temporal dynamics between the two groups must also be different, considering the effects of stimulation decay at varying rates for each model (due to the different stimulation-induced behaviors) for them all to approach equilibrium within 24 hours.

Modelling to optimise the delivery of stimulation
In this section we discuss two studies that address predictive questions on how to optimise closed-loop stimulation [36] and how to force a transition between brain states [39].

"How can we optimise closed-loop stimulation?"
Closed-loop stimulation systems use biomarkers to inform the timing of delivery of stimulation, and have the potential to increase the efficacy of brain stimulation by delivering stimulation protocols when markers of pathology are present, which offers both greater precision and efficiency compared to continuous DBS approaches by mimicking endogenous circuits more closely [96].
There are a number of studies investigating computational models of closed-loop stimulation [97] [98] [36]. Here we focus on the study by Weerasinghe et al [36] who sought to devise a closed-loop DBS strategy that most effectively reduces pathological neural oscillations. To do this, they investigated whether phase, amplitude, or phase and amplitude of neural oscillations provided the best biomarker to inform delivery of closed-loop DBS. In the first part of the study, the authors used a theoretical model to predict changes as a result of stimulation. Then, the authors applied these predictions to empirical data from essential tremor patients, and introduced a novel stimulation strategy.
The model consisted of a system of N Kuramoto oscillators [Box 2 ] where each oscillator was connected to every other oscillator. This model is phenomenological (i.e., descriptive), meaning that each oscillator included in the model does not explicitly represent a particular brain region, and the authors stated that each oscillator represented either a neuron or microcircuit [36]. The parameter k determined the coupling strength in the following equation:

16/38
The third term on the right-hand side describes stimulation. The stimulation of intensity I acts on all oscillators but its effect depends on the phase of each oscillator as through the phase response function Z(θ i ). The stimulation is turned "on or off" with X(t). When X(t) = 1, stimulation is on, and X(t) = 0 means stimulation is off.
In DBS, the response to stimulation can be described by amplitude and/or phase response curves (ARC and PRC, respectively). In the study, ARC and PRC described how amplitude and phase of pathological oscillations (such as tremor, which has been shown to operate within 3 to 8 Hz [99]) varied as a function of the phase at which stimulation was delivered. To capture the relationship between the ARC or PRC and stimulation the authors used the reduction method of Ott and Antonsen [100] [47], which in the limit of large N for an instantaneous pulse and a sinusoidal phase response function, gives analytical expressions that describe how the ARC or PRC change over time.
Weerasinghe et al then simulated a population of N = 3000 oscillators using pulsatile stimulation at ω = 30 Hz and a range of phases. They found that stimulating at phase π resulted in the largest change to the ARC, a result of the strong correlation between the ARC and the negative gradient of the PRC.
Then, the authors also explored how harmonics in the phase response function, Z(θ), influenced the effects of stimulation. They found if Z(θ) contains a dominant first harmonic, and stimulation is applied when the amplitude of oscillation ρ is low, then the effects of stimulation are magnified.
These results provide two predictions which they tested by examining responses to phase-locked DBS in experimental data collected from patients. In this experiment, DBS was applied at a range of phases and the effect of DBS was measured using hand accelerometer data, since the hand tremors of essential tremor patients parallel the neural oscillations that underpin the movement itself (both occurring around 5 Hz) [101]. Though only one out of nine subjects fit their predicted response to theoretically-driven, closed-loop stimulation, it is highly commendable that Weerasinghe et al empirically tested their model-driven predictions.
Research building upon the work of Weerasinghe et al may employ models developed with more consideration for the role of brain regions and circuits underpinning motor disorders, which are often characterised by well-known pathology in the basal ganglia [102]. For example, Meier et al [103] developed a personalised BNM that incorporates a spiking network model of the basal ganglia with a whole-brain model created using The Virtual Brain [104], an online platform that enables researchers to input subject-specific MRI data to inform mean-field oscillatory models of brain function. The authors showed simulated DBS to either the subthalamic nucleus or globus pallidus internus (GPi) restored normative thalamic and cortical activity, providing a proof-of-concept that a multiscale model can be used to develop and test new stimulation paradigms.
Recent modelling work has made advances in developing a mechanistic understanding of how oscillatory synchronisation results in metastable, brain state-like behaviour [117], and how brain states dynamics change when moving from wakefulness to sleep [118] [39] or the response to psychedelics [119]. For example, work by Bansal et al used WCO to develop personalised BNM with the aim of characterising the states that emerge in response to simulated brain stimulation of each node [120]. The authors showed that there is a large variety of states that emerge in response to stimulation, which can be grouped into three categories -coherent states (which exhibit global synchrony), metastable states (which exhibit low levels of synchrony and little network structure), and chimera states (which are dominated by groups of synchronised regions) [120]. Regional stimulation was most likely to evoke chimera states, inviting future work to explore whether prioritising target states exhibiting chimera-like properties provide better stimulation targets than metastable or coherent states.
Here we focus on a study by Deco et al in which a modified SLO system was employed to model the activity of two different brain states (sleep and wakeful resting state) [39]. The authors used their model to predict what stimulation parameters could force a transition between these two states.
First, Deco et al characterised each of the two discrete brain states -sleep and wakefulness -into three metastable substates (patterns of brain network activity that characterise a state) [39]. This was done using the Leading Eigenvector Dynamic Analysis (LEiDA) developed by Cabral et al [121]. LEiDA characterised the temporal evolution of dynamic functional connectivity through calculating the leading eigenvector of a phase coherence matrix.
This simultaneously reduces dimensionality, improves signal to noise ratios, and allows for the smooth characterization of recurring patterns over time. Moreover, the metastable substates identified by LEiDA have been shown to capture functional organisation spanning the coherent, metastable, and chimera categorisations of states observed by Bansal et al [120].
Second, whole-brain models were generated and fitted to the probability distributions of all the substates [Fig 2Aiv]. Similar to [37], the authors parcellated the brain using the AAL atlas, this time into 90 regions, each of which became a node in their network [Fig 2Ai]. Connections between nodes were set by a structural connectivity matrix created via tractography techniques [Fig 2Aii]. The dynamics of each node was described using a modified 18/38 LSO [Fig 2Aiii], written in Cartesian form as Here, x n modelled the BOLD activity for each node. The intrinsic frequencies ω n were 0.04 to 0.07 Hz, in line with the average peak frequency of the BOLD signal in each node. The coupling parameter C np measures the strength of the coupling node p to node n. For a n ≤ 0 the baseline state is a stable fixed point dominated by noise, which the authors describe as "a low-activity noisy state" [39]. At a n = 0 this baseline state loses stability and transitions to a self-sustained oscillatory state. Note this type of transition is a supercritical Hopf bifurcation, as discussed in the Section "Four common construction choices for computational models".
Fitting the model to the empirical data was done in two stages. First, the authors updated the effective connectivity coefficients C np to optimise the fit between the modelled and empirical time-averaged phase coherence matrix of the sleep and wake conditions. Second, the global coupling parameter, G, was tuned to fit the SLO activity to the empirically derived probability of occurrence of the three substates. Previous work by Deco et al showed that the global coupling parameter played a critical role in shaping simulated brain network function to fit empirically derived measures of brain function [122], and recent work by Treibkorn et al [65] evaluated the importance of conduction velocity, global coupling, or graph theoretic properties in shaping modelling dynamics that resemble empirical dynamics. Assessment of the first two free parameters was conducted using a systematic, exhaustive search of all combinations of parameter values, and showed that individuals typically exhibit a narrow range of values for global coupling parameters that generated simulated dynamics close to empirically derived ones, in line with previous work [122]. Thus, by tuning the global coupling parameter, G, Deco et al ensured the simulated data closely emulated properties of the empirical data.
In the final step they demonstrated that they could force a transition between the two states of wake and sleep through stimulation. Stimulation was administered bilaterally to pairs of nodes, so if a node on the left was stimulated, then so was the corresponding node on the right. Stimulation was modelled by changing the bifurcation parameter a n for target nodes. Increasing values of the bifurcation parameter resulted in stimulation-induced oscillations; conversely, decreasing values of the bifurcation parameter modelled noise. This provided a contrast between stimulation parameters that serve to either reinforce or disrupt brain activity.
Furthermore, the authors investigated whether stimulating multiple regions with weaker stimulation more effectively forced a transition between brain states. Their best fit for multi-site stimulation resulted from stimulating 19/38 two pairs of regions.
Though the modelling study by Deco et al [39] was not empirically validated, a study by Chen et al [123] directly compared whether stimulation of the primary motor cortex or multi-site stimulation of the primary motor cortex, premotor cortex, and supplementary motor area increased a measure of the synchrony between cortical and muscular activity. The authors showed that multi-site stimulation resulted in significantly higher and longer-lasting measures of cortico-muscular synchrony, thus corroborating the simulated results from Deco et al that multisite stimulation can more effectively drive a neural system towards a desired state [123] [39].
Steps to move from phenomenological to predictive modelling This work began by introducing different types of brain stimulation, how brain stimulation affects brain function, and how macroscale computational models can represent such relationships. We then described studies that use phenomenological or biophysical models to simulate brain stimulation and its effect on neural function with the aim of answering either a mechanistic or predictive question concerning the effects of brain stimulation. The results from the Muldoon et al study not only suggested that differences in structural connectivity could explain some of the interindividual variability that hounds NIBS studies [27] [124], but also provided a framework for deciding which regions to stimulate based on their controllability [37]. In the study by Giannakakis et al the authors suggested variations in excitability and connectivity may underpin the differences in the brain's immediate response to stimulation [38]. Weerasinghe et al suggested a novel form of closed-loop DBS called "Hybrid DBS". This would deliver high frequency stimulation when endogenous oscillations surpass a threshold amplitude [36]. Finally, results from Deco et al suggest popular one-site or single-channel stimulation setups may not be as effective as those using multi-channel setups concomitantly targeting multiple regions [39].
There is no dispute that the human brain is complex, with trillions of neural synapses packed into 1.41 dmm 3 [125]. Initiatives such as the Virtual Brain [104] and NiMARE [126] aim to capture this complexity in-silico by developing detailed computational models. But such models are limited by computational power and current knowledge of neurological structure and function. Instead, macroscale approximations take a coarse-grain approach, supposing that clusters of billions of neurons act as a single oscillatory unit. These oscillatory models are grounded in physiological truths that have not changed substantially since Wilson and Cowan's seminal paper [40]. Such models are complex enough to exhibit some of the brain's macroscale dynamical patterns, yet simple enough to offer a mechanistic insight into brain dynamics, and to suggest ways to optimise stimulation delivery. Computational models hold real promise as ways to understand both how the brain works and how brain stimulation can be employed with can be used to define the relationship between the neural correlate and outcome response. There are inherent complexities when choosing which neural correlate best captures the experimental spatio-temporal dynamics of the brain, and how it relates to behavioural or clinical outcomes. For example, models simulating brain activity captured with fMRI often aim to fit the functional connectivity between brain areas averaged over time [37]. However, recent work has focused on the temporal nature of brain dynamics by quantifying recurrent, metastable brain states rather than a static representation of network connectivity [39] [121]. Another method to select a neural correlate of clinical or behavioural outcomes is to use a feature extraction method, such as Highly Comparative TimeSeries Analysis (HCTSA) [127] [128]. HCTSA identifies neural features and ranks their ability to predict behavioural or clinical outcome measures. Deciding which neural correlate serves as an optimal target feature and determining how it relates to the behavioural or clinical outcome measure is not straightforward, yet it is a crucial step when developing model-driven targets for brain stimulation. provided can be grouped into two sets, where each set has a different approach to identify the optimal parameters.
The first set of parameter tuning methods (Fig 4i-ii) identifies the optimal parameter values by running the computational models using different parameter configurations, and then determines the combination of parameter values that generated simulated data that most closely resemble the experimental values. The second set identifies the optimal parameter values by searching the parameters space in the training phase of the model, during which the parameter values are updated at each iteration of the training according to specific criteria set by the optimisation approach selected (Fig 4iii-iv].  . We selected Nelder Mead simplex algorithms to represent local search as it does not require the model to be differentiable. Given the fact that some biophyical models, including the WC model, utilise partially differential equations, the Nelder Mead simplex algorithm may be a better local search method than GD. There are several factors that may influence which parameter tuning method is most appropriate, including the size of the parameter space (i.e., the set of values that parameters can acquire), the available computational resources, whether the functions in the models are differentiable, whether the set of potential solutions is continuous or discrete, and whether there are constraints in the values that the parameters can acquire. All these requirements and criteria can make the selection of the appropriate approach difficult. A further complexity in using parameter tuning approaches is that their implementation for the optimisation of computational models of brain function might not be trivial, considering that they are commonly developed and integrated within machine learning pipelines [134].
Therefore, adapting parameter tuning methods for these type of generative models may not be straightforward.
Finally, and most importantly, some machine learning methods (especially those in the second category) are not necessarily expected to identify parameter values within biologically feasible ranges. This is a common concern in NCT studies [37] [139], which have been shown to better capture theoretical controllability, rather than practical controllability [140]. Constraining stimulation parameters within ranges suitable for brain stimulation would ensure 24/38 that model-generated stimulation hypotheses are practical. Yet even this suggestion is not readily implemented, as the relationship between modelled and empirical stimulation parameters is not known. To constrain the simulated stimulation parameters to biologically relevant ranges, we suggest researchers implement an approach such as that described in Phase 1 to identify the response relationship between simulated and empirical brain stimulation.

Box 4: Parameter tuning approaches
Grid Search: The parameter space is divided into different grids [134] [138]. A model is trained for all the combinations of parameters within each grid, and the combination of parameters that leads to a model that generates the closest match between the simulated and empirical target features is identified as optimal. Two notable limitations of grid search are that it is computationally expensive, and it does not use information from past simulations in subsequent applications (e.g. if one parameter value was found to be suboptimal, it may still be used in following simulations if it is in another grid).
Random Search: Unlike the brute force approach of grid searching, in random search combinations of parameters are randomly selected as input to the model [134] [138]. As in grid search, parameters generating a close match between the simulated and target feature are saved. Random search requires fewer computational resources than grid search, yet performance is equitable [138].
Local Search: A loss function is defined as the difference between the experimental and simulated data obtained by the model. The aim of local search algorithms is to find the parameters that minimise this loss function; namely, to reduce to the minimum the difference between the simulated and experimental values. A commonly used local search algorithm is gradient descent (GD) [136]. GD is characterised by an iterative process in which the parameters p are updated gradually at each iteration. To find the minimum of the loss function L(x), the algorithm calculates the gradient at each iteration, which corresponds to the partial derivative of the loss function with reference to each parameter's value L(p). The gradient establishes the direction of the parameter's update, while its magnitude is commonly selected by the user through the learning rate and fine-tuned as hyper-parameter. Two limitations of GD are, however, that the function L(p) needs to be differentiable and that GD tends to get stuck at local minimum. Furthermore, if multiple local minima are present, running the simulation with the same configuration might lead to different results, due to the random initialization of parameter values. An alternative to gradient descent is the Nelder-Mead simplex algorithm, which is a derivative-free local search algorithm [141] [137] [142]. This algorithm represents the search area as a simplex with the number of points determined by the number of dimensions. At each iteration, the point of the simplex with the poorest performance is moved to a new position, resulting in a different simplex. This is repeated until a minimum is found. The algorithm, despite not requiring a derivable function, is computationally expensive at higher numbers of dimensions.
Global Search: Global search distinguishes itself from local search because it finds global minima across the entire parameter space. A commonly used global optimization approach is simulated annealing (SA) [143] [144] [135]. SA is a non-derivative approach that is capable of finding global minima in presence of multiple local minima. At each iteration, SA algorithms move in a random direction along the loss function, and then evaluate the results of the model. If the results improved from the previous iteration (i.e., the loss decreased), then the move is always approved. If not, the move is accepted with a certain probability lower than 1. The probability is regulated by both a parameter defined as Temperature T, which is kept high at the first iteration and then gradually reduced throughout iterations according to an annealing schedule, and the outcome after each move. At high temperatures, the probability of accepting a bad move is higher compared to lower temperatures. The advantages of this methodology are that it doesn't require the function to be derivable, it is capable to find global rather than local minima, and can be used both for discrete and continuous variables. It was successfully applied to fine tune parameters in other computational models [135], and was shown to perform better than Nelder-Mead simplex algorithms [145].

25/38
Conclusion Creating computational models that are predictive, rather than descriptive, would enable a more efficient exploration of the effects of stimulation. In-silico stimulation studies would reduce the substantial time, human, and financial resources associated with a trial-and-error approach [146], but until computational models of the effects of brain stimulation are effective at generating testable stimulation hypothesis, empirical research will remain the gold standard. This work began with the quote "All models are wrong, but some are useful" by George E. P. Box [147].
Yet without empirically testing the conclusions drawn from computational models of the effects of brain stimulation, we cannot answer Prof. Box's subsequent statement -"So the question you need to ask is not 'Is the model true?' (it never is) but 'Is the model good enough for this particular application?' " [148].