MathematicS In Action

Numerical experiments presented in this paper were carried out using the PLAFRIM testbed, being developed under the Inria PlaFRIM development action with support from Bordeaux INP, LABRI and IMB and other entities: Conseil Régional d’Aquitaine, Université de Bordeaux and CNRS (and ANR in accordance to the programme d’investissements d’Avenir,


Introduction
During recent decades, voltage sensitive dye imaging (VSDI) has fruitfully contributed in advancing our knowledge of neuronal circuits and single cells functioning in various brain regions, both in vitro and in vivo [1,15,16,27,31,34].To optically image electrical activity, the preparation under study is stained with a suitable dye, which binds the external surface of cell membranes and acts as transducer transforming variations in membrane potential into changes in absorption or emission of light [15,27].
In the majority of cases, the analysis of VSDI data has so far been restricted to the quantification of membrane potential-dependent variations of fluorescence inside a given region of interest (ROI).Although useful and easy to implement, such analysis misses key information of neuronal computation that in principle is readily obtainable from VSDI data: for instance, (i) how do dynamics of neuronal activity change at time scales relevant for neuronal processes?(ii) are there other parameters during spreading of VSDI signals that could provide additional properties of neuronal circuits?
Several VSDI studies began to address these questions [2,9,26,28,30,34] by developing analytical methods for a deeper interpretation of circuits dynamics.These approaches, however, mainly rely on the analysis of integral quantities [2,9,28,34] and in some cases on ad hoc approaches to determine how the VSDI signal is propagated through the neuronal tissue.For example, Takagaki and collaborators recently employed a flow-detection algorithm [30] to study the propagation of neuronal activity recorded using VSDI.This algorithm exploits a maximum correlation principle, classical in continuum mechanics, to determine spatial flow patterns.For simple propagation patterns, this approach can lead to reasonable results.When, however, the signal spreads in complex threads, a statistical approach displays problems delivering valuable information due to the lack of an underlying model.
Principal component analysis (PCA) is another method classically employed in the analysis of distributed systems [20,26].The main feature of this approach is that it may deliver a small dimensional representation of spatially distributed patterns.A substantial drawback is that the notion of propagating signals is blurred by this empirical modal representation.Wave front propagation of signals is a classical example of these issues: an infinite number of PCA components would be needed to reasonably represent this simple phenomenon.Considering that signal spreading in neuronal networks shares specific features of wave front propagation phenomena, PCA approaches display apparent limited analytical benefit to describe their behavior.
Here we propose a model-based non-statistical objective method for quantifying the distributed information obtained with VSDI, using as experimental model the layer-specific neuronal network dynamics within the CA1 region of the hippocampus.The model originates from the resolution of an optimal mass transportation problem [33] aimed at dissecting propagation patterns from successive time-points of neuronal depolarization spreading.To our knowledge, there is no a connection establishing between optimal transportation and the mathematical models describing stimuli propagation in the neural tissue at the microscopic level.Optimal transportation is considered as a phenomenological model.
The fundamental question behind the optimal mass transportation problem is to determine how to map one image onto the next in an optimal way.In this case, the mapping is optimal in the sense of minimum displacement of elementary signal intensity.We use the Wasserstein distance to measure the optimality and we consider the L 2 -Monge-Kantorovich problem for the analysis.
The nature of the optimal transfer problem is inherently non-linear and may not be solved using linear interpolation techniques.In the simple example of a travelling signal captured in two consecutive images, the two images must be interpolated by an intermediate image to arrive at the underlying transport phenomenon.Linear interpolation techniques are inadequate as they result in an averaged combination of the initial and final signal with no travelling pattern.The Wasserstein distance that is calculated by the optimal mass transportation problem rigorously defines a measure of the difference between several snapshots of a propagation pattern.The dynamic vector fields obtained with this approach determine velocity and overall direction (i.e.divergence) of VSDI signal spreading in a specific region, and can be analyzed by conventional point-wise statistical methods.In addition, taking into account the spatial nature of the timeresolved data, we have the possibility to provide a deeper understanding of the distributed neuronal activity.
Compared to other approaches to detect signal propagation patterns in VSDI, the advantage of the optimal transportation formulation is that it is model-based and no additional parameters are tuned to analyze the data.Furthermore, we show the ease of dissecting the activity and distribution of neuronal signals amongst different sub-layers of a given brain region, which, using classical electrophysiological techniques, is both invasive and time-consuming to achieve.
The outline of the present paper is as follows.In Section 2 we describe the materials and methods used for VSDI experiments.In Section 3 we present the mathematical model based on the optimal transportation theory that we used to analyse the VSDI data.In Section 4 we describe the results obtained from the layer-specific analysis of neuronal dynamics using our model in the CA1 region of the mouse hippocampus.We conclude the paper with a discussion section.

Voltage Sensitive Dye Imaging: materials and methods
In this section, we details all the materials and methods used to extract VSDI data.
Slice preparation and staining with Voltage sensitive dye.Experiments were approved by and carried out according to the local ethical committee of the University of Bordeaux (approval number 501350-A) and the French Ministry of Agriculture and Forestry (authorization number 3306369).8 to 11 weeks-old male C57BL/6-N mice (Janvier, France) were kept with ad libitum access to food and water.After isoflurane anesthesia mice were decapitated and 350 µm thick sagittal slices containing dorsal hippocampus were prepared with a vibratome (VT1200S, Leica, Germany).During this procedure the brain was submerged in oxygenated (95%O 2 /5%CO .Subsequently slices were allowed to recover, at room temperature in the same solution, for at least 30 minutes before starting the dye staining procedure.For dye staining, each slice was incubated for 15 minutes in oxygenated ACSF containing Di-4-ANEPPS (Sigma-Aldrich, France) at a concentration of 16.4 µM (20 mM stock in DMSO, final DMSO < 0.1%).The stained slice was allowed to recover in dye-free ACSF, at room temperature, for at least 45 minutes, before recording.
Optical recording method.Slices were placed in a submerged recording chamber (Membrane Chamber, Scientific Systems Design Inc., Canada) under constant flow of oxygenated ACSF (∼ 2 mL/min), at room temperature.To record neuronal signals with VSDI we used an epifluorescence macroscope (Brainvision, Japan) equipped with the MiCAM02 optical imaging system (MiCAM02 HR; Brainvision, Japan) with a spatial resolution of 33.3 × 37.5 µm (horizontal and vertical, respectively) for each pixel.A stereoscopic microscope (Leica, Germany) was used to visually guide the placement of a concentric bipolar stimulating electrode (FHC Inc., USA) into the proximal (with respect to the CA3) region of the stratum radiatum, to stimulate the Schaffers collateral pathway.For each stimulus, a 200 µs long voltage pulse was applied using an isolated stimulator (DS2A, Digitimer Ltd., United Kingdom).The stimulation intensity in experiments with HU210 with respective control was set at 20 V. Every acquisition consisted of 256 frames sampled every 2.2 milliseconds averaged 15 times.Each experimental point was the mean of four acquisitions interleaved of 20 seconds, averaged by using the utility of the imaging analysis software (Brainvision, Japan).

VSDI data and ROI extractions.
For each averaged data we calculated the fractional change in fluorescence (∆F * F −1 ) using the Brainvision imaging analysis software (Brainvision, Japan) as follows: the intensity of emitted fluorescence prior to stimulation in the first 8 frames was averaged and used as reference intensity (F 0 ).The change in fluorescence [∆F (t) = F (t)−F 0 ] was then normalized by the highest background fluorescence value, and absolute values representing changes respect to background fluorescence were used as the optical signal.From each data file representing each experimental point we extracted 16 frames containing the signal of interest, starting from one frame before stimulus, thus covering all the time-lapse of neuronal depolarization.Exclusively for Figure 4.1 (A) (see Section 4), after calculating ∆F * F −1 , we applied a spatial filter of 5×5 pixels and then isolated the CA1 region by zeroing values of fluorescence outside this region of interest (ROI).A depolarization produces a reduction in fluorescence emitted by Di-4-ANEPPS, while a hyperpolarization an increase; for clarity, ∆F * F −1 values representing depolarization (Figures 4.1 (C) and 4.2) were considered positive.ROIs were manually drawn post-hoc using the Brainvision image analysis-acquisition software (Brainvision, Japan), representative ROIs are shown in Figure 4.1 (B).Great care was taken to match the ROI boundaries with anatomical landmarks, however, the low spatial resolution of our VSDI together with the relatively large pixel size, make an exact subdivision of the CA1 region difficult.As a consequence, the ROI named Radt.Distal (radiatum distal) may contain a limited component of the stratum lacunosum-moleculare, and the ROI named Pyr.Layer (Pyramidal Layer) may include very limited parts of stratum oriens and stratum radiatum.To extract the signal of interest inside each ROI, over the 16 frames of neuronal activity, we used the Brainvision software to zero ∆F * F −1 values outside of ROI boundaries.

Pharmacology.
In experiments with HU210 (Tocris; United Kingdom), we compared separate experiments where we applied either HU210 or the equivalent volume of DMSO (i.e."Vehicle" condition; DMSO < 0.1% v/v final concentration).Drugs were bath-applied for 30 minutes after respective control conditions.

The mathematical model
The optimal transportation problem is increasingly used to model several problems in mechanics, physics, image analysis and other fields, see e.g.[33] and references therein.Because of all these applications, this old topic first introduced by Monge in 1781 has attracted considerable attention these last years.
The model formulation is detailed as follows: let us consider two images of a spreading depolarization signal taken at ∆t one from the other.We consider non-negative signal values defined by ρ 0 (ξ) for the initial frame and ρ 1 (x) for the successive one, where ξ and x are the coordinates of the pixels within each frame.We model the signal spreading from one image to the other by a phenomenological model requiring that the overall squared distance travelled, weighted by signal intensity, is the smallest possible.Let X : Ω 0 → Ω 1 a smooth one-to-one map that takes every point from the first image Ω 0 to the next Ω 1 .The model is then defined introducing the L 2 squared Wasserstein distance to be minimized: inf Let m 0 and m 1 the total intensity of the first and second frame, respectively.The minimization in (3.1) is carried out under the constraint that X(ξ) realizes the transfer of ρ 0 /m 0 onto ρ 1 /m 1 : where c = ln(m 1 /m 0 ) and With respect to the classical optimal transportation problem, we consider here an unbalanced transport to take into the difference of intensity between ρ 0 and ρ 1 .From the biological point of view this means that the firing neuron population is increasing proportionally to its size.Other authors have considered extensions of the optimal transport problem that allow for variation of mass [3,11].The phenomenological model that we propose, hence, relies on the choice of the Wasserstein metric to define the cost of the mapping and on the model of signal growth implicit in the Jacobian equation (3.2), which is an intensity balance equation.This problem is called the exponential growth problem (EGP).The EGP can be reformulated introducing a temporal dimension so that, given Π : where the minimum is taken among all densities ρ(t, x) ≥ 0 and velocity fields v(t, x) ∈ R 2 satisfying the continuity equation: with the initial and final conditions As in [4] we compute the formal optimality conditions of the space-time minimization problem (3.4) satisfying (3.6) and (3.5): where ϕ is the Lagrange multiplier of the constraint (3.6) and (3.5).
Finally, the constrained Euler-Lagrange equations are: and the problem consists in determining the initial velocity u(0, • ) = u 0 .The solution to this problem exists and it is unique as shown in the following.Let ρ1 = ρ 1 (X(ξ))/e c , then the minimization of (3.1) subject to: becomes the classical Monge-Kantorovich problem (MKP), which admits a unique solution [6,32] and the minimizing mapping X * (ξ) can be written as the gradient of some convex potential.
Let us then denote the solution of the MKP by ∇Ψ = X * (ξ) − ξ.From equation (3.5) we have that: and hence The frames relative to the propagation of VSDI-recorded depolarization are such that c is very close to 0 since m 0 and m 1 do not considerably differ.In this limit, the initial velocity u 0 (ξ) ≈ ∇Ψ that is the solution of the equivalent MKP.Finally, this velocity is multiplied by 35.4 µm (the averaged pixel size for both x and y), divided by the actual experimental ∆t (2.2 milliseconds) and converted to meters per second, which is the final unit reported here.
As a result, given a time we can compute Wasserstein distance between two subsequent frames and the initial velocity field of the mapping.In addition, as a mean of interpreting the results, the rate of volume change due to the mapping is computed by integrating the divergence of the velocity field: In this paper, we use a robust and fast numerical scheme based on a Newton approach, proposed in [5] to approximate the solution of the L 2 -Monge-Kantorovich problem.Another efficient algorithm has recently been proposed in [19].We consider in particular the following iterative method that can readily be coded as it implies the solution of a rather simple elliptic PDE, see Algorithm 1.

Validation of the algorithm based on the optimal transportation problem
In order to validate the reliability of the algorithm based on the optimal transportation problem, we analyzed three different simulated datasets: (i) displacement by predetermined distances of a simple Gaussian distribution and (ii) of a VSDI-recorded depolarization, and (iii) the mapping of artificial data with a predefined velocity field.In these conditions, the algorithm should provide vector fields representing Wasserstein distances that coherently reflect the displacements and the mapping of the simulated signals.
7. go to 3 if convergence is not attained; In the first test (Figure 3.1 (A)), we assessed how well the algorithm could describe the translation of a two-dimensional Gaussian in a Cartesian plane.Here a two dimensional Gaussian (ρ 0 ) with x, y coordinates 60, 30 was shifted to give ρ 1 at x, y coordinates 25, 42.From a theoretical point of view, the exact solution of the optimal mass transportation problem coincides with the translation of ρ 0 onto ρ 1 .Indeed, as shown in Figure 3.1 (A), when we computed the Wasserstein distance between the two centers (vector, superimposed blue arrow) we found a value of 37 pixels, which coherently reflects the operated translation.To further test the accuracy of the algorithm, we next analyzed the predetermined motion of a more complex shape, such as a real VSDI-recorded depolarization within the whole CA1 hippocampal region.As shown in Figure 3.1 (B), shifting each value of fluorescence intensity (ρ 0 , left panel) by 1 pixel along the y axis (x, y + 1 coordinates) leads to ρ 1 (right panel) with a corresponding vector field (lower panel) in which every vector far from the border is 1 pixel in length and points in the upward direction, correctly reporting the displacement.We next present a mapping of a predetermined radial initial velocity.We defined a velocity field a priori as the gradient of a given potential (i.e. a solution of the optimal transportation problem, "defined velocity"), which corresponds to a source of given intensity ρ 0 (Figure 3.1 (C), left panel).Using this field, we mapped ρ 0 onto ρ 1 (not shown) and then, with our algorithm, used these two densities to compute the mapping between the two (i.e."solution velocity", Figure 3.1 (C), right panel) a posteriori.As shown in Figure 3.1 (C) (right panel), if we compare the known velocity field defined a priori with the one obtained a posteriori, by optimal displacement, we can appreciate complete overlapping between the two.
As a final validation of our method, below we numerically demonstrate how we obtain corresponding values by calculating the integral of divergence for the velocity fields mentioned above.Let us define a circle domain Ω of center (x c , y c ), radius r, and a velocity field as follows: where η is a real positive number that assures the regularity of v and its integrability over the domain Ω.In this example, we consider η = 0.001, r = 5 and c = 1.We can compute the integral as (3.10): By calculating the integral divergence of the vectorial field obtained a posteriori, we obtain the same number: 123 Altogether, these data demonstrate the reliability of our method in coherently reporting changes in travelled distances, directions and divergence of vectorial fields resulting from the computation of optimal transportation theory.(A) A Gaussian representing a density function ρ 0 with centre coordinates (60, 30) pixels on x and y axis respectively has been shifted 35 pixels to the left and 12 pixels upward to obtain ρ 1 , with x, y coordinates 25, 42.Superimposed blue arrow is the vector (37 pixels in length) reporting the displacement according to our algorithm.(B) Upper panel, a density function ρ 0 representing VSDI-recorded depolarization in the whole CA1 area has been shifted upward by 1 pixel to obtain ρ 1 (heat maps were constructed using a 7 colour lookup table, red denotes the highest depolarization values).Lower panel, the velocity field reports the displacement of each depolarization value in each pixel according to the algorithm.Note that every pixel far from the border has a length of 1 pixel and points upward.(C) The left panel shows a predefined velocity field (superimposed blue arrows) and the associated density function ρ 0 (heat maps constructed using a 7 colour lookup table, red denotes the highest density values).The right panel shows the velocity field obtained after calculation of optimal displacement (red arrows) of ρ 0 to ρ 1 (not shown).Green dashed line represents a circle domain Ω of center (x c , y c ).

Relationship to optical flow
Other approaches can be proposed to compute distributed vector fields from successive images.For instance, a well-known way to determine vector fields is the optical flow method, first introduced by the psychologist James J. Gibson in the 1940s.As an example of this strategy for the analysis of VSDI data, Takagaki and collaborators used a flow-detection algorithm that determines spatial flow between two detectors by searching for the time shift which gives the optimal correlation coefficient between the signal from each detector, within a given correlation window [30].The initial problem formulation in optical flow estimation is based on a simple constant brightness constraint: where ρ is the intensity at the point x at the instant t, dx the displacement in space and dt the time lag.Under the assumption of "small displacements" this constraint can be rewritten as: where v = dx/dt is a velocity field.Coherent with this assumption, the time derivative in this equation is approximated by the difference between the image intensities and the space derivative, by finite differences on the average image.Yet, this equation is underdetermined since there are two unknowns, which are the two components of the velocity field, v. Thus, to estimate the flow, several approaches with additional constraints have been proposed.Among them, the well-known Horn-Schunck method [17] combines density advection with a regularization on the velocity fields, and is based on minimization of the following function: where D is the relevant image domain in R 2 and λ is a regularization parameter (the larger λ is, the smoother the flow).This method depends on: (i) the time lag between consecutive images and (ii) the parameter λ that has to be calibrated.Alternative variations on the optical flow calculation follow comparable concepts (e.g. the Lucas-Kanade method [22]) and they are based on similar regularization terms requiring additional information to solve the problem.In contrast to the Horn-Schunck optical flow method described above, the optimal transportation theory is independent of the tuning of any parameters.More specifically, the velocity field from the solution of the optimal transportation problem is explicitly obtained by direct integration and the result is exact, whereas the Horn-Schunck method provides a velocity field that depends on the parameter λ and the extent of the displacement.
To demonstrate this difference, we calculated the velocity of a rigid translation of a synthetic Gaussian in one space dimension.We considered a Gaussian ρ 0 that is translated by a given ∆x in space ρ 1 over a reference time lag of 1.As shown in Figure 3.2 (A) (middle panel), by considering λ = 0.5 and ∆x = 0.5, the velocity field obtained with the Horn-Schunck optical flow method is reasonably correct, returning a value v = 0.51, within tolerable error of the ground truth v = 0.5 (constant in the area where the Gaussians are centered).Note that with the optimal mass transportation algorithm, we obtain exactly v = 0.5 and there were no parameters to tune.Next we plotted the velocity field for the same translation shown in Figure 3.2 (A) (left panel), but this time the velocity field was determined, using Horn-Schunck method, with λ = 1 (Figure 3.2 (A), right panel).In this scenario we obtained an incorrect velocity of 0.36, thereby demonstrating how the parameter λ must be carefully selected based on previous experience or calibration.In the final simulation (Figure 3.2 (B)), we show that for a relatively large spatial shift (∆x = 1.5) and for λ = 0.5, the numerical velocity obtained with the Horn-Schunck method is variable in space and completely inconsistent compared to the exact value of 1.5, demonstrating the inaccuracy of this method in modelling the optical flow for large image displacements.Yet, solving the optimal transportation problem once again yielded the exact velocity, 1.5.Altogether, these simple simulations show certain limitations of currently used optical flow algorithms, revealing that: (i) the method is not accurate when the displacement between two successive images is too large, and (ii) even in the case of small displacements, λ needs tuning to obtain an accurate result.The choice of λ is a critical step because if it is too small the algorithm does not converge, and too large the approximation is inaccurate.The main advantage of the optimal mass transportation model considered here is that it is self-contained, requiring neither an intrinsic hypothesis on the time lag between two images, nor arbitrary parameters to compute the velocity fields.Compared, however, to the optical flow methods, the mass transportation model is computationally more demanding, yet easily handleable by modern computers.

Results
The 16 frames containing the signal extracted from each ROI were used to calculate the optimal displacement of neuronal depolarization between successive frames, according to the optimal mass transfer problem.All the numerical simulations have been performed using a robust and fast numerical scheme proposed in [5].The script used for the estimation of the optimal displacement was realized using Matlab (MathWorksš, R2016b).

Impact of manipulating excitation and neuromodulation on networks dynamics
within the CA1 strata.After quantitatively analyzing the dynamics of evoked depolarizations, propagating within the hippocampal CA1 circuit under basal conditions, we next studied how manipulation of the circuit can affect the routing of neuronal depolarization while it travels across the CA1 subfields.We analyzed two different conditions: (i) increasing Schaffers collaterals stimulation intensity from 10 to 30 Volts and (ii) application of the synthetic agonist of cannabinoid type 1 receptor (CB1) HU210 at an intermediate stimulation intensity (20 Volts).Each of these manipulations should have a mechanistically different impact on the flow of depolarization through the CA1 network, since increased stimulation intensity is able to recruit more axonal fibers, whereas activation of CB1 receptors is very well-known to regulate neurotransmitters release in the hippocampus [21].We first considered how the two manipulations impact the magnitude of neuronal depolarization.Increasing the stimulation of Schaffers collaterals from 10 to 30 Volts significantly amplified the magnitude of neuronal depolarization in the whole CA1 network throughout the duration of propagation (Figure 4.2 (A)).This effect was present in all the CA1 sub-layers (Figure 4.2 (B)-(E)).Conversely, application of HU210 (1 µM) did not change the overall intensity of neuronal depolarization in the CA1 region as a whole (Figure 4.2 (F)), but it slightly but significantly decreased the signal in the pyramidal layer and particularly during the early/rising phase of signal spreading (Figure 4.2 (H); see time course and inset bar chart).Next, we considered how the two manipulations influence the velocity of spreading of neuronal depolarization.Signals elicited by 30 Volts stimulation reached peak velocity 2.2 ms earlier than those elicited by 10 Volts, and then decayed more slowly.In addition, stimulation at 30 Volts was associated with a significantly lower signal velocity during the middle-end of propagation (Figure 4.3 (A)-(E)).While the impact of stronger stimulation on velocity were similar throughout the CA1 sub-fields,the activation of CB1 receptors by HU210 decreased velocity of the VSDI signal specifically in the stratum oriens and in the pyramidal layer, an influence confined to the early phase of spreading (Figure 4.3 (G) and (H); see time course).Divergence, in the whole CA1, was significantly lower throughout the response when Schaffer's collaterals were stimulated at 30 Volts as compared to 10 Volts (Figure 4.4 (A)).This was true for all phases of the response in all CA1 sub-fields (Figure 4.4 (B)-(E)), with the one exception of the early phase response

Discussion
The spatially-distributed nature of VSDI data requires the development of new methods able to provide valuable information allowing to overcome the mere quantification of local changes in membrane potential-dependent fluorescence emission at time-fixed intervals.These results were obtained using a new approach in the analysis of VSDI data, based on the application of the optimal mass transportation problem centered on Monge-Kantorovich theory, that is a deterministic plan to transfer a quantity of mass from a starting configuration to a final one by minimizing a functional cost [33].In our case, the mass to be transported is the VSDI-recorded depolarization and the functional cost to be minimized is the distance.The advantages of the present method to describe spreading of VSDI signals in brain networks as compared to already existing analytical tools [23,26,28,30,34] rely on the fact that optimal mass transportation models do not require tuning of parameters for the analysis and use unbiased underlying deterministic parameters.With this analysis, we obtained vectorial fields in which each vector represents the shortest distance travelled by neuronal depolarization in each pixel during its propagation in the CA1 hippocampal neuronal networks at milliseconds time resolution.In this study, the combination of VSDI and mathematical analysis based on the optimal transportation theory allowed us to characterize the stimulated spreading of depolarization in basal condition and to observe effects that discriminate between distinct manipulations of neuronal networks, within anatomically defined sub-regions of the hippocampal CA1 at defined time intervals.Analysis of VSDI data from the stimulation of Schaffers collaterals in the CA1 region in basal condition underlines how the signal moves over time at a rather non-uniform speed and tends to converge during its propagation, as divergence of the underlying optimal transportation mapping decreases over time.This is particularly evident from a manipulation that busts network activity such as increased stimulation intensity.This manipulation leads to a slower and less divergent propagation of depolarization signals throughout the CA1, during the whole propagation period.To explain this somehow counter intuitive finding, we may speculate that, at equal time intervals, the increased dendritic spatial summation, provided by the stronger stimulus, leads to a more persistent depolarization along the dendritic field.This is associated with a decreased signal divergence because of the sustained converging dendritic inputs.Indeed, voltage sensitive dye imaging reflects mainly dendritic postsynaptic activity [10].Another important and somehow surprising finding of the present study is that while changes in velocity and divergence induced by increased stimulation intensity are conserved in all CA1 sub-regions, the much less dramatic effects of the activation of CB1 receptors are layer-specific, as significant decreases in fluorescence emission and signal peak velocity appear only in stratum oriens and in pyramidal layer.These phenomena may be explained by the well-defined expression pattern and signaling efficacy of CB1 in the CA1 region of hippocampus.Indeed, CB1 receptors are expressed both in excitatory pyramidal neurons and in inhibitory GABAergic interneurons of the hippocampus [24].Through different mechanisms including both presynaptic and postsynaptic actions, the activation of CB1 receptors generally results in a decrease of neuronal excitability [8].Therefore, application of a CB1 receptor agonist might in principle induce both inhibition of excitation and of inhibition [8], thereby generating opposite effects on the control of synaptic transmission.Interestingly, whereas CB1 receptors are expressed at much higher levels in hippocampal GABAergic interneurons than in glutamatergic pyramidal neurons, their signaling in the latter cell type appears to be much more efficient [29], possibly resulting in a net slight and localized inhibition of neuronal signal spreading.Future studies will determine whether anatomical and functional cell type-specific features of the receptor underline the effect of CB1 agonism, but In the context of the present study, however, it is important to notice that the model we propose does not record only major changes in the   spreading of neuronal signals (such as the one obtained by increasing stimulation intensity), but it is sensitive enough to capture even slight but significant variations in these dynamics.By combining the advantages of VSDI in recording neuronal activity from different region of interests with millisecond time-resolution and the mathematical analysis of the VSDI data, we provide a novel and innovative method for the analysis of VSDI data, which might potentially be used for the extraction of further dynamic information from all imaging techniques.Moreover, similar methods can be applied to a large range of scales, from single cells to neuronal networks in different brain regions, both in slices and in vivo.These data will provide a useful novel tool to further dissect complex phenomena such as spreading of neuronal signals in large brain regions.

Figure 3 . 1 .
Figure 3.1.Validation of the algorithm with artificial data.(A)A Gaussian representing a density function ρ 0 with centre coordinates (60, 30) pixels on x and y axis respectively has been shifted 35 pixels to the left and 12 pixels upward to obtain ρ 1 , with x, y coordinates 25, 42.Superimposed blue arrow is the vector (37 pixels in length) reporting the displacement according to our algorithm.(B) Upper panel, a density function ρ 0 representing VSDI-recorded depolarization in the whole CA1 area has been shifted upward by 1 pixel to obtain ρ 1 (heat maps were constructed using a 7 colour lookup table, red denotes the highest depolarization values).Lower panel, the velocity field reports the displacement of each depolarization value in each pixel according to the algorithm.Note that every pixel far from the border has a length of 1 pixel and points upward.(C) The left panel shows a predefined velocity field (superimposed blue arrows) and the associated density function ρ 0 (heat maps constructed using a 7 colour lookup table, red denotes the highest density values).The right panel shows the velocity field obtained after calculation of optimal displacement (red arrows) of ρ 0 to ρ 1 (not shown).Green dashed line represents a circle domain Ω of center (x c , y c ).

Figure 3 . 2 .
Figure 3.2.Comparison between the optimal mass transportation problem and the Horn-Schunck optical flow method.(A) Left panel, a Gaussian function ρ 0 centered at x = 0 is shifted to x = 0.5(ρ 1 ).Middle panel, velocity field obtained with the Horn-Schunck optical flow method from the translation of the Gaussian function in left panel with λ = 0.5.Right panel, velocity field obtained with the Horn-Schunck optical flow method from the translation of the Gaussian function in left panel with λ = 1.(B) Left panel, a Gaussian function ρ 0 centered at x = 0 is shifted to x = 1.5(ρ 1 ).Right panel, velocity field obtained with the Horn-Schunck optical flow method from the translation of the Gaussian function in left panel with λ = 0.5