Background

Cocaine addiction is a global health issue with limited therapeutic options and a high relapse rate. Attentional bias (AB) towards substance-related cues may be an important factor of relapse. However, it has never been compared in former and current cocaine-dependent patients.

Methods

AB towards cocaine-related words was assessed using an emotional Stroop task in cocaine-dependent patients (CD, N=40), long-term abstinent former cocainedependent patients (ExCD, N = 24; mean abstinence: 2 years) and control subjects (N = 28). Participants had to name the colour of cocaine-related words, neutral words, and colour names. We assessed response times using an automatic voice onset detection method we developed, and we measured AB as the difference in response times between cocaine-related and neutral conditions.

Results

There was an overall group effect on AB towards cocaine, but no group effect on the colour Stroop effect. Two-by-two comparison showed a difference in AB 4 between CD and controls, while ExCD were not different from either. While CD showed a significant AB, consistent with the literature, neither ExCD nor controls showed a significant AB towards cocaine related words. We found no link between AB size and either addiction severity or craving.

Conclusions

Cocaine abstinence was associated with an absence of significant AB towards cocaine-related words which may be interpreted either as absence of AB predicting success in maintaining abstinence, or as AB being able to disappear with long-term cocaine abstinence. Further research is needed in order to distinguish the role of AB in maintaining abstinence.

Introduction

Cocaine addiction is a clinical condition characterised by an excessive intake of cocaine, and a relapsing cycle of intoxication, binging, withdrawal and craving.

Cocaine is the second most used illegal drug in Europe (EMCDDA, 2017) and in France, and a high rate of occasional users fall into cocaine addiction. Cocaine addiction is characterised by the severity of its medical and social consequences [START_REF] Whiteford | Global burden of disease attributable to mental and substance use disorders: Findings from the Global Burden of Disease Study 2010[END_REF]. Indeed, cocaine-dependent patients have a standardised mortality rate 4 to 8 times higher than their group of age peers [START_REF] Degenhardt | Mortality among cocaine users: A systematic review of cohort studies[END_REF]. However, there is to date no available substitution treatment for cocaine addiction [START_REF] Castells | Psychostimulant drugs for cocaine dependence[END_REF], and relapse rates are very high [START_REF] Paliwal | Craving predicts time to cocaine relapse: Further validation of the Now and Brief versions of the cocaine craving questionnaire[END_REF].

One important neurobiological process involved in cocaine addiction is dopaminergic arousal of corticostriatal circuits, and in particular the reward system, which leads to a motor preparation and hyperattentive state towards drug-related cues [START_REF] Franken | Drug craving and addiction: Integrating psychological and neuropsychopharmacological approaches[END_REF]. This arousal could play a causal role in key features of addiction such as drug use despite negative consequences and relapse [START_REF] Pascoli | Sufficiency of Mesolimbic Dopamine Neuron Stimulation for the Progression to Addiction Article Sufficiency of Mesolimbic Dopamine Neuron Stimulation for the Progression to Addiction[END_REF]. Among the cognitive consequences of this arousal is the development of an attentional bias (AB) towards cocaine, which may play an important role in the success or failure of cocaine abstinence. AB is hypothesised to interact with craving, the overwhelming urge to consume the substance, through a loop of mutual reinforcement: the drug cues that patients notice because of AB could heighten their craving, and craving heighten AB [START_REF] Field | Attentional bias in addictive behaviors: A review of its development, causes, and consequences[END_REF].

While the importance of this link has been called into question [START_REF] Field | A Meta-Analytic Investigation of the Relationship Between Attentional Bias and Subjective Craving in Substance Abuse[END_REF], understanding the role of AB in abstinence maintenance would likely allow getting a better picture of the cognitive factors implicated in relapse.

Emotional Stroop tasks are widely used neuropsychological tools for the study of AB in addiction [START_REF] Field | A Meta-Analytic Investigation of the Relationship Between Attentional Bias and Subjective Craving in Substance Abuse[END_REF]. They are based on the classic colour Stroop task, which is used to measure selective attention [START_REF] Wright | What Stroop tasks can tell us about selective attention from childhood to adulthood[END_REF]. In classic Stroop tasks, participants are instructed to name the colour in which colour names are written.

In the cases where a colour name is written in another colour (ie "blue" written in green), this incongruence is associated with a slowing of the colour naming response, which is thought to be caused by the interference between the automatic response (reading the word) and the correct response (naming the ink colour) [START_REF] Macleod | Half a century of research on the Stroop effect: an integrative review[END_REF].

In emotional Stroop tasks, the interference is caused by the presence of emotionally salient words, and the slowing is thus caused by the attention towards the semantic content of these words. The exact mechanisms underlying the emotional Stroop effect are not yet fully understood [START_REF] Cox | The addiction-stroop test: Theoretical considerations and procedural recommendations[END_REF] but the strength of this effect is considered a measure of AB towards a specific semantic category of stimuli, such as a substance in the case of addiction.

Cocaine-related Stroop tasks have shown that cocaine-dependent patients exhibit a higher AB towards cocaine-related stimuli than control participants [START_REF] Copersino | Cocaine craving and attentional bias in cocainedependent schizophrenic patients[END_REF][START_REF] Franken | Selective cognitive processing of drug cues in heroin dependence[END_REF][START_REF] Hester | A consistent attentional bias for drug-related material in active cocaine users across word and picture versions of the emotional Stroop task[END_REF][START_REF] Rosse | Saccadic distractibility in cocaine dependent patients: A preliminary laboratory exploration of the cocaine-OCD hypothesis[END_REF].

Moreover, AB intensity has been linked to both addiction severity and abstinence maintenance in several substance addictions. Higher AB towards the substance is linked to quantity and frequency of substance use [START_REF] Field | Attentional bias in addictive behaviors: A review of its development, causes, and consequences[END_REF], and a higher AB at entry of an inpatient care program has been linked to greater chances of relapse in prospective studies on heroin- [START_REF] Marissen | Attentional bias predicts heroin relapse following treatment[END_REF] and alcoholdependent patients [START_REF] Cox | Alcohol attentional bias as a predictor of alcohol abusers' treatment outcome[END_REF]. Patients who have a documented abstinence for several weeks display lower AB towards heroin [START_REF] Gardini | Decreased drug-cue-induced attentional bias in individuals with treated and untreated drug dependence[END_REF], and alcohol [START_REF] Flaudias | Treatment in hospital for alcohol-dependent patients decreases attentional bias[END_REF] than current users. On the other hand, a similar protocol used to test AB towards smoking words in 24-hour abstinent smokers failed to show a significant difference with ad libitum smokers [START_REF] Munafò | Selective processing of smoking-related cues in current smokers, ex-smokers and never-smokers on the modified Stroop task[END_REF]. However, very few controlled studies have investigated AB among cocaine users and its changes with abstinence. A recent review paper identified only 2 prospective studies on this subject [START_REF] Zhang | A Systematic Review of Attention Biases in Opioid, Cannabis, Stimulant Use Disorders[END_REF]. The first one assessed AB change after 8 weeks of either computer-based cognitive behavioural therapy or counselling. AB was assessed using an emotional Stroop task [START_REF] Devito | Drug Stroop: Mechanisms of response to computerized cognitive behavioral therapy for cocaine dependence in a randomized clinical trial[END_REF], and they found a significant AB towards cocaine-related words pre-treatment, and a drop in AB post-treatment. The second one [START_REF] Mayer | The efficacy of attention bias modification therapy in cocaine use disorders[END_REF] assessed AB, craving and drug use change after 5 sessions of AB modification training (or sham training). AB was assessed through a visual priming task and they found a significant AB pre-treatment and no AB post treatment in both treatment groups.

Because AB is acquired at the onset of cocaine addiction [START_REF] Field | Attentional bias in addictive behaviors: A review of its development, causes, and consequences[END_REF] and drops with short-term abstinence [START_REF] Devito | Drug Stroop: Mechanisms of response to computerized cognitive behavioral therapy for cocaine dependence in a randomized clinical trial[END_REF], one could expect that, as previously demonstrated with heroin and alcohol dependent patients, AB would drop further or disappear in patients with cocaine dependence who maintain a long-term abstinence.

Thus, to further explore this hypothesis, we decided to assess AB towards cocainerelated stimuli with an emotional Stroop task involving cocaine-related words. We compared AB in currently cocaine-dependent patients as well as in formerly cocaine-dependent patients to AB in healthy controls. We recorded addiction severity at the beginning of the task, as it has been reported to influence AB as assessed by the emotional Stroop. We also recorded craving intensity before and after the task.

Methods

Participants

We recruited 40 currently cocaine-dependent patients (CD), 24 formerly cocainedependent patients (ex-CD), and 28 healthy control participants.

Inclusion criteria common to all three groups were being 18 years old or older, being affiliated to the French social security system, and giving informed written consent to participate in the study. Exclusion criteria common to all three groups were colour vision deficit and non-fluency in written or spoken French. Colour vision was assessed based on performance on a training test before the task. All participants spoke and read French fluently.

Patients were recruited among outpatients at an outpatient addiction clinic who had been diagnosed with either current or past cocaine dependence according to the DSM-IV (American Psychiatric [START_REF] Psychiatric Association | Diagnostic and statistical manual of mental disorders: DSM-IV-TR[END_REF]. Diagnosis had been established by their referent psychiatrist at the beginning of treatment through a non-structured clinical interview, and was confirmed by their psychiatrist at the time of the study.

Inclusion criteria for the CD group was current diagnosed cocaine dependence, and declared cocaine use within the past two weeks. Inclusion criteria for the ex-CD group were lifetime diagnosed cocaine dependence, and declared last time of cocaine use being over two months prior.

We chose to rely on self-report for assessing last cocaine use, as it is reliable and easy to implement [START_REF] Brown | Self-reports by alcohol and drug abuse inpatients: factors affecting reliability and validity[END_REF][START_REF] Darke | Self-report among injecting drug users: A review[END_REF]. This self-report was confirmed by patients' referring psychiatrist.

The healthy controls (HC) were recruited through public advertisement on a mailing list of people who volunteer to participate in cognitive science experiments (Risc.cnrs.fr). We excluded participants who had a history of substance abuse (excepting tobacco), neurological or psychiatric disorders through an interview with a psychiatrist. We selected participants on their age, sex, and education level in order to match the patients groups as closely as possible, which somewhat limited us in the number of participants we were able to recruit. The entire process of participating in the study took place in the same day for participants, and it typically lasted between one and two hours.

Two former cocaine addicts withdrew their consent during the task because they felt uncomfortable, although they did not express a rise in their subjective craving.

Task Design

The task was designed using E-Prime 2.0 Standard. We adapted the emotional Stroop test developed for addiction to alcohol by [START_REF] Flaudias | Treatment in hospital for alcohol-dependent patients decreases attentional bias[END_REF]. This tasks consists of three consecutive blocks where participants are asked to name the font colour of words shown on a computer screen. Each block is a different condition: neutral words, colour names and cocaine-related words. We chose not to mix the three conditions but to display stimuli in three distinct, successive blocks, so as to prevent interference that could have been caused by the cocaine-related words [START_REF] Cox | The addiction-stroop test: Theoretical considerations and procedural recommendations[END_REF].

We used two different sets of cocaine-related words, corresponding to the forms of cocaine (cocaine hydrochloride and crack cocaine). Patients were shown the set of words consistent with the form in which they used cocaine, and healthy controls were randomly shown cocaine hydrochloride words.

We chose to use a voice response rather than using button pressing as Flaudias and Llorca (2014) recommend using a vocal response modality for a more natural response and a more pronounced Stroop effect.

Participants were sat about fifty centimetres from a computer screen, asked to focus on a fixation cross, and to name the colour of the words that appear on the screen, regardless of their meaning. All words were randomly displayed in either red, green blue or yellow.

The main variable of interest was the interference caused by the cocaine and colour Stroop effects. This interference was calculated as the difference in reaction times between cocaine and neutral words on the one hand, and colour and neutral words on the other hand.

For reaction times calculations, we considered as usable answers only the correct trials where the first word said by the participant was the correct answer. We therefore excluded from the reaction time calculation trials where participants corrected their answer or started by saying "uh". The number of excluded trials for each group and types of trials can be found in Supplementary Table 1.

For accuracy calculations, were delayed reaction time was not a problem, we included all trials where participants started by giving a correct answer, even if they hesitated before answering.

Choice of words

In order to choose cocaine-related words, we selected potential words with clinicians working with cocaine addicts at the Fernand-Widal hospital. We then showed these words to a group of four cocaine users seeking treatment and asked them to choose the most salient ones and to suggest other words that were not on the list.

The final set of words consisted of four words associated with crack cocaine: "fumer" (to smoke), "pipe", "caillou" (rock) and "crack", and four words associated with cocaine hydrochloride: "sniffer" (to snort), "rail", "ligne" (line) and "coke".

Word frequencies were matched between neutral words and cocaine-related words in order not to overestimate AB towards cocaine. We did so using the Lexique 3.80 lexical database [START_REF] New | Lexique 2: A new French lexical database[END_REF]. We selected: "presser" (to press), "fauteuil" (armchair), "pont" (bridge), "chemise" (shirt). There was no significant difference in frequency between neutral words and cocaine-related words (Kruskal-Wallis χ 2 = 2.58, p-value = 0.28). We were later able to confirm that there was no reaction-time variation between words of the same category (data not shown, available upon request).

Procedure

We did not ask participants to abstain from using cocaine or any other substance prior to the test.

The experiment took place in a quiet room. The task started with two blocks of training. In the first one, participants were presented a series of coloured X (XXXXX) instead of words, for a total of 10 trials.

In the second one, five neutral words were presented twice to participants: "voiture" (car), "livre" (book), "chaussure" (shoe), "route" (road), "chaise" (chair), amounting to 10 trials.

We used four colour names: "bleu" (blue), "rouge" (red), "jaune" (yellow) and "vert" (green), all shown in random, incongruent colours. Participants had three seconds to name each colour, and inter-trial duration was 500 ms.

After the training phase, there were three condition blocks: neutral words, colour names and cocaine-related words were presented in a randomised order. Each word was presented in three different colours, twice for each colour (24 words per condition). Each patient thus named the colour of 92 words during the experiment, including 20 training words.

Patients were shown either cocaine-related words or crack-related words, according to the route of administration that they used most. Healthy controls were shown cocaine-related words.

Patients, but not HC, were asked to rate their craving on a scale from 0 to 10 both just before and just after the task. [START_REF] Gardini | Decreased drug-cue-induced attentional bias in individuals with treated and untreated drug dependence[END_REF] was the only prior study using a drug Stroop task in former and current cocaine users. Based on their effect size, we expected 30 participants in each group to be sufficient to detect the expected effect with a one-sided test.

Data analysis

Power calculation

We were limited in our recruitment by two factors. First, few former cocaine users continue to attend their visits at the addiction clinic. Second, we chose to select healthy controls of a sex, age, and education level similar to those of patients. We chose to perform analysis when we reached 80% of the recruitment goal for all groups.

Accuracy assessment

Responses were manually assessed by listening to the recorded answers. We made two assessments for each trial: whether the answer was correct (i.e. the participant names the correct colour), and whether it was usable for data analysis (i.e. the first word that the participant says is the correct colour).

Reaction time calculation

Reaction times were calculated from the voice response, using the Seewave package for R [START_REF] Sueur | Equipment Review: Seewave , A Free Modular Tool for Sound Analysis and Synthesis[END_REF]. We defined reaction time as the first time when sound intensity was greater than 15% of the maximum sound intensity for the trial.

We ignored sounds that lasted under 100 ms or over 600 ms. We eliminated trials where the detected response time was under 200 ms or over 2000 ms (4% of trials).

Statistical analysis

Because normality assumptions were not always met, we chose to use nonparametric tests.

Patient and HC characteristics are described with frequencies and percentages, mean and standard deviation or median and range, as appropriate. Difference between groups for these characteristics were assessed using Kruskal-Wallis (variance comparison for independent samples), Wilcoxon rank sum (median comparison test for two independent samples), and chi-square tests as appropriate.

We calculated a raw accuracy score on the 72 trials (excluding training) for each participant and a separate accuracy score by condition (colour, neutral, and cocaine words).

Reaction times to accurate response are presented as means and standard deviation by condition.

Colour Stroop interference was calculated as the difference between the mean reaction time for colour words minus mean reaction time for neutral words.

Emotional Stroop interference was calculated as the difference between the mean reaction time for cocaine words minus reaction time for neutral words.

Both interferences were calculated using only response times for usable trials.

Those four measures were described as means and standard deviations and compared between groups. We used Kruskal-Wallis (variance comparison for independent samples), Wilcoxon rank sum (median comparison test for two independent samples) and Jonckheere-Terpstra (similar to Kruskal-Wallis, but the alternative hypothesis assumes an order relation between distributions) tests. The association between clinical factors and those four measures was tested in the two groups of patients with Spearman's correlation or Wilcoxon's rank sum test as appropriate.

Significance threshold was set at p = .05.

Results

Demographic characteristics

The demographic data we collected in the different populations who took part in the experiment is summarised in Table 1. There was no difference between groups on age (Kruskal-Wallis χ² = 3.6, p = 0.2), sex (Kruskal-Wallis χ² = 5.0, p = 0.08), or cognitive functioning (Kruskal-Wallis χ² = 1.1, p = 0.9).

Variable (median

Main analysis

Group comparison showed an overall significant effect on the cocaine interference (Jonckheere-Terpstra JT = 1223, p = 0.03, increasing), and we found a significant difference between CD and HC (Jonckheere-Terpstra JT = 485, p = 0.03, increasing).

The ex-CD group was not significantly different from either CD or HC (Jonckheere-Terpstra JT < 447, p > 0.2).

The CD group showed a significant slowing in the cocaine condition compared to the neutral condition (Wilcoxon W = 565, p = 0.02), while neither ex-CD group (Wilcoxon W = 207, p = 0.1), nor the HC (Wilcoxon W = 334, p = 0.3) did.

All three groups showed a significant slowing in the colour condition compared to the neutral condition, which corresponds to the colour Stroop effect (for the control group: Wilcoxon W = 203, p = 0.002). Group had no significant effect on the colour interference (Jonckheere-Terpstra JT= 1040, p = 0.4, increasing), and there was no significant difference between any two groups on this colour Stroop

346

Healthy controls

Formerly cocainedependent (Ex-CD) There was a significant group effect on overall accuracy (Kruskal-Wallis χ² = 12.5, p = 0.002): the two cocaine groups were not significantly different regarding overall accuracy (Wilcoxon W = 497, p = 0.8), but HC were significantly more accurate than the two other groups (Wilcoxon W > 473, p < 0.005).

Groups were also different in overall response time (Kruskal-Wallis χ² = 12.0, p = 0.003). Similarly, the cocaine groups were not significantly different regarding overall response times (Wilcoxon W = 409, p = 0.3), but HC were significantly faster when compared to either of the two cocaine groups (Wilcoxon W > 209, p < 0.02).

The colour and emotional Stroop effects for each group are plotted in Figure 1. We performed post-hoc analyses, which are reported in the supplementary material. We found no correlation between interferences and any clinical variables, including cocaine severity, and no effect of the task on cocaine craving.

We did find an unexpected order effect for conditions, which is discussed in the supplementary material.

Discussion

The aim of our study was to investigate AB towards cocaine-related words in both currently (CD) and formerly (ex-CD) cocaine-dependent patients, using an emotional Stroop task. Our main hypothesis was that ex-CD would have a lower AB towards cocaine-related words than CD, or even show no AB towards cocaine.

As expected, we found that CD showed a larger AB towards cocaine than controls.

Consistently with the literature, CD show a significant AB towards cocaine-related words, whereas, consistently with our hypothesis, ex-CD showed a non-significant AB, lower than that shown by CD. The control group, HC, shows no significant AB towards cocaine-related words. Additionally, as expected, all groups show a significant colour Stroop effect, and there is no difference between groups on the size of the colour Stroop effect. Our study is the first one to our knowledge to assess cocaine-dependent patients and ex-patients through both a modified drug Stroop task and a colour Stroop task, allowing us to assert that the AB exhibited by participants is actually cocaine-specific, and not the consequence of a general attentional deficit.

The slowing of the CD group on cocaine-related words is coherent with the literature on AB towards various substances [START_REF] Flaudias | Treatment in hospital for alcohol-dependent patients decreases attentional bias[END_REF][START_REF] Gardini | Decreased drug-cue-induced attentional bias in individuals with treated and untreated drug dependence[END_REF]. The 90 ms (i.e. 9%) slowing that we observed is on par with or stronger than other reports (e.g., 39 ms / 4%, from [START_REF] Devito | Drug Stroop: Mechanisms of response to computerized cognitive behavioral therapy for cocaine dependence in a randomized clinical trial[END_REF].

The absence of AB in ex-CD as measured by the cocaine-related Stroop effect is also coherent with the only published longitudinal study using a similar task in a prospective design, showing that an 8-week treatment (outpatient treatment, either counselling or cognitive behavioural therapy) was associated with a significant AB decrease in cocaine-dependent patients who maintain abstinence [START_REF] Devito | Drug Stroop: Mechanisms of response to computerized cognitive behavioral therapy for cocaine dependence in a randomized clinical trial[END_REF].

However, unlike what has been reported with other substances [START_REF] Field | Attentional bias in addictive behaviors: A review of its development, causes, and consequences[END_REF], we found no link between addiction severity and AB size. As our sample size was limited, replication is necessary to confirm this result.

This AB difference could be explained by two different mechanisms, which cannot be distinguished by our experiment: either (1) cocaine-dependent patients with a lower to absent AB have an easier time maintaining abstinence, or (2) the process of maintaining abstinence causes a drop in AB. Hypothesis (1) is in line with findings that show that AB predicts relapse [START_REF] Marissen | Attentional bias predicts heroin relapse following treatment[END_REF] and that training to lower AB can lead to better treatment outcomes in addiction [START_REF] Fadardi | Reversing the sequence: Reducing alcohol consumption by overcoming alcohol attentional bias[END_REF][START_REF] Schoenmakers | Clinical effectiveness of attentional bias modification training in abstinent alcoholic patients[END_REF] -though this effect is not consistent [START_REF] Christiansen | Less than meets the eye: Reappraising the clinical relevance of attentional bias in addiction[END_REF]. However, the fact that maintaining abstinence is associated to a decrease in AB [START_REF] Devito | Drug Stroop: Mechanisms of response to computerized cognitive behavioral therapy for cocaine dependence in a randomized clinical trial[END_REF] gives weight to hypothesis (2).

However, one possible bias in our result could be the presence of a general slowing or attentional deficit in CD and ex-CD: the slowing measured in cocaine-related words colour naming could actually not be specific to cocaine. But indeed, consistently with the literature [START_REF] Hester | A consistent attentional bias for drug-related material in active cocaine users across word and picture versions of the emotional Stroop task[END_REF], we found no difference between groups in the size of the colour Stroop effect: the measured difference in AB is thus not simply due to a general attentional processing difference. However, AB towards cocaine-related words in CD and ex-CD was also correlated with colour Stroop effect size in both the CD and ex-CD groups, which supports the idea, suggested in the literature [START_REF] Compton | Paying attention to emotion: An fMRI investigation of cognitive and emotional Stroop tasks[END_REF], that colour and emotional Stroop effects could have some common basis, such as the involvement of the medial and dorsolateral prefrontal cortex.

Our study is the first study comparing AB towards cocaine in control participants and both current cocaine-dependent patients and long-term abstinent patients, while eliminating the hypothesis of a non-specific attentional effect. It confirms that AB is lower in former-than in current cocaine-dependent patients. This study is a stepping-stone for the design of future prospective studies investigating the possible disappearance of AB with abstinence.

In order to start disentangling these two mechanisms, we are currently recruiting patients for a longitudinal study that will follow them during a 3-month abstinence attempt.

In addition to discussing results themselves, it is important to note that emotional Stroop tasks, despite their widespread use, have important limitations (Ataya et al., 2012a). Their reliability can drop below acceptable levels in some cases, and although we designed our task to minimise this issue by using vocal responses and separating conditions in successive blocks [START_REF] Field | Internal reliability of measures of substance-related cognitive bias[END_REF], replication with other techniques such as eye-tracking [START_REF] Marks | Fixation time is a sensitive measure of cocaine cue attentional bias: Attentional bias and cocaine[END_REF] could be very useful.

The other issue with emotional Stroop task is their specificity: the interference detected in drug Stroop tasks is likely to be influenced not only by attentional bias, but also by other factors such as inhibitory control or cue reactivity (Ataya et al., 2012b). Nevertheless, as Ataya et al. (2012b) point out, these various factors may all play a role in maintaining abstinence, and thus be interesting to measure when trying to understand the dynamics of treatment success.

Finally, Ex-CD reported lower craving than CD before the task, which is coherent with the fact that they successfully avoid using cocaine [START_REF] Preston | Cocaine craving and use during daily life[END_REF].

Craving was non-significantly lower after the task in both groups: we can thus posit that our task does not heighten craving in participants. Several ex-CD participants were distressed by the task and made remarks about the fact that seeing cocainerelated words was unpleasant, and it was therefore important for us to make sure that our task did not have negative consequences on their craving. This distress was not captured by the craving scale, and could be an interesting object to explore in future research. One direction that could be particularly interesting would be the relationship between AB and the frequency of exposure to cocaine-related, although this frequency will be difficult to assess. Indeed, if these variables are negatively correlated, further research around exposure therapy such as the AB modification therapy discussed by [START_REF] Mayer | The efficacy of attention bias modification therapy in cocaine use disorders[END_REF] could help build a path for assisting cocaine-dependent patients in becoming abstinent.

The main limitation of our study is that the number of participants in the three groups is imbalanced, which is due to the difficulty of recruiting former cocaine users in a care setting, and to the difficulty of recruiting healthy controls matched in sex, age, and education level to patients. This may have reduced the power of the study and increased the risk of type 1 error. It is also important to note that women comprised a smaller percentage of the CD group than of either the Ex-CD or the HC group, and that the CD and Ex-CD groups were using psychotropic medications at a very high rate. These differences could contribute to or mask potential group differences.

We can also acknowledge that we did not control for the time of the last cocaine dose or medication or cigarette consumed by patients or healthy controls before the test. Neither did we record the possible withdrawal symptoms in the current cocaine users. Stricter laboratory conditions could be proposed for further studies.

Conclusion

Control participants and formerly cocaine-dependent patients do not have a significant AB towards cocaine, whereas current cocaine addicts display a specific interference effect when assessed with cocaine-related words. This is not a general interference effect, since the colour Stroop effect is observed with the same effect

  CD and ex-CD were recruited within the framework of the Declaration of Helsinki and the ethical guidelines of the Fernand-Widal hospital for the analysis of data already collected during routine care (authorisation 2014-067 given on 15 January 2015 by the CPP (Comité de Protection des Personnes, French regional ethical research committee) and did not receive monetary compensation for their participation. HC were recruited through a physiopathology study (Ethics Committee approval 2012-A01460-43) and received 25€ as compensation for their participation. They all completed the MoCA (Montreal Cognitive Assessment; Nasreddine et al., 2005), which allows screening for mild cognitive impairment. Data describing the severity of the past or current cocaine use were recorded with an ad-hoc questionnaire (age of first use and last use, products used, usage route, dose per day, frequency of use, date of last use), as well as data regarding the lifetime use of other substances. Cocaine craving was assessed with two tools: the OCCS (Obsessive-Compulsive Cocaine Craving Scale ; Vorspan et al., 2012) was used to assess cocaine craving and its consequences on the life of the patients within the last two weeks, and a visual analogue instant scale was used to assess current subjective craving.

  for overall between-groups comparisons for both interference effects.

Figure 1 :

 1 Figure 1: Emotional (cocaine) Stroop effect and colour Stroop effect for the different groups. Stars above each bar denote significance level of the difference between interference sizes and zero; stars on dashes above

Table 1 : Demographic and clinical data.

 1

			Healthy	Ex-CD	CD	P-value
	[range])		controls (HC)	patients (ExCD)	patients (CD)	Overall comparison	Group comparisons
							HC -CD
							HC -ExCD
							ExCD -CD
	N		28	24	40
	Age		46 [24-	45 [24-	41 [26-	ns
			74]	71]	68]
	% women	36	46	20	ns
	% with higher	64	50	38	ns
	education
	% with normal	59	54	60	ns
	MoCA score (≥26)
	%	tobacco	14	88	95	p<10 -10	p<10 -10
	smokers						p<10 -6
							ns
	% current	0	22	51	p<10 -5	p<10 -4
	alcohol use					p<0.05
	disorder						p<0.05
	Alcohol		0.1 [0;2]	0.1	0.2 [0;30] ns
	intake/day		[0;6.4]

Table 2 : Mean Accuracy and Reaction Time (± standard deviation) for the different groups

 2

	% with sedatives	0	42	63	p<10 -5	p<10 -6
	(benzodiazepines					0.001
	and Z-drugs)					ns

We wish to thank Kristel Piani for her help in recruiting patients. Financial support: This work was supported by the Investissements d'Avenir program managed by the ANR (under reference ANR-11-IDEX-0004-02) and by a grant of the FondaMental Foundation. The programme QuitCoc received funding from the French Ministry of Health (PHRC National 2012 AOM12390) and from ERA-NET. It was promoted by the DRCD (Direction de la Recherche Clinique et du Developpement) of the Assistance Publique-Hôpitaux de Paris.


Benjamin Remy
email: benjamin.remy@cea.fr

François Lanusse
email: francois.lanusse@cea.fr

Zaccharie Ramzi
email: zaccharie.ramzi@inria.fr

Jia Liu
email: jialiu@berkeley.edu

Niall Jeffrey
email: niall.jeffrey@phys.ens.fr

Jean-Luc Starck
email: jstarck@cea.fr

Probabilistic Mapping of Dark Matter by Neural Score Matching

The Dark Matter present in the Large-Scale Structure of the Universe is invisible, but its presence can be inferred through the small gravitational lensing effect it has on the images of far away galaxies. By measuring this lensing effect on a large number of galaxies it is possible to reconstruct maps of the Dark Matter distribution on the sky. This, however, represents an extremely challenging inverse problem due to missing data and noise dominated measurements. In this work, we present a novel methodology for addressing such inverse problems by combining elements of Bayesian statistics, analytic physical theory, and a recent class of Deep Generative Models based on Neural Score Matching. This approach allows to do the following: (1) make full use of analytic cosmological theory to constrain the 2pt statistics of the solution, (2) learn from cosmological simulations any differences between this analytic prior and full simulations, and (3) obtain samples from the full Bayesian posterior of the problem for robust Uncertainty Quantification. We present an application of this methodology on the first deep-learning-assisted Dark Matter map reconstruction of the Hubble Space Telescope COSMOS field.

Introduction

Deep Learning methods have proven extremely efficient at solving a wide range of inverse problems such as deconvolution, denoising, super-resolution, inpainting, etc. As a concrete example, deep learning models are now state-of-the-art in MRI [START_REF] Ramzi | Benchmarking MRI reconstruction neural networks on large public datasets[END_REF]. Yet, to apply these methods in the Physical Sciences, two particular preoccupations always remain: How do we robustly quantify uncertainties on a deep learning solution? How do we combine deep learning with existing and trusted physical models? In this work, we present a methodology allowing us to address both of these questions in a consistent Bayesian framework and apply it to a challenging instance of inverse problems from the field of Cosmology: mapping the distribution of Dark Matter in the Large-Scale Structure of the Universe from its weak gravitational lensing effect. Several factors make this problem a perfect example:

-The mapping problem is ill-posed, there is no unique solution and it cannot be solved without making some assumptions on the solution: we want a method that doesn't return a single solution, but a range of solutions covering all possible Dark Matter maps compatible with the observations and our a priori knowledge.

-Analytic cosmological theory can predict some properties of Dark Matter maps, in particular their two-point functions: we want to use as much as possible analytic modeling to constrain these solutions, i.e. we do not want to use deep learning in the regimes where theory applies.

-Numerical cosmological simulations can sample realistic realisations of Dark Matter maps, in all of their complexity: we want to restrict the use of deep learning to only modeling residuals between analytic theory and full simulations.

Mapping the Invisible Dark Matter through its Gravitational Effects

The field of cosmology is investing multiple-billion dollars in space-and ground-based telescopes, with the aim to unveil the fundamental nature of our Universe: dark energy, Dark Matter, cosmic inflation, and Einstein's theory of gravity. The matter distribution-where and how all the matter is positioned throughout space-holds the key to these puzzles. The main difficulty stems from the fact that visible matter, the stars and galaxies, only accounts for 15% of total matter, and the remaining 85% is in the form of invisible Dark Matter. Under Einstein's theory of gravity, the path of a light ray is bent when it passes by massive objects, an effect similar to that of magnifying glasses (and hence the name "gravitational lensing"). Using distant galaxies as our background image, by measuring the distortion of the galaxy shape, we can infer the total matter distribution.

In practice, our task is to infer the matter distribution map, known as a convergence κ map or simply as a mass map, from measurements of individual galaxy ellipticities . These ellipticities are affected by gravitational lensing and can be used as estimators of the lensing signal. From a Bayesian perspective, the shape-to-mass → κ inference problem can be quantified as,

p(κ| , M) ∝ p( |κ, M) p(κ|M), (1)
where M is the cosmological model that encapsulates all the fundamental physics of interest. The likelihood factor p( |κ, M), which describes how light-rays are bent by gravity, how measurements are affected by noise, and accounts for missing observational data. All these factors contribute to making the problem non-invertible, but are all well-understood. The prior distribution however, p(κ|M), which describes the detailed matter distribution under any model M, lacks a full analytic expression. Cosmological theory can provide an analytic Gaussian prior which is accurate but only on large scales.

While a number of approximate methods have been adopted to quantify this prior term on smaller scales [START_REF] Marshall | Maximum-entropy weak lens reconstruction: improved methods and application to data[END_REF][START_REF] Lanusse | High resolution weak lensing mass mapping combining shear and flexion[END_REF][START_REF] Alsing | Hierarchical cosmic shear power spectrum inference[END_REF][START_REF] Alsing | Cosmological parameters, shear maps and power spectra from cfhtlens using bayesian hierarchical inference[END_REF][START_REF] Jeffrey | Improving weak lensing mass map reconstructions using Gaussian and sparsity priors: application to DES SV[END_REF][START_REF] Price | Sparse bayesian massmapping with uncertainties: local credible intervals[END_REF], at present the only accurate physical model of convergence maps for Dark Matter mapping rely on large-scale numerical simulations (presented in [START_REF] Shirasaki | Denoising weak lensing mass maps with deep learning[END_REF][START_REF] Jeffrey | Deep learning dark matter map reconstructions from DES SV weak lensing data[END_REF]), such as the MassiveNuS simulations used in this work [START_REF] Liu | MassiveNuS: Cosmological massive neutrino simulations[END_REF]. As described in [START_REF] Jeffrey | Deep learning dark matter map reconstructions from DES SV weak lensing data[END_REF], forward simulation lets us draw samples from p(κ|M). Our goal in this work is to sample the full posterior distribution p(κ| , M).

Neural Score Matching for Bayesian Inverse Problems

With the fast development of Deep Generative Models, it has recently become possible to learn high dimensional probability distributions from data, opening the prospect of data-driven deep priors. In particular, generative models with explicit likelihoods, like PixelCNNs [START_REF] Salimans | Pixelcnn++: Improving the pixelcnn with discretized logistic mixture likelihood and other modifications[END_REF] and Normalizing Flow [START_REF] Jimenez | Variational inference with normalizing flows[END_REF], can directly learn, evaluate, and sample any p(x) given enough training examples. However, a significant realization is that in the perspective of using such generative models as part of a Bayesian inverse problem such as Equation 1, i.e. alongside a likelihood term, there is no need for a full normalized probability distribution function p(x). Instead, all we need to sample the posterior with modern gradient-based inference techniques such as Hamiltonian Monte-Carlo [HMC, 13] or Black-Box Variational Inference [START_REF] Hoffman | Stochastic variational inference[END_REF], is to have access to the score ∇ x log p(x). Very recent work [START_REF] Song | Generative Modeling by Estimating Gradients of the Data Distribution[END_REF][START_REF] Song | Improved techniques for training score-based generative models[END_REF] has demonstrated that generative models based on learning this score can easily reach GAN quality samples and are both more scalable and easier to train.

Denoising Score Matching As originally identified in [START_REF] Vincent | A connection between scorematching and denoising autoencoders[END_REF] and [START_REF] Alain | What regularized auto-encoders learn from the data generating distribution[END_REF], the score of a given target distribution P can be modeled using a Denoising Auto-Encoder (DAE), i.e. by introducing an auto-encoding function r θ : R n × R → R n trained to reconstruct a true x ∼ P given a noisy version x = x + n with n ∼ N (0, σ 2 I) under an 2 loss. An optimal denoiser r would then be achieved for:

r (x, σ) = x + σ 2 ∇ x log p σ 2 (x) (2)
where p σ 2 = p * N (0, σ 2 ). In other words, the optimal denoiser is closely related to the score we wish to learn and when the noise variance σ 2 tends to zero, should exactly match the score of the target density. In practice to learn this score efficiently, we adopt the noise-conditional Denoising Score Matching (DSM) technique proposed by [START_REF] Lim | AR-DAE: Towards Unbiased Neural Entropy Gradient Estimation[END_REF], and train a model to minimize the loss:

L DSM = E u∼N (0,I) σ∼N (0,s 2 ) u + σr θ (x + σu, σ) 2 2 (3)
In this expression, note that r θ (x + σu, σ) directly models the score ∇ x log p σ 2 and is conditioned on the level of noise σ used to corrupt the input. In practice, the network is a 3-scale U-net with residual blocks composed of convolutions followed by a batch normalisation.

Annealed Hamiltonian-Monte Carlo Sampling Given the noise-conditional neural scores learned with the procedure described above, it is now possible to use a variety of inference methods to access the Bayesian posterior. In this work, we adopt an annealed HMC procedure which provides an efficient way to obtain independent samples from the target posterior despite the high dimensionality of the problem. We consider a tempered version of a target posterior density:

log p σ 2 (x|y) = log p σ 2 (y|x) + log p σ 2 (x) (4)
where σ 2 plays the role of the inverse temperature found in thermodynamic annealing methods, and we recognise the log p σ 2 (x) term modeled by the score network r θ (., σ). At high temperatures, i.e. large σ, the tempered distribution tends to a wide Gaussian, all modes merge together, and an HMC with a diagonal mass matrix can efficiently navigate the posterior distribution. This distribution is then gradually annealed to low temperatures and parallel independent chains progressively moves towards different points in the target distribution. Once we have reached low temperature, we retrieve the last sample from each parallel chain, yielding a batch of independent samples of the posterior distribution. This procedure is similar to the Annealed Langevin diffusion proposed in [START_REF] Song | Generative Modeling by Estimating Gradients of the Data Distribution[END_REF].

Hybrid Analytic and Neural Score Modeling for Mass-Mapping

The method presented in the previous section was based on learning the full score from simulations. However, some properties of the Dark Matter maps are known from analytic cosmology models. We now present how we can directly combine analytic modeling and deep Neural Score Matching. The two-point statistics of weak lensing fields can be accurately predicted from theory over a range of scales, and this forms the basis of most lensing-based cosmological constraints to date. Only on small scales do the theoretical models start to deviate from reality. Since these two-point properties are so well-known, we do not wish to leave it up to the neural network to learn them, and instead The top row shows from left to right, true map from the MassiveNuS simulations, the mask showing the extent of the COSMOS survey, a reconstructed map using a conventional Kaiser-Squires method, and the median of our posterior samples. The bottom row shows multiple samples from the posterior obtained with our procedure.

want to directly impose this prior knowledge on our solution. This prior on κ can be expressed as a Gaussian with a covariance S which becomes diagonal in Fourier space:

p th (κ) = 1 √ det 2πS exp - 1 2 κ † S -1 κ (5)
Using this Gaussian prior alone in a reconstruction would yield Gaussian constrained realisations [START_REF] Hoffman | Constrained Realizations of Gaussian Fields: A Simple Algorithm[END_REF] or the conventional Wiener filter if the maximum a posteriori solution is the target [START_REF] Lahav | Wiener Reconstruction of All-Sky Galaxy Surveys in Spherical Harmonics[END_REF]. But we know the convergence field is not Gaussian, so we aim to learn the difference between the true distribution of convergence maps and a Gaussian model. In our Score Matching framework, we propose the following decomposition of the score of the full prior p(κ):

∇ κ log p(κ) = ∇ κ log p th (κ) + r θ (κ, ∇ κ log p th (κ)) (6)
where r θ is a score network, tasked with modeling the difference between true and Gaussian scores. This network is trained by Denoising Score Matching, simply by adapting Equation 3 to add the theory score to the network output. The main advantage of this approach is that we decrease our reliance on the neural score model to a minimum, lowering the demand on the model complexity.

Results

To demonstrate our methodology, we present the first deep-learning assisted, probabilistic, Dark Matter map reconstruction of the Hubble Space Telescope (HST) Cosmic Evolution Survey (COS-MOS) field [START_REF] Scoville | Cosmos: Hubble space telescope observations[END_REF], using data described in [START_REF] Schrabback | Evidence of the accelerated expansion of the universe from weak lensing tomography with cosmos[END_REF], and basing our prior on the MassiveNuS suite of simulations [START_REF] Liu | MassiveNuS: Cosmological massive neutrino simulations[END_REF]. The simulations contain 1024 3 particles in a (512 h -1 Mpc) 3 box with periodic boundary conditions. 10,000 mass maps are generated by ray-tracing through the simulation boxes, following the physics of gravitational lensing. We make use of the maps for source galaxies at z = 1, which corresponds to the mean redshift of our COSMOS lensing catalog. Along with these maps, we use the corresponding theoretical power spectrum to build the analytic prior p th (κ) for training the score network.

Validation on simulation We first emulate from MassiveNuS simulations a mock COSMOS lensing catalog. We use the actual distribution of galaxies in that survey to create a binary mask representing where lensing is actually measured, and we mock observations by adding realistic levels of noise matching the COSMOS bright galaxy sample of [START_REF] Schrabback | Evidence of the accelerated expansion of the universe from weak lensing tomography with cosmos[END_REF]. Figure 1 illustrates the results of our As can be seen, the posterior is bi-modal in this region, with the cluster present only in some of the samples.

method on this simulated dataset. The bottom row shows 4 independent samples from the posterior. Note that because the data is very noisy, independent posterior samples can exhibit significant variability, except where a strong signal is present. This is for instance the case of the massive cluster visible in the bottom right corner of the field, visible even in the map recovered by a traditional Kaiser-Squires [START_REF] Kaiser | Mapping the Dark Matter with Weak Gravitational Lensing[END_REF] reconstruction. This structure being highly significant, it therefore appears virtually in all posterior samples. Also note how realistic each posterior sample looks compared to the simulation, demonstrating the quality of our score matching model on these 320x320 pixel images.

COSMOS field reconstruction Finally, we present the application of our method to the real COSMOS field, based on the bright galaxy sample described in [START_REF] Schrabback | Evidence of the accelerated expansion of the universe from weak lensing tomography with cosmos[END_REF]. Not only does our posterior mean map represents the highest quality map ever made of the COSMOS field, we have access to the full posterior for interpreting these results. As an example, Figure 2 on the right shows cutouts from independent posterior samples of the central region of the map where a massive cluster is known to be present. We see that the posterior is bi-modal in this region, with the cluster present in some maps and not in others. Measuring the frequency of the cluster appearing in posterior samples would yield a robust quantification of the significance of this structure.

These results demonstrate the benefits of our proposed methodology, which optimally makes use of our physical models and knowledge (through the known likelihood term, theoretical prior on large scales, numerical simulations for small scale non-gaussian prior) and allows us to access a full posterior distribution in high dimension.

Broader Impact

The methodology described here finds many useful applications outside of cosmology. Of particular interest, the same method can be applied to medical imaging (e.g. MRI) and help physicians assess the significance of particular parts of an image before making a diagnostic. We believe this work does not entail any negative consequences or ethical issues.

Figure 1 :

 1 Figure 1: Reconstruction of Dark Matter maps on simulated Hubble Space Telescope lensing measurements.

Figure 2 :

 2 Figure 2: Dark Matter map reconstruction of HST COSMOS survey. We provide both a random sample of the posterior and the posterior mean, to be compared to a conventional Kaiser-Squires reconstruction. Positions of known X-ray clusters and their redshifts are indicated in white. The rightmost plot shows cutouts from independent posterior samples of the central region of the map where a massive cluster is known to be present.As can be seen, the posterior is bi-modal in this region, with the cluster present only in some of the samples.


Zaccharie Ramzi
email: zaccharie.ramzi@inria.fr

Benjamin Remy
email: benjamin.remy@cea.fr

François Lanusse
email: francois.lanusse@cea.fr

Jean-Luc Starck
email: jean-luc.starck@cea.fr

Philippe Ciuciu
email: philippe.ciuciu@cea.fr

Denoising Score-Matching for Uncertainty Quantification in Inverse Problems

Deep neural networks have proven extremely efficient at solving a wide range of inverse problems, but most often the uncertainty on the solution they provide is hard to quantify. In this work, we propose a generic Bayesian framework for solving inverse problems, in which we limit the use of deep neural networks to learning a prior distribution on the signals to recover. We adopt recent denoising score matching techniques to learn this prior from data, and subsequently use it as part of an annealed Hamiltonian Monte-Carlo scheme to sample the full posterior of image inverse problems. We apply this framework to Magnetic Resonance Image (MRI) reconstruction and illustrate how this approach not only yields high quality reconstructions but can also be used to assess the uncertainty on particular features of a reconstructed image.

Introduction

Generative modeling has enjoyed remarkable success in recent years with models such as Generative Adversarial Networks (GANs) [START_REF] Goodfellow | Advances in neural information processing systems[END_REF] reaching extremely high quality results on complex high resolution images [START_REF] Karras | Analyzing and Improving the Image Quality of StyleGAN[END_REF]]. Yet, GANs are still prone to issues including unstable training and mode collapse, i.e. a lack of diversity in generated images. In addition, GANs only provide a convenient way to sample from the learned distribution, they do not give access to the density function itself. Another class of commonly used models, Variational Auto-Encoders (VAEs) [START_REF] Diederik | Auto-encoding variational bayes[END_REF] have also been able to reach high quality sample resolutions [START_REF] Vahdat | NVAE: A Deep Hierarchical Variational Autoencoder[END_REF], but can only provide a bound on the density function of the model.

In some practical applications of these generative models however, the most desirable aspect is not necessarily to be able to sample from the model, or directly have access to its density functions, but have access to gradients of the log density function, i.e. its score. In particular, in this work we focus on the generic problem of performing Bayesian inference to solve an inverse problem (i.e. deconvolution, denoising, inpainting, etc), using a learned generative model as a Bayesian prior. Such problems can typically be solved by gradient-based inference methods such as Variational Inference (VI), and Langevin or Hamiltonian Monte-Carlo (HMC), but all these methods have in common that they only rely on the score of the generative model.

Recent work [Lim et al., 2020, Song and[START_REF] Song | Generative Modeling by Estimating Gradients of the Data Distribution[END_REF] has investigated various approaches to directly train a deep neural network (coined, a score network) to estimate the score of a density function instead of the density itself. One of the most promising and scalable approach relies on denoising score-matching, where the score network is trained using a denoising loss on an augmented dataset where various amounts of Gaussian noise is added to the data. Not only has this approach already been demonstrated to reach convincing results in score estimation, but [START_REF] Song | Improved Techniques for Training Score-Based Generative Models[END_REF] for instance demonstrated that these score networks can be used to sample high resolution natural images through Langevin dynamics.

In this work, we develop the first application of deep denoising score matching to Magnetic Resonance Image (MRI) reconstruction and its use for Uncertainty Quantification (UQ) in the context of imaging inverse problems.

Related Works

Score Matching Score Matching was introduced by [START_REF] Hyvärinen | Estimation of Non-Normalized Statistical Models by Score Matching[END_REF]. It was later revisited in [START_REF] Vincent | A connection between scorematching and denoising autoencoders[END_REF] and [START_REF] Alain | What regularized auto-encoders learn from the data generating distribution[END_REF] in the context of Deep Learning and Denoising Auto-Encoders (DAE) in particular. Recently, [START_REF] Song | Generative Modeling by Estimating Gradients of the Data Distribution[END_REF] applied these findings to images showing the potential of this approach combined with sampling techniques. They refined their work in [START_REF] Song | Improved Techniques for Training Score-Based Generative Models[END_REF] to take it to high dimension and provide some techniques for hyper-parameter setting.

Plug-and-Play Priors The idea of using a denoising model as a prior for inverse problems solving was introduced by Venkatakrishnan et al. [2013] who used non-learned denoisers to replace the proximity operator in the ADMM algorithm [START_REF] Boyd | Alternating Direction Method of Multipliers[END_REF]. More recently, [START_REF] Meinhardt | Learning Proximal Operators : Using Denoising Networks for Regularizing Inverse Imaging Problems[END_REF] have made use of learned denoiser networks in Plug-and-Play inverse problems solving. Annealed HMC sampling can be seen as a fuzzy-version of these approaches with theoretical grounds.

Deep Uncertainty Quantification in Inverse Problems

The idea of using generative models for inverse problems such as data imputation and denoising was for instance suggested in an early VAE work [START_REF] Jimenez Rezende | Stochastic Backpropagation and Approximate Inference in Deep Generative Models[END_REF]. In this approach, a generative network is first trained on high quality data, and inference is then performed in the latent space of the model typically using VI. The main advantage comes from the reduced dimensionality of the parameter space. Recent examples of this approach include [START_REF] Wu | Conditional Inference in Pre-trained Variational Autoencoders via Cross-coding. 5[END_REF][START_REF] Böhm | Uncertainty Quantification with Generative Models[END_REF], and an application to MRI can for instance be found in [START_REF] Edupuganti | Uncertainty Quantification in Deep MRI Reconstruction[END_REF]. An important aspect of this approach is that it relies on the known likelihood of observations. A different approach was proposed in [START_REF] Adler | Deep Bayesian Inversion[END_REF] based on a conditional Wasserstein Generative Adversarial Network (cWGAN). In this formulation, a generative model is trained to sample high quality images conditioned on degraded observations. At test time, independent samples from the posterior distribution induced by the GAN are obtaind by simply sampling different latent space variables. Similar frameworks based on other conditional generative models have alse been proposed, such as [START_REF] Denker | Conditional Normalizing Flows for Low-Dose Computed Tomography Image Reconstruction[END_REF] using conditional Normalizing Flows, or [START_REF] Tonolini | Variational Inference for Computational Imaging Inverse Problems[END_REF] using conditional VAEs. While this conditional generative model approach allows for fast inference, thanks to ancestral sampling, it is worth noting that none of these models include an explicit data consistency step, nor do they make use of a test-time likelihood. As a direct consequence, for instance, these models cannot be used on data with different noise levels as seen during training.

Deep Denoising Score Matching for Posterior Inference

Our main objective in this work is to perform probabilistic inference over a Bayesian model typically represented as p(x|y) ∝ p(y|x)p(x)

(1) where y are some measurements and the posterior p(x|y) is the distribution of possible solutions x compatible with observations and p(x) a prior knowledge. The problem-specific likelihood p(y|x) encodes the forward process of the model and accounts for observational noise, while the prior p(x) encapsulates any a priori information we have on the solution of the problem. Our goal for UQ is to sample solutions x belonging to that posterior. Multiple inference techniques can be leveraged for sampling from this posterior p(x|y), but for high-dimensional problems modern techniques rely on gradient-based methods, including Variational Inference [START_REF] Hoffman | Stochastic variational inference[END_REF] and Langevin Diffusion, or HMC [START_REF] Radford | MCMC using hamiltonian dynamics[END_REF]. All of these techniques have in common that they only require having access to the score ∇ x log p(x) of the target distribution. We illustrate in particular in Appendix B how Metropolis Hastings-calibrated HMC can be conducted with only knowledge of the score. Two terms will enter in the score of the posterior distribution in Equation 1, the score of the likelihood, and the score of the prior. In many problems, such as in the MRI problem presented later in this work, the likelihood score can be derived analytically, only the score of the prior remains unknown but can be learned from data by score matching.

Deep Denoising Score Matching As originally identified by [START_REF] Vincent | A connection between scorematching and denoising autoencoders[END_REF] and [START_REF] Alain | What regularized auto-encoders learn from the data generating distribution[END_REF], the score of a given target distribution P can be modeled using a DAE, i.e. by introducing an auto-encoding function r : R n × R → R n trained to reconstruct under an 2 loss a true x ∼ P given a noisy version x = x + n with n ∼ N (0, σ 2 I). An optimal denoiser r would then be achieved for:

r (x , σ) = x + σ 2 ∇ x log p σ 2 (x )
(2) where p σ 2 = p * N (0, σ 2 ). In other words, the optimal denoiser is closely related to the score we wish to learn and when the noise variance σ 2 tends to zero, should exactly match the score of the target density. In practice to learn this score efficiently, we adopt the residual noise-conditional denoising score matching technique proposed by [START_REF] Lim | AR-DAE: Towards Unbiased Neural Entropy Gradient Estimation[END_REF], and train a model to minimize the following Denoising Score Matching loss:

L DSM = E x∼P E u∼N (0,I) σs∼N (0,s 2 ) u + σ s r θ (x + σ s u, σ s ) 2 2 (3)
In this formulation, the network r θ (x, σ) is now directly modeling the score ∇ x log p σ 2 (x) of the Gaussian-convolved target distribution.

Annealed Hamiltonian-Monte Carlo Sampling Given the noise-conditional neural scores learned with the procedure described above, it is now possible to use a variety of inference methods to access the Bayesian posterior. In this work, we adopt an annealed HMC procedure which provides an efficient way to obtain parallel batches of independent samples from the target posterior despite the high dimensionality of the problem. This is a sampling procedure closely related to the Annealed Langevin Diffusion proposed in [Song andErmon, 2019, 2020], but benefits from the faster Hamiltonian dynamics and Metropolis-Hastings calibration (see Appendix B).

To build our procedure, we consider a Gaussian-convolved version of our target density:

log p σ 2 (x|y) = log p σ 2 (y|x) + log p σ 2 (x) + cst (4)
where σ 2 plays the role of the inverse temperature found in classical annealing. The likelihood can be obtained analytically, and in the case of a Gaussian likelihood takes the following form:

log p σ 2 (y|x) = - x-f (x) 2 2 2(σ 2 n +σ 2 ) + cst,
where σ 2 n is the noise variance in the measurements. As for the prior term, the noise-conditional score network introduced above already models the score of log p σ 2 . This distribution is gradually annealed to low temperatures and the chain progressively moves towards a point in the target distribution.

Application to Bayesian Inverse Problems

The MRI Problem Magnetic Resonance Imaging (MRI) is a non-invasive modality used to probe soft tissues. Compressed sensing, introduced for MRI by [START_REF] Lustig | Sparse MRI: The Application of Compressed Sensing for Rapid MR Imaging Michael[END_REF] is used to reduce its significant acquisition time, and recently, deep learning approaches(see for example [START_REF] Schlemper | A Deep Cascade of Convolutional Neural Networks for MR Image Reconstruction[END_REF][START_REF] Hammernik | Learning a Variational Network for Reconstruction of Accelerated MRI Data[END_REF][START_REF] Pezzotti | An Adaptive Intelligence Algorithm for Undersampled Knee MRI Reconstruction: Application to the 2019 fastMRI Challenge[END_REF]) have been shown to perform extremely well on the reconstruction problem. The idealized reconstruction problem is usually formalized in the following way, for the single coil setting, with uniform acquisition: where y is the acquired Fourier coefficients, also called the k-space data, M Ω is a mask, F is the classical 2D Fourier transform, x is the real anatomical image, and n ∼ N (0, σ 2 n ) is measurement noise (we set σ n = 0.1 in our experiments). The data we consider is the single coil data with the Proton-Density (PD) contrast from the fastMRI dataset introduced in [START_REF] Zbontar | fastMRI: An Open Dataset and Benchmarks for Accelerated MRI[END_REF]. For comparison with deep learning approaches, we use the state-of-the-art Primal-Dual net enhanced with a U-net for image correction (UPDNet) from [START_REF] Ramzi | Benchmarking MRI reconstruction neural networks on large public datasets[END_REF]. The under-sampling was done retrospectively using an acceleration factor of 4 and a random mask as described for knee images in [START_REF] Zbontar | fastMRI: An Open Dataset and Benchmarks for Accelerated MRI[END_REF].

y = M Ω F x + n (5)

Results

We use a simple residual U-net inspired by [START_REF] Ronneberger | U-Net: Convolutional Networks for Biomedical Image Segmentation[END_REF] with ResNet building blocks from [START_REF] He | Deep residual learning for image recognition[END_REF] as our score network. In order to promote the regularity of the learned scores, we regularize the spectral norm of each convolutional layer, and we find that setting the spectral norm to 2 yields the best results. More details about the network architecture and training can be found in appendix A. We then sample from the posterior following our annealed HMC procedure down to relatively low temperatures, and apply one final denoising step on the last sample from the chain, following the Expected Denoised Sample (EDS) scheme detailed in [START_REF] Jolicoeur-Martineau | Adversarial score matching and improved sampling for image generation[END_REF]. Figure 1 compares samples from the MRI posterior to the ground truth (top left), zero-filled image (top second left), and UPDNet reconstruction (top center). We see that although individual samples carry slightly less details than the neural network reconstruction, confidence in any particular part of the image can be gauged by looking for variability or stability across multiple independent posterior samples. We highlight in particular the red region, where posterior samples show significant variability, indicating that this part of the image is poorly constrained by data. In contrast, a direct neural network reconstruction (top center) does not match the ground truth in that region, and does not provide an estimate of uncertainty which may lead a physician to misinterpret the image.

Conclusions and Discussions

We have presented in this work the first (to the best of our knowledge) instance of a framework for Bayesian inverse problems based on Deep Denoising Score Matching and applied to MRI reconstruction and Uncertainty Quantification. We illustrated the merits of this approach on an MRI example where the ability to sample from the full posterior highlights what features present in a reconstruction are not actually constrained by data.

We also showed that the Deep Denoising Score Matching approach can effortlessly be used for dimensions as high as 320 × 320 × 2. It is therefore much more scalable than GAN-or VAE-based models. We note however that finding an optimal HMC annealing schedule and temperature-adaptive step-size proved difficult. In general we did not manage to lower the temperature below a certain level, prompting us to resort to a denoising step on the last samples from our HMC chains. This is a direction for further research.

A MRI experiments A.1 Network architecture

The network is a 3-scale U-net with residual blocks composed of 3 convolutions followed by a batch normalisation. Each batch normalisation is followed by a ReLU non-linearity except the last one. A projection is used for the input of residual blocks whose number of input channels is not the same as that of output channels. Each scale uses 2 residual blocks for the downsampling path and 2 residual blocks for the upsampling path. Downsampling is performed via average pooling and upsampling is performed via up-convolution (as designed in [START_REF] Ronneberger | U-Net: Convolutional Networks for Biomedical Image Segmentation[END_REF]), in order to avoid checkerboard artefacts. We use the following sequence of number of channels: [32,64,128]. In order to deal with MRI images, that are complex in nature, we concatenate the real and imaginary part of the image at the input of the network, forming effectively a 2-channel 320 × 320 image. For input noise level conditioning, we concatenate a noise standard deviation map to the input and in the lower scale of the network. Following the recommendation of [START_REF] Song | Improved Techniques for Training Score-Based Generative Models[END_REF], we also divide the output of the network by the absolute value of the noise power. Finally, we use Spectral Norm regularisation from [START_REF] Yoshida | Spectral Norm Regularization for Improving the Generalizability of Deep Learning[END_REF] in order to make sure the score does not take inconsistent values in low-density regions. The Spectral Norm regularisation indeed forces the spectral norm (maximum eigenvalue) of each layer to a certain value which in turn lowers the Lipschitz constant of the network, preventing it to blow in unseen regions. The influence of the Spectral Norm has been studied in the case of Plug-and-Play approaches by [START_REF] Ryu | Plug-andplay methods provably converge with properly trained denoisers[END_REF].

A.2 Network training

We use the Adam optimizer for network training, with a learning rate of 10 -4 . We add a white Gaussian noise of variance σ 2 s to the images scaled by a factor of 10 6 . σ s is drawn on-the-fly from a Gaussian distribution of variance s = 50. This means that at training time, the standard deviation of the noise can be negative, following the recommendation of [START_REF] Lim | AR-DAE: Towards Unbiased Neural Entropy Gradient Estimation[END_REF] to go from extrapolation to interpolation.

A.3 HMC procedure

The HMC procedure starts from a zero-filled reconstruction of the image with added white Gaussian noise of variance σ init = 100. For the reduction of noise standard deviation at reduction step, we use a factor γ = 0.995, following the recommendation of [START_REF] Song | Improved Techniques for Training Score-Based Generative Models[END_REF]. We take a step size α dependent of the sampling temperature at step i, α =

A.4 Zoom on Reconstruction

The zoom on this part of the reconstruction shown in Figure 2, illustrates a failure of the Neural network reconstruction. It indeed produces a sharp artifact that might hinder the examiner's conclusions. The bayesian sampling allows us to see that this region is not necessarily very sharp.

B Score-Based Hamiltonian Monte Carlo

We used HMC for sampling from the posterior distribution. The sampling procedure follows the leapfrog integrator presented bellow:

m t+ α 2 = m t + α 2 ∇ x log p(x t |y) x t+α = x t + αM -1 m t+ α 2 m t+α = m t+ α 2 + α 2 ∇ x log p(x t+α |y) (6)
where α is the step size, m the momentum and M is a preconditioning matrix taking into account the space metric, but in our case the identity matrix. The gradient of the log posterior is computed as ∇ x log p(x t |y) = ∇ x log p(y|x t ) + ∇ x log p(y|x t ), using the inverse Bayes theorem in equation 1.

It is worth noting that the discretized HMC with final step size is not guaranteed to converge to the desired target distribution. This is the motivation for the Metropolis-Hastings algorithm [START_REF] Metropolis | Equation of state calculations by fast computing machines[END_REF][START_REF] Hastings | Monte carlo sampling methods using Markov chains and their applications[END_REF], which corrects for the discretization errors with a correction step. Concretely, a new step is proposed using Equation 6, but it is accepted with probability of the form min{1, p(x * )q(x n |x * )/p(x n )q(x * |x n )} where x n is the last accepted point in the chain, x * is the proposed update, p is the target density, and q is the proposal corresponding to taking a step of Equation 6. In the context of a score-only approach, we do not have access to p and therefore cannot compute directly this acceptance probability. We can nonetheless estimate from the scores themselves by using a path integral:

log p(x * ) -log p(x n ) = 1 0 ∇ x log p(t * (x * -x n ) + x n ) • (x * -x n )dt (7)
This integral can be evaluated numerically to arbitrary precision (at the cost of additional score evaluations), but we find excellent agreement with the theoretical value on toy model densities like the two-moons distribution with a simple 4 points Simpson integration rule.

Figure 1 :

 1 Figure 1: Bayesian posterior sampling for MRI reconstruction. The top leftmost image is the ground truth image. The top second to the left image is the zero-filled retrospectively under-sampled image F T y. The top third to the left image is the reconstruction of the under-sampled image by the UPDNet. All the other images are denoised samples from the estimated posterior distribution obtained by a tempered HMC. The zero-filling achieves a PSNR of 25.55 dB, each sample 27.63 dB on average, the mean of the samples 30.04 dB and the neural network 32.15 dB. A zoom of the region in the red square is provided in appendix A.4. An animation of posterior samples is available online .

Figure 2 :

 2 Figure 2: Zoom of Bayesian posterior sampling for MRI reconstruction. The ordering is the same as in Figure 1.

Acknowledgments and Disclosure of Funding

We acknowledge the financial support of the Cross-Disciplinary Program on Numerical Simulation of CEA for the project entitled SILICOSMIC. We also acknowledge the French Institute of development and ressources in scientific computing (IDRIS) for their AI program allowing us to use the Jean Zay supercomputer's GPU partitions.

This illustrates that for gradient-based MCMC samplers, which also includes Metropolis Adjusted Langevin Algorithm, these algorithms can fully be implemented with only knowledge of the score.


Zaccharie Ramzi
email: zaccharie.ramzi@inria.fr

Philippe Ciuciu

Jean-Luc Starck

XPDNet for MRI Reconstruction: an Application to the fastMRI 2020 Brain Challenge

We present a modular cross-domain neural network the XPDNet and its application to the MRI reconstruction task. This approach consists in unrolling the PDHG algorithm as well as learning the acceleration scheme between steps. We also adopt state-of-the-art techniques specific to Deep Learning for MRI reconstruction. At the time of writing, this approach is the best performer in PSNR on the fastMRI leaderboards for both knee and brain at acceleration factor 4.

MRI reconstruction context

Magnetic Resonance Imaging (MRI) is a modality used to probe soft tissues. The acceleration of the examination time can be done using Compressed Sensing (CS) [START_REF] Lustig | Sparse MRI: The Application of Compressed Sensing for Rapid MR Imaging Michael[END_REF] and Parallel Imaging (PI). However, to have access the underlying anatomical object x, one must solve the following idealized inverse problem:

M Ω SF x = y (1)
where M Ω is a mask indicating which points in the Fourier space (also called k-space) are sampled, S are the sensitivity maps for each coil, F is the 2D Fourier transform, and y are the coil measurements in the k-space.

This problem is typically solved with optimization algorithms that introduce a regularisation term to solve its indeterminacy.

One point to note, is that S is not generally given (although it can be sometimes roughly pre-acquired).

In this particular instance, it is to be estimated from the data. A classical technique is for example EsPIRIT [START_REF] Uecker | ESPIRiT-An Eigenvalue Approach to Autocalibrating Parallel MRI: Where SENSE meets GRAPPA HHS Public Access[END_REF].

2 Cross-domain networks

General

The general intuition behind cross-domain network is that we are going to alternate the correction between the image space and the measurements space. The key tool for that is the unrolling of optimisation algorithms introduced in [START_REF] Gregor | Learning Fast Approximations of Sparse Coding[END_REF]].

An illustration of what cross-domain networks generally look like is provided in Figure 1.

Preprint. Under review.

Figure 1: General cross-domain networks architecture. Skip and residual connection are omitted for the sake of clarity. y are the under-sampled measurements, in our case the k-space measurements, Ω is the under-sampling scheme, F is the measurement operator, in our case the Fourier Transform, and x is the recovered solution.

Implementation

General backbone The XPDNet is a particular instance of cross-domain networks. It is inspired by the PDNet introduced in [START_REF] Adler | Learned Primal-Dual Reconstruction[END_REF]] by unrolling the PDHG algorithm [START_REF] Chambolle | A First-Order Primal-Dual Algorithm for Convex Problems with Applications to Imaging[END_REF] and applied to MRI in [START_REF] Ramzi | Benchmarking MRI reconstruction neural networks on large public datasets[END_REF]. In particular, a main feature of the PDNet is its ability to learn the acceleration pattern through the unrolled optimisation steps with a buffer, here of size 5.

Image correction network

The plain Convolutional Neural Network (CNN) is replaced by an Multi-scale Wavelet CNN (MWCNN) [START_REF] Liu | Multi-level Wavelet-CNN for Image Restoration[END_REF], but the code1 allows for it to be any denoiser, hence the X. This network was chosen after benchmarking different networks introduced as denoisers and listed in [START_REF] Tian | Deep Learning on Image Denoising: An Overview[END_REF], embedded in the XPDNet but for single-coil reconstruction. Following the insights of [START_REF] Ramzi | Benchmarking MRI reconstruction neural networks on large public datasets[END_REF] we chose to use a smaller image correction network than that presented in the original paper [START_REF] Liu | Multi-level Wavelet-CNN for Image Restoration[END_REF], in order to afford more unrolled iterations in memory.

Additionally, because we use a small batch size, we removed batch normalization layers from the network.

k-space In this challenge, since the data is multicoil, we do not use a k-space correction network which would be very heavy in memory use. However, following the idea of [START_REF] Sriram | End-to-End Variational Networks for Accelerated MRI Reconstruction[END_REF],

we introduce a refinement network for S, initially estimated from the lower frequencies of the retrospectively under-sampled coil measurements. This sensitivity maps refiner is chosen to be a simple U-net [START_REF] Ronneberger | U-Net: Convolutional Networks for Biomedical Image Segmentation[END_REF] like in [START_REF] Sriram | End-to-End Variational Networks for Accelerated MRI Reconstruction[END_REF].

We therefore have 25 unrolled iterations, an MWCNN that has twice as less filters in each scale, a sensitivity maps refiner smaller than that of [START_REF] Sriram | End-to-End Variational Networks for Accelerated MRI Reconstruction[END_REF] and no k-space correction network.

Training details

The loss used for the network training was a compound loss introduced by [START_REF] Pezzotti | An Adaptive Intelligence Algorithm for Undersampled Knee MRI Reconstruction: Application to the 2019 fastMRI Challenge[END_REF], weighted sum of an L1 loss and the Multi-scale SSIM (MSSIM) [START_REF] Wang | Image Quality Assessment: From Error Visibility to Structural Similarity[END_REF].

The optimizer was the Rectified ADAM (RAdam) [START_REF] Liu | On the Variance of the Adaptive Learning Rate and Beyond[END_REF] with default parameters2 .

The training was carried for 100 epochs, each epoch consisting of a pass through all the volumes, with one slice selected at random at each epoch and scaled by a factor of 10 6 as in [START_REF] Ramzi | Benchmarking MRI reconstruction neural networks on large public datasets[END_REF]. The batch size is 1. Masks offset for the equidistant masks3 are sampled on-the-fly. Training is carried out separately for acceleration factors 4 and 8. The networks are then fine-tuned for each contrast for 10 epochs.

With one V100 GPU, the training lasted 1 week for each acceleration contrast.

Conclusion and Discussion

We managed to gather insights from many different works on computer vision and MRI reconstruction to build a module approach. Currently this approach outperforms all others in PSNR and NMSE for both the multicoil knee and brain tracks at the acceleration factor of 4. Furthermore, the modularity of the current architecture allows to use the newest denoising architectures when they become available.

However, the fact that this approach fails to outperform the others on the SSIM metric is to be investigated in further work.

https://github.com/zaccharieramzi/fastmri-reproducible-benchmark

https://www.tensorflow.org/addons/api_docs/python/tfa/optimizers/RectifiedAdam

To see more about the exact nature of the masks: https://github.com/facebookresearch/fastMRI/ issues/54

Acknowledgments and Disclosure of Funding

We acknowledge the financial support of the Cross-Disciplinary Program on Numerical Simulation of CEA, the French Alternative Energies and Atomic Energy Commission for the project entitled SILICOSMIC.

We also acknowledge the French Institute of development and ressources in scientific computing (IDRIS) for their AI program allowing us to use the Jean Zay supercomputers' GPU partitions.