Fighting seizures with seizures: diffusion and stability in neural systems
Erik D. Fagerholm, Chayanin Tangwiriyasakul, Karl J. Friston, Inês R. Violante, Steven Williams, David W. Carmichael, Suejen Perani, Federico E. Turkheimer, Rosalyn J. Moran, Robert Leech, Mark P. Richardson
FFighting seizures with seizures: diffusion and stability in neural systems
Erik D. Fagerholm *, Chayanin Tangwiriyasakul *, Karl J. Friston , Inês R. Violante , Steven Williams , David W. Carmichael , Suejen Perani , Federico E. Turkheimer , Rosalyn J. Moran , Robert Leech , Mark P. Richardson Department of Neuroimaging, King’s College London Department of Basic and Clinical Neuroscience, King's College London Wellcome Centre for Human Neuroimaging, University College London School of Psychology, University of Surrey Developmental Neurosciences, University College London School of Biomedical Engineering & Imaging Sciences, King’s College London UCL Great Ormond Street Institute of Child Health, University College London Centre for Epilepsy, King's College Hospital
Corresponding author: [email protected] * ,§ These authors contributed equally to this work
Abstract
Seizure activity is a ubiquitous and pernicious pathophysiology that, in principle, should yield to mathematical treatments of (neuronal) ensemble dynamics—and therefore interventions on stochastic chaos. A seizure can be characterised as a deviation of neural activity from a stable dynamical regime, i.e. one in which signals fluctuate only within a limited range.
In silico treatments of neural activity are an important tool for understanding how the brain can achieve stability, as well as how pathology can lead to seizures and potential strategies for mitigating instabilities, e.g. via external stimulation. Here, we demonstrate that the (neuronal) state equation used in Dynamic Causal Modelling generalises to a Fokker-Planck formalism when propagation of neuronal activity along structural connections is considered. Using the Jacobian of this generalised state equation, we show that an initially unstable system can be rendered stable via a reduction in diffusivity (i.e., connectivity that disperses neuronal fluctuations). We show, for neural systems prone to epileptic seizures, that such a reduction can be achieved via external stimulation. Specifically, we show that this stimulation should be applied in such a way as to temporarily mirror epileptic activity in the areas adjoining an affected brain region – thus ‘fighting seizures with seizures’. We offer proof of principle using age 2 of 24 simulations based on functional neuroimaging data collected from patients with idiopathic generalised epilepsy, in which we successfully suppress pathological activity in a distinct sub-network. Our hope is that this technique can form the basis for real-time monitoring and intervention devices that are capable of suppressing or even preventing seizures in a non-invasive manner.
Keywords: ensemble dynamics, dynamic causal modelling, epilepsy, chaos control, seizure activity, Fokker-Planck.
Introduction
There is an ongoing interest in treating epilepsy by using brain stimulation
1, 2, 3 , as it allows for direct perturbations of the physiological states of neural systems
4, 5 . However, the three basic questions of when , where and how to stimulate for maximum clinical efficacy remain unanswered
6, 7, 8 . Therefore, there is a pressing need for the development of computational frameworks that can be used to model the effect of brain stimulation on neural populations and for the construction of optimal protocols for therapeutic intervention. Over the last decade, several mathematical models have been proposed to explain the emergence of seizures – primarily with focal epilepsy using electroencephalography (EEG) and electrocorticography (ECoG). For example, Benjamin et al. developed a network-based model to describe a phenomenological model of seizure initiation, while Sinha et al. and Goodfellow et al. developed models to predict neurosurgical outcome. In this study, we propose a model, tested with a combination of EEG and functional magnetic resonance imaging (fMRI) data, in patients with idiopathic generalised epilepsy (IGE). Dynamic causal modelling (DCM) provides a powerful analytical tool in this setting. DCM was originally designed for inferring latent structure from blood-oxygen-level-dependent age 3 of 24 (BOLD) time series, by optimising the parameters of a generative model, such as intrinsic connectivity and external driving inputs. In its simplest mathematical form, DCM rests upon a first-order ordinary differential equation (the neuronal state equation), which describes the ways in which neuronal signals change with respect to time. Here, we show that the neuronal state equation can be generalised to account for the ways in which signals change – not only with respect to time – but also with respect to the connectivity architecture of the network within which the signals are constrained to propagate. We show that this additional structural component conforms to a diffusion process and therefore that the generalised neuronal state equation takes the mathematical form of the Fokker-Planck equation. The latter is a partial differential equation that is used to describe the probabilistic evolution of a system, initially studied in the context of Brownian motion , and increasingly used in the modelling of neural systems
14, 15, 16 . As we will show, it is the diffusive property of the Fokker-Planck equation that facilitates the suppression of activity via gradient modulation in brain regions prone to epileptic seizures. This paper comprises three sections. In the first section, we outline the theoretical basis for subsequent applications by showing that the neuronal state equation takes the functional form of the Fokker-Planck equation, when including structural degrees of freedom (connections) within a network. We provide construct validation that the ensuing model provides a more parsimonious explanation of resting state fMRI. Specifically, we show that time series in both epilepsy patients and healthy control subjects are better modelled by the Fokker-Planck equation, compared with the classic (non-structural) neuronal state equation. In the second section, we show that an initially unstable region can be pushed into a stable regime by virtue of a reduction in network diffusivity. We demonstrate that such a reduction can be achieved by using external stimulation to mirror seizure activity in the area(s) surrounding a pathological brain region. age 4 of 24
In the third section, we report a series of simulations incorporating individualised EEG and fMRI data collected in patients with idiopathic generalised epilepsy. We show promising evidence that electrical stimulation is a viable method for suppressing epileptic seizures.
Methods
The generalised neuronal state equation:
Causal models can be considered as the evolution of states 𝑥 as a static function of themselves. For instance, for the 𝑖 !" region: 𝑥 = 𝑓(𝑥, 𝑣, 𝜃), [1] where 𝑓 is a nonlinear function, 𝑣 are external inputs, and 𝜃 are model parameters – usually interpreted in terms of connectivity or rate constants. DCM extends the assumptions in [1] by considering the ways in which states change with respect to time, so that: 𝑑𝑥 𝑑𝑡 = 𝑓(𝑥, 𝑣, 𝜃), [2] which reduces to [1] in the limit that inputs 𝑣 vary slowly relative to the states 𝑥 . Extrapolating the logical progression from [1] to [2], we can define a generalised neuronal state equation in which states may vary, not just with respect to time as in [2], but with respect to any arbitrary number of dimensions 𝜇 , such that: - 𝑑𝑥 𝑑𝜇 $$%& = 𝑓(𝑥, 𝑣, 𝜃). [3] For example, if we consider the simple three-node network in Fig 1: Figure 1: Three-node network.
The structural connection between nodes 1 and 2 is given by s and the structural connection between nodes 1 and 3 is given by s . We describe the evolution of the first node using [3] as follows: - 𝑑𝑥 & 𝑑𝜇 $’$%& = 𝑓(𝑥, 𝑣, 𝜃). [4]
12 3s s age 5 of 24 We then retain time as the first dimension (𝜇 & = 𝑡) on the left-hand side and assign the remaining two terms to the rate of change of states along the structural dimensions of the network, such that: 𝑑𝑥 & 𝑑𝑡 + 𝜎 1 𝑑𝑥 & 𝑑𝑠 &( + 𝑑𝑥 & 𝑑𝑠 &’ [5] where 𝜎 is a (rate) constant of proportionality with dimensions 𝑇 )& ; and e.g. *+ ! *, !" describes the rate of change of the first node’s activity along the edge connecting the first and second nodes. Intuitively, this kind of extension can be thought of as generalising the dynamics of any one point in the brain to a partial differential equation (PDE) that models the spatiotemporal dynamics (e.g., neuronal field models). However, here, the spatial aspect is generalised to edges for connections of a graph or network. We write the general form of [5] as follows for the 𝑖 !" node in a network of 𝑁 regions: 𝑑𝑥 𝑑𝑡 + 𝜎 - 𝑘 𝑑𝑥 𝑑𝑠 = 𝑓(𝑥, 𝑣, 𝜃), [6] where 𝑘 is the adjacency matrix element connecting the 𝑖 !" and 𝑗 !" regions, the inclusion of which ensures that only nearest neighbours in the network are taken into account. Note that in the context of neural systems, the definition of nearest neighbours may also be taken to mean regions that are structurally connected via white matter tracts. This means that any given region or node can have more neighbours than if we were modelling the cortical sheet as a two-dimensional Markov field. The structural gradients *+ *, in [6] can be discretized by expressing them in terms of the difference in state values (at a given time) at the 𝑖 !" and 𝑗 !" nodes, such that: 𝑑𝑥 𝑑𝑠 → 𝑥 $ − 𝑥 , [7] age 6 of 24 which, together with [6], tells us that: 𝑑𝑥 𝑑𝑡 = 𝑓(𝑥, 𝑣, 𝜃) + 𝜎 - 𝑘 ;𝑥 − 𝑥 $ < -$%& , [8] where we see that in the limit of zero structural gradients ;𝑥 = 𝑥 $ < , the second term on the right-hand side vanishes and [8] reduces to the original DCM neuronal state equation in [2]. In other words, every node behaves in the same way and can be described with an ordinary differential equation. In summary, we obtain a third level of neuronal state equations building on the simple causal models in [1], through the classical DCMs in [2], and ending with the generalised form in [8], with each level reducing to the former in the appropriate limiting cases. DCM and the Fokker-Planck equation:
The second term on the right-hand side of [8] can be re-written as follows: ! 𝑘 !" ! − 𝑥 " & = 𝑥 ! ! 𝑘 !" − ! 𝑘 !" 𝑥 " = 𝑥 ! 𝑑 ! − ! 𝑘 !" 𝑥 " = ! !" 𝑑 ! − 𝑘 !" &𝑥 " = ! 𝑙 !" 𝑥 " , [9] where 𝑑 is the degree of the 𝑖 !" node; 𝛿 is the Kronecker delta function, which equals unity when 𝑖 = 𝑗 , and otherwise equals zero; and 𝑙 is the graph Laplacian matrix element connecting the 𝑖 !" and 𝑗 !" nodes, where the graph Laplacian matrix 𝐿 is defined as the difference between the degree matrix 𝐷 and adjacency matrix 𝐾 , such that 𝐿 = 𝐷 − 𝐾 . Using [9] we can write [8] as follows: 𝑑𝑥 𝑑𝑡 = 𝑓(𝑥, 𝑣, 𝜃) + 𝜎 - 𝑙 𝑥 $-$%& , [10] and as the graph Laplacian is the discretized version of the Laplace operator , the second term on the right-hand side describes a diffusion process, hence lending an interpretation to 𝜎 as a diffusion coefficient. Therefore, [10] takes the form of the discretized Fokker-Planck equation, with a drift term: 𝑓(𝑥, 𝑣, 𝜃) and a diffusion term: 𝜎 ∑ 𝑙 𝑥 $-$%& . age 7 of 24 Linear stability analysis:
In order to model neural time series, we assume first order linear interactions, which means [8] can be written as: 𝑑𝑥 𝑑𝑡 = - 𝑝 𝑥 $ + 𝜎 - 𝑘 ;𝑥 − 𝑥 $ < -$%& + - 𝑞 $ 𝑣 $ + 𝜔 ( , [11] where the matrix element 𝑝 reduces to the DCM intrinsic coupling matrix element 𝑎 in the limiting case of zero structural gradients ;𝑥 = 𝑥 $ < ; the matrix element 𝑞 $ reduces to the DCM extrinsic coupling matrix element 𝑐 in the limiting case of zero structural gradients ;𝑥 = 𝑥 $ < ; and 𝜔 ( are non-Markovian fluctuations in the 𝑖 !" region’s activity . These fluctuations model deviations from the linear flow of states under an adiabatic approximation. In other words, we assume a centre manifold for the dynamics, which are linear and assign fluctuations tangential to the manifold to 𝜔 ( , which decay rapidly and return to the manifold under the centre manifold theorem . To create a DCM of observable timeseries, we can use [11] as a state space model with fast (analytic) fluctuations 𝜔 ( and map the latent states 𝑥 to observable quantities with additive observation noise. Equation [11] is used to model all the neural time series presented in this paper, using standard (variational) routines in the Statistical Parametric Mapping (SPM) software. Extrinsic coupling does not affect resilience to perturbation, as the stability of a linear time invariant (LTI) system is determined by the roots of the characteristic equation (𝑠𝐼 − 𝐽) )& , where the Jacobian 𝐽 comprises the first two terms of [11] only . Retaining these terms, we can re-write [11] as: 𝑑𝑥 𝑑𝑡 = - 𝑝 𝑥 $-$%& + 𝜎 - 𝑘 ;𝑥 − 𝑥 $ < -$%& = - J;𝑝 − 𝜎𝑘 <𝑥 $ + 𝜎𝑘 𝑥 K -$%& , [12] age 8 of 24 which we can write out explicitly for the three-node network in Figure 1 as follows: L𝑥̇ & 𝑥̇ ( 𝑥̇ ’ N = O𝑝 && + 𝜎(𝑘 &( + 𝑘 &’ ) 𝑝 &( − 𝜎𝑘 &( 𝑝 &’ − 𝜎𝑘 &’ 𝑝 (& − 𝜎𝑘 (& 𝑝 (( + 𝜎(𝑘 (& + 𝑘 (’ ) 𝑝 (’ − 𝜎𝑘 (’ 𝑝 ’& − 𝜎𝑘 ’& 𝑝 ’( − 𝜎𝑘 ’( 𝑝 ’’ + 𝜎(𝑘 ’& + 𝑘 ’( )P L𝑥 & 𝑥 ( 𝑥 ’ N. [13] The Jacobian of this generalised DCM can thus be written as follows for a network comprising 𝑁 regions: 𝐽 = ⎣⎢⎢⎢⎢⎢⎡𝑝 && + 𝜎 - 𝑘 &$-$1& ⋯ 𝑝 &- − 𝜎𝑘 &- ⋮ ⋱ ⋮𝑝 -& − 𝜎𝑘 -& ⋯ 𝑝 -- + 𝜎 - 𝑘 -$-$1- ⎦⎥⎥⎥⎥⎥⎤ . [14] Diffusion and stability:
In 1953, Turing showed that an initially stable dynamical system can be rendered unstable by virtue of a diffusion mechanism, allowing for the emergence of spatial inhomogeneities – now known as Turing patterns . Here, we proceed via similar logic, except we begin with the opposite premise and ask the following question: can we push an initially unstable system (such as a neural system prone to seizures) into a stable regime by altering the system’s diffusivity? To answer this, we multiply the diffusion coefficient 𝜎 by a constant 𝛼 : 𝜎 → 𝛼𝜎. [15] By multiplying 𝜎 by some unknown quantity 𝛼 in this way we can determine – in the subsequent linear stability analysis – whether this change increases or decreases diffusivity in order to render an initially unstable system stable. We assume that the system described by [14] is initially unstable, such that the sum of the Real components of its eigenvalues (given by the trace) is positive: 𝑡𝑟𝐽 = 𝑝 && + 𝜎 - 𝑘 &$-$1& + ⋯ + 𝑝 -- + 𝜎 - 𝑘 -$-$1- > 0, [16] age 9 of 24 We then transform [14] via [15] to obtain:
𝐽′ = ⎣⎢⎢⎢⎢⎢⎡𝑝 && + 𝛼𝜎 - 𝑘 &$-$1& ⋯ 𝑝 &- − 𝛼𝜎𝑘 &- ⋮ ⋱ ⋮𝑝 -& − 𝜎𝑘 -& ⋯ 𝑝 -- + 𝜎 - 𝑘 -$-$1- ⎦⎥⎥⎥⎥⎥⎤ , [17] which we assume has been rendered stable due to the altered diffusion coefficient, such that the Real component of the sum of its eigenvalues is now negative: 𝑡𝑟𝐽′ = 𝑝 && + 𝛼𝜎 - 𝑘 &$-$1& + ⋯ + 𝑝 -- + 𝜎 - 𝑘 -$-$1- < 0 , [18] which, together with [16], means that: 𝑡𝑟𝐽′ = 𝑡𝑟𝐽 + 𝜎(𝛼 − 1) - 𝑘 &$-$1& < 0. [19] We then know from [16] that 𝑡𝑟𝐽 > 0 . Furthermore, we know that diffusion coefficients are necessarily positive, i.e. 𝜎 > 0 , given that they play the role of rate constants . Therefore, the only way in which [19] can be satisfied is if 𝛼 is less than unity: (𝛼 − 1) - 𝑘 &$-$1& < 0 ⟹ 𝛼 < 1, [20] i.e. an initially unstable system can be rendered stable by virtue of a reduction in the diffusion coefficient. Diffusivity can be altered via external driving inputs:
In practice it is not possible to change the diffusion coefficient as in [15], due to the fact that it is an intrinsic property of the dynamical system described by [14]. However, as the diffusion coefficient quantifies diffusivity via Fick’s law , we recover an important piece of information from [20]; namely, that if we want to push a system toward stability, we must act so as to decrease diffusivity. age 10 of 24 If we look at the governing equation of motion [11], we note that the diffusion term 𝜎 ∑ 𝑘 ;𝑥 − 𝑥 $ < -$%& comprises three factors: 1) the diffusion coefficient 𝜎 , which we noted above cannot be changed; 2) the adjacency matrix element 𝑘 which, similar to 𝜎 , is intrinsic to the system and is therefore also unchangeable (without resorting to surgical intervention); and finally 3) the gradient ;𝑥 − 𝑥 $ < . It is this last factor that we can influence by applying specific external driving inputs in strategic locations, in order to decrease gradients and thus to decrease diffusivity. For instance, let us consider the same three-node system shown in Figure 1 and assume that the bottom node is prone to instabilities (Fig. 2A). Figure 2: Gradient reduction via stimulation. A)
The seizure-prone node. This displays the activity shown in the red graph therein. B) External driving input is applied in such a way as to allow the two neighbouring nodes to mirror the activity of the seizure-prone node in A). C) Extrinsic coupling to the two regions neighbouring A).
In the case of epilepsy, we therefore propose a form of intervention in which we apply external driving input (Fig. 2B) in such a way as to allow the two neighbouring nodes (Fig. 2C) to mirror the activity in the region that is prone to seizures. In summary, we have derived a generalised dynamic causal model of neuronal activity that can generate empirical timeseries. In what follows, estimate the parameters of the DCM using empirical timeseries from patients with epilepsy. Equipped with these parameters, we can then evaluate the Jacobian and simulate the effects of an intervention that moves the ensemble or population dynamics implicit in [11] from a regime of instability (i.e., seizure activity) to a one of stability.
The seizure network : In the analyses below, we use a network comprising the following regions: frontal mid, frontal mid orbital, precuneus, and thalamus. This network is known to play an important role in generating generalised spike and wave (GSW) discharges and, for simplicity, we will refer to it henceforth as the ‘seizure network’. age 11 of 24 Participants and data acquisition : We analyse data recorded from 15 patients (6 male) with juvenile myoclonic epilepsy (JME) with a mean age of 24.5 years, and 15 age-matched healthy controls (5 male) with a mean age of 25.2 years (see Supplementary Table I). The data were acquired at the Institute of Psychiatry Psychology and Neuroscience (IoPPN), King’s College London. The patients did not have any neurological diagnoses other than epilepsy and had no history of drug or alcohol misuse. The study was approved by the Riverside Research Ethics Committee (12/LO/2005) and all participants signed a written informed consent form prior to the study, according to the declaration of Helsinki (2013). All participants underwent a resting state simultaneous EEG-fMRI and a diffusion tensor imaging (DTI) session with 32 directions and b = 1500 s/mm . A 3T scanner (MR750, GE Healthcare) was used to acquire 300 echo-planer images ( mm, field of view 211 mm, repetition time 2.16 s, echo time 25 ms, flip angle 75 ° , 36 slices, slice thickness 2.5 mm). Bayesian model inversion:
The stochastic differential equation in [11] furnishes a dynamic causal model, where fast fluctuations are assumed to be small. This means that, assuming that the underlying dynamics can be modelled by [11], the model parameters 𝜃 = (𝑝, 𝜎) can be recovered from observations of BOLD signals (Fig. 3A).
Figure 3: Equation of motion. A)
BOLD signal intensity (i.) for an example time course (in seconds) in one patient (blue), together with the time series estimate (red) following Bayesian model inversion. This inversion provides posterior densities over the matrix elements 𝑝 in the drift term, and the diffusion coefficient 𝜎 in the diffusion term, informed by a non-Markovian noise process 𝜔 . B) Example of an adjacency matrix from DTI data, from which we obtain matrix elements 𝑘 in the diffusion term. C) The eight regions of the seizure network, in which the left and right frontal mid regions (red) give rise to unstable activity. Stimulation would then be applied to one or more of the remaining (black) regions in a way that mirrors the activity of the red nodes. age 12 of 24
As the parameter space associated with the eight-region seizure network is larger than can be accommodated by the time points available for each scan, we reduce the number of free parameters via functional connectivity priors
25, 26 . This allows us to constrain the optimization such that we retain ~10 × as many time points as free parameters. The adjacency matrix elements 𝑘 are not included as free parameters as they are known a priori from diffusion imaging (Fig. 3B). We use generalised or variational Bayesian filtering (specifically, Dynamic Expectation Maximisation (DEM)) to a) infer the latent states; b) estimate the parameters; and c) hyperparameters, i.e. the precision components of fluctuations on the states and observation noise. After recovering the posterior densities, hyperparameters, and variational free energies on an individual level for all 30 subjects, we then obtain two group-level models by performing Bayesian model averaging separately across 15 patients and 15 controls. These averaged models are used as the bases for all the forward generative models using external stimulation as a means of reducing diffusivity (Fig. 3C). Data preprocessing : To pre-process the fMRI data, we use the statistical parametric mapping (SPM8, r613) software running on MATLAB (R2016b), together with the FIACH package for R (3.2.2). First, we convert the data from DICOM to Nifti formats. We then delete the first four volumes of each session to avoid magnetic saturation effects. We subsequently re-align all images to the first remaining volume. To correct for possible artefacts, we apply the FIACH toolbox to the BOLD time series. We then normalise all data into standard MNI space with 2 mm isotropic voxels. All images are then spatially smoothed using a Gaussian filter of 8 mm fullwidth at half maximum. The BOLD signal is filtered between 0.04 – 0.07 Hz . This frequency range was chosen to minimise the overlap between the BOLD signal and possible breathing and pulsation artefacts. We parcellate the brain into standard 90 automated anatomical labelling (AAL) regions (excluding the cerebellum) . Finally, we apply principal component analysis (PCA) to the voxel time series within each region, from which we retain the first principal component to summarise the activity in each region . Probabilistic age 13 of 24 tractography is used to pre-process DTI data using the iFOD2 algorithm within the MRtrix software . Streamlines are filtered using SIFT , resulting in streamlines. We estimate a
90 × 90 structural connectivity matrix according to the number of streamlines connecting each pair of regions, normalised by their combined volumes. To reduce inter-subject variability, we then normalise each structural connectivity matrix relative to its maximum value. For each subject, we binarize the structural connectivity matrix according to 18 thresholds between 10% and 95% (in steps of 5%). Equipped with the structural measures 𝑘 and the summaries of regional activity in our seizure network, we estimate the latent states 𝑥 and parameters of the DCM – crucially including the diffusion coefficient 𝜎 (see equation [11]). Our special interest here lies in the diffusion coefficient and its role in mediating dynamical instability that we associate with seizure activity. Results Diffusivity is higher in the patient group:
We find that the patient group has higher diffusivity as compared with the control group across structural adjacency matrix thresholds used to define 𝑘 (see Figure 4). The higher diffusivity in the patient group across all thresholds speaks to the approach of suppressing seizure activity by decreasing diffusivity. Figure 4: Diffusion coefficients.
The diffusion coefficient (𝜎) as a function of structural (DTI) adjacency matrix threshold (%) for patients and controls following Bayesian model averaging.
Stability is lower in the patient group:
We find that the patient group is associated with lower levels of stability – as quantified with the sum of the Real components of the eigenvalues age 14 of 24 – as compared with the control group across structural adjacency matrix thresholds (see Figure 5). The lower stability in the patient group speaks to the approach of aiming to increase stability. All subsequent results shown in this paper are calculated using a threshold of 50% to define structural connectivity (i.e., 𝒌 𝒊𝒋 ). Figure 5: Stability.
The sum of the Real components of the eigenvalues of the Jacobian (∑ ℝ(𝜆)) – measuring intrinsic stability of neural dynamics – as a function of structural (DTI) adjacency matrix threshold (%) for patients and controls following Bayesian model averaging.
Accounting for diffusion improves models:
In order to determine whether we are licensed to include the diffusion term in the equation of motion [11], we first optimise the full model including a non-zero diffusion coefficient 𝜎 (i.e., including diffusion) and subsequently use Bayesian model reduction
34, 35 to estimate the evidence for the reduced model in which 𝜎 = 0 (i.e. excluding diffusion). We specify the reduced model by setting the prior variance over the 𝜎 parameter to zero, where 𝜎 is given a prior mean of zero. We show, in both patients and controls, that the variational free energies and associated probabilities are higher for the full models including diffusion, than in the reduced models excluding diffusion (Figure 6). This Bayesian model comparison provides evidence for our assumption that the diffusion term in the equation of motion [11] is necessary to explain these BOLD time series. Figure 6: Bayesian model reduction. A)
Control group. Approximate lower bound on log model evidence afforded by the free energy (F) following Bayesian model reduction for the reduced model without (w/o) diffusion and the full model with (w) diffusion. B) Probabilities derived from the log evidence in A).
C) & D)
Same layout as A) & B), but for the patient group. age 15 of 24
Ictal onset perturbation:
Following Bayesian model inversion and model averaging in the patient group, we use the resulting parameters to create an in silico seizure network, in which we can perturb the mid-frontal sources with ictal onset activity taken from EEG measurements (Figure 7A & B). We then obtain the response of the mid-frontal regions to this exogenous stimulation, which we see climbs in an uncontrolled manner, due to the associated positive Real eigenvalue (Figure 7C).
Figure 7:
Ictal onset perturbation . A)
Normalised EEG activity from frontal lobe (FP1 & FP2) activity. Ictal activity is shown by the yellow section and ictal onset activity is shown by the red section. B) The seizure network shown in MNI space (left) with the mid-frontal region (left & right) indicated by the red nodes. The mean ictal onset activity (right), corresponding to the red sections in A), is used as the external driving input and is supplied to the red nodes, as indicated by the inward-pointing red arrows. C) The mean response of the mid-frontal region to the stimulus in B), as indicated by the outward-pointing red arrows. Note that the time scale is longer due to the model inversion having been applied to BOLD data.
Seizure suppression:
Using the same setup as in Figure 7, we again run the forward model. However, this time – in addition to the ictal driving stimulus – we supply additional external stimulation to all nodes except the mid-frontal region, in a way that mirrors the ictal onset activity in Figure 7B. We show that, in agreement with our theoretical predictions, reducing activity gradients in this way results in the response of the unstable mid-frontal region being suppressed (Figure 8). age 16 of 24
Figure 8: Simulating seizure suppression. A)
External stimulation applied to all nodes except for the mid frontal region for 1000 forward models ranging from zero stimulation (yellow) to a stimulation profile that perfectly mirrors the ictal onset activity in Figure 7B (black). B) Response of the mid frontal region for the same 1000 forward models in A) with matching colours, i.e. the yellow response corresponds to zero stimulation and the black response corresponds to a stimulation profile that perfectly mirrors the ictal onset activity in Figure 7B.
Note that although these simulations are based on fMRI data collected in patients with epilepsy, these patients did not experience seizures while in the scanner: the ictal onset activity collected with EEG was recorded separately. The output of the model in Figure 7C and Figure 8B should therefore not be viewed as mimicking a seizure-like BOLD response. Instead, what we show here is that, in response to a perturbation in the form of real ictal onset activity (Figure 7B and Figure 8A), an initially unstable system with a diverging response can be suppressed – a result that we suggest could be clinically advantagous in the treatment of epilepsy.
Minimum effective stimulation : In the previous section we stimulated all nodes except for the mid-frontal region (which receives the driving input) to demonstrate proof of principle. However, for treatment purposes, it is clearly better to stimulate as few regions in the seizure network as possible, in order to remain minimally invasive. We therefore proceed by using each of the nodes of the seizure network individually in turn as the stimulation target. This allows us to determine the order of the regions, ranked in terms of the lowest stimulation strength required to suppress the mid-frontal response (Table 1). age 17 of 24
Table 1:
Ranking order in terms of the stimulation strength required to suppress the BOLD response in the mid-frontal region to the same extent as in Figure 8B by stimulating a single region in the seizure network. Stimulation strengths are presented relative to the lowest value (thalamus left), which is assigned a value of unity. N/A values are assigned if it is not possible to suppress the BOLD response in the mid-frontal region by targeting the corresponding single region.
Note that the mid-frontal region is not included in this ranking as it is being supplied with the driving input in the form of ictal onset activity and our technique requires mirroring this activity in neighbouring nodes in the network, in order to decrease gradients and thus diffusivity. Using this in silico stimulation protocol, we find that the thalamus (left) requires the lowest stimulation strength (by a considerable margin) to suppress activity in the unstable mid-frontal region relative to the other six regions in the seizure network.
Discussion
We began by showing that the neuronal state equation from Dynamic Causal Modelling takes the form of the Fokker-Planck equation when generalised to account for the propagation of neuronal activity over structured connections. As the Fokker-Planck equation entails a conservation of probability , the generalised state equation implies a conservation of neuronal activity in the brain – which can be plausibly motivated in terms of conservative aspects of neuronal message passing, such as the balance between excitation and inhibition
37, 38 . This balance lies at the base of functional modes, which are found (particularly in the cortices) to repeat across scales in the brain . The canonical computational units at the most elementary scale take on various ratios of excitatory and inhibitory neurons. However, the current consensus is that the basic unit of the cortical system is the pyramidal interneuron gamma Rank Region Stimulation 1
Thalamus (left) 1 Frontal Mid Orbital (right) 20 Precuneus (left) 76 Precuneus (right) 82 Frontal Mid Orbital (left) N/A Thalamus (right) N/A age 18 of 24 network (PING) . The PING configuration is made up of a pyramidal excitatory neuron (PN) and a fast spiking inhibitory parvalbumin interneuron (IN). One can interpret the diffusion term in equation [11] as an intrinsic mode of interaction between such neuronal units across spatial scales, where the responsible mechanism could be due to: 1) neurotransmission at the microscopic scale of individual neurons, 2) electrical impulses at the mesoscopic scale of microcircuits, and 3) long-range white matter connections between brain regions at a macroscopic scale, as captured by our DTI data. The intervention technique we are proposing is a departure from the status quo, which usually involves the opposite approach – namely, increasing inhibition
41, 42, 43 . Our method was shown to work in the context of forward generative models, parameters of which were optimised from data collected in epilepsy patients. Specifically, we applied stimulation in such a way as to decrease activity gradients between unstable regions and their neighbouring nodes, thereby decreasing diffusivity. It is important to emphasize that phasic stimulation cannot change the system’s long-term stability, as this is determined by intrinsic properties, such as the diffusion coefficient and connectivity. Instead, we are proposing that it should be possible to temporarily suppress pathological regional activity – achieving a ‘quasi-stability’ – sufficiently long enough to attenuate seizure activity. Furthermore, one need not wait for a seizure to be in progress in order to activate stimulation. Rather, it may be advantageous to continuously minimize activity gradients in the areas surrounding a known pathological region (e.g., epileptogenic zone), thereby preventing seizures from occurring in the first place. Overall, the approach proposed here provides a novel framework to address the three fundamental questions of when , where, and how to stimulate the brain in order to suppress pathological activity in the context of epilepsy. In particular, we investigated the impact of stimulation strength in relation to the number of stimulation sites and proposed specific age 19 of 24 stimulation timings and profiles, in such a way as to achieve modulation of activity gradients. As our technique relies upon supplying otherwise healthy regions of the brain with stimulation that mirrors pathological activity, it is clearly clinically advantageous to target as few regions as possible. It is for this reason that we focused on targeting a single region and found that the thalamus required the lowest stimulation strength, in line with the broad literature showing the thalamus to be a key region in seizure generation
44, 45, 46, 47 . However, it is in principle possible to target multiple nodes. The trade-off between the number of target nodes and stimulation strengths required is to be determined going forward in practical applications of this technique on a patient-by-patient basis – informed by the specific pathology of the individual. Stimulation techniques are currently available using both invasive and non-invasive methods. The fast computational processing times associated with our strategy render it compatible with closed-loop approaches, which are increasingly seen as providing the greatest clinical efficacy in delivering personalised therapy . Furthermore, the methodology presented is not limited to applications in epilepsy. For instance, it may be beneficial to use the same gradient reduction technique in the treatment of e.g. disorders associated with cortical spreading depression (CSD) . There are similarities between our approach and coordinated reset strategies, given that our results support targeting several stimulation sites in a spatially and temporally coordinated manner
50, 51 . However, a critical difference in our study is that we demonstrate that abnormal activity can be mitigated by changes in network diffusivity via the counter-intuitive approach of increasing excitation in a strategic manner, rather than by following a phase-resetting mechanism. Specifically, we demonstrated that unstable activity can be suppressed by modulating the neighbourhood of affected brain regions, with the stimulation profile and timing age 20 of 24 chosen in such a way as to mirror the pathological activity – hence ‘fighting seizures with seizures’.
Supplementary Table I:
Demographics of the subjects in this study
References
1. Nitsche MA, Paulus W. Noninvasive brain stimulation protocols in the treatment of epilepsy: current state and perspectives.
Neurotherapeutics , 244-250 (2009). 2. Salanova V. Deep brain stimulation for epilepsy. Epilepsy Behav , 21-24 (2018). 3. Vicario CM, Nitsche MA. Non-invasive brain stimulation for the treatment of brain diseases in childhood and adolescence: state of the art, current limits and future challenges.
Front Syst Neurosci , 94 (2013). 4. Lozano AM , et al. Deep brain stimulation: current challenges and future directions.
Nat Rev Neurol , 148-160 (2019). 5. Polania R, Nitsche MA, Ruff CC. Studying and modifying brain function with non-invasive brain stimulation. Nat Neurosci , 174-187 (2018). 6. Lorenz R , et al. Efficiently searching through large tACS parameter spaces using closed-loop Bayesian optimization.
Brain Stimul , 1484-1489 (2019). 7. Muldoon SF , et al. Stimulation-Based Control of Dynamic Brain Networks.
PLoS Comput Biol , e1005076 (2016). Healthy controls Patients Subject Age Gender Subject Age Gender
H01 24 M P01 16 F H02 22 M P02 15 F H03 24 F P03 17 F H04 20 F P04 26 F H05 26 M P05 16 F H06 25 F P06 20 M H07 23 M P07 25 M H08 23 F P08 37 F H09 22 F P09 28 F H10 20 F P10 30 F H11 28 F P11 35 M H12 17 F P12 21 M H13 44 F P13 22 F H14 39 F P14 39 M H15 21 M P15 20 M age 21 of 24
8. Stiso J , et al.
White Matter Network Architecture Guides Direct Electrical Stimulation through Optimal State Transitions.
Cell Rep , 2554-2566 e2557 (2019). 9. Benjamin O , et al. A phenomenological model of seizure initiation suggests network structure may explain seizure frequency in idiopathic generalised epilepsy.
J Math Neurosci , 1 (2012). 10. Sinha N , et al. Predicting neurosurgical outcomes in focal epilepsy patients using computational modelling.
Brain , 319-332 (2016). 11. Goodfellow M, Rummel C, Abela E, Richardson MP, Schindler K, Terry JR. Estimation of brain network ictogenicity predicts outcome from epilepsy surgery.
Sci Rep , 29215 (2016). 12. Friston KJ, Harrison L, Penny W. Dynamic causal modelling. Neuroimage , 1273-1302 (2003). 13. Feynman RP. The Brownian Movement, The Feynman Lectures on Physics (1964). 14. Breakspear M, Heitmann S, Daffertshofer A. Generative models of cortical oscillations: neurobiological implications of the Kuramoto model.
Front Hum Neurosci , (2010). 15. Deco G, Jirsa VK, Robinson PA, Breakspear M, Friston KJ. The Dynamic Brain: From Spiking Neurons to Neural Masses and Cortical Fields. Plos Computational Biology , (2008). 16. Stephan KE , et al. Charting the landscape of priority problems in psychiatry, part 1: classification and diagnosis.
Lancet Psychiat , 77-83 (2016). 17. Godsil C, Roylse G. Algebraic Graph Theory, Graduate Texts in Mathematics . Springer-Verlag (2001). 18. Li B, Daunizeau J, Stephan KE, Penny W, Hu D, Friston K. Generalised filtering and stochastic DCM for fMRI.
Neuroimage , 442-457 (2011). 19. Carr J. Applications of centre manifold theory . Springer Science & Business Media (2012). 20. Moran RJ, Kiebel SJ, Stephan KE, Reilly RB, Daunizeau J, Friston KJ. A neural mass model of spectral responses in electrophysiology.
Neuroimage , 706-720 (2007). 21. Turing AM. The chemical basis of morphogenesis. 1953. Bull Math Biol , 153-197; discussion 119-152 (1990). 22. Welty JR, Wicks CE, Wilson RE, Rorrer G. Fundamentals of Momentum, Heat, and Mass Transfer . Wiley (2001). 23. Smith WF.
Foundations of Materials Science and Engineering (3rd ed.) . McGraw-Hill (2004). age 22 of 24
24. Vaudano AE , et al.
Causal hierarchy within the thalamo-cortical network in spike and wave discharges.
PLoS One , e6475 (2009). 25. Razi A , et al. Large-scale DCMs for resting-state fMRI.
Netw Neurosci , 222-241 (2017). 26. Seghier ML, Friston KJ. Network discovery with large DCMs. Neuroimage , 181-191 (2013). 27. Friston KJ, Trujillo-Barreto N, Daunizeau J. DEM: a variational treatment of dynamic systems. Neuroimage , 849-885 (2008). 28. Tierney TM , et al. FIACH: A biophysical model for automatic retrospective noise control in fMRI.
Neuroimage , 1009-1020 (2016). 29. Glerean E, Salmi J, Lahnakoski JM, Jaaskelainen IP, Sams M. Functional magnetic resonance imaging phase synchronization as a measure of dynamic functional connectivity.
Brain Connect , 91-101 (2012). 30. Tzourio-Mazoyer N , et al. Automated anatomical labeling of activations in SPM using a macroscopic anatomical parcellation of the MNI MRI single-subject brain.
Neuroimage , 273-289 (2002). 31. Friston KJ, Rotshtein P, Geng JJ, Sterzer P, Henson RN. A critique of functional localisers. Neuroimage , 1077-1087 (2006). 32. Tournier JD , et al. MRtrix3: A fast, flexible and open software framework for medical image processing and visualisation.
Neuroimage , 116137 (2019). 33. Smith RE, Tournier JD, Calamante F, Connelly A. SIFT: Spherical-deconvolution informed filtering of tractograms.
Neuroimage , 298-312 (2013). 34. Friston KJ , et al. Bayesian model reduction and empirical Bayes for group (DCM) studies.
Neuroimage , 413-431 (2016). 35. Rosa MJ, Friston K, Penny W. Post-hoc selection of dynamic causal models.
J Neurosci Methods , 66-78 (2012). 36. Planck M. Über einen Satz der statistischen Dynamik und seine Erweiterung in der Quantentheorie.
Sitzungsberichte der Königlich Preussischen Akademie der Wissenschaften , (1917). 37. Isaacson JS, Scanziani M. How inhibition shapes cortical activity.
Neuron , 231-243 (2011). 38. Xue M, Atallah BV, Scanziani M. Equalizing excitation-inhibition ratios across visual cortical neurons. Nature , 596-600 (2014). 39. Turkheimer FE, Leech R, Expert P, Lord LD, Vernon AC. The brain's code and its canonical computational motifs. From sensory cortex to the default mode network: A multi-scale model of brain function in health and disease.
Neurosci Biobehav Rev , 211-222 (2015). age 23 of 24
40. Borgers C, Kopell N. Synchronization in networks of excitatory and inhibitory neurons with sparse, random connectivity.
Neural Comput , 509-538 (2003). 41. Fritschy JM. Epilepsy, E/I Balance and GABA(A) Receptor Plasticity. Front Mol Neurosci , 5 (2008). 42. Scharfman HE. The neurobiology of epilepsy. Curr Neurol Neurosci Rep , 348-354 (2007). 43. Ziburkus J, Cressman JR, Schiff SJ. Seizures as imbalanced up states: excitatory and inhibitory conductances during seizure-like events. J Neurophysiol , 1296-1306 (2013). 44. Aghakhani Y , et al. fMRI activation during spike and wave discharges in idiopathic generalized epilepsy.
Brain , 1127-1144 (2004). 45. Blumenfeld H. Cellular and network mechanisms of spike-wave seizures.
Epilepsia
46 Suppl 9 , 21-33 (2005). 46. Gotman J, Grova C, Bagshaw A, Kobayashi E, Aghakhani Y, Dubeau F. Generalized epileptic discharges show thalamocortical activation and suspension of the default state of the brain.
Proc Natl Acad Sci U S A , 15236-15240 (2005). 47. Slaght SJ, Leresche N, Deniau JM, Crunelli V, Charpier S. Activity of thalamic reticular neurons during spontaneous genetically determined spike and wave discharges.
J Neurosci , 2323-2334 (2002). 48. Sisterson ND, Wozny TA, Kokkinos V, Constantino A, Richardson RM. Closed-Loop Brain Stimulation for Drug-Resistant Epilepsy: Towards an Evidence-Based Approach to Personalized Medicine. Neurotherapeutics , 119-127 (2019). 49. Lauritzen M, Dreier JP, Fabricius M, Hartings JA, Graf R, Strong AJ. Clinical relevance of cortical spreading depression in neurological disorders: migraine, malignant stroke, subarachnoid and intracranial hemorrhage, and traumatic brain injury. J Cereb Blood Flow Metab , 17-35 (2011). 50. Popovych OV, Tass PA. Control of abnormal synchronization in neurological disorders. Front Neurol , 268 (2014). 51. Tass PA. A model of desynchronizing deep brain stimulation with a demand-controlled coordinated reset of neural subpopulations. Biol Cybern , 81-88 (2003). age 24 of 24 Code availability:
We make all code used to produce the results in this paper available in the following public repository: https://github.com/allavailablepubliccode/Diffusion
Author contributions:
S.P. collected the data under the supervision of M.P.R. and D.C.; C.T. preprocessed the data; all authors designed and performed research, analysed data, and wrote the paper.
Acknowledgements:
We thank all participants and patients for generous involvement in this study, as well as the radiography staff at King’s College London for their excellent technical support in data acquisition. The study was performed at the NIHR ‐ Wellcome Trust King's Clinical Research Facility at King's College Hospital. E.D.F. and R.L were funded by the MRC (Ref: MR/R005370/1); I.R.V. was funded by the Wellcome Trust (Ref: 103045/Z/13/Z) and the BBSRC (Ref: BB/S008314/1); K.J.F. was funded by a Wellcome Principal Research Fellowship (Ref: 088130/Z/09/Z); R.J.M was funded by the Wellcome/EPSRC Centre for Medical Engineering (Ref: WT 203148/Z/16/Z). M.P.R. is supported in part by the National Institute for Health Research (NIHR) Biomedical Research Centre at the South London and Maudsley Hospital NHS Foundation Trust and King's College London, as well as the Engineering and Physical Sciences Research Council (EPSRC) Centre in Predictive Modelling in Healthcare grant number EP/N014391/1, and by the MRC Centre for Neurodevelopmental Disorders grant number MR/N026063/1. This study was also supported by a Medical Research Council (MRC) Programme Grant number MR/K013998/1 and a Wellcome/EPSRC Centre for Medical Engineering Programme Grant number WT/203148/Z/16/Z. The authors would also like to acknowledge support from the Data to Early Diagnosis and Precision Medicine Industrial Strategy Challenge Fund, UK Research and Innovation (UKRI), the National Institute for Health Research (NIHR), the Biomedical Research Centre at South London, the Maudsley NHS Foundation Trust, and King’s College London.