Designing Patient-Specific Optimal Neurostimulation Patterns for Seizure Suppression
Roman A. Sandler, Kunling Geng, Dong Song, Robert E. Hampson, Mark R. Witcher, Sam A. Deadwyler, Theodore W. Berger, Vasilis Z. Marmarelis
DDesigning Patient-Specific Optimal Neurostimulation Patternsfor Seizure Suppression
Roman A. Sandler ∗ , Kunling Geng , Dong Song , Robert E. Hampson , Mark R.Witcher , Sam A. Deadwyler , Theodore W. Berger & Vasilis Z. Marmarelis Department of Physics & Astronomy, University of California, Los Angeles, LosAngeles, CA, USA W. M. Keck Center for Neurophysics, University of California, Los Angeles, LosAngeles, CA, USA Department of Biomedical Engineering, University of Southern California, Los Angeles,CA, USA Department of Physiology & Pharmacology, Wake Forest University, Winston-Salem,NC, USA Department of Neurosurgery, Wake Forest University, Winston-Salem, NC, USAJanuary 15, 2018
Abstract
Neurostimulation is a promising therapy for abating epileptic seizures. However, it is ex-tremely difficult to identify optimal stimulation patterns experimentally. In this study humanrecordings are used to develop a functional 24 neuron network statistical model of hippocam-pal connectivity and dynamics. Spontaneous seizure-like activity is induced in-silico in thisreconstructed neuronal network. The network is then used as a testbed to design and validatea wide range of neurostimulation patterns. Commonly used periodic trains were not able topermanently abate seizures at any frequency. A simulated annealing global optimization al-gorithm was then used to identify an optimal stimulation pattern which successfully abated92% of seizures. Finally, in a fully responsive, or ”closed-loop” neurostimulation paradigm,the optimal stimulation successfully prevented the network from entering the seizure state.We propose that the framework presented here for algorithmically identifying patient-specificneurostimulation patterns can greatly increase the efficacy of neurostimulation devices forseizures.
Epilepsy is a neurological disorder characterized by chronic seizures which affects 1-2% of the USpopulation (Begley et al., 2000). Standard treatments include antiepileptic drugs and resectivesurgery. However, both have major drawbacks. Up to 30% of patients do not respond todrugs; of those who do, many suffer serious side-effects such as nausea, dizziness, drowsiness,and weight-gain (Brodie and Dichter, 1996). Furthermore, Surgery is not an option for manypatients, and when it is, there is a large remission rate within 1-2 years (Engel et al., 2003).In recent years, neurostimulation has emerged as a promising approach to reducing seizures.In 2013, the FDA approved the Neuropace RNS system, the first device for responsive corticalneurostimulation for epilepsy (Sun, Morrell, and Wharen Jr, 2008). However, thus far resultshave shown that neurostimulation provides only palliative relief from seizures rather than a fullcure. For example, the Neuropace device has provided a 60% decrease in median seizure afterthree years (Bergey et al., 2015; Morrell and Halpern, 2016). While these results are impressive,especially considering they were obtained on the most difficult patient category, they are farfrom perfect. For example, 42% of patients did not respond to treatment in the same timeperiod, and no patients were completely seizure free (Bergey et al., 2015). ∗ Corresponding Author: [email protected]. a r X i v : . [ q - b i o . N C ] J a n revious research has shown that the efficacy of neurostimulation could be increased by care-fully designing the temporal pattern of stimulation pulses. Frequency of stimulation (Chkhenkeliand Chkhenkeli, 1997; Cordeiro et al., 2013), periodicity (Wyckhuys et al., 2010; Buffel et al.,2014), and, in the case of multiple electrodes, synchronicity (Good et al., 2009; Nelson et al.,2011; Van Nieuwenhuyse et al., 2014) have all been shown to influence the success of neurostim-ulation. Furthermore, many have argued that stimulation must be custom tailored to the uniqueseizure topology and dynamics of each particular patient (Holt and Netoff, 2014; Mina et al.,2013; Taylor et al., 2015). Presumably, optimal stimulation delivered at the optimal time viaa responsive closed-loop stimulation paradigm could significantly increase therapeutic outcome(Chakravarthy et al., 2009a; Chakravarthy et al., 2007). However, designing optimal neurostim-ulation patterns is extremely challenging in animal models and in human patients. Researcherscannot order seizures ”on-demand” to test a wide range of stimulus patterns. Oftentimes, physi-cians must wait months before learning if a particular stimulus works. Furthermore, researcherscan only stimulate each seizure once, and cannot go ”back in time” to see how the seizurewould have evolved with no stimulation or with different stimulation. Due to these difficultiesthe Neuropace device recommends that physicians keep stimulus frequency fixed at 200 Hz andonly if that fails to adjust the stimulus current. Nonetheless, there is an increasing feeling inthe field that a more principled approach to stimulation design is needed (Nagaraj et al., 2015).In past work we have developed a closed-loop model of the rodent hippocampus and usedthis model to identify the optimal frequency of stimulation needed to reduce network output(Sandler et al., 2015a). However, the above model was limited to 2 reciprocally connectedneurons and thus important features of epilepsy such as population synchrony could not bestudied. Here we use 24 neurons recorded from human hippocampus to reconstruct the pa-tient’s distinct neuronal connectivity and causal dynamics. Spontaneous seizure activity isthen initiated in the reconstructed neuronal network (RNN). Finally, this model is used asan in-silico testbed for designing and testing efficient neurostimulation patterns. The optimalstimulus, obtained via a global optimization algorithm, was found to abate 92% of seizures,and significantly outperformed any other traditional stimulation types. We believe that such apatient-specific algorithmic approach to neurostimulation design can significantly increase theefficacy of neurostimulation devices for epilepsy. One adult patient suffering from pharmacologically refractory temporal lobe epilepsy was surgi-cally implanted with FDA-approved hippocampal electrodes capable of field potential (macro-)and single-unit (micro-) recordings (Ad-Tech Medical Instrumentation Corporation, Racine,WI) for localization of seizures. All procedures were reviewed and approved by the InstitutionalReview Board of Wake Forest Baptist Medical Center, in accordance with National Institutesof Health guidelines. Inclusion in this study was voluntary and consented separately from thesurgical procedure. As the primary contribution of this work is meant to be methodological,data from only a single patient was used to demonstrate proof-of-concept (see discussion).Preoperative planning and intraoperative placement of depth electrodes was performed usinga frameless Brainlab Cranial Navigation System (BrainLab North America, Westchester, IL) toplan and guide electrode entry points, electrode trajectories and target points within the CA3and CA1 subfields of each hippocampus. Electrode localization was confirmed by comparingpreoperative and postoperative MRI in order to identify electrode track and positioning of the”macro” electrode sites with respect to hippocampal morphology. Single unit neural activities(i.e., spike trains) were recorded and isolated using the Blackrock Cervello Elite electrophys-iological recording system with a raw data acquisition frequency of 30 kHz. Electrodes were2xplanted after seizures were localized (14 days). 10 minutes of continuous recordings were usedto estimate all models. Spikes were detected and sorted offline with a 500-5,000 Hz bandpassfilter and discretized using a 2 ms bin.
The aim of the RNN is to use the observed spiking activity to reconstruct in-silico the distincteffective connectivity of the 24 recorded neurons. In other words, for each neuron the RNNattempts to answer which of the other R − t , ˆ y ( t ), was determined by its own past spik-ing activity and the past and present spiking activity of all other N connected CA3 neurons, { x n ( t ) } within a finite memory of M = 100 ms and modeled using a generalized linear model(GLM) with a probit link (see Eq. A2-A4). Each feedforward/feedback filter is either linearor quadratic-nonlinear (2nd order Volterra) (see Eq. A4), as determined by the group reg-ularization algorithm used for model fitting. The feedback component, characterized by thefilter k AR ( τ ), can be intuitively thought of as the afterhyperpotential (AHP) (Spruston andMcBain, 2007) and encapsulates intracellular processes such as the absolute and relative refrac-tory period, slow potassium conductances, and I h conductances. The interneuronal components,characterized by the set of input-output filters { k n } , can be intuited as the waveform of theEPSP from the n th input neuron onto the output neuron. A nonlinear filter is potentially in-cluded to describe interactions between two input pulses, such as paired pulse facilitation anddepression (Song, Marmarelis, and Berger, 2009; Sandler et al., 2015b). More detailed infor-mation on model structure can be found in the appendix. It should be noted that the modelis a ”blackbox”, or entirely based on data, and thus makes no a priori assumptions on thenature of the feedback and interneuronal dynamics. Therefore, the estimated filters include thepreviously listed phenomena as well as more indirect/nonlinear phenomena such as dendriticintegration, spike generation, active membrane conductances, and feedforward interneuronalinhibition (thereby allowing the interactions between two pyramidal cells to be inhibitory).All GLM parameters, which implicitly describe both connectivity and causal dynamics, werefit simultaneously for each neuron using MCP group regularization and a coordinate descentalgorithm (see appendix 4.3.2). It should be noted that because parameters were estimatedfrom spontaneous data rather than from direct perturbations of the network, all parameterestimates may be biased by unobserved inputs (see discussion). A Monte Carlo style shufflingapproach was used to insure all obtained models had significant predictive power and were notsimply overfitting (see appendix 4.3.3). After neuronal connectivity is estimated, network activity can be simulated in a forward mannerusing the estimated coefficients. For each time bin t , the output of the n th neuron, ˜ y n ( t ) ismodeled by an inhomogeneous Bernoulli process (i.e. biased coin flip) whose probability offiring is determined by the past and present spiking activity of other connected neurons:˜ y n ( t ) = ( η n ( t ) , σ ) > u η n ( t ) , σ ) ≤ u (1)Here Φ() is the probit link function (Eq. A2), ˜ η n ( t ) is the linearly weighted combination ofpast and present spiking activity for neuron n (Eq. A3), and u is a standard uniform randomnumber. Thus the neuron spikes when the generated randon number, u , is greater than theneuron’s instantaneous probability of firing, Φ(˜ η n ( t ) , σ ). Equivalently, the neuronal output can3e viewed as being generated by a prethreshold signal, composed of the sum of a deterministiccomponent (Eq. A4) and Gaussian white noise (GWN) of variance σ , followed by a fixedthreshold of 0 which is implicit in probit link formulation (Berger et al., 2012):˜ y n ( t ) = ( η n ( t ) + N (0 , σ ) >
00 ˜ η n ( t ) + N (0 , σ ) ≤ k coefficients (see appendix), the model will generate spikingactivity spontaneously, even without any explicit exogenous input.It was found that when the RNN simulation was allowed to run, it would spontaneously enterdistinct network states for extended time periods. A clustering algorithm was used to identifythe distinct stable dynamical states which emerged in the RNN output (Sasaki, Matsuki, andIkegaya, 2007; Burns et al., 2014). To apply the clustering algorithm, the instantaneous meanfiring rate (MFR) was first computed for each neuron using a 100 ms moving average filter. 5principal components (PCs) were then extracted from the 24 instantaneous MFR vectors (Fig.3c). These 5 PCs were then clustered using a standard k-means algorithm. Finally, for eachstate, the dominant subnetwork of significant neurons was found by identifying the smallest setof neurons which together accounted for >
85% of the cluster energy (see Fig. 3e).
The RNN framework allows binary external stimulation, s ( n, t ) to be applied to neuron (or’electrode’) n at time t by superimposing a spike onto that neuron at that time. Essentially,Eq. 1 is modified to be: ˜ y n ( t ) = ( η n ( t ) , σ ) > u or s ( n, t ) = 10 Φ( η n ( t ) , σ ) ≤ u (3)Intuitively, this model of stimulation allows electrodes to deterministically elicit spikes at precisetimes and thereby ’override’ the neurons’ normal probabilistic mode of spiking defined in Eq.1. After seizure and non-seizure states were identified, a two-stage simulated annealing (SA)algorithm was used to identify the optimal spatiotemporal stimulation pattern, s ∗ ( n, t ), whichwould most quickly and reliably move the network from the seizure state to the non-seizurestate (Kirkpatrick and Vecchi, 1983). Multiple SA streams were run in parallel with 2 differentmodes of stimulation: periodic spiketrains (PTs) and random (Poisson) spiketrains (RTs). Eachelectrode could stimulate at 8 possible parameter values: { OFF, 5 Hz, 20 Hz, 60 Hz, 100Hz, 140 Hz, 180 Hz, 220 Hz } . The electrode parameter determined the periodic stimulationfrequency for PTs and the Poisson rate for RTs. Note that while for PTs, each parameter vectordetermines a unique spatiotemporal pattern, for RTs, each parameter vector only specifies thefiring probability, and could correspond to a near infinite amount of different spatiotemporalpatterns. To further reduce the amount of parameters, the electrodes for neurons which werefunctionally disconnected from the RNN (i.e. had no inputs or outputs) were kept off.The SA algorithm was run on a standard exponential cooling schedule with 120 globaliterations with temperatures logarithmically spaced between T start = 10 and T end = .
01. Eachglobal iteration had 80 local iterations where the temperature was kept constant. At the end ofeach global iteration, the parameter was reset to the minimum of the global iteration (Henderson,Jacobson, and Johnson, 2003).In every SA iteration, the network was initialized in the seizure state. The neurostimulationpattern, s ( n, t ), was designed based on the current parameter values, and applied for 250ms.Since the seizure state was most obviously associated with high firing rates, the SA algorithm4imed to attain the stimulation pattern which most lowered the network MFR in the 100 msafter stimulation ended. Thus the cost function to be minimized was the ratio of stimulus-induced MFR to reference MFR (i.e. the MFR when no stimulation was applied). Note thatthe lowest value the cost function could take was 0, and when it was >
1, the stimulationactually intensified the seizure. At each iteration the same stimulus was applied 3 times, withdifferent random number seeds, and the average MFR ratio over the 3 runs was used as the SAcost function. A sample optimization path can be seen in Fig. 4a.The SA algorithm was run in two consecutive rounds. In both rounds one electrode waschosen whose stimulation parameter value would randomly transition to a neighboring valuewith equal probability; however, the definition of neighboring values was different for bothrounds. In the first round, the neighboring values were defined as the immediately higherfrequency, the immediately lower frequency, and OFF. If the selected electrode was alreadyOFF, all 7 of the frequency values were defined as neighbors (i.e. an OFF electrode wouldtransition to a random frequency with equal probability). Allowing an electrode to turn off ateach iteration encouraged sparsity and facilitated the identification of a subset of electrodes tomore carefully optimize in the second round. In the second round, the OFF electrodes from thefirst round were not adjusted and the remaining electrodes transitioned either to the higher orlower frequency with equal probability.As the SA algorithm does not explicitly penalize nonsparse solutions, a stepwise pruningalgorithm was used to remove superfluous electrodes. During each cycle, the algorithm indi-vidually iterated though all E ”on” electrodes and calculated the cost function if the selectedelectrode was ”turned off”. At the end of each cycle the electrode which caused the greatestdecrease in the cost function was turned off, and the algorithm would continue to the next cycleand consider only the E − This study aims to design a realistic testbed for developing and screening efficient neurostim-ulation patterns for seizure abatement. This testbed, which we have dubbed a reconstructedneuronal network (RNN), is estimated from single unit activity in area CA3 of a human TLEpatient undergoing monitoring for resective surgery (see methods). R = 24 distinct units wereidentified and used for the RNN.The obtained connectivity graph is shown in Fig. 1a. The optimization regularization pathis shown in Fig. 2a,b. Further metrics quantifying the predictive power of the estimated models,including ROC plots and the KS-test are shown in Fig. 2c,d. 18 out of 24 neuronal modelsand 22.8% of all possible connections were found to be significant. The remaining 6 neuronswere functionally isolated from the network: they neither influenced any other neurons or wereinfluenced by them. Of the significant connections, a much larger number than expected werefound to be bidirectionally connected (71.43%, P < . { , , } . Several features of the system can be interpreted from the filters. Note that neurons 8and 14 are connected linearly while neuron 2 is connected with a quadratic filter, thus implyingit exerts some form of short term potentiation. Also note that neuron 14 exerts an entirelyexcitatory influence on neuron 22, while the effect of neuron 8 oscillates between excitation and5 Percent % % Conn% Bidir. 0 5 10 O u t pu t s IO Degree Relationship y=0.69x+1.6 R =0.68. P=0.00026 AB CD
Figure 1: (A) Graph of all identified effective connections. Dashed lines indicate unidirectional connec-tions, while solid lines indicate bidirectional connections. Note that some neurons (like 23) are effectivelydisconnected from the population. (B) Barplot showing % of total possible connections and the % of thosewhich are bidirectional. (C) Positive correlation of the
20 40 60 80
Regularization Strength ~ 1/ o f P a r a m t e r s
20 40 60 80
Regularization Strength ~ 1/ T e s t S e t False Positive Rate (FPR) T r ue P o s i t i v e R a t e ( T P R ) BA C D
Figure 2: (A,B) Plots show ρ , (B) for differentvalues of the regularization parameter, λ . Each plot shows a line for each of the 24 neurons in the RNN.Notice that as in /λ is increased (and thus regularization is weakened), more parameters enter the modeluntil we have a full model. Dots show the optimal λ selected for each neuron. (C) ROC plots for each ofthe 24 neurons, showing model predictive power. Notice that each of the lines are above the dashed blueline (TPR=FPR) which represents a model with no predictive power. (D) vertical KS Plots for each ofthe 24 neurons, along with normalized 95% confidence bounds (Song et al., 2013). As can be seen, mostmodels fall within the bounds. inhibition. Finally, the feedback kernel is composed of initial refractory inhibition, followed byoscillatory bursting activity. Interestingly it oscillates at 5 Hz, which is in the low theta range.This oscillation has been extensively implicated in memory tasks in the hippocampus (Buzsaki,2006; Sandler et al., 2014). 6 .2 Seizure Initiation and Classification Once the effective connectivity and dynamics were estimated, the RNN was allowed to runwithout perturbation and stochastically generate simulated hippocampal CA3 activity, { ˜ y n ( t ) } for all 24 neurons (Pillow et al., 2008). As expected the RNN, which was estimated from normal(nonictal) spiking activity of a TLE patient, generated physiological firing rates with very lowlevels of synchronization (Fig. 3a). In order to induce seizure dynamics, two modificationwere performed. First, the ’in-silico’ neuronal membrane potential was raised by increasing thebaseline firing probability parameter k (see Eq. A3). This mimics the common practice ofinducing seizures experimentally by pharmacologically raising the membrane baseline potential(Fricker, Verheugen, and Miles, 1999; Avoli et al., 2002). Second, the level of stochastic noisedriving the network was reduced by lowering σ . To intuitively understand this modification,note that in the extreme case of σ = 0, the network is entirely deterministic and generates eitherno activity at all oscillates in a fixed limit cycle; alternatively, in the other extreme of σ = ∞ ,the network generates completely random (Poisson) firing. Thus, intuitively, lowering sigmaincreases population control over the neurons and tends to promote the persistent oscillationswhich characterize seizures. It was found that raising the baseline by B = 30% relative to thethreshold and lowering σ to .725 was sufficient to generate spontaneously emerging realisticseizures lasting anywhere between a few seconds to over a minute as seen in real human data(Fig. 3b; Bower et al. (2012) and Truccolo et al. (2014)).In order to gain more intuition about the network dynamics, principal components (PCs)were extracted from the network activity based on instantaneous MFR. Fig. 3c shows thenetwork trajectory through time in PC space. A clustering algorithm was used to identifythe distinct stable dynamical states which emerged in the RNN spiking activity (Fig. 3d-f;Sasaki, Matsuki, and Ikegaya (2007) and Burns et al. (2014)). As can be seen, under theselected { σ, B } parameters, the RNN jumped between only 2 stable states: normal and seizure.The dominant subnetwork of neurons which comprised the seizure cluster is shown in Fig. 3e.These 6 neurons are responsible for most of the spiking activity within the seizure state, andtheir reciprocal connectivity is presumably responsible for maintaining the seizure dynamics.Many neurons maintained their regular firing rate, and 2 neurons even reduced their firing rate.These observations match the heterogeneity of neurons in recorded human seizures (Bower etal., 2012). Finally, the cluster state of the network was identified at each moment in time (Fig.3b,bottom). As can be seen, this method is able to reliably detect when the RNN enters andleaves the seizure state, thus making it a simple seizure detection algorithm (Mormann et al.,2007; Cook et al., 2013; Nagaraj et al., 2015). Once the effective connectivity of the human TLE patient was estimated, and realistic seizureactivity simulated, we aimed to identify a spatiotemporal pattern of electric stimulation whichcould reliably and efficiently induce the network to leave the seizure state. At first glance, thisis a paradoxical task: we want to lower network spiking activity by applying external spikes tothe network. However, several experimental studies have shown that this could be accomplishedusing precisely designed patterned stimulation (Durand and Bikson, 2001; Heck et al., 2014).Our working assumption was that there exist 24 electrodes which could each stimulate a sin-gle neuron without affecting any of the others. A ”pulse” in an electrode at a given time wouldelicit a single contemporaneous spike in the associated neuron at that time. The experimentalimplications of these assumptions will be discussed later. Computationally, however, this in-troduces a highly complex optimization problem. If we only consider whether a given electrodewill be on or off, there are 2 , or over 16 million possibilities. If any of the 24 could take onan arbitrary pattern of spikes over 250 ms, this number would increase to 2 . To make theproblem less formidable, only two modes of stimulation were considered: periodic trains (PTs)7 igure 3: (A) 4 minutes of simulated firing of RNN in physiological conditions. (B) Top: Simulatedfiring of RNN after modifications to induce seizures. Notice the RNN spontaneously enters and leavesthe seizure state. Bottom: the network cluster state through time (see D). (C) Trajectory of RNN activityin (B) within the space of the first 3 principal components. Color indicates time. (D) Results of clusteringalgorithm applied to PC trajectory. Red dots indicate cluster centers. (E) The connectivity between thesubnetwork of neurons which were found to dominate the seizure cluster (yellow). Bigger circles indicatebigger MFRs during seizures. (F) Additional metrics characterizing the two clusters, including proportionof time, MFR, and Fano factor in each state. Note that MFR was scaled to have a maximum of 1 topromote visualization. having a fixed frequency of stimulation, and random, or Poisson, trains (RTs) with differentMFRs (see Fig. 5b). Thus, every electrode was considered a parameter which could take on8 values, including { OF F, Hz, Hz, Hz,
Hz,
Hz,
Hz, Hz } . This allowed atotal of 2 ∗ possible stimulation types.The stimulation length was fixed to 250 ms and a two-stage simulated annealing algorithmwas used to find the optimal parameters which would make the network leave the seizure statemost quickly (see methods). The path of the simulated annealing algorithm is shown in Fig.4a. The optimal stimulation parameters and spatiotemporal pattern are shown in Fig. 4c,d. Itwas found that only 4/24 ( 17%) electrodes needed to be turned on. Of those, 2 electrodes hadfrequencies of 180 Hz, and the remaining two had frequencies of 100 Hz and 220 Hz; furthermore,it was found that periodic stimulation led to better results than Poisson stimulation (Fig. 4b).This confirms experimental evidence showing that high frequency stimulation (HFS) is optimalfor abating seizures (Durand and Bikson, 2001). Interestingly, 2 of the selected electrodes (22and 24) stimulated the epileptic subnetwork (Fig. 3e), while the other two stimulated outsidethe subnetwork, suggesting that direct stimulation of the seizure focus itself may not be themost effective route.By initializing the RNN simulations under identical seizure conditions and using identicalsequences of random numbers, one could compare how a seizure would have evolved underdifferent types of applied stimulation. Fig. 4e,f show how a network initialized in the seizurestate evolves when no stimulation is applied (reference) and when the optimal stimulation inFig. 4d is applied. Additionally, Fig. 4g,h shows the population MFR and PC trajectory inboth simulations. As can be seen from these figures, the optimal stimulation is able induce thenetwork to leave the seizure state in under 250 ms, and more importantly the network stays inthe nonictal state after the 250 ms stimulation ended.Our working premise has assumed that precisely designed independent, or unsynchronized,stimulation across multiple sites could improve responsive neurostimulation. While some workhas supported this hypothesis (Nelson et al., 2011), most DBS studies, due to either experi-mental or theoretical consideration (Durand and Bikson, 2001), have only looked at single site8 igure 4: (A) Evolution of electrode parameters during the second round of simulated annealing. Noticethat some electrodes are consistently kept off as they were not selected during the first round. Whitedots show the relative cost function value for each iteration. (B) Final costs for 6 parallel runs of theSA algorithm using different modes of stimulation: PT,RNB, and RT. Vertical red line shows that thebest results were achieved for PT stimulation and these parameters were the selected ones. (C) Directedconnectivity graph from Fig. 1a, where the optimal electrode frequencies are indicated by the color of thecircles surrounding the neurons. No circles indicates that neuron is not to be stimulated. (D) Rasterplotof optimal spatiotemporal pattern from the SA algorithm. Rasterplots of evolution of a seizure when nostimulation is applied is shown in (E), and when optimal stimulation is applied in (F). Cluster state isshown below both rasterplots (yellow=seizure, blue=normal). (G) Distance from the seizure state clustercenter is shown for both reference and stimulation runs. Notice that stimulation induces the network torapidly move away from the seizure cluster center. This can be seen more clearly in (H) which shows thetrajectory of both runs in PC space (see Fig. 3c. Dotted red lines in (D-G) indicate the beginning andend of stimulation. stimulation where an electric pulse presumably stimulates a large population of neurons simul-taneously (Sun and Morrell, 2014). Furthermore, most studies have used periodic spiketrains(PTs) rather than random spiketrains (RTs) despite a few studies indicating that RTs may besuperior to PTs (Fig. 5a, Wyckhuys et al. (2010) and Van Nieuwenhuyse et al. (2014)). In orderto compare synchronized PT and RT stimulation with independent multi-electrode stimulationwe first attempted to find the optimal stimulation frequency/rate for PT/RT stimulation. Thiswas done by delivering identical patterns of stimulation to all 24 neurons simultaneously for250 ms. The optimal frequency was found by sweeping from 5 Hz to 220 Hz in 50 Monte-Carlo9tyle trials. Once again, in each trial, the RNN was initialized in the seizure state, and identicalrandom numbers were used for each frequency of stimulation in order to allow a fair comparisonof the stimulation frequencies under equivalent conditions (which notably is impossible in reallife). The results are shown in Fig. 5b. Neither synchronized PTs or synchronized RTs, ofany frequency, were found to significantly help in ending seizures. In fact, at many frequenciesthey actually exacerbated the length of seizures. Interestingly however, PTs and RTs at highfrequencies ( >
180 Hz) did temporarily move the network away from the seizure zone (Fig. 6).However, as soon as the stimulation was turned off, the seizure continued. It should be em-phasized however, that these results cannot be generalized beyond the particular patient fromwhom this data was estimated, and high frequency stimulation has been shown to be efficaciousin a large number of patients (Heck et al., 2014). PT Time (ms)
FIT Frequency (Hz)
OFF 50 100 150 200 N on s e i z u r e C l u s t e r D i s t an c e × -4 PTRT
Time (s) % Seizures Aborted
A BC D
Figure 5: (A) Examples of optimal 170Hz PT (top) and 130Hz RT applied equally to all 24 neurons. (B)The optimal PT and RT frequencies were found by comparing seizure length over various frequencies over50 trials. Black lines indicate seizure length in reference (no stimulation) conditions. Line and shadingshow mean ± SEM. (C) Comparison of performance of 6 types of stimulation (see text) over 50 trials.Stimulation was applied for the first 500 ms (vertical red line). Each row shows the network clusterover time (see Fig.2b,d). As can be seen, only the optimal SA stimulation pattern was able to movethe network from the seizure cluster (yellow) to the normal cluster (blue). (D) % of seizures aborted byvarious stimulation modalities.
Finally, the optimal stimulation was compared with other unsynchronized multi-electrodestimulation patterns to insure that the acquired simulated annealing solution is indeed a signifi-cant local, if not global, minima. In each simulation, the optimal stimulation was applied alongwith no stimulation (NON), and various control stimulation patterns including: (1) 200-Hz PT,(2) 200-Hz RT (Fig. 5a,b), (3) random unsynchronized multi-electrode (RM) stimulation havingthe same amount of selected electrodes and total bursts as the SA solution, and (4) stimulationusing the same electrodes as the SA solution, but with mixed frequencies (MF). The results areshown in Fig 5c,d. Again, it can be seen that both types of synchronized stimulation failed toprovide good results. Once again, it can be seen that synchronized stimulation can only tem-porarily move the network out of the seizure state. The mixed-frequency stimulus (MF) wasonly able to abate 5% of seizures confirming that the optimized HFS settings were needed forsuccessful abatement. The random unsynchronized (RM) stimulation, which applied the HFS torandom electrodes, was able to abate 20% of seizures showing that even random unsynchronizedstimulation is, at least in this patient, superior to optimized synchronized stimulation. Finally,the optimal SA solution was able to abort 92% of the seizures. In the four unsuccessful cases, it10 C1 P C xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx Applied StimulationStimulation Output xxxx
AB C
Figure 6: (A) Synchronized 200 Hz periodic stimulus. (B) Network output. Notice that the seizurestemporarily stops during stimulation; however, the seizure returns almost immediately once stimulationis turned off. (C) The same plot in principal-component space. was found that while the stimulation did push the network temporarily out of the seizure state,the seizure restarted after the stimulation was turned off. These results confirm the efficacy ofthe SA algorithm to identify optimal stimulation patterns.
In order to assess the feasibility of model-based responsive neurostimulation, a basic controlstrategy was used (Chakravarthy et al., 2007; Chakravarthy et al., 2009b; Schiff, 2012). Thestrategy consisted of applying the optimal 250 ms neurostimulation pattern (Fig. 4d) in ’real-time’ as soon as a seizure was detected. Causal clustering allowed real-time seizure detection byidentifying when the network shifted from the normal to seizure cluster (Fig. 3). At the end ofthe stimulation, if the network was still in the seizure state, another round of stimulation wasimmediately applied; otherwise, the stimulation therapy ended. Notably, this is the same controlalgorithm currently used in the Neuropace RNS R (cid:13) device (NeuroPace, Inc, 2015). (However,it should be noted that the Neuropace broadens our definition of ’seizure state’ to includeinterictal activity, with the hope that stimulation will prevent the network from ever enteringa full seizure.)Two simulations were performed with and without responsive neurostimulation (Fig. 7).Both simulations were under epileptic { σ, B } parameters and were performed using identi-cal initial conditions and random number generators. As can be seen, in the reference (non-stimulation) case, 4 spontaneous seizures emerges lasting between 5 and 34 seconds. The re-sponsive neurostimulation algorithm was able to detect all these seizures (and the additionalseizures which would have emerged had the network not already been in the seizure state).Furthermore, in most cases a single round of stimulation was able to prevent the network fromgoing into a prolonged seizure state. In one case, at 72 seconds, the first stimulus was un-successful eliminating the seizure, so a second stimulus was applied thereafter and successfullystopped the seizure. In this study single neuron activity from human hippocampus was used to develop a recon-structed neuronal network (RNN) which replicates in-silico the distinctive connectivity andcausal dynamics of the recorded 24 neurons (Fig. 1). the RNN was estimated using a non-parametric/phenomenological approach based entirely on recorded data and which makes few apriori assumptions about the biophysical nature of the network dynamics (Pillow et al., 2008).11 ontrol
Time (s)
Responsive Optimal StimulationAB
Figure 7: (A) 2 minutes of spontaneous network activity under epileptic conditions where 4 seizuresspontaneously emerged. (B) identical network activity in (A), but with responsive neurostimulation,whereby the optimal 250 stimulus was delivered at the times indicated by the red lines. As can be seen,the stimulation was able to avert a prolonged seizure.
The spiking probability of each neuron was estimated using a realistic model incorporating theoutput neuron’s past spiking history and the spiking history of all connected neurons. Groupregularization allowed for the efficient and compact estimation of connectivity and model com-plexity (i.e. whether neuronal interactions are best described by a linear or nonlinear filter).In this study recorded data from a human hippocampus was used to estimate the RNN. Inprinciple the RNN could have been estimated from any simultaneously observed set of pointprocesses, including those produced by artificial spiking neural networks. Human data was cho-sen because it most accurately reflects the conditions which would be met by such an algorithmin practice. Namely, it already inherently contains the desired level of complexity and avoids theneed to make any assumptions on the stochasticity, connectivity, and higher order statistics ofclinical data. Nonetheless, due to (1) the difficulty of obtaining human hippocampal single unitdata and (2) the emphasis of this work on methodology and proof-of-concept only one datasetwas used in the study. A central goal of future work will be to see how well the algorithmperforms on a diverse set of synthetic and recorded datasets and to establish theoretical criteriafor its success.
After the RNN was estimated, seizure dynamics were induced by raising membrane potentialand isolating the network from external noise (Fig. 1), both features which have been impli-cated in initiating physiological seizures (Fricker, Verheugen, and Miles, 1999; Wendling et al.,2003; Warren et al., 2010). Finally, a simulated annealing algorithm was used to design anoptimal stimulation pattern to induce the network to leave the seizure state (Fig. 4). The op-timal stimulation was found to abate 92% of seizures and was successfully used in a responsivestimulation paradigm to prevent seizures from developing in the RNN (Fig. 5-7). These resultslead us to hypothesize that (1) the unique nature of every patient’s seizure focus proscribesany single neurostimulation pattern from being optimal in every patient (2) the distinctive na-ture of the seizure focus can be exploited to algorithmically develop efficient patient-specificneurostimulation patterns.A conceptually similar approach has been used in a previous study where a neural massmodel of the thalamocortical network estimated from patient data was used to explain whyparticular frequencies of stimulation were successful to abate seizures while others were not(Mina et al., 2013). Additionally, such a customized/algorithmic approach has already begun12o be applied to develop neurostimulation patterns for Parkinson’s Disease (PD) (Holt andNetoff, 2014; Grill, Brocker, and Kent, 2014; Brocker et al., 2017; Pe˜na et al., 2017). However,when dealing with PD neurostimulation, one has near instantaneous feedback of the stimulationby assessing its effects on the patient’s tremor (Brocker et al., 2013). The challenges are muchgreater in epilepsy where physicians must oftentimes wait several months before they can accessthe quality of a particular stimulation design due to the infrequency of seizures. Furthermore,physicians cannot ”go back in time” to see how a particular seizure would have evolved hada different stimulation pattern been applied or had no stimulation been applied at all. Thisis particularly important since responsive neurostimulation aims to perturb the network in thepreictal state and thus prevent the seizure from ever occurring; currently, however, devices suchas the Neuropace suffer from a very high false-positive rate. This means that the lack of aseizure following stimulation cannot be used as indicative of its success since in most cases noseizure would have developed regardless of the stimulation pattern. Due to these difficulties,the Neuropace manual recommends that physicians use a 200 Hz periodic stimulus, and if it isunsuccessful to increase the current amplitude (NeuroPace, Inc, 2015). This is despite the factthat the device allows two leads to be independently programmed with a wide range of complexstimuli and a frequency range of 1-333 Hz (Sun, Morrell, and Wharen Jr, 2008). We believe amajor bottleneck in the performance of devices such as the Neuropace, which currently reducesseizures by an impressive but far from perfect 54% (Heck et al., 2014), is not the hardware, butrather the physicians inability to successfully identify optimal stimulation parameters.In this study, model based in-silico neurostimulation optimization is presented as a solu-tion to this vital issue. In this paradigm, a patient-specific model is used as a testbed orhypothesis engine for the design and validation of optimal neurostimulation. This paradigmprovides solutions to many of the experimental issues of neurostimulation design: one mayobtain seizures ”on-demand” by initiating the network in the seizure state. Furthermore, byusing identical sequences of random numbers, one can ”go back in time” and observe how theseizure would have evolved under different applied stimuli. Using this testbed we gained manyinsights into neurostimulation and generated many testable predictions. It was observed thatsynchronized random (Poisson) stimulation provides slight benefits over periodic stimulation -an observation previously observed in the literature both experimentally and in computationalmodels (Wyckhuys et al., 2010; Buffel et al., 2014). However, neither of these stimulus stylesperformed as well as random unsynchronized stimulation over multiple sites, suggesting theneed to conduct more experiments exploiting multiple electrodes for stimulation (Cook et al.,2013; Van Nieuwenhuyse et al., 2014). Furthermore, it was observed that our proposed optimalstimulation disproportionably targeted neurons outside of the epileptic subnetwork (i.e. focus).However, our most important finding was that optimized stimulation over multiple electrodessignificantly outperformed any of the previously mentioned stimulus styles by exploiting theunique connectivity and dynamics of the RNN. Most importantly the optimal neurostimulationpattern identified here can be experimentally validated by applying it to the patient for whichit was estimated on. This ability to experimentally validate our model predictions is lacking inmany of the computational studies exploring neurostimulation.
While the current framework provided very strong computational results, several limitationsneed to be addressed before it can be applied experimentally. The network dynamics wereestimated from spontaneous/observed data and predictive power was used to determine con-nectivity. This Granger-causality approach is biased by unobserved inputs. In this case, neuronswithin CA3 and neurons which project to CA3 (such as those from the entorhinal cortex) arepotential unobserved inputs. Furthermore, our model of electrical stimulation where an elec-trode consistently elicits a single spike uniquely in the neuron that it records, is overly simplistic.Any stimulation will affect at least dozens of surrounding cells (Wei and Grill, 2005; Desai et al.,13014) and have complex temporal dynamics. Additionally, there is evidence that the neuronsan electrode records and those which it stimulates are not equivalent (Histed, Bonin, and Reid,2009). Both of these issues can simultaneously be overcome by experimentally perturbing thenetwork using sequential stimulus pulses across multiple electrodes and analyzing the response.Efficient algorithms are already being developed for how to optimally design such experiments(Lepage, Ching, and Kramer, 2013; Kim et al., 2014)The current model assumes stationarity of neuronal dynamics, while in reality recordedneurons exhibit strong nonstationarities due to synaptic plasticity and electrode drift. Thestationarity assumption was necessary in order to define a tractable optimization problem. Itis unknown to what extent nonstationarities in neuronal dynamics lead to nonstationarities inoptimal stimulation parameters. Future work may aim to address this issue by attempting toidentify optimal neurostimulation patterns in nonstationary networks such as those modeled inRobinson, Berger, and Song (2016).The model presented here relied on single unit activity. In practice, however, recording andsorting stable single units is quite difficult (Gilja et al., 2012; Dhawale et al., 2017). Thus,in future work we shall apply the presented framework for identifying optimal stimulation tocontinuous electrophysiological signals such as ECoG. In this case, each of the specific stepswould be modified, while the overall framework would remain the same. For example, inthe simulated annealing algorithm the current MFR based cost function would need to besubstituted for one which can be applied to continuous signals, such as one based on highfrequency oscillations. Furthermore, due to the difficulty of recording single units during humanseizures, and the difficulty of estimating reliable models from such short data records, thepresent work synthetically induced a seizure. Synthetic seizure dynamics were homogenous,while there is evidence that actual seizures exhibit diversity even within the same patient (Boweret al., 2012; Aarabi and He, 2014). In future work, estimating network dynamics directly fromrecorded seizure data will lead to more conclusive results as less assumptions about seizureinitiation will need to be made. Most importantly, there is a need to validate the obtainedresults experimentally. We imagine that the advocated framework will need to go throughseveral iterations of experimental refinement until the strong computational results achievedhere can be matched in actual animal models or human patients.From a computational perspective, several improvements can be made to the stimulationdesign and control algorithm. While the simulated annealing algorithm considered a very largespace of stimulation possibilities (2 total), several stimulation patterns were not consideredsuch as those which intermixed periodic and random pulse trains and completely arbitrary pulsetrains (Grill Jr and Brocker, 2014; Brocker et al., 2017). Furthermore, results may potentially beimproved by incorporating the relative phases between different pulse trains into the algorithm.Also, a relatively simple control strategy based on the Neuropace device was employed. Theadvantage of this strategy was that the stimulation was independent of the seizure specifics.In the future more sophisticated control algorithms may be employed which emit differentstimulation patterns based on the quality and progression of the specific seizure (Ching, Brown,and Kramer, 2012; Kalitzin et al., 2014; Ehrens, Sritharan, and Sarma, 2015; Tsakalis andIasemidis, 2006; Tsakalis et al., 2006; Jassemidis and Tsakalis, 2012; Schiff, 2012). Overall,while many improvements can be made in the specifics, we believe that the overall frameworkpresented here has the potential to significantly increase the effectiveness of neurostimulationfor epilepsy.
Our speculative and perhaps overly optimistic vision is that in the future epileptic patients willbe implanted with stimulation devices consisting of multiple electrodes which are capable ofboth recording and stimulating (Ryapolova-Webb et al., 2014) and can be independently pro-grammed. Upon implantation, an automatic stimulation algorithm will perturb the network to14stablish safe current levels and to map effective connectivity between the observed areas. Ma-chine learning algorithms will then program the initial stimulation parameters. As is currentlydone in the Neuropace (NeuroPace, Inc, 2015), the device will automatically record all detectedepileptic activity and this data will be uploaded daily to a computer. Then, this data andpatient input will be used offline to analyze the success of yesterday’s stimulation. Finally, areinforcement learning paradigm (Gosavi, 2014), will be used to adjust parameters for the nextday. While admittedly, such a task may seem incredibly difficult to realize we believe that thegrowth of machine learning in the last decade has made this more realistic to accomplish thanever before. Furthermore, since this approach is fundamentally an algorithmic one and does notheavily rely on specific electrode design and placement, it may be backwards compatible withcurrent generation devices and significantly improve their efficacy.
Appendix: Dynamic Connectivity Model and Estimation
A probabilistic model was used to predict the firing probability of a given output CA3 neuronbased on its own spiking history and the past and present spiking activity of all other functionallyconnected CA3 neurons. Thus the probability that a particular neuron, y ( t ) will fire at discretetimebin t is expressed by the probability, ˆ y ( t ):ˆ y ( t ) = P r (cid:16) y ( t ) = 1 | x ( t − τ ) ...x N ( t − τ ) , y ( t − − τ ) (cid:17) = H h x ( t − τ ) ...x N ( t − τ ) , y ( t − − τ ) i (A1)where { x n ( t ) } reflect the N effectively connected spiketrains from CA3, τ reflects the finitememory of the system which ranges from 0 ≤ τ ≤ M , and H [] reflects the mathematical modelwhich is used to describe the dynamical transformation from { x n ( t ) } → y ( t ). The generalizedlinear modeling (GLM) framework was used whereby H [] was decomposed into a linearizedfunction of the inputs, η ( t ), followed by a static nonlinearity, here chosen to be the probit linkfunction (Truccolo et al., 2005; Song et al., 2007):ˆ y ( t ) = Φ (cid:16) η ( t ) , σ (cid:17) = 1 √ πσ Z x −∞ e − ( η ( t ) σ ) (A2) η ( t ), the linearized component of the GLM, takes the form of a nonparametric multiple-inputautoregressive model which describes the dynamical transformation between input and outputspike trains. It consists of a feedforward component, reflecting the effect of the N input cells onthe output cell, and a feedback/autoregressive component reflecting how the cell’s past spikinghistory affects its current probability of spiking. Thus, the output is calculated as: η ( t ) = k |{z} baseline + N X n =1 F [ x n ( t ) , k n ] | {z } interneuronal + F [ y ( t ) , k AR ] | {z } feedback (A3)where F [ x n ( t ) , k n ] models the feedforward effects of input x n ( t ), F [ y ( t ) , k AR ] models feedbackeffects, and k , the constant offset term, models the baseline firing probability.In past studies, feedforward effects were fixed to be either linear (Sandler et al., 2014) ornonlinear (Sandler et al., 2015b; Song et al., 2007), and feedback effects were fixed to be linear(Song et al., 2007). In this study, using the sparse group selection algorithm (see section 4.3.2),it is possible to determine whether particular inputs and feedback effects are best modeled byeither a linear or nonlinear kernel (see Fig. 1d). In this study, linear refers to convolution with alinear filter, k (1) ( τ ), while nonlinear refers to convolution with a quadratic (2nd order Volterra)15lter, k (2) ( τ , τ ) (Marmarelis, 2004; Rajan, Marre, and Tkaˇcik, 2013). Mathematically, theseoperations are respectively defined by Eq.A4a,b: F [ x ( t ) , k (1) ( τ )] = M X τ =0 k (1) ( τ ) x ( t − τ ) (A4a) F [ x ( t ) , k (2) ( τ , τ )] = M X τ =1 M X τ =1 k (2) ( τ , τ ) x ( t − τ ) x ( t − τ ) (A4b)It was found that a memory of 100 ms was sufficient to model the dynamical effects ofmost neurons, and thus M was fixed to 50 (100 ms/2 ms binwidth). In order to reduce theamount of model parameters and thereby increase parameter stability, we applied the Laguerreexpansion technique (LET) to expand the feedforward and feedback filters over L Laguerre basisfunctions (Marmarelis, 2004). Based on previous studies, L = 6 Laguerre basis functions wereused (Sandler et al., 2015a; Song et al., 2007; Song et al., 2009). Correspondingly, the amountof parameters in linear kernels was reduced from M to L (savings of 44 parameters) and in 2ndorder kernels from to M ( M + 1) / L ( L + 1) / α was fixed at 0.542 to reflect this system memory (Marmarelis, 2004). Presumably, only a small portion of the total recorded neurons, R , causally influence anygiven neuron in the reconstructed neuronal network (RNN). This motivates the central task ofidentifying which neurons are effectively connected and which are not (Bullmore and Sporns,2009; Fallani et al., 2014). Most methods which aim to estimate effective connectivity adopt aGranger causality approach whereby neuron A is effectively connected to neuron B only if it canhelp predict when neuron B will spike (Krumin and Shoham, 2010; Kim et al., 2011; Sandleret al., 2014; Zhou et al., 2014). Here a penalized group regression approach was adopted whichimplicitly maximizes predictive power while improving computational efficiency over commonlyused stepwise input selection methods (Song et al., 2013; Song et al., 2016).To proceed with penalized group regression, Eq. A2 is first recast in matrix form: ˆ y = Φ( η ) = Φ( V c ) (A5)where ˆ y and η are the vectors consisting of ˆ y ( t ) and η ( t ) for 1 ≤ t ≤ T , V is the design matrixconsisting of the convolved inputs and, for quadratic kernels, their cross products (Marmarelis,2004), and c is the vector of model parameters to be estimated. In the 24 neuron RNN, thereare 300,000 observations and 649 unknown parameters.The input parameters are divided into 2 R groups consisting of 2 groups for every putativeinput: one group for the L L ( L + 1) / R = 24 recorded neurons, there were 48 groups: 24 1storder kernel groups and 24 2nd order kernel groups. The objective is now to find the optimalparameter vector, c ∗ which minimizes the cost function composed of the sum of the negativelog-loss likelihood and the group regularization term, P ( c ): C ( c ; y , V ) = T X t =1 (cid:16) y ( t ) log ˆ y ( t ) + (1 − y ( t )) log (1 − ˆ y ( t )) (cid:17) + R X g =1 P ( c g ) (A6)where c g is the group parameter vector containing only the parameters within group g . Thefunction of the regularization term is to automatically set to 0 any parameter groups whichare not found to significantly influence the output - thus implicitly estimating sparse functional16onnectivity. Here, the group minimax concave penalty (MCP) regularizer was chosen overthe more conventional group LASSO because (1) it induces much less regularization basedparameter shrinkages (biases) than the latter (2) it leads to much sparser solutions than thelatter (Breheny and Huang, 2014; Zhang, 2010). The MCP regularizer is defined as: P ( c g ; λ, γ ) = ( λ k c g k − k c g k γ k c g k ≤ γλ γλ k c g k > γλ (A7)where λ determines the strength of regularization and γ , which was fixed at 3, determinesthe range over which MCP regularization is applied (Breheny and Huang, 2014). The groupcoordinate descent algorithm outlined in Breheny and Huang (2014) was used to find c ∗ .The optimal λ value was selected using a warm-start regularization path approach using 90logarithmically spaced λ values. At each iteration, the Pearson correlation, ρ was computedbetween the recorded spiketrain, y ( t ), and the estimated spike probabilities, ˆ y ( t ), on a testingset consisting of 20% of data randomly selected from V . ρ was used since (1) it was previouslyfound to be a robust metric of similarity between spiketrains and continuous signals (Sandleret al., 2014) and (2) it led to sparser solutions than the more commonly used cross-entropyerror. Finally, the optimal λ was selected as the largest λ which achieved >
99% of the max ρ value. The regularization path can be seen in Fig. 2.Since any regularization will necessarily bias obtained parameters to varying degrees, allparameters of nonsparse groups were reestimated without sparse groups and without the regu-larization term. Note that due to computational efficiency and simplicity the initial search for λ uses a logit link function (Breheny and Huang, 2014), while the final reestimation was doneusing the probit link of Eq. A2. All computations were done in Matlab using custom codeavailable upon request. A standard 3.2Ghz, 6-core desktop computer was able to estimate the24-neuron RNN in approximately 3 hours. To avoid overfitting, Monte Carlo style simulations were used to select those models whichrepresent significant causal connections between input and output neurons and do not just fitnoise (Sandler et al., 2014). The following procedure was used: in each run the real input wasdivided into 40 blocks and these blocks were randomly permuted with respect to the output. Amodel was then generated between the permuted inputs and the real output, and the Pearsoncorrelation coefficient, ρ i , was obtained as a metric of performance. T = 40 such simulationswere conducted for each output and a set of performance metrics, { ρ i } Ti , was obtained. Then,using Fisher’s transformation, we tested the hypothesis, H , that ρ was within the populationof { ρ i } . If this hypothesis could be rejected at the 99.99% significance level, the model wasdeemed significant. The very conservative threshold ( P < . y ( t ) with the true binary output y ( t ). Receiver Operatorcharacteristic (ROC) curves plot the true positive rate against the false positive rate over theputative range of threshold values for the continuous output, y ( t ) (Fig. 2c, Zanos et al. (2008)).The second metric used was the discrete KS test (Haslinger, Pipa, and Brown, 2010; Song et al.,2013) which compares the ISI distribution of the time rescaled probabilistic estimates with thatof a homogenous Poisson process (Fig. 2d). All model assessment metrics were evaluated onthe testing set. 17 cknowledgements This work was supported by NIH grant P41-EB001978 to the Biomedical Simulations Resourceat the University of Southern California and DARPA contracts N6601-14-C-4016 and N66601-09-C-2081. The authors thank Dr. Daniel E. Couture, and Dr. Gautam Popli, Wake ForestBaptist Medical Center for their contribution of human hippocampal neural recordings.
References [1] Aarabi, Ardalan and Bin He (2014). “Seizure prediction in hippocampal and neocorticalepilepsy using a model-based approach”. In:
Clinical Neurophysiology
Progress in neurobiology
Epilepsia
Neural Systems and Re-habilitation Engineering, IEEE Transactions on
Neurology
Epilepsia
Statistics andComputing , pp. 1–15.[8] Brocker, David T et al. (2013). “Improved efficacy of temporally non-regular deep brainstimulation in Parkinson’s disease”. In:
Experimental neurology
Science Translational Medicine
New EnglandJournal of Medicine
International journal of neural systems
Nature Reviews Neuroscience
Proceedings of the National Academy of Sciences
Rhythms of the Brain . Oxford University Press.[15] Chakravarthy, Niranjan et al. (2007). “Controlling synchronization in a neuron-level pop-ulation model”. In:
International journal of neural systems
Journal of Combinatorial Optimization
Annals of biomedical engineering
PhysicalReview E
Stereotactic and functional neurosurgery
The Lancet Neurology
Epilepsy research
Frontiers in neuroengineering eLife
Proceedings of the IEEE
Frontiersin neuroscience
Neurology
Philosophical Transactions of the RoyalSociety B: Biological Sciences
The Journal of Phys-iology
Nature neuroscience
International journal of neural systems
Simulation-based optimization: parametric optimization techniquesand reinforcement learning . Vol. 55. Springer.[32] Grill Jr, Warren M and David T Brocker (2014).
Non-regular electrical stimulation patternsdesigned with a cost function for treating neurological disorders . US Patent 8,923,981.[33] Grill, Warren M, David T Brocker, and Alexander R Kent (2014).
Devices, systems andmethods for deep brain stimulation parameters . US Patent 20,140,350,634.[34] Haslinger, Robert, Gordon Pipa, and Emery Brown (2010). “Discrete time rescaling the-orem: determining goodness of fit for discrete time statistical models of neural spiking”.In:
Neural computation
Epilepsia
Handbook of metaheuristics . Springer, pp. 287–319.[37] Histed, Mark H, Vincent Bonin, and R Clay Reid (2009). “Direct activation of sparse,distributed populations of cortical neurons by electrical microstimulation”. In:
Neuron
Journal of computational neuroscience
Pacemaker for treating phys-iological system dysfunction . US Patent 8,197,395.[40] Kalitzin, Stiliyan et al. (2014). “Multiple oscillatory states in models of collective neuronaldynamics”. In:
International journal of neural systems
PLoS computational biology
Neural computation .[43] Kirkpatrick, Scott, MP Vecchi, et al. (1983). “Optimization by simmulated annealing”.In: science
Computational intelligence andneuroscience
Journal of computational neuroscience
Nonlinear dynamic modeling of physiological systems . Wiley-Interscience.[47] Mina, Faten et al. (2013). “Modulation of epileptic activity by deep brain stimulation: amodel-based study of frequency-dependent effects”. In:
Frontiers in computational neuro-science
Brain
Neurosurgery clinics of North America
Journal of Clinical Neurophysiology
Epilepsy research
RNS R (cid:13) System User Manual . Accessed: 2015-06-30.[53] Pe˜na, Edgar et al. (2017). “Particle swarm optimization for programming deep brainstimulation arrays”. In:
Journal of Neural Engineering
Nature
Neural computation
Neural Computation .[57] Ryapolova-Webb, Elena et al. (2014). “Chronic cortical and electromyographic record-ings from a fully implantable device: preclinical experience in a nonhuman primate”. In:
Journal of neural engineering
Journal of computationalneuroscience , pp. 1–15.[59] — (2015a). “Hippocampal closed-loop modeling and implications for seizure stimula-tion design”. In:
Journal of neural engineering
Journal of neuroscience methods
The Journal of Neuroscience
Neural control engineering: the emerging intersection betweencontrol theory and neuroscience . MIT Press.[63] Song, Dong, Vasilis Z Marmarelis, and Theodore W Berger (2009). “Parametric and non-parametric modeling of short-term synaptic plasticity. Part I: Computational study”. In:
Journal of computational neuroscience
Biomedical Engineering, IEEE Transactions on
Neural Networks
Journal of computational neuro-science
IEEE Transactions on Neural Systems andRehabilitation Engineering .[68] Spruston, N and C McBain (2007). “Structural and functional properties of hippocampalneurons”. In:
The hippocampus book , pp. 133–201.[69] Sun, Felice T and Martha J Morrell (2014). “Closed-loop Neurostimulation: The ClinicalExperience”. In:
Neurotherapeutics , pp. 1–11.[70] Sun, Felice T, Martha J Morrell, and Robert E Wharen Jr (2008). “Responsive corticalstimulation for the treatment of epilepsy”. In:
Neurotherapeutics
Frontiers in Neuroscience
9, p. 202.[72] Truccolo, Wilson et al. (2005). “A point process framework for relating neural spikingactivity to spiking history, neural ensemble, and extrinsic covariate effects”. In:
Journalof neurophysiology
The Journal of Neuroscience
International Journal of Bifurcation and Chaos
Cybernetics and Systems Analysis
Brain stimulation .[77] Warren, Christopher P et al. (2010). “Synchrony in normal and focal epileptic brain: theseizure onset zone is functionally disconnected”. In:
Journal of neurophysiology
Journalof neural engineering
Brain
Epilepsia
Neural Systems and Rehabilitation Engineering, IEEE Transactionson
The Annals of Statistics , pp. 894–942.[83] Zhou, Douglas et al. (2014). “Granger causality network reconstruction of conductance-based integrate-and-fire neuronal systems”. In: