arXiv:0912.0265v2 [] 22 Jan 2010

Mapping the spatiotemporal dynamics of
calcium signaling in cellular neural
networks using optical flow
Marius Buibas,

1

Diana Yu,1 Krystal Nizar,3 Gabriel A. Silva,1,2,3

1

3

Department of Bioengineering; 2 Department of Ophthalmology;
Neurosciences Graduate Program, University of California, San Diego
9500 Gilman Drive, La Jolla, CA 92093 USA

Abstract- An optical flow gradient algorithm was applied to spontaneously forming networks of neurons and glia in culture imaged by fluorescence optical microscopy in order
to map functional calcium signaling with single pixel resolution. Optical flow estimates
the direction and speed of motion of objects in an image between subsequent frames in
a recorded digital sequence of images (i.e. a movie). Computed vector field outputs by
the algorithm were able to track the spatiotemporal dynamics of calcium signaling patterns. We begin by briefly reviewing the mathematics of the optical flow algorithm, and
then describe how to solve for the displacement vectors and how to measure their reliability. We then compare computed flow vectors with manually estimated vectors for the
progression of a calcium signal recorded from representative astrocyte cultures. Finally,
we applied the algorithm to preparations of primary astrocytes and hippocampal neurons
and to the rMC-1 Muller glial cell line in order to illustrate the capability of the algorithm
for capturing different types of spatiotemporal calcium activity. We discuss the imaging
requirements, parameter selection and threshold selection for reliable measurements, and
offer perspectives on uses of the vector data.
Keywords- Optical flow; Calcium signaling; Calcium imaging; Neural circuits; Neural
networks; Neurons; Glial cells; Astrocytes

Address correspondance to Dr. Gabriel A. Silva, UC San Diego Jacobs Retina Center, 9415
Campus Point Drive, La Jolla, CA 92037-0946 USA. Electronic mail: gsilva@ucsd.edu

1

Buibas et. al.

Mapping neural cell networks with optical flow

Introduction
Calcium signaling is an intermediate step in many of the signaling pathways in neurons
and glial cells and is informative of functional neural activity. In neurons calcium signaling precedes sub threshold and threshold (i.e. action potential) changes in membrane
voltage, and can be used to infer electrophysiology from optical imaging [33, 7, 38, 35].
In astrocyte glial cells it underlies the mechanisms by which these cells communicate in
astrocyte networks and in bi-directional communication with neurons [1, 6, 30]. Relative
changes in cytosolic calcium concentration can be measured using different fluorescence
indicator dyes that can be imaged by optical microscopy in the visual light range, such as
bulk loaded AM esters and genetically encoded calcium indicators [29, 34]. The emitted
fluorescence of indicator dyes change as a function of the relative amount of free calcium
ions individual indicator molecules are able to interact with. Although the relationship between measured fluorescence signals and the calcium levels that produce them is complex
and non-linear, it is assumed that there exists a correlation between measured changes in
emitted fluorescence by indicator molecules and differing cytoslic calcium concentrations.
In this context, the measured fluorescence signal provides a valuable qualitative metric of
changing calcium levels that allow inferences of cell signaling and function. Throughout
the rest of this paper, we will use the terms “calcium signal” or “calcium fluorescence”
to mean a measured calcium indicator fluorescence signal that reflects a relative cytosolic
calcium concentration, as is routinely implied in the literature, even though in practicality
we never know the real, i.e. absolute, free ion concentration that gives rise to the measured
fluorescence signal.
The data collected by a typical experiment records qualitative movies of imaged changes
in calcium fluorescence intensity. One can visualize calcium transients and their relative
positions and durations, but there is no inherent quantitative analysis of the data by the
experiment itself that allows one to derive the dynamics that characterize such signaling
events. For example, things such as propagation speeds and directions (i.e. velocity),
the kinetics of measured waveforms, or analysis that depend on such properties, such as
identifying and mapping the signaling geometry of intercellular calcium waves in networks
of neurons or astrocytes. Measuring and tracing calcium (or other second messenger)
fluorescence signals quantitatively from recorded movies manually is a tedious and labor
intensive process for even small data sets, and involves comparing intensities at different
frames and locations in order to calculate speeds and directions. It is generally not possible
to do so for large data sets that encompass high spatial and temporal resolution detail or
large numbers of cells interacting in a circuit or network.
This can be addressed by analyzing experimental data with a filter algorithm called
optical flow, which can be used to derive quantitative measurements of observed spatiotemporal calcium signals imaged from fluorescence movies. The resultant vector data has a
variety of uses, ranging from deriving basic measurements of signal velocity and direction,
to characterizing and classifying spatiotemporal calcium dynamics between different ex2

Buibas et. al.

Mapping neural cell networks with optical flow

perimental conditions. Optical flow is an imaging technique (i.e. a filter) that calculates
a two-dimensional displacement field between two subsequent frames in a movie, based on
the local spatial and temporal gradients of the two images. The optical flow filter originated in the computer vision field, where it was designed to approximate object motion in
time-ordered image sequences for applications like stereo disparity measurements, motion
estimation, movie encoding and compression, and object segmentation [15]. The algorithm
uses a computed local spatial and temporal gradient to approximate a displacement or
flow vector at each pixel in the image. In both neurons and glial cells cytosolic calcium
concentration changes manifest themselves as transient responses with a rapid increase, i.e.
rising phase, followed by a kinetically slower decaying phase. This is because free calcium
is cytotoxic and therefore kept at nanomolar concentrations in the cytoplasm under normal
conditions. It is only transiently elevated followed quickly by its re-uptake or extrusion.
Temporal changes are typically coupled to spatial changes as a signal propagates through a
cell. Measured fluorescence changes then trace specific paths during periods of observation
(c.f. Fig. 1). Calcium transients start at a particular location, travel in some direction
at a specific speed and terminate at a different location. The typical kinetics of calcium
transients in neural cells are particularly well suited to the computational requirements of
the optical flow algorithm.
We have successfully applied optical flow to calcium fluorescence movies and obtained
displacement vectors that track the spatiotemporal progression of calcium signals. The
filter works for calcium fluorescence data because calcium signals exhibit both spatial
and temporal gradients. The computed vectors provide point estimates of the speed and
direction of signals. Optical flow is ultimately an imaging filter that works on whole
movies, much like edge filters and image segmentation filters are used in static microscopy
[12, 25, 13], and provides a novel and automated way of analyzing the spatiotemporal
dynamics of calcium intracellular signaling in neurons and astrocytes.
We begin by briefly reviewing the mathematics of the optical flow algorithm, describe
how to solve for the displacement vectors, and how to measure their reliability. We then
compare computed flow vectors with manually estimated vectors for the progression of a
calcium signal recorded from representative astrocyte cultures. Finally, we applied the algorithm to preparations of primary astrocytes and hippocampal neurons and to the rMC-1
Muller glial cell line in order to illustrate the capability of the algorithm for capturing
different types of spatiotemporal calcium activity. We discuss the imaging requirements,
parameter selection and threshold selection for reliable measurements, and offer perspectives on uses of the vector data.

Optical Flow Algorithm and Computation
In this section we briefly introduce the concepts and mathematics of optical flow, focusing
in particular on our own implementation of the algorithm to the experimental data that

3

Buibas et. al.

Mapping neural cell networks with optical flow

follows in the Results section. The theory behind the algorithm is well established and the
interested reader is referred to a number of excellent texts on the subject (see for example
[15, 16]). Optical flow is an algorithm that operates at the pixel level and calculates local
displacement or velocity between time ordered image pairs. Optical flow (or equivalently
image flow) is the perceived motion of an object in a field of view (e.g. by the human eye or a
camera), defined as the “flow” or change in space and time of gray values at the image plane.
It is an estimation of the motion field, which is the actual motion of the object in three
dimensional space projected onto the image plane (i.e. what we would like to know). As
long as the frequency of successive frames in an image sequence is shorter than the motion
or displacement of the object of interest (in order to avoid confounding ambiguities in
detecting the components of the motion caused by aperture and correspondence problemssee [15, 16]), the optical flow algorithm is able to track the motion of objects in the field
of view as a function of changing gray scale levels, subject to appropriate constraints and
minimizations. In other words, the algorithm assumes that any changes in gray values are
due to the object moving, and that the irradiance of the object is constant from frame
to frame. (This is actually a weak assumption that is difficult to satisfy since motion
usually causes changes in irradiance, which is why the algorithm is an estimation of the
motion field. In cases where irradiance does not change, the optical flow exactly equals the
motion field.) The algorithm assumes the conservation of gray levels in the field of view
and assumes that any changes in the distribution of gray levels are due to motion. In fact,
the optical flow constraint equation (introduced below) can be derived by analogy from
the continuity equation in fluid dynamics that conserves mass ([16]). By computing the
optical flow for all pixels in a field of view, displacement vectors can be calculated for each
pixel that map where an object moved to from the pixel in the first frame to that in the
second. Intuitively, one can see why the algorithm performs best with objects that have
strong contrasts at boundaries or large signal to background noise ratios. The kinetics of
calcium transient signals display clearly distinguishable rising and decaying phases that
trace specific paths during periods of observation in the form of intracellular calcium waves
(Fig. 1) that are readily detectable by the algorithm.
The underlying assumption for computation is to constrain local temporal gradients
to the product of spatial gradients and displacement vectors. The basic principle of the
algorithm takes as inputs two images and computes a vector for each corresponding pixel
in the images which approximates the displacement of a small window surrounding that
pixel between the two images (Fig. 2). Only intensity values inside the window are used for
computing the pixel displacement value, so the measurement is localized. Adjacent pixels
will have overlapping windows, so their computed vectors will be similar, much like pixels
in a blurred image are similar. Following a mathematical description of the algorithm we
describe the method for its solution and implementation that we used to derive the optical
flow for calcium signals. We also discuss parameters and constraints of relevance to calcium
fluorescence movies.
Consider an arbitrary pixel with gray level intensity I(x, y, t), displaced in the xy plane
4

Buibas et. al.

Mapping neural cell networks with optical flow

3.0s

4.0s

5.0s

5.5s

1000
750
500
250
0

40μm
0.5s

0.6s

0.7s

0.8s

100
80
60
40
20
0

100μm

Figure 1: Selected frames from recorded movies of imaged calcium fluorescence activity in
sparse networks of primary dissociated cortical astrocytes (top) and hippocampal neurons
(bottom- where a single neuron at high magnification is shown). The color coded scale
bars on the right represent fluorescence intensity I in units of ∆I/sec, as a first derivative
of the calcium signal. For neurons the signaling was spontaneous, while for the astrocytes
waves were pharmacologically induced. Fluorescence increases followed a relatively smooth
spatial progression across the frames at the times indicated by the time stamp in the upper
left hand corner of each image. Areas of increasing calcium concentration appear as positive
∆I/sec values, while areas of decreasing calcium concentration appear as negative values,
but at a much smaller magnitudes.

5

Buibas et. al.

Mapping neural cell networks with optical flow

40μm

Figure 2: The optical flow algorithm. A window Ω in the same location in two subsequent
image frames is used to compute a displacement or flow vector (arrow) for the pixel at
the center of the window. Only image intensity values in Ω are used for the calculation.
Vectors are computed for each pixel in an image frame except in border regions where Ω
falls outside of the image. Given the position of the pixel as (x, y) at t seconds, (x, y, t), the
displacement vector defines the motion of the pixel at the subsequent frame at δt seconds
as (x + δx, y + δy, t + δt).

6

Buibas et. al.

Mapping neural cell networks with optical flow

by δx and δy at time δt in an n x n window Ω (Fig. 2). This implies that
I(x, y, t) = I(x + δx, y + δy, t + δt)

(1)

A first order Taylor series approximation of I(x, y, t) by expansion of the right side of
1 results in
I(x + δx, y + δy, t + δt) = I(x, y, t) +

∂I
∂I
∂I
δx +
δy +
δt + higher order terms
∂x
∂y
∂t

(2)

Ignoring higher order terms, which provide negligible contributions, and taking into
consideration equation 1
∂I
∂I
∂I
δx +
δy +
δt = 0
(3)
∂x
∂y
∂t
Dividing by δt
∂I δx ∂I δy ∂I
+
+
=0
∂x δt
∂y δt
∂t

(4a)

∂I
∂I
∂I
ux +
ux +
=0
∂x
∂y
∂t

(4b)

The two spatial and one temporal gradients are defined by
δx
δt

∂I ∂I
∂x , ∂y ,

and

∂I
∂t ,

respec-

δy
δt

and vy =
represent the x and y spatial components of the optical flow
tively. ux =
displacement vector u(x, y) = (ux , vy ). The basic optical flow formulation is to constrain
temporal intensity changes (gradients) to the product of spatial gradients and u(x, y) to
give equation 4. In more compact notation this can be written as
∇I(x, y, t) · u(x, y) +

∂I(x, y, t)
=0
∂t

(5)

Computing optical flow means finding the values of u(x, y) at each location for every time
point that satisfy the above constraint, given the known local image intensity spatial and
temporal gradients.
Two factors establish computability of meaningful non-zero flow vector values. First,
local spatial gradients must be non-zero at the point of interest (x, y, t). There has to
be some image information around the pixel of interest, meaning that neighboring points
have to have different values so that gradients are non-zero. If all pixels in a window
around (x, y, t) have the same intensity values, then spatial gradients are zero and motion
is undetectable by any means. Second, for displacement between subsequent frames to
be computed, there has to be a temporal gradient at (x, y, t), or some change in intensity
between time points. If there is no temporal change in intensity between subsequent time
points, then a value of u(x, y) = 0 satisfies the constraint equation in 5. Both of these
requirements are limitations on the original application of the optical flow when estimating
displacement in natural scenes: objects may have constant intensity in a small window and
7

Buibas et. al.

Mapping neural cell networks with optical flow

still be moving, meaning that motion may occur and the recorded intensity spatial and
temporal gradients equal zero. These limitations are less important when optical flow is
applied to calcium fluorescence movies.
There are many methods for calculating optical flow for recorded movies (see [5] for a
review), and all of them work on digitized movies with discrete pixel values of position and
time, i.e. (x, y, t) ∈ (columns, rows, f rames). We chose the Lucas and Kanade method
because it is conceptually simple and efficient, and flexible in terms of the image processing
steps required for computation [3, 4, 22]. First, computation of the flow vector u(x, y) is
performed on a window or spatial neighborhood Ω of arbitrary size, centered around (x, y),
which is more reliable than a single point estimate at (x, y). Second, a window function
W (x, y) is defined to favor values in the center over those near the edges. The choice
of window size will depend on a variety of factors. It must be large enough to capture
the apparent displacement across frames and small enough to resolve features of interest.
The capture frame rate must be fast enough for displacements to be observable within the
width of the spatial observation window across successive frames. When measuring the
spatiotemporal motion of calcium signals the size of the window Ω, the frame rate, and the
resolution are all deeply tied to the size of the cells or cellular compartments in which the
signal travels. Together, these parameters must be chosen so that the signal is observable
and smooth enough to measure reliably as flow vectors across frames. For example, the
choice of parameter values used to image calcium signals in part of a dendrite or a fine
astrocyte process will necessarily be different than parameter values for broad calcium
signals that fill the soma.The constraint equation is redefined as a weighted least-squares
fit of local first-order constraints to a constant model of a local u(x, y) in each small spatial
neighborhood Ω around the pixel of interest. The goal is to find the value of u(x, y) that
minimizes

X
∂I(x, y, t) 2
W 2 (x, y) ∇I(x, y, t) · u(x, y) +
(6)
∂t
(x,y)∈Ω

The above equation can be rewritten and solved as the linear system:
AT W 2 A · u(x, y) = AT W 2 b

(7)

Where, for neighborhood Ω, consisting of n points centered around the pixel and time of
interest (x, y, t), Ω = {(x1 , y2 , t), (x2 , y2 , t), . . . , (xn , yn , t)}:

A=

∂I
∂x (x1 , y1 , t)

∂I
∂x (x2 , y2 , t)

∂I
∂y (x1 , y1 , t)

∂I
∂y (x2 , y2 , t)

...

...


W = diag W (x1 , y1 ), . . . , W (xn , yn )
iT
h ∂I
∂I
b=−
(x1 , y1 , t), . . . , (xn , yn , t)
∂t
∂t
8

∂I
∂x (xn , yn , t)
∂I
∂y (xn , yn , t)

T


(8)
(9)
(10)

Buibas et. al.

Mapping neural cell networks with optical flow

Ω is usually a square window with sizes typically ranging from 3 x 3 to 15 x 15 or n = 9
to n = 225 points. We have set the weight matrix W to a two dimensional Gaussian with
σ 2 equal to 1/6 of the window width. As an example, for a 5 x 5 or n = 25 point window:


1
6
13
6
1
 6 54 112 54 6 

1 
 13 112 230 112 13 
W =


1000 
6 54 112 54 6 
1
6
13
6
1
Here, the center values in W have a greater contribution to the calculation than the edge
values, favoring gradient values at the pixel of interest.
Solving for the flow vector u(x, y) in equation 7, yields:

−1
u(x, y) = AT W 2 A AT W 2 b

(11)

Equation 11 describes a linear system in matrix form, where the flow vector u at spatial and
time location (x, y, t) is solved from the quantities of A, W , and b, defined
 from the
 spatial
and temporal derivates of n points around (x, y, t). The 2 × 2 matrix AT W 2 A matrix
contains all the image spatial derivatives, and if those values are close to zero, the matrix is
poorly conditioned,
 and flow estimates become unreliable. Ensuring that both eigenvalues
of the AT W 2 A matrix are sufficiently large is a good way to ensure that the matrix is
well conditioned, since a measure of the conditioning number is the ratio of the largest to
the smallest eigenvalue [5]. While this is not the only way to compute conditioning, this
is the test we used for visualization and measurement reliability of computed vectors for
calcium fluorescence data (see Appendix A.3 for more information).
Spatial and temporal derivates were computed using 2 x 2 convolution kernel filters,
where the ∗∗ operator denotes 2-dimensional discrete convolution:


∂I(x, y, t)
1 −1 1
(12)
= I(x, y, t) ∗ ∗
∂x
4 −1 1


∂I(x, y, t)
1 −1 −1
= I(x, y, t) ∗ ∗
(13)
1
1
∂y
4



∂I(x, y, t)
1
1 1 1
=
I(x, y, t + ∆t) − I(x, y, t) ∗ ∗
(14)
∂t
∆t
4 1 1
Here ∆t represents the time between frames or the frame rate 1/∆t. Since the temporal
derivative calculated in 14 forms the basis for the b vector in 11, the frame rate has a linear
effect on the magnitude of the flow vector u.
Optical flow outputs a displacement vector in units of pixels, normalized to the time
difference between the two frames used for computation. When normalized, the vector
takes on velocity units of pixels per frame (for this reason it is called a flow vector).
9

Buibas et. al.

Mapping neural cell networks with optical flow

The conversion to physical units will depend on the spatial resolution of the camera and
microscope, typically expressed in microns per pixel, and the sampling rate for the movie
capture, expressed in frames per second. Spatial resolution is a function of the objectives
used as well as the resolution of the imager and any pixel binning used. The frame rate is
limited at the high end by the camera sampling rate, and at the low end by the minimum
exposure time required to capture a detectable intensity signal. The exposure time may be
reduced by increased gain or pixel binning, but those come at a cost of reduced resolution or
increased noise. The conversion between units of pixels/frame and units of microns/second
is straightforward:
microns
pixels f rames microns
=
·
·
(15)
second
f rame second
pixel
While the optical flow algorithm produces vectors in units of pixels/frame, the analysis
of the data in the Results section below have been converted into physical units of microns/second, using the resolution and capture rate of the recordings given the specifics of
our imaging system.

Results
Comparison between computed and manually estimated flow vectors
We manually estimated flow vectors for 12 images equivalent to 6 seconds of calcium signaling in primary dissociated spinal cord astrocyte cultures (orange arrows in Fig. 3), and
qualitatively compared them to computed optical flow vectors for the same data (green
arrows in Fig. 3; note that only reliable vectors are shown as determined by the eigenvalue
test- c.f. equation 11 and above discussion). Manual estimation required stepping through
frames and approximating roughly how a calcium signal progressed in time, which in this
experimental preparation included intercellular calcium waves that propagated through a
subset of the cell network. The manually traced signals were not the only ones observable in the small movie sequence used, but were chosen to illustrate four representative
signaling paths. Manual estimation was performed in two second intervals, estimating the
incremental spatial progression of a given calcium signal across four frames.
Estimation of the flow or displacement of a cell signal such as calcium between frames
manually like we did for the data in Fig. 3 is a very tedious and labor intensive process,
and can only realistically be done under very sparse conditions where the observer can
clearly delineate the flow of the signal visually. It is nearly impossible to do at the pixel
or small window level. In contrast, optical flow calculates a displacement vector for every
pixel in every frame, operating at a much finer scale and capturing much more detail than
is possible with manual estimates. Nonetheless, in Fig. 3 for the purpose of qualitatively
validating derived optical flow vectors to manually estimated ones, in both cases there was
a clear overlap in vector direction between manual and flow vectors. By contrast, there
were greater differences in vector magnitudes between manual and flow vectors, which is
10

Buibas et. al.

Mapping neural cell networks with optical flow

2.0s

0.0s

4.0s

6.0s

40+m

2.0  sec

4.0  sec

6.0  sec

Figure 3: Comparison between computed optical flow vectors (green) and manually estimated flow vectors (orange) for a primary dissociated spontaneously forming astrocyte
network in culture. While computed vectors were calculated for every pixel and every
frame, manual vectors were estimated every four frames and only trace a few selected signals. Only reliable optical flow vectors are shown, and only one in four vectors in both
horizontal and vertical directions are shown for clarity. Unlike the manual vectors, flow
vectors are only shown for the current frame. The top sequence of panels show vectors
overlaid on extracted frames from the actual movie at the indicated times for the entire
field of view. The bottom panels show the vectors in detail for the 2, 4, and 6 second frames
in order to more clearly assess the qualitative overlap between optical flow computed and
manually estimated results. For optical flow vectors (in green) only vectors that putatively
correspond to manual vectors (in orange) are shown, in contrast to the upper panels which
show all computed vectors (see text). See appendix for details regarding experimental
preparations, imaging, and parameters for calculation.

11

Buibas et. al.

Mapping neural cell networks with optical flow

consistent with the fact that the manual estimates spanned four frames while optical flow
vectors were calculated across adjacent frames. There is also temporal overlap between
the two cases, in the sense that similar displacements were estimated for the same frames
using both approaches. The eigenvalue threshold masked out unreliable vectors, and this
correlated well with calcium activity; only areas of spatial and temporal changes in the
movie produced reliable vectors as assessed visually, which is ultimately the most accurate
estimator of complex motions, but only if given the right conditions (e.g. conditions that
allow the human eye to separate motion). The optical flow algorithm however, is able to
provide reliable quantitative measurements of signaling dynamics at spatial and temporal
scales simply not measurable by qualitative visual inspection or manual estimations of the
data.

Optical flow characterization of intercellular signaling
We applied the optical flow algorithm to typical calcium fluorescence movies of spontaneously forming sparse networks of neural glial cells and neurons in culture, and looked
at the dynamics of intercellular signaling following pharmacological or mechanical stimulation. We purposely chose sparse networks because it facilitates the visual interpretation
of the entire resultant vector field, but the algorithm itself can operate on any data that
displays an appropriate signal. We recorded movies from the rMC-1 Muller glial like cell
line, which mechanistically displays calcium signaling similar to Muller retinal glial cells in
vivo [39], primary dissociated spinal cord astrocytes, and primary dissociated hippocampal
neurons. Intercellular calcium waves in rMC-1 cells and astrocytes were mechanically induced by gently poking an initial cell without penetrating the cell membrane, while calcium
waves in neuronal networks were induced by the localized pharmacological application of
glutamate to one or a small group of cells (see the appendix below for details about experimental preparations and imaging parameters). In particular, intercellular calcium waves
in astrocytes and related anatomically specialized macroglial cells such as Muller cells in
the neural retina or Bergman glia in the cerebellum have been known to occur under experimental conditions for several years now, and have recently been shown in vivo under
both physiological and pathophysiological conditions in different parts of the brain, mediated by intracellular calcium transients that induce paracrine signaling, primarily through
adenosine triphosphate (ATP) [19, 20, 14]. Astrocyte and related macroglial cells engage
in bi-directional chemical signaling with neurons and have the ability to modulate and
directly participate in information processing in the brain, which necessitates more than
just interactions between neurons and almost certainly involves astrocytes somehow. The
functional roles of glial intercellular calcium waves and their contributions to modulating
neuronal information are not yet known, and in fact the dynamics of these signaling events
and the conditions under which they occur are just beginning to be explored.
The key parameter for computing optical flow using the Lucas-Kanade method is the
window size Ω, specified as a square of a given width (see above). It defines the local
12

Buibas et. al.

Mapping neural cell networks with optical flow

neighborhood of pixels along a point of interest that is used to compute the spatial and
temporal gradients required for the calculation. Though not required for computation, a
minimum value for the eigenvalues for the matrix AT W 2 A should be specified to mask
out unreliable measurements. This ensures that only reliable displacement vectors are displayed and used for analysis. Since the intensity values are a function of the experimental
setup, microscope, and camera, the AT W 2 A matrix and its eigenvalues will scale accordingly. The selection of the eigenvalue threshold is thus arbitrary, much like the selections
of the camera gain, exposure time, and other imaging parameters are made to generate
easily visible intensity values (see Appendix A.3 for more information on selecting suitable
eigenvalue thresholds). Table 1 shows the window sizes, eigenvalue thresholds, and capture
frame rates used to calculate the vector fields shown in Fig 4. The displacement vectors
can be converted into velocity by equation 15. The original calcium fluorescence movies
and Matlab code written to implement the optical flow algorithm are freely available by
contacting the corresponding author.
Table 1: Image capture and optical flow parameters for shown figures
Parameter
Frame Capture
Rate (Hz)
Window Size
(pixels at
1.3µm/pixel)
Minimum
Eigenvalue (λ1 , λ2 )
greater than

rMC-1 Cells

Astrocytes

Hippocampal
Neurons

16.4

8

4

11

9

11

11

1.4

0.3

Neuronal cultures displayed derived optical flow vectors along processes as the calcium
signal propagated throughout the network. As expected, computed vectors and the resultant vector field followed the geometry of connected processes (i.e. axons and dendrites)
in the sparse network (Fig. 4a). The pattern of activation in this example proceeded
diagonally from the site of stimulation in the upper left hand corner of the field of view.
Some neurons activated at considerably longer times following the stimulus (i.e. out to
7 or 8 seconds) most likely due to recurrent feedback signaling in the network which can
last several seconds. Note how since only reliable vectors are plotted, as determined by
the eigenvalue test, there are spatial discontinuities in the temporal progression of mapped
signals, which reflect areas where the algorithm could not compute reliable vectors given
the measured data. This may be especially true at lower magnifications as in the example
13

Buibas et. al.

Mapping neural cell networks with optical flow

shown here for comparatively large fields of view that capture many cells. This represents a
challenging task for the algorithm. Nonetheless, both the spatial and temporal progression
of calcium signals are easily visible. The computed data, being in vector form, can complement existing methods like cross-correlation that use only cell body data to establish
relationships between cells for example.
Signal flow patterns were also computed for astrocyte and rMC-1 glial networks (Fig.
4b and c, respectively). Astrocyte signaling showed rapid burst like radial patterns that
was mostly complete by 2 seconds, with some smaller regions of cells activating later as far
out as 6-7 seconds. This is consistent with descriptions of intercellular calcium waves reported previously [8, 26, 2, 31]. rMC-1 cells showed qualitatively similar radial patterns of
activation, with signaling occurring within about 3 seconds following stimulation. However,
unlike the astrocyte response, where there was uniform signaling across the network near
the site of stimulation, rMC-1 cells showed more heterogeneity in spatial activation patterns, with distinct clusters of cells activating and spreading calcium waves. The distances
traveled by the waves in the rMC-1 example roughly agree with previous quantitative characterizations of calcium waves in similar preparations, on average displaying wave distances
of about 60 µm over the first two seconds or so and distances between 50-100 µm over about
4 seconds [39]. It is interesting however that the spatial progression of the calcium signals
in this example was not linear as a function of time, in the sense that cells roughly equidistant from the site of stimulation activated at different times, within about one second for
some versus as late as three seconds for others. The relationship and dynamics between
the spatial versus temporal properties of such waves are difficult at best and usually not
possible to determine by visual inspection of recorded movies alone, and are not captured
by calculations such as the one dimensional signaling speed of a progressing wave front.
Furthermore, speed and distance calculations of neuronal and glial signaling across networks of cells are usually corse approximations computed using low magnification movies
that provide a sufficiently large field of view. In contrast, optical flow provides reliable
single pixel vectors for any sized region of interest that represent very fine grain detailed
descriptions of calcium signal propagation difficult to achieve otherwise. For the astrocyte
data from Fig. 4b, Fig. 5 illustrates the distribution of signaling speeds in µm/second for
65 optical flow vectors for a small 10 x 10 pixel region equivalent to a 13 x 13 µm region
in the field of view (orange box in the figure). Any region size of interest anywhere in the
imaged field that might be of functional interest to the investigator can be similarly characterized. By way of rough comparison, optical flow calculated speeds for calcium signals
in computed window were distributed from 1-10 µm/sec, and are roughly similar to those
reported previously using more approximate methods, in the range of 5 to 10 µm/second
[8, 26, 2, 31]. The bimodal distribution in the figure reflects what is visually apparent
in the source movie: some of the areas in the orange pixel region exhibit spatiotemporal
displacement while others do not, indicating that calcium concentration changes propagate
in specific regions with specific patterns. Manual estimates from the literature typically
look at maximum propagation speeds, as seen in the second peak at about 9 µm/second.
14

Buibas et. al.

40+m/s

(a)

Mapping neural cell networks with optical flow

2.20s

2.40s

2.60s

8
7
6
5

40+m
3.60s

3.80s

4.20s

4
3
2
1

60+m/s

(b)

1.13s

1.38s

1.50s

8
7
6
5

40+m
2.38s

2.88s

3.00s

4
3
2
1

40+m/s

0.48s

0.60s

0.72s

(c)

3.2
2.8
2.4
2

40+m
0.84s

0.96s

1.08s

1.6
1.2
0.8
0.4

Figure 4: Computed optical flow vectors for induced calcium signals in spontaneously
forming in vitro networks of (a) primary hippocampal neurons, (b) primary spinal cord
astrocytes, and (c) the rMC-1 Muller glial-like cell line. Six frames from each representative
recorded movie are shown with the computed vector field superimposed at times indicated
by the time stamps in each frame (left set of six panels). Right panels: Composite temporal
projections of the entire movies. The vector fields show the full spatial progression for the
evolving calcium signals, with time (i.e. temporal progression) color coded by the color
map (in seconds). Plotting the vector fields in this way allows the full spatiotemporal
propagation of derived signals from entire 15
movies to be summarized in a single image.
This facilitates the qualitative visualization and identification of complex dynamic signaling
patterns that would be difficult to detect otherwise, such as for example by simply ”playing
back” the movie.

Buibas et. al.

Mapping neural cell networks with optical flow

20

N=65

Vector count

15
10
5
0

0
5
10
15
20
Signal Speed (microns/second)

Figure 5: Optical flow velocity magnitude distributions for the astrocyte data from Fig.
4b. The flow vector magnitudes for reliable measurements in a 10 x 10 pixel (13 x 13 µm)
region (orange square) are shown as histograms.

Discussion
We describe and show the application of optical flow gradient methods for identifying and
spatiotemporally mapping functional calcium signaling in networks of neurons and glia.
Although we focused on networks of cells here, the method can be equally applied to the
analysis of spatially detailed sub-cellular compartmentalized regions of interest, such as
dendrites or astrocyte processes. The method makes use of the spatial first derivative of
moving objects in a field of view, in this case changes in fluorescence levels of calcium
indicator dyes associated with the free concentration of intracellular calcium, to track their
motion between subsequent frames in an image sequence (i.e. a recorded movie). The
mathematical foundations of optical flow are well established and optical flow algorithms
have been used in a wide variety of fields including applications to cell and molecular biology to track the movement of proteins, vesicles, and even whole cells [24]. In neuroscience
and neural engineering it has been used in electromyography [18, 27] and sensory perception [21, 28], while clinically it has been used to detect seizures in neonatal infants [17],
among other applications. However, the method has not been previously applied to tracking and visualizing calcium signaling and deriving quantitative measurements of calcium
spatiotemporal changes that underlie intracellular and intercellular functional signaling in
neural cells.
Although in this paper we applied the optical flow algorithm to two dimensional flu16

Buibas et. al.

Mapping neural cell networks with optical flow

orescence movies, the algorithm itself can be readily applied to a recorded movie made
up of three dimensional stacks acquired using two-photon microscopy. Work by others is
pushing two photon imaging towards recording real time functional signaling from three
dimensional volumes of active cellular neural networks [11, 10]. If the sampling rate is sufficiently high, optical flow can be computed in three dimensions using a volume instead of
a square window around a pixel to generate a three dimensional displacement vector. The
same constraints on volume size, sampling, and vector reliability metrics in two dimensions
apply to the three dimensional case.
Optical flow methods produce a lot of data, generating a vector for every pixel in every
image pair computed, so further processing, rendering and visualization methods are key to
making quantitative comparisons between experimental setups. Statistical comparisons can
be made from vector values by comparing differences between selected regions in different
preparations; velocity averages for each region can be compared using statistical methods
such as means, standard deviations, and p-values. While vector values from adjacent pixels
are not statistically independent, averaged vector values for a given region of interest may
be used for statistical comparison with another, non-overlapping region.
Another potential use of the vectors is to classify spatiotemporal patterns. Similar to
using a scalar kernel filter to match an image pattern such as an edge or corner, vector
fields themselves can be filtered with a known vector kernel to match a pattern of interest.
This method is called Clifford convolution [9] and has been used to label physical flow
regimes in fluid dynamics applications. By designing a vector field filter and convolving
it with computed optical flow vectors, a scalar map identifying specific patterns of flow
associated with the saptiotemporal dynamics of the measured signal can be constructed in
order to classify regions exhibiting such patterns.
One of the most exciting potential uses of computed flow vectors is in functional network
reconstruction. Borrowing again from the field of fluid mechanics, a dynamic vector field
can be used to reconstruct the path of a hypothetical particle from a given starting point,
tracing out the path that a signal might take between cells, much like a particle in a dynamic
flow field [36, 37]. Geometrically mapped paths of measured signals that originate in an
activating cell and propagate through a network may be very useful for reconstructing
the dynamics of the network. This would complement existing network reconstruction
algorithms which typically rely on temporal data around fixed regions of interest.

Acknowledgments
This work was supported by grant RO1 NS054736 from the National Institute for Neurological Disorders and Stroke (NINDS) at the National Institutes of Health (NIH).

17

Buibas et. al.

Mapping neural cell networks with optical flow

References
[1] C Agulhon, J.-C Platel, B Kolomiets, V Forster, S Picaud, J Brocard, P Faure, and
P Brulet. Bioluminescent imaging of ca2+ activity reveals spatiotemporal dynamics in
glial networks of dark-adapted mouse retina. The Journal Of Physiology, 583(3):945–
958, Jul 2007.
[2] N Allen and B Barres. Signaling between glia and neurons: focus on synaptic plasticity.
Curr Opin Neurobiol, 15(5):542–548, Oct 2005.
[3] S Baker and I Matthews. Lucas-kanade 20 years on: A unifying framework. International Journal of Computer Vision, 56(3):221–255, Jan 2004.
[4] P Baraldi, A Sarti, C Lamberti, A Prandini, and F Sgallari. Evaluation of differential
optical flow techniques on synthesized echo images. Ieee T Bio-Med Eng, 43(3):259–
272, Jan 1996.
[5] JL Barron, DJ Fleet, and SS Beauchemin. Performance of optical flow techniques.
International Journal of Computer Vision, 12(1):43–77, Jan 1994.
[6] M Bennett. A quantitative model of purinergic junctional transmission of calcium
waves in astrocyte networks. Biophysical Journal, 89(4):2235–2250, Oct 2005.
[7] M Canepari, M Djurisic, and D Zecevic. Dendritic signals from rat hippocampal
ca1 pyramidal neurons during coincident pre- and post-synaptic activity: a combined
voltage- and calcium-imaging study. The Journal Of Physiology, 580(2):463–484, Dec
2006.
[8] A Cornell-Bell, S Finkbeiner, and M Cooper. Glutamate induces calcium waves in
cultured astrocytes: long-range glial signaling. Science, Jan 1990.
[9] J Ebling and G Scheuermann. Clifford convolution and pattern matching on vector
fields. Visualization, 2003. VIS 2003. IEEE, pages 193 – 200, Oct 2003.
[10] W Gobel and F Helmchen. In vivo calcium imaging of neural network function.
Physiology, 22(6):358–365, Dec 2007.
[11] Werner Göbel, Björn M Kampa, and Fritjof Helmchen. Imaging cellular network
dynamics in three dimensions using fast 3d laser scanning. Nat Meth, 4(1):73–79, Jan
2007.
[12] Yu Guo, Biao Gong, Sebrina Levesque, Thomas Manfredi, and Ying Sun. Automated
detection and delineation of mitochondria in electron micrographs of human skeletal
muscles. Microsc Res Tech, 63(3):133–9, Feb 2004.

18

Buibas et. al.

Mapping neural cell networks with optical flow

[13] Mahboubeh Hashemi, Marius Buibas, and Gabriel A Silva. Automated detection of
intercellular signaling in astrocyte networks using the converging squares algorithm.
Journal of Neuroscience Methods, 170(2):294–299, Jan 2008.
[14] T Hoogland, B Kuhn, W Göbel, and W Huang. Radially expanding transglial calcium
waves in the intact cerebellum. Proceedings of the National Academy of Sciences, Jan
2009.
[15] B Horn and B Schunck. Determining optical flow(distribution of apparent movement
velocities of image brightness . . . . Techniques and applications of image understanding,
21-23:319–331, Jan 1981.
[16] B Jahne. Digital Image Processing. Springer, ”sixth edition” edition, 2005.
[17] Nicolaos B Karayiannis, Guozhi Tao, James D Frost, Merrill S Wise, Richard A Hrachovy, and Eli M Mizrahi. Automated detection of videotaped neonatal seizures based
on motion segmentation methods. Clin Neurophysiol, 117(7):1585–94, Jul 2006.
[18] M-G Knuttinen, T B Parrish, C Weiss, K S LaBar, D R Gitelman, J M Power, M-M
Mesulam, and J F Disterhoft. Electromyography as a recording system for eyeblink
conditioning with functional magnetic resonance imaging. Neuroimage, 17(2):977–87,
Oct 2002.
[19] KV Kuchibhotla, CR Lattarulo, BT Hyman, and BJ Bacskai. Synchronous hyperactivity and intercellular calcium waves in astrocytes in alzheimer mice. Science Signaling,
323(5918):1211, 2009.
[20] Z. L Kurth-Nelson, A Mishra, and E. A Newman. Spontaneous glial calcium waves in
the retina develop over early adulthood. Journal of Neuroscience, 29(36):11339–11346,
Sep 2009.
[21] Keith Langley and Stephen J Anderson. Subtractive and divisive adaptation in visual
motion computations. Vision Res, 47(5):673–86, Mar 2007.
[22] S Lim, J Apostolopoulos, and A Gamal. Optical flow estimation using temporally
oversampled video. Image Processing, IEEE Transactions on, 14(8):1074 – 1087, Aug
2005.
[23] CL MacDonald, D Yu, M Buibas, and GA Silva. Diffusion modeling of atp signaling
suggests a partially regenerative mechanism underlies astrocyte intercellular calcium
waves. Front. Neuroeng., 1, 2008.
[24] Kota Miura. Tracking movement in cell biology. Adv Biochem Eng Biotechnol, 95:267–
95, Jan 2005.

19

Buibas et. al.

Mapping neural cell networks with optical flow

[25] Eran A Mukamel, Axel Nimmerjahn, and Mark J Schnitzer. Automated analysis of
cellular signals from large-scale calcium imaging data. Neuron, 63(6):747–60, Sep 2009.
[26] E A Newman and K R Zahs.
275(5301):844–7, Feb 1997.

Calcium waves in retinal glial cells.

Science,

[27] Nils Ostlund, Björn Gerdle, and J Stefan Karlsson. Location of innervation zone determined with multichannel surface electromyography using an optical flow technique.
Journal of electromyography and kinesiology : official journal of the International
Society of Electrophysiological Kinesiology, 17(5):549–55, Oct 2007.
[28] C Pagano and G Bingham. Comparing measures of monocular distance perception:
Verbal and reaching errors are not correlated. Journal of Experimental PsychologyHuman Perception, Jan 1998.
[29] R Madelaine Paredes, Julie C Etzler, Lora Talley Watts, Wei Zheng, and James D
Lechleiter. Chemical calcium indicators. Methods, 46(3):143–51, Nov 2008.
[30] Eliana Scemes and Christian Giaume. Astrocyte calcium waves: what they are and
what they do. Glia, 54(7):716–25, Nov 2006.
[31] Eliana Scemes and Christian Giaume. Astrocyte calcium waves: What they are and
what they do. Glia, 54(7):716–725, Nov 2006.
[32] G A Silva, C Feeney, L R Mills, and E Theriault. A novel and rapid method for
culturing pure rat spinal cord astrocytes on untreated glass. J Neurosci Methods,
80(1):75–9, Mar 1998.
[33] D Smetters, A Majewska, and R Yuste. Detecting action potentials in neuronal populations with calcium imaging. Methods, Jan 1999.
[34] L Tian and L Looger. Genetically encoded fluorescent sensors for studying healthy
and diseased nervous systems. Drug discovery today Disease models, 5(1):27–35, Jan
2008.
[35] J Vogelstein, B Watson, and A Packer. Spike inference from calcium imaging using
sequential monte carlo methods. Biophysical Journal, Jan 2009.
[36] D Weiskopf and G Erlebacher. 12 overview of flow visualization. The visualization
handbook, Jan 2005.
[37] D Weiskopf, F Schramm, G Erlebacher, and T Ertl. Particle and texture based spatiotemporal visualization of time-dependent vector fields. Visualization, 2005. VIS
05. IEEE, pages 639 – 646, Oct 2005.

20

Buibas et. al.

Mapping neural cell networks with optical flow

[38] E Yaksi and RW Friedrich. Reconstruction of firing rate changes across neuronal
populations by temporally deconvolved ca2+ imaging. Nat Meth, 3(5):377–383, Jan
2006.
[39] D Yu, M Buibas, SK Chow, IY Lee, Z Singer, and GA Silva. Characterization of
calcium-mediated intracellular and intercellular signaling in the rmc-1 glial cell line.
Cellular and Molecular Bioengineering, 2(1):144–155, 2009.

A
A.1

Appendix
Cell Preparations

rMC-1 glial cells and primary spinal cord astrocyte cultures (the latter dissected and grown
similar to previously described [32, 23] were grown to approximately 80% confluency and
washed twice with Kreb-HEPES buffer (KHB) solution (10 mM HEPES, 4.2 mM NaHCO3 ,
10 mM glucose, 1.18 mM MgSO4 , 7H2O, 1.18 mM KH2 PO4 , 4.69 mM KCL, 118 mM NaCl,
1.29 mM CaCl2 , pH 7.4) and incubated with 5µM Fluo-4AM in KHB for 1 hr at room
temperature. Excess dye was removed by washing twice with KHB and an additional
incubation of 30 min at room temperature was done to equilibrate intracellular dye concentration and ensure complete intracellular hydrolysis. For astrocytes and rMC-1 cells,
calcium transients were initiated by mechanical stimulation of a single cell using a (0.5µm
i.d.) micropipette tip (WPI Inc., Sarasota FL) mounted on a M325 Micrometer Slide Micromanipulator (WPI Inc., Sarasota FL). Comparable data were obtained using adenosine
triphosphate (ATP) pharmacological stimulation.
For hippocampal cultures, dissociated hippocampal neurons from timed-pregnant embryonic day 18 (E18) Sprague-Dawley rats were cultured on glass bottomed tissue culture
dishes coated with poly-D-lysine and laminin (BD Biosciences, San Jose, CA). Cultures
were plated at a cell density of 106 cells/3.8cm2 . Cultures were maintained at 37C in 5%
ambient CO2 . Plating media was composed of basal medial Eagle (Invitrogen, Carlsbad,
CA) with 1X Glutamax, 1000 U/mL penicillin and streptomycin sulfate, 5% FBS, and 1X
N2 supplement. Culture media consisted of Neurobasal (Invitrogen, Carlsbad, CA) with
1X Glutamax, 1000 U/mL penicillin and streptomycin sulfate, 20mM glucose, and 1X B27
supplement. Culture media was supplemented with 10uM Ara-C for 24 hrs at 1DIV to
inhibit overgrowth of glia. All imaging was performed on 3-5DIV.
Bulk loading of hippocampal cell cultures was accomplished via incubation in the dark,
at room temperature, for 30 min in 1µM of the fluorescent Calcium indicator Fluo-4AM in Krebs-HEPES buffer (10mM HEPES, 4.2 mM NaHCO3 , 10mM glucose, 1.18mM
MgSO4 ·7H2 O, 1.18mM KH2 PO4 , 4.69 mM KCl, 118mM NaCl, 1.29 mM CaCl2 , pH 7.4),
followed by 2x 5 min washes in Krebs-HEPES with 100µM sulfinpyrazone. Hydrolysis
was allowed to proceed for an additional 30 min. Stimulation of neurons with glucose was
performed by microinjection of 100uL of 10mM glutamate in PBS from a specified-side

21

Buibas et. al.

Mapping neural cell networks with optical flow

of the culture dish, well outside of the microscope field of view. The fluorescence signal
generated across the monolayer of cells was recorded for 10 sec prior to glutamate injection,
and for 120 sec following injection. Cultured neurons were incubated for 30 min prior to
imaging in Mg2+ -free PBS to induce the synchronization of calcium transients.

A.2

Imaging Setup

Visualization of calcium indicator dye fluorescence was achieved using a 488 nm (FITC)
filter on an Olympus IX81 inverted fluorescence confocal microscope (Olympus Optical,
Tokyo, Japan) that included epifluoresence, confocal, phase, brightfield, and Hoffman differential interference contrast (DIC) modalities. Real-time movie recordings of calcium
transient propagation were acquired with a Hamamatsu ORCA-ER digital camera (Hamamatsu Photonics K.K., Hamamatsu City, Japan) and Image-Pro Plus data acquisition and
morphometric software (version 5.1.0.20, Media Cybernetics, Inc., Silver Spring, MD) or
LabView custom written data acquisition software (ScopeController). All images were
captured with a 10X objective, using a 2x2 binning on the camera, for a resolution on
1.3µm/pixel and a total image size of 612x572 (camera’s maximum resolution is 1224x1144).
Images sampled at frequencies ranging from 2 to 16.4Hz, or 0.5sec to 0.06sec exposure time.

A.3

Reliable Vectors via the Eigenvalue Test

Recall that flow vectors are computed from the linear system in 11. This is a typical linear
system of the form:
M ·u=z
where u is the unknown and z and M are known quantities. The condition number of a
matrix simply describes how a small deviation in the known z translates to an error in
u. A high condition number means the matrix is ill-conditioned, meaning that a small
deviation in z leads to a large deviation in u, making that computation unreliable. One
way to compute the condition number of a matrix is to take the ratio of the largest to
smallest eigenvalue of that matrix:
κ(M ) =

λmax (M )
λmin (M )

Since M is a 2 × 2 matrix, it has 2 eigenvalues so ensuring that both are above a certain
value makes the condition κ value relatively low. The minimum threshold value depends
on the incoming intensity values.
Intensity readings from the CCD camera can take on any number of values, based on
the digitization (8-bit, 12-bit, 16-bit, for example), the exposure time, gain setting on the
camera, and above all the dye loading in the cell preparation. Typically, during observation,
the experimenter manually adjusts gain and exposure time to obtain reasonable intensity
values, typically in the middle of the digitization range.
22

Buibas et. al.

Mapping neural cell networks with optical flow

Eigenvalues for the AT W A used in flow vector calculation typically scale with the range
of recorded intensity values and are calculated for every pixel, producing an eigenvalue
image map. The values chosen in table 1 were manually chosen during examination of the
minimum eigenvalue image for a few representative frames, ensuring that they fell between
ares where we visually detected spatiotemporal changes in intensity and areas were we did
not detect such changes. This is the same process one would undertake when thresholding
a regular monochrome image for the counting of cells: the intensity threshold is set to a
value between the intensity of an area where there is a cell and an area where there is no
cell.

23

