PhD Adam Trickey
email: adam.trickey@bristol.ac.uk

Kathleen Mcginnis

John Gill

Sophie Abgrall

Juan Berenguer

Christoph Wyen

Mojgan Hessamfar

Peter Reiss

Katharina Kusejko

PhD Michael J Silverberg

PhD Suzanne M Ingle

Prof J A C Sterne Phd

Arkaitz Imaz

Ramon Teira

Antonella D'

Arminio Monforte

Robert Zangerle

Prof Jodie L Guest

Vasileios Papastamopoulos

Heidi Crane

Timothy R Sterling

Sophie Grabar

Jonathan A C Sterne

Longitudinal trends in causes of death among adults with HIV on antiretroviral therapy in Europe and North America from 1996 to 2020: a collaboration of cohort studies

des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

Introduction

Before the introduction of combination antiretroviral therapy (ART) in 1996, people with HIV experienced very high mortality rates, mostly due to AIDS. [START_REF] Krentz | Changing mortality rates and causes of death for HIV-infected individuals living in Southern Alberta, Canada from 1984 to 2003[END_REF] Since then, the vast majority of people with HIV in high-income countries are now on ART. 2 The successful suppression of viral replication with ART has led to substantial reductions in the risk of AIDS and death and to corresponding large increases in life expectancy. 3 With reduced AIDS-related mortality, people with HIV are ageing and experiencing an increased mortality due to other age-related causes, such as cardiovascular disease and cancer. [START_REF] Althoff | The shifting age distribution of people with HIV using antiretroviral therapy in the United States[END_REF][START_REF] Marty | Projection of age of individuals living with HIV and time since ART initiation in 2030: estimates for France[END_REF] Substance use and comorbidities such as hepatitis C virus (HCV) infection are more common in people with HIV than in the general population, leading to increased rates of mortality related to these causes. [6][START_REF] Alejos | Overall and cause-specific excess mortality in HIV-positive persons compared with the general population: role of HCV coinfection[END_REF][START_REF] Ingle | Impact of risk factors for specific causes of death in the first and subsequent years of antiretroviral therapy among HIV-infected patients[END_REF] Between 2012 and 2015, following the INSIGHT START Study, [START_REF] Lundgren | Initiation of antiretroviral therapy in early asymptomatic HIV infection[END_REF] international treatment guidelines were changed to recommend that all people with HIV receive ART regardless of disease stage and CD4 cell count. [START_REF] Eholie | Antiretroviral treatment regardless of CD4 count: the universal answer to a contextual question[END_REF] Accordingly, the time between diagnosis and start of ART was reduced. Additionally, more potent and well tolerated ART regimens, particularly second-generation integrase inhibitors, have continued to become available and are now widely used, [START_REF] Vitoria | Current and future priorities for the development of optimal HIV drugs[END_REF] leading to ongoing reductions in AIDS incidence and mortality. [START_REF] Brooks | Integrase inhibitors: after 10 years of experience, is the best yet to come?[END_REF] Policy makers need current data to understand the changing burden of cause-specific mortality among people with HIV, both overall and in subgroups defined by age, sex, and mode of acquisition of HIV.

The Antiretroviral Therapy Cohort Collaboration (ART-CC) [START_REF] May | Cohort profile: Antiretroviral Therapy Cohort Collaboration (ART-CC)[END_REF] has classified causes of death among people with HIV living in Europe and North America since 1996. We investigated longitudinal trends in rates of cause-specific mortality among adult people with HIV who started ART in Europe and North America between 1996 and 2020.

Methods

Study design and population

In this collaborative HIV cohort study, we combined data from 17 European and North American HIV cohort studies participating in the ART-CC for which at least 70% of deaths had been classified by cause. The appendix shows details on the included cohorts (appendix p 1). Ethics committees or institutional review boards approved the individual cohorts, in which standardised data collection methods were used and patients were followed up regularly. Cohorts included information on mortality that was gathered through linkage with vital statistics agencies and hospitals or physician report and on the active follow-up of participants. Eligible participants were aged 16 years or older when starting combination ART (three drugs or more) and had no previous exposure to monotherapy or dual ART regimens. Included participants had a CD4 count and HIV-1 RNA viral load measurement within a window of 1 month before and 1 week after starting ART.

Procedures

We adapted the Cause of Death (CoDe) project protocol [START_REF] Kowalska | The Coding Causes of Death in HIV (CoDe) project initial results and evaluation of methodology[END_REF] to classify causes of death information into a single cause in the HICDEP format. Information on cause of death was recorded either as International Classification of Diseases, Ninth Revision (ICD-9) or Tenth Revision (ICD-10), or free text. If ICD-9 or ICD-10 codes were available, causes of death were classified by a clinician and a computer algorithm. [START_REF] Kowalska | The Coding Causes of Death in HIV (CoDe) project initial results and evaluation of methodology[END_REF] When ICD-9 or ICD-10 codes were not available, two clinicians independently classified each death. Disagreements between clinicians or computer-assigned codes were resolved via panel discussion.

Clinicians classified deaths using summary tables of patient data including ICD-9 or ICD-10 codes or free text for cause of death, patient characteristics at ART initiation (age, sex, transmission risk group, AIDSdefining conditions, and HCV antibody status), use of ART before death, AIDS-defining conditions after starting ART, and latest CD4 cell count within 6 months of death (appendix pp 2-4). Deaths were coded as AIDS-related if there was a serious AIDS-defining condition before death or low CD4 count (<100 cells per μL) within a year of death (18 months if off treatment), and a diagnosis compatible with AIDS as cause of death.

Statistical analysis

Data were split into calendar periods of follow-up (1996-99, 2000-03, 2004-07, 2008-11, 2012-15, and 2016-20). We derived participants' characteristics at the time of starting ART and at the beginning of each subsequent period. The last follow-up date was April, 2020. People were censored at the first of a cohortspecific administrative censoring date, date of death, or date of loss-to-follow-up (defined as a period of at least a year between the person's last visit and the cohortspecific administrative censoring date, with the date of loss-to-follow-up defined as 6 months after this last visit date). When a participant was lost-to-follow-up but then found to have died via linkage to death registry, their follow-up was instead censored at date of death. We did a sensitivity analysis excluding deaths that happened after loss to follow-up. As ART should be continued once started, we did not consider ART discontinuation in the analyses.

We used Poisson regression to estimate cause-specific mortality rate ratios (MRRs) for each follow-up calendar period compared with the 2000-03 period, as well as the average MRR per period assuming that changes across periods were log-linear. The 2000-03 period was chosen as the comparator as all cohorts had data for this period, whereas two cohorts had limited data for the 1996-99 period. The Poisson models accounted for each individual's exposure time (in days) during each calendar year period of follow-up. These periods began on Jan 1, 1996, 2000, 2004, 2008, 2012, and 2016. Estimated MRRs were adjusted for time-updated age (16-29, 30-39, 40-49, 50-59, and ≥60 years), timeupdated CD4 count (0-199, 200-349, 350-499, ≥500 cells per μL, or missing), whether the person was ART-naive at the start of each period, and cohort. To obtain time-updated CD4 counts, we used a window from 365 days before the start of that follow-up period to 30 days after. In the categorised time-updated CD4 count variable, we included a separate missing data category so that people with missing data were not excluded. We estimated MRRs overall and stratified into subgroups based on time-updated age (16-39, 40-59, and ≥60 years), time-updated CD4 count (0-199, 200-349, 350-499, ≥500 cells per μL, and missing), whether ARTnaive at the start of each period (to capture the period at the start of ART when mortality rates are highest), timeupdated previous AIDS-defining event status (no, yes; defined with a HICDEP code), time-updated HCV antibody status (negative or positive), cohort region (Europe or North America), and self-reported sex and HIV acquisition risk group (men who have sex with men [MSM], men who acquired HIV through injecting drug use [IDU], women who acquired HIV through IDU, men who acquired HIV through sex with women, women who acquired HIV through sex with men, and either sex with risk group other or unknown [people in this final category are mostly from the Veterans Aging Study Cohort]). Men and women with unknown or other mode of HIV acquisition were combined into one subgroup for the cause-specific mortality analyses.

Due to data sharing rules regarding cause-of-death information, analyses for one of the cohorts were done separately on the basis of deaths classified using the second computer algorithm. MRRs and their confidence intervals were produced separately for this cohort (in Yale, CT, USA) and for the other cohorts (in Bristol, UK). To combine them, their natural logs were taken and then combined by use of inverse-variance weighted metaanalysis (in Bristol, UK) using the metan command in Stata (version 17.1) on the log-transformed rate ratios and their standard errors. Estimated MRRs for 1996-99 versus 2000-03 exclude this cohort, for which data for the 1996-99 period were unavailable. All other estimated MRRs are based on all cohorts.

In a sensitivity analysis, we removed the two cohorts for which the second computer algorithm was used For the codes for AIDS-defining events see https://hicdep.org/ Wiki/v/10/pt/4/Table/94/ FieldID/1206 instead of clinician input. An additional, post-hoc analysis was done to investigate longitudinal changes in rates of cause-specific deaths among people who acquired HIV through IDU among those living in Europe and those living in North America separately. We estimated per-period MRRs for each cause of death adjusted for time-updated age, CD4 count, whether the person was ART-naive at the start of each period, and cohort. All analyses were done with Stata (version 16.1).

Role of the funding source

The funder of the study had no role in study design, data collection, data analysis, data interpretation, or writing of the report.

Results

16 832 (8•9%) deaths were recorded among 189 301 people with HIV followed up for 1 519 200 person-years. The crude all-cause mortality rate was 11 ) had AIDS when starting ART, and 14 967 (7•9%) had tested positive for HCV antibodies before starting ART. Time-updated CD4 count at the start of a new period was missing for some people, but the proportion of people with missing data was not large (eg, 7•5% for 2016-20). The category with most missing data was HCV antibody testing, but the proportion of people with missing HCV antibody data improved over time from 50•4% in 1996-99 to 11•1% in 2016-20.

The cause of death could not be classified for 3652 (21•7%) of the 16 832 deaths, a proportion that remained quite consistent across the calendar year periods. The most common causes of death were AIDS (4203 deaths; 25•0%), non-AIDS non-hepatitis malignancy (2311 deaths; 13•7%) and cardiovascular or heart-related (1403 deaths; 8•3%). There was a large reduction in the proportion of AIDS-related mortality, from 49% during 1996-99 to 16% during 2016-20, with increases in the proportion of mortality due to cancers from 5% during 1996-99 to 19% during 2016-20 (figure 1). There was a reduction in the proportion of AIDS-related deaths over calendar time for all age groups (figure 2): the proportion of AIDS-related deaths was lowest for the oldest age group who, correspondingly, had the highest proportions of cardiovascular or heartrelated and cancer-related mortality.

Rates of all-cause mortality declined over time (average adjusted MRR [aMRR] per period 0•85, 95% CI 2). The largest average per-period reductions were for AIDS-related mortality and cardiovascular or heart-related mortality.

Rates of all-cause mortality declined more steeply over calendar time in people with HIV aged 16-39 years than in those aged 40-59 years or those aged 60 years or older (table 2). Similar patterns were seen for AIDS-related and cardiovascular or heart-related mortality. For liver mortality and for non-AIDS, non-hepatitis malignancy mortality, the steepest average declines were in the youngest age group, with little evidence of a decline among people with HIV aged 60 years or older. There were larger reductions in all-cause mortality over calendar time among people with HIV who had started ART in the current period than among those who had started ART during a previous period (table 2). This pattern was also seen for most categories of causespecific mortality.

Rates of all-cause mortality declined over calendar time for MSM and for both men and women who acquired HIV through heterosexual sex; the decline was smaller for men who acquired HIV through IDU; and there was an increase for women who acquired HIV through IDU. For women who acquired HIV through IDU, there were increases over calendar time in respiratory-related mortality. For MSM, there was a large decrease over calendar time in rates of cardiovascular or heart-related mortality.

People with HIV with the highest CD4 counts experienced larger average declines in rates of all-cause mortality than those with lower CD4 counts. For many categories of cause-specific mortality, there were also larger average declines in mortality rates among people with HIV with higher CD4 counts than among those with lower CD4 counts. Declines in rates of all-cause mortality over calendar time were larger among people with HIV with no previous AIDS-defining events than for those with previous AIDS-defining events, with similar patterns observed for most causes of death. Rates of all-cause mortality decreased over calendar time in people with HIV without HCV, but the decreases were less steep among people with HIV with HCV.

Declines in all-cause mortality rates were less substantial for people with HIV living in Europe than for those living in North America. For most categories of cause-specific mortality, declines in rates of mortality over calendar time were greater in North America than in Europe. However, rates of substance use-related mortality declined on average among people with HIV in Europe but increased among people with HIV in North America.

In a sensitivity analysis omitting the two cohorts coded by the second algorithm (containing around half of the deaths), the main differences were that there were reductions in rates of CNS-related mortality over calendar time (average aMRR per period 0•75, 95% CI 0•59-0•95) and substance use-related mortality (0•88, 0•79-0•98), which were not seen when including all cohorts (appendix p 6). These differences might possibly have been due to one of these cohorts contributing little data in the 1996-99 period. Results of the sensitivity analysis excluding deaths that happened after loss to follow-up are shown in the appendix (p 7). The average aMRRs per period were similar to those in the analysis including these deaths, except for other causes of death for which the aMRR was 0•87 (0•83-0•91) when excluding the deaths compared with 0•81 (0•77-0•86) when including them.

In post-hoc analysis comparing MRRs for causespecific mortality among people with HIV who acquired HIV through IDU living in Europe with those living in North America, most causes of death had wide confidence intervals, and therefore it was difficult to ascertain differences between people with HIV in Europe and those in North America (appendix p 8).

Discussion

There have been reductions in rates of all-cause mortality and of most major causes of death among people with HIV on ART between 1996 and 2020. For the overall population of people with HIV, reductions in mortality rates over calendar time were seen for each cause of mortality except for CNS, respiratory, and substance use deaths. The biggest reductions were for AIDS-related and cardiovascular or heart-related mortality. The *Adjusted for CD4 cell count category and age group at the start of the period, whether the participant was ART-naive when starting the period, and cohort. †Calendar year period is modelled as a continuous variable; the MRRs can be interpreted as a per-period decrease. reductions in mortality rates over calendar time differed between subgroups of people with HIV. Declines in heart-related or cardiovascular mortality were most marked among MSM, and they might reflect improvements in post-cardiovascular disease care in the same periods among the general population or perhaps the reduced toxicity of ART regimens used in the later studied time periods. For men who acquired HIV through IDU, decreases in mortality rates were smaller than for other groups, whereas mortality rates increased in women who acquired HIV through IDU. For most causes of death, reductions in mortality rates were larger in North America than in Europe. However, there were reductions in substance use-related mortality over calendar time among people with HIV living in Europe and increases in North America.

0•7 0•8 0•9 1•0 1•1
Other published studies have investigated changes in cause-specific mortality over time among people with HIV on ART. A Spanish population-based cohort of people with HIV, including those not on ART, showed that all-cause mortality decreased from 1999 to 2018, with the decrease driven by AIDS-related mortality, while non-AIDS mortality remained stable. [START_REF] Fontela | Trends and causes of mortality in a population-based cohort of HIV-infected adults in Spain: comparison with the general population[END_REF] The largest decrease in mortality was among people with HIV who had had AIDS-defining events. HIV-related mortality decreased among people with HIV without a previous AIDSdefining event, but the number of such deaths was low. [START_REF] Fontela | Trends and causes of mortality in a population-based cohort of HIV-infected adults in Spain: comparison with the general population[END_REF] A study among people with HIV on ART in British Columbia, Canada, observed reductions in all-cause, AIDS-related, cardiovascular, liver-related, and suiciderelated mortality over calendar time from 2001 to 2012, 16 while the HIV Atlanta Veterans Affairs Cohort Study found declines in all-cause, AIDS-related, and non-AIDSrelated mortality rates during the ART era. [START_REF] Vyas | Trends in cause-specific mortality among veterans with HIV: a 35-year (1982-2016) analysis of the HIV Atlanta VA Cohort Study[END_REF] A previous ART-CC study examined changes from 1996 to 2015 in cause-specific mortality by ART start year, which is comparable to the changes in the present study among people with HIV who were ART-naive in each follow-up period. 3 That study also observed reductions over time in AIDS-related mortality, non-AIDS infection, non-AIDS cancers, cardiovascular, suicide or accident-related, and other mortality. 3 Our finding of an increase in substance use-related mortality in North America and a decrease in Europe perhaps reflects an increase in opioid use in the USA, [START_REF] Singh | Opioid epidemic in the United States: empirical trends, and a literature review of social determinants and epidemiological, pain management, and treatment patterns[END_REF] where opioid-related mortality is over ten times higher than in the European Union. [START_REF] Kalkman | Monitoring opioids in Europe: the need for shared definitions and measuring drivers of opioid use and related harms[END_REF] Unlike in the USA, [START_REF] Bradley | Estimated number of people who inject drugs in the United States[END_REF] there have been reductions in the prevalence of injecting drug use in some European countries, [START_REF] Degenhardt | Global prevalence of injecting drug use and sociodemographic characteristics and prevalence of HIV, HBV, and HCV in people who inject drugs: a multistage systematic review[END_REF] although not necessarily opiate use. [START_REF] Kalkman | Monitoring opioids in Europe: the need for shared definitions and measuring drivers of opioid use and related harms[END_REF] Additionally, opiate substitution therapy coverage is higher in western European countries than in North America, which might contribute to this difference. 22 Our finding that rates of respiratory-related mortality increased among people with HCV is likely due to the high prevalence of HCV among people who acquired HIV through IDU, among whom there was also an increase in rates of respiratoryrelated mortality.

The main strengths of this study are the large sample size, geographical diversity, representativeness of the included people with HIV, and the availability of data on cause-specific mortality, classified according to a common protocol. [START_REF] Kowalska | The Coding Causes of Death in HIV (CoDe) project initial results and evaluation of methodology[END_REF] However, causes of death were classified retrospectively and without complete patient histories, so misclassification is more likely than if the causes had been classified on the basis of full medical history. Autopsy, which is becoming less common over time, [START_REF] Weber | Decreasing mortality and changing patterns of causes of death in the Swiss HIV Cohort Study[END_REF] remains the gold standard for classifying causes of death, and clinical classifications might not correlate well with those from autopsy reports. [START_REF] Cox | Autopsy causes of death in HIV-positive individuals in Sub-Saharan Africa and correlation with clinical diagnoses[END_REF][START_REF] Roulson | Discrepancies between clinical and autopsy diagnosis and the value of post mortem histology; a meta-analysis and review[END_REF] For some people with HIV there was information on a number of comorbidities or potential causes of death, which could lead to the death being classified as of unknown cause or unclassifiable. It is possible that changes in how causes of death are captured on death certificates have changed over time. In particular, HIV was previously more likely to be mentioned as a cause on a death certificate because the person who died had HIV, even if it was irrelevant to the death. The CoDe processes and rules were set up to minimise the impact of this, by accounting for recent CD4 counts and AIDS diagnoses rather than relying solely on death certificate information. For our models of cause-specific mortality, we adjusted for time-updated CD4 cell counts and age to capture major changes in HIV prognostic markers for mortality. Data for some established risk factors for chronic diseases, particularly smoking and alcohol consumption, [START_REF] Petoumenos | Smoking, alcohol and illicit drug use effects on survival in HIV-positive persons[END_REF] were not available for all patients and therefore we were unable to estimate the contributions of these risk factors to causespecific mortality. Other potentially important patient characteristics, including socioeconomic factors such as education, poverty, and homelessness, were not routinely collected across the cohorts, so we were unable to include them in our analyses. Data for HIV acquisition risk group were available for all but one cohort, while data for one cohort were unavailable in 1996-99. The cohorts are from various countries with different health systems, population characteristics, and migration patterns, so pooling such data might obscure within-country patterns, so we conducted analyses stratified by region. Loss-tofollow-up was high in some cohorts with our definition of a gap of at least a year between the person's last visit and the cohort-specific administrative censoring date, particularly for hospital-based cohorts in which a large percentage of patients subsequently moved to generalpractitioner-based care (appendix p 1). However, when a participant was lost-to-follow-up and information on death was later available through linkage to registries we censored the patient at the date of death, so that in such situations loss-to-follow-up would not affect ascertainment of deaths. Results of sensitivity analyses excluding deaths after loss to follow-up were similar to the results of the main analysis. Most of our findings are likely to be generalisable to other high-income countries such as Australia and Japan, although some particular trends might be context-specific, for example regarding substance use deaths.

The trends in cause-specific mortality captured here should assist policy makers in targeting improvements in the care of people with HIV towards conditions that are amenable to interventions and have the biggest influence on mortality. The per-period reductions in cause-specific mortality for people with HIV on ART are probably partly due to changes in treatment guidelines meaning that people start ART sooner after diagnosis, [START_REF] Gunthard | Antiretroviral drugs for treatment and prevention of HIV infection in adults: 2016 recommendations of the International Antiviral Society-USA Panel[END_REF] as well as more effective and less toxic regimens becoming available, [START_REF] Vitoria | Current and future priorities for the development of optimal HIV drugs[END_REF] and better care in general for people with HIV. [START_REF] Mugavero | Elements of the HIV care continuum: improving engagement and retention in care[END_REF] Although mortality among people with HIV on ART has decreased, 3 there is still higher mortality among people with HIV than among the general population, due both to the consequences of HIV infection and to a higher prevalence of comorbidities and risk behaviours among people with HIV. [START_REF] Fontela | Trends and causes of mortality in a population-based cohort of HIV-infected adults in Spain: comparison with the general population[END_REF][START_REF] Croxford | Mortality and causes of death in people diagnosed with HIV in the era of highly active antiretroviral therapy compared with the general population: an analysis of a national observational cohort[END_REF][START_REF] Trickey | Parameter estimates for trends and patterns of excess mortality among persons on antiretroviral therapy in high-income European settings[END_REF] Although there were reductions in rates of non-AIDS-related mortality, such as cancer and cardiovascular disease, non-AIDS deaths make up an increasingly large proportion of mortality among people with HIV. Expanding access to prevention, screening, and treat ment of these conditions is required to close the gap of comorbidity prevalence between adults with and those without HIV, and funders should recognise this. Progress in reducing cause-specific mortality has not been evenly spread across subgroups of people with HIV, often with the most marginalised populations experiencing the least benefits. Men who acquired HIV through IDU had the lowest reductions in mortality, while there was some evidence of increases in mortality among women who acquired HIV through IDU. People with histories of substance use conditions have higher rates of homeless ness and other comorbidities [START_REF] Hotton | Trends in homelessness and injection practices among young urban and suburban people who inject drugs: 1997-2017[END_REF][START_REF] Lim | Infectious diseases, comorbidities and outcomes in hospitalized people who inject drugs (PWID) infections in persons who inject drugs[END_REF] and often experience additional barriers and stigma while attempting to access care. [START_REF] Paquette | Stigma at every turn: health services experiences among people who inject drugs[END_REF] This indicates that targeted interventions, such as addressing social determinants of health and bringing comorbidity care to needle and syringe dispensing locations, are required for people who acquired HIV through IDU. That there was a decrease in substance abuse-related mortality in Europe, but an increase in North America, which saw a large increase in opioid use during the study period, [START_REF] Singh | Opioid epidemic in the United States: empirical trends, and a literature review of social determinants and epidemiological, pain management, and treatment patterns[END_REF] shows that contextspecific analyses are also required to understand the epidemiology of mortality among people with HIV in different settings and that the most impactful interventions will probably vary by location.

Contributors

AT and JACS conceived and designed the study. AT and SMI combined, checked, and cleaned the datasets. MJG, SA, JB, CW, MH, and PR were involved in the cross-cohort coding of causes of death. AT conducted the statistical analyses. KMcG performed separate analyses that were then pooled via meta-analysis. AT and JACS drafted the report. All authors contributed to the interpretation of data and critical revisions of the report for important intellectual content. AT and SMI accessed and verified the data underlying this study. AT and JACS were responsible for the decision to submit for publication.

Declaration of interests

AI has received financial compensation for lectures, educational activities, consultancy work, as well as funds for research, from Gilead Sciences, Janssen-Cilag, Merck Sharp & Dohme, and ViiV Healthcare. VP has received honoraria from ad-hoc membership of national HIV advisory boards, Merck, Gilead, and ViiV. PR, through his institution, has received scientific grant support for investigator-initiated studies from Gilead Sciences, Janssen Pharmaceuticals, Merck & Co, and ViiV Healthcare, and has served on scientific advisory boards for Gilead Sciences, ViiV Healthcare, and Merck & Co, honoraria for which were all paid to his institution. JB reports honoraria for advice or public speaking from Gilead, GlaxoSmithKline, Janssen, MSD, and ViiV Healthcare; and grants from Gilead, MSD, and ViiV Healthcare. MJG has received honoraria from ad-hoc membership of national HIV advisory boards, Merck, Gilead, and ViiV. HC has received research grant funding from ViiV Healthcare, National Institutes of Health (NIH), and Agency for Healthcare Research and Quality paid to their institution and sits on the NIH Office of AIDS Research Advisory Council. CW reports honoraria for advice or public speaking from Abbott, Gilead, Janssen, MSD, Pfizer, and ViiV Healthcare. KK, TRS, SMI, SG, SA, RT, RZ, MH, MJS, JACS, AT, JLG, KMcG, and Ad'AM declare no competing interests.

Data sharing

Due to the data sharing agreements between individual cohorts and the Antiretroviral Therapy Cohort Collaboration (ART-CC), the data collected for this study cannot be shared. Data are owned by the individual cohorts and those wishing to access these data should contact the individual cohorts, details of which are given in the appendix.

Figure 1 :Figure 2 :

 12 Figure 1: Proportions of each cause of death category among people with HIV who died, by calendar year period of death

Figure 3 :

 3 Figure 3: Adjusted* cause-specific MRRs per period (1996-99, 2000-03, 2004-07, 2008-11, 2012-15, 2016-2020 †), with 95% CIs CNS category includes Parkinson's and Alzheimer's disease. MRR=mortality rate ratio.*Adjusted for CD4 cell count category and age group at the start of the period, whether the participant was ART-naive when starting the period, and cohort. †Calendar year period is modelled as a continuous variable; the MRRs can be interpreted as a per-period decrease.

table 1 )

 1 •1 deaths . The median age at initiation of ART increased from 35 years (IQR 31-42) during 1996-99 to 38 years (30-48) during 2016-20. Participants' median age increased from 37 years (32-43) on Jan 1, 2000, to 43 years (36-50) on Jan 1, 2008, and 47 years (38-54) on Jan 1, 2016. The proportion of people with HIV with CD4 counts of 0-199 cells per μL dropped from 42•6% at the start of the 1996-99 period to 7•6% at the start of the 2016-20 period. Of the 189 301 people with HIV included in the study, 44 348 (23•4%) were women and 144 953 (76•6%) were men, 36 311 (19•2%

		Overall, when	Number of	Number of participants alive during each calendar period
		starting ART	deaths
		(n=189 301)	(n=16 832)
				1996-99	2000-03	2004-07	2008-11	2012-15	2016-20
				(n=20 448)	(n=47 124)	(n=74 386)	(n=108 946)	(n=139 851)	(n=145 996)
	Age
	16-39 years	103 112 (54•5%)	5728 (34•0%)	13 969 (68•3%)	28 242 (59•9%)	37 463 (50•4%)	46 792 (42•9%)	54 166 (38•7%)	48 140 (33•0%)
	40-59 years	75 731 (40•0%)	8724 (51•8%)	5874 (28•7%)	17 073 (36•2%)	33 049 (44•4%)	54 489 (50•0%)	72 899 (52•1%)	79 989 (54•8%)

www.thelancet.com/hiv Published online January 24, 2024 https://doi.org/10.1016/S2352-3018(23)00272-2

www.thelancet.com/hiv Published online January 24, 2024 https://doi.org/10.1016/S2352-3018(23)00272-2

Acknowledgments

The study was funded by the US National Institute on Alcohol Abuse and Alcoholism.

Funding US National Institute on Alcohol Abuse and Alcoholism.


Marc Briane
email: mbriane@insa-rennes.fr

Some classes of singular ODEs flows

Keywords: ODEs flow, rotation set, invariant measure, first integral Mathematics Subject Classification: 34E05, 34E10, 37C10, 37C40

Some classes of singular ODEs flows

Introduction

In this paper we study some ODEs flows associated with a regular vector field b defined in the torus T d , solutions to

   ∂ t X(t, x) = b Ä X(t, x) ä , t ∈ [0, ∞) X(0, x) = x ∈ T d . (1.1)
A flow X will be said singular if it has one or both of the following properties:

1. The flow X admits only invariant probability measures on T d ( 1 ) which are singular with respect to Lebesgue's measure. For example, consider the Stepanoff flow [START_REF] Stepanoff | Sur une extension du théorème ergodique[END_REF] associated with the vector field defined by b(x) := ρ(x) ξ, where ρ is a non-negative function in C 1 (T d ) with a finite number of roots in T d and σ := 1/ρ ∈ L 1 (T d ) and where ξ ∈ R d has incommensurable coordinates. By virtue of [START_REF] Briane | Asymptotics of ODE's flow on the torus through a singleton condition and a perturbation result. Applications[END_REF]Pro. 5.1]) σ(x)/σ dx is the unique invariant probability measure for this flow which is absolutely continuous with respect to Lebesgue's measure. In contrast, if the function ρ changes sign and has a (at most) countable zero set {ρ = 0}, then by virtue of [5, Pro. 2.1] any invariant probability measure on T d has support in {ρ = 0}, and is thus singular with respect to Lebesgue's measure.

1 2. The flow X has no first integral u (2 ) with Z d -periodic gradient, which is real-analytic. From now one, such a flow will be said not integrable. For example, the so-called ABC flow (Arnold-Beltrami-Childress) introduced by Arnold [START_REF] Arnold | Sur la topologie des écoulements stationnaires des fluides parfait[END_REF] is shown in [START_REF] Ziglin | On the absence of a real-analytic first integral for ABC flow when A=B[END_REF] to be not integrable when A = B.

In view of the previous (apparently) simple examples the singular properties of some given flow are difficult to derive. So, the dynamical systems of classical mechanics (see, e.g. [8, §2.2]) have a smooth invariant probability measure (with a constant density for the Hamiltonian systems) by virtue of Liouville's theorem [8, Thm. 1, §2.2] (see also [START_REF] Briane | A picture of the ODE's flow in the torus: from everywhere or almost-everywhere asymptotics to homogenization of transport equations[END_REF]Pro. 2.2]). In the present case, we need thus to prove that any invariant probability measure µ on T d for the flow (1.1) associated with some vector field b, i.e. solution to the PDE div (µ b) = 0 in T d (1.2)

(see below for an equivalent variational formulation (1.10)), is actually singular with respect to the Lebesgue measure. However, such a solution µ is far to be explicit, which makes the problem quite intricate.

Our aim is to provide different classes of ODEs flows which have only singular invariant probability measures on T d for d = 2, 3, and which are possibly not integrable. First, in Section 2 we prove a general two-dimensional result for a flow (1.1) (see Theorem 2.1) the assumptions of which combine the structure of the vector field b and the one of the Herman rotation set C b (see definition (1.14) below) induced by b. More precisely, we assume that the set C b -which by virtue of the Franks-Misiurewicz theorem (see Theorem 1.1) is known to a closed line segment along a line passing through 0 R 2 -is contained in the line R k for some k ∈ Z 2 \ {0 R 2 }. We also assume that the vector field b satisfies the representation

∀ x ∈ T 2 , b(x) • k ⊥ = α(x • k ⊥ ) + β(x),
with max

x∈T 2 α(x • k ⊥ ) > max T 2 |β|. (1.3)
In other words, (1.3) means that the component b(x)•k ⊥ is dominated by a function which only depends on the variable x • k ⊥ and which changes sign. Then, the flow associated with b has no regular invariant probability measure of the type σ(x) dx, where σ is a positive (in the broad sense of (2.3)) function in L 1 (T 2 ). If in addition the vector field does not vanish in T 2 , then the flow has no real-analytic first integral u such that ∇u is Z 2 -periodic. Such a class is illustrated by the flows of Proposition 2.4 (in which x • k ⊥ = x 1 ) which satisfy the two previous conditions. In particular, this class contains the two-dimensional parametrized Euler flows (see §2.

Notations

• (e 1 , . . . , e d ) is the canonical basis of R d , and 0 R d is the null vector of R d .

• I d is the unit matrix of R d×d , and

R ⊥ := Å 0 -1 1 0 ã .
• A T denotes the transposed matrix of A ∈ R d×d , and R d×d S denotes the subspace of the symmetric matrices in R d×d .

• tr A denotes the trace and det A denotes the determinant of the matrix A ∈ R d×d .

• For ξ = (ξ 1 , ξ 2 ) ∈ R 2 , we denote ξ ⊥ := R ⊥ ξ = (-ξ 2 , ξ 1 ).
• " • " denotes the scalar product and | • | the euclidean norm in R d .

• T d denotes the d-dimensional torus R d /Z d for d ∈ N.
• dx denotes the Lebesgue measure on R d and on T d .

• M (T d ) denotes the set of Radon measures on the torus T d , and M p (T d ) denotes the set of probability measures on T d .

• C k (T d ), respectively W 1,k (T d ), k ∈ N ∪ {∞}, denotes the space of the real-valued functions f ∈ C k (R d ), respectively f ∈ W 1,k (R d ), which are Z d -periodic, i.e. ∀ k ∈ Z d , ∀ x ∈ R d , f (x + k) = f (x). (1.5)
• For any function f ∈ L 1 (T d ), we denote

         f x (x 1 ) := ˆTd-1 f (x 1 , x ) dx , x 1 ∈ T 1 f x 1 (x ) := ˆT1 f (x 1 , x ) dx 1 , x = (x 2 , . . . , x d ) ∈ T d-1 ,
where f :=

ˆTd f (x) dx. (1.6)
For any Borel measure µ on T d and any function f ∈ L 1 (T d , µ), we denote

µ(f ) := ˆTd f (x) µ(dx) (1.7)
which is also extended to vector-valued functions in L 1 (T d , µ) d .

• For two functions f, g defined on a non-empty set E, "f / ≡ g" means that f and g are not identically equal in E.

• C denotes a positive constant which may vary from line to line.

Some recalls of ergodic theory

Invariant measures and rotation sets

Let b a vector field in C 1 (T d ) d . A probability measure µ in M (T d ) is said to be invariant for the flow X associated with the vector field b by (1.1), if one has

∀ t ∈ R, ∀ ψ ∈ C 0 (T d ), ˆT2 ψ Ä X(t,
∀ ϕ ∈ C 1 (T d ), ˆTd b(x) • ∇ϕ(x) µ(dx) = 0. (1.10)
Alternatively, define the set

I b := ¶ µ ∈ M p (T 2 ) : µ invariant for the flow X © , (1.11)
where M p (T d ) is the set of probability measures on T d . Also define the two following non-empty subsets of R d :

• The set of all the limit points of the sequences

Ä X(n, x)/n ä n≥1 for x ∈ T d (denoted by ρ p (b) in [13]) A b := x∈T 2   n≥1 ® X(k, x) k : k ≥ n ´  .
(1.12)

When b does not vanish in T 2 , by virtue of the Peirone result [14, Thm. 3.1] the former definition is reduced to

A b = ® lim t→∞ X(t, x) t : x ∈ T d ´.
(1.13)

• The so-called Herman [START_REF] Herman | Existence et non existence de tores invariants par des difféomorphismes symplectiques[END_REF] rotation set (recall notation (1.7))

C b := ¶ µ(b) : µ ∈ I b © (1.14) which is a compact convex subset of R d .
An implicit consequence of the Misiurewicz-Ziemian results [13, Thm. 2.4, Rem. 2.5] and [START_REF] Misiurewicz | Rotation sets for maps of tori[END_REF]Cor. 2.6] shows that

A b ⊂ C b = conv (A b ).

Herman's rotation set in dimension two

In dimension two [START_REF] Misiurewicz | Rotation sets for maps of tori[END_REF]Thm. 3.4 (c)] shows that actually the set A b (1.12) is a convex set of R 2 , which thus implies that

C b = A b . (1.15)
If in addition b does not vanish in T 2 , we thus deduce from (1.13) that

C b = ¶ µ(b) : µ ∈ I b © = ® lim t→∞ X(t, x) t : x ∈ T 2 ´= A b . (1.16)
This connects closely the invariant probability measures for the flow X (1.1) and the limit points of the flow. Moreover, Herman's rotation set satisfies the celebrated Franks-Misiurewicz ergodic theorem [START_REF] Franks | Rotation sets of toral flows[END_REF] (which actually extends to any continuous flow).

Theorem 1.1 [START_REF] Franks | Rotation sets of toral flows[END_REF]). Let b be a vector field in C 1 (T 2 ) 2 . Then, there exists a vector ξ ∈ R 2 such that the Herman rotation set satisfies C b ⊂ R ξ. More precisely, one of the three following mutually disjoint cases holds:

(i) C b is a unit set.
(ii) C b is a closed line segment of a line passing trough 0 R 2 , and some another point of Q 2 .

(iii) C b is a closed line segment with one end at 0 R 2 , and an irrational slope.

2 A class involving some change of sign

A general two-dimensional result

We have the following result which allows us to exhibit singular two-dimensional ODEs flows.

Theorem 2.1. Let b be a vector field in C 1 (T 2 ) 2 such that the Herman set C b reads as

C b = I k, where I is a closed segment of (0, ∞) and k ∈ Z 2 \ {0 R 2 }. (2.1)
Assume that there exist two functions α ∈ C 0 (T 1 ) and β ∈ C 0 (T 2 ) such that

∀ x ∈ T 2 , b(x) • k ⊥ = α(x • k ⊥ ) + β(x),
with max

x∈T 2 α(x • k ⊥ ) > max T 2 |β|. (2.2)
Then, for any function σ ∈ L 1 (T 2 ) the sections of which satisfy the positivity property for any closed line segment S ⊂ R 2 , ˆS σ > 0 a.e. in T,

where T is any closed line segment orthogonal to S from some end point of S (so that S × T is a closed rectangle of R 2 ), the probability measure σ(x)/σ dx cannot be invariant for the flow (1.1) associated with the vector field b. Moreover, if the vector field b does not vanish in T 2 , there does not exist any real-analytic first integral u for the flow with ∇u ∈ C 1 (T 2 ) 2 . Remark 2.2. In view of the Franks-Misiurewicz theorem, the assumption (2.1) of Theorem 2.1 corresponds in the statement of Theorem 1.1, either to the case (i) when C b is reduced to a commensurable vector in R 2 \ {0 R 2 }, or to the case (ii).

Remark 2.3. It is not a priori evident to provide explicitly a vector field b satisfying the assumptions of Theorem 2.1. Indeed, the same integer vector k has to satisfy both condition (2. Proof of Theorem 2.1. Let b be a vector field in C 1 (T 2 ) 2 satisfying the condition (2.1) with some k ∈ Z 2 \ {0 R 2 }. Assume by contradiction that there exists an invariant probability measure for the flow (1.1), i.e. of the type σ(x) dx where σ is in L 1 (T 2 ) and satisfies the positivity property (2.3). By Liouville's theorem combined with the representation of the divergence free vector fields in R 2 , there exists a current field

ϕ ∈ W 1 loc (R 2 ) with ∇ϕ ∈ L 1 (T 2 ) 2 , such that σ b = ∇ ⊥ ϕ = (∇ϕ) ⊥ in T 2 .
Making the change of variables

x = Ky where K := 1 |k| Ç -k 2 -k 1 k 1 -k 2 å = (K T ) -1 , (2.4
)

which implies that K T k = -|k| e 2 and K T k ⊥ = |k| e 1 .
Then, using (2.2) we get that for a.e. y in the quotient group

R 2 /(K -1 Z 2 ), σ(Ky) b(Ky) • k ⊥ = σ(Ky) α(|k| y 1 ) + σ(Ky) β(Ky) = ∇ ⊥ ϕ(Ky) • k ⊥ = ∇ϕ(Ky) • k. (2.5)
Moreover, by (2.1) we have σ b • k ⊥ = 0. Taking the mean value of (2.5) this yields ∇ϕ(K•) • k = 0. Hence, defining the function ϕ :

Ä x → ϕ(x) -∇ϕ • x ä which belongs to W 1 (T 2 ), we deduce that ∇ϕ • k = ∇ϕ • k -∇ϕ • k = ∇ϕ • k.
Putting this in (2.5) and recalling

K T k = -|k| e 2 , we get that α(|k| y 1 ) + σ(Ky) β(Ky) = ∇ ⊥ ϕ (Ky) • k ⊥ = K ∇ Ä ϕ (Ky) ä • k = -|k| ∇ Ä ϕ (Ky) ä • e 2 = -|k| ∂ y 2 Ä ϕ (Ky) ä .
(

Now, integrating (2.6) with respect to y 2 ∈ R/(k/|k| • Z 2 ) (due to (2.4)) and taking into account the periodicity of ϕ , we obtain that

σ(K•) y 2 (y 1 ) α(|k| y 1 ) + σ(K•) β(K•) y 2 (y 1 ) = 0 for y 1 ∈ R/(k ⊥ /|k| • Z 2 ),
which combined with the positivity property (2.3) of σ, implies that

max x∈T 2 α(x • k ⊥ ) = max y 1 ∈R/(k ⊥ /|k|•Z 2 ) α(|k| y 1 ) ≤ max y 1 ∈R/(k ⊥ /|k|•Z 2 ) σ(K•) β(K•) y 2 (y 1 ) σ(K•) y 2 (y 1 ) ≤ max T 2 |β|.
This contradicts the strict inequality of (2.2) and concludes the first part of the proof. Next, assume that there exists a real-analytic first integral u for the flow, such that the gradient field ∇u is Z 2 -periodic and has a finite number of roots in T 2 . Hence, taking the time derivative in u

Ä X(t, x) ä = u(x), we get that ∀ x ∈ T 2 , ∀ t ∈ [0, ∞), ∇u Ä X(t, x) ä • b Ä X(t, x) ä = 0,
and thus for t = 0, ∇u • b = 0 in T 2 . Now, if the vector field b does not vanish in T 2 , there exists a function σ ∈ C 0 (T 2 ) such that ∇ ⊥ u = σ b in T 2 . So, the continuous function σ does not vanish in R 2 \ {∇u = 0} which by assumption is a connected open set of R 2 . It follows that σ has a constant sign, say positive, in R 2 \ {∇u = 0}. This combined with the continuity of σ yields the positivity property (2.3), which thus contradicts the first part of Theorem 2.1. Therefore, the gradient field ∇u admits at least a cluster point in T 2 , so that u cannot be real-analytic.

The proof of Theorem 2.1 is thus done.

A class of two-dimensional singular flows

We have the following result.

Proposition 2.4. Let b ∈ C 1 (T 2 ) 2 and let α, β ∈ C 1 (T 2 ) be such that the first coordinate of b reads as

b 1 (x) = α(x 1 ) + β(x) for x ∈ T 2 ,
with B := max

T 2 |β| < min Ä -min T 1 α , max T 1 α ä . (2.7)
Then, the two-dimensional flow (1.1) is singular in the sense of Theorem 2.1.

Remark 2.5. It is easy to check that condition (2.7) already prevents Lebesgue's measure from being invariant for the flow associated with the vector field b satisfying (2.7). Indeed, let p 1 , q 1 ∈ T 1 be such that min

T 1 α = α(p 1 ) < 0 and max T 1 α = α(q 1 ) > 0.
Then, we deduce from the inequality of (2.7) that ˆq1

p 1 Å ˆT1 (div b)(x 1 , x 2 ) dx 2 ã dx 1 = Ä α + β x 2 ä (q 1 ) >0 - Ä α + β x 2 ä (p 1 ) <0 > 0.
Therefore, the vector field b of (2.7) cannot be divergence free in T 2 .

Proof of Proposition 2.4. Let us start to prove that for any x 0 = (x 0 1 , x 0 2 ) in T 2 , the first coordinate x 1 := x 1 (•, x 0 ) of the flow X := X(•, x 0 ) (1.1) associated with b, is bounded in [0, ∞). Due to the strict inequality of (2.7) there exists a point p 1 ∈ T 1 such that α(x 0 1 ) = α(p 1 ) < -B. Hence, we have

x 1 -α(x 1 ) = b 1 (X) -α(x 1 ) = β(X) ≤ B < -α(p 1 ) = -α(x 0 1 ) in [0, ∞).
Then, applying the first inequality (2.10) of Lemma 2.6 below to the function x 1 , with f := -α and x 0 := p 1 , we get that

x 1 -x 0 1 = x 1 -x 1 (0) < 1 in [0, ∞).
Similarly, let q 1 be a point in T 1 such that α(x 0 1 ) = α(q 1 ) > B. Then, we have

x 1 -α(x 1 ) = b 1 (X) -α(x 1 ) = β(X) ≥ -B > -α(q 1 ) = -α(x 0 1
) in [0, ∞). Now, applying the second inequality (2.10) of Lemma 2.6 below to the function x 1 , with f := -α and x 0 := q 1 , we get that

x 1 -x 0 1 = x 1 -x 1 (0) > -1 in [0, ∞).
Collecting the two inequalities satisfied by x 1 , we obtain the estimate

x 1 (•, x 0 ) -x 0 1 < 1 in [0, ∞). (2.8)
Next, let us prove that the Herman rotation set satisfies

C b ⊂ R e 2 .
(2.9)

Let µ(b) ∈ C b be associated with an invariant probability measure µ on T 2 . Using successively the invariance property and Fubini's theorem we have for any t > 0,

µ(b) • e 1 = ˆT2 b 1 (x) µ(dx) = 1 t ˆt 0 Å ˆT2 b 1 Ä X(s, x) ä µ(dx) ã ds = ˆT2 Ç 1 t ˆt 0 b 1 Ä X(s, x) ä ds å µ(dx) = ˆT2 x 1 (t, x) -x 1 t µ(dx)
Moreover, taking into account the bound (2.8) Lebesgue's theorem implies that

lim t→∞ ˆT2 x 1 (t, x) -x 1 t µ(dx) = 0,
which proves the inclusion (2.9). Finally, choosing the integer vector k := -e 2 , or equivalently k ⊥ = e 1 , we get that the expression of b 1 in (2.7) agrees with the expression of b • k ⊥ in (2.2). Moreover, it is clear that the inequality satisfied by B in (2.7) implies the inequality in (2.2). Therefore, the vector field b under the assumption (2.7) satisfies the conditions of Theorem 2.1.

We have just proved that the flow associated with any vector field b characterized by the condition (2.7) is singular in the sense of Theorem 2.1.

Lemma 2.6. Let f be a function in W 1,∞ (R/T Z), T > 0 (i.e. f is Lipschitz and periodic of period T in R), let x 0 ∈ R, and let x be a function in C 1 ([0, ∞)) such that f (x(0)) = f (x 0 ). Then, we have the two following implications:

x + f (x) < f (x 0 ) in [0, ∞) ⇒ x -x(0) < T in [0, ∞), x + f (x) > f (x 0 ) in [0, ∞) ⇒ x -x(0) > -T in [0, ∞),
(2.10)

where f (x) := f • x.
The proof of Lemma 2.6 is given in [2, Lem. 2.5].

An extension to dimension three

The two dimensional Proposition 2.4 can be extended partially in dimension three as follows.

Proposition 2.7. Let b ∈ C 1 (T 3 ) 3 and let α, β ∈ C 1 (T 3 ) be such that the first coordinate of b reads as

b 1 (x) = α(x 1 ) + β(x) for x ∈ T 3 ,
with B := max

T 3 |β| < min Ä -min T 1 α , max T 1 α ä . ( 2

.11)

Then, for any non-negative function σ in L 1 (T 3 ) satisfying the weak positivity property

σ x (x 1 ) := ˆT2 σ(x 1 , x ) dx > 0 a.e. x 1 ∈ T 1 , (2.12)
the probability measure σ(x)/σ dx cannot be invariant for the flow (1.1) associated with b.

Proof of Proposition 2.7. First of all, using condition (2.11) and proceeding similarly to the proof Proposition 2.4 with Lemma 2.6, we get that for any x 0 ∈ T 3 , the first coordinate x 1 (•, x 0 ) of the flow (1.1) is bounded in [0, ∞). Hence, it follows immediately that

∀ x ∈ T 3 , lim t→∞ x 1 (t, x) t = 0, (2.13)
which similarly to (2.9) yields

C b ⊂ (R e 2 + R e 3 ). (2.14)
Now, assume by contradiction that there exists an invariant probability measure for the flow (1.1) of the type σ(x) dx where σ is in L 1 (T 3 ) and satisfies the positivity property (2.3). By Liouville's theorem combined with the representation of a divergence free field in dimension free (see, e.g., [START_REF] Girault | Finite element methods for Navier-Stokes equations. Theory and algorithms[END_REF]Thm. 3.4]) there exists a periodic current field ϕ = (ϕ 1 , ϕ 2 , ϕ 3 ) ∈ W 1 (T 3 ) 3 such that

σ b = σ b + ∇ × ϕ in T 3 .
Moreover, by (2.14) we have σ b • e 1 = 0. Then, the representation (2.11) of b 1 gives

σ(x) b 1 (x) = σ(x) α(x 1 ) + σ(x) β(x) = ∂ x 2 ϕ 3 -∂ x 3 ϕ 2 for x ∈ T 3 .
(2.15)

Then, integrating the former equality with respect to x = (x 2 , x 3 ) and taking into account the periodicity of the functions ϕ 2 and ϕ 3 , we obtain that

σ x (x 1 ) α(x 1 ) + σβ x (x 1 ) = 0 a.e. x 1 ∈ T 1 .
This combined with the positivity property (2.12) of σ leads us to max

T 1 |α| ≤ max T 1 σβ x σ x ≤ max T 3
|β|, which contradicts the strict inequality of (2.11), and concludes the proof of Proposition 2.7.

Examples

Gradient flows in any dimension

Consider the flow (1.1) associated with the gradient field b = ∇u ∈ C 1 (T d ) d , where the function u is "affine plus periodic", i.e. the function

Ä x → u(x) -∇u • x ä is in C 2 (T d ) d . Define the T d-1 -mean value function ūx of the function u by ūx (x 1 ) := ˆTd-1 u(x 1 , x ) dx for x 1 ∈ T 1 .
Then, assume that the function ū changes sign, and that the function u is sufficiently close to its mean value ūx in the sense that max

T d ∂ x 1 (u -ūx ) < 1 2 max T d |∂ x 1 u|.
Therefore, due to the equality b 1 = ∂ x 1 u, we easily deduce that the condition (2.11) is fulfilled by the functions α := ∂ x 1 ūx and β := ∂ x 1 (u -ūx ).

The gradient fields of Section 3 below will be used in a completely different manner.

A modified nonlinear pendulum

In [12, §9.2, §9.4, §9.5] the dynamics of the ideal pendulum and of the one with constant forcing are investigated. Here, instead we consider the two-dimensional flow X(•, x) = (ω, θ) solution to the equation

   ∂ t ω(t, x) = ϕ Ä ω(t, x) ä -sin Ä θ(t, x) ä ∂ t θ(t, x) = ψ Ä ω(t, x) ä , with Ä ω(0, x), θ(0, x) ä = x ∈ T 2 , (2.16)
where θ denotes the angular position of the pendulum, ψ(ω) (rather than ω) denotes the angular velocity (which ensures the periodicity condition), and ϕ, ψ are two given functions in C 1 (T 1 ).

If the two following inequalities hold min

T 1 ϕ < -1 and max T 1 ϕ > 1,
then the flow of the modified pendulum (2.16) is singular in the sense of Proposition 2.4.

A two-dimensional Euler flow

We consider the two-dimensional Euler flow X(•, x) for x in the torus R 2 /(2πZ 2 ), solution to the ODEs system

   ∂ t x 1 (t, x) = -A cos Ä x 1 (t, x) ä -B sin Ä x 2 (t, x) ä =: b 1 Ä X(t, x) ä ∂ t x 2 (t, x) = A sin Ä x 1 (t, x) ä + B cos Ä x 2 (t, x) ä =: b 2 Ä X(t, x) ä , (2.17)
where A, B are two fixed non-zero real parameters.

The vector field b defined by (2.17) represents the velocity solution to the steady Euler equation (see, e.g., [START_REF] Wirosoetisno | Nonlinear stability of Euler flows in two-dimensional periodic domains[END_REF])

∂ t ω + ψ 1 ω 2 -ψ 2 ω 1 = 0, (2.18)
where the function ω denotes the (scalar) fluid vorticity, and ψ solution to -∆ψ = ω, denotes the stream function. By definition (2.17) we have

ω(x) = (curl b)(x) = (∂ x 1 b 2 -∂ x 2 b 1 )(x) = A cos x 1 + B cos x 2 and ψ(x) = ω(x).
According to [START_REF] Wirosoetisno | Nonlinear stability of Euler flows in two-dimensional periodic domains[END_REF][START_REF] Constantin | Stratospheric planetary flows from the perspective of the Euler equation on a rotating sphere[END_REF] the study of the stability of the flow (2.17) on the torus is actually relevant for the atmospheric flows.

We have the three following alternatives (see [2, Thm. 2.1] for further details):

• If |A| > |B|, then the functions b 1 , (x → α(x 1 ) := -A cos x 1 ), (x → β(x) := -B sin x 2 )
of the flow (2.17) clearly satisfy condition (2.7). Therefore, the flow (2.17 • If |A| = |B|, then the ODEs system (2.17) is integrable according to [2, Thm. 2.1].

3 A class involving some gradient fields

The main result

We have the following class of singular flows which is based on gradient fields.

the probability measure σ(x)/σ dx cannot be invariant for the flow (1.1) associated b. Moreover, restricting ourselves to dimension two, if the vector field b does not vanish in T 2 , there does not exist any real-analytic first integral v for the flow with ∇v ∈ C 1 (T 2 ) 2 . Remark 3.2. In dimension two the weak positivity on open sets (3.2) is finer than the weak positivity on closed line segments (2.3), and it is also finer than the almost-everywhere positivity on R 2 . The last implication is clear. The first one is deduced by applying Fubini's theorem on any square U of R 2 . Similarly, in dimension three the weak positivity on open sets (3.2) is finer than the weak positivity (2.12). Remark 3.3. In Theorem 3.1 for a given vector field b ∈ C 1 (T d ) d , the non existence of a regular invariant probability measure on T d for the flow (1.1) having a Lebesgue density in L 1 (T d ), does not exclude the existence of a regular invariant measure σ(x) dx with a positive density σ ∈ L 1 loc (R d ), i.e. solution to div (σb) = 0 in R d by virtue of Liouville's theorem. So, by virtue of [6, Thm. 2.15] if b is a non vanishing gradient field in C 1 (T d ) d , such a positive density σ does exist but it is not periodic in general (see [START_REF] Briane | Which electric fields are realizable in conducting materials?[END_REF]Thm. 2.17] for a necessary and sufficient condition). For instance in [6, Exa. 2.12] the two-dimensional gradient field b defined by b := ∇u where u(x • 2π sin(2πx 2 ) e 2 = 4π 2 sin 2 (2πx 2 ) ≥ 0.

) := x 1 -cos (2πx 2 ) for x = (x 1 , x 2 ) ∈ R 2 , ( 3
However, one gets the following explicit regular but non-periodic (with respect to the variable x 1 ) function

σ(x) :=        1 + tan 2 (πx 2 ) e 4π 2 x 1 + e -4π 2 x 1 tan 2 (πx 2 ) if x 2 / ∈ 1 2 + Z e 4π 2 x 1 if x 2 ∈ 1 2 + Z. (3.4
0 = ˆB b(x) • (∇u(x) -∇u) σ(x) dx ≥ c ˆB σ(x) dx > 0,
which leads us to a contradiction. In dimension two, assuming that the vector field b does not vanish in T 2 , the proof of the absence of any first integral with smooth periodic gradient is quite similar to the one of Theorem 2.1. This concludes the proof of Theorem 3.1.

Examples

The first example provides a class of singular gradient flows in any dimension. Note that, since A is symmetric positive definite and inequality (3.8) holds with tr B = 0, the matrix (2A + αB) is also symmetric positive definite. Then, setting X := ∇u and ξ := ∇u, and fixing a symmetric positive definite matrix A and a real number α, by virtue of (3.10) the function F ξ,A,α defined by On the other hand, again using that (2A + αB) is symmetric positive definite, the function F ξ,A,α satisfies min

F ξ,A,α (X) := 1 2 (2A + αB)X • X -(A -α AR ⊥ ) ξ • X for X ∈ R 2 , ( 3
X∈R 2 F ξ,A,α (X) = F ξ,A,α (X ξ,A,α ) = - 1 2 (2A + αB) -1 η • η ≤ 0,
where X ξ,A,α := (2A + αB) -1 η and η := (A -α AR ⊥ ) ξ, (3.13) which implies the equivalences

F ξ,A,α ≥ 0 in R 2 ⇔ η = 0 R 2 ⇔ (A -α AR ⊥ ) ξ = 0 R 2 ⇔ ξ = 0 R 2 . (3.14)
The last equivalence of (3.14) is an immediate consequence of

det (A -α AR ⊥ ) = (1 + α 2 ) det A > 0. (3.15)
Then, we have the following alternative:

• Assume that ∇u = 0 (3.16)

Hence, it follows that if the function u is not a constant and condition (3.8) holds true, the vector field b defined by (3.9) satisfies the non-negativity condition (3.1). Therefore, the flow (1.1) associated with b is singular in the sense of Theorem 3.1.

• Assume that ξ := ∇u = 0 R 2 .

Let K ξ,A,α be the subset of R 2 defined by (recall (3.13))

K ξ,A,α : = ¶ X ∈ R 2 : F ξ,A,α (X) ≤ 0 © = ¶ X ∈ R 2 : (2A + αB) X • X -2 η • X ≤ 0 © .
(3.17

  4.3). Next, we extend the two-dimensional class of singular flows to dimension three through Proposition 2.7, except for the question of integrability which is more delicate in dimension three. Finally, the results of Section 2 are illustrated by various examples in §2.4. Alternatively, in Section 3 we show (see Theorem 3.1) that if the vector field b in C 1 (T d ) d satisfies the non-negativity condition b • (∇u -∇u) ≥ 0 in T d , but not identically null, (1.4) for some periodic gradient field ∇u ∈ C 1 (T d ) d , then the flow (1.1) associated with b has only singular invariant probability measures in any dimension d ≥ 2, and is not integrable in dimension two. This result is illustrated by the class of gradient flows of §3.2.1 in any dimension, and by the class of two-dimensional flows of §3.2.2 which is not restricted to gradient flows.

  2) and the equality set C b = I k. The representation of the Herman rotation set C b (1.14) cannot be simply deduced from the data of the vector field b. However, in §2.2 we provide a class of suitable vector fields b illustrating the theoretical result of Theorem 2.1.

  ) is singular in the sense of Theorem 2.1, and C b ⊂ R e 2 . More precisely, since the vector field b does not vanish in R 2 /(2πZ 2 ), the commensurable case of [14, Thm. 3.1] combined with (1.16) implies that C b = I 2 e 2 where I 2 is a closed interval of (0, ∞).• If |A| < |B|, then permuting the roles of x 1 and x 2 , the flow (1.1) is singular in the sense of Theorem 2.1, and C b = I 1 e 1 where I 1 is a closed interval of (0, ∞).

Theorem 3 . 1 . 0 in T d . ( 3 . 1 )

 31031 Let b be a vector field in C 1 (T d ) d such that there exists a gradient field ∇u in C 0 (T d ) d satisfying the non-negativity condition b • (∇u -∇u) ≥ 0 b • (∇u -∇u) / ≡ Then, for any non-negative function σ ∈ L 1 (T d ) which satisfies the weak positivity property for any non-empty open set U of T d , ˆU σ(x) dx > 0,

. 3 )

 3 satisfies the condition (3.1) with the same gradient field, since b • (∇u -∇u) = Ä e 1 + 2π sin(2πx 2 ) e 2 ä

) Proof of Theorem 3 . 1 .

 31 Let b be a vector field in C 1 (T d ) d satisfying (3.1) with some gradient field ∇u ∈ C 0 (T d ) d . Assume by contradiction that there exists a non-negative function σ ∈ L 1 (T d ) satisfying (3.2) and div (σb) = 0 in T d . Since the function u : (x → u(x) -∇u • x) belongs to C 1 (T d ), the variational formulation (1.10) on T d yields ˆTd b(x) • (∇u(x) -∇u) σ(x) ≥0 dx = ˆTd b(x) • ∇u (x) σ(x) dx = 0, which implies that b • (∇u -∇u) σ = 0 a.e. in T d . Now, since b•(∇u-∇u) is continuous and not identically null in T d , there exists a non-empty open ball B of T d on which it is bounded from below by a constant c > 0. Hence, recalling the weak positivity property (3.2) satisfied by σ, it follows that

. 11 )

 11 provides the representation b • (∇u -∇u) = F ∇u,A,α (∇u) in T 2 . (3.12)

R 2 .

 2 By (3.12) and (3.13) we get the inequalitiesb • ∇u ≥ 0 in T 2 , b • ∇u > 0 if ∇u = 0 R 2 .

)∀ x ∈ T 2 ,

 2 Since the matrix (2A + αB) is symmetric positive definite, K ξ,A,α is a convex compact set of R 2 , which due to (3.14) has a non-empty interior Kξ,A,α . Also note that the boundary ∂K ξ,A,α is the ellipse of equation(2A + αB) (X -λ) • (X -λ) = (2A + αB) λ • λ,centered on the pointλ := (2A + αB) -1 η = (2A + αB) -1 (A -α AR ⊥ ) ξ = 0 R 2 by (3.15),and passing through 0 R 2 . Then, collecting (3.12), (3.13) and (3.17) we get that   b • (∇u -∇u) ≥ 0 if ∇u ∈ R 2 \ K∇u,A,α , b • (∇u -∇u) > 0 if ∇u ∈ R 2 \ K ∇u,A,α = Ø, ∇u(x) ∈ R 2 \ K ∇u,A(x),α(x) ∃ x ∈ T 2 , ∇u(x) ∈ R 2 \ K ∇u,A(x),α(x) ,it follows from (3.18) that the vector field b defined by (3.9) satisfies the non-negativity condition (3.1). Therefore, the flow (1.1) associated with b is singular in the sense of Theorem 3.1.

i.e. a continuous function u satisfying u X(•, x) = u(x) for any x ∈ T d .

Gradient flows in any dimension

Let u be the "affine plus periodic" function defined by (3.5) and assume that

Then, by (3.5) and noting that ∇u = ξ, we get that for any

Therefore, if condition (3.6) is satisfied and the function u is not affine, the gradient field b := ∇u satisfies the non-negativity condition (3.1), so that the flow (1.1) associated with the gradient field b is singular in the sense of Theorem 3.1.

The next example provides a class of singular two-dimensional flows including gradient flows.

A two-dimensional class of singular flows

Let A be a symmetric positive definite matrix-valued function in

) and let B be the matrix-valued function defined by (recall §1.1)

Inequality (3.8) holds in particular for small enough max T 2 |α|, since det A is bounded from below in T 2 by a positive constant.

Let ∇u be a gradient field in C 1 (T 2 ) 2 . Then, we consider the vector field b defined by (3.10)


Claire Lefort
email: claire.lefort@xlim.fr

Emilie Chouzenoux

Jean-Christophe Pesquet

Image reconstruction in multiphoton imaging through multivariate Gaussian fitting

Keywords: Multiphoton Microscopy, Point Spread Function, 3D Gaussian fitting, image restoration, resolution improvement

Multiphoton microscopy (MPM) is a relatively recent tool involved in biological imaging. Its resolution is somewhat limited due to its physical principles, resulting experimentally in a deteriorated planar resolution of about few tenths of micrometers and a reduced axial resolution of a few micrometers. In this communication, we present a numerical approach taking the form of a processing pipeline for image restoration going from the raw images of 0.2 µm diameter fluorescent microbeads to the reconstructed ones. The strategy consists in estimating 3D Point-Spread Function (PSF) shapes in automatically selected volumes of interest, considering a 3D Gaussian profile of intensity. An automatic crop procedure selects areas of a few dozens of individual beads. Our algorithms, based on a variational approach for multivariate Gaussian fitting, allows us to identify evolutions of PSF dimensions along the 3 dimensions of the volume. A proximal alternating optimization method called FIGARO is employed to minimize a cost function serving to quantify the optimality of the process. The difficulty overcome thanks to this novel variational approach for multivariate Gaussian fitting lies in its ability to estimate parametrically Gaussian shapes in a 3D space, which leads to accurate models of the PSF. Thanks to our technique, the PSF measurements show a dimension in (X,Y,Z) of (0.21, 0.27, 1.49) µm. It is the first time to our knowledge that such a direct 3D measurement has been made possible. This process paves the way to a significant improvement of the resolution, perfectly suited to MPM.

INTRODUCTION

The generation of an image requires the use of an instrument in a great majority of situations: a telescope for spatial imaging, a camera for photography, a microscope for microscopy. This intermediary introduces necessarily perturbations inside the resulting image such as blur or noise. Multiphoton microscopy (MPM) is particularly affected by this problem [START_REF] Helmchen | Deep tissue two-photon microscopy[END_REF]. Indeed, this method often involved in biomedical imaging [START_REF] Lefort | A review of biomedical multiphoton microscopy and its laser sources[END_REF] requires the intermediary of (i) a microscope objective, (ii) a tridimensional scanning system, (iii) an infrared laser beam crossing a potentially scattering and absorbing medium. These factors, inherent to the method, deteriorate the image quality and resolution where the image of a lighting point source is a bright spot spread out spatially. This phenomenon is characterized by the Point-Spread Function (PSF) [START_REF] Vo-Dinh | Biomedical Photonics Handbook[END_REF]. In optical microscopy, whatever the solution chosen (confocal, wide-field or multiphoton), the effect of PSF does exist in the 3 dimensions and mainly in the axis of the optical beam, associated with depth deepness imaging and usually designed by axis Z. Consequently, this problematic concerns mainly MPM and less the other solutions which are not adapted to depth imaging.

In the literature, several solutions exist for improving the resolution and reducing the effect of PSF. They are based on sophisticated optical setups [START_REF] Hell | The 2015 super-resolution microscopy roadmap[END_REF], consisting, for example, in shaping the excitation beam. These super-resolution methods allow to reach resolutions at the scale of several tens of nanometers with very nice optical setups, which are also very complex and expensive, resulting in a system far from turnkey commercial solutions.

In this context, the implementation of a numerical method for the computational restoration of multiphoton images and optimization of their resolution appears well-suited: transparent for users, having the ability to be applied in real-time and probably less expensive. In the context of a computational increase of virtual reality, the first step consists in modeling and fitting optimally the profile of the PSF pattern from experimental data. Indeed, the perfect characterization of the PSF of the involved system is essential to perfectly visualize the beam evolution linked to the instrument. This opens the possibility to remove the contribution of the instrument contained into the image in order to keep exclusively the part of the image coming from the samples. In the literature, Gaussian models are considered as good and tractable approximations for PSF modeling in optical microscopy [START_REF] Marim | Improving single particle localization with an empirically calibrated Gaussian kernel[END_REF][START_REF] Zhu | Efficient parallel Levenberg-Marquardt model fitting towards real-time automated parametric imaging microscopy[END_REF][START_REF] Anthony | Image analysis with rapid and accurate two-dimensional Gaussian fitting[END_REF]. However, existing Gaussian fitting techniques remain unsatisfying since (i) they are limited to 1D or 2D fitting exclusively, (ii) they usually assumed zero background (such as in the famous Caruana's approach [START_REF] Caruana | Fast algorithm for the resolution of spectra[END_REF]), (iii) the presence of noise which corrupts the data is quasi never taken into consideration. This limits the applicability of these methods for processing real 3D experimental data, and especially those from MPM. We have recently proposed a mathematical approach aiming at Gaussian fitting multiphoton microbeads. Our mathematical strategy has shown an impressive quality of results and is now in process of publication [START_REF] Tsz-Kit | Optimal Multivariate Gaussian Fitting for PSF Modeling in Two-photon Microscopy[END_REF]. It is based on a variational approach for PSF modeling through multivariate Gaussian fitting allowing to process real 3D datasets highly damaged, in a direct way and without resorting to a combination of 1D or 2D analysis.

In this article we are detailing the steps required for implementing our computational strategy based on a processing pipeline going from the raw images to the reconstructed ones, by estimating 3D PSF shapes in automatically selected volumes of interest and restoring the resulting images. This results in a significant improvement of the resolution. Our strategy performs a trade-off between accuracy and simplicity where a novel variational approach for PSF modeling is proposed through multivariate Gaussian fitting. First of all, the data fitting problem is formulated. Our variational model is proposed and positioned regarding to related works. A proximal alternating optimization method is then introduced, named FIGARO, having a guaranteed convergence in order to find a minimizer of the related cost function. The estimated Gaussian PSF models are then used for restoring a degraded volume. The objective of this publication is to illustrate the validity and the robustness of this multivariate Gaussian fitting approach on real 3D images obtained from a multiphoton microscope. The sample is constituted by fluorescent microbeads placed into a solidified gel and considered as immobile lighting point sources. The resulting images are the input of our processing pipeline for computational restoration. Several relevant images illustrating key steps of the pipeline are presented.

PRESENTATION OF THE EXPERIMENTAL AND COMPUATIONAL STRATEGIES

Experimental setup of multiphoton microscopy

The multiphoton system involved in these experiments is standard and similar to those usually employed in imaging platforms adjacent to labs dealing with life sciences. In our case, a commercial multiphoton microscope (Olympus, BX61WI) has been employed in a routine protocol for two-photon fluorescence imaging allowing the generation of experimental raw dataset of images. The microscope objective is a 25× water immersion (Olympus, XLPLN25×WMP, 1.05 numerical aperture) especially dedicated to multiphoton microscopy. The excitation source is a standard femtosecond titan sapphire laser source, (Chameleon Ultra II, Coherent Inc., 800 nm, 150 fs, 10 nm, 82 MHz, 4 W). Optical densities are preceding the working station input for managing the average power deposited on the sample. An average power of few mW (about 4 mW) is usually enough to generate exploitable raw images in significant depth.

A scanning system, constituted by two galvanometric mirrors, is generating a periodic movement of the excitation beam on all the surface of the microscope objective. It follows a raster scan pattern, resulting in a 2D scanning of a plane of the sample. Each position of the excitation beam on the sample is known. By an analog-to-digital conversion realized by a photomultiplier tube (PMT), the multiphoton process emitted by the samples generates a point on the image whose intensity coded on 4096 grey values is linearly proportional to the intensity of the multiphoton process involved. The detection is synchronized with the scanning system, resulting in a discretization of the detection, presented as a 2D matrix of pixel. Native images presented in this publication often have dimensions of 2048 × 2048 pixels.

The samples imaged are spherical latex microbeads, coated with a fluorophore emitting fluorescence at 515 nm, detected with a PMT coupled with an optical filter 517/45. Specifically dedicated to one-photon fluorescence microscopy (confocal or wide field), these microbeads are absorbing excitation light at 505 nm. In two-photon microscopy, we have chosen the excitation light at 800 nm. A dichroic mirror at 690 nm splits the excitation beam from the laser source and the backfluorescence from the volume of microbeads which is directed to the detection module. The diameter is standardized at 0.2 µm. With such a dimension, the microbeads can be considered as small enough to be assimilated to a spot generating images of the microscope PSF.

Microbeads are spread into gelatin, initially liquid and then solidified in order to have a still bulk volume of microbeads, resulting in a homogeneous distribution of the microbeads along the volume. 3D images of microbeads are systematically recorded, thanks to the inherent ability of multiphoton microscopes to generate optically sectioned images. A planar scanning is generated by galvanometric mirrors and the axial scanning is generated by the motion of the microscope objective with a step mastered with an accuracy of 10 nm.

Theoretical elements of the computational strategy for the restoration of multiphoton images

The retained computational strategyconsists in fitting with Gaussian shapes noisy and blurry bead images represented in a Q-dimensional grid defined by N points. Their coordinates are represented by X=(xn)1≤n≤N . 𝑋 ∈ (ℝ 𝑄 ) 𝑁 and the observed data 𝑦 ∈ ℝ 𝑁 .

∀𝑛 ∈ {1, … , 𝑁}, 𝑦 𝑛 = 𝑎 ̅ + 𝑏 ̅ 𝑝 𝑛 ̅̅̅ + 𝑤 𝑛 (1)
𝑝̅ = (𝑝̅ 𝑛 ) 1≤𝑛≤𝑁 𝑤 = (𝑤 𝑛 ) 1≤𝑛≤𝑁
In this representation of the data, 𝑎 ̅ ∈ ℝ is a background term, 𝑏 ̅ ∈ (0, +∞) a scaling parameter, 𝑤 ∈ ℝ 𝑁 represents the acquisition noise, 𝑝̅ is a noiseless version of the data. Our strategy is based on the estimation of (𝑎 ̅, 𝑏 ̅ , 𝑝̅ ) from the noisy data y assumed to contain a single bead. O strategy stands on the approximation of 𝑝̅ with a multivariate Gaussian shape.

In the situation of our real data from multiphoton microscopy, the resulting images are directly delivered in 3D, resulting in Q=3. At each voxel position xn, each entry 𝑝̅ 𝑛 is assumed to be close to the probability density function named g(u,μ,C) with the inverse covariance matrix 𝐶 ̅ ∈ 𝑆 ++ (𝑄).

With |C| defining the determinant of matrix C. The key point of our numerical strategy results in the fitting of estimates (𝑎 ̂, 𝑏 ̂, 𝑝, µ, ̂𝐶 ̂) of (𝑎 ̅, 𝑏 ̅ , 𝑝̅ , µ ̅, 𝐶 ̅ ), for an optimal description of the data y.

The Gaussian variances in each direction are assumed to be bounded above by some maximal values. We have chosen to introduce some prior knowledge on p and to favor the proximity between p and the Gaussian model parametrized by (μ, C).We propose to measure this closeness by the Kullback-Leibler (KL) divergence [START_REF] Basseville | On entropies, divergences, and mean values[END_REF] from the sought multivariate normal shape to p. The definition of KL divergence in the continuous framework is recalled in the following equation:

(

) 3
where g(x) which represents the function of reference from which p(x) has to be the closest possible. When p(x) is equal to g(x), equation ( 3) is equal to zero.

This approach results in the formulation of a non-convex minimization problem that is solved using a proximal alternating scheme. The resulting algorithm is named "FIGARO". This computational strategy has been implemented by using the software Matlab which is well-suited for managing matrix datasets.

Physical and computational pipeline for data processing applied to images of multiphoton microscopy

On the one hand, our physical strategy for image acquisition from a multiphoton microscope is standard. On the other hand, our computational strategy for Gaussian fitting is totally new. We will illustrate experimentally in the next sections the specific robustness of this method, regarding to the level of degradation of the raw images acquired, a problematic inherent to the multiphoton method of imaging. Indeed, the physical safety of samplesand especially those coming from life sciencesis closely linked to the average powers involved from the Ti: Sa excitation source. Resorting to the weakest average power is of primary interest in imaging applications for life sciences, thus leaving a non-negligible part of perturbations coming mainly from external photons that can be detected experimentally.

Table 1 presents the steps of this computational strategy of restoration of multiphoton images.

Table 1. Pipeline processing of image restoration.

Step

Step The objective is now to illustrate experimentally step by step the result of our computational strategy on raw images obtained from multiphoton microscopy.

ILLUSTRATIONS OF THE PIPELINE FOR THE RESTORATION OF IMAGES OF MULTIPHOTON MICROSCOPY

The mathematical method that we have applied for the computational reconstruction of multiphoton images has been detailed in the previous part. Now, we are going to illustrate each representative step of the computational strategy detailed in Table 1 with images resulting from the application of this strategy.

Illustration of the striking steps of the pipeline

The aim of this part consists in the illustration of the images obtained with the representative steps detailed in Table 1. Indeed, steps 2, 5, 7, 8 and 9 are resulting in visual modification of the images and a gradual increase of their quality, until the restoration of the final image.

These images are presented thanks to the software Image J, an open source especially dedicated to the microscopy community. It is mainly employed for managing the data and for the basic computational processing of images from many kinds of format of images generated by microscope software. Some of the resulting images are presented in 3D where the 3D volume of images is delimited by red lines. The plugin "3D viewer" has been employed for this task. Few images results from a 2D projection either in XY plan stemming from the images directly generated by the microscope or in YZ resulting from a projection of the superimposed stacks.

Step 2: Recording raw images We note here that the resolution in XY is about a half micrometer. In XZ and in YZ, the resolution is about few micrometers. This point is essential in the identification of an important limitation of optical microscopy in general and multiphoton microscopy in particular. Indeed, the interest of optical microscopy lies in its ability to generate images of targets with a minimally invasive process. In other words, the optical beam, in adapted conditions, is able to generate images without a destructive preliminary treatment such as those involved in electron microscopy for example. Optical microscopy is thus an essential solution for obtaining images with a sub micrometric resolution in the XY plan, especially adapted for life sciences. Multiphoton microscopy is particularly interesting thanks to its 2D optical sectioning resulting in the generation of 2D slices at different levels of depth (Figure 1, center and right). But the depth spacing between two consecutive XY stacks is limited by the axial resolution (YZ and XZ) of few micrometers. This has an important consequence on the resulting image, especially in the case of spherical microbeads. Indeed, considering these resolutions linked to the beam spread along the optical axis (Z axis), the image of a sphere appears as an egg-shape, such as a rugby ball. This phenomenon is highly visible on the following images and must be considered and corrected when applied to images of biological samples whose shape in 3D is unknown.

Step 5: search of VOIs Figure 2 illustrates the representation of the beads in the YZ axis from raw extracted data (step 2) and the result of the application of steps 3, 4 and 5 aiming at identifying the quantity of VOIs defined by the user. Here, 40 VOIs have been defined. The resulting VOIs are identified by red rectangles. All the pixels, placed outside the rectangles delimiting the VOIs are equalized to zero.

Please verify that (1) all pages are present, (2) all figures are correct, (3) all fonts and special characters are correct, and (4) all text and figures fit within the red margin lines shown on this review document. Complete formatting information is available at http://SPIE.org/manuscripts Return to the Manage Active Submissions page at http://spie.org/submissions/tasks.aspx and approve or disapprove this submission. Your manuscript will not be published without this approval. Please contact author_help@spie.org with any questions or concerns.

10677 -106 V. A step of binarisation is required to select the VOIs independently from the noise. This step is the cornerstone of the robustness to noise and the reliability of this strategy. The voxels presenting a repartition of intensities close to each other describe the presence of a microbead. The Matlab function "bwconncomp" is involved in this part. This way of process is essential to distinguish the presence of a microbead from the noise having a random repartition which cannot be described by clouds of bright voxels.

The quantity of VOIs can be adjusted by the user. The objective of an adapted adjustment of this quantity lies in the compromise between the thinnest meshing of the PSF in 3 dimensions and the fast delivery of the PSF shape. Indeed, depending on the carrying phase, the PSF can evolve along the depth and a bigger quantity of VOIs is associated with a longer computation time.

Step 7: Crops of raw images

Figure 3 shows few examples, non-exhaustive, of raw data extracted from the selected VOIs. A cloud of bright neighboring pixels spread allowed to recognize the presence of a microbead and to distinguish it from the high level of noise. At this stage, we can clearly see the effect of beam spreading on the image. It results in the 3D representation of an egg-shape and not of a sphere.

Please verify that (1) all pages are present, (2) all figures are correct, (3) all fonts and special characters are correct, and (4) all text and figures fit within the red margin lines shown on this review document. Complete formatting information is available at http://SPIE.org/manuscripts Return to the Manage Active Submissions page at http://spie.org/submissions/tasks.aspx and approve or disapprove this submission. Your manuscript will not be published without this approval. Please contact author_help@spie.org with any questions or concerns.

10677 -106 V. Step

8: Fitting signal with Gaussian model and denoised version of the VOIs

Once the individual volumes are identified and separated from the rest of the total volume, the fitting of the repartition of pixel intensity with a multivariate Gaussian can be performed in reliable and robust conditions. Figure 4 shows the result of denoised version of the three VOIs selected in Figure 3. Starting from the images of Figure 4, the fitting of the pixel intensity by a 3D Gaussian shape has been performed by employing FIGARO algorithm. Our algorithm can deliver a 3D image of the related 3D Gaussian with the multivariate parameters estimated with FIGARO. The resulting 3D images corresponding to the three microbeads chosen for the illustration in Figures 3 and4 are presented in Figure 5.

Please verify that (1) all pages are present, (2) all figures are correct, (3) all fonts and special characters are correct, and (4) all text and figures fit within the red margin lines shown on this review document. Complete formatting information is available at http://SPIE.org/manuscripts Return to the Manage Active Submissions page at http://spie.org/submissions/tasks.aspx and approve or disapprove this submission. Your manuscript will not be published without this approval. Please contact author_help@spie.org with any questions or concerns.

10677 -106 V. Please verify that (1) all pages are present, (2) all figures are correct, (3) all fonts and special characters are correct, and (4) all text and figures fit within the red margin lines shown on this review document. Complete formatting information is available at http://SPIE.org/manuscripts Return to the Manage Active Submissions page at http://spie.org/submissions/tasks.aspx and approve or disapprove this submission. Your manuscript will not be published without this approval. Please contact author_help@spie.org with any questions or concerns.

(

  ∀𝑥 𝑛 ∈ ℝ 𝑄 , ∀µ ∈ ℝ 𝑄 , ∀𝐶 ∈ 𝑆 ++ (𝑄), 𝑔(𝑥 𝑛 , µ, 𝐶) 𝑥 𝑛 -µ) 𝑇 𝐶(𝑥 𝑛 -µ))

8 Fitting 4 Figure 5 FIGARO

 845 Figure 4 Figure 5 FIGARO algorithm with output (𝑎 ̂, 𝑏 ̂, 𝑝, 𝜇, 𝐶 ̂) 9 Reconstruction of images Figure 6 Using OPTIMISM Plugin unsing interpolating PSF

Figure 1

 1 Figure 1 represents the volume of raw images generated by the microscope and presented in the 3D space. Two slices of 2D representation are also shown.

Figure 1 .

 1 Figure 1. Raw images of spherical microbeads obtained with a standard protocol of multiphoton microscopy. From left to right: 3D representation of 230 stacks superimposed in a box with dimensions: 80 × 80 × 23 µm; XY stack 68/230; XY stack 122/230

Figure 2 .

 2 Figure 2. Images of microbeads of 0.2 µm presented in YZ axis with (YZ) dimensions of 80 × 23 µm. Raw images (up) and restored images (down). 40 volumes of interest are selected automatically, including a step of blur reduction.

Figure 3 .

 3 Figure 3. Three examples of 3D VOIs selected resulting in individual crops (Y_crop n°16, 27, 29) from raw data.

Figure 4 .

 4 Figure 4. Selection and crop of the VOIs. Results of the application of steps 6 and 7 on three VOIs selected for the illustration.

Figure 5 .

 5 Figure 5. Images illustrating the 3D Gaussian functions associated to the selected crops presented in the previous Figures3 and 4, thanks to the parameters calculated with FIGARO algorithms.

Figure 6

 6 Figure6summarizes in picture the striking images of steps required for generating a 3D model of the PSF with a Gaussian model on a single microbead. It illustrates the ability of our strategy to define the PSF of the instrument from a cloud of bright pixels thanks to their specific distribution in packets despite the high level of deterioration. This demonstrates experimentally a robustness and reliability yet unequaled in literature, a mandatory step for the restoration of images generated especially with a multiphoton process.

Figure 6 .

 6 Figure 6. Summary of the pipeline for the computational restoration of 3D images of microbeads after the identification of each individual VOI. From left to right: raw data (step 2), raw data with denoising (steps 3, 4 and 5) and reconstruction of microbead with a Gaussian shape (step 8).

  Please verify that[START_REF] Helmchen | Deep tissue two-photon microscopy[END_REF] all pages are present,[START_REF] Lefort | A review of biomedical multiphoton microscopy and its laser sources[END_REF] all figures are correct, (3) all fonts and special characters are correct, and (4) all text and figures fit within the red margin lines shown on this review document. Complete formatting information is available at http://SPIE.org/manuscripts Return to the Manage Active Submissions page at http://spie.org/submissions/tasks.aspx and approve or disapprove this submission. Your manuscript will not be published without this approval. Please contact author_help@spie.org with any questions or concerns.

10677 -106 V. 3 (p.3 of 12) / Color: No / Format: Letter / Date: 4/13/2018 9:38:01 AM SPIE USE: ____ DB Check, ____ Prod Check, Notes:

10677 -106 V. 3 (p.8 of 12) / Color: No / Format: Letter / Date: 4/13/2018 9:38:01 AM SPIE USE: ____ DB Check, ____ Prod Check, Notes:

Starting from these results, a quantitative analysis is possible. The information concerning the positions of VOIs centers in the 3 dimensions, the width of the microbeads image defining the PSF dimension, and the orientation of the bead compared to the optical axis are accessible. Table 2 gathers these values for four microbeads among the 40 that have been selected into the bulk volume. For comparing the relevance of our strategy, we have chosen to compare it with a plugin of Image J, used in routine in imaging platforms involving optical microscopy. MetroloJ is a plugin of Image J especially dedicated to the estimation of PSF relying on a 1D Gaussian fitting along each 3 axis. First, it is not a 3D solution. Then, it requires images of PSF with a good initial contrast and intensity. As illustrated by Table 2, our strategy delivers values of PSF consistent with the expected values of about 0.2 µm. MetroloJ cannot estimate correctly these values due to a too important level of image degradation. The Euler angles (Φ1, Φ2) of microbeads are calculated, characterizing the slope of the main direction of the PSF, which corresponds to the eigenvector of C associated to its largest eigenvalue. Step 9: reconstruction of images

Figure 7 shows the result of the application of the plugin OPTIMISM which aims in the reconstruction of images of microbeads thanks to a step of deconvolution using our PSF model, in a volume where the PSF is assumed to be stationary. Figure 7 illustrates the extrema of the pipeline from raw data to the restored one, involving FIGARO algorithms and OPTIMISM deconvolution. The interest of this strategy is illustrated on images of standardized samples. Indeed, the imaging of known samplesstandardized fluorescent microbeadsallows to check the reliability and robustness of the computational strategy. The use of biological samples could not ensure a reliable result because an unknown shape of the samples cannot guarantee for the identification of consistent results which is the purpose of this presentation.

Exploitation of the results of FIGARO algorithms for a 3D fitting of the PSF

The computational strategy described in Table 1 is reliable and robust to noise. It presents several freedom concerning the quantity of VOIs that can be selected, depending on the level of degradation of the data and depending on the quantity of microbeads that can be identified into the volume. The experimental parameters that result from the fitting of the microbeads by a Gaussian shape thanks to FIGARO algorithms give access to several information concerning the PSF repartition inside the volume: the position of microbeads, their dimensions in the 3 axis and their orientation compared to the axis of the box. This last point was never made possible by processes found in literature, although it gives an important information which is illustrated in Figure 8.

Figure 8 presents a general illustration of the pipeline process that has been applied. A focus on the exploitation of the Gaussian parameters which has been estimated with FIGARO algorithms is now proposed. In Figure 8, the steps of FIGARO algorithms are graphically summarized. Raw data are presented in YZ axis; the repartition of selected and cropped microbeads are then localized at their related position. All the rest of the volume is considered as equal to zero in intensity. A step of denoising inside each of the cropped volumes ensures the reliability of the Gaussian fittings and their robustness to noise. Finally, from the estimation of Gaussian parameters, the orientation of the microbeads into the volume of image is estimated. The 3D representation of Figure 8 gathers the axis of microbeads in a unique center, showing a similar orientation of the microbeads. This last point can reveal a default of orientation of the laser beam or an unplaneity of the target.

CONCLUSION

This paper deals with the problems inherent to the limits of image quality in multiphoton microscopy due to the combined effect of a high sensitivity to noise and deteriorated resolution. This method of biomedical imaging is limited by the presence of blur and noise into the images due to the laser excitation strategy, thus reducing the image contrast and resolution. This results in the degradation of the image quality, lowering the interest of multiphoton microscopy compared to more standard techniques, such as confocal microscopes with their most recent development.

We propose and illustrate here a solution of image restoration based on a computational strategy that we have named "FIGARO". It is composed by 9 steps: 2 for the generation of raw images of multiphoton microscopy and 7 for the numerical analysis and restoration. In order to demonstrate the relevance of our strategy, we have chosen known and standardized samples composed by fluorescence microbeads having sizes standardized at 0.2 µm. Numerically, we have a first stage of preprocessing of the image. It consists in reducing the image blur and noise. Then, a step of localization of clouds of bright pixels is then crucial and determines the robustness of the strategy to external perturbations inside the image. Next, volumes of interest containing a single bead are detected automatically and individually selected and localized. According to the considerations found in literature, we assumed that the intensity profile of the microbead images are Gaussian. The step of Gaussian fitting with our computational strategy of algorithms FIGARO, consists in searching the parameters of a Gaussian having a multivariate shape. The closeness of the parameters of the Gaussian model to the noiseless version of the raw data is evaluated with a Kullback-Leibler divergence. By adding a final step of restoration based on a deconvolution (plugin OPTIMISM from Image J), a complete pipeline processing strategy is demonstrated from raw images up to their restored versions.

The results obtained are presenting a high robustness to noise and blur and show an experimental reliability unequaled in the literature. The bead dimensions calculated are consistent compared to the expected values and the orientation of their axis due to the difference of resolution on the three axis is estimated.

Starting from this results initiating a new computational strategy specifically dedicated to multiphoton microscopy, many other field of improvement of the multiphoton images can benefit from the strategy described by these algorithms. First, switching the excitation source to a one presenting a lower ability of generating multiphoton processes, but interested for other properties such as a wide spectrum, can become more easily conceivable. Then, the application of some of the most recent strategies in terms of image reconstruction and resolution improvement becomes possible. These are a few examples of the future research tracks.