MD Auguste Herissay

MD Laurent Siproudhis

E Le Balc'h

MD Alexandre Merlini L'heritier

MD Marie Dewitte

MD Timothée Wallenhorst

MD Guillaume Bouguen

Dr Charlène Brochard
email: charlene.brochard@chu-rennes.fr

Eric Le Balc'h Md

COMBINED STRATEGIES FOLLOWING SURGICAL DRAINAGE FOR PERIANAL FISTULISING CROHN'S DISEASE: FAILURE RATES AND PROGNOSTIC FACTORS

Keywords: perianal fistulising Crohn's disease, fistulotomy, immunosuppressants

HAL is

INTRODUCTION

Perianal Crohn's disease (CD) includes non-fistulising(fissures, ulcers and strictures) and fistulising (fistulas and abscesses) lesions. Perianal disease generally is associated with a more aggressive and disabling CD phenotype [START_REF] Ingle | The natural history of perianal Crohn's disease[END_REF][START_REF] Schwartz | The natural history of fistulizing Crohn's disease in Olmsted County, Minnesota[END_REF]. Moreover, patients with perianal CD complain of anal pain, discharge, faecal incontinence and dyschezia; emotional symptoms and impairment of social and sexual functioning significantly alter their quality of life [START_REF] Mahadev | Quality of life in perianal Crohn's disease: what do patients consider important?[END_REF]. Perianal fistulising CD is frequent since it occurs either before or after the diagnosis of CD in 14%-37% of cases [START_REF] Schwartz | The natural history of fistulizing Crohn's disease in Olmsted County, Minnesota[END_REF][START_REF] Tang | Are perineal and luminal fistulas associated in Crohn's disease? A population-based study[END_REF][START_REF] Lapidus | Crohn's disease in Stockholm County during 1990-2001: an epidemiological update[END_REF].

The treatment of perianal fistulising CD remains a challenge in clinical practice. The optimal management of perianal fistulising CD requires the combination of perianal surgeries with systemic treatments. During the past two decades, tumour necrosis factor antagonists (anti-TNFα) have revolutionized the approach to treating fistulising CD that is refractory to standard medications [START_REF] Bouguen | Longterm outcome of perianal fistulizing Crohn's disease treated with infliximab[END_REF][START_REF] Present | Infliximab for the treatment of fistulas in patients with Crohn's disease[END_REF]. At present, anti-TNFα antagonists are recommended as the first line of treatment for perianal fistulising CD [START_REF] Gecse | A global consensus on the classification, diagnosis and multidisciplinary treatment of perianal fistulising Crohn's disease[END_REF]. However, less than half of patients experience sustained remission after one year of follow-up [START_REF] Sands | Infliximab maintenance therapy for fistulizing Crohn's disease[END_REF].

Following the initial drainage, the surgical management of perianal fistulising CD remains complex and debated. Several types of surgeries have been proposed to achieve the closure of the fistula, with a focus on sphincter-sparing techniques (SST) (rectal advancement flap [START_REF] Bessi | Advancement flap procedure in Crohn and non-Crohn perineal fistulas: a simple surgical approach[END_REF], glue, plug [START_REF] Senéjoux | Fistula Plug in Fistulising Ano-Perineal Crohn's Disease: a Randomised Controlled Trial[END_REF], ligation of the intersphincteric fistula (LIFT), radiofrequency, laser, cell injection strains [START_REF] Panés | Long-term Efficacy and Safety of Stem Cell Therapy (Cx601) for Complex Perianal Fistulas in Patients With Crohn's Disease[END_REF], seton removal or fistulotomy [START_REF] Bouchard | How to manage anal ulcerations and anorectal stenosis in Crohn's disease: algorithm-based decision making : French National Working Group Consensus 2018[END_REF]. However, these techniques are associated with failure and side effects.

Accepted Article

The aims of this study were (i) to describe the failure rate of therapeutic interventions after drainage of the fistula tract and (ii) to determine the factors associated with failure to optimize the medico-surgical strategies.

PATIENTS AND METHODS

Patients

All prospective medical data of patients with perianal CD were recorded between January 2005 and August 2018 in the database namely, fondamentum, and were reviewed. The database fondamentum is secured (nominative access code) with a CNIL (Commission nationale de l'informatique et des libertés) declaration (CNIL n°1412467). The database is filled in at consultation, during hospitalization and some items are constrained. The data recorded in a secure database were as follows: sex, height, weight, body mass index (BMI), age and luminal CD phenotype at diagnosis according to the Montreal classification [START_REF] Satsangi | The Montreal classification of inflammatory bowel disease: controversies, consensus, and implications[END_REF], extraintestinal manifestations, smoking habits, past and actual treatments (steroids, immunosuppressants and anti-TNFα) and past surgical history (anal fistula surgery, ileal and/or colonic surgery and stoma).

Immunosuppressants included azathioprine, methotrexate, and 6-mercaptopurine.

For the present study, the inclusion criteria were (i) adult age (i.e., aged over 18 years); (ii) CD diagnosis confirmed on clinical, radiological, endoscopic and/or histologic evidence; (iii) undergone both surgical drainage of anal fistula and SST, seton removal or fistulotomy; and (iv) follow-up of at least 12 months. The exclusion criteria were (i) anal neoplasia, (ii) CD diagnosis uncertain, and (iii) anastomotic fistula or enteroanal fistula. Patients with proctectomy, stoma or permanent loose setons were not analysed (Figure 1).

Accepted Article

The study methodology is detailed in the Supplementary Material. The characteristics of the fistula were assessed on the days of i) surgical drainage(s) and ii) surgery for fistula tract management or iii) seton removal. Anatomical classification of perianal CD according to the Cardiff-Hughes classification (15) was used to classify the tract: supralevator fistula, fistula location (anterior or not/high or low), horseshoe fistula, anovaginal fistula, number of fistula tract(s), inter/trans-or suprasphincteric tract, abscess and associated proctitis. Supralevator fistulas were defined as fistulas with supralevator extension. Complex fistulas were defined as supralevator fistulas, deep trans-sphincteric fistulas, fistulas with abscess or extension with the primary tract anovaginal fistulas or horseshoe fistulas. Rectal inflammation was evaluated by rectoscopy. Additionally, CD activity was assessed with the Harvey-Bradshaw Index (16) (HBI), and perianal CD activity was assessed with the Perianal Disease Activity Index (17) (PDAI).

Therapeutic procedures were recorded between the first surgical drainage and the treatment of the tract, including medical treatments (introduction, optimization, switch of immunosuppressant/biological therapy) and surgeries (fistula drainage, intestinal resection). The Harvey Bradshaw Index (HBI), Cardiff classification and PDAI index were used to identify persistent disease; biologics and immunosuppressants were reviewed and adapted according to tertiary centre consensus and staff. According to the HBI, CD is considered to be in clinical remission with an HBI score below 4, mild to moderately active with an HBI score ranging from 4

to 12 and severely active with an HBI score above 12 [START_REF] Harvey | A simple index of Crohn's-disease activity[END_REF]. The PDAI index and clinical assessment were used to identify fistulas that were insufficiently drained. Collections were surgically drained prior to biologic adaptations.

Therapeutic strategies and follow-up

Fistula tract management (FTM) was defined by the day of SST or fistulotomy or seton removal.

The choice of therapeutic strategy was made by the operator according to the anatomical characteristics of the fistula and the patient's field. SST included rectal flap advancement, ligation of the intersphincteric fistula tract (LIFT), fistula plug, fibrin glue injection, Ovesco clip or radiofrequency ablation.

Accepted Article

The evaluation data in the follow-up were extracted from the clinical database. The following data were recorded: occurrence of an abscess, purulent discharge, external fistulous orifice, and discontinuation of immunosuppressant/biological therapy.

The failure of the procedure (SST or fistulotomy or seton removal) was defined as the occurrence of at least one of the following items: abscess recurrence, purulent discharge from the tract, visible external opening, or drainage surgery [START_REF] Bessi | Advancement flap procedure in Crohn and non-Crohn perineal fistulas: a simple surgical approach[END_REF].

Follow-up was defined as the duration between the day of the visit for fistula tract management and the last visit for patients for whom the procedure failed or the last visit.

Statistical analysis

Data are expressed as the the median and percentile (interquartile range: 25% and 75%) or as the

Ethical considerations

The study protocol was approved by the Comité de Protection des Personnes du Sud-Ouest et Outre-mer III. ClinicalTrials.gov ID: NCT03878498.

RESULTS

Study population

Among 212 CD patients who underwent anal fistulas and/or abscesses drainage over a thirteenyear period, 169 patients underwent SST, fistulotomy or seton removal and were included (Figure 1). The baseline characteristics of the population are depicted in

Management between drainage and fistula tract management

A total of 24/169 patients (14%) experienced at least 2 fistula drainages, including 5 (3%) who underwent at least 3 fistula drainages (including the first drainage) between the first drainage and the fistula tract management. Most patients underwent fewer than 2 surgical procedures for fistula drainage (145/169, 86%).

The drainage of anal fistula induced a change in Crohn's treatments for 137/169 patients (81%):

focusing on biologics, anti TNF alpha was introduced in 84/169 (50%), and immunosuppressants were introduced in 93/169 (55%). Anti-TNF alpha was switched or optimized in 13/122 (10%)

patients between drainage and fistula tract management.

Fistula tract management

The median duration between drainage and fistula tract management was 7. and 5 years, respectively. These rates tended to be higher than those of patients who had a fistulotomy (p=0.075) (Figure 3B) and were comparable to those of patients who had SST (p=0.900) (Figure 3C). The cumulative failure rates according to the surgical option are detailed in Supplementary material Table 2.

Factors associated with treatment failure

Figure 4 and Table 3 illustrate items significantly associated with higher failure rates (univariate analyses; p <0.1), including the absence of ulceration (p=0.082), ileal and/or colonic resection (p=0.060), the need for at least 2 fistula drainages (p=0.0009), anterior fistula (p=0.047),

Accepted Article

supralevator extension fistula tract management (p=0.0001) and the absence/discontinuation of immunosuppressants after fistula tract management (p=0.037). In contrast, Montreal classification, gender, smoking, fistula characteristics at drainage, medical treatments at drainage and medical treatment modifications between drainage and fistula tract management tract management showed no impact.

Anterior fistula (HR = 2.52 [1.13-5.61], p=0.024) or supralevator extension (HR = 20.78 [3.38-127.80] p=0.001) at fistula tract management and the absence/discontinuation of immunosuppressants after fistula tract management (HR = 3.74 [1.11-12.5], p=0.032) were significantly associated with higher failure rate in the multivariate analysis model (Table 3).

DISCUSSION

The present work quantified the rate of failure of fistula tract management and assessed the factors associated with failure in perianal fistulising CD. Approximately one-third of patients failed to maintain a sustained benefit. Analyses showed that anterior fistula, supralevator extension and the absence/discontinuation of immunosuppressants after fistula tract management were significantly associated with a less favourable outcome. Fistulotomy was a therapeutic strategy that significantly increased the sustained healing rate.

The main strengths of this work are the sample size, the systematic assessment of perianal CD in a referral centre for CD and the duration of the follow-up. The data were recorded in a prospective database using recommended classifications and validated scales [START_REF] Harvey | A simple index of Crohn's-disease activity[END_REF][START_REF] Irvine | Usual therapy improves perianal Crohn's disease as measured by a new disease activity index. McMaster IBD Study Group[END_REF][START_REF] Satsangi | The Montreal classification of inflammatory bowel disease: controversies, consensus, and implications[END_REF]. Disease characteristics of CD and perianal CD and both medical and surgical strategies were assessed. All events were prospectively recorded. The main judgement criterion was based on failure because it is easier to evaluate and is absolutely objective.

However, our study results should be interpreted with caution for several reasons. The data were prospectively recorded but the analysis was retrospective. The recruitment of patients might have resulted in selection bias: this study was performed in a tertiary referral centre, the number of patients treated with anti-TNFα was high and fistulas were severe. Magnetic resonance imaging (MRI) and systematic assessment of faecal continence were not done for all consecutive patients so these data could not be analysed. Finally, patients with stoma or permanent loose seton were

Accepted Article

not included since we have analyzed the failure rate of therapeutic interventions after drainage of the fistula tract. The stoma and permanent loose seton may be transient but may also be associated with an unfavourable outcome: it may have led to an underestimation of the failure rate.

One of our objectives was to compare the rate of failure of management of drained fistulas.

Overall, the failure rate was 30% at 3 years, which is in accordance with the data in the literature [START_REF] Bouguen | Longterm outcome of perianal fistulizing Crohn's disease treated with infliximab[END_REF][START_REF] Sands | Infliximab maintenance therapy for fistulizing Crohn's disease[END_REF]). It appears that patients who had an SST or seton removal have a nearly 2-fold risk of failure compared to those who have a fistulotomy. This risk is comparable in the short and long term.

These data are in line with previous reports [START_REF] Williams | Fistula-in-ano in Crohn's disease. Results of aggressive surgical treatment[END_REF][START_REF] Papaconstantinou | Surgical Management of Fistula-in-ano Among Patients With Crohn's Disease: Analysis of Outcomes After Fistulotomy or Seton Placement-Single-Center Experience[END_REF]. Some words of caution are needed. Indeed, fistulotomy was performed only when fistulas were simple. Fistulas can be considered simple, even in the case of CD, in case of single fistulas, unbranched, with low tract and without supralevator extension to avoid continence disorders. If these conditions are met, fistulotomy may be an option with respect to the continence status and risk of incontinence.

Importantly, anatomical features of the anal fistula including anterior fistula and supralevator extension negatively impacted the outcome of perianal fistulising CD. These data are in line with previous studies [START_REF] Brochard | Anal Fistulas in Severe Perineal Crohn's Disease: Mri Assessment in the Determination of Long-Term Healing Rates[END_REF]. Some anatomical features of fistula can be modified by surgical drainage techniques and lowering of the tract. The goals of fistula surgery are to achieve efficient drainage and simplify fistula tracts. Lowering the tract is an important strategy in managing high/complex fistulas. Therefore, effective fistula drainage which may include several interventions is needed to avoid SST failure in particular rectal flap advancement [START_REF] Bessi | Advancement flap procedure in Crohn and non-Crohn perineal fistulas: a simple surgical approach[END_REF]. Furthermore, luminal disease was controlled and no patient had active proctitis. The failure rate of supralevator extension is so high that a complementary strategy including seton removal is questionable. Loose seton could be preferred to avoid abscess recurrence in such a situation.

In contrast to previous works [START_REF] Present | Infliximab for the treatment of fistulas in patients with Crohn's disease[END_REF][START_REF] Sands | Infliximab maintenance therapy for fistulizing Crohn's disease[END_REF], this study does not show an association between introduction or adaptation of anti-TNFα and a better healing rate. One hypothesis may be that the modalities of anti-TNFα treatment were heterogeneous due to the retrospective cohort (no long-term treatment and fewer combination therapies at the beginning of the study period). Another hypothesis may be that most patients underwent both surgical and biologic treatments concomitantly as compared with previous studies. Nevertheless, the absence/discontinuation of immunosuppressants after fistula tract management was significantly associated with a less favourable outcome. These data are in accordance with some studies [START_REF] Schröder | Combining infliximab and methotrexate in fistulizing Crohn's disease resistant or intolerant to azathioprine[END_REF][START_REF] Roumeguère | Combined approach with infliximab, surgery, and methotrexate in severe fistulizing anoperineal Crohn's disease: results from a prospective study[END_REF][START_REF] Topstad | Combined seton placement, infliximab infusion, and maintenance immunosuppressives improve healing rate in fistulizing anorectal Crohn's disease: a single center experience[END_REF] suggesting that immunosuppressants would be

Accepted Article

effective in perianal fistulising CD. Medical treatment remains essential to maintain the healing of perianal fistulising CD; this work put into question whether immunosuppressants should be maintained.

Prospective evaluation of faecal incontinence with validated questionnaires is needed to better quantify the benefit/risk ratio of the surgical strategy with a special focus on fistulotomy.

CONCLUSION

Combined surgical and medical strategies after surgical drainage in perianal fistulising CD lead to failure rate of 36% at 5 years. Fistulotomy is effective in treating perianal fistulising CD in a lowrisk situations of incontinence; this option should not always be ruled out in principle.

Immunosuppressants may have an important role in maintaining control of perianal fistulising CD.

  number and percentage of the cohort. The nonparametric Wilcoxon signed-rank test for continuous variables and nonparametric Pearson or Fisher exact tests for categorical variables were used. The judgement criteria were the occurrence of failure after the fistula tract management. The cumulative probabilities of failure were estimated using the Kaplan-Meier method with 95% confidence intervals [CIs]. To identify factors associated with failure, a univariate analysis was first performed using the log-rank test. The results are shown as hazard ratios (HRs) with 95% confidence intervals [CIs]. Statistical analyses were performed using JMP Pro 13.0.0 software (SAS institute).

Table 1 .

 1 Only 10 patients (6%) BMI>=30 kg/m 2 ). The duration of CD was less than one year at first drainage in 74/169 patients (44%), and 31/169 (18%) had a past history of anorectal stricture or anal ulceration. Most patients received medical treatment before drainage: 62/169 patients (37%) were or had been treated with immunosuppressant monotherapy, 20/169 (12%) were or had been treated with anti-TNFα only and 19 /169 (11%) were or had been treated with combination therapy.

	Accepted Article

Accepted Article

Supplementary material:

(1) Figure : Description of the study methodology

(2) Supplementary material Table 1: Fistula characteristics at fistula tract management according to therapeutic strategy

(3) Supplementary material Table 2: Failure rates according to the surgical option


M J Page

M Symeonidis

J D Vieira

B Altieri

A Amblard

V Arumugam

H Aussel

T Babbedge

A Blain

J Bock

A Boselli

V Buat

N Castro-Rodríguez

A Cava

P Chanial

D L Clements

A Conley

L Conversi

A Cooray

C D Dowell

E N Dubois

J S Dunlop

E Dwek

S Dye

S Eales

D Elbaz

D Farrah

M Fox

A Franceschini

W Gear

J Glenn

M Griffin

M Halpern

E Hatziminaoglou

E Ibar

K Isaak

R J Ivison

G Lagache

L Levenson

N Lu

S Madden

B Maffei

G Mainetti

L Marchetti

H T Nguyen

B O'halloran

S J Oliver

A Omont

P Panuzzo

A Papageorgiou

C P Pearson

I Pérez-Fournon

M Pohlen

J I Rawlings

D Rigopoulou

L Riguccini

D Rizzo

G Rodighiero

I G Roseboom

M Rowan-Robinson

M Sánchez Portal

B Schulz

Douglas Scott

N Seymour

D L Shupe

A J Smith

J A Stevens

M Trichas

K E Tugwell

M Vaccari

I Valtchanov

M Viero

L Vigroux

L Wang

R Ward

G Wright

C K Xu

M Zemcov

The suppression of star formation by powerful active galactic nuclei

The old, red stars which constitute the bulges of galaxies, and the massive black holes at their centres, are the relics of a period in cosmic history when galaxies formed stars at remarkable rates and active galactic nuclei (AGN) shone brightly from accretion onto black holes. It is widely suspected, but unproven, that the tight correlation in mass of the black hole and stellar components 1 results from the AGN quenching the surrounding star formation as it approaches its peak luminosity 2;3;4 . X-rays trace emission from AGN unambiguously [START_REF] Brandt | Deep Extragalactic X-Ray Surveys[END_REF] , while powerful star-forming galaxies are usually dust-obscured and are brightest at infrared-submillimetre wavelengths [START_REF] Sanders | Luminous Infrared Galaxies[END_REF] . Here we report observations in the submillimetre and X-ray which show that rapid star formation was common in the host galaxies of AGN when the Universe was 2-6 Gyrs old, but that the most vigorous star formation is not observed around black holes above an X-ray luminosity of 10 44 erg s -1 . This suppression of star formation in the host galaxies of powerful AGN is a key prediction of models in which the AGN drives a powerful outflow 7;8;9 , expelling the interstellar medium of its host galaxy and transforming the galaxy's properties in a brief period of cosmic time.

Measuring star formation in galaxies containing powerful active galactic nuclei has long been a problem, because the radiation from the AGN outshines that from star formation in almost all wavebands. Of all parts of the electromagnetic spectrum, the far-infrared to millimetre waveband offers the best opportunity to measure star formation in galaxies hosting AGN because, in contrast to strongly star-forming galaxies, AGN emit comparatively little radiation at these wavelengths [START_REF] Hatziminaoglou | HerMES: Far infrared properties of known AGN in the HerMES fields[END_REF] . The combination of deep X-ray and submillimetre observations therefore offers the best prospects for studying the association of star formation and accretion during the 1 < z < 3 epoch when star formation and black hole growth in massive galaxies were at their most vigorous.

The X-ray catalogue of the Chandra Deep Field North (hereafter CDF-N) derives from a series of observations made with the Chandra X-ray observatory with a total of 2Ms exposure time [START_REF] Alexander | The Chandra Deep Field North Survey. XIII. 2 Ms Point-Source Catalogs[END_REF] . We restrict the sample to those sources detected in the most penetrating 2-8 keV band to minimize the influence of obscuration on our results, and we further limit the sample to those sources (64%) for which spectroscopic redshifts are available in the literature 12;13 . 2-8 keV luminosities were calculated assuming that AGN X-ray spectra are power laws of the form 14 S ν ∝ ν -0.9 ; the luminosities are not corrected for absorption intrinsic to the AGN or their host galaxies. In order to restrict the X-ray sample to AGN we have discarded any sources with 2-8 keV luminosity L X < 10 42 erg s -1 . SPIRE 15 observations of the CDF-N were carried out in October 2009 as part of the HerMES programme [START_REF] Oliver | HerMES: SPIRE galaxy number counts at 250, 350, and 500 µm[END_REF] . Maps and source catalogues at 250, 350 and 500 µm were constructed [START_REF] Smith | HerMES: point source catalogues from deep Herschel-SPIRE observations[END_REF] . At the depth of the SPIRE maps, the dominant source of uncertainty in the maps is confusion noise due to the high sky density of sources. For cross-matching with the Chandra source catalogue [START_REF] Alexander | The Chandra Deep Field North Survey. XIII. 2 Ms Point-Source Catalogs[END_REF] we chose the 250 µm catalogue, which has the most precise positions, and we used only sources with 250 µm flux densities > 18 mJy, which corresponds to a 1. The expected level of spurious associations between X-ray and 250 µm sources was calculated from the sky density of 250 µm sources in 10 to 30 arcsecond-radius annular regions around the X-ray source positions, and is reported in Table 1.

The distribution of CDF-N AGN in the redshift -luminosity plane is shown in Figure 1, and shows a striking trend of 250 µm detectability with X-ray luminosity: of the 24 AGN detected at 250 µm, none of them have L X > 10 44 erg s -1 . The redshift range between 1 and 3 is of most interest, because it corresponds to the epoch in which powerful AGN accreted most of their black hole mass and present-day massive galaxies formed most of their stars.

Within this redshift range, Fig. 1 shows that 11 out of 44 AGN (25 +8 -7 %) with 10 43 erg s -1 < L X < 10 44 erg s -1 are detected at 250 µm, while none of the 21 objects with L X > 10 44 erg s -1 are detected. The difference in detection rates has a significance of > 99% according to a single-tail Fisher's exact test. We have considered the effects that incompleteness in the spectroscopic redshifts, or absorption of the X-ray flux by gas and dust, might have on our results. We find that the systematic non-detection of the powerful AGN is robust against both effects, although X-ray absorption does appear to be a common property of the 250 µm-detected AGN. We have also verified the low 250 µm detection rate of AGN with L X > 10 44 erg s -1 using the Extended Chandra Deep Field South field, finding that of 49 such sources with 1 < z < 3, only 1 is detected at 250 µm.

Infrared spectral energy distributions for the 250 µm-detected AGN were constructed by combining the SPIRE photometry with Spitzer 3.6 -160 µm photometry. X-ray and infrared properties of the 11 250 µm-detected AGN with 1 < z < 3 and L X in the range 10 43 -10 44 erg s -1 are given in Table 2. In most cases the AGN contributes less than 10% to the infrared luminosity. The best-fit infrared luminosities lie between 4 × 10 erg s -1 < L X < 10 44 erg s -1 . The first column gives the ID number of the source in the X-ray catalogue [START_REF] Alexander | The Chandra Deep Field North Survey. XIII. 2 Ms Point-Source Catalogs[END_REF] . The second column gives the redshift, and the third column gives the logarithm of the 2-8 keV X-ray luminosity. The fourth column gives the log of the column density of absorbing gas (units of hydrogen atoms per cm 2 ) implied by the ratio of 2-8 keV to 0.5-2 keV X-rays; a dash indicates no evidence for photoelectric absorption in X-rays, and an asterix is used where the lower limit to the column density is unconstrained. The fifth column gives the correction to log L X to account for the absorption. The sixth column gives the 8-1000 µm IR luminosity. The seventh column gives the maximum likely contribution of an AGN to the infrared luminosity, and the eighth column gives the range of star formation rate implied by the infrared luminosity, where the upper and lower limits correspond to zero AGN contribution and the maximum AGN contribution to the IR luminosity, respectively. Photometry for the spectral energy distributions was extracted from Spitzer and SPIRE images using the X-ray and 24 µm catalogue positions as priors [START_REF] Roseboom | The Herschel Multi-Tiered Extragalactic Survey: source extraction and crossidentifications in confusion-dominated SPIRE images[END_REF] . Total 8-1000 µm infrared luminosities were then determined by fitting templates [START_REF] Siebenmorgen | Dust in starburst nuclei and ULIRGs. SED models for observers[END_REF] to the spectral energy distributions [START_REF] Symeonidis | The link between SCUBA and Spitzer: cold galaxies at z ≤ 1[END_REF] . Upper limits to the AGN contribution to the infrared luminosities were obtained by normalising an AGN template in the mid infrared [START_REF] Seymour | HerMES: SPIRE emission from radio-selected active galactic nuclei[END_REF] .

ID Redshift Log L X (erg s -1 ) Log N H (cm -2 ) N H corr Log L IR (L ⊙ ) AGN (%) SFR (M ⊙ yr -1 ) 35
10 13 L ⊙ , implying star formation rates between 50 and 1750 M ⊙ per year [START_REF] Kennicutt | The global Schmidt law in star-forming Galaxies[END_REF] .

We performed a stacking analysis for the 1 < z < 3 AGN to probe below the confusion limit of the SPIRE images. We split the sample into five bins of L X from 10 43 to 10 45 erg s -1 and determined the average star formation rates of AGN in each bin. The results are shown in Figure 2. In the redshift range 1 < z < 3, the mean star formation rate in AGN with L X of 10 43 -10 44 erg s -1 is 214±25 M ⊙ per year, compared to a mean star formation rate for AGN with L X > 10 44 erg s -1 of 65±18 M ⊙ per year. These averages are independent of the SPIRE 250 µm detection limit because they are obtained from a stack of all sources within a given range of L X , whether detected at 250 µm or not.

At redshifts of 1-3, the X-ray luminosity of 10 44 erg s -1 , which divides the regions of 250 µm detection and non-detection in Fig. 1, corresponds approximately to the knee in the luminosity function of AGN [START_REF] Ebrero | The XMM-Newton serendipitous survey. VI. The X-ray luminosity function[END_REF] . The steep shape of the luminosity function at L X > 10 44 erg s -1 implies that this part of the luminosity function is dominated by objects which are at the peak of their accretion rates. Our observations therefore imply that the most prodigious episodes of star formation are common in the host galaxies of 1 < z < 3 AGN, but avoid powerful AGN in which accretion is at its peak.

Figure 2: Star formation rates derived from averaged far-IR luminosities of 1 < z < 3 AGN. We converted the 250, 350 and 500 µm flux densities for each source into an equivalent 8-1000µm luminosity by fitting a grey-body curve, with a temperature of 30 K in the rest-frame of the source, an emissivity index of β = 1.6, and a powerlaw extension to the Wien side 29 and multiplying by 4πD 2 L where D L is the luminosity distance. Fluctuations in the map sometimes scatter the fluxes of undetected sources to negative values, which translate to negative luminosities when multiplied by 4πD 2 L . Such negative solutions for individual AGN were allowed so as not to produce an artificial positive bias in the averages. The luminosities were averaged in 5 bins in L X , which were chosen to include a similar number of AGN in each bin. The average luminosities were then converted to star formation rates [START_REF] Kennicutt | The global Schmidt law in star-forming Galaxies[END_REF] . AGN which are individually detected at 250 µm are included in the averages shown in bold black, but have been excluded from the averages which are shown in grey, to show the contribution that these sources make to the average star formation rates. The grey points have been offset horizontally from the bold black points for clarity. Error bars correspond to 68% confidence and were determined by bootstrap resampling, with a 7% systematic error added in quadrature to account for the calibration error on SPIRE photometry. This systematic non-coincidence of the peak periods of star formation and accretion implies a palpable interaction between the two processes, and provides a powerful discriminator for the form of AGN feedback which is responsible for terminating star formation in the host galaxy. Two families of feedback models have been proposed, widely referred to as 'quasar mode' and 'radio mode' [START_REF] Croton | The many lives of active galactic nuclei: cooling flows, black holes and the luminosities and colours of galaxies[END_REF] . In quasar-mode feedback, a luminous AGN generates a powerful wind which terminates star formation by driving the interstellar medium from the surrounding host galaxy. In radio-mode feedback, star formation is suppressed because collimated jets of relativistic particles emitted by a radiatively-inefficient AGN prevent gas in the surrounding hot halo from cooling, thereby starving the galaxy of cool gas from which to form stars.

Radio-mode feedback is commonly invoked in semi-analytical models to limit galaxy masses and luminosities 21;20 . In these models, black holes grow through luminous accretion episodes and black hole mergers. The correlation between black hole and bulge mass comes from assuming that a fixed fraction of the gas is accreted by the nucleus during each star forming episode that results from a galaxy merger or disc instability, and hence star formation and accretion rate should be correlated over the full range of luminosity. Our observations are therefore inconsistent with models in which AGN influence their host galaxies only through radio mode feedback 21;20 . In contrast, models of galaxy formation in which quasar-mode feedback is responsible for terminating the star formation 22;8;7;9 , and which have received some observational support recently 23;24 , predict that the AGN luminosity peaks later than the star formation rate, and thus are consistent with our observations. These models also predict that residual star formation, at the level of a few tens of per cent of the peak, will continue during the period in which the AGN luminosity is at its maximum, consistent with our stacked results which show that on average AGN with L X > 10 44 erg s -1 are still forming stars at approximately 65 M ⊙ per year.

Our observations do not discriminate between models invoking major mergers 8 or accretion of gas into a massive halo [START_REF] Granato | A physical model for the coevolution of QSOs and their spheroidal hosts[END_REF] as the trigger for the intense star formation. After the interstellar medium has been driven out by the luminous AGN and the AGN itself becomes starved of fuel, radio-mode feedback is the most credible agent by which further star formation is inhibited.

Source number in X-ray catalogue [START_REF] Trouille | The OPTX Project. I. The Flux and Redshift Catalogs for the CLANS, CLASXS, and CDF-N Fields[END_REF] S1: Spectroscopic redshifts which do not come from our main source [START_REF] Barger | A Highly Complete Spectroscopic Survey of the GOODS-N Field[END_REF] and references to the literature in which these redshifts are published.

221, a photometric redshift of 1.9 is found elsewhere in the literature [START_REF] Donley | Spitzer Power-Law Active Galactic Nucleus Candidates in the Chandra Deep Field-North[END_REF] , in which case it has L X < 10 44 erg s -1 in keeping with the other 250 µm-detected AGN. Source 107 is within 6 arcsec of the same 250 µm source as a brighter X-ray source (source 109), which has a spectroscopic redshift of 2.58. The brighter X-ray source is closer to the 250 µm source position than source 107, and source extraction using X-ray and 24 µm positions as priors [START_REF] Roseboom | The Herschel Multi-Tiered Extragalactic Survey: source extraction and crossidentifications in confusion-dominated SPIRE images[END_REF] does not recover significant 250 µm emission for source 107, leaving some doubt about the reality of the association. Source 107 has a photometric redshift of 1.15, but could be physically associated with source 109 at z = 2.58 given that they are only separated by 2.7 arcsec. It can be seen from Fig. S1 that whichever of the two redshifts is adopted for source 107, it has L X between 10 43 erg s -1 and 10 44 erg s -1 , and hence if the association of source 107 and the 250 µm source is real, it too has L X < 10 44 erg s -1 , in keeping with the other 250 µmdetected AGN. Finally, we come to source 406, which is the brightest of the six X-ray sources. From the track in Fig. S1, it appears the most likely of the six to be found within the 1 < z < 3, 10 44 erg s -1 < L X < 10 45 erg s -1 region of our survey. This source is extremely faint in the optical, with an R-band magnitude of 27.8 [START_REF] Barger | Optical and Infrared Properties of the 2 Ms Chandra Deep Field North X-Ray Sources[END_REF] , and there are no estimates of photometric redshift available in the literature. The association of this X-ray source with the 250 µm source must be regarded with some doubt, however, because it is located just 5.5 arcsec to the south of a brighter 24 µm and radio source which is closer to the 250 µm position. Source extraction using -9 ) % 1 < z < 3, 10 44 < L X < 10 45 erg s -1 25 0 0.3 < 4 (< 11) % Table S2: 250 µm detection statistics in the various parts of parameter space when AGN with photometric redshifts are included as well as those with spectroscopic redshifts. The first row corresponds to the full AGN luminosity and redshift parameter space in the CDF-N, while the second and third rows correspond to the regions enclosed within the blue dashed lines in Fig. 1.

X-ray and 24 µm positions as priors [START_REF] Roseboom | The Herschel Multi-Tiered Extragalactic Survey: source extraction and crossidentifications in confusion-dominated SPIRE images[END_REF] detects source 406 only at the 3σ level. Nonetheless, this source does represent the single example of an X-ray source, potentially associated with a 250 µm source, with a reasonable likelihood of lying within the 1 < z < 3, L X > 10 44 erg s -1 region. Equally, it may fall within the 1 < z < 3, L X < 10 44 erg s -1 region. To summarise, it appears that the sources without spectroscopic redshifts are likely to harbour at most one object within the redshift range 1 < z < 3, which has L X > 10 44 erg s -1 .

Effect of X-ray absorption

X-rays in the 2-8 keV band can penetrate considerable column densities of material, but column densities > 10 23 cm -2 can lead to significant attenuation of the X-ray flux. Conceivably, the different 250 µm detection rates for AGN above and below 10 44 erg s -1 could be a consequence of X-ray absorption, if 250 µm-luminous AGN are systematically more heavily absorbed than AGN which are not detected at 250 µm. X-rays in the 0.5-2

keV band are absorbed more easily than those in the 2-8 keV band which was used to select our X-ray sample, and hence the hardness ratio of 2-8 keV / 0.5-2 keV X-rays is a crude but effective estimator of the degree of X-ray absorption.

For each 2-8 keV X-ray source, we used the Portable Interactive Multi-Mission Simulator (PIMMS; http://heasarc.nasa.gov/docs/software/tools/pimms.html) to convert from hardness ratios to spectral indices and column densites. The intrinsic spectra of faint AGN without absorption are power laws of the form S ν ∝ ν -α with a mean α = 0.9, and an rms scatter in α of 0.2 [START_REF] Griffin | The Herschel-SPIRE instrument and its in-flight performance[END_REF] . We took a hardness ratio corresponding to α = 0.5 (2σ from the mean) as the threshold beyond which we considered sources as absorbed. For absorbed sources, we used PIMMS to estimate the column density, assuming an intrinsic spectral index α = 0.9. PIMMS was then used to calculate the expected degree of attenuation to the observed 2-8 keV X-rays from the column density estimate. The X-ray luminosities were then adjusted accordingly. Where only a lower limit to the hardness ratio is available (due to the faintness of the object in the 0.5-2 keV X-ray band), the luminosity correction corresponding to the lower limit is applied. The fraction of absorbed AGN is known to depend on luminosity [START_REF] Ceca | The cosmological properties of AGN in the XMM-Newton Hard Bright Survey[END_REF] such that absorption is less common in the most luminous AGN. In light of this interdependency between luminosity and absorption, we have also examined the 250 µm detection rate of AGN with 1 < z < 3 and 10 43 erg s -1 < L X < 10 45 erg s -1 as a function of X-ray absorption rather than X-ray luminosity. Dividing the sample into absorbed and unabsorbed sources using the hardness-ratio criterion given above, we have 9 absorbed and 2 unabsorbed AGN detected at 250 µm, and 27 each of absorbed and unabsorbed sources which are not detected at 250 µm. Applying a Fisher's exact test, we find that the difference in 250 µm detection rate between absorbed and unabsorbed sources is significant at a little less than 95%. We can reduce the influence of the intrinsic anti-correlation between absorption and luminosity within the sample under test by excluding the sources with X-ray luminosities above 10 44 erg s -1 . We then have 9 absorbed and 2 unabsorbed sources detected at 250 µm, with 19 absorbed and 14 unabsorbed sources undetected at 250 µm. Applying the Fisher's test, the statistical significance at which the 250 µm detection rate depends on X-ray absorption drops to 90%. Hence, whereas we can demonstrate that the 250 µm detection rate depends on X-ray luminosity, we can not state with such confidence that the 250 µm detection rate depends on X-ray absorption.

Verification in the Extended Chandra Deep Field South

We have verified that the main result of this work, the paucity of 250 µm-detections amongst X-ray-luminous AGN, is reproduced in the Extended Chandra Deep Field South (ECDF-S), which is a larger field than the CDF-N, containing a larger sample of AGN with L X > 10 44 erg s -1 . We used the published X-ray catalogue [START_REF] Lehmer | The Extended Chandra Deep Field-South Survey: Chandra Point-Source Catalogs[END_REF] and spectroscopic redshifts [START_REF] Silverman | The Extended Chandra Deep Field-South Survey: Optical Spectroscopy of Faint X-ray Sources with the VLT and Keck[END_REF] . The SPIRE observations of the ECDF-S are shallower than those of CDF-N, but nonetheless are limited primarily by confusion. As with CDF-N, we have adopted a detection threshold of The images are arranged in increasing wavelength from left to right and decreasing X-ray luminosity from bottom to top. (b) Star formation rates obtained from the stacked fluxes and median redshifts (grey) compared to the star formation rates obtained from stacking the luminosities as in Fig. 2 (black). The blue, green and red points show the star formation rates obtained in a similar fashion to the black points, but in these cases using only one of the SPIRE bands to determine the luminosities prior to averaging. Blue, green and red correspond to 250, 350 and 500µm bands. The points have been offset slightly in the horizontal direction to aid clarity.

Figure 1 :

 1 Figure1: Redshifts and 2-8 keV X-ray luminosities of AGN in the CDF-N. The luminosities have been Kcorrected assuming a spectrum S ν ∝ ν -0.9 and are not corrected for intrinsic absorption. The blue dashed rectangles delimit the luminosity decades above and below 10 44 erg s -1 in the 1 < z < 3 redshift range. Error bars correspond to 68% confidence.

Fig. S2 shows

  Fig.S2shows how Fig.1would be affected by the absorption corrections. More than half of the 250 µmdetected AGN increase in luminosity when the absorption corrections are applied, but on the whole these corrections are quite modest, with the majority amounting to less than 40% change in luminosity. None of the 250 µm-detected AGN move into the 1 < z < 3, 10 44 erg s -1 < L X < 10 45 erg s -1 region. Although several of the 250 µm-detected AGN move close to the 10 44 erg s -1 boundary, and the absorption corrections for two of the 250 µm-detected AGN are lower limits such that they may have intrinsic luminosities > 10 44 erg s -1 , the essential characteristic of Fig.1, the low 250 µm detection rate above 10 44 erg s -1 , persists.

  18 mJy at 250 µm, and cross-correlated the X-ray and submillimetre catalogues with a matching radius of 6 arcseconds. The 250 µm detections in the ECDF-S (z, L) plane are shown in FigureS3, analogous to Fig.1of the main article. Of 49 sources with 1 < z < 3 and L X > 10 44 erg s -1 , only 1 is detected at 250 µm, confirming the low detection rate obtained in the CDF-N.

Figure S4 :

 S4 FigureS4: Infrared spectral energy distributions of 250 µm-detected AGN with 1 < z < 3 and 10 43 erg s -1 < L X < 10 44 erg s -1 , based on SPIRE and Spitzer data. We included 850 µm and 1.3 mm photometry from ground-based telescopes where available41;42 . The black datapoints are the measured photometry, and the red lines show the star-forming model templates[START_REF] Siebenmorgen | Dust in starburst nuclei and ULIRGs. SED models for observers[END_REF] which best fit the data at 24 µm and longward. The black lines show AGN templates[START_REF]Atlas of quasar energy distributions[END_REF] , normalised in the mid-IR to estimate the AGN contribution[START_REF] Seymour | HerMES: SPIRE emission from radio-selected active galactic nuclei[END_REF] . Each panel is labelled with the identity of the object in the X-ray catalogue[START_REF] Trouille | The OPTX Project. I. The Flux and Redshift Catalogs for the CLANS, CLASXS, and CDF-N Fields[END_REF] .

Fig. S5 ,

 S5 Fig. S5, and are consistent with each other as well as those shown in Fig. 2.

Table 1 :

 1 250 µm detection statistics in the various parts of parameter space. The first row corresponds to the entire sample of secure AGN in the CDF-N, while the second and third rows correspond to the regions enclosed within the blue dashed lines in Fig.1. Confidence intervals on the fraction of AGN associated with 250 µm sources are given at 68%, with 95% intervals enclosed in brackets. It should be noted that there is one case of two AGN being associated with the same 250 µm source. The two AGN have very similar spectroscopic redshifts, and both have X-ray luminosities between 10 43 and 10 44 erg s -1 . Although the two AGN cannot be resolved at 250 µm, source extraction using X-ray and 24 µm positions as priors[START_REF] Roseboom | The Herschel Multi-Tiered Extragalactic Survey: source extraction and crossidentifications in confusion-dominated SPIRE images[END_REF] indicates that both AGN are 5σ sources at 250 µm.signal/noise ratio >3 when the effects of confusion are included 17 . X-ray sources were matched to 250 µm sources within 6 arcseconds, corresponding to approximately 95% confidence in the 250 µm positions. The detection statistics are given in Table

	Region of (z,L) parameter space	Number Number of AGN Expected number of Fraction of AGN
		of AGN	associated with spurious associations	associated with
			250 µm sources		250 µm sources
	all z, 1 < z < 3, 10 43 < L X < 10 44 erg s -1 10 42 < L X < 10 45 erg s -1	176 44	24 11	2.1 0.5	14±3 ( +6 -5 ) % 25 +8 -7 ( +15 -12 ) %
	1 < z < 3, 10 44 < L X < 10 45 erg s -1	21	0	0.2	< 5 (< 13) %

Table 2 :

 2 Properties of SPIRE-detected AGN. Data are given for AGN with 1 < z < 3 and 10

	109 135 158 190 331 366 368 384 455 500	2.203 2.580 2.466 1.013 2.015 1.253 1.970 1.996 1.021 1.168 1.990	43.59 +0.08 -0.11 43.58 +0.07 -0.09 43.69 +0.07 -0.09 43.05 +0.04 -0.05 43.81 +0.04 -0.04 43.48 +0.03 -0.04 43.11 +0.08 -0.12 43.42 +0.07 -0.09 43.49 +0.03 -0.03 43.96 +0.02 -0.02 43.89 +0.10 -0.14	23.6 +0.1 -0.2 23.4 +0.1 - * > 24.0 23.1 +0.1 -0.1 23.6 +0.1 -0.1 -23.4 +0.1 -0.2 > 23.8 23.4 +0.1 -0.1 -23.2 +0.2 - *	0.15 +0.05 -0.07 0.09 +0.03 -0.09 > 0.30 0.15 +0.02 -0.02 0.16 +0.02 -0.02 -0.11 +0.03 -0.04 > 0.26 0.25 +0.02 -0.03 -0.07 +0.03 -0.07	12.70 ± 0.03 13.01 ± 0.05 12.81 ± 0.05 12.29 ± 0.09 12.88 ± 0.03 12.51 ± 0.07 12.84 ± 0.05 12.41 ± 0.06 11.55 ± 0.17 11.97 ± 0.04 12.62 ± 0.03	12 5 4 4 21 5 3 3 11 30 4	750 -850 1660 -1750 1060 -1110 320 -330 1030 -1300 530 -550 1140 -1170 420 -430 50 -60 110 -160 690 -710

  [START_REF]Atlas of quasar energy distributions[END_REF]

Acknowledgments SPIRE has been developed by a consortium of institutes led by Cardiff Univ. (UK) and including Univ. Lethbridge (Canada); NAOC (China); CEA, LAM (France); IFSI, Univ. Padua (Italy); IAC (Spain); Stockholm Observatory (Sweden); Imperial College London, RAL, UCL-MSSL, UKATC, Univ. Sussex 10 (UK); Caltech, JPL, NHSC, Univ. Colorado (USA). This development has been supported by national funding agencies: CSA (Canada); NAOC (China); CEA, CNES, CNRS (France); ASI (Italy); MCINN (Spain); SNSB (Sweden); STFC, UKSA (UK); and NASA (USA).

Supplementary Information is linked to the online version of the paper at www.nature.com/nature Author Contributions This paper represents the combined work of the HerMES collaboration, the SPIRE Instrument Team's extragalactic survey. M. Page planned the study, and wrote the draft version of the paper.

M. Symeonidis fitted models to the spectral energy distributions of the sources and J.D. Vieira performed the stacking analysis. All other coauthors of this paper contributed extensively and equally by their varied contributions to the SPIRE instrument, Herschel mission, analysis of SPIRE and HerMES data, planning of HerMES observations and scientific support of HerMES, and by commenting on this manuscript as part of an internal review process.

Author Information Reprints and permissions information is available at www.nature.com/reprints. The authors declare that they have no competing financial interests. Correspondence and requests for materials should be addressed to M.Page (mjp@mssl.ucl.ac.uk).

Supplementary Information

Spectroscopic redshifts of CDF-N X-ray sources

We have taken a recent compilation [START_REF] Barger | A Highly Complete Spectroscopic Survey of the GOODS-N Field[END_REF] as our main source of spectroscopic redshifts. These redshifts are given to two decimal places, but the majority of them originate from earlier studies, so we have in most cases obtained more precise redshifts by tracing them back to the original publications. Where there is a significant discrepancy (> 2%) between a spectroscopic redshift in our main source and that given in earlier literature, we have adopted the redshift from our main source [START_REF] Barger | A Highly Complete Spectroscopic Survey of the GOODS-N Field[END_REF] . We have utilised a number of spectroscopic redshifts which are not listed there. The redshifts and literature references for these cases are given in Table S1.

Effect of redshift incompleteness

While we have restricted our study to objects within the CDF-N which have optical spectroscopic identifications and redshifts to ensure fidelity in redshifts and luminosities, we must examine whether the limited spectroscopic completeness (64%) has any effect on our results. A large fraction of the X-ray sources within the CDF-N have photometric redshifts (i.e. redshifts estimated by comparing photometry in a number of bands to that expected from galaxy template spectra as a function of redshift) in the same compilation that we use for our spectroscopic redshifts [START_REF] Barger | A Highly Complete Spectroscopic Survey of the GOODS-N Field[END_REF] . Although these photometric redshifts are of low accuracy and reliability compared to the spectroscopic redshifts, they raise the redshift completeness of the 2-8 keV X-ray sources to 87%, and so offer a good test as to whether the spectroscopic completeness has a significant influence on our results. Table S2 gives statistics equivalent to those shown in Table 1 of the main paper, but with the photometric redshifts included. As found for the spectroscopic-only sample, no AGN more luminous than 10 44 erg s -1 are associated with 250 µm sources.

We can be further reassured that our main result is not affected by the spectroscopic incompleteness by looking specifically at the identification statistics of those 2-8 keV X-ray sources which are associated with 250 µm sources. A total of 33 2-8 keV X-ray sources are associated with 250 µm sources, of which 27 are spectroscopically identified. The spectroscopically-identified fraction of these sources is thus 82%, higher than the 2-8 keV source population as a whole. For the 6 sources without spectroscopic redshifts, Fig. S1 shows the tracks of X-ray luminosity as a function of redshift. We can see that with respect to the 1 < z < 3, 10 43 erg s -1 < L X < 10 45 erg s -1 region which is well sampled by the X-ray and submillimetre surveys, sources 31, 37 and 431 from the X-ray catalogue [START_REF] Trouille | The OPTX Project. I. The Flux and Redshift Catalogs for the CLANS, CLASXS, and CDF-N Fields[END_REF] have L X < 10 44 erg s -1 for all 1 < z < 3. The tracks for sources 107, 221 and 406 all cross the upper rectangle corresponding to 1 < z < 3, 10 44 erg s -1 < L X < 10 45 erg s -1 , and so should be examined in more detail. The track for source 221 only enters this rectangle for a tiny redshift interval (2.94 < z < 3.0). Although our main redshift compilation [START_REF] Barger | A Highly Complete Spectroscopic Survey of the GOODS-N Field[END_REF] does not give a photometric redshift for source Infrared spectral energy distributions of 250 µm-detected AGN

The spectral energy distributions and best fit models for the AGN associated with 250 µm sources with 1 < z < 3 and 10 43 erg s -1 < L X < 10 44 erg s -1 , are shown in Fig. S4 Obtaining the average star formation rates through stacking

The average star formation rates shown in Fig. 2 were obtained by averaging the far-infrared luminosities of the AGN to account for the different redshifts of the different AGN. It is more common in stacking analyses to average the fluxes of the objects, and then to use the average flux and the median redshift of the sample to infer an average luminosity. Here we examine how the results would change if we used this latter method rather than stacking directly in luminosity. Fig. S5 shows the star formation rates obtained by the two different methods. The two methods give consistent results, and in both cases the average star formation rates are significantly higher for AGN with log L X < 10 44 erg s -1 than for AGN with log L X > 10 44 erg s -1 . We have also examined how the star formation rates would differ if only one of the SPIRE bands were used, instead of all three, to determine the luminosities. For each of the SPIRE bands, we derived a luminosity for each source by normalising a 30 K grey-body in the rest frame to the SPIRE flux and multiplying by 4πD 2 L , before averaging the luminosities in each bin. The star formation rates so derived from the three SPIRE bands are shown in


Anais Barateau
email: anais.barateau@gmail.com

Renaud De Crevoisier

Axel Largent

Eugenia Mylona

Nicolas Perichon

Joël Castelli

Enrique Chajon

Oscar Acosta

Antoine Simon

J-C Nunes

Anaïs Barateau

Jean- Claude Nunes

Caroline Lafond

Eugène Centre

Marquis

Comparison of CBCT-based dose calculation methods in head and neck cancer radiotherapy: from Hounsfield unit to density calibration curve to deep learning

Keywords: Deep learning, CBCT dose calculation, head and neck cancer

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

Introduction

Recent developments in modulated treatment plans and image-guided radiotherapy (IGRT) devices have allowed more precise and targeted head and neck (H&N) treatments, with improved sparing of organs at risk while covering target volumes. However, during H&N radiotherapy, patients are typically subject to anatomical variations such as tumor shrinkage or weight loss [START_REF] Brouwer | Identifying patients who may benefit from adaptive radiotherapy: Does the literature on anatomic and dosimetric changes in head and neck organs at risk during radiotherapy provide information to help?[END_REF][START_REF] Castadot | Adaptive Radiotherapy of Head and Neck Cancer[END_REF] .

These variations can induce discrepancies between planned and delivered doses. Acquisitions of daily or weekly 3D kV cone-beam computed tomography (CBCT) images can be used to quantify dose discrepancies related to anatomical variations. CBCT images can be used to calculate the delivered "dose of the day" [START_REF] Veiga | Toward adaptive radiotherapy for head and neck patients: Feasibility study on using CT-to-CBCT deformable registration for "dose of the day" calculations[END_REF] to viewable and critical structures such as the spinal cord and parotid glands (PGs). However, dose calculation from CBCT is challenging not only owing to the "poor" image quality, but also to the limited field of view (FOV) and the inconsistency of the Hounsfield units (HUs). Moreover, CBCT dose calculation often lacks ground truth (dose computed from the reference CT acquired at the same time as CBCT) for quantifying the uncertainties.

Several methods to perform CBCT dose calculation have been proposed: i) calibration curve between the HU and densities (HU-D curve) [START_REF] Richter | Investigation of the usability of conebeam CT data sets for dose calculation[END_REF][START_REF] Fotina | Feasibility of CBCTbased dose calculation: Comparative analysis of HU adjustment techniques[END_REF][START_REF] Yoo | Dosimetric feasibility of cone-beam CT-based treatment planning compared to CTbased treatment planning[END_REF] , ii) density assignment method (DAM) [START_REF] Fotina | Feasibility of CBCTbased dose calculation: Comparative analysis of HU adjustment techniques[END_REF][START_REF] Dunlop | Comparison of CT number calibration techniques for CBCT-based dose calculation[END_REF][START_REF] Van Zijtveld | Correction of conebeam CT values using a planning CT for derivation of the "dose of the day[END_REF][START_REF] Barateau | A density assignment method for dose monitoring in headand-neck radiotherapy[END_REF][START_REF] Giacometti | An evaluation of techniques for dose calculation on cone beam computed tomography[END_REF] , iii) deformable image registration (DIR) between CT and CBCT [START_REF] Veiga | Toward adaptive radiotherapy for head and neck patients: Feasibility study on using CT-to-CBCT deformable registration for "dose of the day" calculations[END_REF][START_REF] Giacometti | An evaluation of techniques for dose calculation on cone beam computed tomography[END_REF][START_REF] Marchant | Accuracy of radiotherapy dose calculations based on cone-beam CT: comparison of deformable registration and image correction based methods[END_REF][START_REF] Kurz | Investigating deformable image registration and scatter correction for CBCT-based dose calculation in adaptive IMPT: CBCT correction to enable IMPT dose calculation[END_REF] and iv) machine learning to generate a pseudo-CT (pCT) [START_REF] Li | A preliminary study of using a deep convolution neural network to generate synthesized CT images based on CBCT for adaptive radiotherapy of nasopharyngeal carcinoma[END_REF][START_REF] Chen | Synthetic CT generation from CBCT images via deep learning[END_REF][START_REF] Liang | Generating Synthesized Computed Tomography (CT) from Cone-Beam Computed Tomography (CBCT) using CycleGAN for Adaptive Radiation Therapy[END_REF] . i) The HU-D curve established from a CBCT image can be used to convert CBCT HUs to densities for dose calculation. This curve can be defined with either an "adapted" phantom [START_REF] Yoo | Dosimetric feasibility of cone-beam CT-based treatment planning compared to CTbased treatment planning[END_REF] (according to anatomical localization) or patient CBCT images [START_REF] Richter | Investigation of the usability of conebeam CT data sets for dose calculation[END_REF][START_REF] Fotina | Feasibility of CBCTbased dose calculation: Comparative analysis of HU adjustment techniques[END_REF] . Although these methods are straightforward, they are sensitive to CBCT artefacts [START_REF] Schulze | Artefacts in CBCT: a review[END_REF] and patient scattering. ii)

The density assignment method (also known as the bulk density method) involves segmenting an image into two to six tissue classes (e.g., soft tissues, air and bones) before assigning density to each class. Nevertheless, this method is dependent of structure segmentation and provide an image with homogeneous tissues. iii) By deforming CT to CBCT, a "deformed" CT is generated and can be used for dose calculation. CT-CBCT DIR can be difficult owing to intrinsic CBCT limitations, such as noise, low contrast, and reduced FOV. Moreover, the registration step is complex owing to large anatomical variations, such as tumor shrinkage or weight loss [START_REF] Rigaud | Deformable image registration for radiation therapy: principle, methods, applications and evaluation[END_REF] . iv) Machine learning methods are based on patches or deep learning (DL), to generate a pCT (i.e., synthetic images) from CBCT. Machine learning methods require a large training cohort and most of these methods

Accepted Article

require co-registered data for the training step. The main advantages and drawbacks of the HU-D curve method, DAM, DIR method, and DLM are summarized in the two first columns of Table 1.

Deep learning methods (DLMs) are models composed of multiple processing layers that learn multiscale representations of data through multiple levels of abstraction [START_REF] Lecun | Deep learning[END_REF] . These methods have recently been introduced in radiotherapy for various applications, such as image segmentation, image reconstruction, image registration, treatment planning, and radiomics [START_REF] Meyer | Survey on deep learning for radiotherapy[END_REF][START_REF] Higaki | Improvement of image quality at CT and MRI using deep learning[END_REF][START_REF] Alkadi | A Deep Learning-Based Approach for the Detection and Localization of Prostate Cancer in T2 Magnetic Resonance Images[END_REF][START_REF] Laukamp | Fully automated detection and segmentation of meningiomas using deep learning on routine multiparametric MRI[END_REF][START_REF] Liang | Deep-learning-based detection and segmentation of organs at risk in nasopharyngeal carcinoma computed tomographic images for radiotherapy planning[END_REF][START_REF] Nyflot | Deep learning for patientspecific quality assurance: Identifying errors in radiotherapy delivery by radiomic analysis of gamma images with convolutional neural networks[END_REF][START_REF] Sahiner | Deep learning in medical imaging and radiation therapy[END_REF] . DLMs have been primarily proposed for pCT generation from magnetic resonance imaging (MRI) [START_REF] Gong | Attenuation correction for brain PET imaging using deep neural network based on Dixon and ZTE MR images[END_REF][START_REF] Ladefoged | Deep Learning Based Attenuation Correction of PET/MRI in Pediatric Brain Tumor Patients: Evaluation in a Clinical Setting[END_REF][START_REF] Leynes | Zero-Echo-Time and Dixon Deep Pseudo-CT (ZeDD CT): Direct Generation of Pseudo-CT Images for Pelvic PET/MRI Attenuation Correction Using Deep Convolutional Neural Networks with Multiparametric MRI[END_REF][START_REF] Torrado-Carvajal | Dixon-VIBE Deep Learning (DIVIDE) Pseudo-CT Synthesis for Pelvis PET/MR Attenuation Correction[END_REF][START_REF] Kläser | Deep Boosted Regression for MR to CT Synthesis[END_REF][START_REF] Nyholm | MR and CT data with multiobserver delineations of organs in the pelvic area-Part of the Gold Atlas project[END_REF] . They are particularly appealing owing to their fast computation time. One of the first DLMs for pCT generation was based on the U-Net architecture. More recently, DLMs using generative adversarial networks (GANs) have been proposed, offering the theoretical advantage of providing more realistic pCTs by obtaining an adversarial feedback from a discriminator network [START_REF] Yi | Generative Adversarial Network in Medical Imaging: A Review[END_REF] . Some studies have recently proposed DLMs for pCT generation from CBCT, mainly for scatter correction [START_REF] Maier | Real-time scatter estimation for medical CT using the deep scatter estimation: Method and robustness analysis with respect to different anatomies, dose levels, tube voltages, and data truncation[END_REF][START_REF] Nomura | Projection-domain scatter correction for cone beam computed tomography using a residual convolutional neural network[END_REF] . Other studies proposed the DLM for pCT generation from CBCT in prostate [START_REF] Landry | Comparing Unet training with three different datasets to correct CBCT images for prostate radiotherapy dose calculations[END_REF][START_REF] Hansen | ScatterNet: A convolutional neural network for cone-beam CT intensity correction[END_REF][START_REF] Kurz | CBCT correction using a cycle-consistent generative adversarial network and unpaired training to enable photon and proton dose calculation[END_REF] , pancreas [START_REF] Liu | CBCT-based Synthetic CT Generation using Deep-attention CycleGAN for Pancreatic Adaptive Radiotherapy[END_REF] , and H&N [START_REF] Li | A preliminary study of using a deep convolution neural network to generate synthesized CT images based on CBCT for adaptive radiotherapy of nasopharyngeal carcinoma[END_REF][START_REF] Chen | Synthetic CT generation from CBCT images via deep learning[END_REF][START_REF] Liang | Generating Synthesized Computed Tomography (CT) from Cone-Beam Computed Tomography (CBCT) using CycleGAN for Adaptive Radiation Therapy[END_REF] , for dose calculation. H&N studies have been performed using either U-Net or cycleGAN architectures to generate pCT from CBCT.

To the best of our knowledge, no DLM to generate pCT from H&N CBCT has been compared with other methods for dose calculation from CBCT. In H&N radiotherapy, previous studies showed dose differences of less than 3% for the DAM [START_REF] Fotina | Feasibility of CBCTbased dose calculation: Comparative analysis of HU adjustment techniques[END_REF][START_REF] Dunlop | Comparison of CT number calibration techniques for CBCT-based dose calculation[END_REF][START_REF] Van Zijtveld | Correction of conebeam CT values using a planning CT for derivation of the "dose of the day[END_REF][START_REF] Barateau | A density assignment method for dose monitoring in headand-neck radiotherapy[END_REF][START_REF] Giacometti | An evaluation of techniques for dose calculation on cone beam computed tomography[END_REF] , less than 2% for the DIR method [START_REF] Veiga | Toward adaptive radiotherapy for head and neck patients: Feasibility study on using CT-to-CBCT deformable registration for "dose of the day" calculations[END_REF][START_REF] Giacometti | An evaluation of techniques for dose calculation on cone beam computed tomography[END_REF][START_REF] Marchant | Accuracy of radiotherapy dose calculations based on cone-beam CT: comparison of deformable registration and image correction based methods[END_REF] , and less than 1% for the DLM.

This study aims to evaluate the accuracy of the DL method for CBCT H&N dose calculation and to compare this method with the HU-D curve, density assignment, and DIR methods.

Materials and methods

Patient data

Forty four patients with locally advanced oropharyngeal carcinomas were retrospectively selected.

Image acquisition

All patients had a planning CT scan, and 14 patients had weekly CT scans. All CT images were acquired on a BigBore (Philips) scanner, with a 2 mm slice thickness from the vertex to the carina.

Patient positioning was assured using a personalized thermoplastic head and shoulder mask with

Accepted Article

five fixation points. The treatment isocenter was set at a reproducible bony place, between the C3-C4 vertebrae.

Weekly CBCTs were acquired with an XVI (Elekta) on a VERSAHD linac (Elekta). The settings for the H&N CBCT acquisition were as follows: 120 kV, 0.4 mAs per frame, 660 frames, M20 FOV, and 2 mm slice thickness. For the14 patients who had weekly CT scans, the weekly CBCTs were acquired for a time close to the CTs.

Delineation and dose calculation

Structure delineation of the planning CTs was performed by the same radiation oncologist according to the European recommendations [START_REF] Brouwer | CT-based delineation of organs at risk in the head and neck region: DAHANCA, EORTC, GORTEC, HKNPCSG, NCIC CTG, NCRI, NRG Oncology and TROG consensus guidelines[END_REF] . Positron emission tomography and MRI coregistration were used for tumor delineation. The gross tumor volume (GTV) corresponded to the primary tumor and involved lymph nodes. Three target volumes were generated for each patient.

A clinical target volume receiving 70 Gy (CTV 70 , 35 fractions) was equal to the GTV plus a 5 mm 3D margin, adjusted to exclude any air cavities and bone mass that indicate no evidence of tumor invasion. CTV 63 (receiving 63 Gy in 35 fractions) corresponded to the high-risk area of the microscopic spread, whereas CTV 56 (receiving 56 Gy in 35 fractions) corresponded to the low-risk subclinical area. The GTV, CTV 63 , CTV 56 , and all the organs at risk, in particular the PGs and spinal cord, were manually delineated on each CT slice. The planning target volumes (PTV 70Gy , PTV 63Gy , and PTV 56Gy ) were generated by adding a 5 mm 3D margin around the CTVs. For all patients, the total prescribed dose was 70 Gy (35 fractions) delivered by a simultaneous integrated boost (70-63-56 Gy) using VMAT technique. The dose calculation was performed with Pinnacle v.9.10 (Philips) treatment planning system (TPS) using the collapsed cone convolution algorithm and a dose grid resolution of 3 mm. GORTEC recommendations were used for organs at risk limitation doses [START_REF] Toledano | Intensity-modulated radiotherapy in head and neck cancer: Results of the prospective study GORTEC 2004-03[END_REF] . The dose constraints for the contralateral PG were a mean dose lower than 30 Gy and a median dose lower than 26 Gy. For the ipsilateral PG, the mean dose was as low as possible. For the spinal cord, the maximum dose was lower than 45 Gy. The minimum PTV volume coverage by the 95% isodose was 95%.

On the pCT images generated by the evaluated methods, structures (target volumes, PGs and spinal cord) were propagated from the planning CTs, and the dose distributions were calculated after copying the beam parameters.

Accepted Article

Intra-patient CBCT to CT registration

To minimize the effects of anatomic variations between CT and CBCT, each CBCT was registered to its corresponding CT initial using a rigid registration, followed by a non-rigid registration (NiftyReg). This CT initial was considered as the reference (CT ref ).

Deep learning method

The GAN DLM architecture was composed of two networks: a generator (G) and a discriminator (D), which were trained in competition with each other (Fig. 1).

Generator network

The generator network aims to provide a pCT from each patient's CBCT. In this study, the generator network used a 2D architecture similar to the U-Net DLM proposed by Han [START_REF] Han | MR-based synthetic CT generation using a deep convolutional neural network method[END_REF] . A singlescale perceptual loss function was used to train this network. This loss function mimics the human visual system to compare the CT and pCT images using similar features as opposed to only the intensities [START_REF] Nie | Medical Image Synthesis with Context-Aware Generative Adversarial Networks[END_REF][START_REF] Yang | Low-Dose CT Image Denoising Using a Generative Adversarial Network With Wasserstein Distance and Perceptual Loss[END_REF] . The features inside the CT and pCT images were computed using the Visual Geometry Group (VGG) 16 network, which is used for perceptual loss computation and appears relevant for various tasks (e. g., image deblurring, super-resolution, and computer vision) [START_REF] Yang | Low-Dose CT Image Denoising Using a Generative Adversarial Network With Wasserstein Distance and Perceptual Loss[END_REF][START_REF] Simonyan | Very Deep Convolutional Networks for Large-Scale Image Recognition[END_REF] . The perceptual loss function of the generator (L G ) is defined as follows:

𝐿 𝐺 (𝐼, 𝐶) = ‖𝑉𝐺𝐺(𝐶) -𝑉𝐺𝐺(𝐺(𝐼))‖ 2 2
where is the output of the 7 th VGG16 convolutional layer, the CBCT, and the 𝑉𝐺𝐺 𝐼 𝐶 corresponding CT.

Discriminator network

The discriminator network aims to classify the pCT image as a real or fake CT. Hence, the output of this network is a probability value ranging between 0 and 1 depending on whether the pCT appears to be fake or real. The architecture was composed of six convolutional layers and one fully connected layer. Each convolutional layer was followed by batch normalization and Leaky-ReLu activation functions. The number of filters for these six layers were 8, 16, 32, 64, 64, and 64.

Accepted Article

The filter size was 3 × 3 (stride = 2) for the first four layers and 1 × 1 (stride = 1) for the remaining layers. The fully connected layer had one filter followed by a sigmoid activation function.

The loss function of the discriminator (L D ) was a binary cross entropy expressed as 𝐿 𝐷 (𝐺(𝐼), 𝐶)

where is the pCT computed by the

= -∑ 𝑛 𝑖 = 1 𝐶 𝑖 𝑙𝑜𝑔(𝐺(𝐼) 𝑖 ) + (1 -𝐶 𝑖 )log (1 -𝐺(𝐼) 𝑖 ) ,

𝐺(𝐼)

generator from the target CBCT , and is the number of voxels inside and images.

𝐼 𝑛 𝐶 𝐼

The generator and discriminator losses were combined to form the following adversarial loss:

, where is the discriminator loss,

𝐿 𝑎𝑑𝑣𝑒𝑟𝑠𝑎𝑟𝑖𝑎𝑙 (𝐼, 𝐶) = 𝜆 1 𝐿 𝐷 (𝐼, 𝐶) + 𝜆 2 𝐿 𝐺 (𝐼, 𝐶) 𝐿 𝐷 (𝐼, 𝐶) 𝐿 𝐺 (𝐼, 𝐶)
is the generator loss, and and are the weights of the discriminator and generator losses, 𝜆 1 𝜆 2 respectively. First, the discriminator was trained using the discriminator loss, followed by generator training using the fully adversarial loss. These training steps were performed iteratively until 300 epochs when the discriminator could not accurately determine if the pCTs provided by the generator were real or false CTs.

Training data

The DLM was trained using 30 anatomically paired (CT-CBCT) data: axial 2D slices of the training CT and CBCT images (7600 slices). Data augmentation was performed to artificially increase the size of the training cohort. It was conducted by randomly applying affine registrations (by NiftyReg) on the slices (translated from -5% to 5% per axis, rotated from -10° to +10°, sheared from -10° to 10°). A mini-batch size of four slices and 300 epochs was considered. The network parameters were optimized using the Adam algorithm [START_REF] Kingma | A Method for Stochastic Optimization[END_REF] . The parameters of this algorithm were as follows: , , and For the GAN, the weights of the 𝛼 = 1 × 10 -4 𝛽 1 = 0.9 𝛽 2 = 0.9.

discriminator and generator loss functions were: and , respectively. The convergence

𝜆 1 = 5 𝜆 2 = 1
curves of the GAN generator and discriminator are presented in Fig. S-1.

The DLM was implemented in Python using Keras [START_REF]Contribute to Keras-Team/Keras Development by Creating an Account on GitHub[END_REF] . The training computation time for the networks was approximately 72 h using an Nvidia GTX 1070 TI 8 GB GPU.

Other CBCT-based dose calculation methods

Accepted Article

The three other assessed methods are described following.

i. HU to density curve (HU-D curve) from phantom CBCT image

This method involves establishing an HU-D curve of the CIRS 062 phantom with clinical H&N CBCT protocol acquisition. Only the central part of the phantom was used to mimic the dimensions of the patient head. This phantom contained eight heterogeneous inserts and is considered appropriate for CBCT dose calculation owing to its size and chemical composition [START_REF] Guan | Dose calculation accuracy using cone-beam CT (CBCT) for pelvic adaptive radiotherapy[END_REF] .

An HU-D H&N curve was then implemented into the TPS. iii.

CT-CBCT deformable image registration (DIR)

A deformed image was created with automatic DIR between the CBCT (fixed image) and planning CT (moving image) images using Admire software (research v.3.3.1, Elekta). A study was performed to assess the DIR accuracy of Admire (Appendix 1).

Endpoints

Imaging and dosimetric endpoints were considered for the 14 patients (total of 64 CT-CBCTs).

Imaging endpoints

A 𝐻𝑈 𝐶𝑇𝑟𝑒𝑓 (𝑖) -𝐻𝑈 𝑝𝐶𝑇 (𝑖) of voxels.

Dosimetric endpoints

Owing to the limited size of the CBCT FOV, the water equivalent density was assigned inside the CT body contour when no CBCT or pCT information was available to perform an evaluation of the method. The accuracy of the methods was first evaluated by computing the dose and systematic dose discrepancies. The DVH calculated from the CT ref was subtracted from the DVH calculated from the pCT. The dose discrepancy was defined by the mean absolute dose and the systematic dose discrepancy by the mean DVH differences. The results were reported for the GORTEC reference DVH points (V 95% for the PTVs and D 2% for the spinal cord and mandible), and D mean to the PGs was considered. A spatial dose evaluation was finally conducted by performing 3D gamma analyses (local, 2%/2 mm, low dose threshold: 30%) using the dose distributions from the CT ref and pCTs.

Statistical analysis

Wilcoxon signed-rank tests were performed to compare the endpoints. For the DVH comparisons based on the pCT generation methods, a nonparametric permutation test was performed [START_REF] Ross | Chapter 14-Nonparametric Hypotheses Tests[END_REF] to control the presence of false positives in multiple statistical tests (10 cGy DVH bin-wise). In this case, 1000 permutations were performed, where for each permutation i, randomly selected DVHs were swapped (CTref <->pCT) and the average difference was computed for each dose-bin. For each permuted sample and the original sample, the average difference was normalized to the standard deviation computed over all the 1000 permutations, and the maximum observed difference was selected as the test statistic (TS). A distribution of TS across all the permuted samples (TS i,max ) was obtained and compared with that one from the observed sample (TS max ). The adjusted p-value was then computed as the probability of having a TS max greater than the TS i,max at the significance level of 5% (p ≤ 0.05). The corresponding percentile over the distribution of all the TS i,max provides a threshold value that determines the dose DVH bins where a statistically significant dose difference arises. Unlike bin-wise tests, a permutation test provides a single

Accepted Article

number that summarizes the discrepancy in the DVH between the two groups, rather than the discrepancy of a particular bin, which therefore accounts for multiple comparisons. The mathematical formulation of the permutation test has been provided by Chen et al. [START_REF] Chen | Multiple comparisons permutation test for image based data mining in radiotherapy[END_REF] . The test allowed us to report a robust bin-wise comparison based on the DVH value of each method as well as to compare the lowest MAE among all the methods with the MAE of each method and the ME of each method with a null distribution.

Results

Imaging endpoints and calculation time

Examples of CBCT, CT ref , and pCT generated by each method are illustrated in Fig. 2. The mean calculation time to generate one pCT was 30 s for the GAN DLM.

Dosimetric endpoints

Fig. 3 shows the mean DVHs for the CT ref and each method, by volume of interest. The DVHs of all methods differed significantly from that of CT ref .

Table 3 shows the D mean differences of DVH points for the PTVs, spinal cord, and mandible and the D mean to the PGs. The density assignment as well as DIR and DL methods provided dose discrepancies lower than 0.6% for the PTVs, PGs, and spinal cord. The lowest absolute mean D max difference of the spinal cord was obtained using the DIR method (17 ± 20 cGy). The lowest absolute mean D mean differences of the PGs were obtained using the DLM. The mean differences differed significantly from the null distribution for the PTV, spinal cord, and mandible DVH points in the four methods. The mean D mean differences of the PGs did not significantly differ from the null distribution in the four methods. Fig. 4 shows the dose discrepancies of each method along the DVHs, by volume of interest. The method with the

Accepted Article

highest discrepancies was the HU-D curve method. The density assignment and DIR methods were less significantly different from the DLM compared with the HU-D curve method. Fig. S-2 shows the systematic dose discrepancies (ME) of each method along the DVHs, by volume of interest. All methods showed significant differences with a null distribution.

Table 4 shows the mean gamma and gamma pass-rate values calculated from the CT ref and pCT dose distributions for each method. The density assignment, DIR, and DL methods provided mean gamma pass-rates that exceeded 97.9%. The gamma pass-rates and mean gamma results of the HU-D curve method, DAM, and DIR method differed significantly from those of the DLM. Fig.

S-3 shows examples of gamma maps in the coronal plane for one patient.

The last two columns of Table 1 are a qualitative analysis of the image and dose endpoints among the four compared methods.

Discussion

The aim of this study was to compare a DLM (GAN) with three other methods (HU-D curve method, DAM, and DIR method) to perform dose calculation from H&N CBCT. Compared with the CT ref , the DAM, DIR method, and DLM provided low dose discrepancies, thereby rendering them clinically acceptable for CBCT-based dose calculation (Fig. 3). Regarding dose accuracy as well as calculation time (< 30 s), the DIR method and DLM appeared to be the most attractive methods (Table 1).

The main disadvantages of CBCT for dose calculation is HU inconsistency. Three studies have proposed comparing several CBCT dose calculation methods (scatter correction, DAM, DIR) but not with the DLM [START_REF] Dunlop | Comparison of CT number calibration techniques for CBCT-based dose calculation[END_REF][START_REF] Giacometti | An evaluation of techniques for dose calculation on cone beam computed tomography[END_REF][START_REF] Marchant | Accuracy of radiotherapy dose calculations based on cone-beam CT: comparison of deformable registration and image correction based methods[END_REF] . We conducted a comparison study using the DLM and three other methods (well-known methods in the literature) for CBCT dose calculation. This comparison study was performed using imaging and dose endpoints. Table 1 summarizes the benefits, drawbacks and qualitative comparison of the four investigated methods for H&N CBCT dose

Accepted Article

calculation. The DLM provided the lowest MAE (Table 2) in the entire body (82.4 HU), soft tissue (69.2 HU), and bones (207.6 HU). However, the DLM did not show the lowest dose discrepancies. The relationship between image and dose uncertainties is particularly complex. In fact, the dose uncertainty in a voxel is not only related to the image uncertainty (HU) in the considered voxel, but also to several parameters, such as the image uncertainty of the surrounding voxels or dose gradients. For the bone tissues, the systematic error of the methods (HU-D curve method, DAM, and DLM) showed an underestimation of HU values (Table 2), whereas dose underestimation was recorded for the spinal cord (surrounded by bone) and mandible (Table 3).

The density assignment and DIR methods provided better or similar dose results to the DLM (Table 3, Fig. 4, Fig. S-2, and Table 4). The method with the worst dose discrepancy (Table 3, Table 4, Fig. 4), except for the mandible, was the method using the HU-D conversion curve. This simple method appeared to be insufficient owing to artefacts and patient scattering [START_REF] Schröder | Evaluating the impact of cone-beam computed tomography scatter mitigation strategies on radiotherapy dose calculation accuracy[END_REF][START_REF] Schröder | Technical Note: Long-term stability of Hounsfield unit calibration for cone-beam computed tomography[END_REF] . Using this method, the main dose discrepancies were obtained at the shoulder level (Fig. S-3), where scatter artefacts are more present. The PTV 56 Gy was the most affected structure because of this issue (Fig. 2).

We proposed a DLM to generate pCT from CBCT to perform dose calculation. Such methods have already been used for pCT generation from MRI in the context of an MRI-only workflow [START_REF] Wang | Medical Imaging Synthesis using Deep Learning and its Clinical Applications: A Review[END_REF] .

The two recent MRI H&N studies using conditional GAN [START_REF] Qi | Multi-sequence MR image-based synthetic CT generation using a generative adversarial network for head and neck MRI-only radiotherapy[END_REF] or U-Net 54 provided an MAE of approximately 70-75 HU in the entire body contour. However, the aim of pCT generation from CBCT images is different. In fact, pCT generated from CBCT are used to monitor delivered doses or to estimate the cumulative delivered dose during the treatment course in the context of doseguided adaptive radiotherapy. Studies using DL for pCT generation from CBCT are scarce for brain [START_REF] Harms | Paired cycle-GAN-based image correction for quantitative cone-beam computed tomography[END_REF] , H&N [START_REF] Li | A preliminary study of using a deep convolution neural network to generate synthesized CT images based on CBCT for adaptive radiotherapy of nasopharyngeal carcinoma[END_REF][START_REF] Chen | Synthetic CT generation from CBCT images via deep learning[END_REF][START_REF] Liang | Generating Synthesized Computed Tomography (CT) from Cone-Beam Computed Tomography (CBCT) using CycleGAN for Adaptive Radiation Therapy[END_REF] , pancreas [START_REF] Liu | CBCT-based Synthetic CT Generation using Deep-attention CycleGAN for Pancreatic Adaptive Radiotherapy[END_REF] or prostate cancer [START_REF] Landry | Comparing Unet training with three different datasets to correct CBCT images for prostate radiotherapy dose calculations[END_REF][START_REF] Hansen | ScatterNet: A convolutional neural network for cone-beam CT intensity correction[END_REF][START_REF] Kurz | CBCT correction using a cycle-consistent generative adversarial network and unpaired training to enable photon and proton dose calculation[END_REF][START_REF] Harms | Paired cycle-GAN-based image correction for quantitative cone-beam computed tomography[END_REF][START_REF] Kida | Visual enhancement of Cone-beam CT by use of CycleGAN[END_REF] . The studies involved an imaging analysis (pCT versus reference CT), but only half of them evaluated the dose accuracy. Among the three H&N studies using DL for pCT generation from CBCT [START_REF] Li | A preliminary study of using a deep convolution neural network to generate synthesized CT images based on CBCT for adaptive radiotherapy of nasopharyngeal carcinoma[END_REF][START_REF] Chen | Synthetic CT generation from CBCT images via deep learning[END_REF][START_REF] Liang | Generating Synthesized Computed Tomography (CT) from Cone-Beam Computed Tomography (CBCT) using CycleGAN for Adaptive Radiation Therapy[END_REF] , one involved training a U-Net neural network on 50 co-registered CBCT/CT images and performing a test based on data from 10 patients [START_REF] Li | A preliminary study of using a deep convolution neural network to generate synthesized CT images based on CBCT for adaptive radiotherapy of nasopharyngeal carcinoma[END_REF] . MAEs from 6 to 27 HU were obtained for the pCTs generated from Varian CBCTs. The dose differences were lower than 1%, and the mean gamma pass-rate (global, dose difference: 1%, DTA: 1 mm) was 95.5%. The second study involved training a cycle-consistent generative adversarial network (cycleGAN) on 81 CBCTs [START_REF] Liang | Generating Synthesized Computed Tomography (CT) from Cone-Beam Computed Tomography (CBCT) using CycleGAN for Adaptive Radiation Therapy[END_REF] . The evaluation performed on data from 20 patients provided a mean MAE of 29.9 HU. The mean gamma pass-rate (3D, dose difference: 1%,

Accepted Article

DTA: 1 mm) was 96.3%. The cycleGAN DLM involved two GANs: one to generate pCT from CBCT and a second GAN to generate pCBCT from pCT (the output of the first GAN). These two dual GANs learn simultaneously to improve theoretically the efficiency. The cycleGAN DLM does not require paired data (no challenging intra-patient registration). Therefore, the amount of training data can largely be increased. However, the cycleGAN requires more complex training compared with the GAN because a larger number of parameters require optimization. The third study involved U-Net training on 37 patients to generate pCT from CBCT [START_REF] Chen | Synthetic CT generation from CBCT images via deep learning[END_REF] , where only image analysis was performed. Furthermore, the MAE of 19 HUs for seven patients was obtained. We obtained a higher MAE with a larger number of patients and Elekta CBCT acquisitions. It is noteworthy that the comparison can only be indirect because the patient data differed in each study (CBCT imaging device, FOV size, tumor location, etc.).

Our study has some limitations. First, we trained the DL network with data from only 30 patients.

Second, before the learning process, non-rigid registration was used to obtain the same H&N anatomy between CBCT and CT ref , with the uncertainties depending of the deformable image registration algorithm. Third, the limited size of the CBCT FOV is still an issue. In this study, we added water equivalent density to override it and evaluate only the methods. Moreover, no scatter correction [START_REF] Rührnschopf | A general framework and review of scatter correction methods in xray cone-beam computerized tomography. Part 1: Scatter compensation approaches: Scatter compensation approaches[END_REF] was applied on CBCT images. Such image correction in a pre-processing step could decrease image and dose discrepancies. For image analysis, we did not consider metrics other than HU differences (MAE and ME) because the focus was on dose calculation. Image quality metrics such as Peak Signal-to-Noise Ratio (PSNR), Normalized Mutual Information (NMI), Structural SIMilarity (SSIM) or Visual Information Fidelity (VIF) should be implemented for further image analysis. For the DLM, our GAN parameters (e.g., loss function, VGG layer, discriminator weight, and mini-batch size) should be optimized to improve the image and dose accuracies. Moreover, the DL network was trained with 2D axial slices and not with fully 3D images because of memory GPU limitations. Finally, DL architectures such as the cycleGAN, which may have enabled some intra-individual co-registration issues to be overcome, could be further investigated.

Conclusions

Accepted Article

Four methods of H&N CBCT-based dose calculation were compared in this study. The use of an HU-D curve from CBCT was the simplest but the least accurate method. The DAM, use of DIR, and GAN DLM provided similar dose discrepancies. The DIR method and DLM appeared to be the most attractive methods in terms of dose accuracy as well as calculation time. Using one of such a method with CBCT images would enable dose monitoring in the PGs during the H&N treatment course to trigger replanning.

Figure captions

Accepted Article

The systematic dose discrepancy was defined as the mean DVH differences between the reference CT and the pCTs generated by each method. Permutation tests were performed to compare the DVH differences of each method to a null distribution. Significant differences (p ≤ 0.05) are displayed at the top of each figure using the symbol *. In the two last columns, "+" represents the qualitative comparison of the four methods in terms of image and dose endpoints. The more the number of "+", the more accurate is the method. The dose endpoints are the DVH and gamma results.

Accepted Article

  ii.Density assignment method (DAM)Anautomatic thresholding was performed inside the patient body contour for all the CBCT and planning CT images based on HU values to create bone and air cavity classes. Manual corrections of the segmentation were necessary in case of dental artefacts and shoulder areas in the CBCT images. A soft tissue class was then created by the subtraction of air and bone classes within the patient body contour. The HU values of each class from the mean values of the planning CT were assigned to the corresponding classes of the CBCTs for each patient. Hence, three classes were obtained with the following values: -730 HU for the air cavity, 0 HU for the soft tissue, and 550 HU for the bone.

  voxel-wise comparison of the HU between CT ref and pCT was performed. The mean absolute error (MAE) and the mean error (ME) were calculated between the CT ref and pCT obtained from

Fig. 1 .

 1 Fig. 1. Deep learning method for pseudo-CT generation from CBCT The training of the generative adversarial network (GAN) comprises two competing multilayer networks: the generator and the discriminator. The discriminator aims to distinguish the real image (ground truth) from the realistic fake image (pCT) produced by the generator. The input data of the generator are CBCT and CT images that provide pCTs. The discriminator classifies these pCTs as real or fake CTs until the discriminator cannot determine whether the pCT is a real CT. In the testing step, for a new test patient, the CBCT goes through the trained network to obtain the corresponding pCT.

Fig. 2 .Fig. 3 .

 23 Fig. 2. CBCT, reference CT, and pseudo-CT images from the deep learning method

Fig. 4 .

 4 Fig. 4. Dose discrepancies for CBCT-based dose calculation methods along the entire DVH for the PTV, PGs, and spinal cord The dose discrepancy is defined as the mean absolute DVH differences between the reference CT and the pCT corresponding to each method. Permutation tests were performed to compare the absolute DVH differences of the DLM to those of the other methods. Significant differences (p ≤ 0.05) are displayed at the top of each figure with *.

Fig. S- 1 .

 1 Fig. S-1. Convergence of generator and discriminator loss functions on the training cohort

Fig. S- 3 .

 3 Fig. S-3. Illustration of pseudo-CTs, dose distributions, and gamma maps for one patient in coronal planeThe dose distributions are displayed for each of the four investigated methods. The dose differences are illustrated in terms of gamma analysis. 3D gamma pass-rate values corresponding to each method are indicated for this patient.Accepted Article

Table 2

 2

lists the imaging endpoints for the entire body contour, soft tissue, and bone. In the entire body, the MAE and ME were 266.6 and 208.9 HU, 113.2 and 14.2 HU, 95.5 and -36.6 HU, and 82.4 and 17.1 HU for the HU-D curve method, DAM, DIR method, and DLM, respectively. The MAE results of the DLM differed significantly from those of other methods (p < 0.05). The ME results of the four methods differed significantly from the null distribution.

Table 2 . Imaging endpoints comparing the reference CT to the pseudo-CTs obtained by each method MAE

 2 : mean absolute error of HU values defined as the mean difference (in absolute value) of HU values per voxel between the reference CT and pseudo-CT; ME: mean error, defined as the mean difference of HU values per voxel between the reference CT and the pseudo-CT.The imaging endpoint values are expressed as mean ± standard deviation.The Wilcoxon test was used to first compare the MAE of the DLM to those of the other methods and then to compare the ME of the methods to a null distribution. Significant differences (p ≤ 0.05) are displayed using the symbol *.

	Entire body Soft tissue only Bone only Accepted Article	Endpoints MAE (HU) ME (HU) MAE (HU) ME (HU) MAE (HU) ME (HU)	HU-D curve method 266.6* ± 25.8 208.9* ± 36.1 260.9* ± 25.9 206.4* ± 35.7 344.6* ± 64.7 253.1* ± 127.2	Density assignment method 113.2* ± 6.7 14.6* ± 11.2 91.2* ± 6.4 9.6* ± 11.0 388.1* ± 38.3 84.7* ± 77.6	DIR method 95.5* ± 21.2 -36.6* ± 14.7 85.0* ± 20.2 -33.3* ± 14.7 226.2* ± 44.4 -77.4* ± 28.6	Deep learning method 82.4 ± 10.6 17.1* ± 19.9 69.2 ± 15.3 12.9* ± 19.4 207.6 ± 41.8 64.9* ± 56.1

Acknowledgments: With financial support from ITMO Cancer AVIESAN (Alliance Nationale pour les Sciences de la Vie et de la Santé/ National Alliance for Life Sciences & Health) within the framework of the Cancer Plan. ADMIRE software was supplied by Elekta AB, Stockholm, Sweden as part of a research agreement.

Conflict of Interest Statement:

The authors have no relevant conflicts of interest to disclose.

Accepted Article

Accepted Article