Roland Topor

The nocturnal condition Luc Gwiazdzinski

According to the Genesis, "God saw the light, that it was good: and God divided the light from darkness. And God called the light Day, and the darkness he called Night. And the evening and the morning after were the first day." Just as it happens with the human organism, the city has always had its existence defined by the day/ night cadency. As a period of obscurity symbolized by the curfew, by the stoppage of all activities and closing of the city gates, the night has for a very long time been considered a time of social rest and return to the intimacy of the private sphere. As a forgotten dimension of the city, the night inspired poets in their search for freedom, served as refuge for scoundrels and disquieted power-structures, which have set out to control it. Mayors, technicians and researches have for quite some time imagined, lead or managed cities as if they only functioned 16 hour per day, Monday to Saturday, oblivious to the city's nightlife. Still today, the night usually is absent from reflections related to urbanism, development and territorial management, or merely limited to inconvenient aspects or merely functional ones, such as lighting. Only artists, writers or poets have had an interest in the urban night. For centuries the 'other side' of our cities has remained a space-time largely disinvested by human activity, a "Finisterre"1 human ambitions have turned against, a last internal frontier that begins now to instigate the appetites.

But times have changed. There is a life after day in our cities where the non-day sometimes makes up two-thirds of the workday. Fortunately for night-birds, the conquest of the night has in the meantime begun. There is a life after day in our cities where the non-day sometimes makes up two-thirds of the workday. Freed from natural problems, our cities liven up under the influence of evermore desynchronized modes of life, of reduced working days or of new lighting and communication technologies. Insofar as society profoundly redefines its nycthemeron 2 , the city reconfigures itself. Human activities progressively open themselves up to the night and compose a new space of work and leisure forcing us to consider the spaces and times of the city under a new guise. Colonized by typically day-related activities, the night is under new pressures. Not always the city which sleeps (and works) and the city which has fun make a good couple. Jobs involving atypical shifts, specific remuneration for night jobs, electricity outages, night-related commerce, lights, white nights, transport, late opening of retail, pollution (light and noise) or urban violence: between insecurity and freedom, for the better or for the worst, the night is definitely here... And our organizations aren't ready for it.

As a central field of tension in our society, a space of ideas and last-standing frontier of the city, the 'night' has plenty to say to the day'. The night should be open to scientific investigation, to the future, to creativity. The time has arrived to explore the nocturnal archipelago, decode the urban night, discover its actors, its limits, its centers, its margins, its rhythms. In face of pressures resulting from the continuous temporality of the economy and of the system, we should invest in the night, explore the last-standing frontier in the city so as to anticipate conflicts and imagine possible futures of the Urbs and Civitas. We need to place the night in the political agenda and launch a grand debate in favor of more beautiful cities than the ones we live in during the day.

A SPACE-TIME IN MUTATION

In his constant search to free himself for natural rhythms, Man has little by little colonized the night.

A progressive colonization by daytime activities.

Human activities progressively open themselves up to the night and compose a new space of work and leisure. With the passing of the centuries, political affirmation as well technical developments with regard to lighting have allowed for the extension of typically day-related activities into the night and the emergence of a nocturnal public space. Under the pressure of the continuous time of the economy and the system, day-related activities increasingly spread into the night, leading to the becoming-day of nocturnal space as the final phase in the artificialization of the city and the death of the original bipolar alternation between day and night.

Between insecurity and freedom, entrepreneurs take over our nights under the risk of losing themselves and leading us to the same fate: marketers of the future who, however, do not promise us the moon. The limits of the 'day' have been surpassed. Natural rhythms have gradually ceased. Daylight saving time allows us to take advantage of public space until later.

Light has with time taken over urban public space, partly wiping out the threatening darkness of our nights and allowing for the continuity of daytime activities. At least in cities, lamp posts turn on before the stars do. Lighting passes from security to pleasure. Sounds, lights, the luminous decorations of buildings, multiply as they swathe the cities of the world with the same coating of yellow light. 'Light designers' sculpt the night, giving to our cities a nocturnal identity. At night, the urban labyrinth rearranges itself. Another city enters the stage with its lights, its decoration, its new actors. The same and yet another city.

At night, the urban labyrinth rearranges itself. Another city enters the stage with its lights, its decoration, its new actors. Industries function continuously in order to maximize their equipments, and, in the majority of sectors, nocturnal work has become trivial. More and more service-providers adhere to the 7 day a week, 24h per day, publicity catch-line, making it a commonplace. Numerous activities alternate their advantage of nighttime. Stores offer more and more nocturnal events. All over, the supply of night-related leisure expands. In a few years, the night has become a separate economic sector. Thematic nights bring in revenue and the nocturnal calendar thickens with the 'Night of the Museums', the 'European Night of Science' in Berlin and other cities, the 'White Nights' in Paris, Rome, Brussels, Montreal and other places.

All over Europe, legislation on the opening hours for night-related commerce is slackening. In Great Britain, pubs can remain open after 23h. In Germany, metro stations are gradually becoming 'oasis of nocturnal services'. In Poland, many places remain open round the clock. In England, only the opening times of supermarkets are predetermined. In Belgium, a law legalized the famous night shops. In Denmark, shop-owners can now freely choose their working hours during the week. In Spain, there is talk of total liberalization. In Italy, each commune decides its own working hours, while in Portugal shops can stay open until midnight.

Festive nights begin later and later, for the desperation of club-owners. For quite a while now media curfews have ended, with radio and television stations working round the clock. We sleep an hour less than our grandparents and we go to bed later.

A nocturnal archipelago. In the majority of cities, public transportation stops, taxis circulate with more expensive rates and in lesser quantity, the waiting times (while they still run) of buses is steeper, making the cost of accessing urban space necessarily rise. As a consequence of these pressures, the night, defined as an interruption of these activities, is a space-time amounting to three hours: between 1h30 and 4h30. Vanguards also progress in space discontinuously: sectors dedicated to nocturnal leisure have sprung in the hearts of old city centers as well as in competing peripheral areas which have gradually included multiplex businesses and clubs. Automatic points of distribution (for tickets, drinks, etc.) are installed in every corner while conduits of international flux (streets, railroads, airports) cut across metropolitan centers, feeding nocturnal nomads with oasis of permanent services (service stations, retail in train stations or airports...). At the night's core, the city transforms itself into a fortress, accessible only to those on foot or possessing private means of transport. In the majority of cities, public transportation stops, taxis circulate with more expensive rates and in lesser quantity, the waiting times (while they still run) of buses is steeper, making the cost of accessing urban space necessarily rise. It is the image of the archipelago which comes to mind, imposing itself while imagining the geography of urban nightlife insecurity and limited freedom. If the night isn't as dangerous as it is generally said to be, neither is it the space of freedom and of encounter dreamt by the poets. The range of services offered at night reduces. The majority of public spaces close. The limitation or simple absence of public means of transport is a barrier. The night is expensive and free collective spaces are scarce. There is a tendency to segregate and box the publics according to age, sex, taste. Each, therefore, in his own bar, club, and that is ok, what the night is meant to be. Light and media manipulate us, and just as if we were moths, we flee the stigmatized peripheries to dry out our wings in the neon lights of cities.

PRESSURES WHICH HEIGHTEN

The pressures which heighten upon the night crystallize fundamental economic, political and social interests.

A critical field of tensions. In the backdrop, dominant classes perk up to suppress the night. Colonized by light and day-related activities, traversed by users living increasingly alternative rhythms, the night has transformed itself into a critical field of tensions. The continuous time of the economy and the system is opposed to the cyclical rhythm of our bodies and cities. World time clashes with local time. Conflicts multiply between individuals, groups and neighborhoods of the city which sleeps, works, and the ones of the city which enjoys itself. There are rebellions against light-pollution, the fact that it has killed the magic of our nights, depriving us of the stars as a free spectacle. In the city-centers, conflict emerges between residents worried about tranquility and people frequenting bars and clubs, themselves symbols of the surging of a nocturnal public space. Prostitution prospers. In peripheral neighborhoods, the burning of vehicles occurs between 22h and 1h, right when the natural functioning of urban societies (commerce, social centers, police, etc.) has disappeared and the streets are taken over by male adolescents. In many European agglomerations, let us also point out that there is opposition by people living close to airports to the implementation of a 24h public system of transport.

Alarming conflicts. Notwithstanding the usual clichés, it is true that urban nightlife questions our capacity of living together. Permanent paradox, the city, place of reunion, of differentiation, is by its very nature the theatre of conflicts between neighborhoods, activities, groups and individuals of a 'polychronic' society. During the night, these conflicts, more visible than during the day, allow us to understand the night as a frontline, and this in the sense of demarcating the line delimiting the spatial and temporal conquests of typically day-related activities. We flee more at night. Dark or lit, alive or aseptic, dangerous or policed, centralized or not, warm or repulsive, agitated or calm, the night has invited itself into the party we have usually thought as pertaining to the day. These evolutions instigate us to change our approach to these agglomerations, to cities, and to do so in terms of social time, rhythms. They oblige us to occupy and explore urban nightlife.

Investigations which are needed. Researchers, public agencies and citizens, should occupy this space-time so as to anticipate conflicts between individuals, groups or neighborhoods, and together imagine the contours of a new urban life. It is time to anticipate the predictable developments of nocturnal activities in order to think a global, 24h, management of the city. Researchers, public agencies and citizens, should occupy this space-time so as to anticipate conflicts between individuals, groups or neighborhoods, and together imagine the contours of a new urban life. It is specifically along these spatial and temporal margins, along the peripheries, that the future of our society stages itself. It is necessary to populate and liven up the night while facing the coercive temptations related to security.

Difficulties. « The Latin word Nox designates a period of time in which the sun disappears in the horizon » (Dictionnaire Zedler). This definition is no longer apt to describe the night of our cities. It is necessary to go further, even if the night doesn't unveil itself that easily. In many languages, eight nocturnal hours almost seem like a negation of the day, an invitation for us to follow suit. Noite (Não-oito) Nuit (Non-huit), Night (No-Eight), Nacht (Nein-Acht); Notte; Noche: the night nurtures paradoxes and seems to have a dark pleasure is erasing clues. Those who want to capture it, know there is always a game of craps out there which makes the night evade any possible equation.

Paradoxes. One may expect the night to bring counsel, but, in our language, it is the day that enlightens positively. Always. At night, colors disappear, one distinguishes badly forms and contours. However, the night discloses man and the city. Seen from satellites or on a plane, it is the night, the millions of urban lights, which reveal the presence of man on earth. It was at night, guided by the stars, that maritime civilizations set off to discover the world, and Columbus discovered the New World. It was at night, from the height of a promontory, that the organization of cities was better understood. Last paradox: if the fairy of electricity killed the night, she gave birth to a magical couple: the city and its night. Without light there isn't urban nightlife. And what a show, a spectacle, it is! Last paradox: if the fairy of electricity killed the night, she gave birth to a magical couple: the city and its night. Without light there isn't urban nightlife. And what a show, a spectacle, it is!

A DEBATE TO BE INITIATED

In face of the explosion of social times and generalized social acceleration, urban nightlife is surely the last terrain wherein we may enunciate, 'we', 'making the city, society or nation'. In the incessant to and fro between individualization and socialization, me and us, the night still allows for resynchronizations. She continues to be a possible moment of pause, for the encounter and discovery of the other. We can take refuge in the night, resist and question anew our capacity of living together. In this temporarily autonomous zone, we can still define a moment of standstill and think together. Field of tensions and creation with the capacity of re-enchanting our cities and our lives. Faced with the promise of a 'commercialized' night, the questions of security and the temptations of a city which doesn't sleep, we can propose a night both responsive and pacific. It is a question of sustainable development for welcoming cities, both for their nights or for their days. The night -and time in general -is one of the questions pertaining to public policy and whose responsibility is totally transversal to society as a whole. It is an opportunity. Do we want to both maintain our traditional rhythms and allow ourselves to be in a society which doesn't sleep, as a synonym of comfort for some, of hell for others? Without public debate and without conciliation, we risk creating new tensions, reinforcing inequalities among those who take advantage of urban nightlife and the others obliged to work in order to keep the city going. Hiding these questions or allotting them to the private sphere, we allow the market economy to enforce its laws on the less socially favored, bearing in so doing the risk of seeing a set of isolated decisions spark off new conflicts and inequalities.

(

  *) Luc Gwiazdzinski was born in Lorraine, in the northeast of France. PhD in geography, professor and researcher in the Joseph Fourier University of Grenoble and in the Polytechnic of Milan, he is member of the lab Pacte (UMR 5194 CNRS which is associated to MOTU [Milan] and EIREST [Paris, Sorbonne]). European specialist on cities (questions of mobility, social temporalities and urban nightlife) and responsible for conferences and research programs, he also has published extensively on these thematics: La ville 24h/24, 2003, Aube; La nuit dernière frontière de la ville, 2005, Aube; Si la ville m'était contée (ed.), 2005, Eyrolles; La nuit en questions (ed.), 2005, Aube et Nuits d'Europe, Pour des villes accessibles et hospitalières, UTBM, 2007. Citer l'article: GWIAZDZINSKI L., 2014, « The nocturnal condition », Night Manifesto. Seeking Citizenshup 24h. Colaboratória. -São Paulo: Invisíveis Produções, pp. 53-70 Contact: luc.gwiazdzinski@ujf-grenoble.fr

During the pre-Christian era they peregrinated Fisterra and, for many centuries, it was considered as the westernmost spot of land of the known world and -as such -the end of the world. Until today, Cape Finisterre is the final destiny for many pilgrims on the Way of St. James.

The nychthemeron is a physiological unit of time defining, for man and for the majority of animals, a period of wakefulness and a period of sleep, under the cadency of day and night.


Birgitta Dresp

Psychophysical measures of two-dimensional curvature

We investigated 'curvedness' of 2-D line stimuli using Stevens' method of magnitude estimation, where observers had to assign a number to a given stimulus. This number was to reflect the subjective estimate of the magnitude of curvature a given curve conveyed. Planar geometry permits generating curved lines in the 2-D plane with consistent variations in critical curve parameters, such as sagitta, chord-length, and height-to-width ratio. The model that yields optimal psychometric functions of curvature was found to be a combination of halves of vertically and horizontally oriented ellipses derived from concentric circles. A linear function accounts for perceived curvature as a function of the sagitta and chordlength of such curves, and a power function for perceived curvature as a function of the curves' height-towidth ratio. Perceptual estimates of curvature in planar images are reliably predicted by the three geometrical parameters tested.

Introduction

In classical psychophysics, our immediate subjective responses to sensory data, commonly called 'sensations', are considered to reflect the expression of brain mechanisms processing spatial or temporal characteristics of physical objects. Fechner (1860), Wundt (1874), Thurstone (1927), and subsequently Stevens (1957) used psychophysical scales as measures for the variations in our subjective perceptual responses to variations in the physical characteristics of objects and events. The psychometric functions obtained through such measures are to account for internal representations of the loudness, brightness, size, taste, smell, warmth, tactile quality or duration of physical events, also called stimuli. The psychophysiological mechanisms such psychometric functions are believed to reflect are sometimes referred to in terms of 'psychophysical laws', or 'laws of sensation'. A well-known example of a 'psychophysical law' is Steven's power law, which claims that the relationship between the physical magnitude of a stimulus and its perceived strength or intensity is universally accounted for by a power function. Stevens' law is written in terms of:

 I = k I b
where I is the physical magnitude of the stimulus, d the subjective value (sensation) associated with it, k a constant of proportionality and b the exponent. According to Stevens (1957), k and b would depend on the kind of stimulus presented in the psychophysical experiment. In its classic form, the power law is written in terms of:

y = a x b
Stevens' law, adapted for the subjective magnitudes associated with the intensity of sounds, lights, tastes or other sensory data, is based on power fits to average magnitude estimation data from multiple observers without considering inter-individual differences, or analyzing individual psychometric functions. This somehow seems to imply that Stevens considered his power law in terms of a unique law of sensation, or the reflection of a universal psycho-physiological mechanism common to all individuals. Critical minds, however, have questioned since whether the idea of a 'unique law of sensation' stands ground or not (e.g. Green & Luce, 1974). Also, Stevens assumed that individuals can perceive and report numbers reflecting physical intensities or ratios in an absolute and veridical manner independent of contextual influence (see also Narens, 1996), which poses a further problem.

Is it realistic to assume that perceptual judgements can be made at any given moment of time without any internal reference to prior experience or immediate context? Power laws are among the most frequent scaling laws to account for a myriad of phenomena in biology, physics, or economics. Power laws are asymptotic distributions which, in contrast with other asymptotic laws, produce curves that approach but never meet or cross the theoretical asymptote, and thus tend toward infinity. Assuming that sensations elicited by physical intensities of stimuli obey a power law therefore implies that the underlying mechanism justifies, or can at least be reconciled with, the concept of infinite progression.

How does this take into account the fact that receptors of any given sensory modality are known to saturate at some intensity of stimulation? Or, how would a supposedly universal, infinitely asymptotic power law fare in the light of visual discrimination at the extreme end of the contrast sensitivity function (CSF), which shows that the human visual system cannot distinguish spatial frequencies higher than 60 cycles per degree of visual angle (e.g. Mannos & Sakrison, 1974) because the number of photoreceptors in the human visual system is insufficient? In other words, how can we reconcile the functional hypothesis of the infinitely progressing power law with all that is known about the capacity limits of sensory and cognitive information processing? Some critical investigations into the power law as a universal function for internal representations of sensory have been reported. Certain authors, for example, examined the validity of Stevens' general assumption that magnitude estimates scale with physical or numerical values, and found a dependence on the assumed nature of the internal scale, with evidence pointing toward relative scales or ratio scales rather than a veridical representation of absolute numbers (e.g. Narens, 1996). Also, the power law does not necessarily provide the best fit when individual rather than average data are considered (Green & Luce, 1974). Our contribution here is based on an analysis of individual psychometric functions from a magnitude estimation study combined with numerical simulations where different mathematical functions to account for variations in visual sensations of curvedness in the twodimensional plane as a function of the geometric parameters of the stimuli are compared.

The perception of curves appears to be largely determined by a statistical relationship between their geometric properties and physical object attributes (e.g. Attneave, 1954, Ogilive & Daicar, 1967, Morgan, 2005[START_REF] Dresp | Which geometric model for the perceived curvature of 2-D shape contours?[END_REF]. In a recent study, we [START_REF] Dresp | Which geometric model for the perceived curvature of 2-D shape contours?[END_REF] investigated the geometric representations underlying visual sensations of curvature, or curvedness, produced by stimuli represented in the two-dimensional plane. The latter consisted of arcs corresponding to the lower and upper halves of concentric circles and ellipses derived mathematically through planar projection by affinity with the circles. Such a descriptive geometric model permits generating curved lines in the 2-D plane with consistent variations in critical curve parameters, such as sagitta, chord-length, and height-to-width ratio. The model that was found to yield the best psychometric functions for visual sensations of curvedness was a combination of halves of vertically and horizontally oriented ellipses derived from concentric circles (see Figure 1). These elliptic arcs represented a sub-set of curves with optimal variations in height-to-width ratio, which is defined in terms of the ratio between the sagitta and the chord-length of a given curve.

The halves of concentric circles of varying diameter but with constant height-to-width ratio produced a poor psychometric function, represented by a flat data curve within the middle region of the subjective scale, which ranged from 0 to 10. The average subjective magnitudes produced by elliptic curves, in contrast, increased exponentially with their heightto-width ratio and linearly with their sagittae across the full range of the subjective scale (see [START_REF] Dresp | Which geometric model for the perceived curvature of 2-D shape contours?[END_REF]. The best psychometric fits were obtained with the height-to-width ratio.

Thus, visual sensations of 2-D curvature were shown to be largely size-invariant and best predicted by a geometric parameter that conveys spatial information relative to the visual area covered by a curve, which is given by the height-to-width ratio.

On the basis of these previous findings, we investigated the psychometric functions describing individual visual sensations of 2-D curvedness elicited by halves of vertically and horizontally oriented ellipses. A specific geometric transform was applied to generate elliptic arcs (see Figure 1 and the corresponding text in Stimuli) within the limits of spatial resolution of a computer screen. In such elliptic arcs, chord-length and sagitta are related to the heightto-width-ratio by exponential functions; sagitta is related to chord-length by a linear function.

The mathematical functions describing relations between the geometric curve parameters considered here are represented by the graphs in Figure 2.

Material and methods

In a psychophysical experiment investigating visual sensations relative to the 'curvedness' of 2-D line stimuli, we used Stevens' method of magnitude estimation, where observers had to assign a number to a given stimulus. Here, this number was to reflect the subjective estimate of the "amount", or magnitude, of curvature a given curve conveyed. The curvature of a visual object in terms of a perceptual quality is referred to as 'curvedness'.

Observers

Fourteen observers (seven women and seven men) with normal or corrected-to-normal vision participated in the experiments.

Stimuli

The stimuli for the experiments were generated by a computer graphic design software tool (AUTOCAD) on an IBM computer (Pentium III) equipped with a standard colour screen with a display resolution of 1024 × 768 pixels. Curves with 'positive' (upward) and 'negative' (downward) curvature in the two-dimensional plane were generated from arcs of 11 ellipses derived from concentric circles with varying diameter through planar projection by affinity (see Figure 1). The rationale of this procedure for generating curves is described in detail elsewhere [START_REF] Dresp | Which geometric model for the perceived curvature of 2-D shape contours?[END_REF]. The diameter of the concentric circles varied between 1 and 9 degrees of visual angle, leading to projected ellipses the arcs of which had sagittae (sagitta being the distance between the chord and the parallel line tangent to it) within a range of variations between 0.5 and 4.5 degrees of visual angle. Chord-lengths varied between 1 and 9 degrees of visual angle. The height-to-width ratio, defined as the ratio between the sagitta and the chord-length of a given curve, varied between 0.05 and 4.5. All stimuli were defined by a curved line of the width of a single pixel and a constant length of 9 degrees of visual angle.

The luminance of the curves (40 cd/m 2 ) was determined by means of a photometer (Minolta LS 110) used for calibrating grey levels (R-G-B combinations) of the computer screen. The dark background of the screen on which the curves were presented had a constant luminance (2 cd/m 2 ).

Procedure

Observers were seated comfortably in a semi-dark room in front of the computer screen. A perfectly straight, white horizontal line with zero curvature was presented on a computer screen at the beginning of the trials. This line, of the same length as the lines used for generating the curves, served as standard and was assigned the modulus 'zero'. Starting from there, observers were told that they were going to view a series of curves, one at a time, and were asked to produce a number between 0 and 10 that was to reflect the magnitude of the curvedness they spontaneously perceived when a given stimulus came up on the screen.

None of the 14 observers had ever participated in a psychophysical experiment or knew what a psychophysical scale was. The curves were presented in random order, and for each observer, a different random sequence of stimuli was generated. The duration of a presentation or trial was not limited, but observers were encouraged to make their decisions as rapidly as possible. Typing the 'enter' key triggered the presentation of the next stimulus.

Each curve was shown only once to a given observer within a single series of trials.

Results

The magnitudes of curvedness recorded from each observer in response to the different height-to-width ratios of the 22 curves were analyzed. Individual psychometric functions of subjective curve magnitude as a function of height-to-width ratio, sagitta and chordlength were plotted. "Positive" and "negative" curves with identical height-to-width ratios produced identical or very similar magnitudes, as could be expected from previous data reported by [START_REF] Dresp | Which geometric model for the perceived curvature of 2-D shape contours?[END_REF]. The individual data, averaged over the curve orientation factor; are shown in the graphs represented in the graphs in Figures 3,4 and 5. To account for variations in visual sensations of curvedness as a function of the different geometric parameters of the curves, different mathematical functions were fitted to the data, including Stevens' Power Law and the exponential-rise-to-maximum function. While the power law represents an asymptotic mechanism which continues to progress towards infinity, the exponential-rise-to-maximum law describes a hypothetical mechanism that, once the critical asymptotic level attained, levels out flat without progression toward infinity. The exponentialrise-to-maximum function is written in terms of:

 y = a (1-exp(-b x))
Stevens' psychophysical power law can be expressed in terms of:

 I = k (1-exp(-b I))

Discussion

The data show that planar geometry permits generating curved lines in the 2-D plane with consistent variations in critical curve parameters, such as sagitta, chord-length, and height-towidth ratio. The model that yields optimal psychometric functions for visual sensations of curvedness is a combination of halves of vertically and horizontally oriented ellipses derived from concentric circles (see Figure 1). These elliptic arcs represent a sub-set of curves with variations in height-to-width ratio, defined in terms of the ratio between the sagitta and the chord-length of a given curve. Psychometric functions such as Steven's power law predict perceived 2D curvature as a function of the height-to-width ratio. Linear functions predict perception as a function of the sagitta or chordlengths of the curves. We conclude that perceptual estimates of curvature in planar images are reliably predicted by these three geometric parameters.

Figure captions

Figure 1 Vertically (a) and horizontally (b) oriented ellipses were obtained from concentric circles through a geometric transform called planar projection by affinity (see also [START_REF] Dresp | Which geometric model for the perceived curvature of 2-D shape contours?[END_REF]. In Cartesian space, an ellipse may be defined as the projected image of two concentric circles. In the two examples given here in (a) and (b), images (x, y) = (b x , a y ) of the principal circle C( 0 , a ) and images (x, y) = ((a/b) x, y) of the secondary circle C( 0,b ) generate ellipses through planar projection by affinity with the two circles. In our study, the upper and lower halves of eleven such ellipses (c), derived from concentric circles with varying diameter, were taken to generate two sets of eleven curves with 'positive' (d) and 'negative' (e) curvature in the two-dimensional plane. The sagitta of a curve described by one half of a given ellipse is identical to the radius (r) of the principal circle the full ellipse is derived from; the chordlength of such a curve is identical to twice the radius (r x 2) of the secondary circle the full ellipse is derived from.

Figure 2

 2 Figure 2 Arcs of ellipses, generated here within the limits of spatial resolution of the computer screen, possess geometric parameters which are linked by specific mathematical functions. The sagittae and chord-lengths of elliptic arcs are related to their height and width-ratios by exponential functions (see upper and middle graphs). The sagittae are linked to the chordlengths by a linear function (see lower graph).

Figure 3

 3 Figure 3Individual data (13 graphs) showing psychometric functions of subjective curvature as a function of the height-to-width ratio.

Figure 4

 4 Figure 4 Individual data (13 graphs) showing psychometric functions of subjective curvature as a function of the sagitta.

Figure 5

 5 Figure 5 Individual data (13 graphs) showing psychometric functions of subjective curvature as a function of the chordlength.


Fabian Pedregosa
email: fabian.pedregosa@inria.fr.

Michael Eickenberg

Philippe Ciuciu

Bertrand Thirion

Alexandre Gramfort

Data-driven HRF estimation for encoding and decoding models

Keywords: Functional MRI (fMRI), Hemodynamic response function (HRF), machine learning, optimization, BOLD, Finite inpulse response (FIR), decoding, encoding

come

Introduction

The use of machine learning techniques to predict the cognitive state of a subject from their functional MRI (fMRI) data recorded during task performance has become a popular analysis approach for neuroimaging studies over the last decade [START_REF] Cox | Functional magnetic resonance imaging (fMRI) brain reading: detecting and classifying distributed patterns of fMRI activity in human visual cortex[END_REF][START_REF] Haynes | Decoding mental states from brain activity in humans[END_REF]. It is now commonly referred to as brain reading or decoding. In this setting, the BOLD signal is used to predict the task or stimulus that the subject was performing. Although it is possible to perform decoding directly on raw BOLD signal [START_REF] Miranda | Dynamic discrimination analysis: a spatial-temporal SVM[END_REF][START_REF] Miyawaki | Visual image reconstruction from human brain activity using a combination of multiscale local image decoders[END_REF], the common approach in fast event-related designs consists in extracting the activation coefficients (beta-maps) from the BOLD signal to perform the decoding analysis on these estimates. Similarly, in the voxel-based encoding models [START_REF] Kay | Identifying natural images from human brain activity[END_REF][START_REF] Naselaris | Encoding and decoding in fMRI[END_REF], the activation coefficients are extracted from the BOLD signal, this time to learn a model to predict the BOLD response in a given voxel, based on a given representation of the stimuli. In addition, a third approach, known as representational similarity analysis or RSA [START_REF] Kriegeskorte | Representational similarity analysis-connecting the branches of systems neuroscience[END_REF] takes as input the activation coefficients. In this case a comparison is made between the similarity observed in the activation coefficients, quantified by a correlation measure, and the similarity between the stimuli, quantified by a similarity measure defined from the experimental setting.

These activation coefficients are computed by means of the General Linear Model (GLM) [START_REF] Friston | Statistical parametric maps in functional imaging : A general linear approach[END_REF]. While this approach has been successfully used in a wide range of studies, it does suffer from limitations [START_REF] Poline | The general linear model and fMRI: does love last forever?[END_REF]. For instance, the GLM commonly relies on a data-independent canonical form of the hemodynamic response function (HRF) to estimate the activation coefficient. However it is known [START_REF] Handwerker | Variation of BOLD hemodynamic responses across subjects and brain regions and their effects on statistical analyses[END_REF][START_REF] Badillo | Group-level impacts of within-and between-subject hemodynamic variability in fMRI[END_REF] that the shape of this response function can vary substantially across subjects and brain regions. This suggests that an adaptive modeling of this response function should improve the accuracy of subsequent analysis.

To overcome the aforementioned limitation, Finite Impulse Response (FIR) models have been proposed within the GLM framework [START_REF] Dale | Optimal experimental design for event-related fMRI[END_REF][START_REF] Glover | Deconvolution of impulse response in event-related BOLD fMRI[END_REF]. These models do not assume any particular shape for the HRF and amount to estimating a large number of parameters in order to identify it. While the FIRbased modeling makes it possible to estimate the activation coefficient and the HRF simultaneously, the increased flexibility has a cost. The estimator is less robust and prone to overfitting, i.e. to generalize badly to unseen data. In general, FIR models are most appropriate for studies focused on the characterization of the shape of the hemodynamic response, and not for studies that are primarily focused on detecting activation (Poldrack et al., 2011, Chapter 5).

Several strategies aiming at reducing the number of degrees of freedom of the FIR model -and thus at limiting the risk of overfitting -have been proposed. One possibility is to constrain the shape of the HRF to be a linear combination of a small number of basis functions. A common choice of basis is formed by three elements consisting of a reference HRF as well as its time and dispersion derivatives [START_REF] Friston | Statistical parametric maps in functional imaging : A general linear approach[END_REF], although it is also possible to compute a basis set that spans a desired function space [START_REF] Woolrich | Constrained linear basis sets for HRF modelling using variational bayes[END_REF].

More generally, one can also define a parametric model of the HRF and estimate the parameters that best fit this function [START_REF] Lindquist | Validity and power in hemodynamic response modeling: A comparison study and a new approach[END_REF]. However, in this case the estimated HRF may no longer be a linear function of the input parameters.

Sensitivity to noise and overfitting can also be reduced through regularization. For example, temporal regularization has been used in the smooth FIR [START_REF] Goutte | Modeling the haemodynamic response in fMRI using smooth FIR filters[END_REF][START_REF] Ciuciu | Unsupervised robust nonparametric estimation of the hemodynamic response function for any fMRI experiment[END_REF][START_REF] Casanova | The impact of temporal regularization on estimates of the BOLD hemodynamic response function: a comparative analysis[END_REF] to favor solutions with small second order time derivative. These approaches require the setting of one or several hyperparameters, at the voxel or potentially at the parcel level (if several voxels in a pre-defined parcel are assumed to share some aspects of the HRF timecourse). Even if efficient techniques such as generalized cross-validation [START_REF] Golub | Generalized cross-validation as a method for choosing a good ridge parameter[END_REF] can be used to choose the regularization parameters, these methods are inherently more costly than basis-constrained methods. Basis-constrained methods also require setting the number of basis elements; however, this parameter is not continuous (as in the case of regularized methods), and in practice only few values are explored: for example the 3-element basis set formed by a reference HRF plus derivatives and the FIR model. This paper focuses on basis-constrained regularization of the HRF to avoid dealing with hyperparameter selection with the goal of remaining computationally attractive. A different approach to increase robustness of the estimates consists in linking the estimated HRFs across a predefined brain parcel, taking advantage of the spatially dependent nature of fMRI [START_REF] Wang | Multiscale adaptive smoothing models for the hemodynamic response function in fMRI[END_REF]. However, hemodynamically-informed parcellations [START_REF] Chaari | Hemodynamic-informed parcellation of fMRI data in a joint detection estimation framework[END_REF]Badillo et al., 2013a) rely on the computation of a large number of estimations at the voxel or sub-parcel level. In this setting, the development of voxel-wise estimation procedures is complementary to the development of parcellation methods in that more robust estimation methods at the voxel level would naturally translate into more robust parcellation methods. In this paper we focus on voxel-wise estimation methods.

We propose a method for the simultaneous estimation of HRF and activation coefficients based on lowrank modeling. Within this model, and as in [START_REF] Makni | Bayesian deconvolution of fMRI data using bilinear dynamical systems[END_REF][START_REF] Kay | Identifying natural images from human brain activity[END_REF][START_REF] Vincent | Spatially adaptive mixture modeling for analysis of fMRI time series[END_REF][START_REF] Degras | A hierarchical model for simultaneous detection and estimation in multi-subject fMRI studies[END_REF], the HRF is constrained to be equal across the different conditions, yet permitting it to be different across voxels. Unlike previous works, we formulate this model as a constrained least squares problem, where the vector of coefficients is constrained to lie within the space of rank one matrices. We formulate the model within the framework of smooth optimization and use quasi-Newton methods to find the vector of estimates. This model was briefly presented in the conference paper [START_REF] Pedregosa | HRF estimation improves sensitivity of fMRI encoding and decoding models[END_REF].

Here we provide more experimental validation and a more detailed presentation of the method. We also added results using a GLM with separate designs [START_REF] Mumford | Deconvolving BOLD activation in event-related designs for multivoxel pattern classification analyses[END_REF]. Ten alternative approaches are now compared on two publicly available datasets. The solver has also been significantly improved to scale to full brain data.

The contributions of this paper are two-fold. First, we quantify the importance of HRF estimation in encoding and decoding models. While the benefit of data-driven estimates of the HRF have already been reported in the case of decoding [START_REF] Turner | Spatiotemporal activity estimation for multivoxel pattern analysis with rapid event-related designs[END_REF] and encoding approaches [START_REF] Vu | Encoding and decoding V1 fMRI responses to natural images with sparse nonparametric models[END_REF], we here provide a comprehensive comparison of models. Second, we evaluate a method called GLM with Rank-1 constraint (R1-GLM) that improves encoding and decoding scores over state-of-the-art methods while remaining computationally tractable on a full brain volume. We propose an efficient algorithm for this method and discuss practical issues such as initialization. Finally, we provide access to an open source software implementation of the methods discussed in this paper.

Notation:

• and • ∞ denote the Euclidean and infinity norm for vectors. We use lowercase boldface letter to denote vectors and uppercase boldface letter to denote matrices. I denotes the identity matrix, 1 n denotes the vector of ones of size n, ⊗ denotes the Kronecker product and vec(A) denotes the concatenation of the columns of a matrix A into a single column vector. A † denotes the Moore-Penrose pseudoinverse. Given the vectors {a 1 , . . . , a k } with a i ∈ R n for each 1 ≤ i ≤ k, we will use the notation [a 1 , . . . , a k ] ∈ R n×k to represents the columnwise concatenation of the k vectors into a matrix of size n × k.

We will use Matlab-style colon notation to denote slices of an array, that is x(1 : k) will denote the first k elements of x.

Methods

In this section we describe different methods for extracting the HRF and activation coefficients from BOLD signals. We will refer to each different stimulus as condition and we will call trial a unique presentation of a given stimulus. We will denote by k the total number of stimuli, y ∈ R n the BOLD signal at a single voxel and n the total number of images acquired.

The General Linear Model

The original GLM model [START_REF] Friston | Statistical parametric maps in functional imaging : A general linear approach[END_REF] makes the assumption that the hemodynamic response is a linear transformation of the underlying neuronal signal. We define the n×k-matrix X GLM as the columnwise stacking of different regressors, each one defined as the convolution of a reference HRF [START_REF] Boynton | Linear Systems Analysis of Functional Magnetic Resonance Imaging in[END_REF][START_REF] Glover | Deconvolution of impulse response in event-related BOLD fMRI[END_REF] with the stimulus onsets for the given condition. In this work we used as reference HRF the one provided by the software SPM 8 [START_REF] Friston | Statistical parametric mapping: The analysis of functional brain images: The analysis of functional brain images[END_REF]. Assuming additive white noise, n ≥ k and X GLM to be full rank, the vector of estimates is given by βGLM = X † GLM y, where βGLM is a vector of size k representing the amplitude of each one of the conditions in a given voxel.

A popular modification of this setting consists in extending the GLM design matrix with the temporal and width derivatives of the reference HRF. This basis, formed by the reference HRF and its derivatives with respect to time and width parameters, will be used throughout this work. We will refer to it as the 3HRF basis. In this case, each one of the basis elements is convolved with the stimulus onsets of each condition, obtaining a design matrix of size n × 3k. This way, for each condition, we estimate the form of the HRF as a sum of basis functions that correspond to the first order Taylor expansion of the parametrization of the response function. Another basis set that will be used is the Finite Impulse Response (FIR) set. This basis set spans the complete ambient vector space (of dimension corresponding to the length of the impulse response) and it is thus a flexible model for capturing the HRF shape. It consists of the canonical unit vectors (also known as stick function) for the given duration of the estimated HRF. Other basis functions such as FMRIB's Linear Optimal Basis Sets [START_REF] Woolrich | Constrained linear basis sets for HRF modelling using variational bayes[END_REF] are equally possible but were not considered in this work.

More generally, one can extend this approach to any set of basis functions. Given the matrix formed by the stacking of d basis elements B = [b 1 , b 2 , . . . , b d ], the design matrix X B is formed by successively stacking the regressors obtained by convolving each of the basis elements with the stimulus onsets of each condition. This results in a matrix of size n × dk and under the aforementioned conditions the vector of estimates is given by βB = X † B y. In this case, βB is no longer a vector of size k: it has length k × d instead and can no longer be interpreted as the amplitude of the activation. One possibility to recover the trial-by-trial reponse amplitude is to select the parameters from a single time point as done by some of the models considered in [START_REF] Mumford | Deconvolving BOLD activation in event-related designs for multivoxel pattern classification analyses[END_REF], however this procedure assumes that the peak BOLD response is located at that time point. Another possibility is to construct the estimated HRF and take as amplitude coefficient the peak amplitude of this estimated HRF. This is the approach that we have used in this paper.

GLM with rank constraint

In the basis-constrained GLM model, the HRF estimation is performed independently for each condition. This method works reliably whenever the number of conditions is small, but in experimental designs with a large number of conditions it performs poorly due to the limited conditioning of the problem and the increasing variance of the estimates.

At a given voxel, it is expected that for similar stimuli the estimated HRF are also similar [START_REF] Henson | Detecting latency differences in event-related BOLD responses: application to words versus nonwords and initial versus repeated face presentations[END_REF]. Hence, a natural idea is to promote a common HRF across the various stimuli (given that they are sufficiently similar), which should result in more robust estimates [START_REF] Makni | Bayesian deconvolution of fMRI data using bilinear dynamical systems[END_REF][START_REF] Vincent | Spatially adaptive mixture modeling for analysis of fMRI time series[END_REF].

In this work we consider a model in which a common HRF is shared across the different stimuli. Besides the estimation of the HRF, a unique coefficient is obtained per column of our event matrix. This amounts to the estimation of k + d free parameters instead of k × d as in the standard basis-constrained GLM setting.

The novelty of our method stems from the observation that the formulation of the GLM model with a common HRF across conditions translates to a rank constraint on the vector of estimates. This assumption amounts to enforcing the vector of estimates to be of the form

β B = [hβ 1 , hβ 2 , • • • , hβ k ]
for some HRF h ∈ R d and a vector of coefficients β ∈ R k . More compactly, this can be written as β B = vec(hβ T ). This can be seen as a constraint on the vector of coefficients to be the vectorization of a rank-one matrix, hence the name Rank-1 GLM (R1-GLM).

In this model, the coefficients have no longer a closed form expressions, but can be estimated by minimizing the mean squared error of a bilinear model. Given X B and y as before, Z ∈ R n×q a matrix of nuisance parameters such as drift regressors, we define F R1 (h, β, ω, X B , y, Z) = 1 2 y -X B vec(hβ T ) -Zω 2 to be the objective function to be minimized. The optimization problem reads: ĥ, β, ω = arg min h,β,ω

F R1 (h, β, ω, X B , y, Z) subject to Bh ∞ = 1 and Bh, h ref > 0 , (1)
The norm constraint is added to avoid the scale ambiguity between h and β and the sign is chosen so that the estimated HRF correlates positively with a given reference HRF h ref . Otherwise the signs of the HRF and β can be simultaneously flipped without changing the value of the cost function. Within its feasible set, the optimization problem is smooth and is convex with respect to h, β and ω, however it is not jointly convex in variables h, β and ω.

From a practical point of view this formulation has a number of advantages. First, in contrast with the GLM without rank-1 constraint the estimated coefficients are already factored into the estimated HRF and the activation coefficients. That is, once the estimation of the model parameters from Eq. ( 1) is obtained, β is a vector of size k and ĥ is a vector of size d that can be both used in subsequent analysis, while in models without rank-1 constraint only the vector of coefficients (equivalent to vec(hβ T ) in rank-1 constrained models) of size k × d is estimated. In the latter case, the estimated HRF and the beta-maps still have to be extracted from this vector by methods such as normalization by the peak of the HRF, averaging or projecting to the set of Rank-1 matrices.

Second, it is readily adapted to prediction on unseen trials. While for classical (non rank-1 models) the HRF estimation is performed per condition with no HRF associated with unseen conditions, in this setting, because the estimated HRF is linked and equal across conditions it is natural to use this estimate on unseen conditions. This setting occurs often in encoding models where prediction on unseen trials is part of the cross-validation procedure.

This model can also be extended to a parametric HRF model. That is, given the hemodynamic response defined as a function h : R d1 → R d of some parameters α, we can formulate the analogous model of Eq. ( 1)

as an optimization over the parameters α and β with the design matrix X FIR given by the convolution of the event matrix with the FIR basis:

α, β, ω = arg min α,β,ω F R1 (h(α), β, ω, X FIR , y, Z) subject to h(α) ∞ = 1 and h(α), h ref > 0 (2)
In section 2.4 we will discuss optimization strategies for both models.

Extension to separate designs

An extension to the classical GLM that improves the estimation with correlated designs was proposed in [START_REF] Mumford | Deconvolving BOLD activation in event-related designs for multivoxel pattern classification analyses[END_REF]. In this setting, each voxel is modeled as a linear combination of two regressors in a design matrix X GLM . The first one is the regressor associated with a given condition and the second one is the sum of all other regressors. This results in k design matrices, one for each condition. The estimate for a given condition is given by the first element in the two-dimensional array X Si † y, where X Si is the design matrix for condition i. We will denote this model GLM with separate designs (GLMS). It has been reported

to find a better estimate in rapid event designs leading to a boost in accuracy for decoding tasks [START_REF] Mumford | Deconvolving BOLD activation in event-related designs for multivoxel pattern classification analyses[END_REF][START_REF] Schoenmakers | Linear reconstruction of perceived images from human brain activity[END_REF][START_REF] Lei | A mixed L2 norm regularized HRF estimation method for rapid event-related fMRI experiments[END_REF].

This approach was further extended in [START_REF] Turner | Spatiotemporal activity estimation for multivoxel pattern analysis with rapid event-related designs[END_REF] to include FIR basis instead of the predefined canonical function. Here we employ it in the more general setting of a predefined basis set. Given a set of basis functions we construct the design matrix for condition i as the columnwise concatenation of two matrices X 0 BSi and X 1 BSi . X 0 BSi is given by the columns associated with the current condition in the GLM matrix and X 1

BSi is the sum of all other columns. In this case, the vector of estimates is given by the first d vectors of X † BSi y. See [START_REF] Turner | Spatiotemporal activity estimation for multivoxel pattern analysis with rapid event-related designs[END_REF] for a more complete description of the matrices X 0 BSi and X 1 BSi . It is possible to use the same rank-1 constraint as before in the setting of separate designs, linking the HRF across conditions. We will refer to this model as Rank-1 GLM with separate designs (R1-GLMS). In this case the objective function has the form

F R1-S (h, β, ω, r, X B , y, Z) = 1 2 k i y-β i X 0 BSi h-r i X 1 BSi h-Zω 2
, where r ∈ R d is a vector representing the activation of all events except the event of interest and will not be used in subsequent analyses. We can compute the vector of estimates β as the solution to the optimization problem β, ω, ĥ, r = arg min

h,β,ω,r F R1-S (h, β, ω, r, X B , y, Z) subject to Bh ∞ = 1 and Bh, h ref > 0 (3)

Optimization

For the estimation of rank-1 models on a full brain volume, a model is estimate at each voxel separately.

Since a typical brain volume contains more than 40,000 voxels, the efficiency of the estimation at a single voxel is of great importance. In this section we will detail an efficient procedure based on quasi-Newton methods for the estimation of R1-GLM and R1-GLMS models on a given voxel.

One approach to minimize (1) is to alternate the minimization with respect to the variables β, h and ω. By recalling the Kronecker product identities (Horn and Johnson, 1991, Chapter 4.3), and using the identity vec(hβ T ) = β ⊗ h we can rewrite the objective function (1) to be minimized as:

1 2 y -X B (β ⊗ h) -Zω 2 = 1 2 y -X B (I ⊗ h)β -Zω 2 = 1 2 y -X B (β ⊗ I)h -Zω 2 . ( 4
)
Updating h, β or ω sequentially thus amounts to solving a (constrained) least squares problem at each iteration. A similar procedure is detailed in [START_REF] Degras | A hierarchical model for simultaneous detection and estimation in multi-subject fMRI studies[END_REF]. However, this approach requires computing the matrices X B (β ⊗ I) and X B (I ⊗ h) at each iteration, which are typically dense, resulting in a high computational cost per iteration. Note also that the optimization problem is not jointly convex in variables h, β, ω, therefore we cannot apply convergence guarantees from convex analysis.

We rather propose a more efficient approach by optimizing both variables jointly. We define a global variable z as the concatenation of (h, β, ω) into a single vector, z = vec([h, β, ω]), and cast the problem as an optimization with respect to this new variable. Generic solvers for numerical optimization [START_REF] Nocedal | Numerical optimization, series in operations research and financial engineering[END_REF] can then be used. The solvers that we will consider take as input an objective function and its gradient. In this case, the partial derivatives with respect to variable z can be written as

∂F R1 /∂z = vec([∂F R1 /∂h, ∂F R1 /∂β, ∂F R1 /∂ω])
, whose expression can be easily derived using the aforementioned Kronecker product identities:

               ∂F R1 ∂h = -(β T ⊗ I)X T (y -X vec(hβ T ) -Zω) ∂F R1 ∂β = -(I ⊗ h T )X T (y -X vec(hβ T ) -Zω) ∂F R1 ∂ω = -Z T (y -X vec(hβ T ) -Zω)
If instead a parametric model of the HRF is used as in Eq. ( 2), the equivalent partial derivatives can be easily computed by the chain rule.

For the sake of efficiency, it is essential to avoid evaluating the Kronecker products naively, but rather reformulate them using the above mentioned Kronecker identities. For example, the matrix M = X(I ⊗ h)

should not be computed explicitly but should rather be stored as a linear operator such that when applied to a vector a ∈ R k it computes M (a) = X(a ⊗ h), avoiding thus the explicit computation of I ⊗ β.

Similar equations can be derived for the rank-1 model with separate designs of Eq. ( 3) (R1-GLMS), in which case the variable z is defined as the concatenation of (h, β, ω, r), i.e. z = vec ([h, β, ω, r]). The gradient of F R1-S with respect to z can be computed as ∂F R1-S /∂z = vec([∂F R1-S /∂h, ∂F R1-S /∂β, ∂F R1-S /∂ω, F R1-S /∂r]).

The partial derivatives read:

               ∂F ∂h = k i -(X 0 BSi β i -X 1 BSi r i ) T (y -β i X 0 BSi h -w i X 1 BSi h) ∂F ∂βi = -(X 0 BSi h) T (y -β i X 0 BSi h -w i X 1 BSi h) ∂F ∂ωi = -Z T (y -β i X 0 BSi h -w i X 1 BSi h) ∂F ∂ri = -(X 1 BSi h) T (y -β i X 0 BSi h -w i X 1 BSi h)
A good initialization plays a crucial role in the convergence of any iterative algorithm. Furthermore, for non-convex problems a good initialization prevents the algorithm from converging to undesired local minima. We have used as initialization for the R1-GLM and R1-GLMS models the solution given by the GLM with separate designs (GLMS). Since the GLM with separate designs scales linearly in the number of voxels, this significantly reduces computation time whenever an important number of voxels is considered.

Whenever the design matrix X B has more rows than columns (as is the case in both datasets we consider with B the 3HRF basis), it is possible to find an orthogonal transformation that significantly speeds up the computation of the Rank-1 model. Let Q, R be the "thin" QR decomposition of X B ∈ R n×dk , that is, QR = X B with Q ∈ R n×dk an orthogonal matrix and R ∈ R dk×dk a triangular matrix. Because of the invariance of the Euclidean norm to orthogonal transformations, the change of variable X B ← Q T X B , y ← Q T y yields a Rank-1 model in Eq. ( 1) with equivalent solutions. This reduces the size of the design matrix to a square triangular matrix of size dk × dk (instead of n × dk) and reduces the explained variable y to a vector of size kd (instead of n). After this change of variable, the convergence of the Rank-1 model is significantly faster due to the faster computation of the objective function and its partial derivatives. We have observed that the total running time of the algorithm can be reduced by 30% using this transformation.

Some numerical solvers such as L-BFGS-B [START_REF] Liu | On the limited memory BFGS method for large scale optimization[END_REF] require the constraints to be given as box constraints. While our original problem includes an equality constraint we can easily adapt it to use convex box constraints instead. We replace the equality constraint Bh ∞ = 1 by the convex inequality constraint Bh ∞ ≤ 1, which is equivalent to the box constraint -1 ≤ (Bh) i ≤ 1 supported by the above solver. However, this change of constraint allows solutions in which h can be arbitrarily close to zero. To avoid such degenerate cases we add the smooth term -Bh 2 2 to the cost function. Since there is a free scale parameter between h and β, this does not bias the problem, but forces Bh to lie as far as possible from the origin (thus saturating the box constraints). Once a descent direction has been found by the L-BFGS-B method we perform a line search procedure to determine the step length. The line-search procedure was implemented to satisfy the strong Wolfe conditions [START_REF] Nocedal | Numerical optimization, series in operations research and financial engineering[END_REF]. Finally, when the optimization algorithm has converged to a stationary point, we rescale the solution setting to ensure that the equality constraint. This still leaves a sign ambiguity between the estimated HRF and the associated beta-maps. To make these parameters identifiable, the sign of the estimated HRF will be chosen so that these correlate positively with the reference HRF.

We have compared several first-order (Conjugate Gradient), quasi-Newton (L-BFGS) and Newton methods on this problems and found that in general quasi-Newton methods performed best in terms of computation time. In our implementation, we adopt the L-BFGS-B as the default solver.

In Algorithm 1 we describe an algorithm based on L-BFGS that can be used to optimize R1-GLM and R1-GLMS models (a reference implementation for the Python language is described in subsection Software).

Variable r is only used for the R1-GLMS method and its use is denoted within parenthesis, i.e. (, r), so that for the R1-GLM it can simply be ignored.

The full estimation of the R1-GLM model with 3HRF basis for one subject of the dataset described in section Dataset 2: decoding of potential gain levels (16 × 3 conditions, 720 time points, 41, 622 voxels) took 14 minutes in a 8-cores Intel Xeon 2.67GHz machine. The total running time for the 17 subjects was less than four hours.

Software

We provide a software implementation of all the models discussed in this section in the freely available (BSD licensed) pure-Python package hrf estimation5 .

Algorithm 1 Optimization of R1-GLM and R1-GLMS models Input: Given initial points

β 0 ∈ R k , h 0 ∈ R d , ω 0 ∈ R q (, r 0 ∈ R k ), convergence tolerance ǫ > 0, inverse Hessian approximation H 0 . Output: β m , h m 1: (Optional): Compute the QR decomposition of X B , QR = X B , and replace X B ← Q T X B , y ← Q T y 2: Initialization. Set m ← 0, z ← vec([h 0 , β 0 , ω 0 (, r 0 )]) 3: while ∇f > ǫ do 4: Compute search direction. Set p m ← -H m ∇f (h m , β m , ω m (, r m )
) by means of the L-BFGS algorithm.

5:

Set z m+1 = z m + γ m p m , where γ m is computed from a line search procedure subject to the box constraints h m ∞ ≤ 1.

6:

m ← m + 1 7: end while

8: Extract R1-GLM(S) parameters from z m . Set h m ← z m (1 : d), β m ← z m (d + 1 : m + d) 9:
Normalize and set sign so that the estimated HRF is positively correlated with a reference HRF:

q m ← h m ∞ sign(h T m h ref ), h m ← h m /q m , β m ← β m q m 3. Data description
With the aim of making the results in this paper easily reproducible, we have chosen two freely available datasets to validate our approach and to compare different HRF modeling techniques.

Dataset 1: encoding of visual information

The first dataset we will consider is described in [START_REF] Kay | Identifying natural images from human brain activity[END_REF][START_REF] Naselaris | Bayesian reconstruction of natural images from human brain activity[END_REF][START_REF] Kay | fMRI of human visual areas in response to natural images[END_REF].

It contains BOLD fMRI responses in human subjects viewing natural images. As in [START_REF] Kay | Identifying natural images from human brain activity[END_REF] run of the first session for that subject. Across-session alignment was performed manually. Additionally, data were temporally interpolated to account for slice-timing differences. See [START_REF] Kay | Identifying natural images from human brain activity[END_REF] for further preprocessing details.

We performed local detrending using a Savitzky-Golay filter [START_REF] Savitzky | Smoothing and differentiation of data by simplified least squares procedures[END_REF] Ridge regression (regularization parameter chosen by Generalized Cross-Validation [START_REF] Golub | Generalized cross-validation as a method for choosing a good ridge parameter[END_REF])

was then used to learn a predictor of voxel activity on the training set. By using this encoding model and the estimated HRF it is possible to predict the BOLD signal for the 70 images in the test set (20 % of the original session). We emphasize that learning the HRF on the training set instead of on the full dataset is necessary to avoid overfitting while assessing the quality of the estimated HRF by any HRF-learning method:

otherwise, the estimation of the HRF may incorporate specificities of the test set leading to artificially higher scores.

In a first step, we perform the image identification task from [START_REF] Kay | Identifying natural images from human brain activity[END_REF]. From the training set we estimate the activation coefficients that will be used to compute the activation maps. We use an encoding model using Gabor filters that predicts the activation coefficient from the training stimuli. From the stimuli in the validation set we predict the activation coefficients that we then use to identify the correct image.

The predicted image is the one yielding the highest correlation with the measured activity. This procedure mimics the one presented in (Kay et al., 2008, Supplementary material).

In a second step, we report score as the Pearson correlation between the measurements and the predicted BOLD signal on left out data. The prediction of BOLD signal on the test set is performed from conditions that were not present in the train set. In order to do this, an HRF for these conditions is necessary. As highlighted in the methods section, the construction of an HRF for these conditions is ambiguous for non Rank-1 methods that perform HRF estimation on the different stimuli. In these cases we chose to use the mean HRF across conditions as the HRF for unseen conditions. Finally, linear predictions on the left out fold were compared to the measured BOLD signals.

Dataset 2: decoding of potential gain levels

The second dataset described in [START_REF] Tom | The neural basis of loss aversion in decision-making under risk[END_REF] is a gambling task where each of the 17 subjects was asked to accept or reject gambles that offered a 50/50 chance of gaining or losing money. The magnitude of the potential gain and loss was independently varied across 16 levels between trials. Each gamble has an amount of potential gains and potential losses that can be used as class label. In this experiment, we only considered gain levels. This leads to the challenge of predicting or decoding the gain level from brain images.

The dataset is publicly available from http://openfmri.org under the name mixed-gambles task dataset.

The data preprocessing included slice timing, motion correction, coregistration to the anatomical images, tissue segmentation, normalization to MNI space and was performed using the SPM 8 software through the Pypreprocess6 interface.

For all subjects three runs were recorded, each consisting of 240 images with a repetition time (TR) of 2 seconds and a stimulus presentation at every 4 seconds. In order to perform HRF estimation on more data than what is available on a single run, we performed the estimation on the three runs simultaneously. This assumes HRF consistency across runs, which was obtained by concatenating the data from the three runs and creating a block-diagonal design matrix correspondingly (each block is the design of one run).

After training a regression model on 90% of the data, we predict the gain level on the remaining 10%.

As a performance measure we use Kendall tau rank correlation coefficient [START_REF] Kendall | A new measure of rank correlation[END_REF] between the true gain levels and the predicted levels, which is a measure for the orderings of the data. We argue that this evaluation metric is better suited than a regression loss for this task because of the discrete and ordered nature of the labels. Also, this loss is less sensible to shrinkage of the prediction that might occur when penalizing a regression model [START_REF] Bekhti | Decoding perceptual thresholds from MEG/EEG[END_REF]. The Kendall tau coefficient always lies within the interval [-1, 1], with 1 being perfect agreement between the two rankings and -1 perfect disagreement. Chance level lies at zero. This metric was previously proposed for fMRI decoding with ordered labels in [START_REF] Doyle | Multivariate decoding of brain images using ordinal regression[END_REF].

Results

In order to compare the different methods discussed previously, we ran the same encoding and decoding studies while varying the estimation method for the activation coefficients (beta-maps). The methods we considered are standard GLM (denoted GLM), GLM with separate designs (GLMS), Rank-1 GLM (R1-GLM) and Rank-1 GLM with separate designs (R1-GLMS). For all these models we consider different basis sets for estimating the HRF: a set of three elements formed by the reference HRF and its time and dispersion derivative, a FIR basis set (of size 20 in the first dataset and of size 10 in the second dataset) formed by the canonical vectors and the single basis set formed by the reference HRF (denoted "fixed HRF"), which in this case is the HRF used by the SPM 8 software.

It should be reminded that the focus of this study is not the study of the HRF in itself (such as variability across subjects, tasks or regions) but instead its possible impact on the accuracy of encoding and decoding paradigms. For this reason we report encoding and decoding scores but we do not investigate any of the possible HRF variability factors.

Dataset 1: encoding of visual information

In the original study, 500 voxels were used to perform image identification. These voxels were selected as the voxels with the highest correlation with the true BOLD signal on left-out data using a (classical) GLM with the reference HRF. These voxels are therefore not the ones naturally benefiting the most from HRF estimation.

We first present the scores obtained in the image identification task for different variants of the GLM.

This can be seen in Figure 1. The displayed score is the count of correctly identified images over the total number of images (chance level is therefore at 1/120). The identification algorithm here only uses the betamaps obtained from the train and validation set. This makes the estimation of the HRF an intermediate result in this model. However, we expect that a correct estimation of the HRF directly translates into a better estimation of the activation coefficients in the sense of being able to acheive higher predictive accuracy.

Our results are consistent with this hypothesis and in this task the rank-one (R1) and glm-separate (GLMS) models outperform the classical GLM model. The benefits range from 0.9% for R1-GLM in subject 2 to 8.2% for the same method and subject 1. It is worth noticing that methods with FIR basis obtain a higher score than methods using the 3HRF basis.

In order to test whether this increase is statistically significant we performed the following statistical test. The success of recovering the correct image can be modeled as a binomial distribution, with p A being be the probability of recovering the correct image with method A and p B be the probability of recovering the correct image with method B. We define the null hypothesis H 0 as the statement that both probabilities are equal, H 0 : p A = p B , and the alternate hypothesis that both probabilities and not equal, H 1 : p 1 = p 2 (this test is sometimes known as the binomial proportion test [START_REF] Röhmel | Unconditional non-asymptotic one-sided tests for independent binomial proportions when the interest lies in showing non-inferiority and/or superiority[END_REF]). The score test statistic for the one-tailed test is T = (p A -p B )/ p(1 -p) 2 n , where p = (p A + p B )/2 and n is the number of repetitions, in this case n = 120. This statistic is normally distributed for large n. The p-value associated with this statistical test when comparing every model (by order of performance) with the model "GLM with with fixed HRF" is (0.10, 0.10, 0.15, 0.19, 0.21, 0.26, 0.5, 0.5, 0.82, 0.81) for the first subject and (0.18, 0.18, 0.25, 0.34, 0.34, 0.44, 0.5, 0.5, 0.86, 0.93) for the second.

We will now use a different metric for evaluating the performance of the encoding model. This metric is the Pearson correlation between the BOLD predicted by the encoding model and the true BOLD signal, averaged across voxels. We will compute the this metric on a left-out session, which results in five scores for each method, corresponding to each of the cross-validation folds. Given two methods, a Wilcoxon signedrank test can be used on these cross-validation scores to assess whether the score obtained by the two methods are significantly different. This way, irrespective of the variance across voxels, which is inherent to the study, we can reliably assess the relative ranking of the different models. In Figure 2 we show the scores for each method (averaged across sessions) and the p-value corresponding the Wilcoxon test between a given method and the previous one by order of performance.

GLMS GLM

Figure 1: Image identification score (higher is better) on two different subjects from the first dataset. The metric counts the number of correctly identified images over the total number of images (chance level is 1/120 ≈ 0.008). This metric is less sensitive to the shape of the HRF than the voxel-wise encoding score. The benefits range from 0.9% points to 8.2% points across R1-constrained methods and subjects. The highest score is achieved by a R1-GLM method with a FIR basis set for subject 1 and by a R1-GLMS with FIR basis for subject 2. Average Correlation Score, subject 2 R1-GLMS R1-GLM

GLM with FIR basis

GLMS GLM

Figure 2: Average correlation score (higher is better) on two different subjects from the first dataset. The average correlation score is the Pearson correlation between the predicted BOLD and the true BOLD signal on left-out session, averaged across voxels and sessions. Methods that perform constrained HRF estimation significantly outperform methods that use a fixed reference HRF. As for the image identification performance, the best performing method for subject 1 is the R1-GLM, while for subject 2 it is the R1-GLMS model, both with FIR basis. In underlined typography is the GLM with a fixed HRF which is the method used by default in most software distributions. A Wilcoxon signed-rank test is performed between each method and the next one in the ordered result list by considering the leave-one-session out cross-validation scores for each method. We report p-values to assess whether the score differences are statistically significant.

We observed in Figure 2 that methods that learn the HRF together with some sort of regularization (be it Rank-1 constraint or induced by separate designs) perform noticeably better than methods that perform unconstrained HRF estimation, highlighting the importance of a robust estimation of the HRF as opposed to a free estimation as performed by the standard GLM model with FIR basis. This suggests that R1 and GLMS methods permit including FIR basis sets while minimizing the risk of overfitting inherent to the classical GLM model.

We also observed that models using the GLM with separate designs from [START_REF] Mumford | Deconvolving BOLD activation in event-related designs for multivoxel pattern classification analyses[END_REF] perform significantly better on this dataset than the standard design, which is consistent with the purpose of these models. It improves estimation in highly correlated designs. The best performing model for both subjects in this task is the R1-GLMS with FIR basis, followed by the R1-GLM with FIR basis model for subject 1 and GLMS with FIR basis for subject 2. The difference between both models (Wilcoxon signed-rank test) was significant with a p-value < 10 -6 . Since the results for both subjects are similar, we will only use subject 1 for the rest of the figures.

To further inspect the results, we investigated the estimation and encoding scores at the voxel level. This provides some valuable information. For example, parameters such as time-to-peak, width and undershoot of the estimated HRF can be used to characterize the mis-modeling of a reference HRF for the current study. Also, a voxel-wise comparison of the different methods can be used to identify which voxels exhibit a greater improvement for a given method. In the upper part of Figure 3 we show the HRF estimated on the first subject by our best performing method (the Rank-1 with separate designs and FIR basis). For comparison we also present two commonly used reference HRFs: one used in the software SPM and one defined in (Glover, 1999, auditory study) and used by software such as NiPy 7 and fmristat 8 . Because the HRF estimation will fail on voxels for which there is not enough signal, we only show the estimated HRF for voxels for which the encoding score is above the mean encoding score. In this plot the time-to-peak of the estimated HRF is color coded. One can observe a substantial variability in the time to peak, confirming the existence of a non-negligeable variability of the estimated HRFs, even within a single subject and a single task. In particular, we found that only 50% of the estimated HRFs on the full brain volume peaked between 4.5 and 5.5 seconds.

In the lower part of Figure 3 we can see a scatter plot in which the coordinates of each point are the encoding scores with two different methods. The first coordinate (X-axis) is given by the score using a canonical GLM whilst the second coordinate (Y-axis) corresponds to the Rank-1 separate with FIR basis.

Points above the black diagonal exhibit a higher score with our method than with a canonical GLM. As previously, the color represents the time to peak of the estimated HRF. From this plot we can see that voxels that have a low correlation score using a canonical GLM do not gain significant improvement by using a 7 http://nipy.org 8 http://www.math.mcgill.ca/keith/fmristat/ The metric is Pearson correlation. Points above the black diagonal correspond to voxels that exhibit a higher score with the R1-GLMS method than with a standard GLM.

Rank-1 Separate FIR model instead. However, voxels that already exhibit a sufficiently high correlation score using a canonical GLM (> 0.05) see a significant increase in performance when estimated using our method.

These results suggest as a strategy to limit the computational cost of learning the HRF on an encoding study to perform first a standard GLM (or GLMS) on the full volume and then perform HRF estimation only on the best performing voxels.

The methods that we have considered for HRF estimation can be subdivided according to the design matrices they use (standard or separate) and the basis they use to generate the estimated HRF (3HRF and FIR). We now focus on the performance gains of each of these individual components. In the upper part of Figure 4 we consider the top-performing model, the Rank-1 GLMS, and compare the performance of two different basis sets: FIR with 20 elements in the Y-axis and the reference HRF plus its time and dispersion derivatives (3HRF) in the X-axis. The abundance of points above the diagonal demonstrates the superiority of the FIR basis on this dataset. The color trend in this plot suggests that the score improvement of the FIR basis with respect to the 3HRF basis becomes more pronounced as the time-to-peak of the estimated HRF deviates from the reference HRF (peak at 5s), which can be explained by observing that the 3HRF basis corresponds to a local model around the time-to-peak. In the bottom part of this figure we compare the different design matrices (standard or separate). Here we can see the voxel-wise encoding score for two Rank-1 models with FIR basis and different design matrices: separate design on the Y-axis and classical design on the X-axis. Although both models give similar results, a Wilcoxon signed-rank test on the leaveone-session-out cross-validation score confirmed the superiority of the separate designs model in this dataset with p-value < 10 -3 .

In Figure 5 we can see the voxel-wise encoding score on a single acquisition slice. In the upper column, the score is plotted on each voxel and thresholded at a value of 0.045, which would correspond to a p-value < 0.05 for testing non-correlation assuming each signal is normally distributed, while in the bottom row the 0.055 contour (p-value < 0.001) for the same data is shown as a green line. Here it can be seen how the top performing voxels follow the gray matter. A possible hypothesis to explain the increase of the encoding score between the method R1-GLMS with FIR basis and the same method with 3HRF basis could be related either to the shape of the HRF deviating more from a canonical shape in lateral visual areas or to the higher signal-to-noise ratio often found in the visual cortex when compared to lateral visual areas.

Dataset 2: decoding of potential gain levels

The mean decoding score was computed over 50 random splittings of the data, with a test set of size 10%. The decoding regression model consisted of univariate feature selection (ANOVA) followed by a Ridge regression classifier as implemented in scikit-learn [START_REF] Pedregosa | Scikit-learn : Machine learning in python[END_REF]. Both parameters, number of voxels and amount of ℓ 2 regularization in Ridge regression, were chosen by cross-validation. Average Decoding Score R1-GLMS R1-GLM

GLMS GLM

Figure 6: Averaged decoding score for the different method considered (higher is better) on the second dataset. The metric is Kendall tau. Methods that perform constrained HRF estimation significantly outperform methods that use a fixed (reference) HRF. In particular, the best performing method is the R1-GLM with 3HRF basis, followed by the R1-GLMS with 3HRF basis.

In underlined typography is the GLM with a fixed HRF which is the method used by default in most software distributions. As in Figure 2, a Wilcoxon signed-rank test is performed and the p-value reported between a given method and the next method in the ordered result list to assess whether the difference in score is significant.

The mean score for the 10 models considered can be seen in Figure 6. Similarly to how we assessed superiority of a given method in encoding, we will say that a given method outperforms another if the paired difference of both scores (this time across folds) is significantly greater than zero. This is computed by performing a Wilcoxon signed rank test across voxels. For this reason we report p-values together with the mean score in Figure 6.

As was the case in encoding, Rank-1 constrained methods obtain the highest scores. In this case however, methods with 3HRF basis outperform methods using FIR basis. This can be explained by factors such as smaller sample size of each of the runs, smaller number of trials in the dataset and experimental design.

Discussion

We have compared different HRF modeling techniques and examined their generalization score on two different datasets: one in which the main task was an encoding task and one in which it was a decoding task. We compared 10 different methods that share a common formulation within the context of the General Linear Model. This includes models with canonical and separate designs, with and without HRF estimation constrained by a basis set, and with and without rank-1 constraint. We have focused on voxel-independent models of the HRF, possibly constrained by a basis set, and have omitted for efficiency reasons other possible models such as Bayesian models [START_REF] Marrelec | Robust bayesian estimation of the hemodynamic response function in event-related BOLD fMRI using basic physiological information[END_REF][START_REF] Ciuciu | Unsupervised robust nonparametric estimation of the hemodynamic response function for any fMRI experiment[END_REF][START_REF] Makni | Joint detection-estimation of brain activity in functional MRI: a Multichannel Deconvolution solution[END_REF] and regularized methods [START_REF] Goutte | Modeling the haemodynamic response in fMRI using smooth FIR filters[END_REF][START_REF] Casanova | The impact of temporal regularization on estimates of the BOLD hemodynamic response function: a comparative analysis[END_REF].

Other models such as spatial models [START_REF] Vincent | Spatially adaptive mixture modeling for analysis of fMRI time series[END_REF], and multi-subject methods [START_REF] Zhang | Nonparametric inference of the hemodynamic response using multi-subject fMRI data[END_REF][START_REF] Zhang | A semi-parametric model of the hemodynamic response for multi-subject fMRI data[END_REF] that adaptively learn the HRF across several subjects are outside the scope of this work. The latter models are more relevant in the case of standard group studies and second level analysis.

Our first dataset consists of an encoding study and revealed that it is possible to boost the encoding score by appropriately modeling the HRF. We used two different metrics to assess the quality of our estimates.

The first metric is the fraction of correctly identified images by an encoding model. For this we computed the activation coefficients on both the training and validation dataset. We then learned a predictive model of the activation coefficients from the stimuli. This was used to identify a novel image from a set of 120 potential images from which the activation coefficients were previously computed. The benefits range from 0.9% points to 8.2% points across R1-constrained methods and subjects. The best-performing model in this task is the R1-GLM with FIR basis. The second metric is the Pearson correlation. By considering the voxel-wise score on a full brain volume we observed that the increase in performance obtained by estimating the HRF was not homogeneous across voxels and more important for voxels that already exhibited a good score with a classical design (GLM) and a fixed HRF. The best-performing method is the Rank-1 with separate designs (R1-GLMS) and FIR basis model, providing a significant improvement over the second best-performing model. We also found substantial variability of the shape in the estimated HRF within a single subject and a single task.

The second dataset consists of a decoding task and the results confirmed that constrained (rank-1) estimation of the HRF also increased the decoding score of a classifier. The metric here is Kendall tau.

However, in this case the best performing basis was no longer FIR basis consisting of ten elements but the three elements 3HRF basis (HRF and derivatives) instead, which can be explained by factors such as differences in acquisition parameters, signal-to-noise ratio or by the regions involved in the task.

A higher performance increase was observed when considering the correlation score within the encoding model. This higher sensitivity to a correct (or incorrect) estimation of the HRF can be explained by the fact that the estimation of the HRF is used to generate the BOLD signal on the test set. The metric is the correlation between the generated signal and the BOLD signal. It is thus natural to expect that a correct estimation of the HRF has a higher impact on the results.

In the decoding setup, activation coefficients (beta-map) are computed but the evaluation metric is the accuracy at predicting the stimulus type. The validation metric used for decoding is less sensitive to the HRF estimation procedure than the correlation metric from the encoding study, although it allowed us to observe a statistically significant improvement.

Conclusion

We have presented a method for the joint estimation of HRF and activation coefficients within the GLM framework. Based on ideas from previous literature [START_REF] Makni | Bayesian deconvolution of fMRI data using bilinear dynamical systems[END_REF][START_REF] Vincent | Spatially adaptive mixture modeling for analysis of fMRI time series[END_REF] we assume the HRF to be equal across conditions but variable across voxels. Unlike previous work, we cast our model as an optimization problem and propose an efficient algorithm based on quasi-Newton methods. We also extend this approach to the setting of GLM with separate designs.

We quantify the improvement in terms of generalization score in both encoding and decoding settings.

Our results show that the rank-1 constrained method (R1-GLM and R1-GLMS) outperforms competing methods in both encoding and decoding settings.

  , we performed prediction of BOLD signal following the visual presentation of natural images and compared it against the measured fMRI BOLD signal. As the procedure consists of predicting the fMRI data from stimuli descriptors, it is an encoding model. This dataset is publicly available from http://crcns.org Two subjects viewed 1750 training images, each presented twice, and 120 validation images, each presented 10 times, while fixating a central cross. Images were flashed 3 times per second (200 ms on-off-onoff-on) for one second every 4 seconds, leading to a rapid event-related design. The data were acquired in 5 scanner sessions on 5 different days, each comprising 5 runs of 70 training images -each image being presented twice within the run-and 2 runs of validation images showing 12 images, 10 times each. The images were recorded from the occipital cortex at a spatial resolution of 2mm×2mm×2.5mm and a temporal resolution of 1 second. Every brain volume for each subject has been aligned to the first volume of the first

  with a polynomial of degree 4 and a window length of 91 TR. The activation coefficients (beta-map) and HRF were extracted from the training set by means of the different methods we would like to compare. The training set consisted of 80% of the original session (4 out of 5 runs). This resulted in estimated coefficients (beta-map) for each of the 70 × 4 images in the training set. We proceed to train the encoding model. The stimuli are handled as local image contrasts, that are represented by spatially smoothed Gabor pyramid transform modulus with 2 orientations and 4 scales.

Average

  = *< 0.05, **< 10 -3 , ***< 10 -6 = *< 0.05, **< 10 -3 , ***< 10 -6

Figure 3 :

 3 Figure 3: Top: HRF estimated by the R1-GLMS method on voxels for which the encoding score was above the mean encoding score (first dataset), color coded according to the time to peak of the estimated HRFs. The difference in the estimated HRFs suggests a substantial variability at the voxel level within a single subject and a single task. Bottom: voxel-wise encoding score for the best performing method (R1-GLMS with FIR basis) versus a standard GLM (GLM with fixed HRF) across voxels.

Figure 4 :Figure 5 :

 45 Figure 4: Voxel-wise encoding score for different models that perform HRF estimation (first dataset). As in figure 3, color codes for the time to peak of the estimated HRF at the given voxel. Top: two Rank-1 separate design models with different basis functions: FIR with 20 elements in the Y-axis and the reference HRF with its time and dispersion derivatives (3HRF)in the X-axis. The color trend in this plot suggests that the score improvement of the FIR basis with respect to the 3HRF becomes more pronounced as the time-to-peak of the estimated HRF deviates from the reference HRF (peak at 5s). This can be explained by taking into account that the 3HRF basis is a local model of the HRF around the peak time of the canonical HRF. Bottom: voxel-wise encoding score for two Rank-1 models with FIR basis and different design matrices: separate design on the Y-axis and classical design on the X-axis. Although both models give similar results, a Wilcoxon signed-rank test on the leave-one-session-out cross-validation score (averaged across voxels) confirmed the superiority of the separate designs model in this dataset with p-value < 10 -3 .

https://pypi.python.org/pypi/hrf estimation

https://github.com/neurospin/pypreprocess

Acknowledgements This work was supported by grants IRMGroup ANR-10-BLAN-0126-02 and Brain-Pedia ANR-10-JCJC 1408-01. We would like to thank our colleagues Ronald Phlypo and Gael Varoquaux for fruitful discussions.