Searching for Conservation Laws in Brain Dynamics—bold Flux and Source Imaging

Blood-oxygen-level-dependent (BOLD) imaging is the most important noninvasive tool to map human brain function. It relies on local blood-flow changes controlled by neurovascular coupling effects, usually in response to some cognitive or perceptual task. In this contribution we ask if the spatiotemporal dynamics of the BOLD signal can be modeled by a conservation law. In analogy to the description of physical laws, which often can be derived from some underlying conservation law, identification of conservation laws in the brain could lead to new models for the functional organization of the brain. Our model is independent of the nature of the conservation law, but we discuss possible hints and motivations for conservation laws. For example, globally limited blood supply and local competition between brain regions for blood might restrict the large scale BOLD signal in certain ways that could be observable. One proposed selective pressure for the evolution of such conservation laws is the closed volume of the skull limiting the expansion of brain tissue by increases in blood volume. These ideas are demonstrated on a mental motor imagery fMRI experiment, in which functional brain activation was mapped in a group of volunteers imagining themselves swimming. In order to search for local conservation laws during this complex cognitive process, we derived maps of quantities resulting from spatial interaction of the BOLD amplitudes. Specifically, we mapped fluxes and sources of the BOLD signal, terms that would appear in a description by a continuity equation. Whereas we cannot present final answers with the particular analysis of this particular experiment, some results seem to be non-trivial. For example, we found that 3690 during task the group BOLD flux covered more widespread regions than identified by conventional BOLD mapping and was always increasing during task. It is our hope that these results motivate more work towards the search for conservation laws in neuroimaging experiments or at least towards imaging procedures based on spatial interactions of signals. The payoff could be new models for the dynamics of the healthy brain or more sensitive clinical imaging approaches, respectively.


Introduction
Understanding functional mechanisms underlying brain dynamics is essential for clinical applications such as recovery from stroke, traumatic brain injury, and disorders of consciousness.Neuronal activity in response to perceptual stimuli or cognitive tasks can be mapped by using the blood-oxygen-level-dependent (BOLD) response as a proxy signal in functional MRI experiments.Functional MRI (fMRI) studies have provided us with detailed insights into local functional organization of the brain.Recently, fMRI also has allowed for mapping of brain networks, mostly during unconstrained tasks or "rest".Whereas physiological models exist that describe the local relationship of neuronal processing, hemodynamics, and the observed signal [1][2][3], which can be numerically validated [4,5], the physiology underlying the spatial relationships among signals is much less understood [6][7][8][9].Nevertheless, many recent studies relate BOLD resting state networks to functional networks of the brain, for example by utilizing information theoretical concepts or other notions from physics such as small world networks [10][11][12][13][14].
Despite these efforts to link BOLD signals of the brain to physical concepts, it seems that one fundamental way of physical modeling has not been applied to functional brain imaging yet: The concept of conserved quantities and conservation laws.Most laws of modern physics can be derived from conservation laws, which again follow from fundamental properties of the world such as space and time homogeneity.As the brain becomes more and more understood as a neurovascular system, and as functional imaging technologies can depict brain function with reasonable high spatial resolution, in this contribution we ask: Does the brain described on spatiotemporal scales accessible with fMRI experiments reveal conservation laws of observable quantities that might assist us in modeling brain dynamics?Although we will not be able to provide final answers within this manuscript, we hope that the proposed research methodology will inspire more applications in the quest for a better understanding of spatiotemporal brain dynamics.
As a first step towards this goal, we performed an fMRI experiment on a group of nine healthy subjects, in which the response to a complex mental motor imagery task [15] was mapped.We then quantified the spatial relationship between local BOLD signals in order to define BOLD fluxes, sources and sinks.Conservation laws then might show up as a balance between sources and sinks, for example.
The manuscript is organized as follows: We briefly review the importance of conservation laws and mention possible hints for conservation laws affecting brain dynamics (Section 2).We then describe a model based on continuity equations and associated quantities such as fluxes and sources (Section 3).These quantities are then measured in an fMRI imaging experiment (Section 4).Although this experiment does not provide definite answers yet, it provides hints for a more global organization of brain function not seen before in conventional BOLD mapping (Section 5).A discussion concludes the paper (Sections 6).

The Importance of Conservation Laws
We are aiming at a description of the spatiotemporal dynamics of the BOLD signal by physical principles.The reason is that physical models of complex systems are often generalized models encompassing a whole class of specific systems.For example, many different pattern-forming systems can be described by generalized equations independent of the specifics of the system's components [16].Similarly, conserved quantities such as energy, mass, or momentum are applicable for a wide range of systems.In this sense, we are not aiming at describing the spatiotemporal dynamics of the brain by specific physiological parameters (although this is of course a desirable goal, too), but on a phenomenological level by general physical laws.
In particular, conservation laws state that certain system properties do not change as the system evolves, or, if they change, how they are balanced by sources and sinks.For example, the total energy is conserved for many systems during time evolution, but not necessarily the entropy [17,18].In particular for systems out of equilibrium there are models for the entropy that take its change into account by assigning it to source terms [19].
Here we are interested in conservation laws, including non-conserved quantities described by sources and sinks, for two reasons: (i) There exists a quite general way of expressing conservation laws in spatiotemporal systems without the need to specify the specific nature of the quantity: Continuity equations explicitly relate the local change of a physical quantity to nonlocal influences, namely flow or flux terms, and to source and sink terms.(ii) On a certain level of description, the brain can be considered as a flow-equilibrium system, which naturally gives rise to a possible description by continuity equations.

Possible Hints for Conservation Laws Ruling the Brain
It is well known that the brain consumes, relative to its contribution to body mass, an over-proportional amount of energy [20].During normal brain operation, blood continuously transports the main source of energy, glucose and oxygen, into the brain, and waste is exported in the same way.On average, the inflow of blood into the brain must equal the outflow.This already gives rise to a mass conservation law that can principally be expressed by continuity equations.Similar mass conservation laws will hold for smaller vessels down to the capillary level [21,22].In addition, cerebral autoregulation ensures that average energy supply is kept approximately in equilibrium.Unlike a muscle that consumes less energy during rest than during work, the brain seems to consume a considerable amount of energy even in its "resting" condition.Fluctuations from this background activity, which can be picked up with fMRI, for example, are thought of consuming only a small percentage of this energy [23].Therefore, seen as a physical system, even the working brain can then be considered on average to be in a flow-equilibrium.In addition, tight energy constraints on the brain created a highly optimized organ with probably little energy to waste [24].It stands to reason that flow-equilibrium together with optimized energy usage might give rise to conserved quantities even trickling down to local levels of description.
The simplest sign of conserved quantities observable in BOLD signals would be a negative response that on average offsets the positive responses.However, in many experimental tasks including the one used in this paper, negative responses are not seen consistently and therefore do not persist on a group level analysis over subjects.Local negative responses do exist however and have been consistently found in occipital cortex in humans [25], adjacent to regions with positive response.Negative responses have been related to decreased neuronal activation [26][27][28] but also alternatives such as reallocation of resources and hemodynamic effects [29][30][31] have been proposed.Further, negative responses have been observed in high-resolution imaging of cortical layers [32], in pathological conditions [33][34][35], to noxious stimuli [36], and, again adjacent to positive responses, in indirect response to optogenetic stimulation [5,37,38].Notwithstanding the still debated origins of negative BOLD [39,40], these and other observations of negative BOLD show that there might be mechanisms of neuronal, vascular, or neurovascular organization that could involve yet to be determined conserved quantities.
There is another hint towards an organization of global brain dynamics that might give rise to conservation laws: During rest, BOLD amplitude baseline fluctuations exhibit a well-defined spatiotemporal dynamics [6,7].The most prominent resting state network, the "default mode" was originally discovered in resting PET data [41] but also reveals itself in resting state as well as task-dependent BOLD [8,9].There are task-positive as well as task-negative responses of areas belonging to the default mode network during attention demanding tasks or goal directed behavior [41,42].A biological selective pressure for the evolution of conservation laws controlling the distribution and regulation of changes in blood flow in the brain is the closed volume of the skull somewhat limiting the expansion of brain tissue by increases in blood volume [43].An active law-like reallocation of blood volume based on the regularities of correlations of components of the specific neuronal networks identified in studies of the resting state [44][45][46] may help keep intracranial blood volume relatively constant over the wide range of blood pressures seen in the flat part of the cerebral autoregulation curve.
As another argument for the possible importance of conservation laws in brain dynamics we cite very recent work that found conservation laws ruling at completely different time scales: biological brain evolution, brain growth, and ageing, demonstrated on cortical minicolumns [47].These findings concern structural changes of the brain, but brain growth and ageing also can be seen as very slow dynamical changes.
Finally, from cognitive-behavioral studies it becomes more and more appreciated that attention is directed at only one (or at most two [48]) goals at a time and that "multitasking" [49] of the brain probably is really a shift of attention between multiple cognitively demanding tasks.This phenomenon shows that some parts of the brain, even if thought to be attributed to different tasks, probably cannot function independently at the same time.One might be tempted to speculate that the phenomenon of attention focusing follows from constraints on the total computational power of the neuronal network of the brain, or in other words, a conservation law.

Continuity Equations
We are aiming at an interpretation of data in terms of continuity equations, which are tightly connected to the concept of conserved quantities.The continuity equations used in fluid dynamics, electromagnetism, and many other physical disciplines relate the time change of a quantity to fluxes and sources of that quantity.Continuity equations are either differential or integral equations, which are related to each other by the Gauss theorem.Here we use the differential notation because differential quantitates better suit the local estimation approach that we are adopting for the analysis of fMRI data.However, integral equations might be advantageous for further analysis procedures not described here, for example to measure the fluxes into or out of particular brain regions.
Our main assumption is that the BOLD amplitude represents a density or concentration of something that can drive a gradient flux, and that we can neglect advective flux on this level of description.If the observed BOLD amplitude in an imaging voxel located at position r = (x, y, z) and time t is denoted by ρ(r,t), the corresponding gradient flux j(r,t) = (j x (r,t), j y (r,t), j z (r,t)) can be written as: The constant D relates the magnitude of the flux to the concentration gradient grad ρ(r,t) = (∂ρ(r,t)/∂x, ∂ρ(r,t)/∂y, ∂ρ(r,t)/∂z)), and is assumed constant in space and time.In the neurovascular coupling balloon model of Buxton and others [2,3,[50][51][52], the BOLD signal is a nonlinear function of blood deoxyhemoglobin content and blood volume.As a working hypothesis we therefore assume that blood volume or hemoglobin concentrations can influence neighboring regions in the way as described by Equation ( 1), either directly or mediated by neuronal activity.Without specifying the underlying physiological mechanisms, here we just use our working hypothesis in order to test for a continuity equation of the BOLD signal that might have to be modified in further searches for conservation laws.The continuity equation: where div j(r,t) = ∂j x (r,t)/∂x + ∂j y (r,t)/∂y + ∂j z (r,t)/∂z is the divergence of the flux j(r,t) and σ(r,t) a source term, can now be written as: with the Laplace operator In an fMRI block design experiment, one typically has two alternating states, periods of task execution followed by control periods of rest.Denoting the task state with "1" and the control state with "0", from Equation (3) follows: ∂ρ 0 (r,t)/∂t = D Δ ρ 0 (r,t) + σ 0 (r,t) and ∂ρ 1 (r,t)/∂t = D Δ ρ 1 (r,t) + σ 1 (r,t) ( Finally, the BOLD response is approximated as a block function proportional to a characteristic function consisting of "0"s and "1"s indicating periods of rest and task, respectively: During rest, ρ 0 (r,t) = ρ 0 (r), during task, ρ 1 (r,t) = ρ 1 (r).In other words, the BOLD amplitude is modeled as constant over time (but not space) during task and control periods and its explicit time derivatives vanish.Therefore, we end up with expressions for the BOLD sources tailored to a simple block design fMRI experiment: In other words, the Laplacians of the BOLD signal for task and control blocks provide estimates for the sources, up to a constant.It is this variable then that is regressed against the characteristic function in order to estimate an effect, instead of the BOLD amplitude itself.

Estimation of Fluxes and Sources from BOLD Data
To estimate sources of the BOLD signal, as per Equation ( 5) the Laplacians of the signals need to be computed for each voxel and time point.The Laplacian can be approximated by a finite difference scheme that involves nearest neighbor voxels.Here we adopt a scheme that is accurate to second order in the time interval [53] in order to avoid having to use voxels further away from the voxel of interest, which could blur the estimate.Estimation of the fluxes j(r,t), up to the constant D, works equivalently by computing the numerical gradient for second order accuracy (Equation ( 1)).However, the three components of j(r,t) each on their own do not contain much information.Therefore, in the application we test for the hypothesis that the Euclidean norms of j 0 (r,t) and j 1 (r,t) do not differ.This hypothesis does not account for changes of direction, unfortunately, and there is still room for improvement.The details of the statistical procedures as applied to experimental data are provided in Section 4.

Generalizations
If the explicit time dependence of the BOLD amplitude ρ(r,t) has to be taken into account, for example in the case of slow changes that cannot be described by a block function, or when there are no defined alternating task and rest periods such as in resting state fMRI, Equation ( 5) does not hold.In this case, ∂ρ(r,t)/∂t needs to be estimated from data, too, and used as an additional regressor.For example, in a two-variable regression model with the variables ∂ρ(r,t)/∂t and Δ ρ(r,t), where D is then the fit coefficient (see Equation ( 3)), the sources σ(r,t) would have to be defined differently.It depends on the specific situation if these are valid assumption and is not topic of this paper.

Subjects
Experiments were performed on nine healthy adult subjects (37 ± 8 years, one female) and one 23-year-old man who had had a severe traumatic brain injury about 4 years before scanning.All experiments were approved by the Institutional Review Board of Weill Cornell Medical College, and informed consent was obtained from the subjects or their proxies.

MRI
Data were acquired on a 3.0 T General Electric MRI scanner (healthy subjects) and on a 3.0 T Siemens TIM Trio MRI scanner (patient).For fMRI, a GE-EPI sequence was used (TR = 2 s, TE = 30 ms, flip angle = 70°, voxel size = 3.75 × 3.75 × 5 mm on the GE and 3.75 × 3.75 × 4 mm on the Siemens scanner).Before the scan, the healthy subjects were instructed to imagine themselves swimming, starting and stopping with pre-recorded cues received over MRI compatible headphones.The patient was instructed to imagine himself swinging a tennis racket with his right hand, and cues were delivered over in-room loudspeakers.In the interim, all subjects were required to think of nothing in particular.Eight blocks of 16 s rest alternating with 16 s task were used.Data acquisition consisted of 128 scans with a total time of 4:16 min.In addition, anatomical high-resolution scans of the whole head were acquired.

fMRI Data Analysis
The following preprocessing procedures were applied using SPM5 [10]: motion correction, slice timing correction, and spatial normalization to an MNI template [icbm_avg_152_t1_tal_lin.mnc, MNI, Canada (ICBM, NIH P-20 project)], including interpolation to a voxel size of 2 × 2 × 2 mm.Then, using Matlab, the data was normalized by subtracting the mean over all voxels and time points and then dividing by the standard deviation over all voxels and time points.Note that we did not normalize each voxel time series by its own mean and standard deviation as this would render the spatial derivatives, computed later on, useless after averaging.Each voxel time series was then high-pass filtered with a cutoff frequency of 0.01 Hz.A characteristic function consisting of "0"s and "1"s indicating blocks of eight samples was shifted two samples against the functional paradigm to account for hemodynamic delay.We did not use a more detailed hemodynamic model in order to keep the analysis comparable between conventional BOLD amplitude mapping and our approach of spatial interaction of amplitudes mapping.
(1) Conventional BOLD amplitude mapping: The linear correlation coefficient between the characteristic function and each preprocessed voxel time series was converted into z-values by Fisher's z-transform [67].This linear least-squares procedure models the case that the relationship between data and characteristic function is linear with unexplained variance expressed through normally distributed residuals.
(2) Spatial interaction of BOLD amplitude mapping: For each voxel, the linear correlation coefficient between the characteristic function and the flux norm was computed and converted to z-values.A positive z-value means that during the task the absolute value of the flux as per Equation (1) increases.The procedure was repeated for the source estimate as per Equation (5).Again, any unexplained variance is modeled by normally distributed residuals.
For the group analysis, the z-maps were spatially smoothed with SPM5 with an 8 mm full-widthhalf-mean Gaussian kernel [68].Subsequently, a mixed effects group analysis of all nine scans was performed on those z-maps, in which the z-values were used as the fit coefficients.This quite conservative approach includes inter-subject variance and allows for inference on the population level.The null hypothesis tested was that this group does not activate on average, using a one-sample t-test within a summary statistics approach, as typical in a mixed (or random) effects analysis in SPM [69].The occurrence of false positives in the group SPMs of the flux was controlled by using a false-discovery rate (FDR) of q = 0.04, based on the data masked by the MNI template to exclude most extracerebral voxels, and computed by xjView (http://www.alivelearn.net/xjview); in the flux map, q = 0.04 corresponds to an uncorrected p < 0.005, which was used in all other maps.Also, clusters smaller than 27 voxels were discarded.Fixing the threshold to p < 0.005 (t > 3.36) was necessary for two reasons: (1). to make different maps comparable as the FDR critically depends on the volume of activation.
(2).It was not possible for example to fix q = 0.05 as otherwise the clear and strong BOLD amplitude map observed, concordant with literature results in motor imagery [15,[54][55][56][57]59,61,63], would not have shown up in this conservative analysis, although visible in all individual subjects.Maximal intensity projections or "glass brains" and SPMs were finally displayed by xjView, too.
For the case study subject in the discussion section, no smoothing was applied and the threshold for p < 0.005 was set to |t| > 2.81 under the assumption of a z-distribution due to the larger number of degrees of freedom now available.

BOLD Amplitude Mapping
A conventional BOLD amplitude group map revealed positive activation bilaterally in one cluster covering mostly premotor and supplementary motor areas in BA6 (p < 0.005, t = 3.36).A glass brain and t-map are shown in Figure 1(a).Additionally, there was significant positive bilateral BOLD response in the superior and inferior parietal gyri, and a small cluster within the right putamen.No negative BOLD response was observed.

BOLD Flux Mapping
The first map based on a spatial analysis of the BOLD signal, the BOLD flux map (Figure 1(b)), revealed more widespread positive values with considerable overlap to premotor cortex and smaller overlap to areas located towards the interhemispheric fissure (SMA).In addition, there was congruency with superior and inferior parietal gyri.Notable additional increased flux during task was observed in the mid and anterior cingulate gyrus, the calcarine sulcus, cuneus and precuneus, and bilaterally in the superior temporal gyrus.In addition, there were larger clusters, now bilateral, in the putamen.Again, no negative response was observed.

BOLD Source Mapping
The second map based on a spatial analysis of the BOLD signal, the BOLD source map (Figure 1(c)), revealed both positive and negative values (see also Figure 2(a), right panel).The largest positive cluster (i.e., BOLD source) was located bilaterally in the putamen, and the second largest cluster in the SMA.Comparatively small clusters included locations in the left precentral gyrus, the right and left superior frontal gyrus, and the left superior and inferior parietal gyrus.Negative values (i.e., BOLD sink) were observed in small parts of the SMA but were located more towards adjacent white matter regions (Figure 2(a), right).Additional very small regions of negative values included the left postcentral gyrus, the precuneus, right Heschl's gyrus, and left amygdala.These small regions might be more likely to be false positives as we did not strictly correct for multiple comparison specific to source mapping but used the same p-values as for BOLD flux mapping, or leakage from the auditory stimuli.

Implications towards Conservation Laws Affecting Neuronal Dynamics
To summarize the findings: The conventional BOLD amplitude map is in agreement with previously published fMRI motor imagery experiments [15,[54][55][56][57][58][59][60]62,64,70].BOLD flux maps were more widespread than BOLD amplitude maps, but there were also regions in BOLD amplitude responses without corresponding increases in BOLD flux in the SMA.BOLD sources were mostly restricted to the putamen and SMA, whereas BOLD sinks were only sparsely distributed and small.
Taken together, it is not too evident that the sources and sinks can account for all the increased fluxes during task observed for example in the cingulate cortex.On the other hand, there were promising hints towards a continuity equation: Detectable sources in the SMA surrounded by flux, and also sources in the putamen.Overall, the observations do not seem to point strongly towards or against our main hypothesis of conserved quantities governed by a continuity equation.
There could be a number of reasons, including: First, it might be an interpretation difficulty as there is no experience with how the flux, source and sink patterns actually should appear in the brain.Second, the sought processes might be sub-threshold and not showing in this quite conservative statistical approach used.Also, the use of derivatives increases noise and might mask the sought effects.By using derivatives, we are actually doing the opposite of the conventional smoothing approach of fMRI data analysis, where smoothing is used in order to reduce the noise by averaging [10,71,72].Third, model inadequacy; it might be that some model assumptions, such as the homogeneity of the constant D, are not valid approximations, for example due to inhomogeneities in vascularity, neuronal density, etc.. Also, the model neglects energy storage within the brain, for example in the form of ionic gradients and releasable transmitter, such that models based on instantaneous differential balances of sources and sinks are too limited to fully embrace the local energy balance [73][74][75].In other words, homogeneous time-invariant model approach might be too simplistic.Fourth, the BOLD amplitude might not be the right proxy variable to use in such an analysis, or the BOLD amplitude as a nonlinear combination of other factors might mask the relevance of these factors.Therefore, other parameters such as cerebral blood flow, volume, or metabolic rate of oxygen might be more suitable.Fifth, even if there are conserved quantities, the description by a continuity equation might not be right or optimal.
Finally, it might be that conservation laws just play not a major role in the dynamics of the brain on this spatiotemporal scale.

BOLD Flux Imaging and Model-Free Approaches
Our finding that the flux appears to extend well beyond areas of primary relationship to the task in a group analysis is reminiscent of recent findings that if sufficient data is measured and analyzed in a model-free way, most of the brain appears to be involved even in simple tasks [76].Similar to that finding, our model is quite versatile as it depends on spatial correlations among voxels only, thereby dismissing information of the BOLD amplitudes themselves.More work is needed to relate BOLD flux imaging to those approaches.In addition, flux maps were derived under the assumption that flux values are distributed normally; however, due to their computation as vector norms flux distributions are skewed and it has not been investigated here if statistical tests that account for this will yield different results than the tests used here to compute flux maps.

Putamen
We found a strong BOLD source signal in particular in the putamen, meaning that during task signal fluctuations in putaminal areas become more homogeneous.This has been demonstrated recently in a similar experiment using local spatial synchronization fMRI (lssfMRI), where it was shown that during the task voxels appear more synchronized than during the control period [62].There, it has been argued that lssfMRI might be sensitive to spatial neuronal synchronization, specifically on spatial scales accessible to BOLD fMRI.A recent investigation of such spatial synchronization measured with electrodes in the striatum of macaques [77] might support this notion: These authors found a functional specialization of the putamen in the sense that pairs of putaminal medium spiny neurons were synchronized (by means of correlation) during a classical conditioning task.In contrast, neurons in the caudate and ventral striatum were not.Specifically, the electrodes covered a region with a spatial extension accessible to fMRI.

Relation to Resting State fMRI
The here presented approach of BOLD flux and source mapping specifically probes for spatial dependencies (computed as spatial derivatives) between nearest neighbor voxels, disregarding the BOLD amplitude itself.It therefore differs somewhat from correlation-based or independent component resting state analyses [6,7,9,78,79] in which correlations between BOLD amplitudes are measured.In addition, this approach is task dependent and looks at the differences between task and control states [80].A positive flux response here means that the norm of the BOLD gradient increases during task.A positive BOLD source signal means that the second spatial derivative decreases, indicating an increased homogeneity during the task period.

Differential Operators, Smoothing, and Nonlinear Transformations
The Laplacian is widely used in image analysis for blob detection.Blobs are areas of local homogeneity in pixel values.The here proposed source/sink estimation is based on applications of the Laplacian, too, before correlating the resulting time series voxelwise with the characteristic function of the task.From image analysis research it is known that the particular definition of the Laplacian actually pre-determines the local sensitivity, or the blob size to be detected, and it has been suggested to apply multi-scale approaches [81] to get a better picture of the blob distribution.In our application, the Laplacian has been defined by involving only nearest neighbors of the voxel of interest, thereby setting the scale or bandwidth accordingly.This can be generalized by using other spatial derivative estimators [82].This way, the results presented here might change and more work is in order to define the optimal definition of Laplace and gradient operators.
If applied to pre-smoothed data, derivatives partially revert smoothing, if the bandwidths of the differential operators and the smoothing operator are of comparable value.Then, the source maps as defined here would not necessarily provide much new information.Therefore, it is essential that no smoothing of the data is performed beforehand.
Another issue not addressed so far is how the nonlinear transformation to MNI space and the involved interpolation to an isotropic voxel size distort the derivative estimates.For example, it might be worthwhile to compute sources and fluxes before rather than after those transformations, and then transform the resulting maps to MNI space for the subsequent group analysis.
Finally, the local anatomy of the brain given by folded cortical layers etc. might dictate better local spaces to work with than the Cartesian grid approach used by us.

Possible Clinical Applications
Once conservation laws are established for normal, healthy brain dynamics, the hope is that they would inform our understanding of altered brain function.In pathological conditions one might then use conservation laws derived from imaging to detect how the brain copes with these conditions, which might then again enable imaging modalities with better discriminatory power to diagnose disease.Whereas this is a distant goal, already we found in our group study with healthy subjects that BOLD flux imaging was comprising larger areas than conventional BOLD amplitude maps, as well as larger significance measured by t-values (Figures 1(a,b)).
In order to provide an example of an individual clinical subject analysis, rather than a group study, here we consider a subject in the minimally conscious state (MCS).The MCS is a recovery step often following the vegetative state after severe brain injuries [83].Functional MRI in MCS can be used to locate residual cognitive capacities, for example by having subjects follow simple commands via mental motor imagery, providing data sets as used here [55,84].These experiments can assist diagnostics in order to distinguish MCS from other conditions [85] such as the locked-in syndrome [86], in which subjects are conscious but the ability to communicate is severely limited to absent.They also might assist clinical applications such as neural control engineering or planning for brain computer interfaces [87][88][89][90][91].In our previous MCS neuroimaging studies, we found that overall BOLD amplitude response was considerably weaker than in healthy control subjects, even when there was unambiguous response to certain tasks the patients were still able to perform.We also found dissociations between bedside examinations and functional imaging, which require further explanations [64].Modeling these disorders of consciousness goes hand in hand with advanced imaging technologies [43,92]; for example, the mesocircuit model for MCS predicts specific structural and dynamic changes occurring during recovery [93].In addition, cerebral autoregulation itself might provide markers for the important discrimination between the minimally conscious and vegetative state via observable effects of neurometabolic coupling [94].It would be very interesting to see if imaging conservation laws can shed more light on these models, or if at least it could provide us with more sensitive neuroimaging modalities.In order to demonstrate this point, BOLD amplitude mapping is compared with BOLD flux and source mapping in a single subject who was diagnosed with MCS (Figure 3).The task consisted of imagining swinging a tennis racket with the right hand.Whereas the amplitude map shows activation that is partially in agreement with the expected result for this motor imagery task, as in control subjects the flux and source maps appear entirely different to the amplitude map.In this case, there is only partial agreement with the control group results; for specifics, see caption of Figure 3.However, the brain injuries and accompanying morphologic changes make their exact location difficult to assess.In any case, the information obtained from functional scanning of subjects with disorders of consciousness generally is still quite limited, and any further information that can be extracted from this data could be very valuable.

Conclusions
Starting with the hypothesis that competition between brain regions affects the BOLD signal, and that this can be observed in the form of some conservation laws, we modeled the spatial dependences of the BOLD signal and its change with a cognitive task with a continuity equation.Flux and source terms in the equation where estimated from data by using spatial derivatives applied to the data.We demonstrated this procedure on a group of human subjects participating in a mental motor imagery fMRI experiment.We could not find clear evidence for a conservation law but introduced BOLD flux and source imaging as markers that might be sensitive to spatial communication between brain regions.The quest for finding the rule of conservation laws in the brain just started and hopefully will lead to new models for the complex functional organization of the brain.

Figure 1 .
Figure 1.BOLD amplitude, flux, and source maps.Comparison of BOLD amplitude (a), flux (b) and source (c) maps in a group of nine healthy subjects imagining themselves swimming while being imaged with an MRI scanner.Glass brains (left column) show maximum intensity projections of significant clusters.Negative-t clusters appear in c) only and are shown in a lighter shade.The three orthogonal cuts per brain (right column) show color coded t-values.The voxelwise thresholds are set to p < 0.005 (|t| > 3.36) in a two-sided mixed-effects t-test over the nine subjects, testing against the null hypothesis of zero activation.This corresponds to a false-discovery rate in the flux map of q = 0.04.

Figure 2 .
Figure 2. Sources and sinks, demonstrated on the supplementary motor area.(a) BOLD amplitude (left), flux (middle) and source map (right) in parts of the SMA, for the same group of subjects as shown in Figure 1.The source map shows regions with positive and negative values, i.e., sources and sinks, respectively.Color denotes t-values scaled in the same way as in the corresponding color bars shown in Figure 1 (|t| > 3.36, p < 0.005).(b) Vector plots to visualize the direction of the normalized flux vectors in the region indicated by the boxes in (a).In the left image, vectors are overlaid to the flux in gray scale denoting t-scores as indicated in scale on right of that panel; in the right image, vectors are overlaid to the sources.In the latter image, it is clearly discernible how the flux emanates from source regions (positive values, bright) and then proceeds to sink regions (negative values, dark).The tips of vectors are indicated by diamonds.

Figure 3 .
Figure 3. Application to a subject in the minimally conscious state.Images are from a 23-year-old man with a severe traumatic brain injury about four years before scanning.The anatomical image (see labels in figure) shows some cerebral atrophy as well as some image artefacts caused by a shunt.The BOLD amplitude map shows activation in SMA and premotor areas, among others.The BOLD flux map shows flux increase with task primarily in superior medial frontal cortex.The BOLD source map shows sources and sinks in the left caudate (|t| > 2.81, p < 0.005).