Imaging Brain Activation Streams from Optical Flow Computation on 2-Riemannian Manifolds Julien Lefèvre1 , Guillaume Obozinski2 , and Sylvain Baillet1 1

Cognitive Neuroscience and Brain Imaging Laboratory, CNRS UPR640–LENA, Université Pierre et Marie Curie–Paris6, Paris, F-75013, France {julien.lefevre,sylvain.baillet}@chups.jussieu.fr, 2 Department of Statistics, University of California, Berkeley, CA 94720-3860, USA .

Abstract. We report on mathematical methods for the exploration of spatiotemporal dynamics of Magneto- and Electro-Encephalography (MEG / EEG) surface data and/or of the corresponding brain activity at the cortical level, with high temporal resolution. In this regard, we describe how the framework and numerical computation of the optical flow — a classical tool for motion analysis in computer vision — can be extended to non-flat 2-dimensional surfaces such as the scalp and the cortical mantle. We prove the concept and mathematical well-posedness of such an extension through regularizing constraints on the estimated velocity field, and discuss the quantitative evaluation of the optical flow. The method is illustrated by simulations and analysis of brain image sequences from a ball-catching paradigm.

1

Introduction

Magnetoencephalography (MEG) and electroencephalography (EEG) respectively measure magnetic fields and electric potentials on the scalp surface, which provides investigation of neural processes with exquisite time resolution within the millisecond range. Estimation techniques of MEG/EEG generators have been considerably improving recently [5], thereby making electromagnetic brain mapping become a true functional imaging modality. However, most studies report on the classical brain mapping questions of ‘Where?’ and When?’ specific brain processes have occurred, but rarely address ‘How?’ these latter might be embedded in space and time altogether. The theory of electrographic objects [18] was the first computational approach with the objective of deciphering the ‘Rosetta Stone’ of brain language at the macroscopic scale. The concept of spatiotemporal elements, that were assumed to structurally sustain neural activity, was later reformulated by Lehmann in terms of microstates as ‘building blocks of mentation’ [14]. This concept emerged from empirical observations of the time-evolving topographies of scalp potential maps in EEG, which could be described as the succession of episodes of relativelystable spatial configurations. Multiple techniques have been proposed to detect

II

and characterize microstates and they have been applied in a great number of paradigms, from cognitive to clinical experimental neuroscience (e.g. [22]). However, they all refer to some static and geometrical analysis of EEG scalp patterns — through PCA analysis of time segments [17], dissimilarity measures on successive EEG topographies [24] — and have never been adapted to the analysis of brain image sequences so far. Our contribution suggests a new framework to investigate the spatiotemporal dynamics of brain activations in terms of the estimation and analysis of their displacement field. We first prove in Section 2 the concept and work at the theoretical aspects of the computation of optical flow on a 2-dimensional surface, through the generalization of existing variational formulations. In Section 3, we evaluate the consistency of the estimates and run realistic numerical simulations. Finally, the method is illustrated from experimental data in the context of a ball-catching experimental investigation with MEG.

2

Velocity field of neural activity

In multiple applications, analysis of dynamical phenomena through the computation of a velocity vector field has contributed to the description and extraction of informational contents about the processes involved (see e.g. [7],[9]). Such approaches have barely been suggested for the analysis of structured patterns within brain signals and image sequences with high temporal resolution (e.g. using MEG/EEG or optical imaging). In [13], estimation of velocity fields was restricted to 2D images of narrow-band scalp EEG measures in the α range (i.e. typically within [7, 14] Hz), with limited quantitative analysis. Our approach to the computation of the velocity vector field descends from optical flow techniques as introduced originally by Horn & Schunk [12]. These techniques have been demonstrating efficiency in the analysis of video sequences for about the last two decades (see e.g. [8, 15] which review a selection of computational methods associated with sound performance evaluation). The computation of vector flow is generally driven by basic assumptions which postulate conservation of brightness of moving objects. These restrictive hypotheses may not fit rigourously the exact nature of phenomena but have proven to yield commensurate estimations of vector fields provided they are valid locally in time and space [8]. The exquisite time-resolution of MEG/EEG images is generally compatible with these assumptions, as brain responses unfold to a large extent with substantial spatiotemporal smoothness. In the context of brain imaging though, we are facing the issue of distributed intensity variations in 3D. Detection can be restricted to the cortical surface as a first approach, hence recent surface flattening tools could be applied [23] prior to using classical 2D flow estimators. However, such a transform entails multiple kinds of limitations: the required topological cuts in the closed surface of the brain induce linking problems at boundaries; moreover, local distortions

III

of angles and distances are problematic when it comes to estimate the local orientation of the flow. Here we introduce a formalism based on differential geometry to extend the computation of optical flow equations on non-flat surfaces (see [10] for an introduction), hence on the folded geometry of the brain. 2.1

Optical flow on a non-flat domain

Let us define M, a 2-Riemannian manifold representing the imaging support (i.e. the scalp or cortical surfaces), parameterized by the local coordinates system φ : p ∈ M 7→ (x1 , x2 ) ∈ R2 . We introduce a scalar quantity defined in time on a 2-dimensional surface — e.g. brain activity from scalp MEG/EEG topographies or cortical activation maps — as a function I(p, t) ∈ M, where (p, t) ∈ M × R. We note eα = ∂xα p := ∂α p, the canonical basis of the tangent space Tp M at a S point p of the manifold, and T M = p Tp M the tangent bundle of M. M is equipped with a Riemannian metric, meaning that at each point p of manifold M, there exists a positive-definite form: gp : Tp M × Tp M → R, ¡ ¢ which is differentiable with respect to p. We later denote (gp )α,β = gp eα , eβ . A natural choice for gp is the restriction of the Euclidian metric to Tp M, which we have adopted for subsequent computations. For concision purposes, we will now only refer to gp as g. As in classical computation approaches to optical flow, we now assume that the activity of a point moving on a curve c(t) in M is constant along time. The condition dI = 0 yields : ∂t I + Dc(t) I(˙c) = 0, (1) where Dc(t) I is the differential of I at point c(t) applied to c˙ = V = (V 1 , V 2 ), the unknown vector field. We express the linear application Dc(t) I as a scalar product and introduce ∇M I, the gradient of I which is defined as the vector field satisfying at each point p : ¡ ¢ ∀V ∈ Tp M, g(∇M I, V) = Dp I V . (1) can thereby be transformed into an optical-flow type of equation: ∂t I + g(V, ∇M I) = 0.

(2)

We note that (2) takes the same form as general conservation laws defined on manifolds in [19]. Here, only the component of the flow V in the direction of the gradient is accessible to estimation. This corresponds to the well-known aperture problem [12], which requires additional constraints on the flow to yield a unique solution. This approach classically reduces to minimizing an energy functional such as in [12]: ¶2 Z Z µ ∂I + g(V, ∇M I) dµM + λ C(V)dµM , (3) E(V) = ∂t M M

IV

p where dµM is a volume form on M, for which we suggest det(gα,β )dx1 dx2 as a natural choice. The first term is a measure of fit of the optical flow model to the data, while the second one acts as a spatial regularizer of the flow. The scalar parameter λ tunes the respective contribution of these two terms in the net energy cost E(V). Here we use the smoothness term from [12], which can be expressed as a Frobenius norm: C(V) = Tr(t ∇V.∇V), (4) where

X ¡ ¢β ∇V α = ∂α V β + Γβαγ V γ γ

is the covariant derivative of V, a generalization of vectorial ∂α V β is P βgradient. γ the classical Euclidian expression of the gradient, and γ Γαγ V reflects local deformations of the tangent space basis since the Christoffel symbols Γβαγ are the coordinates of ∂β eα along eγ . This rather complex expression ensures the tensoriality property of V, i.e. invariance with parametrization changes. This constraint will tend to generate a regularized vector field with small spatial derivatives, that is a field with weak local variations. Such a regularization scheme may be problematic in situations where spatial discontinuities occur in the image sequences. For example, in the case of a moving object on a static background, the severe velocity discontinuities around the object contours are eventually blurred in the regularized flow field. In the case of brain activations revealed by MEG/EEG though, spatial patterns are naturally smooth thus we adopt the basic Horn & Schunk regularization scheme (see [25] for a taxonomy of other possible terms). 2.2

Variational formulation

Variational formulation of 2D-optical flow equation has been first proposed by Schnörr in [20]. The advantage of such formulation is twofold. Theoretically, it ensures the problem is well-posed; that is there exists a unique solution in a specific and convenient function space e.g. a Sobolev spacen [20], or a space of functions with bounded variations [4]. Numerically, it allows to solve the problem on discrete irregular surfacic tessellations and to yield discrete solutions belonging to the chosen function space. We demonstrate these assertions in the case of Horn & Schunk isotropic smoothness priors, but the general framework remains the same for Nagel’s anisotropic image-driven regularization approach [16]. Considering M, we need to define a working space of vector fields Γ 1 (M) on which functional E(V) will be minimized. Let us first denote the Sobolev space H 1 (M) defined in [11] as the completion of C 1 (M) (the space of differentiable functions on the manifold) with respect to k . kH 1 derived from the following scalar product : Z Z < u, v >H 1 = uv dµM + g(∇u, ∇v) dµM . M

M

V

We chose a space of vector fields in which the coordinates of each element are located in a classical Sobolev space: o n P2 Γ 1 (M) = V : M → T M / V = α=1 V α eα , V α ∈ H 1 (M) , with the scalar product : Z < U, V >Γ 1 (M) =

Z

M

g(U, V) dµM +

M

Tr(t ∇U∇V) dµM .

E(V) can be simplified from (3) as a combination of the following constant, linear and bilinear forms : Z ¡ ¢2 K(t) = ∂t I dµM , M Z f (U) = − g(U, ∇M I)∂t I dµM , M Z Z a(U, V) = g(U, ∇M I)g(V, ∇M I)dµM +λ Tr(t ∇U∇V) dµM . M M | {z } | {z } Fit to data

Regularizing term

Minimizing E(V) on Γ 1 (M) is then equivalent to the following problem : ¡ ¢ min a(V, V) − 2f (V) + K(t) . (5) 1 V∈Γ (M)

Lax-Milgram theorem ensures unicity of the solution with the following assumptions [3]: 1. a and f are continuous forms; 2. Γ 1 (M) is complete, the bilinear form a(., .) is symmetric and coercive (elliptic), that is, there exists a constant C such that : ∀ V ∈ Γ 1 (M), a(V, V) ≥ C k V k2Γ 1 (M) . Moreover, the solution V to (5) satisfies: a(V, U) = f (U), ∀ U ∈ Γ 1 (M).

(6)

Continuity of f and a are straightforward. Completeness of Γ 1 (M) is ensured because any Cauchy sequence has components in H 1 (M) which are also Cauchy sequences since k . kH 1 is bounded by k . kΓ 1 (M) . Proof of coercivity can be adapted — analogously to flat domains [20] — thanks to isothermal coordinates. Indeed, the Korn–Lichtenstein theorem (1914) allows to find a system of coordinates for which the two basis vectors of tangent space are orthogonal. In this basis, calculus are similar to those in Euclidian case by introducing a multiplicative coefficient equal to the norm of the basis vectors. The coercivity and therefore well-posedness requires only a similar assumption about linear independency of the two components of the gradient ∇M I.

VI

2.3

FEM computation of the optical flow

Now that we have proven the well-posedness of the regularized optical flow problem on a manifold M, we derive computational methods from the variational formulation, which are defined on a tessellation Mh which approximates the manifold. Let us consider the vector space of continuous piecewise affine vector fields on Mh which belong to the tangent space at each node of the tessellation. A convenient basis is: ¡ ¢ © ª wα,i = w(i)eα (i) for 1 ≤ i ≤ Card Mh , α ∈ 1, 2 , where w(i) stands for the continuous piecewise affine function which is 1 at node i and 0 at all other triangle nodes, and eα (i) is a basis of tangent space at node i. The variational formulation in (6) yields a classical linear system: © ª ∀j ∀β ∈ 1, 2 ,

Card(Mh )

X

2 X

i=1

α=1

a(wα,i , wβ,j )Vα,i = f (wβ,j ),

(7)

where Vα,i are the components of the velocity field V in the basis wα,i . Note that a(wα,i , wβ,j ) and f (wβ,j ) can be explicitly computed with first-order finite elements by estimating the integrals on each triangle T of the tessellation and summing the different contributions. Practically, ∇M I is obtained on each triangle T = [i, j, k] from the linear interpolation: ∇M I ∼ I(i)∇T w(i) + I(j)∇T w(j) + I(k)∇T w(k)., with ∇T w(i) =

hi , k hi k2

where hi is the height of triangle T from vertex i. Let us define Pn (i) as the projection operator onto the local tangent plane, which is obtained at node i by estimating the local normal n as the sum of normals of each triangle containing i. For each i, eα (i) is chosen as a basis of the kernel of Pn (i). The Christoffel symbols Γkij vanish in our discretization with first-order finite elements since the eα (i) have no variations on each triangle.

3

Simulations

We address the quantitative and qualitative evaluation of optical flow with simple and illustrative simulations on a selection of surfaces.

VII

3.1

Angular error estimation via parallel transport

Most of optical flow techniques are evaluated with quantitative comparison beˆ This criterion tween the true velocity field V and the estimated optical flow V. can be reduced as the evaluation of an angular error (AE) between two vectors, taking into account theirprelative amplitude. We introduce ||.|| = g(., .) + 1 and adapt the angular error estimate originally introduced in [6]: Ã ! ˆ +1 g(V, V) ˆ AE(V, V) = arccos . ˆ ||V||.||V|| Nevertheless, there remains the nontrivial issue of defining a velocity field on a manifold. During the uniform translation of a uniformly-colored patch for instance, the velocity vectors of moving points are not parallel in the Euclidian sense anymore. Hence we extend the definition of parallelism via the notion of parallel transport. Let us consider the centroid G of a translating illuminated patch with given velocity VG . In order to calculate the velocity at any point A of the patch, we transport VG along γ, a geodesic curve joining G to A (γ(0) = G, γ(1) = A). Practically we have to solve the differential equation : ∇γ(u) Y(u) = 0, ˙ with ∇γ(u) Y(u) = ˙

X

¡ ¢j γ(t) ˙ i ∇Y i , Y(0) = VG .

i

The velocity at point A is obtained taking Y(1). In the case of spherical geometry, the parallel transport along a geodesic reduces to a simple rotation whose axis is orthogonal to the plane containing the geodesic3 (Fig. 1). The implementation of general parallel transport has not been addressed since it reaches far beyond the scope of this study. The angular error provides a simple evaluation index as well as a quantitative criterion to adjust the regularizing parameter λ in (3), which was fixed to λ = 0.1 in the rest of the study. 3.2

Results on synthetic data

Two types of synthetic data were created. They approximate typical situations encountered when dealing with EEG or MEG images evolving on the surface of the head or on the cortical surface. The first set of simulations illustrates the emergence and fading of activity within on a single region of the brain (see Fig. 2). Even if this situation infringes the hypothesis of intensity conservation across time 1, the radial structure of optical flow gives an indication about how the system is evolving, spreading then collapsing. 3

Note the parallel transport cannot be defined for diametrally-opposite points.

VIII

Fig. 1. Evaluation of computation of the vector flow: a patch with Gaussian distribution of intensity travels on a sphere. The outcome of the optical flow computation is shown in green and compared with the true velocity field in red at two time instants. Mean angular error AE is given in degrees.

The second kind of synthetic data simulates propagations of activity across distant brain regions. Fig. 1 shows the translation of a Gaussian patch of activity and compares the true displacement field with its optical-flow estimation. At each step of time, we indicate the mean angular error as defined in Section 3.1 for each point of the patch. A similar type of unfolding patterns of activation is shown Fig. 3 on real brain surface geometry. An approximation of a Gaussian patch propagates from a rather flat domain of the cortical manifold and travels down into a sulcal fold. We represent the velocity of the patch centroid and the mean optical flow projected on a plane containing the true displacement. We speculate the irregular angular error could be improved with the method exposed in Section 3.1.

4 4.1

Application to the investigation of spatiotemporal dynamics of MEG signals The electromagnetic brain imaging problem

Magneto and electroencephalography stem from similar physical principles since they are directly related to the electromagnetical activity of neurons. Magnetic fields (and similarly electric potentials) are sampled on s sensors, M1 , . . . , Ms , as a linear combination of p sources signals S1 , . . . , Sp which can be written as M = AS, where A is the gain matrix from the forward model.

IX

Fig. 2. Simulating local emergence and waning of brain activity with resulting flow. Top row illustrates the progressive emergence followed by fading of a 33 cm2 region in the posterior medio-temporal brain area. The entire process unfolds within 100 time samples. Bottom row displays the corresponding distributions of vector flow; initially diverging from (left) then converging to (right) the center of the activation zone.

Electromagnetic brain mapping consists of the estimation of MEG/EEG sources S from scalp measures M . However, this inverse problem is underdetermined since there are far more possible cortical sources than sensors. s is typically on the order of a few hundreds while p amounts to about 10000 elemental sources constrained onto the surface of the cortical manifold, which was extracted from MRI image sequences [1]. Inverse modelling can therefore be approached as in many other image reconstruction applications, i.e. though the introduction of priors in addition to data. Here we used a weighted-minimum norm estimate (WMNE) of source amplitudes and its implementation in the BrainStorm software [2]. 4.2

Evaluation on experimental data

We applied optical flow computation to magnetic evoked-fields in a ball-catching paradigm [21]. The subjects had to catch a free-falling tennis ball which fall was initiated at time t = 0 ms. The second experimental condition (‘NoCatch’) consisted for the subject in only looking at the ball falling without catching it. Fig. 4 shows a rapid decrease of cerebellar activity around 40 ms after the ball hits the subject’s hand (371 ± 7 ms). We can notice the convergent and

X

Fig. 3. Simulation of brain activation propagating at the surface of the cortex. Top row, from left to right: displacement field of an activation patch is translating along a predefined path during 100 steps of time. Bottom row: the mean vector field and true displacement are shown in green and blue, respectively, with indication of instantaneous angular error. For clarity purposes, only one brain hemisphere is shown within the scalp surface.

radial structure of the flow. This phenomenon was not found in the NoCatch condition where no motor program was required at the time of ball-impact.

5

Conclusion

This article introduced an extension of computational approaches to optical flow on non-flat domains. The framework of Riemannian geometry allows to adapt a variational formulation for this ill-posed problem and to derive evaluation tools. We suggest new applications of the quantization of the displacement field to spatiotemporal data in MEG and EEG and to the question of neural information directionality. Results from simulations indicate the flow has satisfactory behavior in terms of spatial structure and angular errors for the application in question. Encouraging preliminary results have been presented from real experimental data. Ongoing developments consist in relating measures from this computational approach to their physiological origins. Particularly we intend to use these new indices for data mining and visualization, which could offer local and global descriptors of the brain dynamic.

References [1] http://brainvisa.info/.

XI

Fig. 4. MEG Activities in the catch (top) and no-catch (middle) experimental conditions. The corresponding velocity field in the catch condition is shown in green (bottom). The shrinking of activities specific to the catch condition is clearly elucidated by the convergent structure of the flow.

[2] http://neuroimage.usc.edu/brainstorm. [3] G. Allaire. Analyse numérique et optimisation. Editions de l’Ecole Polytechnique, first edition, 2005. [4] G. Aubert, R. Deriche, and P. Kornprobst. Computing optical flow via variational techniques. SIAM Journal of Applied Mathematics, 60(1):156–182, 1999. [5] S. Baillet, J.C. Mosher, and Leahy R.M. Electromagnetic brain mapping. IEEE Signal Processing Magazine, november 2001. [6] J.L. Barron, D.J. Fleet, and S.S. Beauchemin. Performance of optical flow techniques. International Journal of Computer Vision, 12:43–77, 1994. [7] J.L. Barron and A. Liptay. Measuring 3-d plant growth using optical flow. Bioimaging, 5:82–86, 1997. [8] S.S. Beauchemin and J.L. Barron. The computation of optical flow. ACM Computing Surveys, 27(3):433–467, 1995. [9] T. Corpetti, D. Heitz, G. Arroyo, E. Mémin, and A. Santa-Cruz. Fluid experimental flow estimation based on an optical-flow scheme. Experiments in fluids, 40(1):80–97, 2006. [10] M.P. Do Carmo. Riemannian Geometry. Birkhäuser, 1993.

XII [11] O. Druet, E. Hebey, and Frédéric Robert. Blow-up theory for elliptic PDEs in Riemannian geometry, chapter Background Material, pages 1–12. Princeton University Press, Princeton, N.J., 2004. [12] B.K.P. Horn and B.G. Schunck. Determining optical flow. Artificial Intelligence, 17:185–204, 1981. [13] T. Inouye, K. Shinosaki, S. Toi, Y. Matsumoto, and N. Hosaka. Potential flow of alpha- activity in the human electroencephalogram. Neurosci Lett, 187:29–32, 1995. [14] D. Lehmann. Quantitative and Topological EEG and MEG Analysis (Proceedings, Third Hans Berger Symposium, Jena 1996). Druckhaus Maier, Jena and Erlangen, Germany (1997), chapter From EEG waves to brain maps and to microstates of conscious mentation., pages 139–149. Witte, H.H., Zwiener, U., Schack, B. and Doering, A. (eds.), 1997. [15] H. Liu, T. Hong, M. Herman, T. Camus, and R. Chellapa. Accuracy vs. efficiency trade-off in optical flow algorithms. Comput. Vision Image Understand., 72:271– 286, 1998. [16] H.-H. Nagel. On the estimation of optical flow : relations between different approaches and some new results. Artificial Intelligence, 33:299–324, 1987. [17] R.D. Pascual-Marqui, C.M. Michel, and D. Lehmann. Segmentation of brain electrical activity into microstates: Model estimation and validation. IEEE Trans. Biomed. Eng., 42:658–665, 1995. [18] A. Remond and B. Renault. The theory of electrographic objects. Rev Electroencephalogr Neurophysiol Clin., 2(3):241–56, july-september 1972. [19] J.A. Rossmanith, D.S. Bale, and R.J. LeVeque. A wave propagation algorithm for hyperbolic systems on curved manifolds. J. Comput. Phys., 199:631–662, 2004. [20] C. Schnörr. Determining optical flow for irregular domains by minimizing quadratic functionals of a certain class. Int. J. Computer Vision, 6(1):25–38, 1991. [21] P. Senot, M. Zago, F. Lacquaniti, and J. McIntyre. Anticipating the effects of gravity when intercepting moving objects: Differentiating up and down based on nonvisual cues. J Neurophysiol, 94:4471–4480, 2005. [22] W.K. Strik, R. Chiaramonti, G.C. Muscas, M. Paganini, T.J. Mueller, A.J. Fallgatter, A. Versari, and R. Zappoli. Decreased eeg microstate duration and anteriorisation of the brain fields in mild and moderate dementia of the alzheimer type. Psychiatry Res., 75:183–191, 1997. [23] D.C. Van Essen, H.A. Drury, S. Joshi, and M.I. Miller. Functional and structural mapping of human cerebral cortex: Solutions are in the surfaces. Proc Natl. Acad. Sci. USA, 95:788–795, 1998. [24] J. Wackermann, D. Lehmann, C.M. Michel, and W.K. Strik. Adaptive segmentation of spontaneous eeg map series into spatially defined microstates. Int J Psychophysiol., 14(3):269–83, May 1993. [25] J. Weickert and C. Schnörr. A theoretical framework for convex regularizers in pde-based computation of image motion. International Journal of Computer Vision, 45(3):245–264, December 2001.

Acknowledgements : This project is supported by the French Ministry of Research (ACI New Applications of Mathematics). Thanks to Patrice Senot (Laboratoire de Physiologie de la Perception et de l’Action, Collège de France & CNRS, Paris) for providing the data. Thanks to Stéphane Dellacherie (Commissariat à l’Énergie Atomique) for discussions about mathematical formulation.