The role that choice of model plays in predictions for epilepsy surgery

Mathematical modelling has been widely used to predict the effects of perturbations to brain networks. An important example is epilepsy surgery, where the perturbation in question is the removal of brain tissue in order to render the patient free of seizures. Different dynamical models have been proposed to represent transitions to ictal states in this context. However, our choice of which mathematical model to use to address this question relies on making assumptions regarding the mechanism that defines the transition from background to the seizure state. Since these mechanisms are unknown, it is important to understand how predictions from alternative dynamical descriptions compare. Herein we evaluate to what extent three different dynamical models provide consistent predictions for the effect of removing nodes from networks. We show that for small, directed, connected networks the three considered models provide consistent predictions. For larger networks, predictions are shown to be less consistent. However consistency is higher in networks that have sufficiently large differences in ictogenicity between nodes. We further demonstrate that heterogeneity in ictogenicity across nodes correlates with variability in the number of connections for each node.

Mathematical models are increasingly being used to study the emergence of large-scale, spatiotemporal brain dynamics, which are typically recorded using techniques such as electroencephalography (EEG) or functional magnetic resonance imaging (fMRI) [1][2][3] . A common approach is to formulate models based upon large-scale brain networks in order to understand how dynamics are shaped by network connectivity and intrinsic node properties [4][5][6] . There is an extensive literature for this kind of modelling [7][8][9][10][11] , but a popular approach is to use ordinary or stochastic differential equations to model the temporal evolution of nodes (i.e. regions of brain tissue) in combination with an estimate of network structure and coupling equations, which define how nodes interact with one another 3,[12][13][14][15] . Thus, physiological mechanisms that are incorporated into these models include the presence (or absence) and weight of large-scale connections. The choice of model for nodes can be broadly split into two categories. So-called "physiological" models incorporate intrinsic node mechanisms that are derived from the properties of large regions of brain tissue 1,9,16,17 . On the other hand, "phenomenological" models do not explicitly model the physiological mechanisms of node dynamics but represent pertinent features of brain dynamics using more abstract or canonical forms. Nevertheless, they retain mechanisms relating to network connectivity [18][19][20][21][22] . Both approaches have been widely used to study the emergence of large-scale brain dynamics, either to better understand healthy brain functioning or the effects of disruptions associated with neurological conditions like epilepsy 8,[23][24][25][26][27][28][29] .
In addition to a fundamental understanding of spontaneous brain dynamics, models of the response of the brain to perturbations are also crucial in order to better understand sensory processing and responses to treatment 30,31 . For example, in order to design optimal treatments for neurological and neuropsychiatric disorders, we should seek to use models to understand the effects of treatment perturbations on healthy and abnormal brain dynamics 29 . A prototypical example is epilepsy, which affects approximately 50 million people worldwide 32 . Around a third of people with epilepsy do not respond to antiepileptic drugs or other treatments 33 , and for these people epilepsy surgery can be the only way to eliminate or mitigate seizures. However, long-term post-operative seizure freedom is only achieved in around 50% of patients 34 . In addition, many patients are not referred to surgery due to difficulties in identifying the epileptogenic focus 35 . Mathematical models of seizures have the potential to help us improve our understanding of epilepsy and provide quantitative prognoses for treatment outcome 36,37 . www.nature.com/scientificreports www.nature.com/scientificreports/ which converts net post-synaptic potentials into afferent firing rates, and R i is the coupling term, given by where α P is the global coupling strength and M is the adjacency matrix of the network representing the interaction between nodes (brain regions). The parameter values used to solve Eq. (1) and their biological interpretation are given in Table 1. The output of the model is given by (y 3 − y 5 − y 7 ), which corresponds to the membrane potential of pyramidal cells, resulting from the interactions between three populations of interneurons, one excitatory and two inhibitory 1 . The choice of parameters in this model places the system near a saddle node on invariant circle (SNIC) bifurcation 52 , such that transitions from a fixed point (background state) to a high amplitude oscillation (epileptiform dynamics) can arise due to noise. From a dynamical viewpoint, the value of the excitability parameter p P quantifies the distance from the bifurcation point and therefore contributes to the propensity of the system to transition from the background to epileptiform dynamics. This propensity is also influenced by the input from other nodes, which is quantified here by the coupling strength α P . Theta model. As described above, the dynamic mechanism underlying transitions between states in our implementation of the Physiological model is a SNIC bifurcation. The normal form of this bifurcation is given by the Ermentrout-Kopell canonical model, also known as "Theta-Neuron" or simply "Theta" model 48,53 . Originally proposed to describe neuron firing, this model has been used to represent large-scale neural masses embedded in networks 42 . When network connectivity is incorporated, this model can be described as 42 : ν 0 ,e 0 ,r Parameters of the sigmoid function ν 0 = 6 mV, e 0 = 2.5 s −1 , r = 0.56 mV −1 In these equations, θ i is the phase of node i, ξ is Gaussian noise with zero mean and standard deviation σ = 8, M is the adjacency matrix, θ s is the steady state of a single (uncoupled) node, given by and the parameters α T and p T represent the coupling strength and node excitability, respectively. For uncoupled nodes and in the absence of noise, the SNIC bifurcation takes place at p T = 0. When p T < 0 the stable attractor is the steady state θ s , while for p T > 0 a limit cycle emerges and is the only stable attractor of the system. The phenomenological Theta model presents some important advantages when compared to the physiological model described in the previous Subsection. The computational cost can be significantly lower, and the reduction in dimensionality and number of independent parameters make it more straightforward to understand and control the influence of connectivity and excitability in the emergence of epileptiform dynamics 42 .
Bistable model. A phenomenological model that is commonly used to simulate transitions between healthy and epileptiform dynamics in brain networks is a modified version of the normal form of the subcritical Hopf bifurcation 22,40,43,54 , which we refer to here as the "bistable model". The formulation of this model for a single node in the deterministic case is given by  where z is a complex variable that defines the two states associated with normal and pathological activity (stable fixed point and stable limit cycle, respectively), p B is the excitability parameter, and ω controls the oscillation frequency (we use ω = 20 as in Benjamin et al. 54 ). This model has a globally stable fixed point for p ≤ −1 and a globally stable limit cycle for p B ≥ 0. In the interval −1 ≤ p B ≤ 0, the stable fixed point (z = 0) and the stable limit cycle coexist. In this work, we represent the interaction between nodes through an additive coupling in the real part of z, as used in Petkov et al. 40 and Sinha et al. 43 . For z i (t) = x i (t) + iy i (t), the model is described by M is the adjacency matrix, α B is the coupling coefficient and ξ represents the Gaussian noise with zero mean and standard deviation σ = 1.85 (β = 0.01) 54 .
Quantification of state-switching dynamics. The aim of epilepsy surgery is to remove regions of brain tissue such that the brain can no longer generate seizures. As mentioned in the Introduction, Goodfellow et al. 25 introduced a model-based framework to quantify the ictogenicity of networks in terms of their propensity to generate recurrent state-switching dynamics, and thereby also quantify the reduction in seizures that would result from a model surgery. Specifically, Brain Network Ictogenicity (BNI) was used to quantify the propensity of a network to transition to the ictal state 40,45 . BNI can be evaluated in practice in different ways, but a useful method is to consider it as the average proportion of time the nodes of the network spend in the ictal state, compared to a reference period of time: Average time nodes spend in the ictal state Duration of the reference time (9) The definition above is most useful when the model of interest displays spontaneous, recurrent transitions between the healthy and ictal state, and is therefore applied in this study to the Physiological and Theta models. In the use of the Bistable model, an alternative calculation of BNI has been considered, which is to use the escape time from steady state to the limit cycle. Specifically, for the case of the Bistable model, we initiate the dynamics in the z = 0 fixed point (background state) and define the BNI as: Average time nodes take to transit to the ictal state Duration of the reference time (10) When placed into networks, the Physiological, Theta and Bistable models display spontaneous transitions between different types of solutions for a subset of model parameters 20,25,42 . For a given network, it has been shown that these state-switching dynamics can be found by varying the global connectivity strength or the intrinsic excitability of nodes 20,25,42 . In this study, we seek to obtain a more comprehensive quantification of the presence of different states in each model. To do this, we consider the dynamics of each model over a range of values of both the connectivity strength and the intrinsic "excitability" of each node (i.e. its proximity to a bifurcation, see Fig. 1). Specifically, we define a range of values of p and α. For a given network, we simulate each model, for each www.nature.com/scientificreports www.nature.com/scientificreports/ parameter value over a regular grid of 192 × 192 points using an Euler-Maruyama scheme with fixed time step, optimized individually for each model (between 1 × 10 −3 and 5 × 10 −3 ) and calculate BNI according to Eq. (9) (or (10), for the Bistable model), for each set of parameters. The result is a two-dimensional map of BNI values (see Fig. 1). This approach allows a more extensive quantification of changes in ictogenicity compared to previous studies 25,42,43 , taking into account the influence of the excitability and connectivity parameters at the same time.
Quantification of the effect of node removals. The definitions of BNI above provide a starting point to quantify the effect of removing nodes from networks, in a way that is pertinent for epilepsy surgery. Specifically, we define node ictogenicity (NI) as the change in BNI when a node is removed from the network: The subindices 0 and i refer to the complete network (before the removal of any node) and the network after the removal of node i, respectively. In previous work 25 , the reference state BNI 0 was defined to be a state in which the unperturbed (complete) network spent half of the time in epilieptiform dynamics (BNI 0 = 0.5). This value was originally chosen in order to try to optimize detection of changes in BNI. For a network in which the intrinsic parameters of every node are equal, such a state can be achieved by increasing the global connectivity strength from zero, until a state in which BNI 0 = 0.5 is attained. Here, we use a different and more comprehensive approach, computing statistics over the 2-dimensional map of BNI to calculate the change in dynamics as a result of the removal of a node. Specifically, NI is then calculated as the percentage reduction in average BNI over the whole map as follows: where the brackets 〈〉 represent the average over all simulations over the α × p diagram.
The greater the reduction in epileptiform activity after the removal of node i, the greater the value of NI i . The ictogenicity of a node can be negative if the removal of this node results in an increase in the average spiking (or a reduction in the escape time, for the Bistable model) of the remaining network.
The ranges of coupling and excitability for the three dynamical models are given in Table 2. The size of the windows in parameter space where the BNI is calculated for each model might influence the accuracy in capturing changes in the BNI, so we aimed to define regions that are large enough to include parameter values already considered in the literature, as well as to include whole regions of interest for the single node cases (like the bistability region for the Bistable model), but at the same time small enough to avoid considering wide regions where small or no changes in BNI are seen, consequently "diluting" the effects we wish to observe.
Comparison of perturbation effects. For each network realization and choice of model, the distribution of NI was calculated according to Eq. 12 (i.e. {NI i }, where i indexes nodes in the network). In order to compare how these distributions differ for different choices of model within the same network realization, we took into account the ranking of NI across nodes as well as their relative values. To do this, we use a weighted Kendall rank correlation measure 42,55 , defined by: www.nature.com/scientificreports www.nature.com/scientificreports/ where P (Q) is the number of pairs of nodes ranked in the same (inverse) order by both models. In order to quantify the differences between models in a more precise way, the sums in P and Q were weighted by the term where i and j refer to the nodes being ranked, and A and B represent the different models under comparison. The weighted Kendall rank is a number in the range [−1, 1], where 1 represents total agreement and −1 represents total disagreement (rank in the inverse order) between the two compared ranks. It is important to point out that the weighted Kendall rank does not take into account nuances of the shape of the NI distributions being compared beyond the rank ordering of NI values. However, given that here we are mainly interested in identifying the most ictogenic nodes, the rank is our measure of interest.
In this work we analyze binary directed networks with three, four and nineteen nodes. We consider all 13 3-node and 199 4-node networks that are connected and nonisomorphic (the latter were obtained using the software Nauty 56 ). In addition, 125 19-node random networks were selected by varying the probability of node connectivity between 0 and 1, choosing connected and nonisomorphic networks with an approximately uniform distribution of number of edges in the interval between 18 and 342 (minimum and maximum number of edges for 19-node directed connected networks).

Results
Networks with three and four nodes. In order to systematically analyze the role that the choice of model has on the effect of node-removals, we first studied small networks in which the full set of distinct topologies could be elucidated. There are 13 3-node networks and 199 4-node networks that are directed, connected and topologically distinct (nonisomorphic). For each of these networks, we calculated the distribution of NI (see Methods) for each of the three models of interest.
As shown in Fig. 2, the resulting NI distributions for 3-node networks can be divided in four groups. For the first network, shown on top of Fig. 2, the three nodes have quite distinct ictogenicity, where node 1 seems to act as a "controller" of the dynamics resulting from the interaction between the other two nodes, in a way that when this node is removed, the ictogenicity of the remaining network increases. For the second group, one of the nodes (node 2) clearly has a larger NI than the other nodes. For the networks in this group, node 2 has a larger degree (the sum of in-and out-degrees) than the other nodes, as well as an equal number of connections to nodes 1 and 3 (Fig. 2, second row). For the third group (Fig. 2, third row), two of the nodes have relatively high NI (nodes 1 and 3). These networks are characterized by having only one pair of nodes that are mutually connected. Finally, for the fourth group (Fig. 2, bottom row), all nodes have a similar value of NI. For each network in the fourth group, all nodes have the same number of connections. In each case we observe that the distribution of NI is equivalent for the three models, in terms of the ranking of nodes. The similarity between model predictions for 3-node networks show that, for these networks, the models are characterizing the networks in the same way regarding the contribution of each node to the generation of epileptiform dynamics. The NI diagrams presented in Fig. 2 are representative examples of each of the network groups. For the individual diagrams associated with all 13 networks see Supporting Information.
An analysis of NI distributions in 4-node networks is presented in Fig. 3(a). Since the number of networks in this case is large, we studied the average weighted Kendall rank (〈τ〉, see Methods) across all 199 topologically distinct networks. Figure 3(a) demonstrates that 〈τ〉 = 1 for all pairwise comparisons of the three models, indicating that, as in the case of 3-node networks, the response to node-removal perturbations is equivalent among the considered models.
Networks with a larger number of nodes. In order to better understand the influence that the choice of model has on the effect of node removal for larger networks, which we might consider more in line with those derived from clinical data 57 , we examined the effect of this type of perturbation in a set of 19-node networks. Since the number of topologically distinct networks with 19 nodes is too large to systematically explore, given the cost of computing NI, we quantified NI for 125 sampled random networks, taking into account a broad spectrum of topologies by stratifying our sampling according to the number of edges (see Methods). The average Kendall rank 〈τ〉 calculated over all 125 networks is presented in Fig. 3(b). It can be seen that the Physiological and Theta models show highest concordance in ranking, having an average τ (〈τ〉) of 0.84. The Physiological and Bistable models have 〈τ〉 = 0.72, and the concordance between the Theta and Bistable models is lowest, with 〈τ〉 = 0.56. Thus in larger networks, the choice of dynamics of the individual nodes plays a more significant role than in smaller networks.
In order to investigate the reasons for this, we examined the distribution of NI across networks and sought to assess how this may be associated with similarities in perturbation effects across models. The results in Fig. 2 highlight that some networks give rise to heterogeneous NI distributions, so we asked whether the agreement between the NI distribution of models is related to the extent of this heterogeneity. To quantify NI heterogeneity we calculated, for each network and each model, the difference between the highest and lowest values of NI (which we denote as ΔNI).  www.nature.com/scientificreports www.nature.com/scientificreports/ In Fig. 4, we demonstrate how the agreement between the three models (as quantified by 〈τ〉) varies as a function of the heterogeneity in the NI distribution of each model (as quantified by ΔNI). In this analysis, we normalize the values of ΔNI across all 125 networks to values between 0 and 1, independently for each model. We then considered the networks with normalized ΔNI > 0.05 and split them into four equally sized bins according to their value of ΔNI. We used each model in the comparison pair as the reference for this, so that we could mitigate potential differences in scaling of NI for different models. For all of the three pair-wise comparisons it is clear that the more heterogeneous the NI distribution is, the better the two models being compared agree in ranking  www.nature.com/scientificreports www.nature.com/scientificreports/ NI, independently from which model's ΔNI is taken as a reference. For example, the average value of τ is highest, and close to a value of one, when one only considers networks in the top quartile of ΔNI. We further considered whether there is increased concordance in model predictions in subsets of the nodes. To this end we focused on networks with τ < 1 and calculated τ when only the 5 or 10 most ictogenic nodes are considered. Figure 4 of the Supplementary Material shows that 〈τ〉 is higher when only the 5 nodes with largest NI are considered for comparisons involving the Bistable model. However, the value of 〈τ〉 itself is still small in this case.

Model Excitability Coupling
Since heterogeneity in NI is an important determinant of concordance between models, we sought to better understand whether discrepancies between NI distributions occurred for networks with particular topologies. Figure 5(a,b) show how ΔNI relates to the number of edges in each network and the normalized standard deviation of outdegree, for all four models in the case of 19 node networks. Figure 5(a) shows that the greatest differences in NI are observed for relatively sparse networks (low number of edges). As the networks become progressively more dense, the ictogenicity of its nodes become more homogeneous (lower ΔNI). Figure 5(b) indicates the presence of a nonlinear correlation between inhomogeneities in node degree (as quantified by the normalized standard deviation of outdegree) and inhomogeneities in NI. We quantified this correlation and found Spearman correlation coefficients between σ out /<k out > and ΔNI of 0.979 (Physiological), 0.961 (Theta) and 0.972 (Bistable) (p < 0.001 for all models). Thus, to summarize, the effect of perturbations (node removal) is dependent on the choice of model in 19-node networks, but only when NI is homogeneous. In contrast, networks with heterogeneous NI, for example those with a large difference between the minimum and maximum values of NI across nodes (as quantified by ΔNI), display good concordance in the effect of perturbations across different node dynamics. Furthermore, we find that networks with homogeneous NI distributions correspond to networks with low normalized degree variance and greater number of edges.

Discussion
In this study we systematically characterized the effects that different choices of dynamical description have on predictions made using network models of epilepsy surgery. We considered three alternative models that have been used to understand the emergence of seizures in large-scale brain networks to represent the dynamics of nodes 1,10,18,25,40,42,43 . We used the measure NI to quantify the change in epileptiform dynamics that occurs upon  www.nature.com/scientificreports www.nature.com/scientificreports/ removal of a node from the network. Thus, the NI is a prediction for the effect of epilepsy surgery. Our analysis points to the importance of heterogeneity of networks as a determinant of the importance of choice of model. By this we mean the extent to which properties of each node (for example the NI, or the degree) varies across different nodes of a network. We find that when networks are sufficiently heterogeneous, the three models we considered give equivalent rankings of nodes in terms of what happens to epileptiform dynamics when they are removed. For example, for all possible 13 3-node and 199 4-node networks, all three models produce the same ranking of nodes, in terms of their NI (Figs 2 and 3(a)). Therefore, for these small networks our results show that the choice of model has no influence on which nodes would be predicted to be targets for epilepsy surgery. As we consider larger networks, with 19 nodes, the choice of model dynamics can begin to cause deviations in predictions from each model, i.e. affecting the relative ranking of nodes in terms of their NI values. However, we find that networks with a sufficiently heterogeneous distribution of NI, i.e., where at least some nodes contribute very differently to the emergence of epileptiform dynamics, the level of agreement between predictions from different models improves considerably, yielding high values of the Kendall rank for all models (Fig. 4).
We further show that networks in which nodes differ significantly in NI (i.e. high ΔNI) display greater variability in node degree (Fig. 5). Since we demonstrated predictions from the three models are in better agreement when the NI distribution is heterogeneous, this suggests that networks that contain some highly connected nodes (compared to the average), will present an NI distribution which is robust to the choice of model used to calculate it. This result corroborates and extends the findings of Lopes et al. 42 , where the authors show that the Physiological and Theta models present good agreement in ranking nodes according to their NI values for random networks of 15, 30 and 50 nodes, and higher agreement in 50 nodes scale-free networks (which are highly heterogeneous networks).
The contribution of each individual node to the emergence of epileptiform dynamics is, in general, dependent on the size of the network. The larger the network is, the smaller contribution we expect each node to make to the emergence of epileptiform dynamics, which can lead to a more homogeneous distribution of NI. If differences in NI between nodes are very small their rank ordering could be more susceptible to noise, and therefore less robust. This is a contributing factor to why we observe reduced concordance between predictions of models when NI is homogeneous (low ΔNI). From a practical perspective, we are predominantly interested in identifying which brain regions (nodes) are the most ictogenic for the purpose of recommending these for surgical resection. Such regions would be more readily distinguishable if the NI distribution is sufficiently heterogeneous. Our results suggest that predictions from the three considered models are concordant when there are nodes that are clearly more ictogenic than others within the network. Recent studies demonstrate that in fact both functional 42 and structural 58 brain networks derived from epilepsy patients display such network heterogeneity (as demonstrated, for example, by the presence of hub nodes), thus suggesting that predictions from the three models should be concordant when applied to real brain networks.
It is important to point out that the distribution of NI for any given model emerges from an interplay between the node dynamics and the network structure. One can envisage extremes in which the network structure plays a very large or very small role in the emergent dynamics. Networks with homogeneous connectivity (i.e. low standard deviation in degree, as shown in Fig. 5) present a case in which the connectivity of the network does little to facilitate differences in the dynamics of different nodes. Therefore, the distribution of NI is homogeneous (as shown in Fig. 5) and the largest determinant of the dynamics of the network is the choice of model for the nodes. On the other hand, networks with heterogeneous connectivity are examples for which the network structure itself places constraints on the emergent dynamics of nodes since, for example, there will be nodes that can heavily influence the dynamics of other nodes due to the nature of their connections. In this case, the network topology, rather than the choice of model, plays the biggest role in determining the dynamics of the network.
We also show that ΔNI decreases as the number of edges in networks increases. This would suggest that a greater level of agreement between predictions of models is obtained for sparse networks. Networks with fewer edges have heterogeneous, rather than homogeneous, connectivity (the extreme case of homogeneous connectivity has full connectivity, i.e. the maximum number of edges). In clinical applications that have used functional connectivity networks, thresholding or surrogate methods are typically used to focus on the most significant functional connections 25,59,60 , which can yield sparse networks. This is mainly justified to avoid spurious connections due to indirect correlations or random effects. Hence, we might expect large-scale brain networks derived from clinical data to yield predictions that are robust to the choice of model.
In terms of the methodology we employed, an advance in the current work is the more comprehensive approach we used to calculate BNI and NI, where changes in dynamics due to both changing node excitability and network connectivity are taken into account. Such an approach avoids setting specific arbitrary choices for parameters. Nevertheless, this method depends on the choice of boundaries for these parameters. It is impossible to define a priori an optimal parameter window in which heterogeneity in the NI distribution is maximized for all models, and a systematic search for such a window would be an extremely demanding numerical task. Our strategy was to balance the definition of boundaries in order to include all the dynamical changes of interest (i.e. both "healthy" and epileptiform dynamics). In the future, machine learning algorithms and other advanced statistical tools could be used to search for windows that optimize NI heterogeneity 50 . In addition, it is interesting to note that, in spite of the Bistable model being analyzed using a different definition of BNI (see Methods), the NI distributions closely matched those obtained with the Physiological and Theta models. This shows that the obtained results are not only robust across different models, but also across alternative definitions of ictogenicity.
Despite various attempts to derive mathematical models for the generation of seizures 1,18,29,38 there is still a lack of understanding regarding which dynamical description should be used to represent node dynamics in models of ictogenic networks. Indeed, it is likely that different dynamic descriptions are appropriate for different patients. In order to develop a deep understanding about the fundamental mechanisms generating the dynamics observed in real brain networks, and also to optimize these methods for individuals, the understanding of how (2019) 9:7351 | https://doi.org/10.1038/s41598-019-43871-7 www.nature.com/scientificreports www.nature.com/scientificreports/ alternative descriptions compare is paramount. The present work addresses this by clarifying under what circumstances alternative dynamical models yield different predictions for the effect of removing nodes from a network, as quantified by NI.
Comparisons between different models of node dynamics have been performed in other contexts. Messe et al. 61 compared seven different computational models, showing that the simplest model considered (simultaneous autoregressive model) performed better than the other more sophisticated dynamical descriptions in predicting functional connectivity from structural connectivity using DWI and MRI. In the wider context of systems biology, many efforts have been recently undertaken to simplify complex models of biochemical reaction networks 62,63 . Santolini and Barabasi 64 analyzed 87 biological models and showed that "Dynamics-Agnostic Network Models", a framework based exclusively on network topology, can provide 65-80% accuracy in predicting the impact of specific perturbation patterns when compared to the complete biochemical model. This result highlights the importance of understanding the relationship between network topology and model dynamics. In our study we focussed on node-removal perturbations and demonstrated that whether choice of model is important depends on network topology. In practice, if a network under consideration is sufficiently heterogeneous, our analysis shows that we could use any of the three models studied, and the criterion used to choose between the three could be computational efficiency. Thus the lower dimensional theta or bistable models would be preferable.
Furthermore, a comprehensive understanding of the effects of alternative dynamical descriptions on the spontaneous activity of large-scale networks can help provide additional insights into other more general modelling frameworks used to study responses to perturbations. Computational modelling techniques have been used to optimize targeting in Deep Brain Stimulation (DBS) 65 , to explain changes in brain rhythms induced by Transcranial Magnetic Stimulation (TMS) 66 , and to estimate the influence of electrode displacement on Transcranial Direct Current Stimulation (tDCS) 67 , to mention a few examples. These perturbation techniques are becoming increasingly popular in the treatment of several neurological and neuropsychological disorders, including epilepsy [68][69][70] . We hope the present work can serve to highlight the importance of understanding the influence of using alternative dynamical descriptions in predictive modelling, and can serve as a basis for this type of systematic analysis.