Locally embedded presages of global network bursts
Satohiro Tajimaa,b,c,1, Takeshi Mitad, Douglas J. Bakkume, Hirokazu Takahashid,f, and Taro Toyoizumica
Department of Basic Neuroscience, University of Geneva, Centre Médical Universitaire, Genève 1211, Switzerland; bPrecursory Research for Embryonic
Science and Technology (PRESTO), Japan Science and Technology Agency, Kawaguchi, Saitama 332-0012, Japan; cRIKEN Brain Science Institute, Wako,
Saitama 351-0198, Japan; dGraduate School of Information Science and Technology, The University of Tokyo, Bunkyo-ku, Tokyo 113-8656, Japan;eDepartment of Biosystems Science and Engineering, ETH Zurich, Basel 4058, Switzerland; and fResearch Center for Advanced Science and Technology, The
University of Tokyo, Meguro-ku, Tokyo 153-8904, Japanburst| network | prediction | nonlinear dynamics | spontaneous activityT

he collective firing dynamics of neural population have been
related to emergent functions and dysfunctions in the brain.
Unraveling the structure of dynamics emerging from the collective behavior of neural networks is an ultimate issue in studies
of complex systems (1). The simplest and most profound example of collective dynamics is synchronous burst: the simultaneous
occurrence of closely spaced action potentials (“bursts”) across alarge number of neurons in the network. Within neural circuits
in vivo, the synchronous bursts are believed to serve importantroles in information storage and transmission (2, 3) as well as indisease states including epileptic seizures (4–6). Intriguingly,
neurons cultured in vitro also display spontaneous synchronousbursts akin to those observed in vivo (7–9). This indicates thatthe collective bursting is a universal behavior of neural networks
in a wide variety of preparations, for which new approaches tounderstanding its foundations are much needed (1, 10, 11).
Classic biophysical models of the synchronized bursts assume
broadly diverging synaptic interaction that propagates activity ofa single cell to the remaining neural population (4, 12–14). Thisview has been partially supported by experiments demonstratingthat stimulating a single neuron can change the large-scale dynamics of neural populations (15–17). These models and experiments suggest a close link between single neuron activity andthe global state of a neural population. On the other hand, theprevious stochastic models do not specify when a global burst
occurs; forecasting a spontaneous occurrence of synchronousburst from the preceding individual neural activities in nonbursting states remains a highly challenging issue, particularly inreal nervous systems which feature the heterogeneity of local–
global coupling (e.g., nonuniform coupling strength between
individual neural activity and the overall firing of the population)(18). Indeed, recent studies demonstrate potentially diverse dynamical mechanisms accompanied by complex but reproducible
patterns of activity sequences leading to the synchronous bursts

www.pnas.org/cgi/doi/10.1073/pnas.1705981114(11, 19), implying the existence of nonlinear deterministic
mechanisms around the burst initiations.
The present paper investigates the deterministic aspects of burst
initiation. In particular, we question whether individual neuron
dynamic in nonbursting periods can predict an upcoming synchronous burst in cultured neural population, which provides a simplest
model system of self-organized neurons. To cope with the heterogeneous nature of local–global coupling, we introduce an modelfree forecasting method based on a nonlinear state-space reconstruction technique. The method extends the data-driven nonlinear
forecasting techniques (20–25) to characterize the global–local
correspondence in network dynamics, and reveals the deterministic
relationships among individual neuron traces and the global state.
Using this method, we report that dynamic of only one neuron canrobustly predict the upcoming synchronous burst in the neural
population, at high signal-to-noise ratio. Surprisingly, the burst
predictability with an appropriately chosen neuron even outperforms that with the global (i.e., population mean-field) activity ofneural population, demonstrating that macroscopic fluctuationof neural population is better predicted by the microscopic dynamicof a specific single cell, rather than the macroscopic state itself. Thisphenomenon is explained based on the heterogeneous causal interactions among neurons. The present findings demonstrate theeffectiveness of the nonlinear state-space reconstruction techniques
for analyzing the causal coupling in in vitro system dynamics.
Results
Spontaneous Synchronous Bursts in Cultured Neurons. Rat corticalneurons were grown on high-dense complementary metal–oxide–
semiconductor (CMOS)-based multielectrode arrays (MEAs) for
SignificanceThis paper reports an approach to detecting the “early warnings” of upcoming global state transitions in a network basedon its local dynamics, demonstrating that seemingly stochastic
global events can be predicted by local deterministic dynamics.
Based on the method using a nonlinear state-space reconstruction, we show that, surprisingly, dynamics of individual
neurons can robustly predict the upcoming synchronous burst in
the neural population at high signal-to-noise ratios, which evenoutperform the predictions based on population activity. We
explain this apparently counterintuitive property by the network
structures realizing in the nonbursting period, which is supported by a manipulative experiment and analyses. These results
reveal basic properties of the bursting network dynamics.
Author contributions: S.T., T.M., H.T., and T.T. designed research; T.M., D.J.B., and H.T.
contributed new reagents/analytic tools; S.T. and T.T. analyzed data; and S.T., T.M., H.T.,
and T.T. wrote the paper.
The authors declare no conflict of interest.
This article is a PNAS Direct Submission.
Freely available online through the PNAS open access option.
1

To whom correspondence should be addressed. Email: satohiro.tajima@gmail.com.

This article contains supporting information online at www.pnas.org/lookup/suppl/doi:10.
1073/pnas.1705981114/-/DCSupplemental.

PNAS Early Edition | 1 of 6

NEUROSCIENCE

Spontaneous, synchronous bursting of neural population is a widely
observed phenomenon in nervous networks, which is considered
important for functions and dysfunctions of the brain. However, how
the global synchrony across a large number of neurons emerges froman initially nonbursting network state is not fully understood. In thisstudy, we develop a state-space reconstruction method combinedwith high-resolution recordings of cultured neurons. This method extracts deterministic signatures of upcoming global bursts in “local”
dynamics of individual neurons during nonbursting periods. We findthat local information within a single-cell time series can comparewith or even outperform the global mean-field activity for predicting
future global bursts. Moreover, the intercell variability in the burst
predictability is found to reflect the network structure realized in the
nonbursting periods. These findings suggest that deterministic local
dynamics can predict seemingly stochastic global events in selforganized networks, implying the potential applications of the present methodology to detecting locally concentrated early warnings of
spontaneous seizure occurrence in the brain.

APPLIED
MATHEMATICS

Edited by Terrence J. Sejnowski, Salk Institute for Biological Studies, La Jolla, CA, and approved July 28, 2017 (received for review April 14, 2017)

A

C
RecordCell ID

20

Stimulation

Neuron

40
60
80

500μm

Excitatory cells
Inhibitory cells

1
0.8
0.6
0.4

0

-1

Chip #

-2
-3
-4
0

1
2
3
Log ISI(normalized)“Burst”

B

Global mean (normalized)

D“Burst” 14
4
24 0
2427
4
14 0
1444
37

ElectrodeLog fraction of time bins

20 s

50μm

0.2
0
20 s

Fig. 1. Global bursts in cultured neural populations. (A) CMOS-based recording apparatus. (B) Spatial distribution of cells recorded in a representative preparation, Chip 1440. (C) Spontaneous activity of neurons in thesame preparation, Chip 1440. The color indicates the spike count within each5-ms time bin (white: 0 spikes; black: >5 spikes). The abscissa represents thetime. The majority of neurons fired together during network bursts (columnsof dark-colored bins). (D) The temporal evolution of global mean field(normalized population firing rate) computed from C. The horizontal dashedline indicates the threshold for burst detection for this preparation. (Inset)
The distributions of interspike intervals (ISIs) in five different preparations(solid lines) and the burst-detection thresholds for the individual preparations (dashed lines). The horizontal axis shows the normalized log-ISI (thelogarithm of inverse global activity; Methods).

12–41 d in vitro (Fig. 1A) (26, 27). The devices allow us to simultaneously stimulate and record from multiple neurons lying on thearray surface at a high spatiotemporal resolution (Methods). Each
recording yielded simultaneous traces of 60–98 cells at a temporal
resolution of 2 kHz. Fig. 1B shows the cell distribution recordedwith a representative array preparation. In our preparations ofsparse culturing, the cell bodies were well isolated, and the action
potential of each cell was identified accurately based on its spatial
location and waveform. Spike action potentials were detected using
standard LimAda method (28).
We first analyzed the spontaneous activities of the recorded
neural population. During the spontaneous firing, the neurons
showed intermittent synchronous bursts, which were interleavedwith silent periods (Fig. 1C). To describe the macroscopic network
dynamics, we defined population activity by the mean firing rate
across all recorded neurons (Fig. 1D). The synchronous bursting
period was defined by time bins in which the population activity
exceeded a threshold determined for each preparation, according
to the previously established method based on interspike interval
distributions (27) (Fig. 1D, Inset). The timing of synchronous burst
occurrence did not show any clear periodicity. Duration of bursts
varied from about 200 to 800 ms. The different preparations
showed distinct complex patterns of collective activity dynamics
when visualized with standard dimensionality reduction techniques such as principal component analysis (Fig. S1).
Predicting Global Network Bursts. To capture the potentially complex temporal structures in neural activities around bursting dynamics, we developed a model-free method that is able to quantify
the general relationships between the network and single-element
dynamics in high-dimensional nonlinear systems. The present protocol is based on the state-space reconstruction using the delayembedding theorems in deterministic dynamical systems (29, 30).
2 of 6 | www.pnas.org/cgi/doi/10.1073/pnas.1705981114

The method can be interpreted as an application of the previous
nonlinear forecasting methods, known as “simplex projection” (20)and “convergence cross-mapping (CCM)” (31), but here appliedin a cross-scale domain with time delays and targeting at a specific global event. It extends/integrates the previous methods inthree aspects. First, the method focuses particularly on the relationships between a global state variable (e.g., the mean-field
variable) and local variables. Second, the method quantifies thepredictability of the global variable based on temporally distantdynamics of each local variable. Third, it constructs a detector ofa specific global event (e.g., the spontaneous synchronous burst)
through bidirectional mappings between the global and local
variables. Because “targeting” a specific global event is a prominent
feature of our present approach, we will sometimes mention themethod parsimoniously as “event-targeted CCM (eCCM)” in what
follows. Note that the idea of mapping variables under time delays
is also found in previous studies (23). In addition, the general
problem of how the nonlinear signals in time series are affected byscale has also been discussed in other contexts (20, 24, 25). Below
we summarize our basic protocol for the eCCM analysis.
Fig. 2 illustrates the analytic protocol (further details areprovided in Methods). We first reconstruct the dynamics of
population activity and each single-neuron activity, respectively,
in delay-coordinate state spaces (Fig. 2 A and B):
bdt = bt , bt-τ . . . , bt-ðd-1Þτ ,

[1]



xdi,t = xi,t , xi,t-τ . . . , xi,t-ðd-1Þτ ,[2]where bt, and xi,t are the population mean-field and neuron i’s
activity at time t, respectively; d denotes the embedding dimensions (the number of delay coordinates); τ is the unit delay size.
In this study, the local time series xi,t was defined by each neuron
i’s spike train smoothed with a Gaussian kernel, and bt was defined by the mean of xi,t over all of the simultaneously recorded
neurons (Methods). For convenience, we refer to btd and xi,td as
the “global” and “local” state trajectories, respectively. We measure the accuracy of synchronous-burst prediction by seeking“presages” of global bursts in each local state trajectories, basedon the accuracy of forecasting using a nearest-neighbor model(Fig. 2A; see also the figure legend for a step-by-step descriptionof the protocol). The rationale behind this protocol is as follows:
if neuron i has enough information to predict the future synchronous burst occurring after a time span −Δt (−Δt > 0), neuron i’s
state up to time −Δt before the targeted burst should be alreadydifferentiated from the states without any future bursting, forming a cluster of states that deviate from those predicting nonbursting periods in the reconstructed state space. This early
warning is detected by finding the bidirectional maps between
the targeted global event (i.e., the synchronous burst) and the
local states (i.e., single-neuron dynamics) as illustrated in Fig.
2A. Similarly, we define the accuracy of “postdiction” of bursts(i.e., detecting a burst event based on neural activities after its
occurrence) with the same protocol as that for prediction except
for using time span Δt with a positive value. We iterated thisprocedure for all of the neurons to characterize the burst predictabilities in individual neurons. In addition, we quantified the “selfpredictability” of the mean-field activity by replicating the above
procedure based on its own (global) state trajectory bt, instead of a
local state trajectory in individual neurons xi,t (Fig. 2B).
Using this method, we investigated the predictability of synchronous bursts based on the individual neural activities or the
population mean field. In particular, we addressed two keyquestions: (i) how accurately can each single neuron predict the
future synchronous burst events, and (ii) if a subset of neurons
predicts synchronous bursts better than others, how are they
related to the underlying neural network structure?
Tajima et al.

B Forecasting by the global signal

A Forecasting by a local signal(local-global prediction)(global-global prediction)

Local state
Global state(the mean field) (a single cell)x1(t)Pick a global burst state.

(ii) Trace-back map

b(t)

(ii)

(ii)

Trace-back the time along
the local-state trajectory.

(i)

Trace-back the time along
the global-state trajectory.

(iii) Finding neighbors(iii)(iv) Local-to-global map

Tajima et al.

0
-1000

0

0
1000

Robustness
1

0.5
0
0

Time from burst peak (-Δt) [ms]

0
-1000

0

0.5
0
1000

D1

1

0.5

0.5

0
-1000

0

Prediction Postdiction
Time from burst peak (-Δt) [ms]

0
1000

1

F Superiority to firing rate
Fraction of cells

0.5

1

0.5

0.5
0
-1 0

G
Relative success rate
(single cell - mean field)

C1

0.5

Success rate(even trials)

1437
1440
1444
2427
2440

1

0
-1 -0.50
0.5
1
Relative success rate
(eCCM - FR)
Single cell vs. mean field

0.5

eCCM
FR

0
-0.5
-1
0

APPLIED
MATHEMATICS

0.5

Fo
r
Fo ecasrec t bas y lot b ca
yg l
lob
al

0.5

Success rate
(odd trials)

1

100 200 300 400
Success rate rank

Fig. 3. Presages of global bursts within single neurons. (A) Burst detection
with the eCCM based on a single representative neuron. The black trace withthe gray shade represents the median and quartiles of the estimated global^ t , in individual trials around burst peaks. The burst-detection
states, bthreshold (dashed green line) was defined by the 95th percentile of theestimated global activity outside the bursting period. The red curve shows
the fraction (“hit rate”) of successfully detected bursts in all of the burst
trials. The curves are smoothed with 250-ms boxcar kernel for the visualization. Δt takes a negative value when predicting and a positive value when“postdicting” burst events. (B) Burst detection with the eCCM based on themean-field activity in the population. (C) Burst detection with the momentary firing rate (FR) of the same representative single neuron as the one
shown in A. (D) Burst detection with the momentary mean-field firing rate.
(E) Burst predictability is a robust property of individual neurons. Dividing
the data into two halves (the odd and even trials of spontaneous burst)
shows that variability in burst-detection success rate is highly consistent overtime. Each marker corresponds to a neuron. The filled markers represent theresults using the mean-field activity. The success rate was averaged over therange of time span −200 ms < Δt < −50 ms. (F) Burst detection with the eCCM ismore accurate than that based on the momentary firing rate. The comparison of
success rates in burst detection with the eCCM and with the firing-rate-based
method. (Inset) The control analysis in which we compared the success rates ofthe eCCM with that of the firing-rate-based method using multiple time bins, by
matching the number of bins used in the eCCM. The gray bars show the fractionof cells in all of the preparations. The arrowheads indicate the medians. The
colored lines show the fractions in the individual preparations. (G) A single cell
can outperform the mean field at predicting the spontaneous network bursts,
particularly when we use the information in the temporal patterns of neuralactivity rather than the momentary activity. The light- and dark-blue curves,
respectively, show the success rates in burst detection with single neurons relative to that with the mean field, for each of the eCCM and the firing-ratebased analysis. For each method, the cells were sorted based on the success rate.

PNAS Early Edition | 3 of 6

NEUROSCIENCE

0

0
1000

Predicted mean-field(with eCCM)

0
-1000

1

Detection success rate

Single-Neuron Dynamics Can Predict Global Network Bursts. First, we
asked how accurately individual neurons predict (or postdict) the
occurrence of synchronous bursts, and how their predictions
compare with that based on the population mean field, which
defines synchronous burst. We quantified the rate of successful
synchronous-burst predictions in each single neuron with keeping the false-positive rate at 5% (Methods). Interestingly, in somecells (e.g., Fig. 3A), the rates of successful burst predictions
based on the eCCM analysis could be relatively high (>50% of
trials) even when the network appeared to be silent in terms of

0.5

Detection success rate

Fig. 2. eCCM analysis relates the global-network events to the preceding
local states through the state-space reconstruction. (A) Forecasting by a local
signal. In each box, the left and right trajectories depict the delay-coordinate
reconstruction (see the text) of the mean field (global state) and of a singleneuron activity (local state), respectively. A representative neuron (Cell 1) isshown. Dark-blue dots: the time points of peak global activities during the
bursts; magenta dots: the temporally preceding states from the individualglobal-activity peaks (here, shift time = −100 ms); magenta dots: theneighboring states within the reconstructed state space for Cell 1 (10 nearest
neighbors for each burst peak); cyan cross: the output of the present protocol (i.e., the predicted mean field). A single neuron’s ability at predicting
each burst event is quantified in five steps: (i) selecting a time point t on a
local-state trajectory (say, of neuron i) that corresponds to a peak population activity found in the global-state trajectory (which is a target
bursting state to be predicted), (ii) tracing back the local-state trajectory fora given time span, Δt, (iii) collecting “neighbor” time points corresponding
to states nearby t − Δt in the local-state trajectory (avoiding points temporally too near), (iv) mapping the neighbor time points back into theglobal-state trajectory, and (v) forwarding the time with Δt on the globalstate trajectory to obtain the global-state prediction (the cyan cross). Whenthe global-state prediction deviates from the nonbursting-state distribution
by having a significantly large value, we defined it as a successful detectionof burst (refer to Methods for formal descriptions). The accuracy of burst
prediction is higher when the traced-back states (the red dots) deviate morefrom the ones in the nonbursting period in the local trajectory. Although the
combination of tracing time in the retrograde direction with Δt (step ii) andmutual mapping between global and local states (steps i and iv) are newly
introduced features in this study, step v alone could be regarded as a variantof the previously proposed forecasting methods based on simplex projection(20). On the other hand, if steps iii–v alone are applied (i.e., without any
constraint of the targeted event, unlike our current protocol with steps i andii), it reduces to the CCM algorithm including time delays (23, 31) but applied
to the cross-scale predictions. (B) Forecasting by a global signal. The same as
in A, except that the global state itself is used instead of the local state.

0.5

Mean field
B(global-global
E
prediction)

FR (normalized)

(v)

1

Chip#

Trace-forward the timealong the global-state
trajectory to obtainthe forecast.

1

Detection success rate

Single cellA(local-global
prediction)(iv)

Detection success rate(v) Trace-forward map (v)Trace-forward the timealong the global-state
trajectory to obtainthe forecast.

(iii & iv)

De
thr tecti
es on
ho
ld

Identify the corresponding
global state for each of
those points.

Pick the data points near
the traced-back statein the global-state space.

Predicted mean-field(with eCCM)

Pick the data points near
the traced-back statein the local-state space.

FR (normalized)b(t)x1(t+τ)b(t+τ)(i)b(t+τ)(i) Global-to-local map
Pick a global burst state,
and identify thecorresponding
local state.

Global state(the mean field)

the global mean-field firing rate (Fig. 3D, −1,000 ∼ −500 ms).
Corresponding to this, the success rate of burst detection using asingle-neuron dynamic in the eCCM could be even higher thanthat with the global mean-field dynamics (Fig. 3B). The predictability of burst occurrence varied among neurons, suggesting
the heterogeneity of burst predictability across neurons. However, the accuracy of predictions in individual neurons was highly
consistent among different sets of burst trials within each preparation (ρ = 0.79 ∼ 0.93, P < 2 × 10−23, Spearman’s rank correlation between the first and the second halves of trials; Fig.
3E), demonstrating that this cell-to-cell variability reflects a robust property of each neuron. Across all of the preparations,about 1/3 of neurons were found to predict synchronous bursts
robustly better than population mean field over some prediction
span, Δt (Fig. 3G). Note that the prediction by mean field is notequal to the average performance of single neurons. How asingle neuron can better predict the mean field than the mean
field itself does is a nontrivial issue.
The fact that a single neuron can compare or even outperformthe population mean at predicting the mean-field dynamics
may appear surprising, considering that the single-cell activities
fluctuate due to spiking much more than the mean-field activity(e.g., compare the black traces in Fig. 3 C and D). However, it isexplained by the delay-embedding theorems (29, 30) that the
global state information could be reconstructed from observationof time series in a single variable participating in the dynamicalsystem (Fig. S2). The present result suggests that such a mathematical property can be demonstrated in a biological system.
Importantly, to reconstruct the global state required, we need toobserve multiple time points rather than the momentary state in
the time series of the observed variable (which is why the theorem is called “delay-embedding” theorem). It suggests that the
information about the global state is carried by the temporal
patterns including multiple time points, rather than momentarysnapshots of activity in neurons.
To explore the importance of temporal patterns of activities in
predicting bursts compared with the momentary firing rates, we
compared the obtained prediction accuracy with an analogous
index based on momentary firing rates, xi,t or bt (the red traces inFig. 3 C and D, respectively; see also Methods). In single neurons,
the eCCM tended to provide more accurate predictions than the
momentary-activity-based method (P < 5 × 10−14 in all of theindividual preparations, sign test with paired samples; Fig. 3F).
This superiority of the eCCM to the momentary firing rate is not
likely to be due to the difference in the sampled spike numbers
because it remained even if we lengthened the bin size forcomputing activity so as to have the same total duration as usedin the eCCM analysis [P < 6 × 10−3 in all of the preparations
except for Chip 2427 (which had much lower firing rate thanothers), sign test with paired samples; Fig. 3F, Inset]. Finally,
when we used the momentary firing rates instead of the eCCM
based on temporal sequence, the fraction of cells whose burst
detectability outperformed that of the mean field reduced from
∼1/3 to ∼1/5 (Fig. 3G). Together, these results suggest that thetemporal patterns of neural activity, not only momentary activity,
are crucial for early detection of synchronous bursts.
Burst Predictability Reflects Network Structures Realized in Nonbursting
Periods. Next, we examined whether the predictability of synchro-nous bursts corresponds to any specific network structure definedby synaptic interactions. The structural correspondences of burst
predictability were investigated in three steps (Fig. 4): we firstreanalyzed the nonbursting spontaneous activities to infer the directed causal network structure; next, the inferred causal networks
were verified by comparing them to synaptic interactions measured
by an additional electrical stimulation experiment with the same
preparations; finally, we related each neuron’s network connectivity with its burst prediction accuracy.
In the causal network analysis, we applied a previously proposed method that is capable of detecting nonlinear coupling(31) (Methods). Theoretically, the method detects causal interactions including indirect ones (32). Applying this method to allof the neuron pairs showing nonbursting spontaneous activities
yielded a matrix that represents pairwise, directed interactions
within each preparation (Fig. 4A).
The relevance of the inferred causal couplings to synaptic interactions among neurons was verified using an independent
dataset in which synaptic interactions were directly measured (Fig.
4 B and C). The present CMOS-based MEA system allowed us to
stimulate the neural tissue local to given electrodes at submillisecond accuracy, directly activating a single neuron’s soma at a hightemporal precision with few artifact (26) (Fig. S3). A subset ofneurons was electrically stimulated while the evoked activities in
4 of 6 | www.pnas.org/cgi/doi/10.1073/pnas.1705981114

A

B

C

D

Fig. 4. Network structures in nonbursting periods explain the local burst
predictability. (A) Causal network analysis. The network depicts the directed
pairwise interactions among neurons. The filled (orange) and open (blue)
markers represent the putative excitatory and inhibitory cells, respectively.
The thickness of arrows indicates the absolute strength of causal coupling(“causality”). For clarity of illustration, only the top 5% strong couplings areshown, and the nodes are distributed by multidimensional scaling with the
distances defined to be inversely proportional to the causality. The figure
shows the result for Chip 1440. (B) Synaptically connected neuron pairs show
larger causal interaction. The bar labels show the time span to be used forthe causality analysis. Bars in dark- and light blue: average cross-embedding
values for cells that showed short-latent (dark, latency <10 ms) and long(dark, latency> 10 ms) stimulus-triggered spike increase. The average causality index within time span [−100, 100] ms dissociated the postsynaptic cells
from other cells (P < 0.005, Wilcoxon sign rank test, independent samples).
The error bar shows the SEM across neurons. (C) Identification of synaptic
connectivity by electrical stimulation experiment. The figure shows the peristimulus time histograms of three representative neurons (from top tobottom, putative direct, indirect, and nonpost cells, respectively). The asterisks indicate the earliest significant increase (P < 0.05) of spike count,compared with the spike-count distribution in the prestimulus period(from −500 to −100 ms). The figure shows the spike histogram smoothed with aboxcar kernel of 5-ms width for visualization purpose. (D) The relationshipbetween burst predictability and the interaction strength in different neuron
types (putative excitatory and inhibitory neurons). The neurons were classified
into two groups depending on the strength of input or output causal couplings:
the top half of neurons was labeled as “strong,” whereas the bottom half waslabeled as weak. The error bar shows the SEM across neurons.

other neurons were monitored simultaneously. According to theevoked response latencies, we identified the short-latency (<10 ms)
and long-latency (≥10 ms) downstream cells for each stimulated
neuron (Fig. 4C), between which the latter was expected more
likely to include multisynaptic interactions and thus weaker causal
effects. As expected, we confirmed that the causality inferred basedon nonbursting spontaneous data consistently reflected the identified short- and long-latency interactions (Fig. 4B). Specifically, theinferred causality differentiated the short- and long-latency downstream cells from the remaining population (P = 0.0001, sign test,
data pooled across all of the preparations). This result also provides
an empirical validation for the causality inference with nonlinear
state-space reconstructions such as CCM (31), by directly relating
the causality extracted from the dynamics to the one defined by the
manipulation experiment. Based on these confirmations, we used
the inferred causality as proxies for the effective synaptic interactions in the network during the nonbursting periods.
Finally, we compared these causal network structures and theburst predictabilities, finding that the burst predictability in eachTajima et al.

Tajima et al.

Methods
Cell Culturing, Recording, and Causal Network Analysis. All experimental
protocol was approved by the Committee on the Ethics of Animal Experiments at the Research Center for Advanced Science and Technology, theUniversity of Tokyo (Permit number: RAC130106). The cell culturing, recording, stimulation, and burst detection were based on the protocols
reported previously (26, 27). The causal networks in the nonbursting period
were estimated using the CCM (31). The details of protocols are described in
SI Methods.
Burst Predictability Measured by Event-Targeted CCM. The preset methodborrows the technique used in CCM to quantify an early warning of a specificglobal event, rather than the average interactions among the variables
throughout long observations. The method extends/integrates the original
CCM and the related methods by (i) designing the analytic protocol to assess
the predictability of a specific event (in our case, the spontaneous synchronous burst), (ii) relating the global state to local components while theoriginal protocol only quantified the relationships between local variables,and (iii) quantifying on the predictability of temporally distant global statebased on local variable dynamics. This is done by finding a mapping between
the system’s global state bt and each neuron xt−Δt including temporal lag Δt,instead of simultaneous prediction among individual neurons. For burst
detection, we set t as the time corresponding to the peak of global state
during a burst period t ∈ ft1 , . . . , tn g, where n is the number of burst trials in
the experiment. We first reconstruct the state-space trajectory for neuron x
and the normalized global state b in the (randomized) delay coordinates, xtd =(xt, xt-τ . . ., xt-(d-1)τ) and btd = (bt, bt-τ . . ., bt-(d-1)τ). It is mathematically shown thatboth of these trajectories reproduce the topological structure (the proximity ofdata points) of attractor dynamics that they participate, with which a (near-)
future state of each variable is accurately predicted from the current state of

PNAS Early Edition | 5 of 6

APPLIED
MATHEMATICS

Discussion
In this study, we have focused on the relationship between local
and global neural properties by analyzing how single-cell activity
predicts the population mean-field dynamics. This could also beinterpreted as an implementation of generalized Takensʼ theorem with general observation functions of a dynamical system(33), subject to limited data size and noise (34). The method
revealed the heterogeneity of neurons also in the predictabilityof upcoming occurrence of a synchronous burst of a neural
population, which is temporally distant on the order of 100 ms.
The bursts were predicted accurately by dynamics of some single
neurons. Remarkably, such neuron even outperformed the
population mean-field dynamics itself in terms of prediction
accuracy, particularly when we extracted the information embedded in the temporal patterns of neural activity rather than the
momentary firing rates. From causal network analysis combinedwith the electrical stimulation experiment, the heterogeneity of
burst predictability was explained by the structures of neuralnetworks: the “hub”-like neurons having stronger synaptic interactions with other cells can better predict the upcoming global
network burst than others. It means that, even without observing
any burst events or measuring burst predictabilities in eachneuron, we can predict to some extent which neurons should be
used to forecast the bursts by looking at the causal network. In
addition, the present results based on the cultured neurons indicate
the heterogeneous relationships between single neurons and population activity are basic properties observed in the nervous network that is self-organized without external stimulus. We also
found that the excitatory neurons in the network hub better predict
the bursts than the inhibitory neurons. This could reflect thefunctional asymmetry between the cell types: the activity of specificexcitatory neurons may be sufficient for a burst, thus making themgood predictors. In contrast, inactivity among inhibitory neurons
may be a necessary but not a sufficient condition for a bursting
event and thus less useful as predictors.
Classic models of the synchronized population burst assume
diverging interaction that propagates activity of a cell to theother cells (4, 12–14). A variant of this model features “avalanche”-like bursting sequences, which models synchronous burstwith noise being amplified and propagated through the network(9, 11, 35–39). These models, focusing mainly on the stochastic
aspects of burst initiation and propagation, have found successes
in explaining a variety of statistical properties of bursting neuralnetworks. On the other hand, it has been challenging to predict
the occurrence of synchronous burst, partially because they canbe triggered by complex and heterogeneous mechanisms (40).
The present study focuses on deterministic aspects of the neural
population dynamics, with which the burst occurrence is predictable to the extent based on a preceding single-neuron activity. The existence of specific neurons predicting the upcoming
burst is also consistent with the previous idea of “leader neurons” that fire at the beginning of bursts (8, 41). Notably, however, the predictor of the global burst is not identical to theincreased activities in those single neurons as expected from previous studies (15, 37, 41), because the magnitude of firing rate
alone is not sufficient to explain the highly accurate predictability
based on the eCCM analysis (Fig. 3E). Rather, the earliest presages are likely to be embedded in the dynamical patterns of thoseneurons. On the other hand, not everything is predictable in a
deterministic manner before the burst occurrence. For example,
the trial-to-trial fluctuation in burst size was not predicted frompreburst dynamics (ρ = −0.289 ∼ −0.077, P > 0.4, Spearman’s rankcorrelation between the odd and even trials), suggesting that the
burst size is determined by stochastic mechanisms (11, 39).
Finally, the predictability of upcoming burst occurrence basedon single-cell dynamics suggests a potentially effective methodology to capture the early warnings of population-state transition
in a variety of networks, including epileptic seizures in the brain(5, 42), rumor propagation in social networks (43–45), or epidemic spreading in human networks (46, 47). In clinical purposes, for example, accurate prediction and characterization ofpathologically synchronized neural firing is an important issue (5,42). Although nonlinear features of field potential sequence
have been proposed as useful markers of epileptic state (42, 48–51), the predictability of global dynamics based on single-neuron
activity has not been thoroughly studied. In particular, the present study suggests that single neuron could even outperform thedirect observation of global state itself in predicting the globalstate dynamics. In principle, the same framework could be applied to study other network events than the synchronous burst,
such as the changes in more complex correlation patterns amongneurons. The success of the present state-space reconstructionbased method implies a new approach to detect the early warningsof transitions in global-network state based on local rather thanglobal features. However, it should be noted with an emphasis thatthe in vivo neural dynamics are generally much more complex; notnecessarily all of the in vivo dynamics show the same properties as
what are observed in vitro. Nonetheless, our analytic protocol canbe applied to in vivo experiments with a wide variety of systems. If
some network events differed from others regarding the local early
warnings detected via state-space reconstructions, it can imply
potentially distinct properties of underlying dynamical systems.

NEUROSCIENCE

neuron displayed a positive correlation to the average inferredcausality between it and the other neurons (Fig. 4 D and E). We
also classified putative excitatory and inhibitory neurons basedon the spike waveforms (SI Methods, Fig. S4). Interestingly, theexcitatory neurons predicted bursts better than the inhibitory
neurons when we compared the neurons showing strong causality(P < 0.002, both for the causality for input and output connections, rank sum test), but such cell-type dependence was unclear
in the neurons having weak interactions (P > 0.5, Fig. 4D). These
results suggest that the stronger network interaction via excitatory cells underlies the higher burst predictability. The above
results together led us to conclude that the cell-to-cell variability
in burst predictability is likely to reflect the heterogeneous network structures that shape nonbursting activity dynamics.

other variables, if the dynamics (approximately) follows deterministic mechanisms (29–32). Practically, however, the accuracy of prediction is limited byvarious factors including data size, process/observation noises, as well as theuncertainty due to asymmetric causal interactions.
^ t based on a single neuron x with time
The forecast for the global state b
lag Δt is given by


^t xd
b
t−Δt =

X
t′s. t. x dt′−Δt ∈Bðx dt−Δt Þw x dt′−Δt − x dt−Δt  bt′ ,with the k-nearest-neighbor set Bðx dt−Δt Þ of x dt−Δt , in the same manner as in
CCM. In the data analysis we used k = 4 and the same weight function,w x ddt′−Δt − x t−Δt


 =


 
exp −x dt′−Δt − x dt−Δt 
.
  d
Σt′s.t.x d ∈Bðx d Þ exp −x t′−Δt − x dt−Δt 
t′−Δt
t−Δt

threshold used for defining the burst period (for details of burst period
definition, see SI Methods). The burst predictability in each cell was quantified by the rate of trials in which the upcoming bursts were successfully
predicted, where the success rate was summarized as an average over rangeof time span −200 ms < Δt < −50 ms.
Similarly, the forecast based on the history of the global state itself is givenby the same protocol except for using b instead of x:^ t bdbt−Δt =

X





w bdt′−Δt − bdt−Δt  bt′ .

t′s. t. bdt′−Δt ∈Bðbdt−Δt ÞThis quantifies the self-predictability of burst based on the global state itself,
through the success rate computed based on the corresponding threshold θb,^ t ðbd
the 95th percentile of null-hypothesis distribution Pðb
Null
tNull −Δt ÞÞ.

The success rate in burst detection was defined by the ratio between the^ t > θx , where θx is the criterion for detecting
number of burst trials such that b
an upcoming burst based on neuron x’s activity sequence. In the present
study, we defined θx by the 95th percentile in the null-hypothesis distribu^ t ðx dtion Pðb
Null
tNull −Δt ÞÞ, where tNull is randomly selected time outside thebursting period. Note that this detection criterion is different from the

ACKNOWLEDGMENTS. The authors thank Dr. Urs Frey and Prof. Andreas
Hierlemann for providing CMOS-based MEA with their technical support. This
work was supported by Japan Science and Technology Agency PRESTO Grant
JPMJPR16E6 (to S.T.); Asahi Glass Foundation (H.T.); Kayamori Foundation ofInformation and Science Advancement (H.T.); KAKENHI Grants 26242040,
16H01604, and 16K14191 (H.T.); Brain/MINDS from Japan Agency for Medical
Research and Development (T.T.); and RIKEN Brain Science Institute (T.T.).

1. Chialvo DR (2010) Emergent complex neural dynamics. Nat Phys 6:744–750.
2. Lisman JE (1997) Bursts as a unit of neural information: Making unreliable synapses
reliable. Trends Neurosci 20:38–43.
3. Izhikevich EM, Desai NS, Walcott EC, Hoppensteadt FC (2003) Bursts as a unit of neural
information: Selective communication via resonance. Trends Neurosci 26:161–167.
4. Traub RD, Wong RK (1982) Cellular mechanism of neuronal synchronization in epilepsy. Science 216:745–747.
5. McCormick DA, Contreras D (2001) On the cellular and network bases of epileptic
seizures. Annu Rev Physiol 63:815–846.
6. Staley KJ, Dudek FE (2006) Interictal spikes and epileptogenesis. Epilepsy Curr 6:
199–202.
7. Feinerman O, Segal M, Moses E (2007) Identification and dynamics of spontaneous
burst initiation zones in unidimensional neuronal cultures. J Neurophysiol 97:
2937–2948.
8. Eckmann JP, Jacobi S, Marom S, Moses E, Zbinden C (2008) Leader neurons in population bursts of 2D living neural networks. New J Phys 10:015011.
9. Yang H, Shew WL, Roy R, Plenz D (2012) Maximal variability of phase synchrony incortical networks with neuronal avalanches. J Neurosci 32:1061–1072.
10. Petermann T, et al. (2009) Spontaneous cortical activity in awake monkeys composedof neuronal avalanches. Proc Natl Acad Sci USA 106:15921–15926.
11. Beggs JM, Plenz D (2004) Neuronal avalanches are diverse and precise activity patterns that are stable for many hours in cortical slice cultures. J Neurosci 24:5216–5229.
12. Johnston D, Brown T (1981) Giant synaptic potential hypothesis for epileptiform activity. Science 211:294–297.
13. Bush PC, Douglas RJ (1991) Synchronization of bursting action potential discharge in amodel network of neocortical neurons. Neural Comput 3:19–30.
14. Pinsky PF, Rinzel J (1994) Intrinsic and network rhythmogenesis in a reduced Traubmodel for CA3 neurons. J Comput Neurosci 1:39–60.
15. Li C-YT, Poo M-M, Dan Y (2009) Burst spiking of a single cortical neuron modifies
global brain state. Science 324:643–646.
16. Miles R, Wong RK (1983) Single neurones can initiate synchronized population discharge in the hippocampus. Nature 306:371–373.
17. Takahashi H, et al. (2012) Light-addressed single-neuron stimulation in dissociatedneuronal cultures with sparse expression of ChR2. Biosystems 107:106–112.
18. Okun M, et al. (2015) Diverse coupling of neurons to populations in sensory cortex.
Nature 521:511–515.
19. Nowotny T, Rabinovich M (2007) Dynamical origin of independent spiking andbursting activity in neural microcircuits. Phys Rev Lett 98:128106.
20. Sugihara G, May RM (1990) Nonlinear forecasting as a way of distinguishing chaosfrom measurement error in time series. Nature 344:734–741.
21. Ye H, et al. (2015) Equation-free mechanistic ecosystem forecasting using empiricaldynamic modeling. Proc Natl Acad Sci USA 112:E1569–E1576.
22. Ye H, Sugihara G (2016) Information leverage in interconnected ecosystems: Overcoming the curse of dimensionality. Science 353:922–925.
23. Ye H, Deyle ER, Gilarranz LJ, Sugihara G (2015) Distinguishing time-delayed causal
interactions using convergent cross mapping. Sci Rep 5:14750.
24. Sugihara G, et al. (1999) Residual delay maps unveil global patterns of atmospheric nonlinearity and produce improved local forecasts. Proc Natl Acad Sci USA 96:14210–14215.
25. Sugihara G, Grenfell B, May RM (1990) Distinguishing error from chaos in ecologicaltime series. Philos Trans R Soc Lond B Biol Sci 330:235–251.
26. Bakkum DJ, et al. (2013) Tracking axonal action potential propagation on a highdensity microelectrode array across hundreds of sites. Nat Commun 4:2181.
27. Bakkum DJ, et al. (2013) Parameters for burst detection. Front Comput Neurosci 7:
193, 1–12.
28. Wagenaar D, Demarse TB, Potter SM (2005) MeaBench: A toolset for multi-electrodedata acquisition and on-line analysis. 2nd International IEEE EMBS Conference on
Neural Engineering (IEEE, Arlington, VA), pp 518–521.

29. Takens F (1981) Detecting strange attractors in fluid turbulence. Dynamical Systems
and Turbulence, Lecture Notes in Mathematics, eds Rand DA, Young L-S (Springer,
Berlin), pp 366–381.
30. Sauer T, Yorke JA, Casdagli M (1991) Embedology. J Stat Phys 65:579–616.
31. Sugihara G, et al. (2012) Detecting causality in complex ecosystems. Science 338:496–500.
32. Tajima S, Yanagawa T, Fujii N, Toyoizumi T (2015) Untangling brain-wide dynamics inconsciousness by cross-embedding. PLOS Comput Biol 11:e1004537.
33. Deyle ER, Sugihara G (2011) Generalized theorems for nonlinear state space reconstruction. PLoS One 6:e18295.
34. Casdagli M, Eubank S, Farmer JDD, Gibson J (1991) State space reconstruction in thepresence of noise. Physica D 51:52–98.
35. Orlandi JG, Soriano J, Alvarez-Lacalle E, Teller S, Casademunt J (2013) Noise focusing
and the emergence of coherent activity in neuronal cultures. Nat Phys 9:582–590.
36. Shew WL, Yang H, Yu S, Roy R, Plenz D (2011) Information capacity and transmission are
maximized in balanced cortical networks with neuronal avalanches. J Neurosci 31:55–63.
37. Hahn G, et al. (2010) Neuronal avalanches in spontaneous activity in vivo. J Neurophysiol
104:3312–3322.
38. Shriki O, et al. (2013) Neuronal avalanches in the resting MEG of the human brain.
J Neurosci 33:7079–7090.
39. Beggs JM, Plenz D (2003) Neuronal avalanches in neocortical circuits. J Neurosci 23:
11167–11177.
40. Pinto DJ, Patrick SL, Huang WC, Connors BW (2005) Initiation, propagation, andtermination of epileptiform activity in rodent neocortex in vitro involve distinct
mechanisms. J Neurosci 25:8131–8140.
41. Eytan D, Marom S (2006) Dynamics and effective topology underlying synchronization in networks of cortical neurons. J Neurosci 26:8465–8476.
42. Elger CE, Lehnertz K (1998) Seizure prediction by non-linear time series analysis ofbrain electrical activity. Eur J Neurosci 10:786–789.
43. Boyd D, Golder S, Lotan G (2010) Tweet, tweet, retweet: Conversational aspects ofretweeting on twitter. 43rd Hawaii International Conference on System Sciences(Institute of Electrical and Electronics Engineers, New York). Available at http://
ieeexplore.ieee.org/abstract/document/5428313/. Accessed January 4, 2013.
44. Ishii A, et al. (2012) The “hit” phenomenon: A mathematical model of human dynamics interactions as a stochastic process. New J Phys 14:63018.
45. Zhou J, Liu Z, Li B (2007) Influence of network structure on rumor propagation. Phys
Lett A 368:458–463.
46. Gross T, D’Lima CJD, Blasius B (2006) Epidemic dynamics on an adaptive network. Phys
Rev Lett 96:208701.
47. Pastor-Satorras R, Vespignani A (2001) Epidemic spreading in scale-free networks.
Phys Rev Lett 86:3200–3203.
48. Burns SP, et al. (2014) Network dynamics of the brain and influence of the epileptic
seizure onset zone. Proc Natl Acad Sci USA 111:E5321–E5330.
49. Wang Y, Goodfellow M, Taylor PN, Baier G (2012) Phase space approach for modelingof epileptic dynamics. Phys Rev E Stat Nonlin Soft Matter Phys 85:061918.
50. Stam CJ (2005) Nonlinear dynamical analysis of EEG and MEG: Review of an emerging
field. Clin Neurophysiol 116:2266–2301.
51. Lehnertz K, Elger CE (1995) Spatio-temporal dynamics of the primary epileptogenic area
in temporal lobe epilepsy characterized by neuronal complexity loss. Electroencephalogr
Clin Neurophysiol 95:108–117.
52. Frey U, et al. (2010) Switch-matrix-based high-density microelectrode array in CMOS
technology. IEEE J Solid-State Circuits 45:467–482.
53. Livi P, Heer F, Frey U, Bakkum DJ, Hierlemann A (2010) Compact voltage and current stimulation buffer for high-density microelectrode arrays. IEEE Trans Biomed Circuits Syst 4:372–378.
54. Lorenz EN (1969) Atmospheric predictability as revealed by naturally occurring analogues. J Atmos Sci 26:636–646.
55. Kennel MB, Abarbanel HD (2002) False neighbors and false strands: A reliable minimum embedding dimension algorithm. Phys Rev E Stat Nonlin Soft Matter Phys 66:
026209.

6 of 6 | www.pnas.org/cgi/doi/10.1073/pnas.1705981114

Tajima et al.

Supporting Information
Tajima et al. 10.1073/pnas.1705981114
SI Methods
Cell Culturing. All experimental protocol was approved by theCommittee on the Ethics of Animal Experiments at the ResearchCenter for Advanced Science and Technology, the University ofTokyo (Permit number: RAC130106).
A protocol for cell culturing had been reported previously (26).
Briefly, E18 Wistar rat cortices were dissociated using trypsin andmechanical trituration. Next, 20–40 k/μL neurons and glia wereseeded over an area of ∼12 mm2 on top of the CMOS chip. Layers ofpoly (ethyleneimine) followed by laminin were used to adhere cells.
Plating media consisted of Neurobasal-B27 supplemented with 10%
horse serum and 0.5 mM GlutaMAX during the first 24 h. Growth
media consisted of DMEM supplemented with 10% horse serum,
0.5 mM GlutaMAX, and 1 mM sodium pyruvate. Experiments wereconducted inside an incubator to control of environmental conditions(36 °C and 5% CO2).
CMOS-Based Recording and Stimulation of Network Activity. Cultured neuron activities were simultaneously recorded with highdensity MEA as described before (26, 27). Cortical networks weregrown over 11,011-electrode CMOS-based MEAs (52, 53), which
provide enough spatial and temporal resolution to detect action
potentials from any neuron lying on the array: 1.8 × 2.0-mm2 area
containing 8.2 × 5.8-μm2 electrodes with 17.8-μm pitch, sampled
at 20 kHz. Subsets of 126 electrodes can be read out (andstimulated) at one time, and electrode selection can be reconfigured within a few milliseconds. To identify the locations ofneurons growing over the array, a sequence of about 100 recording configurations were scanned across the whole array
while recording spontaneous activity. Locations of active somata
were identified because action potentials could usually be detected from multiple nearby electrodes. Electrode selection wasthen reconfigured such that a single electrode was assigned toeach identified soma, and spontaneous activities were measured
simultaneously from all of these cells. The putative neuron types(excitatory or inhibitory) were estimated based on spike shapes(Fig. S3). We could recode from 93, 47, 98, 92, and 53 neurons in
Chip 1437, 1440, 1444, 2427, and 2440, respectively.
Microstimulation-elicited network activities were then investigated to characterize a pairwise synaptic strength between
neurons. An adequate stimulating electrode was explored suchthat a single target neuron was directly activated through axonal
stimulation and that the directly evoked spikes were exclusively
measured from the target cell. The directly evoked spikes could be
easily distinguished from postsynaptic activations because they
were very reliable (i.e., 100 spikes elicited out of 100 stimulation
trials) and exhibited a small temporal jitter (Fig. S4). The microstimulation was a single, positive–negative, biphasic pulse with acharge-balanced amplitude of ±300–900 mV and a duration of
200 μs/phase, and was delivered 100 times every 3 s.
Burst Detection. Burst in spontaneous activity was detectedmodifying a protocol that was previously established by the authors’ group (27). The previous study proposes to determine thethreshold for burst detection based on the distribution of ISI of
consecutive spikes (27). In the bursting neurons, the ISI distribution typically has a bimodal structure whose valley can be used
as an objective criterion for burst detection. In this study, we firstcomputed sequences of population firing rate that were normalized such that it ranges from 0 to 1. We defined this as thesequence of global state bt, where t is the time. This was to use acommon burst-detection threshold across different preparations
Tajima et al. www.pnas.org/cgi/content/short/1705981114

of neuron cultures, which generally vary in terms of the absolute
firing rate. We used 5-ms bin for the firing rate computation. We
next derived the distributions of inverse firing rate over the bins.
This yielded an ISI distributions for each preparation, in which we
confirmed that all of them had bimodal shapes (Fig. 1D, Inset). For
each preparation, we selected the smallest ISI providing the valleyof distribution as the burst-detection threshold. The burst periods
were determined by the at least three consecutive sequences of binsthat have average ISIs under this threshold, identifying 26 bursts on
average for each preparation. We defined the peak burst timing byselecting the center of the bin having the smallest average ISI (i.e.,
the largest global state) within each burst period.
Estimation of Synaptic Connectivity from Electrical Stimulation.

Synaptic connectivity was estimated based on the evoked responses during electrical stimulation experiment. To reliably
stimulate a single neuron, we searched a stimulation site around atarget neuron in each MEA that could elicit an action potential
exclusively at the target neuron. After the careful selection of
stimulation sites, we could evoke the action potentials at almost
100% probability for each single stimulation, with very little jittersin the timings of action potentials in the stimulated neurons. We
first aligned the spike raster to the timings of stimulation. We nextcomputed the sequence of firing rate xi,t in a way described above.
The firing rates before the stimulation (−500 ms < t < −100 ms)
were used to produce the reference probability distribution ofeach neuron’s firing rate, P(xi). Next, we computed the probability, P(xi > xi,t), for each time bin centered at t in a shortpoststimulation period (2.5 ms < t < 15 ms); the responses
during 0 ms < t < 2.5 ms were omitted to eliminate the artifactsof electrical stimulation. We defined the smallest time t thatsatisfies P(xi > xi,t) < 0.01 in the poststimulation period as thelatency of neuron i. The neurons that had at least one time binsatisfying this condition were defined as “downstream” cells of
the stimulated neuron; the neurons that did not have any suchtime bin were defined as “non-downstream” cells.
According to the evoked response latencies, we identified theshort- (<10 ms) and long-latency (≥10 ms) downstream cells foreach stimulated neuron. The multisynaptic downstream cells fora stimulated neuron comprised relatively small fraction n (17%on average) of the entire population. This indicates that they aresubsets of neurons receiving effectively strong input from thestimulated cell via relatively a small number of path length, although the interaction with longer path length, which was
not detected here, could include the larger fraction of the cell
population.
Causal Network Analysis. Pairwise causal interaction among neurons was estimated based on a variant of CCM. CCM was developed by Sugihara et al. (31) as an extension of nonlinear
forecasting method of time sequence based on nearest-neighbormodels (20, 54). The method is capable of detecting relatively
weak causal interactions in deterministic dynamical systems(which can include some stochastic components) if the system andobservations are deterministic, the embedding dimension is
sufficiently large, and data size is sufficiently large for the given
embedding dimension. A variant of CCM that can be used forsystems in which the variables have heterogeneous time scales(such as in neural system) was developed by the authors (32).
Suppose that we want to know the interaction from neuron
x to another neuron y. We first reconstruct the state-space
trajectories of each neuron x in the delay coordinates,
1 of 4xdmax= ðxt , xt−τ , . . . , xt−ðdmax−1Þτ Þ, where dmax represents the maxt
imum number of dimensions (number of delay coordinates) tobe considered, t is the time point, and τ is the unit delay length.
We used τ = 5 ms dmax = 8 in this study. In the embedding-based
analyses, we convolved the firing-rate sequence of each neuron
with a Gaussian kernel that has half-width-of-half-height of
25 ms, and normalized such that each neuron trace ranges from0 to 1. To avoid a known vulnerability of the state-space reconstruction to the time-scale heterogeneity (55), we projected
delay vector xdmax
to a randomized coordinate space by multit
plying a square random matrix, R, from the left, to obtain atransformed vector: xdmax= R xdmax. A d-dimensional delay vectt
tor xdt is constructed by selecting the first dð≤dmax Þ components
of xdmax. The causal interaction form neuron x to another neuronty is quantified based on the correlation coefficient, ρyt ,^y ðxdt Þ, between the true (yt) and forecast (^yt ðxdt Þ) signals, where ^yt xdt =X
t′s. t. xdt′ ∈Bðxdt Þw xdt′ − xdt j yt′ ,with the k-nearest-neighbor set Bðxdt Þ of xdt in the delay-coordinate
space. We set k=4, weight wðjxdt′ − xdt jÞ = ½expð−jxdt′ − xdt jÞ=
Σt′s.t.xd ∈Bðxdt Þ expð−jxdt′ − xdt jÞ, and jxdt′ − xdt j as the square distancet′
between xdt′ and xdt in the data analysis. Note that ρyt ,^y ðxdt Þ = 1
means perfect prediction. We selected the optimal dimension d =d* so as to maximize the prediction performance, ρyt ,^y ðxdt Þ. Typicalvalue of d* was distributed from 2 to 4 in the current data.

Fig. S1. Population activity dynamics visualized with principal component analysis (PCA). The trajectories of population activity were plotted within the spaceof top three principal components (PCs one to three), by applying PCA to the multineuron time series for each preparation. The results for three representative
preparations are shown.

ac Delay-coordinates

Global stateb

Local dynamics

Fig. S2. Delay-embedding theorem (known as Takens’ theorem). Suppose that we have the time series of three variables, ðxðtÞ, yðtÞ, zðtÞÞ, generated by some
_ = f ðxðtÞ, yðtÞ, zðtÞÞ, yðtÞ_ = gðxðtÞ, yðtÞ, zðtÞÞ, and zðtÞ_ = hðxðtÞ, yðtÞ, zðtÞÞ. (A) The evolution of the global state ðxðtÞ, yðtÞ, zðtÞÞ is repdifferential equations:xðtÞ
resented by a trajectory in the space of those three variables. (B) Consider observing the temporal sequence of a single variable, e.g., xðtÞ. (C) The attractor
topology in the original global-state space is fully recovered in a delay coordinate of the observed variable, ðxðtÞ, xðt − τÞ, xðt − τÞÞ with an arbitrary unit delay, τ.
Namely, we can construct a smooth one-to-one map from the reconstructed attractor to the original attractor.

Tajima et al. www.pnas.org/cgi/content/short/1705981114

2 of 4

a
Pre

Post

b
100 µV
10 ms

Trial #c 100

1

Spike #

d 100

0

0

10

10
20
20
30 0
Time after stimulation [ms]30

Fig. S3. Microstimulation-based estimation of synaptic connectivity in a pairwise manner. (A) Axonal stimulation on an arbitrary neuron elicited bidirectional
action potential propagation. (B) Raw data of neural responses at a putative presynaptic neuron and postsynaptic neuron. Data from 100 trials are superimposed. (C) Raster plot of B. Antidromic, direct action potentials exhibited precise temporal responses, while orthodromic, synaptic action potentials wereelicited stochastically with significant temporal jitters. (D) Poststimulus spike histograms of C.

Tajima et al. www.pnas.org/cgi/content/short/1705981114

3 of 4aInhibitory

MAP2 / GABA

Excitatory

50µm

200 µV
1ms

Tpp

b 20
Count of neurons

16
12
8
4
0

0

0.2

0.4

0.6

0.8

1.0

1.2

1.4

1.6

Tpp [ms]
Fig. S4. Identification of excitatory and inhibitory neurons. (A) Representative neurons in immunostaining with MAP2 and GABA. Action potentials below theinsets were obtained at white rectangles, putatively from a neighboring neuron in a circle. The peak-to-peak time, Tpp, was defined as time duration fromnegative peak to positive peak of action potential. (B) Histogram of Tpp. Excitatory neurons (green) had larger Tpp than inhibitory neurons (magenta).
K-means method to Tpp was used to separate excitatory and inhibitory neurons.

Tajima et al. www.pnas.org/cgi/content/short/1705981114

4 of 4