The interaction between neural populations: Additive versus diffusive coupling

Models of networks of populations of neurons commonly assume that the interactions between neural populations are via additive or diffusive coupling. When using the additive coupling, a population’s activity is affected by the sum of the activities of neighbouring populations. In contrast, when using the diffusive coupling a neural population is affected by the sum of the differences between its activity and the activity of its neighbours. These two coupling functions have been used interchangeably for similar applications. Here, we show that the choice of coupling can lead to strikingly different brain network dynamics. We focus on a model of seizure transitions that has been used both with additive and diffusive coupling in the literature. We consider networks with two and three nodes, and large random and scale-free networks with 64 nodes. We further assess functional networks inferred from magnetoencephalography (MEG) from people with epilepsy and healthy controls. To characterize the seizure dynamics on these networks, we use the escape time, the brain network ictogenicity (BNI) and the node ictogenicity (NI), which are measures of the network’s global and local ability to generate seizures. Our main result is that the level of ictogenicity of a network is strongly dependent on the coupling function. We find that people with epilepsy have higher additive BNI than controls, as hypothesized, while the diffusive BNI provides the opposite result. Moreover, individual nodes that are more likely to drive seizures with one type of coupling are more likely to prevent seizures with the other coupling function. Our results on the MEG networks and evidence from the literature suggest that the additive coupling may be a better modelling choice than the diffusive coupling, at least for BNI and NI studies. Thus, we highlight the need to motivate and validate the choice of coupling in future studies. Author summary Most models of brain dynamics assume that distinct brain regions interact in either an additive or a diffusive way. With additive coupling, each brain region sums incoming signals. In contrast, with diffusive coupling, each region sums the differences between its own signal and incoming signals. Although they are different, these two couplings have been used for very similar applications, particularly within models of epilepsy. Here we assessed the effect of this choice on seizure behaviour. Using a model of seizures and both artificial and real brain networks, we showed that the coupling choice can lead to very different seizure dynamics. We found that networks that are more prone to seizures using one coupling, are less prone to them using the other. Likewise, individual brain regions that are more likely to drive seizures when using additive coupling, are more likely to prevent them when using diffusive coupling. Using real brain networks, we found that the additive coupling predicted higher seizure propensity in people with epilepsy compared to healthy controls, whereas the diffusive coupling did not. Our results highlight the need to justify the choice of coupling used and show that the additive coupling may be a better option in some applications.

Modelling large-scale brain activity is key to better understanding macroscopic brain 2 dynamics [1]. Merging such models and experimental data enables posing and testing 3 hypotheses about brain function and dysfunction [1,2]. There are two main classes of 4 large-scale brain dynamic models, neural field models and brain network models 5 (BNM) [1]. Neural field models treat the cortex as a continuous medium, whereas 6 BNMs discretise the cortex into nodes. A node may typically represent a local 7 population of excitatory and inhibitory neurons, whose activity may be modelled using 8 a neural mass model such as the Wilson-Cowan model [3]. An ensemble of such coupled 9 nodes is a BNM. The BNMs are particularly suited to study the role of brain network 10 connectivity in shaping healthy and pathological dynamics as they can readily 11 incorporate a brain connectome into a brain dynamics modelling framework [4]. For 12 example, Hansen et al. [5] used a BNM and structural brain connectivity to simulate 13 functional connectivity dynamics. Goodfellow et al. [6] used a BNM and functional 14 connectivity derived from intracranial electroencephalography (EEG) to simulate and 15 predict the outcome of epilepsy surgery. Demirtaş et al. [7] used a BNM to investigate 16 the mechanisms responsible for connectivity changes in Alzheimer's disease. Given the 17 potential of BNMs to bring mechanistic understanding into the field of network 18 neuroscience, it is important to be aware of its assumptions and choices.
(NI) [6,14]. The escape time quantifies the average time taken to transit from a resting 53 state to a seizure state; the BNI measures the likelihood of a network to generate 54 seizures; and the NI quantifies the contribution of single nodes to the network's seizure 55 propensity. We first apply all three measures to quantify the behaviour of artificial 56 networks consisting of two, three and 64 nodes with additive and diffusive coupling. We 57 then test whether the two couplings provide similar results in terms of BNI when 58 applied to functional brain networks inferred from resting-state 59 magnetoencephalography (MEG) with the aim of distinguishing individuals with 60 juvenile myoclonic epilepsy (JME) and healthy controls. Concordant results in terms of 61 escape time, BNI and NI between models using additive and diffusive coupling would 62 imply that the choice of coupling is inconsequential with little impact on the predictions 63 relevant to epilepsy; whereas discordant results would ask for careful consideration when 64 choosing the coupling. To assess seizure-like dynamics with a BNM using both additive and diffusive coupling, 68 we consider a commonly used phenomenological model of seizure transitions that is 69 based on the normal form of the subcritical Hopf bifurcation [11-13, 16, 17, 19]. In this 70 model, each network node can be represented by a noisy bi-stable oscillator, where a 71 stable fixed point coexists with a stable limit cycle. Fluctuations around the fixed point 72 represent resting dynamics, whereas large oscillations around the limit cycle correspond 73 to seizure dynamics. Transitions between the two states are driven by noise and the 74 influence of other nodes in the network. 75 The network dynamics is described by the following system of stochastic differential 76 equations: where z k (t) is a complex variable that describes the dynamics of node k 78 (k = 1, 2, . . . , N , and N is the number of nodes). The function f (z) that defines the 79 activity of a single node is, The parameter ν controls the stability of the node and ω defines the frequency of the 81 oscillations that the node may display depending on ν. At ν < 0, the origin is an 82 unstable fixed point and the node oscillates around a stable limit cycle. At ν = 0, the 83 unstable point becomes stable in a subcritical Hopf bifurcation. For 0 < ν < 1, the node 84 is bi-stable with a stable fixed point at the origin and a stable limit cycle separated by 85 an unstable limit cycle. At ν = 1, the stable and unstable limit cycles meet each other 86 in a saddle-node bifurcation. In this paper, we use ν = 0.2 and we fix ω = 20 in line 87 with previous studies [11,17,19,20] (except in the two-node networks, where we use 88 ω = 0, as described below). Each node has an independent (identically distributed) 89 complex white noise process W k (t) ∈ C, the strength of which is governed by the noise 90 amplitude α > 0. 91 The interaction between nodes is determined by the adjacency matrix A jk and the coupling function g(z k , z j ). Here we use three coupling functions: g a (z k , z j ) = γz j (4) g m (z k , z j ) = β(z j − z k ) + γz j (5) corresponding to diffusive coupling (3), additive coupling (4), and a combination of the 92 two which we will call 'mixed' coupling (5). Note that if β = 0, then 93 g m (z k , z j ) = g a (z k , z j ) = γz j , whereas if γ = 0, then 94 g m (z k , z j ) = g d (z k , z j ) = β(z j − z k ). The diffusive coupling was used for example by 95 Benjamin et al. [11], Terry et al. [21], Hebbink et al. [16], Creaser et al. [18,19], and 96 Junges et al. [17]. The additive coupling was used for example by Petkov et al. [12], 97 Sinha et al. [13], and Junges et al. [20]. To the best of our knowledge, the mixed 98 coupling has never been considered. 99 Quantifiers of seizure transitions 100 To quantify and compare the effect of the three coupling functions on the behaviour of 101 the BNM we use the following measures. The first is based on escape time theory, and 102 has previously been used to classify the behaviour of motif networks of this BNM with 103 diffusive coupling [11,18,19]. The second is brain network ictogenicity which is suitable 104 for comparing the propensity of different networks to generate seizure dynamics 105 (ictogenicity) [6,12,14,20]. The third is node ictogenicity, a quantity that assesses the 106 contribution of each node to the network's ictogenicity [6,14,20].

107
Escape time 108 We will first characterise the transition to seizure dynamics in two coupled nodes. To 109 this end, we consider the mean time taken for nodes to transition from the stable fixed 110 point to the stable oscillatory seizure state. For one node, we define the escape time λ 111 as the moment t at which the amplitude of its activity z crosses a given threshold. happen on a much longer timescale. The escape time λ is a random variable and so we 118 define the mean escape time as T = E[λ]. The mean escape time of a single node (1) for 119 k = 1 has been fully characterised using Eyring-Kramers' escape time theory in [11,19]. 120 For two coupled nodes we define the first escape time as the time the first node 121 transitions to the oscillatory state, which can be either node. The second escape time is 122 then the time that it takes the other node to transition after the first one has escaped. 123 As before these are random variables and so we can define the mean first escape time 124 and mean second escape time; here we will refer to these means as the first escape and 125 November 25, 2021 4/33 second escape. The first and second escape times in a two node system with diffusive 126 coupling were characterised in [19]. Here, we will extend this work to compare the effect 127 of all three coupling functions on the escape times. 128 We identify qualitatively different regimes of escape time behaviour. By converting 129 the BNM with two nodes into polar coordinates z k (t) = R k (t) exp[ıθ k (t)] using Itô's 130 formula, we can restrict our attention to the amplitude dynamics of the system. Note 131 that the oscillatory phase of the periodic orbit does not affect the escape times as we 132 consider the case where the phase difference between nodes is zero. The steady states of 133 the amplitude dynamics are the "transition states" where either none, one or both 134 nodes have escaped. The number and stability of these states change with the strength 135 of the coupling. We perform bifurcation analysis on the transition states as we vary 136 each of the coupling strength parameters β and γ, using specialist continuation software 137 AUTO-07P [22]. We consider here only the symmetrically coupled case (bi-directional 138 coupling) for which we can consider the amplitude dynamics as a potential system and 139 identify the potential landscape V [19]. 140 We numerically compute the mean escape times for two bidirectionally coupled 141 nodes with each of the three coupling functions using custom code written in MATLAB. 142 For each coupling strength we compute 1000 simulations of the two node model using 143 the stochastic Euler-Maruyama method with step size h = 10 −3 and initial conditions 144 z k (0) = 0. To leading order the escape times do not depend on the choice of escape 145 threshold provided it lies beyond the unstable limit cycles, and in the following we fix it 146 to be a node amplitude of 0.5 in line with [11,19]. For each simulation we identify the 147 first and second escape time, then take the means over the 1000 simulations. The 148 escape times do not depend on ω and so we set this to zero in our simulations.

149
Brain network ictogenicity 150 To characterise seizure-like dynamics in networks, we use the concept of brain network 151 ictogenicity (BNI) [6,12,14,20]. The BNI quantifies the likelihood of a network to 152 generate seizures and corresponds to the average time that each network node spends in 153 the seizure state. Since the initial conditions are such that all nodes start in the resting 154 state and once they transition to the seizure state, they do not return to the resting 155 state, then the BNI is formulated as where M is a sufficiently long simulation time and λ k is the escape time of the k th node. 157 The higher the fraction λ k /M is, the longer the node k takes to seize. To quantify the contribution of each node to the network's ability to generate seizures, 169 we use the concept of node ictogenicity (NI) [6,14,23]. The NI (k) measures the relative 170 difference in BNI upon removing node k from a network: where BNI pre is the BNI prior to node removal, and BNI (k) post is the BNI after the 172 removal of node k. Note that to compute BNI (k) post , the coupling term in (1) is not 173 normalised by N , but by N + 1 (i.e., the size of the network before node removal), so 174 that the effective coupling strength is kept fixed. As in previous studies, we set the 175 coupling strength parameters such that BNI pre = 0.5 [6,14,23] (except in the three node 176 networks, where we use the network BNI computed as above), and the same parameters 177 are used to compute BNI  To compare different NI distributions obtained using the different coupling functions, 185 we consider the weighted Kendall's rank correlation τ [14,[23][24][25], where P is the number of pairs of nodes with the same order in two rankings (e.g., pairs 187 of NI values using the additive coupling that are ordered in the same way as pairs of NI 188 values using the diffusive coupling), and Q is the number of pairs in reverse order. To with BNMs applied to functional brain networks [6,13,14]. We considered random and 202 scale-free networks, both directed and undirected [14,26,27]. We generated random 203 networks using the Brain Connectivity Toolbox [28]. To build undirected scale-free 204 networks with degree distribution P (x) ∝ x −a , we used the static model [29] and a = 3. 205 Finally, we employed the Barabási-Albert algorithm to construct directed scale-free 206 networks [30]. We considered networks with mean degree c = 4 and c = 8. In the case of 207 directed networks, we used mean in-degree c in equal to the mean out-degree c out , 208 c in = c out = c. We discarded networks with disconnected components and analysed 10 209 network realisations per network topology. Thus, we studied 80 networks.

MEG functional networks 211
To verify how our results generalise to real-world brain networks we compared additive 212 and diffusive coupling in terms of BNI on MEG functional networks. We used 213 resting-state MEG functional networks from people with JME and healthy controls.

214
This MEG dataset was previously used to demonstrate that the BNI framework was 215 capable of differentiating the two groups of individuals [31]. Subsequently, the underlying sources were inferred using a linear constrained minimum 234 variance (LCMV) beamformer on a 6-mm template with a local-spheres forward model 235 in Fieldtrip [32]. The source signals were then mapped into the 90 brain regions of the 236 Automated Anatomical Label (AAL) atlas [33].

237
To obtain MEG functional networks, the source reconstructed data were divided into 238 10, non-overlapping, 20 s segments. A functional network was computed from each 239 segment using a surrogate-corrected amplitude envelope correlation (AEC) with 240 orthogonalised signals [31,33]. Thus, we considered 10 MEG functional networks per  Due to the symmetry of the system the diagrams plotted against R 2 are identical. For 253 each coupling type the transition states undergo saddle node and pitchfork bifurcations 254 as the coupling strength increases. We follow these bifurcation points in both β and γ. 255 The resulting two dimensional bifurcation diagram shows that these bifurcations   (5), when the coupling is weak, γ, β < SN 1 , the system behaves as if 289 uncoupled and neither type of coupling dominates. When γ is small (γ < 10 −1 ), the 290 first escape times follow the same pattern as for the diffusive only coupling and for large 291 β the diffusive coupling dominates. Conversely, when β is small (β < 10 −1 ), the first 292 escape times follow the same pattern as for the additive only coupling and for large γ 293 the additive coupling dominates.
Mean times for the first node to escape (First Escape, row 1) and second node to escape (Second Escape, row 2). Column (a) shows mean escape times against β for different values of γ; Column (b) shows mean escape times against γ for various β values. A legend is given for each column. Column (c) plots the mean escape times on the (γ, β)-plane with the bifurcation curves overlaid.
The diffusive dominated coupling is characterised by the area of very large escape 295 times in yellow for high β and low γ. However, for large β and γ > SN 2 the first escape 296 time no longer depends on β and additive coupling dominates. This is illustrated by the 297 almost flat lines for γ = 0.2 and 1 in panel (a1) and the coalescence of all the lines in 298 panel (b1) for γ > 0.1. The mean second escape times show a decreasing trend with 299 increasing coupling strength for both the additive and diffusive only cases. We note that 300 around the pitchfork bifurcation (PF) in the diffusive case (large β) the second escape 301 times become noise dominated.

302
Taken together, these escape times show that as coupling strength is increased the 303 state where only one node has escaped disappears and the behaviour of the nodes 304 synchronises. In other words, as soon as one node escapes the other immediately follows. 305 This is because when either of the coupling strengths are large the input from connected 306 nodes dominates the dynamics. The key difference is the time that it takes the first 307 node to escape, which is fundamentally different depending on whether additive or 308 diffusive coupling dominates the system.

309
Three-node networks 310 To consider the effect of network structure (topology) on the escape times of the 311 network and its nodes we consider the BNI and NI of three-node motifs. For simplicity, 312 in this section we compare only additive and diffusive coupling. Figure 3 shows the BNI 313 computed via (6) and NI computed via (7) for all non-isomorphic three-node networks. 314 For each network, we observe that the BNI is higher for the additive coupling than for 315 the diffusive coupling. Furthermore, we observe that networks with higher number of 316 connections tend to have higher BNI for the additive coupling, but lower BNI for the 317 diffusive coupling. These results show that, networks with more additive connections 318 are self excitatory and have more nodes with low escape times, whereas networks with 319 more diffusive connections are more self inhibitory and have more nodes with low escape 320 times. It also indicates that increasing the number of connections has a similar effect as 321 increasing the coupling strength. Thus, the discrepancy in BNI between the two 322 couplings tends to be greater in networks with more connections. 323 Figure 3(c) shows that whilst the NI is generally positive for the additive coupling, it 324 is usually negative for the diffusive coupling. This implies that node activities drive 325 seizures in the network if the coupling is additive, but tend to prevent seizures if the 326 coupling is diffusive. The higher the number of connections of the node, the stronger its 327 ability of driving (preventing) seizures if the coupling is additive (diffusive). We note 328 that as expected for symmetric networks, where the removal of one node is topologically 329 equivalent to the removal of one of the two other nodes, the NI for each node is the 330 same, i.e., within error bars (see e.g. network 8). Overall, we note that the NI 331 distributions are very different, in some cases opposite, for each type of coupling.

332
Moreover, we show the node with the highest NI is different for each network depending 333 on the coupling function.

334
Networks 335 We now turn our attention to larger networks. In this section we compare the role of the different coupling functions on the transient dynamics of the bi-stable model in networks with 64 nodes using the concepts of BNI and NI. To aid our interpretation of the following numerical results we note the relationship between the coupling function and node degree. The mixed coupling function (5) can be rewritten as  where d k is the in-degree of node k. For simplicity of interpretation, we assume that 336 z k (t) is a positive variable, such as the amplitude R k (or as the node output in the 337 theta model [14] described in the supplementary material). As above we consider only β 338 and γ positive. With this set up, the first term in (9) promotes the increase of activity 339 z k , whereas the second term may only suppress it. Therefore, the additive coupling 340 models excitation, whereas the diffusive coupling may model both excitation or 341 inhibition depending on node activities and network structure.

342
Substituting (9) into (1) gives (10) This shows that inhibition is node dependent, being proportional to a node's own 344 activity and number of in-connections. This effect will lead to very different node 345 behaviour in networks where the in-degree is highly heterogeneous between nodes.

346
The difference between additive and diffusive coupling is particularly distinct in 347 all-to-all networks (A jk = 1 for all j = k). In these networks, all nodes are topologically 348 equivalent with in-degree d k = (N − 1), and if the network is sufficiently large, then 349 their activity is on average the same. As a consequence, before any node escapes, the , is approximately zero. In contrast, the additive 351 term, γ j =k A jk z j is approximately equal to γ(N − 1)z k . This suggests that the 352 diffusive coupling tends to have a weak influence on the dynamics of resting 353 well-connected networks, whereas the additive coupling term tends to be stronger as the 354 number of connections increase.

355
Below we consider large random and scale-free networks, both directed and 356 undirected. Whilst random networks are fairly homogeneous with regards to the degree 357 distribution, scale-free networks are highly heterogeneous, having some highly connected 358 nodes [26].

359
Brain network ictogenicity 360 We first focus on the BNI, the network's propensity to generate seizure activity, across 361 different coupling functions and network structures. Figure 4 shows the BNI as a 362 function of the coupling strength γ when using the additive coupling (4). We chose a 363 range of γ such that we could observe the greatest overall possible variation in BNI, and 364 considered five levels of noise α to show its impact on BNI. We observe that in all 365 considered network structures, the BNI grows monotonically with γ. This result means 366 that in all these types of networks, the stronger the connection strength between nodes, 367 the more likely the network is to generate seizure activity. This result is in agreement 368 with our observations in the two-nodes motifs where the first and second escape times 369 decrease as γ increases (see Fig. 2), as well as the BNI results in the three-node 370 networks (see Fig. 3). We also observe that the higher α is, the higher the BNI is. Both 371 the coupling and noise terms are positive, and at larger values the nodes are more likely 372 to transit to the seizure state, hence increasing the BNI. The BNI dependence with γ is 373 qualitatively similar across the four types of network topologies, although we note that 374 in directed scale-free networks the growth in BNI is not as steep as in the other  This result is to be expected: higher mean degree implies on average stronger influence 384 from the coupling term and so the BNI reaches it maximum at lower values of γ.  These results suggest that increasing the mean degree is effectively similar to increasing 416 the coupling strength (regardless of the type of coupling). BNI as a function of diffusive coupling strength β. Each color corresponds to a different level of noise α: green is α = 0.01, orange is α = 0.03, purple is α = 0.05, brown is α = 0.1, and pink is α = 1. All other parameters are the same as in Fig. 4. Figure 6 shows the BNI as a function of the coupling strength γ (with fixed (β) and 418 β (with fixed γ) when using the mixed coupling with each of the four network types. As 419 in the additive coupling case, the BNI grows monotonically with increasing γ, and, as in 420 some of the diffusive coupling cases, the BNI decreases monotonically with increasing β 421 in all of the networks considered. These curves can be considered as cross sections of a 422 generalised BNI surface in a γ − β plane (or a hypersurface if we consider a range of α 423 values). Such surface would contain the curves observed in Figs. 4, 5, and 6 for a given 424 network and a fixed α. Interestingly, Fig. 6(b) suggests that a sufficiently strong 425 additive component γ can prevent the diffusive coupling component of suppressing the 426 BNI. Note that at α = 0.05, the diffusive coupling was capable of reducing the BNI to 427 zero in the case of the undirected networks, see the purple curves in Figs. 5(a) and (c). 428 We also observe higher agreement between the BNI curves across different network BNI as a function of (a) γ(β = 10) and (b) β(γ = 5) using the mixed coupling. Each color corresponds to a different network topology: blue corresponds to undirected random networks, red to directed random networks, green to undirected scale-free networks (overlaps the blue), and orange to directed scale-free networks. Each curve corresponds to a different network realisation. We fix α = 0.05 and mean degree c = 4.

Node ictogenicity 434
The NI quantifies the contribution of each node to the overall network's propensity to 435 generate seizures. Identifying the nodes with the highest contribution (i.e., highest NI) 436 can be useful to inform epilepsy surgery [6,14,34]. Figure 7 shows representative NI 437 distributions for each of the four network topologies, for a given level of noise. These 438 representative NI distributions show how the NI values depend on the coupling function. 439 We observe that the NI when using the additive coupling has a much larger range, with the additive coupling, meaning that the removal of nodes contributes to an overall 444 reduction in BNI, whereas the majority of nodes with diffusive coupling have negative 445 NI values, which implies that removing nodes can increase the network's ability to 446 generate seizure activity. We find that for the diffusive coupling nodes with the highest 447 degree are more likely to have the lowest NI. From (9) we observe that inhibition is 448 node dependent, being proportional to a node's own activity and number of 449 in-connections. In contrast, the NI is highest in highly connected nodes when using the 450 additive coupling because such nodes are more likely to be both excited into the seizure 451 state and to be capable of exciting their neighbours.  To better compare NI distributions from the different coupling functions, we 458 computed the weighted Kendall correlation rank τ . Figure 8 shows that the NI 459 distributions from additive and diffusive couplings are not just different, they actually 460 tend to rank the nodes in opposite order. Nodes with the highest NI in the additive Representative NI distributions of (a) undirected random, (b) directed random, (c) undirected scale-free, and (d) directed scale-free networks using the different couplings. The blue squares represent the NI computed using the additive coupling, the red triangles correspond to the diffusive coupling and the green circles correspond to the mixed coupling. The nodes were sorted such that the NI grows monotonically for the additive coupling. The error bars represent the standard error across 1000 realisations. We used α = 0.005 for the additive coupling, α = 0.03 for the diffusive coupling, and α = 0.01 for the mixed coupling. In all three coupling cases, parameters γ and β were chosen such that BNI pre = 0.5; for the mixed coupling, we fixed β = 10 and chose γ. All networks had mean degree c = 4.
coupling are likely to be the nodes with the lowest NI in the diffusive coupling, and vice 462 versa, as observed in the three-node networks (see Fig. 3(c)). Additionally, as observed 463 in Fig. 7, the NI orderings of the additive and mixed couplings are in almost perfect showing that the NI is related to the number of connections that each node has. We 474 find a positive correlation between NI and node degree in the additive and mixed 475 couplings (see S3 Fig and S5 Fig),  presumably range between the two extremes depending on parameters.

BNI of MEG functional networks 480
To further assess the impact of choosing either additive or diffusive coupling in studies 481 that aim to investigate the emergence of seizures on real-world brain networks, we 482 computed the BNI of 26 MEG functional networks from people with JME and 26 from 483 healthy controls using both additive and diffusive coupling. Figure 9 shows that the 484 BNI based on additive and diffusive couplings rank the individuals in reverse order.

485
Individuals with the highest 'additive BNI' have the lowest 'diffusive BNI', and vice 486 versa. Furthermore, while most individuals have similar diffusive BNI, they are well 487 distinguished in terms of additive BNI. Finally, we observe that people with JME have 488 on average higher additive BNI than controls (Mann-Whitney U test, p = 0.0062), but 489 lower diffusive BNI than controls (p = 0.0026). It is important to note that higher BNI 490 in people with JME relative to controls is the expected, given that higher BNI is 491 assumed to characterize a brain network with a higher propensity to produce seizures. 492 Therefore, only the additive coupling provides the hypothesized BNI distinction between 493 the two groups. We performed the same comparison for other γ and β values and found 494 similar results (not shown).

496
BNMs are useful tools to understand the role of brain network structure on healthy and 497 pathological brain dynamics. These models make assumptions about how brain regions 498 behave and how they interact. Here, we investigated two main modelling choices of 499 interaction, the additive coupling and the diffusive coupling, aiming to understand their 500 role on simulated brain dynamics. We focused on a bi-stable model of seizure 501 transitions, a model that has been used in the literature with both additive and 502 diffusive couplings [11-13, 16, 17]. We further considered a mixed coupling, combining particularly at large coupling values, which explains the different dynamical behaviours 512 in terms of escape times. In three-node networks, we found that generally networks with 513 additive coupling have higher BNI than networks with diffusive coupling, and the 514 difference in BNI increases in networks with more connections. Furthermore, nodes with 515 higher NI with the additive coupling had the lowest NI with the diffusive coupling. Our 516 simulations in larger networks with 64 nodes further supported these observations.

517
Finally, we applied the BNI to MEG functional networks from people with JME and 518 healthy controls using both additive and diffusive couplings. We found that people with 519 JME had a higher BNI based on the additive coupling than the controls, and we 520 observed the opposite using the diffusive coupling.

521
Our observations that the BNI and NI are different depending on whether we use 522 additive or diffusive coupling have consequences for studies aiming to apply these 523 measures to interrogate functional networks obtained from clinical data, as illustrated 524 with our BNI results on the MEG networks. The BNI has been used to differentiate the 525 functional networks of healthy people from people with epilepsy inferred from 526 resting-state data [11,12,31]. The hypothesis is that people with epilepsy have a higher 527 enduring propensity to generate seizures than healthy individuals, and that this 528 propensity can be assessed from their resting-state functional networks. The BNI is 529 meant to assess this enduring feature of the epileptic brain and, therefore, it is 530 hypothesized that people with epilepsy have a higher BNI than healthy controls.

531
However, our findings based on artificial networks suggest that whether the BNI is 532 higher in one group relative to the other depends on the choice of coupling function.

533
Group A may have higher BNI than group B if we choose the additive coupling, 534 whereas group B may have higher BNI than group A if we choose the diffusive coupling. 535 We confirmed this expectation by computing the BNI on the MEG functional networks 536 using the two couplings. We observed that the BNI using the diffusive coupling 537 provided almost the opposite BNI ranking of individuals compared to the BNI using the 538 additive coupling. In particular, we found that only the BNI based on the additive 539 coupling distinguished the JME group with higher BNI values than the healthy group. 540 Together with previous evidence [12,31], this result suggests that the additive coupling 541 may be a better modeling choice than the diffusive coupling, for the purpose 542 characterizing ictogenic brain networks with BNI. 543 The NI has been used to model epilepsy surgery and to make predictions about the 544 epileptogenic zone [6,14,16,17,34]. Nodes with the highest NI are taken as predictors of 545 the epileptogenic zone. The fact that the NI distribution strongly depends on the 546 coupling choice implies that predictions about the epileptogenic zone also depend on 547 this choice. The nodes with the highest NI in the additive coupling are likely to be the 548 nodes with the lowest NI in the diffusive coupling. Such disagreement between additive 549 and diffusive couplings highlights the need of finding which coupling choice is most 550 appropriate to model the brain's ictogenicity. Future work may attempt to answer this 551 question by fitting the model with the mixed coupling function to electrophysiological 552 recordings using a search over the full (γ, β) parameter space. A more detailed study 553 may even consider weighted networks in which different connections are characterised by 554 different γ and β values. However, we note that most studies that have used the 555 bi-stable model or other models of ictogenicity to investigate data have used the 556 additive coupling [6,[12][13][14]34]. Their promising results suggest that the additive 557 coupling, or perhaps the mixed coupling, may be more appropriate than the diffusive 558 coupling for such investigations. 559 We highlight that our simulations based on the additive and diffusive couplings 560 provide different understandings about the role of single node dynamics and network 561 structure on ictogenicity. In the case of the additive coupling, ictogenicity can result 562 from pathological single nodes and/or pathological brain structure. In this context, 563 'pathological single nodes' are nodes that are able to generate seizure activity even in 564 isolation. Within the phenomenological model framework, such nodes have been 565 modelled by making either ν or α node specific [16,21]. With additive coupling, 566 pathological single nodes may cause seizures due to the excitatory nature of all 567 interactions. However, pathological single nodes are not necessary for a network to 568 generate seizures. The network ictogenicity may result from the network's structure.

569
For example, highly connected nodes are likely to be prone to generate seizures, which 570 in turn makes the whole network prone to seizures. Given the excitatory nature of the 571 additive coupling, the network structure can only enhance the ictogenicity of individual 572 nodes. In contrast, in the case of diffusive coupling, ictogenicity is the result of both 573 pathological nodes and pathological brain structure. A network without pathological 574 nodes, i.e., nodes capable of generating seizures in isolation, cannot generate seizures. In 575 our simulations, we had to choose a sufficiently large level of noise α so that the BNI 576 curves had a non-zero BNI at β = 0, otherwise, the BNI would be zero at all values of β. 577 On the other hand, the existence of pathological nodes does not guarantee network network structure is such that enables their pathological activity. As a consequence, our 582 results suggest that if the additive coupling is a better model of large-scale brain 583 interactions, then knowledge about network structure may be sufficient to assess brain 584 ictogenicity because pathological nodes are not necessarily required to drive seizures. 585 On the other hand, if the diffusive coupling is a better approximation of large-scale 586 brain interactions relevant for epilepsy, then knowledge about network structure is 587 insufficient to assess brain ictogenicity. Unfortunately, evaluating whether single nodes 588 are pathological remains a challenge.

589
The results presented in the main text using the bi-stable model of seizure coupling [20]. Also, the relation between ictogenicity and the number of connections 599 uncovered in our analysis is also in agreement with previous findings using both 600 additive [14] and diffusive coupling [17]. 601 We emphasise that our results should be broadly relevant for studies using BNMs biophysical models of large-scale brain interactions [36].

614
Here we compared the impact of using additive or diffusive coupling on the dynamics of 615 a BNM relevant for epilepsy. We showed that the two couplings are not interchangeable. 616 On the contrary, the dynamics on the networks are different and the predictions of node 617 and network relative ictogenicity are often opposite. We used the two coupling 618 frameworks to assess resting-state functional networks inferred from MEG from people 619 with JME and healthy controls and found opposing results in terms of network's 620 propensity to generate seizures. The additive coupling provided the hypothesized result 621 of higher ictogenic propensity on brain networks from people with JME relative to 622 networks from controls. Thus, our results and evidence from the literature suggest that 623 the additive coupling may be a better modeling choice than the diffusive coupling, at  degree is 0.69 and the average p-value is 0.040. All parameters are the same as in Fig. 7. 670 S1 Appendix. Theta model.

671
In the main text we compare additive, diffusive and mixed couplings within the 672 bi-stable model presented in the Methods. Here, we expand our comparison to a 673 different model of ictogenicity, the theta model [14,23,31].

674
The theta model can be used as a phenomenological BNM of seizure dynamics [14]. 675 The phase θ k characterises the activity at node k and it is described by the differential 676 equation, 677 θ k = (1 − cos(θ k )) + (1 + cos(θ k ))I k (t) (11) where I k (t) is an input current to the node. At I k < 0, the node is at a fixed stable 678 phase θ (s) k (a 'resting state'). At I k = 0, there is a saddle-node on invariant circle 679 (SNIC) bifurcation. At I k > 0, the node oscillates (the 'seizure state'). The current I k (t) 680 is defined by where I 0 is the excitability of node k, ξ k (t) are noisy inputs, and the third term on the right hand side is the coupling term. N is the number of nodes, A jk is the adjacency matrix, and G(θ k , θ j , θ (s) k ) is the coupling function. We consider two coupling functions: where G a corresponds to additive coupling, and G d to diffusive coupling. G a has been 682 considered in the literature [14,23,31,34], whereas G d has not.

683
The phase θ (s) k is a stable fixed point obtained fromθ k = 0 in the absence of noise 684 and node interactions (ξ (k) (t) = 0 and γ = β = 0), At I 0 > 0 there is no stable fixed point and we use θ The noise term is independent across nodes and it is Gaussian distributed with zero 687 mean and variance σ 2 . Following previous studies [14,23,31,34], we use I 0 = −1.2 and 688 σ = 0.6 in the case of the additive coupling. We consider σ = 1.44 in the case of the 689 diffusive coupling. As in the bi-stable model, we had to consider a higher level of noise 690 in the diffusive coupling than in the additive coupling because otherwise the network 691 would always be in the resting state regardless of β.

692
To compare additive and diffusive couplings using this model, we computed the BNI 693 and NI as we did in the main text for the bi-stable model. We also used the same 694 networks (see Section ) so that to enable comparisons between the theta and bi-stable 695 models. We integrated the stochastic equations (11) using the Euler-Maruyama method 696 with step size h = 10 −2 .

697
In the theta model, the BNI is defined as follows where t (k) sz is the time that node k spends in the seizure state during a total simulation 699 time M . We used M = 4 × 10 6 as in Ref. [23,31]; see Lopes et al. [14] for more details 700 on the calculation of t increasing the mean degree of the networks has a similar effect to increasing γ or β. A 705 difference between these results and those observed with the bi-stable model is that in 706 the case of the diffusive coupling in directed networks, we found a local minimum in   Briefly, the NI correlates (anti-correlates) with node degree in the additive (diffusive) 734 coupling (except when using additive coupling in directed scale-free networks). The 735 weighted Kendall correlation rank τ is close to −1 in all networks except directed 736 scale-free networks, showing that the NI from additive and diffusive couplings rank 737 nodes in reverse order. The directed scale-free networks is presumably an exception due 738 to its highly heterogeneous nature in terms of degree distribution. We speculate that 739 the role of in-degree and out-degree is different in the two coupling cases, but not the 740 'reverse', as it seems to be the case when considering undirected networks, where the 741 in-degree is equal to out-degree.