PhD Mathijs Raemaekers
email: m.a.h.l.l.raemaekers@uu.nl

Martin J M Lankheet

Sanne Moorman

Zoe Kourtzi

Richard J A Van Wezel

Directional anisotropy of motion responses in retinotopic cortex

Keywords: fMRI, motion, visual motion, MT, occipital fMRI, Motion, Vision, visual motion, MT, Occipital

teaching and research institutions in France or abroad, or from public or private research centers.

Introduction:

Visual motion processing is a fundamental capability of the visual system of humans and animals. Motion sensitive neurons can be found throughout the cortex including in visual areas V1, V2, V3 and the middle temporal visual area (MT+) [START_REF] Albright | Direction and orientation selectivity of neurons in visual area MT of the macaque[END_REF]Orban et al. 1986; Felleman and Van Essen 1987; [START_REF] Levitt | Receptive fields and functional architecture of macaque V2[END_REF]. A motion sensitive neuron selectively responds to a preferred direction and speed of motion, and motion direction is encoded in the pattern of activity across a pool of differently tuned neurons. It is generally thought that direction tuning of a neuron is directly related to its orientation tuning, or to that of its inputs, because only motion orthogonal to the preferred orientation can be signaled unambiguously.

Cells in area MT of the macaque indeed show such a close correspondence between orientation tuning and direction tuning [START_REF] Albright | Direction and orientation selectivity of neurons in visual area MT of the macaque[END_REF].

Recently, a functional magnetic resonance imaging (fMRI) study observed a large response bias for radial orientations as compared to tangential orientations in early visual cortex in humans and monkeys [START_REF] Sasaki | The radial bias: a different slant on visual orientation sensitivity in human and nonhuman primates[END_REF]). This finding leads to the expectation that motion sensitivity in early visual cortex should be anisotropic as well, with a bias towards tangential motion relative to the fovea. The possible presence of directional biases are further supported by human psychophysical studies [START_REF] Georgeson | Apparent foveofugal drift of counterphase gratings[END_REF][START_REF] Ball | Human vision favors centrifugal motion[END_REF][START_REF] Edwards | Asymmetries in the sensitivity to motion in depth: a centripetal bias[END_REF][START_REF] Raymond | Directional anisotropy of motion sensitivity across the visual field[END_REF][START_REF] Giaschi | The role of cortical area V5/MT+ in speedtuned directional anisotropies in global motion perception[END_REF]) and electrophysiological recordings in monkeys in MT [START_REF] Albright | Centrifugal directional bias in the middle temporal visual area (MT) of the macaque[END_REF]), parietal cortex [START_REF] Steinmetz | Functional properties of parietal visual neurons: radial organization of directionalities within the visual field[END_REF], and the frontal eye fields (FEF) [START_REF] Xiao | Radial motion bias in macaque frontal eye field[END_REF].

Directional anisotropies of the motion response are not necessarily based on local inhomogeneities in cells that are tuned for a specific motion direction. E.g. it is known that visual information can be integrated across hypercolumns by means of horizontal connections within early visual cortex. Horizontal connections extend between cells with a similar tuning in adjacent locations within the visual field [START_REF] Ts'o Dy | Relationships between horizontal interactions and functional architecture in cat striate cortex as revealed by cross-correlation analysis[END_REF]).

For motion perception, such cells could facilitate motion detection through recruitment effects along the path of motion (van Doorn and Koenderink 1984). If anisotropies are related to horizontal connections, then the length of the motion stimulus would have a strong effect on directional biases. The presence of any systematic directional biases would have important consequences for recent fMRI studies have used multivoxel pattern analysis of visual areas to predict perceived motion direction [START_REF] Kamitani | Decoding seen and attended motion directions from activity in the human visual cortex[END_REF]; Serences and Boynton 2007a; Serences and Boynton 2007b). Suggestions have been made that these decoding strategies exploit subvoxel inhomogeneities in direction sensitivity [START_REF] Kamitani | Decoding the visual and subjective contents of the human brain[END_REF][START_REF] Kamitani | Decoding seen and attended motion directions from activity in the human visual cortex[END_REF]. This notion implicitly assumes equal responses for all motion directions across the visual field. Large scale anisotropies in the motion response would have implications for this interpretation. To test whether early visual cortical areas change their response with motion direction, we measured the BOLD response with fMRI across the retinotopic visual field representations of V1, V2, V3 and MT+, in human subjects observing translating coherent motion for a full complement of motion directions.

The current study demonstrates that indeed systematic motion direction biases exist in early visual cortex in humans. Activation enhancements were observed for radial motion directions and the relative magnitude of the bias for centrifugal and centripetal motion changed as a function of visual eccentricity. This result is in conflict with a tangential motion direction bias that would be expected on the basis of a radial orientation bias. An additional experiment was performed to resolve this apparent conflict in results by testing the notion that directional biases are caused by anisotropies in the modulation of neuronal activation along the path of motion. The results from the additional experiment revealed that the observed motion direction biases most likely originate from anisotropies in long range horizontal connections within visual cortex.

Methods:

Subjects:

12 healthy subjects that where recruited from Utrecht University participated in the experiment. All subjects gave informed consent for participation (approved by the local ethics committees of Utrecht University or the University of Birmingham). 9 subjects performed the fMRI experiment at Utrecht University, and 3 subjects performed the experiment at the University of Birmingham with concurrent measurement of eye movements. 10 subjects performed the additional experiment.

Stimuli:

For task presentation we used a desktop PC, a projection screen, and a video-projector system. All stimuli were programmed in C++ software (Bjarne Stroustrup, 1983, Bell Laboratories, USA). The start of each series of stimuli was triggered by the scanner. During all stimuli, there was a red central fixation dot (radius of 0.08° visual angle) that was surrounded by a circular aperture (radius of 0.4° visual angle). Subjects were requested to maintain fixated on the fixation dot regardless of the presented stimuli. The average luminance of the entire screen was constant during all stimuli and was 42.2 cd/m 2 .

For eccentricity mapping, we used an expanding ring with a maximum eccentricity of 7.5° visual angle (Figure 1A). After the ring was fully expanded, it returned to its minimum eccentricity (0.4° visual angle).

The width of the ring was 1/5 th of the maximum stimulus radius. There was 1 series (900 images) with 8 cycles of 54000 ms (100 images) and there was a blank period with only the fixation dot during the first and last 27000 ms (50 images) of the series. For polar angle mapping, we used a rotating wedge (45° circular angle) that extended to a maximum eccentricity of 7.5° visual angle (Figure 1B). There was 1 series (900 images) with 8 full clockwise rotations that lasted 54000 ms (100 images) each. The screen was The main stimulus was a random dot pattern with dots that translated at a constant speed of 4.81°/s within a circular aperture with a radius of 6.08° visual angle (Figure 1C). The direction of translation changed 360° in cycles of 54000 ms (100 images). During one series the direction of translation rotated clockwise, in the other series it rotated counter clockwise. There were 8 full rotations in translation direction during each series of 900 images. The dots were static during the first and last 27000 ms (50 images) of each series.

The control stimulus consisted of 2 random dot patterns, one in the left, and one in the right hemifield, within a circular aperture with a radius of 6.08° visual angle (Figure 1D). In one condition, the dots of the two patterns were moving outwards, parallel to the horizontal meridian (always 4.81°/s). In the other condition, they were moving inwards, parallel to the horizontal meridian. There were 8 blocks of 19980 ms (37 images) of each condition which were orderly alternated and separated by 12420 ms (23 images) periods where all dots were static. Responses were measured along the horizontal and vertical meridian.

The control stimulus is less likely to induce shifts in spatial attention or eye movements than the main stimulus.

The stimulus of the additional experiment was largely similar to the stimulus of the main experiment. The important difference was an occlusion of part of the random dot pattern with grey bars with a width of 0.71° each, and a 1.43° gap between bars (Figure 1E). The random dot pattern was moving in blocks of 16200 ms (30 images) with a constant speed (4.81°/s) and direction. The direction of motion could be upward, downward, leftward, or rightward and blocks of motion were always interrupted with 16200 ms (30 images) periods where the dots were static. In the first session (960 images) the direction of motion was always parallel to the bars, in the second session (960 images) the direction of motion was always orthogonal to the bars. Furthermore, the radius of the stimulus was increased to 7.50°, as results from the main experiment indicated that there were no spatial shifts in representation of the stimulus (see further below). During the additional experiment, there was an attentional control task within the aperture at central fixation. A white cross (0.25° visual angle) was presented every 1000 ms for a duration of 300 ms.

During 25% of the presentations, one of the bars of the cross became an arrow (up, down, left, right) by the placement of a small rectangle (0.05° visual angle) at one of the extremities. Subjects had to press the button on a button box (up, down, left, right) that corresponded to the indicated direction of the arrow.

Random dot patterns contained an equal number of white and black dots (0.14° visual angle). On average, there were 4300 dots visible at a particular moment in time. When dots reached the border of the pattern, they were replotted at the other extremity.

Eye movements:

Eye movements were measured at 60 Hz during the main stimulus in 3 out of 12 subjects and using the ASL EyeTrac 6 MR-compatible video based eye tracker (Applied Science Laboratories, Bedford, USA) at the University of Birmingham. Blinks and low frequency drifts (<0.00926 Hz. ~ 2 cycles of the main stimulus) were removed. At each time point during the stimulus, eye position and displacement were recalculated in terms of eye position and displacement relative to the motion direction of the stimulus.

Saccades were detected by identifying any sample in which the velocity exceeded 30°/s.

Statistical Analysis:

All preprocessing steps were done using SPM5 (http://www.fil.ion.ucl.ac.uk/spm/). After realignment, the functional images were coregistered and resliced to the FA27 volume, using the first functional volume as a source. The T1 structural image was also coregistered to the FA27-image, thereby providing spatial alignment between the structural image and the functional volumes. Low frequency noise in the fMRI data was modeled and removed from the data using a GLM and a design matrix containing the mean of each image, and 8 cosine functions forming a high-pass filter with a cutoff at 8.2*10 -1 Hz. For polar angle and eccentricity mapping, a vector was created that represented cyclic activation during presentation of wedges and rings (7200 ms activation every 54000 ms) and was convolved with a hemodynamic response function [START_REF] Friston | Characterizing evoked hemodynamics with fMRI[END_REF]. The cross correlation between the fMRI data and this vector was calculated for every voxel and for 100 lags (0-99; i.e. every image within a cycle) and the peak cross correlation determined the receptive field location of the voxel in polar angle and eccentricity. Polar angle and eccentricity measures of voxels were interpolated to 8 steps (45° circular angle) and 5 steps (1.42° visual angle) respectively. The activated voxels (p<0.05; Bonferroni corrected) formed a visual field representation that consisted of 40 segments (8 x 5), which was further subdivided in V1, V2, V3, and MT+ (see further below). For each subject, the average BOLD response was calculated in these segments during cycles of the main stimulus and the second control stimulus. Average BOLD responses were shifted 7 images back in time to account for the delay of the hemodynamic response. Subsequently, responses for clockwise and counterclockwise changing translation direction were averaged and rescaled to 'translation direction' over the x-axis (0°-360°), and 'percent signal change' over the y-axis. Effects of translation direction on the BOLD response were analyzed for each segment using repeated measures General Linear Models (100 layers ~ number of images in each cycle) and results were Bonferroni corrected for the total number of tested segments in each comparison (p<0.05).

Average BOLD responses during cycles of the control stimulus and the additional experiment were also generated for each segment. The amplitude of the responses were estimated by fitting a BOLD response model of block activation followed by a rest period [START_REF] Friston | Characterizing evoked hemodynamics with fMRI[END_REF].

Segmentation of retinotopic areas:

The T1 image was corrected for intensity inhomogeneities using the segmentation utility in SPM5 [START_REF] Ashburner | Unified segmentation[END_REF]. The bias corrected T1 images were then imported in the Computerized Anatomical Reconstruction and Editing Tool Kit (CARET) [START_REF] Van Essen | An integrated software suite for surface-based analyses of cerebral cortex[END_REF]). T1 images were resliced to a 1 mm isotropic resolution, manually placed in Talairach orientation, and subdivided in left and right hemisphere. All subsequent steps were done per hemisphere. The intensity of the grey/white matter border was determined, followed by automatic extraction of the cortex. A white matter segment was For identification of MT+, the voxel time-series of the control stimulus were entered into a GLM, with a one factor representing BOLD activation during inward motion, and one factor representing BOLD activation during outward motion. A volume of t-values was generated for the contrast representing activation during both inward and outward motion. MT+ was defined as the bilateral clusters of significant voxels in the ascending limb of the inferior temporal sulcus (p<0.05; Bonferroni corrected)

Results

:

Main Experiment:

We presented translating motion that slowly changed direction (Figure 1C), and measured the BOLD response across the visual field representations of V1, V2, V3, and MT+, which were all subdivided in 40 segments (8 polar angles * 5 eccentricities). All except 3 of the total of 120 segments in V1, V2 & V3 demonstrated a systematic group-wise effect on brain activation (F (99,1089) >1.60; p<0.05). Thus, there were systematic anisotropies in V1, V2, and V3. However, there were no significant effects anywhere in MT+.

To further investigate the nature of the anisotropies, we averaged over eccentricities, and found a significant effect of motion direction in all 8 polar angle segments of V1, V2 and V3 (F (99,1089) >1.51; p<0.05), but not in MT+. The mean BOLD responses for each motion direction in the 8 polar angle segments of V1, V2, and V3 are depicted in figure 2. It can be seen that the direction bias in the BOLD response corresponds to enhanced radial motion sensitivity. The direction bias regularly rotates with polar angle relative to the fovea. These results show that the motion direction relative to the fovea is an important For investigating effects of eccentricity on this bias, the BOLD responses per segment where first rescaled to motion direction relative to polar angle over the x-axis (0º=centripetal motion, -90º/90º=tangential motion, -180º/180º=centrifugal motion), and subsequently averaged over the 8 polar angles per eccentricity.

The average responses per eccentricity and visual area are depicted in figure 3. There were significant effects of motion direction for all eccentricities of V1, V2 & V3 (F (99,1089) >1.47; p<0.05), and these effects consisted of either a centrifugal and/or centripetal response biases. Again, no significant effects were observed in MT+.

To estimate how the amplitudes of the centrifugal and centripetal response biases changed with eccentricity in the retinotopic areas, the responses in each subject were fitted with a linear combination of two peaks (based on a cosine function) using a multiple regression analysis. The amplitude fits of the centrifugal and centripetal response are depicted in figure 4. Effects of eccentricity (5 layers; 5 steps of 1.42°), motion direction (2 layers; centrifugal/centripetal motion), and retinotopic area (3 layers; V1, V2 & V3) were analyzed using a repeated measures general linear model. The test revealed a main effect of motion direction (F (1,11) =30.99; p<0.001), and an interaction effect between eccentricity and motion direction (F (4,44) =51.33; p<0.001). Centripetal biases were larger than centrifugal biases, and this difference between centrifugal and centripetal biases depended on eccentricity. Between 0.40° and 1.82°, and between 4.66° and 7.50°, there was a combination of a centrifugal and centripetal bias. Between 1.82° and 4.66°, there was a large centripetal bias and nearly absent centrifugal bias.

Control Experiment:

As the full-field stimulus (Figure 1C) that we used could have induced eye movements or shifts in spatial attention, we investigated whether the same effects were also present during the control stimuli (Figure 1D and 1E). In the control stimulus the visual field was bisected across the vertical meridian and we showed opponently moving random dot patterns in the two halves. Together with the absence of motion in the central 0.40° degrees this minimizes consistent effects of eye movements or shifts in attention. We investigated the BOLD responses in the segments along the horizontal meridian during this stimulus (Figure 5). Amplitude fits of the BOLD response were entered into a repeated measures GLM (motion direction * eccentricity * retinotopic area; similar as for the main stimulus). The results followed the same pattern as for the main stimulus. This consisted of larger responses for centrifugal motion at the lowest eccentricity, and larger responses for centripetal biases at higher eccentricities. The differences between centrifugal and centripetal motion were somewhat larger during the control stimulus than during the main stimulus (88% in V1, 35% in V2, and 19% in V3). Again, there were no differences in MT+.

Additional Experiment:

In the additional experiment, the motion stimulus was interrupted either parallel to the path of motion or orthogonal to the path of motion (1E). BOLD responses were measured along the horizontal and vertical meridians for motion parallel to the bars and motion orthogonal to the bars. BOLD responses were classified as centrifugal, centripetal, or tangential, dependent on the position in the visual field relative to the motion direction. BOLD responses were averaged over retinotopic areas (V1, V2 & V3), as the results of the main experiment and control experiment demonstrated no qualitative differences in directional biases between visual areas. Amplitude fits of the BOLD response were entered into a repeated measures GLM (motion direction * eccentricity * condition (parallel/orthogonal interruption). There was a significant interaction effect between motion direction and condition (F (2,72) =17.44; p<0.001) and a significant interaction effect between motion direction, eccentricity, and condition (F (8,72) =3.24; p=0.003). These effects were best described as the presence of a centripetal motion direction bias when the stimulus was interrupted parallel to the path of motion, but absent direction biases when the stimulus was interrupted orthogonally to the path of motion (except for the lowest eccentricity, where there were no biases in both conditions).

Furthermore, the performance on the attentional control task during the additional experiment was not different for the stimulus with interruptions parallel or orthogonal to the path of motion considering reaction times (mean RT parallel ± SD =782 ± 99 ms; mean RT orthogonal ± SD = 776 ± 101 ms; paired F o r P e e r R e v i e w t (7) =1.20; p=.27) or the proportion of correct responses (mean correct parallel ± SD = 0.86 ± 0.11; mean correct orthogonal ± SD = 0.87 ± 0.09; paired t (7) = -0.70; p=0.50). Subjects reported that they were fully engaged in optimally performing the task. Behavioral data of two subjects was lost due to a defective buttonbox.

Eye movement control:

In addition to the split-field stimulus as a control we also measured eye-movements for three subjects during the main experiment. Eye movement recordings showed no deviations in fixation position with changing motion direction (mean deviation from central fixation was 0.04° visual angle, in the direction of '+14° of polar angle relative to the motion direction'). In addition, no saccades (see Methods) were detected in all three subjects for the entire duration of the stimulus.

To further analyze if there was a systematic relationship between motion direction and direction of microsaccades or smooth pursuit (e.g. as a part of optokinetic nystagmus), the absolute displacement parallel to the direction of motion and perpendicular to the direction of motion was calculated between each sample. Nystagmus e.g. would on average induce a larger amount of eye displacement parallel to the path of motion compared to orthogonal to the path of motion. This analysis revealed a mean displacement during the stimulus of 0.12°/sample parallel to the motion direction, and 0.12°/sample perpendicular to the direction of motion, indicating no relation between direction of eye movements and motion direction. The three subjects with simultaneous eye movement acquisition had fMRI results that were largely similar to the group mean (Figure 4)

Discussion:

We found strong evidence for directional anisotropies of motion responses across the visual field representations of V1, V2 and V3 when observing whole field coherent translating motion. In all three retinotopic areas, the amplitude of the motion response depended on the motion direction relative to polar angle and the eccentricity relative to the fovea. All biases consisted of enhanced responses for centrifugal and centripetal motion, relative to tangential motion. The centripetal bias was present for all except very low eccentricities. The centrifugal bias was present for low eccentricities, and for eccentricities beyond ±4.5º. The response to centripetal motion was larger than the response to centrifugal motion, except for low eccentricities. The results from the control experiment, where the left and right visual field contained opponent motion directions, confirmed the results from the main experiment. We did not observe any significant motion direction biases in MT+. Our additional experiment showed that these directional biases did not occur when the motion stimulus was interrupted orthogonally to the path of motion.

On first sight, the current finding of a radial motion bias is unexpected. Motion perception of gratings and random dot patterns are both thought to be based on orientation detection in space and time [START_REF] Skottun | On the directional selectivity of cells in the visual cortex to drifting dot patterns[END_REF]). As it is known that the preferred motion directions of cells is orthogonal to their orientation preference [START_REF] Albright | Direction and orientation selectivity of neurons in visual area MT of the macaque[END_REF], the previously reported radial orientation bias should therefore predict a tangential motion bias [START_REF] Sasaki | The radial bias: a different slant on visual orientation sensitivity in human and nonhuman primates[END_REF]). In our study we found the opposite. However, the additional experiment shows that it is highly unlikely that the directional biases are related to local anisotropies in orientation sensitivity. If the currently observed directional biases were based on classical receptive field effects, then the additional stimulus would have produced a directional bias whether the stimulus was interrupted parallel or orthogonal to the direction of motion. Even when the stimulus was interrupted orthogonally, the length of each motion path was 1.43º, thereby still far exceeding receptive field sizes in the central 15º of the visual fields in V1 and V2 [START_REF] Dumoulin | Population receptive field estimates in human visual cortex[END_REF]. Directional biases however completely disappeared when the length of the motion path was 1.43º.

Anisotropies in extraclassical receptive field effects are thus producing the bias, which means either long range horizontal (lateral) connections within the visual areas, or interareal feedback connections. Although fMRI does not have the temporal resolution to completely resolve this issue, we believe that the data from the additional experiment indicate that anisotropies in horizontal connections are the most likely candidate.

Interareal feedback signals would originate from areas upstream containing motion sensitive neurons with large receptive field sizes. Neurons that are responding to motion in a large portion of the visual field would be largely indifferent to whether motion is interrupted parallel or orthogonal to the path of motion, as the total amount of motion information across large portions of the visual field remains the same. Anisotropies were however only observed for interruptions parallel to the path of motion. Horizontal Their role in visual motion processing has been far less investigated. These lateral connections could play a role in motion information integration across regions beyond the classical receptive field sizes, which could contribute to solving the aperture problem. Guo et al. found activation modulation in V1 neurons when the stimulation within the classical receptive field was identical, but the global direction of the stimulus movement was different [START_REF] Guo | Primary visual cortex neurons that contribute to resolve the aperture problem[END_REF]. Such connections could also play a role in defining shapes of objects that are defined by motion. The current study would predict that the strength of modulations through horizontal connections would be different along the path of motion compared to orthogonal to the path of motion. Moreover, our study suggests that modulations along the path of motion are anisotropic across the visual field.

A couple of factors could have confounded the current results. Although eye movements are unlikely to have accounted for the data, the motion stimulus could have induced covert shifts in spatial attention. Spatial attention is known to locally enhance the visual response at the focus of attention [START_REF] Tootell | The retinotopy of visual spatial attention[END_REF][START_REF] Jack | Separate modulations of human V1 associated with spatial attention and task structure[END_REF]). The single translating random dot pattern may have induced drifts in the focus of spatial attention, which could account for the currently observed data. However, there are several arguments that oppose such an explanation. Biases were also present when subjects were fully engaged in a in an attentional control task. Furthermore, attentional shifts would most likely enhance the signal in the direction of motion as a result of attentional tracking, and not opposite to the direction of motion as was observed in this study. In addition, during the control stimulus where two opponently moving random dot patterns were presented simultaneously, there were similar differences between centrifugal and centripetal motion, although such a stimulus would be far less effective in inducing attentional shifts. Although it is technically possible that this stimulus would have induced two separate attentional spotlights, divided attention roughly halves the amount of activation per spotlight [START_REF] Mcmains | Multiple spotlights of attentional selection in human visual cortex[END_REF]. This would result in lower biases for the control stimulus compared to the main stimulus, but we observed the opposite. Furthermore, the additional experiment produced only biases when the stimulus was interrupted parallel to the path of motion, whereas both the stimulus with orthogonal and parallel interruptions would to some extent be effective in inducing attentional shifts.

Psychophysical studies have extensively addressed direction biases in motion sensitivity. These studies have predominantly tested and found differences between centrifugal and centripetal motion. While some studies have reported increased sensitivity for centripetal motion [START_REF] Edwards | Asymmetries in the sensitivity to motion in depth: a centripetal bias[END_REF][START_REF] Raymond | Directional anisotropy of motion sensitivity across the visual field[END_REF][START_REF] Giaschi | The role of cortical area V5/MT+ in speedtuned directional anisotropies in global motion perception[END_REF]), other studies have reported a centrifugal bias [START_REF] Scott | Directional asymmetry of motion aftereffect[END_REF][START_REF] Georgeson | Apparent foveofugal drift of counterphase gratings[END_REF] Harris 1978; Edwards and Badcock 1993; [START_REF] Lewis | Bias to experience approaching motion in a three-dimensional virtual environment[END_REF]. However, see also van de Grind et al., who observed a tangential motion bias at eccentricities beyond 24º when measuring detection thresholds for different motion directions across the visual field (van de [START_REF] Van De Grind | Inhomogeneity and anisotropies for motion detection in the monocular visual field of human observers[END_REF]). This apparent conflict in results may originate from differences in the eccentricity at which the motion direction biases are measured. Centripetal biases have been reported for stimuli that were presented at relatively low eccentricities and the centripetal bias decreases with increasing eccentricity [START_REF] Edwards | Asymmetries in the sensitivity to motion in depth: a centripetal bias[END_REF].

Contrary to the centripetal bias, the centrifugal bias increased with increasing eccentricity [START_REF] Ball | Human vision favors centrifugal motion[END_REF]. The larger centripetal bias in this study was observed at relatively low eccentricities, which is thus in line with previous psychophysical studies.

Surprisingly, in our study we did not observe any effects of motion direction on activation in MT+ in spite of it's well known involvement in perception of coherent motion [START_REF] Braddick | Form and motion coherence activate independent, but not dorsal/ventral segregated, networks in the human brain[END_REF][START_REF] Braddick | Brain areas sensitive to coherent visual motion[END_REF]). Contrary to our findings, Giaschi et al. [START_REF] Giaschi | The role of cortical area V5/MT+ in speedtuned directional anisotropies in global motion perception[END_REF]) measured activation differences in MT+ during horizontal and vertical centrifugal and centripetal motion and observed a centripetal response bias. This contradiction could be the result of the difficulties that are present when defining the retinotopy of MT+ [START_REF] Smith | Sensitivity to optic flow in human cortical areas MT and MST[END_REF]). However, the control stimulus that was used in the current experiment was roughly similar to a stimulus that was used by Giaschi et al., and measurement of centrifugal vs. centripetal response biases with this stimulus does not depend on the retinotopy of MT+. Still, our control stimulus did not produce a centripetal response bias. A more likely other explanation is that the motion stimuli of Giaschi et al. were presented at 25% and 75% coherence level respectively, compared to 100% in our study. The centripetal motion bias that they observed in MT+ was only significant for stimuli near the detection threshold (25% coherence). Detection of motion at low coherence levels may require attention to motion which enhances activation levels in MT+. Thus direction biases could very well be present in MT+, but not for fully coherent motion stimuli as in our experiment. For early visual cortex, we would expect no qualitatively different results for incomplete coherence as these areas are thought to be less influenced by attention. Furthermore, incomplete coherence (and also limited lifetime) of dots does not truly interrupt a motion stimulus. An uninterrupted motion stimulus will produce the directional bias.

The currently observed directional anisotropy of the motion response has implications for the interpretation of earlier studies in which perceived motion directions were predicted with multivoxel pattern analysis Tong have argued that the use of multivoxel pattern analysis provides the ability to detect activation differences as a result of random fluctuations in the number of columns within a voxel that are tuned for a particular motion direction [START_REF] Kamitani | Decoding the visual and subjective contents of the human brain[END_REF][START_REF] Kamitani | Decoding seen and attended motion directions from activity in the human visual cortex[END_REF]. Although the acquired form of our data did not allow for multivoxel pattern analysis, our results suggest that large scale anisotropies across the visual field may have contributed to decoding in these studies. This view is further supported by the finding that discrimination with support vector machines was worse for opposing motion directions, and worse in MT+ [START_REF] Kamitani | Decoding seen and attended motion directions from activity in the human visual cortex[END_REF]. Both observations are in line with the notion that the currently observed anisotropies have contributed to the results of previous decoding studies, in addition to possible subvoxel inhomogeneities. The removal of these global effects requires a full correction per presented motion direction for the mean activation within each eccentricity by polar angle segment of the visual field.

In conclusion, we found consistent anisotropies of the BOLD response across the visual field representations of V1, V2 & V3, which are most likely related to anisotropies in horizontal connections.

Horizontal connections could thus make a substantial contribution to the fMRI signal in visual cortex. It There was an effect of motion direction on the BOLD response in all these segments (p<0.05; Bonferroni corrected). The minimum in each graph corresponds to the lowest response within each segment. Gridlines in the radial plots indicate 0.5 % signal change relative to the minimum BOLD response. All average responses could be fitted by a combination of a centrifugal and centripetal response bias. Note that these plots are not meant to represent actual direction tuning curves of the segments. Due to the dispersion of the BOLD response relative to the neural response, underlying neural response biases may be more specifically tuned. In addition, the plotted BOLD response is not absolute, but relative to the minimum response in each segment. In MT+ (not displayed) there were no significant effects of motion direction in any of the 8 polar angle segment.

Human Brain Mapping

  At both locations scanning was performed on a Philips Achieva 3T scanner (Philips Medical Systems, Best, the Netherlands) with a Quasar Dual gradient set. For functional images, a navigated 3D-PRESTO pulse sequence was used (van Gelderen et al. 1995; Ramsey et al. 1998). The acquisition parameters were: TR=30 ms (time between 2 subsequent RF pulses); effective TE=43.87 ms; FOV(anterior-posterior, inferior-superior, right-left)= 65*200*160 mm; flip angle=10 degrees; matrix: 26*80*64 slices; voxel size 2.5 mm isotropic; 8 channel head coil; SENSE factors=2.0 (left-right) and 1.8 (anterior-posterior). A new volume was acquired every 540 ms, and encompassed the posterior 65 mm of the brain. A T1 weighted structural image of the whole brain (voxel resolution=0.875*0.875*1.00 mm; FOV=168*224*160 mm) was acquired at the end of the functional series. Immediately before the T1 image, an additional PRESTO image of the same volume of brain tissue was acquired with a high flip-angle (27 degrees, FA27) for the image coregistration routine (see below).

  first and last 27000 ms (50 images) of the series, except for the central fixation dot. Both the rotating wedge and the expanding ring contained a checkerboard pattern with white and black squares that reversed contrast every 125 ms.

  automatically corrected for topological errors. Remaining topological errors were removed manually. A surface reconstruction was generated and inflated. Several cuts were applied on the inflated surface, amongst others along the calcarine fissure and the medial wall. The surface was flattened and geometric distortions were reduced. Results of the polar angle mapping were mapped to the anatomical surface using the 'enclosing voxel' algorithm. Retinotopic areas V1, V2, and V3 were manually segmented by drawing borders along the reversals in the change of the polar angle representation. The resulting flat segments were converted back to volumetric format and used as ROIs in further analysis.

  of the amplitude of the BOLD response. In areas V1, V2 and V3 BOLD responses were enhanced for centrifugal and centripetal motion, relative to tangential motion. In area MT we found no significant variation of the BOLD response as a function of motion direction.

  better suited to generate modulations that are specifically dependent on the length of the path of motion. Horizontal connections in visual cortex have been implicated in the facilitation of the detection of spatial contours based on signals from local orientation detectors (Gilbert and Wiesel 1983; Schmidt et al. 1997).

(

  Kamitani and Tong 2006; Serences and Boynton 2007a; Serences and Boynton 2007b). Kamitani and

Figure 1 :

 1 Figure 1: Schematics of the presented stimuli: Colored arrows and lines are for illustration purposes only. A: An expanding circle for eccentricity mapping. B: A clockwise rotating wedge for polar angle mapping. C: A moving random dot pattern which slowly changed translation direction in cycles of 54000 ms. During the first functional series the translation direction of dots changed clockwise (as indicated by the dotted green arrows), in the second series it changed counter clockwise D: Two random dot patterns, one in each hemifield (the dotted green line marks the border). During one condition, dots were moving outwards (blue arrows), during the other condition they were moving inwards (red arrows), parallel to the horizontal plane. E: A random dot pattern that was regularly interrupted. In one condition the interruptions were parallel to the path of motion (blue arrow), in the other condition they were orthogonal to the path of motion (red arrow).

Figure 2 :

 2 Figure 2: Mean (n=12) BOLD responses during 360° cycles of motion directions for the 8 polar angle segments of V1, V2 & V3. There was an effect of motion direction on the BOLD response in all these segments

Figure 3 :

 3 Figure 3:

Figure 4 :

 4 Figure 4: Average amplitudes (n=12) of the centrifugal and centripetal response biases as a function of eccentricity in areas V1, V2 & V3 as estimated by fitting cosine functions. Bars indicate standard error of the mean. Dotted lines represent the mean effects of the 3 subjects that performed the task with simultaneous acquisition of eye movement data.

Figure 5 :

 5 Figure 5: Average BOLD responses during presentation of the control stimulus (Figure 1D) for different eccentricities in the segments along the horizontal meridian in V1, V2, V3 & MT+ (n=12). During the first 19980 ms of each series, the two random dot patterns where either moving inward or outward. During the last 12420 ms, the random dot patterns were static. Error bars indicate standard error of the mean. The eccentricity between 1.82° and 3.24° of MT+ was not represented in all subjects, resulting in unreliable estimates.

Figure 6 :

 6 Figure 6: Average BOLD responses (n=10) during presentation of the additional stimulus (Figure 1E) for different eccentricities for motion parallel to the bars (6A) and motion orthogonal to the bars (6B). Responses are the average over all meridians (horizontal & vertical) and visual areas (V1, V2, V3). During the first 16800 ms

John Wiley & Sons, Inc.Human Brain Mapping


Julien Bect

A unifying formulation of the Fokker-Planck-Kolmogorov equation for general stochastic hybrid systems 6,66

Keywords: Stochastic hybrid systems, Stochastic system with jumps, Markov processes, Fokker-Planck equation

A general formulation of the Fokker-Planck-Kolmogorov (FPK) equation for stochastic hybrid systems is presented, within the framework of Generalized Stochastic Hybrid Systems (GSHSs). The FPK equation describes the time evolution of the probability law of the hybrid state. Our derivation is based on the concept of mean jump intensity, which is related to both the usual stochastic intensity (in the case of spontaneous jumps) and the notion of probability current (in the case of forced jumps). This work unifies all previously known instances of the FPK equation for stochastic hybrid systems, and provides GSHS practitioners with a tool to derive the correct evolution equation for the probability law of the state in any given example.

Introduction

Among all continuous-time stochastic models of (nonlinear) dynamical systems, those with the Markov property are especially appealing because of their numerous nice properties. In particular, they come equipped with a pair of operator semigroups, the so-called backward and forward semigroups, which are the analytical keys to many practical problems involving Markov processes. When the system is determined by a stochastic differential equation, these semigroups are generated by Partial Differential Equations (PDE) -respectively the backward and forward Kolmogorov equations. The forward Kolmogorov PDE, also known as the Fokker-Planck equation, rules the time evolution t → µ t , where µ t is the probability distribution of the state X t of the system at time t. This paper deals with the generalization of this Fokker-Planck-Kolmogorov (FPK) equation to the framework of General Stochastic Hybrid Systems (GSHSs) recently proposed by Bujorianu and Lygeros [START_REF] Bujorianu | General stochastic hybrid systems: modelling and optimal control[END_REF][START_REF] Bujorianu | Toward a general theory of stochastic hybrid systems[END_REF].

The GSHS framework encompasses nearly all continuous-time Markov models arising in practical applications, including piecewise deterministic Markov processes [START_REF] Davis | Piecewise deterministic Markov processes: a general class of nondiffusion stochastic models[END_REF][START_REF] Davis | Markov Models and Optimization[END_REF], switching diffusions [START_REF] Ghosh | Optimal control of switching diffusions with application to flexible manufacturing systems[END_REF][START_REF] Ghosh | Ergodic control of switching diffusions[END_REF] and the stochastic hybrid systems of Hu et al. [START_REF] Hu | Towards a theory of stochastic hybrid systems[END_REF]. The reader is referred to [START_REF] Bujorianu | A stochastic hybrid system modeling framework[END_REF][START_REF] Pola | Stochastic hybrid models: an overview with application to air traffic management[END_REF] for a detailed overview of these classes of models with a view towards applications in Air Traffic Management (ATM). Two kinds of jumps are allowed in a GSHS: spontaneous jumps, defined by a state-dependent stochastic intensity λ(X t ), and forced jumps triggered by a so-called guard set G. Generalized FPK equations have been given in the literature, in the case of spontaneous jumps, for several classes of models; see Gardiner [START_REF] Gardiner | Handbook of Stochastic Methods for Physics[END_REF], Kontorovich and Lyandres [START_REF] Kontorovich | Dynamic systems with random structure: an approach to the generation of nonstationary stochastic processes[END_REF], Krystul et al. [START_REF] Krystul | Risk decomposition and assessment methods[END_REF] and Hespanha [START_REF] Hespanha | A model for stochastic hybrid systems with application to communication networks[END_REF] for instance. The case of forced jumps is harder to analyze, at the FPK level, because no stochastic intensity exists for these jumps. Until recently, the only results available in the literature were dealing with one-dimensional models; see Feller [START_REF] Feller | The parabolic differential equations and the associated semi-groups of transformations[END_REF][START_REF] Feller | Diffusion processes in one dimension[END_REF] and Malhamé and Chong [START_REF] Malhamé | Electric load model synthesis by diffusion approximation of a high-order hybrid-state stochastic system[END_REF]. These results have been extended to a class of multi-dimensional models by Bect et al. [START_REF] Bect | Generalized Fokker-Planck equation for piecewise-diffusion processes with boundary hitting resets[END_REF].

The main contribution of this paper is a general formulation of the FPK equation for GSHSs. It is based on the concept of mean jump intensity, which conveniently substitutes for the stochastic intensity when the latter does not exist. This equation unifies all previously known instances of the FPK equation for stochastic hybrid systems, and provides GSHS practitioners with a tool to derive the correct evolution equation for the probability law of the state in any given example. The results presented in this paper are extracted from the PhD thesis of the author [START_REF] Bect | Processus de Markov diffusif par morceaux: outils analytiques et numériques[END_REF].

The paper is organized as follows. Section 2 introduces our notations for the GSHS formalism, together with various assumptions that will be needed in what follows. Section 3 defines the crucial concept of mean jump intensity, which is used in Section 4 to derive our unified measure-valued formulation of the generalized FPK equation for GSHSs. Section 5 shows that the measure-valued equation of Section 4 yields an evolution equation with associated boundary conditions in the case where a piecewise smooth exists. Section 6 provides several examples, showing that the generalized FPK equation allows to recover all known instances of the FPK equation for stochastic hybrid systems. Finally, Section 7 concludes the paper and provides directions for future work.

Remark 1. The stochastic processes that we call GSHSs, following the papers by Bujorianu and Lygeros [START_REF] Bujorianu | General stochastic hybrid systems: modelling and optimal control[END_REF][START_REF] Bujorianu | Toward a general theory of stochastic hybrid systems[END_REF], are also called GSHPs in [START_REF] Bujorianu | A stochastic hybrid system modeling framework[END_REF] -where GSHP stands for "General Stochastic Hybrid Process". Note also that the terms GSHS / GSHP are used by [START_REF] Everdij | Hybrid Petri nets with diffusion that have into-mappings with generalised stochastic hybrid processes[END_REF] to make a clear distinction between the formal data defining the process and the process itself. We shall stick to the "GSHS" terminology in this paper.

General Stochastic Hybrid Systems

From the probabilistic point of view, the object of interest in the GSHS formalism is a continuous-time strong Markov process X = (X t ) t≥0 , with values in a metric space E 0 . It is defined on a filtered space (Ω, A, F ), equipped with a system P x ; x ∈ E 0 of probability measures on (Ω, A), with the property that X starts from x under P x for all x ∈ E 0 . As usual in the theory of Markov processes, E x denotes the expectation operator corresponding to P x . The reader is referred to [START_REF] Blumenthal | Markov Processes and Potential Theory[END_REF][START_REF] Dynkin | Markov Processes[END_REF][START_REF] Sharpe | General Theory of Markov Processes[END_REF] for background information on continuous-time Markov processes.

It is assumed that, for each ω ∈ Ω, the samplepath t → X t (ω) is right-continuous, has left limits X - t (ω) in the completion E of E 0 , and has a finite number of jumps, denoted by N t (ω), on the interval (0; t] for all t ≥ 0. The last condition can be seen as a "pathwise non-Zenoness" requirement. We will denote by τ k the k th jump time, with τ k = + ∞ if there is less than k jumps.

The hybrid state space

The (completed) state space of the model is assumed to have a hybrid structure: E = ∪ q∈Q {q} × E q , where Q is a finite or countable set, and each E q is either the closure of some connected open subset D q ⊂ R nq (n q ≥ 1) or a singleton (in which case we set n q = 0). The state at time t can therefore be written as a pair X t = (Q t , Z t ), where Q t ∈ Q and Z t ∈ E Qt . We denote by Q d = q ∈ Q n q = 0 the set of all "purely discrete" modes, and by E d = ∪ q∈Q d {q} × E q the corresponding subset of E. The usual definitions for smooth maps and vector fields extend without difficulty to such a hybrid structure (see Appendix A for details).

The state space E is regarded as the disjoint sum of the sets E q , q ∈ Q, and endowed with the disjoint union topology 1 . We denote by E the Borel σ-algebra, and by E c the subset of all relatively compact Γ ∈ E. Moreover, we define a "volume measure" on E by the relation

m(Γ) = q ∈Q d m q (Γ ∩ E q ) + x∈E d δ x (Γ) , Γ ∈ E, (1)
where m q is the n q -dimensional Lebesgue measure on E q and δ x the Dirac mass at x. (Note that E q ⊂ R nq has been tacitly identified with {q} × E q ⊂ E.) Let ∂E q be the boundary of E q in R nq , with the convention that ∂E q = ∅ when n q = 0. We define the boundary ∂E of the state space by the relation ∂E = ∪ q∈Q {q} × ∂E q , and the guard set by

G = E \ E 0 . It is not required that G = ∂E.
Notations. Let µ : E → R be a (signed) measure, K : E × E → R a kernel and ϕ : E → R a measurable function. The following notations will be used throughout the paper, assuming the integrals exist: (µK)(dy) = µ(dx) K(x, dy), (Kϕ)(x) = K(x, dy) ϕ(y) and µϕ = µ(dx) ϕ(x).

Stochastic differential equation with jumps

The process X is assumed to be driven by an Itô stochastic differential equation between its jumps: there exist r + 1 smooth vector fields f l and a r-dimensional Wiener process B such that, in mode q ∈ Q \ Q d ,

dZ t = f 0 (q, Z t ) dt + r l=1 f l (q, Z t ) dB l t . (2)
In other words, for all ϕ ∈ C2 (E), X satisfies the following generalized Itô formula

ϕ(X t ) -ϕ(X 0 ) = t 0 (Lϕ)(X s ) ds + r l=1 t 0 (f l ϕ)(X s ) dB l s + 0<τ k ≤t ϕ(X τ k ) -ϕ(X - τ k ) , (3)
where L is the differential generator associated with (2), i.e.

L = i f i 0 ∂ ∂z i + 1 2 i,j r l=1 f i l f j l ∂ 2 ∂z i ∂z j . ( 4
)
We make the following smoothness assumptions:

Assumption 2.
The drift f 0 is of class C 1 , and the other vector fields f l , 1 ≤ l ≤ r, are of class C 2 .

Remark 3. It would be possible to slightly generalize the model by considering a mode-dependent number r q of Wiener processes. All the results of the paper would still hold with the same proofs. We choose to use the same number r of Wiener processes in each mode for the sake of notational simplicity. Note that this is consistent with the most recent definitions of GSHS [START_REF] Bujorianu | Toward a general theory of stochastic hybrid systems[END_REF][START_REF] Everdij | Hybrid Petri nets with diffusion that have into-mappings with generalised stochastic hybrid processes[END_REF], but not with the one in [START_REF] Bujorianu | General stochastic hybrid systems: modelling and optimal control[END_REF] (which uses a mode-dependent number of noise processes).

Two different kinds of jumps

We assume that there exists a Markov kernel K from E to E 0 and a measurable locally bounded function λ : E 0 → R + , such that the following Lévy system identity holds for all x ∈ E 0 , t ≥ 0, and for all measurable ϕ : E×E 0 → R + :

E x 0<τ k ≤t ϕ(X - τ k , X τ k ) = E x t 0 (Kϕ)(X - s ) dH s (5)
where (Kϕ)(y) = E 0 K(y, dy ′ ) ϕ(y, y ′ ) and H is the predictable increasing process defined by

H t = t 0 λ(X s ) ds + τ k ≤t 1 X - τ k ∈G . (6)
The first part corresponds to spontaneous jumps, triggered "randomly in time" with a stochastic intensity λ(X t ), while the other part corresponds to forced jumps, triggered when X hits the guard set G.

Remark 4. The terms "spontaneous" and "forced" seem to have been coined by Bujorianu et al. [START_REF] Bujorianu | A stochastic hybrid system modeling framework[END_REF]. They are closely related to the probabilistic notions of predictability and total inaccessibility for stopping times [see, e.g., 32, chapter VI, § §12-18], but we shall not discuss this point further in this paper.

Remark 5. The pair (K, H) is a Lévy system for the process X in the sense of Walsh and Weil [35, definition 6.1].

Most authors require that H be continuous in the definition of a Lévy system, thereby disallowing predictable jumps.

Mean jump intensity

From now on, we assume that some initial probability law µ 0 has been chosen, with µ 0 (G) = 0 since the process cannot start from G. All expectations will be taken, without further mention, with respect to the probability P µ0 = µ 0 (dx)P x . It is assumed that E(N t ) < + ∞. This is a usual requirement for stochastic hybrid processes 2 , which is clearly stronger than piecewise-continuity of the samplepaths. Its being satisfied depends not only on the dynamics of the system but also on the initial probability law µ 0 .

Definition and connection with the usual stochastic intensity

In order to introduce the main concept of this section, let us define a positive measure R on E × (0; +∞) by

R (A) = E k≥1 1 A X - τ k , τ k . (7)
For any Γ ∈ E, the quantity R (Γ × (0; t]) is the expected number of jumps starting from Γ during the time interval (0; t]. The measure R is in general unbounded, but its restriction to E × (0; t] is bounded for all t ≥ 0 because E(N t ) < + ∞. Definition 6. Suppose that there exists a mapping r : t → r t , from [0; +∞) to the set of all positive bounded measures on E, such that, for all Γ ∈ E,

a) t → r t (Γ) is measurable, b) for all t ≥ 0, R (Γ × (0; t]) = t 0 r s (Γ) ds.
Then r is called the mean jump intensity of the process X (started with the initial law µ 0 ).

Let us split R into the sum of two measures R 0 and R G , corresponding respectively to the spontaneous and forced jumps of the process. Then, using the Lévy system identity, it is easy to see that a mean jump intensity r 0 always exist for the spontaneous part R 0 : it is given by

r 0 t (Γ) = E λ(X t ) 1 Xt∈Γ = Γ λ(x) µ t (dx) . (8)
In other words: for spontaneous jumps, a mean jump intensity always exists, and it is the expectation of the stochastic jump intensity λ(X t ) on the event {X t ∈ Γ}. Forced jumps are more problematic. The Lévy system identity is powerless here, since no stochastic intensity exists (because forced jumps are predictable). All hope is not lost, though: a simple example will be presented in the next subsection, proving that a mean jump intensity can exist anyway. This is fortunate, since the existence of a mean jump intensity will be an essential ingredient for our unified formulation of the generalized FPK equation. See subsection 6.2 for further details on that issue.

Where µ 0 comes into play: an illustrative example

Consider the following hybrid dynamics on E = [0; 1]: the state X t moves to the right at constant speed v > 0 as long as it is in E 0 = [0; 1), and jumps instantaneously to 0 as soon as it hits the guard G = {1} (i.e., the reset kernel is such that

K(1, • ) = δ 0 ).
If we take µ 0 = δ 0 for the initial law, then the process jumps from

1 to 0 each time t is a multiple of 1/v, i.e. τ k = k/v and X - τ k = 1 almost surely.
There is therefore no mean jump intensity in this case, since R = k≥1 δ (1, k/v) . Now take µ 0 to be the uniform probability on [0; 1] (which is, incidentally, the only stationary probability law of the process). Then

R Γ × (0; t] = δ 1 (Γ) 1 0 argmax k≥1 k -x v ≤ t dx (9)
= δ 1 (Γ)

1 0 ⌈vt + x⌉ dx (10) = vt δ 1 (Γ) , (11)
where ⌈vt + x⌉ is the smaller integer greater or equal to vt + x. Therefore the mean jump intensity exists in this case, and is equal to v δ 1 (it is of course time-independent, since µ 0 is stationary). In particular, the global mean jump intensity is r t (E) = v.

Generalized FPK equation

A weak form of the FPK equation

Taking expectations in (3), the following generalized Dynkin formula is obtained: for all compactly supported ϕ ∈ C 2 (E) and all t ≥ 0,

E {ϕ(X t ) -ϕ(X 0 )} = E t 0 (Lϕ)(X s ) ds + E 0<τ k ≤t ϕ(X τ k ) -ϕ(X - τ k ) . ( 12
)
Let us assume the existence of a mean jump intensity r t at all times. Then ( 12) can be rewritten as

(µ t -µ 0 ) ϕ = t 0 µ s (Lϕ) ds + t 0 r s (K -I)ϕ ds , (13)
where µ t is the law of X t and I is the "identity kernel" on E, i.e. the kernel defined by I(y, dy ′ ) = δ y (dy ′ ). Formally differentiating (13) yields

µ ′ t = L * µ t + r t (K -I) , (14)
where t → µ ′ t is the time derivative of t → µ t (in a sense to be specified later), and L * the "distributional adjoint" of L, defined over the set M c (E) of all signed Radon measures3 on E by

(L * ν) (ϕ) = ν (Lϕ) = E (Lϕ)(x) ν(dx) , ∀ν ∈ M c (E), ∀ϕ ∈ C 2 c (E) . (15)
As a consequence of Assumption 2, the result L * ν of applying L * to a Radon measure ν is, in general, a second-order distribution. It is important to note that, because the state space E has a boundary ∂E, the operator L * is not a simple second-order partial differential operator -it also includes "boundary terms". Equation ( 14) begins like the usual Fokker-Planck equation for diffusion processes (µ ′ t = L * µ t ) and ends with an additional term that accounts for the jumps of the process. Definition 7. We will say that t → µ t is a solution in the weak sense of the generalized FPK equation for the GSHS if a) there exists a mean jump intensity t → r t , b) there exists a mapping t → µ ′ t , from [0; +∞) to M c (E), such that t → µ t (Γ) is absolutely continuous with a.e.-derivative t → µ ′ t (Γ), for all Γ ∈ E c , c) L * µ t is a Radon measure for all t ≥ 0, d) equation ( 14) holds as an equality between Radon measures, i.e.

µ ′ t (Γ) = (L * µ t )(Γ) + r t (K -I)(Γ) for all t ≥ 0 and all Γ ∈ E c .
Such a weak form of the FPK equation is the price to pay for a unified treatment of both kind of jumps. Conditions 7.a and 7.b can be seen as smoothness requirements with respect to the time variable, and 7.c with respect to the space variables.

"Physical" interpretation

The usual FPK equation admits a well-known physical interpretation as a conservation equation for the "probability mass" [see, e.g., 17]. Indeed, assuming the existence of a smooth pdf p ∈ C 2,1 (E × R + ), the equation µ ′ t = L * µ t can be rewritten as a conservation equation ∂p t /∂t + div(j t ) = 0, with the probability current j t defined by

j i t = f i 0 p t - 1 2 j ∂(a ij p t ) ∂z j , a ij = r l=1 f i l f j l . (16)
The additional "jump term", in the generalized FPK equation, admits a nice physical interpretation as well. To see this, let us rewrite it as the difference of two bounded positive measure: r t (K -I) = r src t -r t , where r src t = r t K. Therefore r t and r src t behave respectively as a sink and a source in the generalized FPK equation: for each Γ ∈ E, r t (Γ) dt is the probability mass leaving the set Γ during dt, because of the jumps of the process, while r src t (Γ) dt is the probability mass entering Γ.

These two measures are in fact connected by the reset kernel K(x, dy). In particular, the relation r t (E) = r src t (E) holds at all times t ≥ 0, ensuring that the total probability mass is conserved. Moreover, introducing the measures W t (dx, dy) = r t (dx)K(x, dy), we have r t = W (•, dx), r src t = W (dx, •) and the generalized FPK equation can be rewritten more symmetrically as

µ ′ t = L * µ t + (W t (dx, •) -W t (•, dx)) . ( 17
)
It appears clearly, under this form, as a generalization of the differential Chapman-Kolmogorov formula of Gardiner [17, equation 3.4.22] -which only allows spontaneous jumps.

Sufficient conditions for the existence of a weak solution

The main result of this paper shows that the various requirements of definition 7 are not independent. We denote by |ν| the total variation measure of a Radon measure ν, which is finite on E c . We shall say that a function t → ν t from [0; ∞) to M c (E) is right-continuous (resp. locally integrable) if t → ν t ϕ is right-continuous (resp. locally integrable) for all bounded measurable ϕ : E → R.

Theorem 8. Consider the following assumptions:

a) there exists a mean jump intensity r,

such that t → r t is right-continuous, b) t → µ t is differentiable in the sense of 7.b, t → µ ′ t is right-continuous and t → |µ ′ t | locally integrable, c) L * µ t is a Radon measure for all t ≥ 0, t → L * µ t is right-continuous and t → |L * µ t | is locally integrable.
If any two of these assumptions hold, then the third holds as well and t → µ t is a solution in the weak sense of the generalized FPK equation.

The proof of this theorem is given in Appendix B. We will not try to give general conditions under which assumptions 8.a-8.c are satisfied, since such conditions would inevitably be, in the general setting of this paper, very complicated (involving the initial law µ 0 , the vector fields g of the stochastic differential equation, the geometry of the state space E and the reset kernel K).

The case when a piecewise smooth pdf exists

Equation ( 14) is an evolution equation for the measure-valued function t → µ t . In many situations of practical interest, the measures µ t admit a pdf p t , with respect to the volume measure m on E. In this section we show that, if the function p : (x, t) → p t (x) is -at least piecewise -smooth, then equation ( 14) simultaneously gives birth to an evolution equation for t → p t and to static relations that hold for all t ≥ 0 (the so-called "boundary conditions", although the name is not entirely appropriate here).

Assumptions about the guard and the boundary

Turning equation ( 14) into an evolution equation for the pdf ultimately boils down to playing with "integration by parts" formulas, for judiciously chosen test functions. To do so, we shall need additional assumptions concerning the topological regularity of the guard set and the smoothness of the boundary.

Assumption 9.

The guard set G is a regular closed subset of ∂E (i.e., G is a closed set and it is equal to the closure of its interior in ∂E). Assumption 10. For each q ∈ Q such that n q ≥ 2, the domain E q is C 2 -manifold with corners.

See Lee [26, chapter 14] for basic definitions and results concerning manifolds with corners. Assumption 10 is sufficient for the divergence theorem to hold (see Appendix A for a precise statement). The divergence theorem is a multi-dimensional generalization of the "integration by parts" formula, and will be the key tool to compute L * ν for Radon measures with a smooth density.

We denote by s q the surface measure on ∂E q , and define the surface measure s on ∂E by

s = q∈Q nq≥2 s q + q∈Q nq=1 x∈∂Eq δ x . (18)
We further denote by n the outward-pointing unit normal vector on ∂E, which is well-defined s-almost everywhere on ∂E. Since the process X is allowed to start on ∂E \ G, which is a subset of E 0 (see Section 2), the vector fields have to satisfy the following conditions (on the smooth part of ∂E \ G, hence s-almost everywhere):

f 0 , n ≤ 0 , and f l , n = 0 , 1 ≤ l ≤ r . (19)
Otherwise, for any (q, x) ∈ ∂E \ G, the solution of equation ( 2) would leave the domain "instantaneously" (i.e. almost surely in any time neighborhood of 0).

Connecting the mean intensity of forced jumps with the probability current (local result)

Let G 0 denote the subset of the guard set G where at least one of the "noise" vector fields is not tangent to the boundary, i.e. G 0 = {x ∈ G, ∃l ∈ {1, . . . , r}, f l , n = 0}. The following results relates the mean intensity of forced jumps with the probability current j t defined by equation [START_REF] Feller | Diffusion processes in one dimension[END_REF].

Proposition 11. Assume that the measures µ t admit a pdf p t = p( • , t) for all t ≥ 0 on some open subset U ⊂ E, with p ∈ C 2,1 (U × R + ). Define the outward probability current j out t = j t , n on U ∩ G. Then, for all t ≥ 0,

a) j out t ≥ 0 and r G t (Γ) = Γ j out t ds is the mean intensity of forced jumps on U ∩ G, b) the pdf p t vanishes on U ∩ G 0 .
See Appendix C for the proof. This proposition provides two important conclusions concerning forced jumps. The first one is that, when a smooth pdf exists in a neighborhood of the guard set, the mean intensity of forced jumps (which appears in the FPK equation) is equal to the outward flow of the probability current. This is consistent with the physical interpretation of the probability current: j t ds dt is the probability mass ("number of particles") escaping from the domain through ds during dt.

The second conclusion is that the familiar "absorbing boundary" condition p t = 0 holds on the guard set as soon as one of the "noise" vector fields is active in the normal direction. Note that the pdf does not vanish on the boundary in the example of subsection 3.2, which is a piecewise deterministic process with forced jumps. A "physical" explanation of absorbing boundaries, in the spirit of subsection 4.2, can be found in [START_REF] Malhamé | Electric load model synthesis by diffusion approximation of a high-order hybrid-state stochastic system[END_REF] and also, more recently, in [START_REF] Lubashevsky | Boundary singularities and boundary conditions for the Fokker-Planck equations[END_REF].

Evolution equation for the pdf and "boundary" conditions (global result)

The local result of subsection 5.2 will now be used to obtain a general formulation of the FPK equation ( 14) in terms of a probability density function, when one exists and is smooth enough. Let H ⊂ E 0 \ E d be a closed set of m-measure zero -typically, H will be a closed hypersurface in applications. Note that U = E \ H is an open neighborhood of the boundary ∂E. Assume now that the following holds: Assumption 12.

a) µ t admits a pdf p t with respect to m, on the whole state space, for all t ≥ 0, b) p ∈ C 2,1 (U × R + ), with ∂p ∂t and F p locally integrable on E × R + .

Then, it follows from the proof of Proposition 11 (see Appendix C, equation (48)) that

(L * µ t ) (Γ) = Γ F p t dm + ∂E∩Γ j out ds , ∀Γ ∈ E c such that Γ ⊂ U , ( 20
)
where F is the formal adjoint of L, i.e., the differential operator defined by

F : q → - i ∂ f i 0 q ∂z i + 1 2 i,j ∂ 2 a ij q ∂z i ∂z j . (21)
The (possible) lack of differentiability of p t on H therefore translates into the fact that the Radon measures β t ,

β t (Γ) = Γ F p t dm + ∂E∩Γ j out ds -(L * µ t ) (Γ) . ( 22
)
do not vanish in general. This, in turn, is closely related to the existence of a non-vanishing m-singular part (see Appendix D for a definition) in the source term r src t = r t K, as stated by the following result.

Theorem 13. Let Assumption 12 hold. Then the conditions 8.a-8.c of Theorem 8 are satisfied, and the following evolution equation holds on E 0 \ H, for all t ≥ 0: Finally, the m-singular part (r t K) ⊥ of r t K is supported by the set H ∪ (∂E \ G) and satisfies the following "conservation equations":

∂p ∂t = F p t + d(r t K) dm -λ p t . (23
c) (r t K) ⊥ = β t ≥ 0 on H, d) (r t K) ⊥ = - • j out t ds ≥ 0 on ∂E \ G.
See Appendix E for the proof.

Examples

A class of models with spontaneous jumps

Our first series of examples covers a large family of models without forced jumps (G = ∅). The reset kernel K is assumed to satisfy the following assumption: Assumption 14. There exists a kernel K * on E such that m(dx) K(x, dy) = m(dy) K * (y, dx) .

(

and the following evolution equation holds:

∂p t ∂t = L * p t + K * (λ p t ) -λ p t . (26)
Assumption 14 holds for several classes of models known in the literature: pure jump processes with an absolutely continuous reset kernel, the switching diffusions of Ghosh et al. [START_REF] Ghosh | Ergodic control of switching diffusions[END_REF][START_REF] Ghosh | Optimal control of switching diffusions with application to flexible manufacturing systems[END_REF] and also the SHS of Hespanha [START_REF] Hespanha | A model for stochastic hybrid systems with application to communication networks[END_REF].

Example 16. Pure jump processes occur when L = 0, i.e. when there is no continuous dynamics. We consider here the case where K is absolutely continuous: K(x, dy) = k(x, y) m(dy). For instance, if the amplitude of the jumps is independent of the pre-jump state and distributed according to the pdf ρ, then k(x, y) = ρ(y -x). In this case Assumption 14 holds with K * (x, dy) = k(y, x) m(dy). Introducing the function γ(x, y) = λ(x)k(x, y), equation 26 turns into the well-known master equation [17, eq. 3.5.2]:

∂p ∂t (y, t) = γ(x, y)p(x, t) -γ(y, x)p(y, t) m(dx) . ( 27
)
In particular, when all modes are purely discrete (n q = 0), this is just the usual forward Kolmogorov equation for a continuous-time Markov chain.

Example 17. In the case of switching diffusions, the state space is of the form E = Q × R n (with Q a countable set and n ≥ 1) and the reset kernel of the form

K (q, z), • = q ′ =q π qq ′ (z) δ (q ′ ,z) , (28)
where π(z) = (π qq ′ (z)) is a stochastic matrix for all z ∈ R n . Assumption 14 is fulfilled with K * defined by

K * (q, z), • = q ′ =q π q ′ q (z) δ (q ′ ,z) . ( 29
)
Equation 15 becomes in this case the familiar generalized FPK equation for switching diffusion processes [see, e.g., [START_REF] Kontorovich | Dynamic systems with random structure: an approach to the generation of nonstationary stochastic processes[END_REF][START_REF] Krystul | Risk decomposition and assessment methods[END_REF]: for all x = (q, z) ∈ E and t ≥ 0,

∂p ∂t (x, t) = (L * p t )(x) + q ′ =q λ q ′ q (z) p t (q ′ , z) -λ(x) p t (x) , (30)
where λ q ′ q (z) = λ(q ′ , z) π q ′ q (z).

Example 18. The SHS of Hespanha [START_REF] Hespanha | A model for stochastic hybrid systems with application to communication networks[END_REF] are also defined on E = Q × R n , but this time the post-jump state X τ k is determined by applying a reset map Ψ : E → E 0 to the pre-jump state X - τ k , Ψ being chosen randomly in a finite of reset maps Ψ k . The reset kernel can therefore be written as

K(x, •) = k π k (x) δ Ψ k (x) , (31)
with π k (x) the probability of choosing the reset map Ψ k given that X - τ k = x. Provided that the functions Ψ k are local C 1 -diffeomorphisms, the kernel K fulfills Assumption 14 with

K * (x, •) = k y∈Ψ -1 k ({x}) π k (y) J k (y) -1 δ y , (32)
where J k (y) is the Jacobian determinant of Ψ k at y. Therefore, introducing a stochastic intensity λ k = λ ̺ k for each one of the reset maps, we recover thanks to Corollary 15 the generalized FPK equation given by Hespanha [20, p. 1364]:

∂p ∂t (x, t) = (L * p t )(x) + k y∈Ψ -1 k ({x}) λ k p t |J k | (y) -(λ k p t )(x) . ( 33
)
H ∂E Ψ E 0 Figure 1:
State space for the model of subsection 6.2.

A class of models with forced jumps

The measure-valued formulation of the generalized FPK equation ( 14) paves the way for an easier proof of some recent results [START_REF] Bect | Generalized Fokker-Planck equation for piecewise-diffusion processes with boundary hitting resets[END_REF], concerning GSHS with forced jumps and deterministic resets. A typical example of this class of process is the thermostat model of Malhamé and Chong [START_REF] Malhamé | Electric load model synthesis by diffusion approximation of a high-order hybrid-state stochastic system[END_REF], which has been extended to several dimensions in [START_REF] Bect | Generalized Fokker-Planck equation for piecewise-diffusion processes with boundary hitting resets[END_REF]. We consider the class of GSHS models satisfying the following assumptions.

Assumption 19.

a) The model only has forced jumps (λ = 0) with deterministic resets, i.e. there exists a map

Ψ : G → E 0 such that K(x, • ) = δ Ψ(x) for all x ∈ G.
b) All modes have the same dimension n q = n, the guard set is the whole boundary (G = ∂E) and is of class C 1

(in particular, it has no corners).

c) H = Ψ(G) is a C 1 hypersurface, closed in E, and Ψ is a C 1 -diffeomorphism from G to H.
The state space of this model is depicted on figure 1. The assumption that G = ∂E is only here for the sake of simplicity and could easily be relaxed. Boundaries with corners and piecewise smooth reset maps could be considered as well. The model considered in [START_REF] Bect | Generalized Fokker-Planck equation for piecewise-diffusion processes with boundary hitting resets[END_REF] also includes purely discrete modes (i.e., n q ∈ {0, n} for each q ∈ Q), which cause no real additional difficulty. The measure K(x, • ) is supported by H for all x, which implies that the source term r t K is also supported by H, hence is m-singular. Therefore, even if the diffusion is non-degenerate (i.e. the diffusion matrix (a ij ) is uniformly positive definite), we know from subsection 5.3 that the pdf p t will not be smooth on H. Accordingly, we make the following smoothness assumption for the measures µ t : Assumption 20.

a) µ t admits a pdf p t with respect to m on the whole state space, for all t ≥ 0, b) p ∈ C 2,1 ((E \ H) × R + ), p t and ∇p t have at most a jump discontinuity (discontinuity of the first kind) on H.

Then Assumption 12 holds, which allows Theorem 13 to be applied. Moreover, the result of subsection C.1 holds on each component

C of E \ H -i.e., for all ϕ ∈ C 2 c (E), (L * µ t ) (ϕ |C ) = C ϕ F p t dm + ∂C ϕ j t , n ∂C ds + 1 2 r l=1 ∂C f l ϕ p t f l , n ∂C ds ,
where n ∂C is the outward-pointing unit vector on ∂C. Summing over the components and using Proposition 11.b yields

(L * µ t ) ϕ = E F p t ϕ dm + G j out t ϕ ds - H j in t ϕ ds - 1 2 r l=1 H f l ϕ f l , n ab p b t -p a t ds H , (34)
where n ab is the unit normal vector on H oriented from side a to side b and

j in t = j (b) t -j (a)
t , n ab . The superscripts a/b indicate the value of a discontinuous function on the corresponding side of H (but none of these quantities actually depend on the chosen labelling of the sides). The last term on the right-hand side vanishes, because L * µ t is a Radon measure by Theorem 8, whereas this term only involves the first-order derivatives of ϕ (through f l ϕ). Therefore, p t is continuous on the set H 0 = {x ∈ H, ∃l ∈ {1, . . . , r}, f l , n ab = 0} and the measures β t of Theorem 13 are given by

β t (Γ) = H∩Γ j in t ds H . (35)
The conclusions of Theorem 13 can then be summarized as follows:

a) The usual Fokker-Planck equation, ∂p t /∂t = F p t , holds on each component of E 0 \ H.

b) The conservation probability current through the reset map is ensured by the relation

j out t = d(s H • Ψ) ds j in t • Ψ , (36)
where

s H • Ψ is the pushforward (image measure) of s H by Ψ -1 .
c) The absorbing boundary condition, p t = 0, holds on G 0 .

d) The density p t is continuous on H 0 .

Conclusions

A general measure-theoretic formulation of the Fokker-Planck-Kolmogorov equation has been presented, in the modern framework of GSHSs. This formulation is new, and should allow GSHS practitioners to get a better understanding of how a given system behaves in terms of evolution of the probability mass in the state space.

Technical tools have been provided, in order to derive the explicit form of the evolution equation when a probability density function exists and satisfy sufficient regularity conditions. In particular, it has been shown that the general FPK equation allows to recover all previously known instances of the FPK equation for stochastic hybrid systems.

Of course, an important issue is now to provide sufficient conditions for the existence of a "smooth enough" probability density function and the uniqueness of the solution to the generalized FPK equation. The literature already provides such conditions for processes defined by stochastic differential equations and (in some cases) switching diffusions; see, e.g., [START_REF] Bichteler | Malliavin Calculus for Processes with Jumps[END_REF][START_REF] Ikeda | Stochastic differential equations and diffusion processes[END_REF][START_REF] Bris | Existence and uniqueness of solutions to Fokker-Planck type equations with irregular coefficients[END_REF] and the references therein. Extending these results to other types of GSHSs, including models with forced jumps, is an important perspective for future work.

In addition to providing a better understanding of GSHSs, the FPK equation is also a powerful tool for the analysis of low-dimensional systems, for which an approximate solution can be obtained using numerical methods (for instance finite volume methods). It is especially useful for the computation of the stationary distribution, as shown in [START_REF] Bect | On the Fokker-Planck equation for stochastic hybrid systems: application to a wind turbine model[END_REF] a nontrivial three-dimensional model of a wind turbine. This type of application of the FPK equation to the analysis of GSHSs relies on the availability of software components allowing an easy implementation of efficient numerical methods, in the spirit of Mitchell's Level Set Toolbox [START_REF] Mitchell | The flexible, extensible and efficient toolbox of level set methods[END_REF][START_REF] Mitchell | A toolbox of Hamilton-Jacobi solvers for analysis of nondeterministic continuous and hybrid systems[END_REF] for Hamilton-Jacobi equations. The development of such a toolbox is another important direction for future work.

A. Smooth maps and vector fields on E

The following definitions are natural extensions to the hybrid state space E of the usual definition on subsets of R n .

A map ϕ :

E → R is said to be k-times continuously differentiable on E -in short, ϕ ∈ C k (E) -if ϕ q = ϕ(q, • ) is C k on E q in the usual sense for all q ∈ Q \ Q d , i.e. if
there is an open subset U of R nq such that E q ⊂ U and ϕ extends to a C k map on U .

A vector field g on E is defined as a first-order differential operator with respect to the continuous variables. Its action on a continuously differentiable function ϕ ∈ C 1 (E) will be denoted by gϕ, where

(gϕ)(q, z) = nq i=1 g i (q, z) dϕ dz i (q, z) on E \ E d , 0 on E d . ( 37
)
The number of "components" of g depends on the mode q. To simplify the notations, we shall agree that the indices i and j always correspond to summations on the number of continuous variables, and drop the explicit dependence on q.

For instance, the definition of gϕ on E \ E d can be rewritten as gϕ = i g i ∂ϕ ∂z i . A vector field is said to be k-times continuously differentiable on E if g i (q, • ) is C k on E q in the usual sense for all q ∈ Q \ Q d and all i ∈ {1, . . . , n q }. Finally, under Assumption 10, the following version of the divergence theorem holds (since each component of the state space is a C 2 -manifold with corners) : Theorem 21 (see, e.g., [START_REF] Lee | Introduction to smooth manifolds[END_REF]). Let Assumption 10 hold. Then, for all compactly supported C 1 vector field f on E,

E div(f ) dm = ∂E f , n ds . (38)

B. Proof of Theorem 8

Let C 2 c (E) denote the set of all compactly supported ϕ ∈ C 2 (E). The following lemma is an easy consequence of the smoothness of the vector fields:

Lemma 22. For all ϕ ∈ C 2 (E), t → t 0 (L * µ s )(ϕ)
ds is differentiable on the right, with the right-continuous derivative t → (L * µ t )(ϕ).

In the sequel, "right-continuous" is abbreviated as "rc". ⋄ Assume that both 8.a and 8.b hold. Then each term of ( 13) has a t-derivative on the right. Differentiating both sides proves that [START_REF] Everdij | Hybrid Petri nets with diffusion that have into-mappings with generalised stochastic hybrid processes[END_REF]

µ ′ t ϕ = (µ t -µ 0 )ϕ , ∀t ≥ 0 , ∀ϕ ∈ C 2 c (E) . (40)
Moreover, for all Γ ∈ E c ,

t 0 |µ ′ s | (Γ) ds ≤ t 0 |L * µ s | (Γ) ds + 2 E N t ≤ +∞ , (41)
which shows that t → |µ ′ s | is locally integrable. Therefore, using standard approximation techniques and a monotone class argument, it can be proved that (40) still holds for ϕ = 1 Γ , Γ ∈ E c , i.e. that t → µ ′ t is the "derivative" of t → µ t in the sense of definition 7.b.

⋄ Finally, assume that 8.b and 8.c hold. Then, for all ϕ ∈ C 2 c (E), equation ( 13) can be rewritten as

G×]0;t] ϕ(x) R G (dx, ds) -(L * µ s )(dx)ds = E 0 ×]0;t] ϕ(x) (R G K)(dx, ds) -ξ s (dx) ds , (42)
where

ξ s = µ ′ s -L * µ s (E 0 ∩ • ) -r 0 (K-I).
The measures R G and r 0 have been defined in subsection 3.1. Clearly, ξ t ∈ M c (E) and t → ξ t is locally integrable. Using once more standard approximation techniques, one can prove that (42) still holds when ϕ = 1 Γ , with Γ a compact subset of G. In this case the right-hand side vanishes, yielding

R G (Γ×]0; t]) = t 0 (L * µ s )(Γ) ds . (43)
Moreover, since t → R G (Γ×]0; t]) is increasing and t → (L * µ t )(Γ) is rc, we have (L * µ t )(Γ) ≥ 0 for all t ≥ 0. This allows to extend (42) to all Γ ∈ E c , using a monotone class argument, thus proving the existence of a mean jump intensity r G t = (L * µ s )(G ∩ • ) for the forced jumps.

C. Proof of Proposition 11

Since p is of class C 4)) and j t (equation ( 16)), we find that

Lϕ p t =   i f i 0 ∂ϕ ∂z i + 1 2 i,j a ij ∂ 2 ϕ ∂z i ∂z j   p t = j t ϕ + 1 2 i,j ∂(a ij p t ) ∂z j ∂ϕ ∂z i + i,j a ij p t ∂ 2 ϕ ∂z i ∂z j = j t ϕ + 1 2 i,j ∂ ∂z j a ij p t ∂ϕ ∂z i = j t ϕ + 1 2 r l=1 i,j ∂ ∂z j f i l f j l p t ∂ϕ ∂z i = j t ϕ + 1 2 r l=1 div (f l ϕ p t f l ) . (44)
Moreover, using the product rule for the divergence operator and the fact that F p t = -div(j t ), which is a direct consequence of equations ( 16) and ( 21), we get

j t ϕ = div(ϕj t ) -ϕ div(j t ) = div(ϕj t ) + ϕ F p t . (45)
Finally, equations ( 44) and (45) together with the divergence theorem yield :

(L * µ t ) (ϕ) = E Lϕ dµ t = E Lϕ p t dm = E div(ϕj t ) + ϕ F p t + 1 2 r l=1 div (f l ϕ p t f l ) dm = ∂E ϕj t , n ds + E ϕ F p t dm + 1 2 r l=1 ∂E f l ϕ p t f l , n ds = E ϕ F p t dm + ∂E ϕj out t ds + 1 2 r l=1 ∂E f l ϕ p t f l , n ds . (46)

C.2. Proof of assertion 11.b

Let V = (∂E) smooth ∩ U , where (∂E) smooth denote the smooth part of the boundary.

V is an open subset of ∂E. For each η ∈ C 2 c (V ), there exists a sequence of functions ϕ n ∈ C 2 c (U ), with their support in a fixed compact set, such that ϕ n = 0 and ∂ϕ n /∂n = η on V , and ϕ n → 0 uniformly. Equation (46) holds for each n. Taking limits with respect to n on both sides, and using the fact that (L * µ t ) (ϕ n ) → 0 (since L * µ t is a Radon measure on U ), we find that l ∂E η p t f l , n 2 ds = 0. Therefore, x → p(x, t) l f l , n 2 vanishes s-almost everywhere on V , hence everywhere on ∂E by continuity. This proves assertion 11.b since l f l , n 2 > 0 on G 0 .

C.3. Proof of assertion 11.a

It is now proved that, for each x ∈ ∂E, either p t (x) = 0 or f l , n x = 0 for all l ∈ {1, . . . , r}. As a consequence, for all ϕ ∈ C 2 c (U ), the last term of equation (46) vanishes :

(L * µ t ) (ϕ) = E F p t ϕ
Finally, for all t ≥ 0 and all Γ ∈ E c such that Γ ⊂ G∩U , we have µ ′ t (Γ) = 0 and Γ F p t dm = 0, both because m(Γ) = 0, and (r t K)(Γ) = 0 because K is a kernel from E to E 0 . Therefore, as a consequence of (49), r G t (Γ) = Γ j out t ds for all Γ ∈ E c such that Γ ⊂ G ∩ U . The outward current j out t is thus positive s-almost everywhere on U , hence everywhere by continuity. This proves assertion 11.a.

D. The Lebesgue-Radon-Nikodym theorem

This appendix recalls a fundamental result of measure theory, which is used in the statement and proof of Theorem 13. The reader is referred to, e.g., [10, chapter 4], for basic definitions and terminology not recalled here. In this section, (E, E) denotes any measurable space -not necessarily the hybrid state space as defined in subsection 2.1.

Theorem 23 (Lebesgue-Radon-Nikodym). Let ν 1 , ν 2 be positive σ-finite measures on (E, E). Then there is a unique positive σ-finite measure ν ⊥ 1 and a unique (up to ν 2 -everywhere equality) measurable function f such that a) ν ⊥ 1 and ν 2 are mutually singular, i.e. ∃A ∈ E, ν ⊥ 1 (A) = 0 and ν 2 (E \ A) = 0, b) ν 1 has the following decomposition

ν 1 (A) = A f dν 2 + ν ⊥ 1 (A), ∀A ∈ E . (50)
The measures A → A f dν 2 and ν ⊥ 1 are respectively called the absolutely continuous part and the singular part of ν 1 with respect to ν 2 . Equation (50) is called the Lebesgue-Radon-Nikodym decomposition of ν 1 with respect to ν 2 . The function f is called the Radon-Nikodym derivative of ν 1 with respect to ν 2 , is usually denoted by dν1 dν2 .

E. Proof of Theorem 13

The assumption 8.b of Theorem 8 holds with µ ′ (Γ) = Γ ∂p ∂t , since ∂p ∂t exists m-almost everywhere and is locally bounded. Moreover, as in the proof of Proposition 11, the existence of r G t on G ∩ U = G follows from the fact that 8.b and 8.c hold on U . We have thus proved that the assumptions 8.a (existence of r t ) and 8.b of Theorem 8 (existence of µ ′ ) hold on the whole state space, which implies that 8.c and the generalized FPK equation ( 14) hold as well.

As a consequence of equation ( 48 The first line of the system (53)-( 54) is precisely [START_REF] Kontorovich | Dynamic systems with random structure: an approach to the generation of nonstationary stochastic processes[END_REF], and the second one readily splits into 13.c and 13.d by considering the terms that are supported respectively by H and ∂E \ G. The proof of Theorem 13 is thus complete.

)

  Moreover, according to Proposition 11, a) r G t (Γ) = Γ∩G j out t ds is the mean intensity of forced jumps, b) and the absorbing boundary condition, p t = 0, holds on G 0

  If there exists a pdf p ∈ C 2,1 (E × R + ), then the measures r t and r src

					t are absolutely continuous with
	respect to m,	dr t dm	= λ p t ,	dr src t dm	= K

)

(

We do not assume that K * is a Markov kernel, i.e. that K * (y, •) is a probability measure.) The following result is an easy consequence of Theorem 8: Corollary 15. * (λ p t ) ,

  holds for all t ≥ 0, hence that L * µ t is a Radon measure and that t → L * µ t is rc. Moreover, integrating the inequality |L * µ t | ≤ |µ ′ t | + 2r t yields that, for all Γ ∈ E c , Therefore t → |L * µ s | is locally integrable, which proves 8.c. ⋄ Assume now that 8.a and 8.c hold, and set µ ′ t = L * µ t + r t (K -I), for all t ≥ 0. Clearly, µ ′ t is a Radon measure, t → µ ′ t is rc and

	0	t	|L * µ s | (Γ) ds ≤	0	t	|µ ′ s | (Γ) ds + 2 E N t ≤ +∞ .	(39)

t 0

  dm + ∈ E c such that Γ ⊂ U . Substituting the result into the generalized FPK equation (14) yields

	for all Γ µ ′ t (Γ) =
	∂E	ϕj out t	ds ,	(47)
	∂E∩Γ	j out t	ds	(48)

and therefore the Radon measure L * µ t can be rewritten as

(L * µ t ) (Γ) = Γ F p t dm + Γ F p t dm + ∂E∩Γ j out t ds + (r t K)(Γ) -r t (Γ) .

  ∈ E c , since m(E ⊂ U ) = 0 and ∂E ⊂ U . Equation (22) thus simplifies into β t (Γ) = -(L * µ t ) (Γ ∩ H),which proves that the measures β t are supported by H.According to (51), the generalized FPK equation (14) can be decomposed as + (r t K) ⊥ (dx) , r t (dx) = λ(x) p t (x) m(dx) + r G t (dx) .

	where we have used that
				µ ′ t (dx) =	∂p ∂t	(x, t) m(dx) ,
				(r t K)(dx) =	d(r t K) dm	(x) m(dx)
									∂E∩Γ∩U	j out t	ds =	Γ	F p t dm +	∂E∩Γ	j out t	ds	(51)
	for all Γ µ ′ t (Γ) =	Γ	F p t dm +	∂E∩Γ	j out t	ds + (L (52)
	Uniqueness of the Lebesgue-Radon-Nikodym decomposition with respect to m yields
	∂p ∂t	= F p t +	d(r t K) dm	-λ p t ,	(53)
	0 =								(54)

), we have

(L * µ t ) (Γ ∩ U ) = Γ∩U F p t dm + * µ t ) (Γ ∩ H) + r t (K -I)(Γ) . ∂E∩Γ j out t ds + (L * µ t ) (Γ ∩ H) + (r t K) ⊥ (Γ) -r G t (Γ) ,

A shorter version of this paper was presented at the 17th IFAC World Congress (IFAC'08) in Seoul, Korea[START_REF] Bect | A unifying formulation of the Fokker-Planck-Kolmogorov equation for general stochastic hybrid systems[END_REF].66 The results presented in this paper come from the PhD thesis of the author[START_REF] Bect | Processus de Markov diffusif par morceaux: outils analytiques et numériques[END_REF], under the supervision of Pr. Gilles Fleury and Dr. Hana Baili.Email addresses: julien.bect@supelec.fr (Julien Bect)Paper published by Elsevier in Nonlinear Analysis: Hybrid Systems, 4(2):357-370, may 2010. This is an author-generated post-print version.

which is (here) locally compact, separable and completely metrizable

See, e.g., Davis[START_REF] Davis | Piecewise deterministic Markov processes: a general class of nondiffusion stochastic models[END_REF] or Bujorianu and Lygeros[START_REF] Bujorianu | General stochastic hybrid systems: modelling and optimal control[END_REF].

In this paper, a "Radon measure" will always be a signed Radon measure, in other words a distribution of order zero; see Rudin[START_REF] Rudin | Function Analysis[END_REF] for the basic definitions and properties of distributions. Any signed Radon measure ν can be written as the difference ν = ν + -ν -of two positive Radon measures (i.e. locally finite measures); see, e.g., Cohn[10, chapter 7] for more information on the connection between the functional analytic and the measure-theoretic point of view.