Pierre-Yves Boelle
email: py.boelle@gmail.com

Dominique Debray
email: dominique.debray@aphp.fr

Loic Guillot
email: loic.guillot@inserm.fr

Annick Clement
email: annick.clement@aphp.fr

Harriet Corvol
email: harriet.corvol@aphp.fr

Cystic Fibrosis Liver Disease: Outcomes and risk factors in a large cohort of patients

Keywords: cystic fibrosis, liver disease, cirrhosis, portal hypertension, meconium ileus Cystic Fibrosis Transmembrane conductance Regulator Ultrasonography CT: Computed Tomography MRI: Magnetic Resonance Imaging

Background: Cystic fibrosis related liver disease (CFLD) is a common symptom in patients with cystic fibrosis (CF). However, its prevalence, risk factors, and evolution are unclear. We used a large database of patients with CF to investigate the incidence of CFLD and related risk factors.

Methods:

We retrospectively analyzed 3,328 CF patients with pancreatic insufficiency born after 1985 and recruited into the French CF Modifier Gene Study since 2004. We determined liver status, age at CFLD and severe CFLD onset, sex, CFTR genotype, history of meconium ileus, treatment with ursodeoxycholic acid (UDCA), and respiratory and nutritional status.

Results:

The incidence of CFLD increased by approximately 1% every year, reaching 32.2% at age 25. The incidence of severe CFLD increased only after the age of 5, reaching 10% by age 30. Risk factors for CFLD and severe CFLD were male sex, CFTR F508del homozygosity, and history of meconium ileus. Increasingly precocious initiation of UDCA treatment did not change the in CF

patients. Finally, patients with severe CFLD had worse lung function and nutritional status than other CF patients.

Conclusion:

CFLD not only occurs during childhood as often as believed, but also occurs later in the lifetime of patients with CF. Male sex, CFTR F508del homozygosity, and history of meconium ileus are independent risk factors for CFLD development. Earlier use of UDCA over the last 20 years has not changed the incidence of severe CFLD, a fact that should be considered when prescribing UDCA in the young children, along with its possible adverse effects.

Cystic fibrosis (CF) is the most common, severe, autosomal recessive genetic disease in Caucasians. It is caused by mutations in the gene encoding the cystic fibrosis transmembrane conductance regulator (CFTR), a chloride channel expressed in epithelial cells throughout the body [START_REF] Riordan | Identification of the cystic fibrosis gene: cloning and characterization of complementary DNA[END_REF]. The disease affects several organs such as the lungs, pancreas, intestine, and liver. More than 2000 mutations in the CFTR gene have been described, the most frequent being F508del. The CFTR genotype strongly influences pancreatic function, which is either deficient (PI for pancreatic insufficiency), or normal (PS for pancreatic sufficiency). It is recognized that patients carrying two severe CFTR mutations, also called "CF-causing mutations", have a classical form of CF associated with PI, whereas others have a milder form of disease associated with PS [START_REF] Corvol | Translating the genetics of cystic fibrosis to personalized medicine[END_REF].

Cystic fibrosis-related liver disease (CFLD) includes a wide range of hepatobiliary abnormalities. Focal biliary cirrhosis is the most clinically relevant CFLD, since extension of the initially focal fibrogenic process may lead to multilobular biliary cirrhosis with subsequent portal hypertension and related complications [START_REF] Debray | Best practice guidance for the diagnosis and management of cystic fibrosis-associated liver disease[END_REF]. Indeed, multilobular cirrhosis is recognized to have a significant impact on morbidity, and accounts for ~2.5% of mortality in patients with CF (third cause after respiratory failure and transplantation-related complications) [START_REF] Rowland | Outcome in patients with cystic fibrosis liver disease[END_REF]. The prevalence of CFLD remains controversial, with estimates ranging from 2 to 68% in young patients with CF, due to the lack of a consistent definition [START_REF] Rowland | Outcome in patients with cystic fibrosis liver disease[END_REF][START_REF] Chryssostalis | Liver disease in adult patients with cystic fibrosis: a frequent and independent prognostic factor associated with death or lung transplantation[END_REF]. Identification of patients with CF who are at risk of developing cirrhosis and regular screening for early detection of liver involvement is important as therapeutic intervention is likely to be more effective in patients with early liver disease [START_REF] Kobelska-Dubiel | Liver disease in cystic fibrosis[END_REF]. However, evidence of CFLD is usually subclinical until the disease is advanced, and, therefore, often under-diagnosed. Some risk factors have been recognized as being associated with CFLD, such as severe CFTR genotypes and PI, both of which are dependent upon one another. However, others are still controversial such as male sex and meconium ileus, a severe neonatal intestinal obstruction occurring in ~15% of patients with CF [START_REF] Rowland | Outcome in patients with cystic fibrosis liver disease[END_REF].

In this study, we aimed to investigate the evolution of liver disease and its related risks factors in a large cohort of patients with CF, using longitudinal clinical data collected in the French CF Gene Modifier Study. CFLD was homogeneously defined in the entire cohort, in accordance with the European best practice guidance [START_REF] Debray | Best practice guidance for the diagnosis and management of cystic fibrosis-associated liver disease[END_REF].

Methods

Since 2001, all patients with CF in France are evaluated, at least once a year, by one of the 45 French hospital-based CF centers according to the national CF care recommendations (https://www.hassante.fr/portail/jcms/c_2792719/fr/mucoviscidose) [START_REF] Has | Protocole national de diagnostic et de soins pour une maladie rare[END_REF]. Neonatal CF screening was generalized in France in 2002.

Patients

Patients with CF treated in 38 participating CF centers between January 2004 and January 2017 were enrolled in the French CF Modifier Gene Study. As of January 1, 2017, 4,798 patients with CF had been included (corresponding to ~80% of all French patients with CF ( 8)). Longitudinal data were obtained from electronic medical records or abstracted from the patients' paper records, retrospectively before 2004, and prospectively after January 2004. We analyzed patients born after 1985 with PI (n=3,328) since patients with PS are known to have milder disease (2) (see Supplementary Figure S1 for the flowchart, and Supplementary Information for pancreatic status definition). Patients born before 1985 were excluded to limit selection biases due to an overrepresentation of patients with milder disease in those surviving longer (see Supplementary Information). The following patient characteristics were analyzed: sex, CFTR genotype (described as homozygous for the CFTR F508del mutation, heterozygous for this mutation, or others), and history of meconium ileus (Table 1). Information on liver status, measurements of the forced expiratory volume in one second (FEV 1 ), and body mass index (BMI) were collected at each visit. When no history on CFLD was available in the patient's record, he/she was considered not to be affected. Age at onset of liver disease and severe liver disease were determined (see definition below). Age at ursodeoxycholic acid (UDCA) treatment initiation was also recorded.

CFLD definition

CFLD was defined according to the European best practice guidance by Debray et al. (3) when at least two of the following characteristics were present: 1) abnormal physical examination: hepatomegaly and/or splenomegaly; 2) abnormalities of liver function tests defined as an increase of transaminase (alanine aminotransferase and/or aspartate aminotransferase) and/or gammaglutamyl transpeptidase levels above the upper normal limits (see Supplementary information); 3) ultrasonographic (US) evidence of liver involvement (heterogeneous echogenicity, irregular margins, or nodularity), portal hypertension (splenomegaly, increased thickness of the lesser omentum, spontaneous splenorenal anastomosis, large collateral veins, or ascites), or biliary abnormalities (bile duct dilatation). Patients with cirrhosis, diagnosed by US, computed tomography (CT) and/or magnetic resonance imaging (MRI), and/or portal hypertension (splenomegaly, hypersplenism (platelets <150,000 10 9 /L and white blood cells (WBC) <3000 10 9 /L), and/or spontaneous portosystemic shunts on US), and/or esophageal varices were classified as having "severe CFLD" [START_REF] Debray | Cystic Fibrosis-related Liver Disease: Research Challenges and Future Perspectives[END_REF].

Statistical analyses

Descriptive statistics used mean ± standard deviation (SD), percentages, and hazard ratios (HR) with 95% confidence intervals (CI), as appropriate. The cumulative distribution of CFLD with age was computed using the non-parametric maximum likelihood estimator approach for interval censored data [START_REF] Fay | Rank invariant tests for interval censored data under the grouped continuous model[END_REF]. All patients were considered to be at risk of CFLD since birth. Patients who did not have a CFLD diagnosis before January 2017 were censored at the last visit. In the other patients, the date of diagnosis was computed as the first date when the diagnostic criteria were met based on medical record data. When this date was not precisely known, we considered that it was "interval censored" between birth and the age of the first report (i.e., CFLD onset could occur at any time in this range -this was the case for 142 out of 605 patients with CFLD [24%]). We defined onset of severe CFLD as the first-time cirrhosis and/or portal hypertension and/or esophageal varices were reported, and applied the same approach regarding uncertainty to this date (the date of severe CFLD was interval censored in 19 out of 175 patients [11%]). Comparisons between the cumulative incidences curves were done using the log-rank test adapted to interval censored data [START_REF] Fay | Rank invariant tests for interval censored data under the grouped continuous model[END_REF]. Factors linked to age at CFLD and the time interval between CFLD and severe CFLD were tested using a Coxregression model adapted to interval censored data [START_REF] Pan | Extending the Iterative Convex Minorant Algorithm to the Cox Model for Interval-Censored Data[END_REF]. Confidence intervals were computed by the bootstrap method. We also performed a sensitivity analysis including only patients prospectively enrolled since 2004 (see Supplementary Information).

Age at treatment with UDCA was analyzed as above, with patients not under treatment censored at the age of their last visit. We also tested the impact of treatment on the occurrence of severe CFLD. We used instrumental variable analysis to estimate the true causal association between treatment and disease [START_REF] Kjaersgaard | Instrumental Variable Method for Time-to-Event Data Using a Pseudo-Observation Approach[END_REF]. Indeed, protopathic bias may affect this analysis if UDCA is preferentially prescribed in those who are more likely to develop the disease [START_REF] Bosco | A most stubborn bias: no adjustment method fully resolves confounding by indication in observational studies[END_REF], leading to reverse causality with a higher incidence of CFLD in patients under treatment. A statistical instrument must satisfy: 1) association with UDCA treatment, 2) changing the risk of severe CFLD only through UDCA prescription (i.e. no direct link between instrument and outcome), and 3) not sharing a common cause with treatment (see Supplementary Information). For each analysis, the patients were divided in two groups (those born before and after 1995, in early/late prescribing centers) (see Supplementary Information). The effect of UDCA on severe CFLD incidence was computed as the restricted (at age 20) mean difference in age at severe CFLD onset, using pseudo-observations for survival and the generalized method of moments [START_REF] Kjaersgaard | Instrumental Variable Method for Time-to-Event Data Using a Pseudo-Observation Approach[END_REF].

We computed percent-predicted FEV 1 [START_REF] Quanjer | Multi-ethnic reference values for spirometry for the 3-95-yr age range: the global lung function 2012 equations[END_REF] and BMI Z-scores (WHO2007) with respect to those of healthy populations, and CF-specific percentiles and Z-scores for BMI and FEV 1 using previously published methodology [START_REF] Boelle | Reference percentiles for FEV(1) and BMI in European children and adults with cystic fibrosis[END_REF][START_REF] Kulich | Diseasespecific reference equations for lung function in patients with cystic fibrosis[END_REF]. BMI and FEV 1 Z-scores were averaged from measurements taken in the last 3 years for FEV 1 and the last 2 years for BMI, and compared according to CFLD stage (no CFLD, CFLD, and severe CFLD) using the Kruskal-Wallis test.

Results

CFLD cumulative incidence and risk factors

At the time of the study, 18% of the patients had CFLD and 5% had severe CFLD. The incidence of CFLD increased by approximately 1% every year from birth, reaching 32.2% (CI95%: 29.7, 35.2) by age 25 and leveling out thereafter (Table 2, Figures 1A &S4). The most frequent factors defining CFLD were the joint presence of clinical and biochemical abnormalities (50% of the cases), clinical and ultrasonographic (US) abnormalities (26%) and biochemical and US abnormalities (24%).

The cumulative incidences of individual items defining CFLD (clinical, biochemical, and US abnormalities) are reported in the Supplementary information (Table S2).

The incidence of CFLD was higher in male patients [HR=1.15; 95% CI=0.99-1.36], in CFTR F508del homozygous patients (HR=1.17; 95% CI=1.00-1.37] and in those with a past history of meconium ileus (HR=1.66; 95% CI=1.36-2.01) (Table 3 & Supplementary Figure S5). The CFLD cumulative incidence at 20 years was 32% in patients with meconium ileus vs. 21% in others, 24% in males vs. 21% in females, and 25% in CFTR F508del homozygous vs. 20% in other CFTR genotypes. A multivariable analysis showed that these factors were independently associated with time to CFLD occurrence. Male CFTR F508del homozygous patients with meconium ileus (~4% of the population)

were the most at risk of developing CFLD with a 33% cumulative incidence at age 20, while the least at risk were non-CFTR F508del homozygous females without meconium ileus (~21% of the population) who reached 17% cumulative incidence at the same age (Supplementary Figure S6). The sensitivity analysis of patients prospectively enrolled since 2004 showed similar results (see Supplementary Information, Figure S8 and Table S6).

Severe CFLD and risk factors

Severe CFLD complications were rare up to the age of 5 (Table 2). Overall, the incidence of severe CFLD increased by approximately 0.4% per year after the age of 5, reaching 10.2% (95% CI: 8.5-12.9%) by age 25 (Table 2, Figures 1A &S4). Severe CFLD was diagnosed by cirrhosis (56%), portal hypertension (41%) or esophageal varices (3%). Risk factors for severe CFLD were similar to those for CFLD (Table 3). At age 20, the incidence of severe CFLD was higher in males (8.5% vs. 5.4%), in patients homozygous for the CFTR F508del mutation (8% vs. 6%), and in those affected by meconium ileus (8% vs. 7%) (Supplementary Figure S5).

Treatment with UDCA

Treatment with UDCA increased with age (Figure 1A). At the age of 30, 38% of all patients were under treatment. At the time of the study, 64% of patients with known CFLD were receiving UDCA treatment. Treatment with UDCA had been initiated in the absence of CFLD in 83% of those treated. UDCA prescriptions changed with time and according to CF center. Patients born after 1995 received UDCA earlier than those born before (Figure 1B): at 10 years of age, 7% of patients born between 1986 and 1995 were treated with UDCA compared with 23% of patients born after 1995. CF centers could also be split by early or late UDCA prescribers (Figure 1C). Yet, irrespective of the split, the cumulative incidence of severe CFLD remained the same between groups (Figure 1B &1C). The difference in restricted mean age at severe CFLD onset, according to UDCA treatment, was 2 ± 3.5 years (P=0.56 for comparison to 0) with birth cohort as the instrument, 0 ± 1.0 (P=0.98) with CF center as the instrument, and 0 ± 0.8 (P=0.95) with both (see Supplementary Information for a detailed description of this analysis).

BMI and FEV 1 change according to CFLD and severe CFLD

At the time of the study, BMI and FEV 1 measurements were available for most patients (Table 4). As expected, patients had impaired lung function, with percent predicted FEV 1 (ppFEV 1 ) of 75% of the reference population in the non-CFLD group, 71% in those with CFLD, and 63% in those with severe CFLD (P<0.001). With reference to CF patients, those without CFLD had slightly better lung function (CF-specific FEV 1 percentile 53%) and those with severe CFLD slightly worse (CF-specific FEV 1 percentile 47%) (P=0.04). BMI measurements showed the same trend, with CF patients having an altered nutritional status compared with the normal population, and increasingly worse performance in those affected with CFLD and severe CFLD (P<0.001 and P=0.01 for BMI Z-score and CF-specific BMI percentile respectively). Similar results were obtained in a subset of patients exactly matched for age and sex (see Supplementary Information).

Discussion

The French CF Modifier Gene Study provided an unprecedentedly large database of 3,328 PI patients with CF born after 1985, enabling us to study the incidence of CFLD and severe CFLD. We observed that CFLD increased with age up to 32% by age 30, and increased with independent risk factors such as male sex, CFTR F508del homozygosity, and a history of meconium ileus at birth. We showed that severe CFLD was rare before age 5, increased to 10.2% by age 30 with risk factors similar to CFLD, and was not modified by UDCA treatment. Interestingly, we also found that liver disease was associated with worse lung function and nutritional status.

While CFLD is often believed to develop during childhood, we observed that the incidence rates continued to increase in young adults. CFLD incidence increased between 1 and 2% every year from birth up to the age of 25, reaching 32.2% at this age, at which it a plateaued, which may be because of fewer participants in our study. So far, few studies have reported on CFLD incidence after infancy. In a cohort of 241 young patients with CF followed up for a mean duration of 9.8 years (5 months to 20.6 years), CFLD was mostly diagnosed before 12 years of age but not thereafter [START_REF] Lamireau | Epidemiology of liver disease in cystic fibrosis: a longitudinal study[END_REF]. In a cohort of 177 patients, Colombo et al. observed a higher incidence of CFLD than that in the present study (1.8% per year overall), but a plateau was reached at 20 years for male patients and at 10 years for female patients [START_REF] Colombo | Liver disease in cystic fibrosis: A prospective study on incidence, risk factors, and outcome[END_REF]. These two studies included a limited number of patients, mostly children, whereas we were able to analyze a very large cohort of patients with CF, including children and adults, that is likely more representative of all patients with CF under clinical care. Furthermore, in these two studies, patients were only followed up over a short time, meaning that limited information was available to estimate CFLD incidence over the patients' lifetime. In accordance with our findings, a recent study reported that adult-onset CFLD was likely underestimated [START_REF] Koh | Adult-onset cystic fibrosis liver disease: Diagnosis and characterization of an underappreciated entity[END_REF]. The improved survival of patients with CF in the last decade may also explain the change in the incidence profiles of CFLD [START_REF] Cohen-Cymberknoh | Managing cystic fibrosis: strategies that increase life expectancy and improve quality of life[END_REF].

Severe CFLD complications were rare before 5 years of age, but increased afterwards by 0.3% per year, reaching 10.2% by the age of 30. The prevalence of severe CFLD was higher in male patients, in patients homozygous for the CFTR F508del mutation and with a history of meconium ileus at birth. An international study based on 561 CF patients with severe CFLD (cirrhosis and portal hypertension) showed similar results [START_REF] Stonebraker | Features of Severe Liver Disease With Portal Hypertension in Patients With Cystic Fibrosis[END_REF]. The complications were mostly cirrhosis and portal hypertension, as shown by others [START_REF] Rowland | Outcome in patients with cystic fibrosis liver disease[END_REF][START_REF] Lewindon | Importance of hepatic fibrosis in cystic fibrosis and the predictive value of liver biopsy[END_REF][START_REF] Witters | Noncirrhotic presinusoidal portal hypertension is common in cystic fibrosis-associated liver disease[END_REF][START_REF] Witters | Liver disease in cystic fibrosis presents as non-cirrhotic portal hypertension[END_REF]. As recently described, we also observed a subset of patients with portal hypertension but no cirrhosis [START_REF] Witters | Liver disease in cystic fibrosis presents as non-cirrhotic portal hypertension[END_REF][START_REF] Hillaire | Cystic fibrosis liver disease in adults: Limits of noninvasive tests of fibrosis[END_REF]. Differences in biochemical characteristics in patients with CFLD and severe CFLD were consistent with aggravation of the disease (Supplementary Table S4).

We identified three independent risk factors associated with CFLD and severe CFLD: male sex, CFTR F508del homozygosity, and history of meconium ileus at birth. One third of the patients at most risk (CFTR F508del homozygous male patients with a history of meconium ileus) developed CFLD by the age of 20, while only 16% of the lower-risk factor group (non-CFTR F508del homozygous females without meconium ileus) had developed CFLD at the same age. While several studies reported similar associations between CFLD and severe CFTR genotype [START_REF] Lamireau | Epidemiology of liver disease in cystic fibrosis: a longitudinal study[END_REF][START_REF] Colombo | Liver disease in cystic fibrosis: A prospective study on incidence, risk factors, and outcome[END_REF][START_REF] Efrati | Liver cirrhosis and portal hypertension in cystic fibrosis[END_REF][START_REF] Slieker | Risk factors for the development of cystic fibrosis related liver disease[END_REF], male sex [START_REF] Colombo | Liver disease in cystic fibrosis: A prospective study on incidence, risk factors, and outcome[END_REF][START_REF] Efrati | Liver cirrhosis and portal hypertension in cystic fibrosis[END_REF][START_REF] Slieker | Risk factors for the development of cystic fibrosis related liver disease[END_REF][START_REF] Bartlett | Genetic modifiers of liver disease in cystic fibrosis[END_REF][START_REF] Scott-Jupp | Prevalence of liver disease in cystic fibrosis[END_REF], and meconium ileus [START_REF] Chryssostalis | Liver disease in adult patients with cystic fibrosis: a frequent and independent prognostic factor associated with death or lung transplantation[END_REF][START_REF] Lamireau | Epidemiology of liver disease in cystic fibrosis: a longitudinal study[END_REF][START_REF] Colombo | Liver disease in cystic fibrosis: A prospective study on incidence, risk factors, and outcome[END_REF][START_REF] Efrati | Liver cirrhosis and portal hypertension in cystic fibrosis[END_REF][START_REF] Slieker | Risk factors for the development of cystic fibrosis related liver disease[END_REF][START_REF] Colombo | Analysis of risk factors for the development of liver disease associated with cystic fibrosis[END_REF], some discrepancies exist [START_REF] Feranchak | Cholangiocyte biology and cystic fibrosis liver disease[END_REF][START_REF] Lindblad | Natural history of liver disease in cystic fibrosis[END_REF]. So far, the literature is inconclusive about the impact of CFLD on CF evolution. While some studies have shown that patients with CFLD were likely to have a more severe CF phenotype with altered nutrition and lung function statuses [START_REF] Chryssostalis | Liver disease in adult patients with cystic fibrosis: a frequent and independent prognostic factor associated with death or lung transplantation[END_REF][START_REF] Corbett | Cystic fibrosisassociated liver disease: a population-based study[END_REF][START_REF] Minicucci | Liver disease as risk factor for cystic fibrosis-related diabetes development[END_REF][START_REF] Rowland | Liver disease in cystic fibrosis[END_REF][START_REF] Konstan | Risk factors for rate of decline in forced expiratory volume in one second in children and adolescents with cystic fibrosis[END_REF], others did not observe any association [START_REF] Lamireau | Epidemiology of liver disease in cystic fibrosis: a longitudinal study[END_REF][START_REF] Colombo | Liver disease in cystic fibrosis: A prospective study on incidence, risk factors, and outcome[END_REF][START_REF] Lindblad | Natural history of liver disease in cystic fibrosis[END_REF][START_REF] Gooding | Variceal hemorrhage and cystic fibrosis: outcomes and implications for liver transplantation[END_REF][START_REF] Wilschanski | Clinical and genetic risk factors for cystic fibrosis-related liver disease[END_REF]. We were able to confirm in this large cohort that patients with CFLD had more severe lung disease, with the worst sufferers being patients with severe CFLD complications. Nutritional status showed the same trend. One particularity of this study was the use of CF-specific percentiles and Z-scores for BMI and FEV 1 that allowed referencing of the patients with CF with respect to their peers and, thus, peer to peer comparisons [START_REF] Boelle | Reference percentiles for FEV(1) and BMI in European children and adults with cystic fibrosis[END_REF][START_REF] Kulich | Diseasespecific reference equations for lung function in patients with cystic fibrosis[END_REF]. Finally, we did not observe less severe CFLD as a result of precocious use of UDCA in time or between centers. UDCA efficacy in CFLD remains controversial [START_REF] Cheng | Ursodeoxycholic acid for cystic fibrosis-related liver disease[END_REF]. In sharp contrast with pulmonary and nutritional CF complications, where new treatments have been shown to improve survival and quality of life, no drugs have demonstrated effect on CFLD. Although a recent Cochrane systematic review concluded that "There is currently insufficient evidence to justify [UDCA] routine use in cystic fibrosis" [START_REF] Cheng | Ursodeoxycholic acid for cystic fibrosis-related liver disease[END_REF], several guidelines for CFLD patients' care still recommend its use [START_REF] Has | Protocole national de diagnostic et de soins pour une maladie rare[END_REF][START_REF] Debray | Cystic Fibrosis-related Liver Disease: Research Challenges and Future Perspectives[END_REF]. High doses of UDCA treatment have even been detrimental to patients with primary sclerosing cholangitis, as its biotransformation in the colon produced lithocholic acid, a secondary hydrophobic bile acid with potential toxicity [START_REF] Sinakos | Bile acid changes after high-dose ursodeoxycholic acid treatment in primary sclerosing cholangitis: Relation to disease progression[END_REF][START_REF] Lindor | High-dose ursodeoxycholic acid for the treatment of primary sclerosing cholangitis[END_REF]. However, this was not seen in CF patients in whom the prescribed doses of UDCA are much lower [START_REF] Colombo | Long-Term Ursodeoxycholic Acid Therapy Does Not Alter Lithocholic Acid Levels in Patients with Cystic Fibrosis with Associated Liver Disease[END_REF]. In our cohort, UDCA treatment has been initiated earlier than in CF patients over the last decades, even before the onset of CFLD. In turn, one could expect the occurrence of severe CFLD to either decrease or to be delayed in patients treated with UDCA earlier in life. The same trend would be expected in centers where treatment is given earlier in life.

Analyzing the impact of preventive treatments requires taking into account protopathic bias, whereby treatment is preferentially prescribed to those who are at risk of developing the disease. This could cause CFLD incidences to be larger in those treated, reversing the cause and effect in the conclusion [START_REF] Bosco | A most stubborn bias: no adjustment method fully resolves confounding by indication in observational studies[END_REF]. To avoid this bias, we utilized an instrumental analysis, with birth cohorts and early/late prescribing centers as instruments. We found no changes between early and late treatment with UDCA, based on birth periods or on centers. These results were expected in the case of a limited effect, or in the absence, thereof, of UDCA treatment in the prevention of CFLD complications. However we acknowledge that a true effect of UDCA could be masked for the following reasons: first, if the signs of severe CFLD were recognized increasingly earlier in the most recent cohort, this could, by chance, balance out a delay in severe CFLD symptoms induced by UDCA.

It is, however, unlikely that severe CFLD symptoms would be missed in CF patients born since the 1980s, and this does not explain the absence of a difference between early-and late-treating CF centers. Second, and more subtly, the observed differences in age at UDCA treatment initiation could be irrelevant, if all patients at risk of CFLD were treated early enough in the end. One would only observe the severe CFLD cases of patients who would have progressed to severe CFLD disease irrespective of UDCA treatment, which does not rule out the possibility that other treated patients were actually protected in both cohorts. In that case, we would only have shown the absence of an effect for early versus later UDCA treatment for severe CFLD.

Other limitations in our analysis mainly concern its design, where data for the years before 2004 were obtained by retrospective examination of medical records. On the basis of CF mortality rate, we have estimated that between 2 and 7% of CF patients born between 1986 and 2004 were not included in our study, which could lead to selection by survival (see Supplementary Information). It is also possible that CF clinicians grew more aware of CFLD with time, especially after the study by [START_REF] Colombo | Liver disease in cystic fibrosis: A prospective study on incidence, risk factors, and outcome[END_REF] wherein the diagnostic criteria of CFLD were formalized [START_REF] Colombo | Liver disease in cystic fibrosis: A prospective study on incidence, risk factors, and outcome[END_REF].

Analysis of CFLD incidence in the post-2004 cohort hints that the disease is reported earlier, with a cumulative incidence of 13% [10.5, 15.7] at 10 years of age versus 10% [8.7, 11.0] in the whole cohort (see Figure S8 and Table S6). As mentioned above, this is less likely to be the case for severe CFLD, because the clinical signs are less dependent on interpretation (i.e. cirrhosis, portal hypertension, esophageal varices) and its incidence has not changed over the last 20 years. This is the reason we focused on severe CFLD in the UDCA treatment analysis. Yet, our results will have to be verified by future studies based on prospectively collected data alone.

To conclude, we observed, in this large CF cohort, a high incidence of CFLD. While liver disease is thought to develop in the pediatric age in patients with CF, we were able to show that its incidence continuously increased over time, with a consequential rate of progression to severity.

CFLD and severe CFLD were associated with male sex, CFTR F508del homozygosity, and history of meconium ileus at birth, as well as with worse lung function and nutritional status. Finally, the absence of effect of UDCA treatment on the incidence of severe liver disease is an important finding, as UDCA is commonly prescribed in young patients. In the future, it will be critical that potential therapeutic agents be evaluated in well-designed randomized clinical studies, to ensure that the patients most likely to benefit from the treatment are identified.

Figure 1 .

 1 Figure 1. A) Cumulative incidence of cystic fibrosis related liver disease (CFLD), severe CFLD and

Acknowledgements

The authors would like to thank the clinical research assistants involved in the data collection, namely, Pauline Touche, Tania Kombi, and Julie Mesinele. We would like to especially thank the patients and families who participated in this study and the contributing French CF Centers (see the list of investigators in Supplementary information).

Ethics

The study was approved by the French ethical committee (CPP n°2004/15), and the information collection was approved by the Commission Nationale de L'informatique et des Libertés (n°04.404). Informed consent in writing was obtained from each patient and/or guardian.


Vasilij Gomonaj

Hervé Toulhoat

Selective Oxidation of Methane to Formaldehyde Catalyzed by Phosphates: Kinetic Description by Bond strengths and specific total acidities

Keywords: Heterogeneous catalysis, methane oxidation, formaldehyde, kinetic modeling, specific total acidity descriptor, metal-oxygen bond strength descriptor, DFT, Bond Order Conservation principle

A kinetic model predicts the experimental activity and selectivity patterns obtained in oxidation of methane for divided unsupported orthophosphates (M II 3(PO4)2, M III PO4,)

and pyrophosphates (M II 2P2O7, M IV P2O7) catalysts spanning extended ranges of M-O bond strength (calculated from first principles) and specific total acidity H 0 . A volcano curve is obtained when total methane oxidation rates are plotted against M-O bond strength. M IV P2O7 Pyrophosphates are the most acidic catalysts (H 0 >10 µmol.m -2 ) and remarkably selective to formaldehyde (>90%). For this subset of catalysts, experimental oxidation rates deviate however from predictions in correlation with their acidity, suggesting a significant deactivation phenomenon.

H 0 is shown to be described, and therefore predictable for these phosphates, by a structural parameter, the average P-Oa bond length, Oa being a bridging atom between M and P.

Differences in acidity are interpreted as a consequence of the compression or elongation of the Oa-P bond inside the M-Oa-P bridges.

Introduction

Direct conversion of CH4 with the help of oxidants got much more attention, because it is thermodynamically much easier than under non-oxidative conditions [1].

According to the importance of process and scale of production, production of formaldehyde by direct methane oxidation is on the first places. It is connected with the constant increase in the formaldehyde needs for the processes of organic synthesis, as one important monomer for production of plastics, synthetic fibers and other materials.

In early research of methane oxidation catalysts, metallic platinum and palladium were generally preferred, while silver, gold, nickel and copper were also studied. Among simple oxides, Со3О4, NiO, MnO2 showed the highest activity, Cr2O3, Fe2O3, CuO were moderately active, and ZnO, TiO2, V2O5, MoO3, WO3 were the less active. However, the products of methane oxidation on all oxide catalysts were СО2 and Н2О, while only SiO2 oxide in the form of quartz was discovered selective in methane oxidation into formaldehyde in 1963 [START_REF] Gomonaj | The Role of the Vessel Surface in the Methane Oxidation reaction[END_REF].

Taking into account the high selectivity of quartz at methane oxidation, experiments were carried out to study aluminum silicates as catalysts of selective methane oxidation. For that purpose, Н + ions were exchanged with various cations of valency 2 to 4. Ions exchange from Н + to Na + in aluminosilicate led to full loss of formaldehyde selectivity, indicating that sodium ions are catalytic poisons [START_REF] Gomonaj | Heterogeneously Homogeneous Oxidation of Methane on Aluminosilicates[END_REF][START_REF] Gomonaj | Investigation of the Catalytic Properties of Aluminosilicates in the Methane Oxidation Reaction[END_REF].

Oxides on the basis of cations pairs Al-Si, Al-P, B-P, Al-B, turned out the most selective methane oxidation catalysts in formaldehyde [START_REF] Gomonaj | Self-oscillations in the Heterogeneous Catalytic Reaction of Oxidation of Methane with Oxygen[END_REF][START_REF] Gomonaj | Method for Obtaining Formaldehyde[END_REF][START_REF] Gomonaj | Regularities in the Selection of Catalysts for the Selective Oxidation of Normal Alkanes[END_REF][START_REF] Borko | Relationship between Catalytic Properties and Surface Reducibility of Catalysts in the Partial Oxidation of Methane by Molecular Oxygen[END_REF][START_REF] Gomonaj | Regularities in the Formation of Oxygen-Containing Products and C2-Hydrocarbons from Methane on Acidic and Basic Catalysts[END_REF].

On the majority of catalysts tested in the selective oxidation of methane under conditions of a flow reactor, the yield of formaldehyde (the amount of CH4 (in mol%) converted into CH2O) in a single run of the reaction mixture through the reactor was even lower than 1% [START_REF] Bobrova | Direct Catalytic Oxidation of Methane to Formaldehyde: New Investigation Opportunities Provided by an Improved Flow Circulation Method[END_REF]. This parameter was 2-4% for catalytic systems recognized to be the most promising systems. In the majority of cases, the selectivity of methane conversion into formaldehyde above 50% was observed only at very low (less than 2%) degrees of methane conversion. CH3O " * radicals, were detected by freezing free radicals from a reaction mixture of methane oxidation on SiO2 in a CO2 matrix and a correlation was found between the yield of СН2О and the concentration of HO " * and CH3O " * peroxide radicals, which were detected using EPR spectroscopy, in the reaction vessel volume.

Lately, a lot of publications have appeared dealing with catalysts preparations with different elements and oxides supported on SiO2, on which high selectivity of formaldehyde can be obtained. These works are examined in detail in the thesis work by H. Launay [START_REF] Launay | Développement de Catalyseurs à Base d'Oxyde de Molybdène ou de Vanadium Supporté sur Silice pour l'Oxydation Ménagée du Méthane en Formaldéhyde[END_REF] and related papers ( [START_REF] Nguyen | Study of new catalysts based on vanadium oxide supported on mesoporous silica for the partial oxidation of methane to formaldehyde: Catalytic properties and reaction mechanism[END_REF], [START_REF] Launay | Vanadium species in new catalysts for the selective oxidation of methane to formaldehyde: Specificity and molecular structure dynamics with water[END_REF])

Investigations aimed at elucidating more precisely the nature of the sites responsible for the high rate of generation of formaldehyde on SiO2 and V2O5-SiO2 catalysts were made in [START_REF] Parmaliana | Novel High Activity Catalysts for Partial Oxidation of Methane to Formaldehyde[END_REF].

The formation of formaldehyde as a partial oxidation product over silica was associated with high activities. The results indicated that not all silica samples are equally active, "pyrolytic"

or "fumed" silicas being very much less active than "precipitated" silicas: even the bare silica shows appreciable activity. Vanadium oxide-silica catalysts can perform the partial oxidation of methane to formaldehyde with extremely high activities: the yield of formaldehyde can reach a value of 800 g.kg cat -1 .h -1 ; bare silica also shows appreciable values up to 300 g.kg cat -1 .h -1 .

To obtain information about the mechanism of methane oxidation over surface, experiments were carried out in which were studied the influence of different factors on selectivity [START_REF] Nasser | In situ DRIFT study of nonoxidative methane reaction on Mo/SnO2 catalyst[END_REF]. The nature of intermediate compounds in the catalytic process was detected by recording the DRIFT spectra of surface species formed upon contact of CH4 with Mo/SnO2. The spectra exhibit a discrete band at 1680 cm -1 together with the band at 1190 cm -1 g(H-CO). These bands are typical for n(C=O) stretching in coordinatively bonded formaldehyde on Lewis acid sites. Accordingly, this suggests that the first detectable adsorbed species is formaldehyde.

Consequently, the generation of formaldehyde upon CH4 decomposition may indicate that the Mo/SnO2 has a high concentration of Lewis acid sites. However, the desorption of formaldehyde may be due to its low concentration and/or due to extracting lattice oxygen to form additionally CO and CO2.

In a recent review of "catalytic oxidation of light alkanes (C1-C4) by heteropoly compounds", Sun et al. extensively address many practical as well as fundamental aspects [START_REF] Sun | Catalytic Oxidation of Light Alkanes (C1-C4) by Heteropoly Compounds[END_REF].

They conclude on the difficulty to understand the variety of activities that can be obtained with heteropolyanions catalysts (HPA) in view of the complexity and multiplicity of factors involved, so that "To fully make a deep structure activity relationship in this kind of catalysis, it is necessary to decompose elementary steps and study them in a systematic way".

In that spirit we report here a systematic comparison of divided unsupported phosphates (M II 3(PO4)2, M III PO4,) and pyrophosphates (M II 2P2O7, M IV P2O7) catalytic properties in methane oxidation, showing large variations in oxidation activity and selectivities towards formaldhehyde, CO and CO2. We interpret these patterns through a microkinetic model according to which three descriptors suffice to predict the activity and selectivity of a given phosphate catalyst: the M-O bond strength EMO, where M is the hetero-element, the P-Oa bond strength EPOa, where Oa is a bridging oxygen between M and P, and the experimental specific total acidity H 0 . The EMO/POa descriptors are calculated from first principles [START_REF] Toulhoat | Kinetic interpretation of catalytic activity patterns based on theoretical chemical descriptors[END_REF] on unit-cells of the corresponding bulk phosphate crystal. H 0 turns out to be correlated with the computed average P-Oa distances in M-Oa-P bridges. This approach allows to underline the crucial role of surface acidity in the origin of selectivity to formaldehyde.

Materials and Methods

a) Catalysts preparation and characterization

Catalysts were prepared by precipitating hydroxides by ammonia from solutions of nitrates.

Phosphates were then precipitated from hydroxides by addition of orthophosphoric or pyrophosphoric acid. Precipitated phosphates were collected on a Büchner filter, washed from ions NH4 + and NO3 -, and the filter cake dried at 393K. The dried up powdery catalysts were tableted then calcined at 973K. Tablets were further divided into scraps, which were finally screened to collect the 3-4 mm fraction. These catalysts were characterized by X-ray diffractometry, BET method, and for acidity.

Typical X-Ray diffractograms and BET areas for the various structures obtained are reported in Supporting Information SI-1 as well as the corresponding technical details. BET areas were found to vary in the range 6-90 m 2 .g -1 .

Acidities were determined according to a method derived from that described by Benesi ([18], [START_REF] Benesi | Acidity of catalyst surfaces. II. Amine titration using Hammett indicators[END_REF]), namely titration by n-butylamine of weighed dried powders suspensions in dry benzene, in presence of Hammett indicators in the pKa interval -5.6 to +6.8. We report acidities as 𝐻 & in µmol.m -2 , the integral amount of acid sites measured in the pKa interval, normalized by the BET surface area.

b) Catalytic experiments

Catalytic experiments have been performed using a specifically designed batch reactor, provided with an external recycle pump and a liquid product condenser placed downstream of the reactor and maintained at 10-15°C, which traps the oxygenated products and prevents their further oxidation. This setup is schematised in Supporting Information SI-2. The feed rate of reaction mixture is set between 3 and 15 Nl.hr -1 , while the recycle flow-rate is maintained at 150 Nl.hr -1 so as to ensure the condition of equivalence with a perfectly stirred differential reactor. In these conditions all methane, passed through catalyst, can be completely transformed into products of reaction [START_REF] Gomonaj | Physicochemical and catalytic properties of the phosphate system Co-Ni[END_REF]. The adequacy of this particular setup to compare selective oxidation catalysts has been independently demonstrated [START_REF] Bobrova | Direct Catalytic Oxidation of Methane to Formaldehyde: New Investigation Opportunities Provided by an Improved Flow Circulation Method[END_REF]. Reagents and products of reaction were analyzed by gas chromatography, which used TCD and FID. Analysis of CH4, HCOOH, CO, CO2 and other possible products were carried out with the help of a column filled by Porapak Q. Carbon balance was usually in the range of 1%.

c) DFT calculations

The "Yin-Yang" algorithm has been described in reference [START_REF] Toulhoat | Kinetic interpretation of catalytic activity patterns based on theoretical chemical descriptors[END_REF]. In the variant used here, absolute values of total energies per unit cell are computed for:

1. the bulk compound MiXj, 𝐸 ()*+(,*-, 2. the same unit cell as for the bulk with target atom M removed, 𝐸 ()* , 3. the same unit cell as for the bulk with all atoms except M removed, 𝐸 (,*-.

With 𝑛 the number per unit cell of nearest neighbour atoms X to target atom M, we define the bond energy descriptor 𝐸 /0 as:

𝐸 /0 =
𝐸 /0 can therefore be viewed as a "rebonding" energy of X to its complement in the original bulk structure MiXj.

Calculations were performed through the MedeA interface [START_REF]MedeA: Materials Exploration and Design Analysis software from Materials Design Inc[END_REF]. Structures of the reference phosphates were recovered from the Pearson and ICSD crystallographic databases [START_REF]Pearson's Crystal Structure Database is maintained an published by the American society of Metals[END_REF] through the InfomaticA module in MedeA, which was also helpful in analyzing these structures for the determination of 𝑛 as well as distances and angles, and providing them as inputs of ab initio calculations. Total energy calculations were performed with the VASP version 5.3 software [START_REF] Kresse | Efficient iterative schemes for ab initio total-energy calculations using a plane-wave basis set[END_REF] within the density functional theory, using projected augmented wavefunctions, the PBE functional [START_REF] Perdew | Generalized Gradient Approximation Made Simple[END_REF] in the generalized gradient approximation, Spin Orbit Coupling option, and periodic boundary conditions. The unit cell parameters and atomic positions of the bulk compounds were optimized in order to minimize the total energy 𝐸 ()*+(,*-under the approximations made in VASP, but these geometric parameters were kept unchanged to compute 𝐸 ()* and 𝐸 (,*-. All series of calculations were submitted through the PREDIBOND® module [25] within MedeA, which implements the Yin-Yang algorithm. and 5 pyrophosphates (Mn II 2P2O7, Ti IV P2O7, Si IV P2O7, Ge IV P2O7, Sn IV P2O7). Combining densities and stoichiometries of unit-cells, we report also the number of P atoms per unit volume of bulk material, which is revealed to be almost constant at 1.5x10 -5 P.nm -3 in average, with a mere 13% standard deviation. From this result, we expect superficial densities of P atoms, and therefore of P-Oa bridges, to be of the same order for all the divided unsupported phosphates in our set of catalysts. Assuming that exposed catalytic surfaces have unit-cells scaling like VUC 2/3 where VUC is the unit-cell volume of the bulk material, we estimate this average site density to be of the order of 1.26x10 19 sites.m -2 M-Oa, P-Oa, and P-Ob bond strengths were computed according to the Yin-Yang method as described above, starting from the computed minimum energy unit-cells. The values are reported in Table 2, where structures are listed in the same order as for Table 1. Table SI However, as mentioned above, from a structural argument we expect superficial densities of P-Oa-M bridges to be almost identical for the whole set of divided unsupported phosphates considered here. We expect these superficial bridges to react in presence of traces of water (e.g. product of the total oxidation of CH4), the weakest M-Oa bond being hydrolyzed preferentially according to *P-Oa-M* + H2O = *P-OaH + *M-OH where * denotes a superficial atom.

Results

Therefore, it can be deduced that 𝐻 & actually compares acid strengths for similar acid sites densities, and essentially reduces itself in this context to an intensive quantity.

Computed P-Oa bond lengths are similarly anticorrelated to 𝐻 & , although with less regularity, and might as well be used as a predictive descriptor of 𝐻 & for phosphates.

(refs given) (VO)2P2O7 (refs given) as well as phosphorus-treated Silica, (refs given). In all such

cases this sharp band is surface-sensitive, being perturbed upon adsorption, and has been assigned to the OH stretching of free surface POH groups".

Although we cannot provide IR spectra for the whole set of our phosphate catalysts, we present in Supporting Information the data we have for CrPO4, Ni3(PO4)2, Co3(PO4)2, FePO4, AlPO4 and BPO4: They systematically exhibit similar bands in the vicinity of 3600 cm -1 revealing the existence of surface OH groups on these solids. Details of the derivation are provided in Supporting Information section SI-4.

b-1) Kinetic model of selective methane oxidation

To carry out the catalytic process of methane oxidation, conditions were chosen under which homogeneous stages of formaldehyde decomposition and oxidation are practically excluded (the usage of hardening, filling of free volumes of reactor with matter and other), so in these conditions a heterogeneous mechanism dominates. Table 3 presents the assumed nine steps of such a mechanism, as proposed in [START_REF] Gomonaj | Aluminosilicates -selective catalysts for the oxidation of n-alkanes[END_REF]. We consider these steps a priori reversible.

Table 3: Catalytic elementary steps of methane selective oxidation: Z: free Lewis acidic sites on the catalyst's surface (either M or P cations exposed in surface and coordinatively unsaturated); ZO, ZOHCH3, ZOHCH2, ZOHCH2O .... sites occupied by adsorbed species, respectively -О, -ОНСН3, -ОНСН2, -ОНСН2О ....

Step № Reaction equation

1 O2 (g) + 2Z Û 2ZO 2 CH4 (g) + ZO Û ZOHCH3 3 ZOHCH3 + ZO Û ZOHCH2+ ZOH 4 ZOHCH2 + ZO Û ZOHCH2O + Z 5 ZOHCH2O Û CH2O + ZOH 6 ZOHCH2O+ ZO Û ZHCOOH + ZOH 7 ZHCOOH Û CO (g) + H2O (g) + Z 8 ZHCOOH + ZO Û CO2 (g) + H2O (g) + 2Z 9 2ZOH Û H2O (g) + ZO + Z
Step 1 is the dissociative chemisorption of dioxygen occurring on a superficial cationic site exhibiting a Lewis acidic character. Since all phosphate catalysts may express both M and P cations in their outer surface, our notation describes formally two simultaneous networks of surface reactions contributing in parallel to the gaseous products with different energetic descriptors, but following the same mechanisms. Competing steps 5 and 6 are crucial since the latter is the entrance to oxidative degradation steps 7 and 8. Activation of step 5, i.e.

formaldehyde desorption, relative to step 6, will therefore orient the selectivity towards mild oxidation into formaldehyde instead of total oxidation.

We have reported in Supporting Information the detailed analysis of our kinetic model, leaving in this section only the key results. Our main assumptions involve i) a Langmuir-Hinshelwood mechanism for both types of Z sites, determining the surface coverages of competing species; ii) entrance channel of carbon at the catalyst's surface by a carbidic insertion of methane into a ZO bond yielding HO-Z-CH3 surface species; iii) neglecting coverages of species of order equal or greater than two in coverage by ZO; iv) neglecting the rates of reverse elementary surface reactions before those of the corresponding direct reactions since experiments were conducted in a mode simulating a differential reactor with continuous removal of products and therefore surface intermediates must be far from equilibrium.

For a given catalytic surface denoted by subscript j, we obtain the following equations for the total rate of methane oxidation 𝑤 JK L ,N and selectivities 𝑆 JK P @,N , 𝑆 J@,N and 𝑆 J@ P ,N towards carbonaceous products CH2O, CO and CO2 :

𝑤 JK L ,N = 𝑘 ",N 𝜃 S@KJKT,N (2)
With:

𝜃 S@KJKT,N = ;𝑃 @ P 𝐾 W,N = W " ⁄ 𝑃 JK L 𝐾 ",N Y1 + ;𝑃 @ P 𝐾 W , 𝑗= W " ⁄ ;1 + 𝑃 JK L 𝐾 ",N =] +" (3)
𝑆 JK P @,N = ^_,; ^_,`< ^a,`bcd,`= (4)
𝑆 J@,N = ^e,`^a,`bcd,`

;^_,N<^a ,`bcd,`= ;^e ,`< ^f,`bcd,`=

𝑆 J@ P ,N = ^a,`^f,`bcd,`P ;^_,N<^a ,`bcd,`= ;^e ,`< ^f,`bcd,`=

and:

𝜃 S@,N = ;𝑃 @ P 𝐾 W,N = W " ⁄ Y1 + ;𝑃 @ P 𝐾 W,N = W " ⁄ ;1 + 𝑃 JK L 𝐾 ",N =] +W [START_REF] Gomonaj | Regularities in the Selection of Catalysts for the Selective Oxidation of Normal Alkanes[END_REF] where 𝐾 W,N and 𝐾 ",N denote equilibrium adsorption constants of O2 and CH4 from gas phase on catalyst surface j , and 𝑘 g,N denote forward rate constants with subscript l = 2,5,6,7,8 relating to the step number, and subscript j to the catalyst.

Equilibrium adsorption constants can be expressed as function of adsorption enthalpies and entropies with:

𝐾 ),N = W ? 6
h 𝑒𝑥𝑝 l+mK 6,h <nmo 6,h p qn ⁄ [START_REF] Borko | Relationship between Catalytic Properties and Surface Reducibility of Catalysts in the Partial Oxidation of Methane by Molecular Oxygen[END_REF] and rate constants 𝑘 g,N can be expressed according to Eyring [START_REF] Eyring | The Activated Complex in Chemical Reactions[END_REF]:

𝑘 g,N = ^rs n 𝑒𝑥𝑝 t+mK u,± (n)<wmo u,± (n)x qn y (9)
where 𝑘 z and ℎ are the Boltzmann and Planck constants respectively, and Δ𝐻 g,N ± (𝑇) and Δ𝑆 g,N ± (𝑇) are respectively the free enthalpies and entropies of activation of step l for catalyst j.

As explained in Supporting Information, we postulate linear continuous relationships between enthalpies and bond energies for i =1,2 and l =2,6,7,8:

Δ𝐻 ) & = 𝐸 ) & -𝛽 ) 𝐸 0@ ( 10
)
Δ𝐻 g ± = 𝐸 g ± -𝛼 g 𝐸 0@ (11)
while for forward step 5 we use by analogy the descriptor of catalyst acidity 𝐻 & :

Δ𝐻 € ± = 𝐸 € ± -𝛼 € 𝐻 & (12)
Where 𝛼 € is positive and has the appropriate dimension kJ.m 2 .µmol -2 .

We justify this approach as follows: forward step 5 is the desorption of CH2O from the Z center of a ZOH site, as pictured by Combining equations ( 2) to ( 12), we note that the oxidation rate 𝑤 JK L and selectivities in the three gaseous carbonaceous products are completely determined by parameters

Δ𝑆 ),N & , Δ𝑆 g,N ± , 𝐸 ) &
, 𝐸 g ± , 𝛽 ) , 𝛼 g , for i = 1, 2 and l = 2, 5, 6, 7, 8 as long as catalysts can be described by 𝐸 0@ , for X = P, M, and 𝐻 & . We further assume:

Δ𝑆 ),N & = ∆𝑆 ) & ( 13
)
Δ𝑆 g,N ± = ∆𝑆 ± [START_REF] Parmaliana | Novel High Activity Catalysts for Partial Oxidation of Methane to Formaldehyde[END_REF] and in order to express the expected increasing affinity of O2 for Z sites and decreased affinity of CH4 for ZO sites as 𝐸 0@ increases, we set 𝛽 W = -1 and 𝛽 " = 1.

Therefore 15 parameters remain to be determined through a fitting procedure:

∆𝑆 W & , ∆𝑆 " & , 𝐸 W & , 𝐸 "
& , ∆𝑆 ± , 𝐸 g ± and 𝛼 g for l = 2, 5, 6, 7, 8.

Since we consider X=M and X=P separately, we sum up predicted rates for both kinds of sites and add moreover a catalyst independent thermal rate:

𝑤 JK L = 𝑊 JK L / + 𝑊 JK L ? +𝑊 JK L ns = 𝑘 " / 𝜃 /@KJKT + 𝑘 " ? 𝜃 ?@KJKT +𝑊 JK L ns [START_REF] Nasser | In situ DRIFT study of nonoxidative methane reaction on Mo/SnO2 catalyst[END_REF] where the superscripts M and P refer to the choice of descriptors 𝐸 /@ or 𝐸 ?@ in equations SI- [START_REF] Bobrova | Direct Catalytic Oxidation of Methane to Formaldehyde: New Investigation Opportunities Provided by an Improved Flow Circulation Method[END_REF] and [START_REF] Launay | Développement de Catalyseurs à Base d'Oxyde de Molybdène ou de Vanadium Supporté sur Silice pour l'Oxydation Ménagée du Méthane en Formaldéhyde[END_REF]. The thermal rate, a 16 th free parameter, sets a lower limit to the rate in absence of catalyst, and imparts more physical realism by allowing a smooth convergence of the predicted residual rate for vanishing or very high values of the descriptors 𝐸 /@ or 𝐸 ?@ .

Accordingly, predicted selectivities become:

𝑆 ƒ = " … † L ‡ ˆ… † L 𝑆 ƒ / + " … † L ‰ ˆ… † L 𝑆 ƒ ? (16)
For 𝐴 = 𝐶𝐻 " 𝑂, 𝐶𝑂, 𝐶𝑂 "

We describe our two steps fitting procedure in Supporting Information, and report the best fitting 15 parameters in Table SI-4-1.

Resulting from this optimization, best fitted forward rates constants 𝑘 g,N / for l = 2, 5, 6, 7, and the twelve phosphate catalysts considered are plotted against 𝐸 /@ in Figure 2: while forward step 2 appears limiting, forward steps 5, 6, 7 and 8 exhibit higher rate constants by to 13 orders of magnitude. This is consistent with one of the key assumptions underlying our derivation of the kinetic model.

𝑘 €,N is according to our model independent of 𝐸 /@ but linearly dependent of 𝐻 & , overcomes 𝑘 •,N for the 4 most acidic catalysts, while both constants are very close in the case of BPO4.

Plotting rates versus 𝐸 0@ , we clearly obtain a volcano curve, as displayed in Figure 3. Except for the four pyrophosphates, for which predicted rates appear well below experimental rates.

By contrast, predicted selectivities match experiments very well for the complete set of 12 divided unsupported phosphate catalysts, as shown in Figures 4 and5.

Coming back to the outlying data in Figure 3 corresponding to the four pyrophosphates M IV P2O7, the nice linear relationship obtained when plotting in a log scale for these catalysts the ratios of predicted over measured oxidation rates versus total acidity function 𝐻 & (see Figure 6), suggests that for this set of highly acidic catalysts a significant deactivation has occurred along the catalytic test, for instance due to the deposition of carbonaceous species as expected when acidic catalysts are exposed to hydrocarbons at high temperatures. This important point must be further investigated, since if such a deactivation could be avoided or limited, it might provide a way to develop CH2O selective catalysts with greatly improved productivities. For instance, we find here a productivity in CH2O of 10.3 g.kg -1 .h -1 for the catalyst SnP2O7 which would amount to 262 g.kg -1 .h -1 if the theoretical activity was recovered, not so far from the state of the art. (𝐸 0@ computed at the DFT level). Deviations of the experimental rates from predicted rates for strongly acidic pyrophosphates MP2O7 (brown symbols) are the larger the higher the acidity. However, we expect a non-linear relationship to hold between bond energies and bond lengths for a homologous family of compounds like phosphates: such a potential energy function can be represented by the well-known Morse potential:

𝐸(𝑟) = 𝐷 • l31 -𝑒 '('+" h ) > " -1p (17)
On Figure 7, we have plotted EPOa (Table 2) against d1aC (Table 1) together with a Morse potential best fitting the data. We compare DFT computed values of bond energies and bond lengths in order to ensure consistency. We evaluate the relative error on EPOa to 5%, and to 0.5 % on d1aC. The 2xstandard deviation from the best fitting Morse potential is 4 %. This potential has a minimum for d1aC= 1,53 Å, very close to the value of 1,535 Å reported for the P-Oa bond distance in H3PO4, which can be considered as representative for the unconstrained phosphate ion.

In this way, we understand that in pyrophosphates M IV P2O7 the P-Oa bonds are compressed, while they are elongated in Mn2P2O7 and orthophosphates, relative to the unconstrained PO4 3- ion. Constraints acting on P-Oa bonds are the resultant of several factors brought in by heteroions M, like their oxidation state, ionic radii, electronic configuration, which contribute to determine the free energy minimum corresponding to the actual crystal structures. In a compressed bond, where electrons are more confined, the bonding electronic density at the critical mid-point must increase, therefore the bond order, so as to compensate for the increase repulsion of the closer nuclei. It is the reverse in an elongated bond, relative to the unconstrained same bond. Finally, in order to understand the connection shown on Figure 1 between 𝐻 & and the P-Oa bond length, we invoke the Bond Order Conservation (BOC) principle [START_REF] Shustorovich | Chemisorption phenomena: Analytic modeling based on perturbation theory and bond-order conservation[END_REF], [START_REF] Van Santen | On Shustorovich's Bond-Order conservation method as applied to chemisorption[END_REF]: indeed, Brönsted acid hydroxyls will be necessarily expressed on the phosphate catalysts surfaces at Oa* surface atoms inherited from bulk bridging Oa atoms, and will result as discussed above from the dissociative chemisorption of water, implying hydrolysis of the weakest bonds, which are always the M-Oa bonds. From that, we expect that the more compressed, or higher order the Oa-P bond in the bulk, the weakest and lowest order the corresponding H-Oa* bond at the surface. Bond orders of Oa-P bonds estimated according to equations ( 1) and ( 2) given by van Santen in [START_REF] Van Santen | On Shustorovich's Bond-Order conservation method as applied to chemisorption[END_REF] are reported in Table 4. The total acidity 𝐻 & and bond orders A or B are indeed fairly correlated, as shown on figure 8.

Ramis, Busca et al. [START_REF] Ramis | Surface Acidity of the Layered Pyrophosphates of Quadrivalent Ti, Zr, Ge, and Sn and Their Activity in Some Acid-catalysed Reactions[END_REF] also raised the important question of the origin of variations in acid strength of surface P-OH groups in metal phosphates, and tentatively addressed this question in terms of nature of the chemical bonds involved, as it should make sense: "It then seems reasonable to propose that the Brönsted-acid strength of the surface POH groups increases with increasing covalency of the (PO)-M bond, which may cause an increase in the ionic character of the POH bonds. However, increased covalency of the metal-phosphate bonds decreases the cationic character of the bonded metal atoms, so decreasing their Lewis-acid strength when they are exposed on the external crystal surface." Their statement implicitly relates surface and bulk bonding properties. The latter are reflected in the bulk structures. We had independently tested a possible correlation between the Pauling electronegativity of elements M and the measured (surface) acidity of our metal phosphates but found it rather unconvincing (R 2 = 0,17). Our approach based on a DFT evaluation of bond strengths and from that on bond orders provides an alternative interpretation.

A confirmation of this ansatz will however necessitate a further study, out of scope of the present report, involving for instance DFT simulations of realistic models of phosphate catalysts surfaces.

Conclusions

We have reported the catalytic properties towards oxidation of methane at 873 K of an ensemble of twelve unsupported phosphate catalysts of general formulae M II 3(PO4)2, M III PO4, M II 2P2O7, and M IV P2O7, exhibiting a wide range of acidity.

Oxidation products collected from a simulated continuous stirred tank reactor based on a high recycle flow differential reactor, were formaldehyde CH2O, H2O, and the carbon oxides CO and CO2.

Selectivity towards CH2O is strongly correlated to the specific total acidity 𝐻 & , measured using Hammett indicators according to the Benesi method. The activity and selectivity patterns are shown to be accurately predicted thanks to a microkinetic model involving three descriptors to represent the catalysts: 𝐻 & , and the M-Oa and P-Oa bond strengths 𝐸 /@, and 𝐸 ?@, , as computed by DFT according to the "Yin-Yang" method, where Oa is a bridging O atom between ions M and P in the perfect crystal unit-cell.

The activity of catalysts, except that of the four highly acidic M IV P2O7 pyrophosphates, follow a volcano curve when plotted against the descriptors 𝐸 /@, or 𝐸 ?@, . Outlying pyrophosphates experimental activities however deviate from the model prediction by amounts strictly correlated to the corresponding 𝐻 & values, suggesting they were affected by a deactivation phenomenon linked to their high acidities, e.g. irreversible coverage of surface O vacancies (Z sites) by unsaturated hydrocarbons products of acid-catalyzed side reactions.

Acidities 𝐻 & are shown to be closely anticorrelated to P-Oa bond lengths, as measured in the bulk unit-cells of the corresponding cristallized phosphates. We interpret this finding in terms of Brönsted acid sites exclusively associated to surface P-Oa bonds, the surface density of the latter reflecting the nearly constant volumic density of P atoms across the set of phosphates considered. The extensive quantity 𝐻 & can hence coincide with an intensive structural quantity.

Being able to fit the plot of computed P-Oa bond energies versus P-Oa bond lengths with a Morse potential, we analyze the observed variation of P-Oa bond lengths in the set of phosphates studied in terms of compression or elongations relative to that in the free orthophosphate ion, the latter being characterized by a bond length close to the Morse potential minimum.

Further, we propose that the observed anticorrelation (Figure 1) is a manifestation of the Bond Order Conservation principle, since protons are expected to interact all the less with surface Oa atoms, as the latter are engaged in a P-Oa bond of higher order, more confining to their valence electrons.

We expect our findings will be useful in the quest of new catalysts allowing high productivities in valuable mild oxidation products of methane, providing computational and experimental descriptors in the broad family of phosphates.

Analogous concepts might apply for other families of catalysts like sulfates, nitrates, borates, aluminates and silicates. Importantly, all used catalysts were again characterized, showing no change with respect to the freshly calcined solids.

Table of Contents Graphic

We present hereafter for each phosphate catalyst representative of a particular structure (i.e.

space group): 1) experimental X Ray diffraction patterns; 2) simulated diffractions patterns for the matching reference structure; 3) a list of experimental and calculated reticular distances and the corresponding relative intensities for the main diffraction peaks; 4) the corresponding correlation diagrams between distances, with the regression line (dotted line, with equation and squared coefficient of correlation in inset).

Mn P 2 O 7 C 2/m

Mn Step 9 represents the associative desorption of water from two vicinal ZOH sites.

In what follows, we will distinguish steps 1, 2 and 9 involving reversible exchanges with the gas phase from steps 3, 4, 6, which describe the successive oxidation steps of a chemisorbed methane molecule at the surface, and from steps 5, 7 and 8 the desorption to gas phase of the final products formaldehyde, CO, CO 2 and H 2 O. Since the experiments involve a differential reactor with continuous removal of products from the gas phase recycled at high flow-rate at fixed composition, we will consider steps 3, 4, 6, 5, 7 and 8 as very far from equilibrium, i.e.

negligible rates of the reverse reactions.

Steps 1, 2 and 9 contribute to determine surface coverages by O, OH and CH 4 (and its products of oxidation in adsorbed state) as ZO, ZOH and ZOHCH3 sites competing for free Z sites according to a Langmuir-Hinshelwood scheme. The chemisorption of CH 4 on ZO can be considered as well as following the Eley-Rideal mechanism, however it also possible to consider it formally as the chemisorption on a Z site of methanol produced from CH 4 in the gas phase.

With [Z°] the total concentration of surface sites, the balance of surface sites in equilibrium with the gas phase may be written: eq. SI-4_5

𝑍 > =
Eliminating further [ZO]:

1 = H H U 1 + 𝑃 : J 𝐾 L L A 1 + 𝑃 𝐶𝐻 4 𝐾 2 + 𝑃 : J 𝐾 L 𝑃 𝐻 2 𝑂 S T 1 2
Assuming 𝐾 L ~𝐾W and since the recycle gas flow is dried 𝑃 : J ≫ 𝑃 M J : we end up with tha approximation:

1 ≈ 𝑤 N: J = 𝑘 j 𝜃 H: 𝜃 HMN::

M = ] J ] d ] f ^_R a ^_RQ`Qa ] b c] d ^_R ] e c] f ^_R eq. SI-4_16
While the total rate of methane oxidation is:

𝑤 NM P = 𝑘 A 𝜃 H:MNMO eq. SI-4_17
The selectivity towards formation of a given product is determined by the ratio of its formation to the overall rate of reaction. Dividing equations SI-4_14 to SI-4_16 by equation SI-4_17 we obtain the adequate fractional selectivities (summing up to 1):

𝑆 NM J : = ] b ] b c] d ^_R eq. SI-4_18 𝑆 N: = ] e ] d ^_R ] b c] d ^_R ] e c] f ^_R eq. SI-4_19 𝑆 N: J = ] d ] f ^_R J ] b c] d ^_R ] e c] f ^_R
𝑘 z,l = ] { | } 𝑒𝑥𝑝 tuv ~,w ± x yx eq. SI-4_22
Where 𝑘 € and ℎ are the Boltzmann and Planck constants respectively, and Δ𝐺 z,l ± 𝑇 is the free energy of activation of step l for catalyst j.

In order to test the ability of this kinetic model to predict the patterns of activity and selectivities obtained from our methane oxidation experiments catalyzed by the set of 12 solid phosphates described in the present report, the next step was to introduce quantitative dependencies of free energies Δ𝐺 &,l > 𝑇 and Δ𝐺 z,l ± 𝑇 with respect to numerical descriptors of the catalysts. For the latter, the natural choice was bond energies 𝐸 7:,l and 𝐸 9:;,l computed by DFT, and experimental total acidities 𝐻 l > as presented and tabulated in the main text.

Decomposing free energies of adsorption and activation into their entropic Δ𝑆 &,l > Δ𝑆 z,l ± , and enthalpic Δ𝐻 &,l > , Δ𝐻 z,l ± components assumed temperature independent, we postulate linear continuous relationships between enthalpies and bond energies:

Δ𝐻 & > = 𝐸 & > -𝛽 & 𝐸 …: eq. SI-4_23
Δ𝐻 z ± = 𝐸 z ± -𝛼 z 𝐸 …: eq. SI-4_24

Where 𝐸 …: is the bond energy, and in our case X=M or X=P, since we consider M and P cations at the surface as contributing potential sites. Energies 𝐸 & > , 𝐸 z ± and adimensional coefficients 𝛽 & , 𝛼 z are parameters to be determined in the fitting protocol described below.

Equations SI-4_23 are usually referred to as compensation relationships, and equations SI-4_24

as Brönsted-Evans-Polanyi (BEP) relationships.

The selectivity in formaldehyde is, according to this model, crucially determined by the competition between forward steps 5 and 6. Δ𝐻 h ± = 𝐸 h ± -𝛼 h 𝐻 > eq. SI-4_25

Where 𝛼 h is positive and has the appropriate dimension kJ.m 2 .µmol -2 .

Combining equations S-4_7 to S-4_25, we note that the oxidation rate 𝑤 NM P and selectivities in the three gaseous carbonaceous products are completely determined by parameters Δ𝑆 &,l > , Δ𝑆 z,l ± , 𝐸 & > , 𝐸 z ± , 𝛽 & , 𝛼 z , for i = 1, 2 and l = 2, 5, 6, 7, 8 as long as catalysts can be described by 𝐸 …: , for X = P, M, and 𝐻 > . We further assume:

Δ𝑆 &,l > = ∆𝑆 & > eq. SI-4_26

Δ𝑆 z,l ± = ∆𝑆 ± eq. SI-4_27

and in order to express the expected increasing affinity of O 2 for Z sites and decreased affinity of CH 4 for ZO sites as 𝐸 …: increases, we set 𝛽 L = -1 and 𝛽 A = 1.

Therefore 15 parameters remain to be determined through a fitting procedure: ∆𝑆 L > , ∆𝑆 A > , 𝐸 L > , 𝐸 A > , ∆𝑆 ± , 𝐸 z ± and 𝛼 z for l = 2, 5,

Since we consider X=M and X=P separately, we sum up predicted rates for both kinds of sites and add moreover a catalyst independent thermal rate: eq. SI-4_29

𝑤
For 𝐴 = 𝐶𝐻 A 𝑂, 𝐶𝑂, 𝐶𝑂 A

We have adopted a two steps fitting procedure. In the first step an objective function is defined as the root mean square deviation (rmsd) of predicted from experimental Neperian logarithms

and discussion a) Characterization of catalysts a- 1 )

 1 Experimental structural properties Twelve catalysts are compared in the present study, including 7 divided unsupported orthophosphates (Ni II 3(PO4)2, Co II 3(PO4)2, Mg II 3(PO4)2, Fe III PO4, Cr III PO4, Al III PO4, B III PO4)

Figure 1 :- 2 b)

 12 Figure 1: Correlation between the specific total acidity of phosphate catalysts and the average

  Fig. SI-5-3, therefore breaking a Z-C bond. CH2O is a very polar molecule, and in the adsorbed state its carbonyl group will undergo an electrophilic attack of the terminal oxygen by any vicinal acidic proton (e.g. H d+ ….O d-=P) and a nucleophilic attack of the carbon by the conjugated base (e.g. O -=P) ,thus weakening the Z-CH2O bond. Therefore, the more acidic the catalyst, the easier the desorption of formaldehyde and this explains why the total acidity 𝐻 & can be taken as a proximal descriptor of Z-CH2O bonds weaknesses.

Figure 2 : 1 Figure 3 :

 213 Figure 2: Plot against 𝐸 /@ of the best fitted forward rates constants 𝑘 g,N / for l = 2, 5, 6, 7, 8

Figure 4 :Figure 5 :Figure 6 :

 456 Figure 4: Comparison of predicted and experimental selectivities for the complete set of 12

Figure 7 :

 7 Figure 7: Plot of computed P-Oa bonds strengths (Table2) versus the corresponding computed

Figure 8 :Mn 2

 82 Figure 8: Plot showing the correlation obtained between experimental total acidity 𝐻 & and

For catalysts Mn 2 P 2 O 7 ,

 227 CrPO 4 , Co 3 (PO 4 ) 2 , Ni 3 (PO 4 ) 2 , FePO 4 , AlPO 4 , BPO 4 , SiP 2 O 7 , Powder diffraction spectra were recorded in the 2q = 20-50° Bragg angle interval thanks to an automated DRON-UM1 diffractometer equipped with a Cobalt anticathode and an iron filter. Inter-reticular distances d i were computed according to the Wulff-Bragg formula 𝑛𝜆 = 2𝑑 & sin 𝜃 & with 𝜆𝑚 = 1,7909 Å the Co Ka incident radiation wavelength. For catalysts TiP 2 O 7 , GeP 2 O 7 , SnP 2 O 7 , and Mg 3 (PO 4 ) 2 powder diffraction spectra were recorded in the 2q = 10-60° Bragg angle interval using a more recent diffractometer equipped with a Copper anticathode ( 𝜆𝑚 = 1,5418 Å ). The diffraction peaks intensities were normalized between 0 and 100 arbitrary units.The powder diffraction spectra were recorded on all phosphate samples firstly after drying the freshly prepared precipitates at 393 K, and secondly after calcination at 973K. Dried samples exhibited very diffuse diffraction lines, indicative of poorly crystallized or amorphous solids, while calcined phosphates exhibited thin lines with well-defined peak positions, indicative of well crystallized solids. The corresponding crystal structures were determined by comparison with tabulated reference spectra and moreover with simulated spectra as presented below.Simulations were performed thanks to the VESTA program (1) on the basis of the published reference structures as cited in the main text.

  ] d ^_R ] e c] f ^_R eq. SI-4_13 And therefore the rates of production of the final gaseous products of methane oxidation CH 2 O, CO and CO 2 : 𝑤 NM J : = 𝑘 h 𝜃 H:MNMA: = ] b ] J ^_R ^_RQ`Qa ] b c] d ^_R eq. SI-4_14 𝑤 N: = 𝑘 i 𝜃 HMN::M = ] e ] J ] d ^_R J ^_RQ`Qa ] b c] d ^_R ] e c] f ^_R eq. SI-4_15

  Forward step 5 is the desorption of CH 2 O from the Z center of a ZOH site, as pictured below by Fig. SI-5-3, therefore breaking a Z-C bond. CH 2 O is a very polar molecule, and in the adsorbed state its carbonyl group will undergo an electrophilic attack of the terminal oxygen by any vicinal acidic proton (e.g. H d+ ….O d-=P) and a nucleophilic attack of the carbon by the conjugated base (e.g. O -=P) , thus weakening the Z-CH 2 O bond. Therefore, the more acidic the catalyst, the easier the desorption of formaldehyde and this explains why the total acidity 𝐻 > can be taken as a proximal descriptor of Z-CH 2 O bonds weaknesses. By analogy with BEP relationships like equations SI-4_24, we introduce for forward step 5:

Figure SI- 5 - 3 :Figure SI- 5 - 4 :

 5354 Figure SI-5-3: atomistic model of ZOHCH2O surface species produced from step 6 in the case of a tri-periodic 2 layers slab of SnP 2 O 7 exposing two (200) surfaces (the z axis is parallel to the original 200 vector of the SnP 2 O 7 Pa-3 cubic unit-cell). The "bottom" surface is fully covered by ºP-OH and ºM-OH hydroxyl groups. The "top" surface shows arrays of HO-Sn-CH 2 O surface species (noted ºM(O=CH 2 )-OH), POH species (i.e. ZOH with Z=P, noted ºP-OH) and coordinatively unsaturated surface Sn cations (i.e. Z with Z=Sn, noted ºM). Sn centers in blue, P in magenta, O in red, Carbon in grey, Hydrogen in white. Perspective view.

Table 1

 1

	In the case of orthophosphates, all O atoms

reports the space groups and a choice of geometric features obtained by inspection in real space of the unitcells of reference materials matching the observed powder X-ray diffraction spectra of our catalysts as presented in Supporting Information.

Table 1 :

 1 (continued)

	Solid	SG	Ref	d1aC	d2aC	qaC	10 2 .P/VC
				Å	Å	°	DFT Å -3
	Mn2P2O7	C2/m	[26]	1.552	2.20	128.5	1.54
	CrPO4	Cmcm	[27]	1.557	2.02	128.3	1.56
	Co3(PO4)2	P-1	[28]	1.557	2.09	124.5	1.23
	Ni3(PO4)2	P1	[29]	1.566	2.13	123.1	1.41
	FePO4	P3121	[30]	1.548	1.88	141.8	1.12
	Mg3(PO4)2	P2_1/c	[31]	1.553	2.00	146.0	1.22
	AlPO4	P3221	[32]	1.540	1.75	145.8	1.22
	BPO4	I-4	[33]	1.543	1.47	134.5	1.52
	TiP2O7	Pa-3	[34]	1.520	1.94	167.1	1.45
	SiP2O7	Pa-3	[35]	1.511	1.76	158.7	1.73
	GeP2O7	Pa-3	[36]	1.514	1.87	157.3	1.60
	SnP2O7	Pa-3	[37]	1.514	2.03	161.1	1.38
	a-2) Computed geometries and bond strengths
	Starting from reference structures unit-cells, computed geometries were obtained from first-
	principles according to the method described previously. Experimental and computed minimum
	energy unit-cells differ slightly, as expected from the level of approximation involved in our
	DFT calculations. Computed values are also reported in Table 1, so as to allow easier
	comparisons with experimental values.

Table 2 :

 2 Computed bond energies in kJ.mol -1 for DFT optimized reference structures of

	-3-1 in

Ni 3 (PO 4 ) 2 Co 3 (PO 4 ) 2 Mn 2 P 2 O 7 CrPO 4 FePO 4 Mg 3 (PO 4 ) 2 AlPO 4 BPO 4 TiP 2 O 7 SiP 2 O 7 GeP 2 O 7 SnP 2 O 7

	1,00E+03
								Theoretical rate
								W(EMO)+W(EPO)
								Experiments
								MP2O7
	1,00E+02							Theory W(EMO)
								Theory W(EPO)
	1,00E+01
	1,00E+00
	0	50	100	150	200	250	300	350	400	450	500

E MO/PO (kJ.

mol -1 ) Rate of CH 4 oxidation 10 9 x mol.m -2 .s -1

Table 1 )

 1 . Error bars indicate estimated relative errors. The blue line is the Morse potential (equation 17) with parameters r0 = 1,53 Å, a = 20,97 Å -1 and De = 491,2 kJ.mol -1 .

	Mn 2 P

2 O 7 CrPO 4 Co 3 (PO 4 ) 2 Ni 3 (PO 4 ) 2 FePO 4 Mg 3 (PO 4 ) 2 AlPO 4 BPO 4 TiP 2 O 7 SiP 2 O 7 GeP 2 O 7 SnP 2 O 7

	-350
	-370
	-390
	-410
	-430
	-450
	-470
	-490
	-510
	-530
	1,5	1,51	1,52	1,53	1,54	1,55	1,56

Table 4 :

 4 bond orders of P-Oa bonds estimated according to equations (1) (BO A, parameter a fitted to 0,0474) and (2) (BO B) of reference[START_REF] Van Santen | On Shustorovich's Bond-Order conservation method as applied to chemisorption[END_REF], using bond lengths computed by DFT, and the optimal Morse potential shown on Figure5.

	Solid	BO A BO B
	Mn2P2O7	0,631 0,633
	CrPO4	0,571	0,572
	Co3(PO4)2	0,569 0,571
	Ni3(PO4)2	0,472 0,474
	FePO4	0,688 0,690
	Mg3(PO4)2	0,618 0,620
	AlPO4	0,815 0,816
	BPO4	0,768 0,769
	TiP2O7	1,248 1,246
	SiP2O7	1,484 1,480
	GeP2O7	1,401	1,399
	SnP2O7	1,401	1,399

Mn 2 P 2 O 7 C 2/m Calculated XRD diagram

	Exp d (A°) exp I CrPO 4 Cmcm CrPO 4 Main lines Co 3 (PO 4 ) 2 Main Lines Ni 3 P 2 O 8 Main lines	Calc d (A°) Calc I
	5,2 4,45 Exp d (A°) Exp d(A°) Exp d(A°) 3,1 2,95 2,6 2,37 2,17 2,07 1,85	4,30 4,33 5,82 3,88 3,51 3,24 2,58 2,48 2,31 2,15 2,03 1,975 4,09 3,81 3,85 3,56 3,44 3,34 3,21 2,82 3 2,68 2,9 2,45 2,78 2,30 2,6 2,24	8 exp. I exp I exp I 6 100 74 45 100 50 20 1 18 18 10 41 95 51 50 100 54 5 5 37 30 63 30 22 100 100 23 24 13 46 3 51 20 33 7 83	Mn 16 9 CrPO 4 Cmcm Experimental XRD diagram 2 P 2 O 7 Experimental Diagram 5,16 Calc. d (A°) Calc. I Calc d (A°) Calc I Calc d (A°) Calc I 4,43 4,30 5,83 34 4,34 3,08 100 2,94 53 2,62 18 2,37 4 2,17 18 2,07 16 1,85 15 3,88 3,52 3,28 2,59 2,49 2,31 2,15 2,03 1,97 3,78 54 4,07 3,66 8 3,85 3,43 100 3,44 2,91 27 3,22 2,74 45 3,01 2,49 86 2,93 2,30 11 2,80 2,62 2,23 12
			1,938 2,51	20 33		1,94 2,51
			1,845 2,48	4 50		1,85 2,44
	3 4 4 5 6 5 6 7	1	y = 0,9964x R² = 0,99987 3,5 Mn 2 P 2 O 7 exp. vs calc. reticular distances 1,5 2 2,5 3 1,531 12 1,53 1,502 5 1,50 1,87 3 1,83 13 1,84 1,87 1,76 28 1,681 2 1,637 20 1,64 1,621 5 1,62 2,03 10 1,93 7 1,93 2,05 R² = 0,99788 1,68 2,1 2 2,10 y = 0,9898x 1,76 2,23 7 2,23 1,828 6 1,84 2,3 3 2,30	4	4,5	5
	1 2 2 3		1,443 1,63 1,473 1,69 1,484 1,74	7 10 CrPO 4 Cmcm Calculated XRD diagram 1,44 1,63 16 1,47 17 1,69 4 1,49 7 1,77
	0 1	0	1 1,433 1,59 1,383 1,56	10 3 37 10	2	3 1,59 1,43 1,38 1,56		4	5	6
	0	0	1,337 1,319 1,29 1,54 1,51 1,48 1	Mn2P2O7 6 6 10 10 7 7 2	1,34 1,32 1,29 1,54 1,51 1,48	3	Linéaire (Mn2P2O7) 4	5	6	7
			1,26 1,45	14 10		1,28 1,45
			1,202 1,39	12 7		1,20 1,39
			1,37	7		1,37
		0	1		2	3		4	5	6	7	8
		P 2 O 7 1	1,5		2	2,5		3	3,5	4	4,5	5
	Main lines
								S4 S5 S6

Ni 3 P 2 O 8 expt. vs calc. reticular distances

	FePO 4 Main lines Mg 3 (PO 4 ) 2 P2_1/c
	Exp d(A°)		exp I	Calc d (A°) Calc I
	100	4,31	45	4,42		29,3
		3,85	9,6	3,77		0,6
		3,49	6,2	3,48		43,6
		3,35	100	3,48	100,0 Mg 3 (PO 4 ) 2 P 2_1/c Calculated Cu K⍺ XRD diagram
		2,95	41,4	2,87		0,5
		1,905	36	1,89		16,0
			FePO 4 expt. vs calc. reticular distances
	5
				y = 0,9953x
	4			R² = 0,98834
	3 0
	10	15	20	25	30	35	40	45	50	55	60
	2
				Mg 3 (
	1
	0
	0,00		1,00	2,00	3,00		4,00		5,00
	500
	0
	10	15	20	25	30	35	40	45	50	55	60
					S15

PO 4 ) 2 P 2_1/c Experimental Cu K⍺ XRD diagram

	Mg 3 (PO 4 ) 2 Main lines
	Exp d (A°)	exp. I	Calc. d (A°) Calc. I
	5,82	20	5,58	10
	4,37	5	4,36	22
	4,14	12	4,12	23
	4,10	8	4,08	25
	3,86	17	3,84	83
	3,45	17	3,44	100
	2,43	10	2,42	25
	2,11	100	2,13	22
		Mg 3 (PO 4 ) 2 expt. vs calc. reticular distances
	7,00
	6,00			y = 1,0137x
	5,00			R² = 0,99589
	4,00
	3,00
	2,00
	1,00
	0,00
	0,00	1,00	2,00	3,00	4,00	5,00		6,00
	0
	0	1	2	3	4	5	6	7	8

AlPO 4 P3221 Calculated XRD diagram AlPO 4 Experimental Diagram

	AlPO 4 Main lines BPO 4 I-4 BPO 4 Main lines TiP 2 O 7 Main lines
	Exp d (A°) Exp d (A°) Exp d (A°)	exp. I exp. I exp. I	Calc. d (A°) Calc. I Calc. d (A°) Calc. I Calc. d (A°) Calc. I
	100		3,34 2,15 3,63 4,52 3,32 3,91	72 76 100 11 5 100		3,37 2,14 3,65 4,50 3,33 3,90	100 10 100 8 6 100
	90		2,09 3,06 3,50	5 79		2,10 3,09 3,49	1,5 10 76
			1,82 2,25 3,20	41 5 57		1,83 2,26 3,18	15 30 72
	80		1,64 2,17 2,76	2 12		1,64 2,19 2,76	0,5 AlPO 4 P3121 Calculated XRD diagram 1 15
	70		1,35 1,97 2,36 1,86 2,18	2 18 10 8		1,38 1,98 2,35 1,88 2,16	16 3 14 10 9
	60		1,92	6		1,89	11
	50		1,76	AlPO 4 expt. vs calc. reticular distances BPO 4 expt. vs calc. reticular distances 8 1,74 12
	4,00 4,00 40 3,00 20 30 3,00 5				y = 0,9939x y = 0,9938x R² = 0,99953 R² = 0,99987 y = 0,9961x
	10 4					R² = 0,99994
	2,00 0 2,00
	3	0	1		2	3	4	5	6	7	8
	1,00 1,00
	2
	0,00 0,00 1	0,00 0,00		1,00 1,00	2,00 BPO 4 Experimental XRD diagram 3,00 2,00 3,00	4,00 4,00
	0
	0,00		1,00		2,00	3,00		4,00		5,
	0	1		2	3	4	5	6	7	8
						S18 S19 S20

00 TiP 2 O 7 expt. vs calc. reticular distances

	SiP 2 O 7 Main lines GeP 2 O 7 Main lines SnP 2 O 7 Pa-3 SnP 2 O 7 Main lines
	Exp d (A°) Exp d (A°) Exp d (A°)	exp. I exp. I exp. I	Calc. d (A°) Calc. I Calc. d (A°) Calc. I Calc. d (A°) Calc. I
	1200 1400 2⍬	I	3,69 3,33 3,02 4,04 800 4,54 3,99 3,60 700 3,89 3,48	100 100 100 81 50 72 100 100 27		3,74 3,34 3,05 4,34 4,56 4,34 3,76 3,95 3,53	67 100 88 5 44 3 100 100 SnP 2 O 7 Experimental XRD diagram ( Cu K⍺ ) Sn 24
				2,63 3,14 600 3,17	20 33 19		2,64 GeP 2 O 7 Experimental XRD diagram( Cu K⍺ ) 7 3,07 25 3,22 23
	1000		2,27 2,14 2,88 2,84 2,75 500 2,36	60 20 41 57 25 22		2,25 2,16 2,66 2,51 2,79 2,38		3 2 46 5 35 42
				2,07 2,31 400 2,25	60 36 9		2,07 2,27 2,28		10 18 8
				1,81 2,11 1,75	50 24 6		1,81 2,17 1,76		18 19 12
				1,67 1,89 300	40 12		1,67 2,09		11 4
				2,01 1,84 200	15 17 SnP 2 O	2,01 1,82		6 11
			5	1,71 0 1,79 100	16 SiP 2 O 7 expt. vs calc. reticular distances 1,73 26 y = 1,0103x 16 1,77 3
	3 100 5 4 0	10 3 4	15 1,63	R² = 0,99949 30 35 GeP 2 O 7 expt. vs calc. reticular distances 20 25 40 45 2 Q y = 0,9943x 15 1,64 11 20 30 40 R² = 0,99973	50 50	55	60	60
	2 4		2				GeP 2 O 7 Pa-3 Calculated XRD diagram R² = 0,95809 y = 0,9787x
									SnP 2 O 7 Pa-3 Calculated XRD diagram
	1 3		1
	2		0	0	1		2			3		4	5
	0
	0,00 1	1,00			2,00			3,00		4,00
	0
	0	10	15	20	25	30	35	40	45	50	55	60
			0		1		2		3			4	5
				10	15	20	25	30	35	40	45	50	55	60
								S27

7 expt. vs calc. reticular distances BET Areas

  The following table provides the BET surface areas measured on the considered phosphates after calcination at 973 K:

	SI-3 Supplementary Tables
	Solid	SG	BET area
			m 2 .g -1
	Mn 2 P 2 O 7	C2/m	6
	CrPO 4	Cmcm	65
	Co 3 (PO 4 ) 2	P-1	91
	Ni 3 (PO 4 ) 2	P1	29
	FePO 4	P3121	17
	Mg 3 (PO 4 ) 2	P2_1/c	46
	AlPO 4	P3221	15
	BPO 4	I-4	12
	TiP 2 O 7	Pa-3	27
	SiP 2 O 7	Pa-3	35
	GeP 2 O 7	Pa-3	42
	SnP 2 O 7	Pa-3	76

Table SI -

 SI 3-1: Bader charges in electronic units for DFT optimized reference structures of phosphate bulk catalysts. 𝑄 7 : averaged Bader charge on M atoms; 𝑄 78 : formal charge of M ion; 𝑄 9 : averaged Bader charge on P atoms; 𝑄 :; : averaged Bader charge on Oa atoms; averaged 𝑄 :< Bader charge on Ob atoms. Formal charges for P and O are +5 and -2 respectively.Table SI-3-2: Experimental specific total acidities 𝐻 > in µmol.m -2 , catalytic oxidation rates w in 10 9 mol.m -2 .s -1 and selectivities in moles % of oxidation products for phosphate catalysts at 873 K, atmospheric pressure and molar ratio O 2 /CH 4 = 0.5. Acidities and selectivities in formaldehyde are correlated.Table SI-3-3: Comparison of experimental and simulated catalytic oxidation rates w in 10 9 mol.m -2 .s -1 and selectivities in moles % of oxidation products for phosphate catalysts at 873 K, atmospheric pressure and molar ratio O 2 /CH 4 = 0.5 . Descriptors E Moa (computed by DFT, see

	Solid	𝐻 > w	%𝐶𝐻 A 𝑂 %CO %𝐶𝑂 A
	Mn 2 P 2 O 7	0.09 185 1	15	84
	Solid CrPO 4 Mn 2 P 2 O 7 Co 3 (PO 4 ) 2 0.83 35.6 6 𝑄 7 𝑄 78 𝑄 9 0.35 76.3 3 1.533 2 3.687 -1.478 -1.568 𝑄 :; 19 78 𝑄 :< 23 71 CrPO 4 1.815 3 Ni 3 (PO 4 ) 2 1.30 25.4 10 26 64 3.652 -1.368 Co 3 (PO 4 ) 2 2 2 1.312 2 FePO 4 2.10 16.4 15 33 52 3.656 -1.402 Ni 3 (PO 4 ) 2 1.236 2 Mg 3 (PO 4 ) 2 3.30 11.8 24 36 40 3.616 -1.368 FePO 4 1.703 3 AlPO 4 5.00 8.76 38 35 27 3.703 -1.351 Mg 3 (PO 4 ) 2 1.697 2 BPO 4 6.90 6.28 51 32 17 3.694 -1.558 AlPO 4 2.492 3 TiP 2 O 7 10.6 4.85 79 22 9 3.756 -1.562 BPO 4 2.380 3 SiP 2 O 7 13.6 3.55 83 12 5 3.756 -1.534 TiP 2 O 7 2.380 4 GeP 2 O 7 15.2 2.94 86 11 3 3.751 -1.391 -1.538 SiP 2 O 7 3.262 4 3.772 -1.543 -1.548 SnP 2 O 7 16.5 2.45 96 3 1
	GeP 2 O 7	2.640 4		3.776 -1.438 -1.564
	SnP 2 O 7	2.692 4		3.756 -1.441 -1.556

Table 2 )

 2 and 𝐻 > (experimentally measured, see Table4) are inputs in the kinetic model.As presented in the main text, we assume that methane oxidation can be decomposed in the following elementary steps for each phosphate catalyst:

	Solid	w exp w sim %𝐶𝐻 A 𝑂 %𝐶𝐻 A 𝑂𝑠𝑖𝑚
	Mn 2 P 2 O 7	185 185.0 1	4.1
	CrPO 4	76.3 94.4 3	4.8
	Co 3 (PO 4 ) 2 35.6 61.8 6	7.4
	Ni 3 (PO 4 ) 2 25.4 25.3 10	11.3
	FePO 4	16.4 25.6 15	13.8
	Mg 3 (PO 4 ) 2 11.8 11.6	24	23.7
	AlPO 4	8.76 8.41 38	37.1
	BPO 4	6.28 6.18 51	56.6
	TiP 2 O 7	4.85 8.99 79	77.5
	SiP 2 O 7	3.55 14.6 83	86.8
	GeP 2 O 7	2.94 53,4 86	87.7
	SnP 2 O 7	2.45 62,4 96	90.7
	Solid	%CO %COsim %𝐶𝑂 A %𝐶𝑂 A 𝑠𝑖𝑚
	Mn 2 P 2 O 7	15	16.8	84	79.1
	CrPO 4	19	21.5	78	71.4
	Co 3 (PO 4 ) 2 23	21.1	71	71.5
	Ni 3 (PO 4 ) 2 26	23.7	64	65.0
	FePO 4	33	33.9	52	52.3
	Mg 3 (PO 4 ) 2 36	36.5	40	39.8
	AlPO 4	35	34.9	27	28.0
	BPO 4	32	27.3	17	16.1
	TiP 2 O 7	22	11.4	9	11.1
	SiP 2 O 7	12	5.9	5	7.7
	GeP 2 O 7	11	3.2	3	9.1
	SnP 2 O 7	3	2.3	1	7,0

  further assume: i) that forward step 2, the entrance channel for carbonaceous species to the adsorbed state at the catalyst surface is the overall limiting step; ii) that the comparatively very fast reaction rates for forward steps 3 to 8 are of first order with respect to the respective reacting species. It follows that fractional coverages by other carbonaceous species than ZOHCH3 remain always small fractions of 𝜃 H:MNMO which can be considered as a bundled coverage by carbon in the completion for available surface sites. Denoting 𝑘 & the rate constant for step 𝑖, and since the rates of consecutive reverse steps 3, 4 and 6 are considered negligible before the corresponding forward rates, we obtain from simple material balances:

	Next, we
	H H U 1 + 𝑃 : J 𝐾 L	L A 1 + 𝑃 𝐶𝐻 4 𝐾 2	eq. SI-4_6
	Introducing fractional coverages 𝜃 H: H U ,	H:MNMO H U	,	H:M H U ,
	respectively, we get :
					ZL
	𝜃 H = 1 + 𝑃 : J 𝐾 L	L A 1 + 𝑃 𝐶𝐻 4 𝐾 2	= 𝐷 ZL	eq. SI-4_7
	𝜃 H: = 𝜃 H 𝑃 : J 𝐾 L	L A	eq. SI-4_8
	𝜃 H:MNMO = 𝜃 H 𝜃 H: 𝑃 𝐶𝐻 4 𝐾 2	eq. SI-4_9
	𝜃 H:M = 𝜃 H 𝜃 H:	𝑃 𝐻 2 𝑂 S T	1 2	eq. SI-4_10

H , 𝜃 H: , 𝜃 H:MNMO , 𝜃 H:M for H H U ,

  NM P = 𝑊 NM P7 + 𝑊 NM P9 + 𝑊 NM P |} = 𝑘 A 7 𝜃 7:MNMO + 𝑘 A 9 𝜃 9:MNMO + 𝑊 NM P |} eq. SI-4_28where the superscripts M and P refer to the choice of descriptors 𝐸 7: or 𝐸 9: in equations SI-4_23 and SI-4_24. The thermal rate, a 16 th free parameter, sets a lower limit to the rate in absence of catalyst, and imparts more physical realism by allowing a smooth convergence of the predicted residual rate for vanishing or very high values of the descriptors 𝐸 7: or 𝐸 9: .

	Accordingly, predicted selectivities become:
	𝑆 ‰ =	Š `QP ‹ OE `QP	𝑆 ‰ 7 +	• Š `QP OE `QP	𝑆 ‰ 9

Table SI -

 SI 4-1: Best fitted parameters of the kinetic model. Other inputs are operating conditions, average site density for phosphate catalysts 1.25 10 19 site.m -2 ; T = 873 K, [CH 4 ]= 10 -2 mol.l -1 ; [O 2 ]= 5.10 -3 mol.l -1 and Boltzmann and Planck fundamental constants k B and h

	Δ𝑆 L >	73.5	J.mol -1 .K -1
	Δ𝑆 A >	-18,19 201.04	J.mol -1 .K -1
	-18,19 Δ𝑆 ± -56.78	J.mol -1 .K -1
	𝐸 L >	-18,19 0.00	kJ.mol -1
	𝐸 A >	29.96 kJ.mol -1
	𝐸 A ±	kJ.mol -1 179.12 kJ.mol -1
	𝛼 A -0.16 -kJ.mol -1 𝐸 h ± 4.29 kJ.mol -1 kJ.mol -1 𝛼 h 0.32 kJ.m 2 .µmol -2
	𝐸 • ±	kJ.mol -1 15.62 kJ.mol -1
	𝛼 • -0.04 -
	𝐸 i ±	kJ.mol -1 52.21 kJ.mol -1
	𝛼 i -0.29 -
	𝐸 j ±	kJ.mol -1 28.98 kJ.mol -1
	𝛼 j -0.36 -
	𝑊 NM P |}	2.00	kJ.mol -1 10 9 .mol.m -2 .s -1

ASSOCIATED CONTENT Supporting Information. One pdf file including four sections: SI-1: Catalysts characterization (X-Ray diffraction, BET areas, IR spectra of phosphate catalysts); SI-2 : Bench-Scale Catalysts testing equipment; SI-3 : Supplementary Tables; SI-4 : Kinetic Model; SI-5 : Molecular models of key surface intermediates along the catalytic pathways.

Author Information

IR Spectra of catalysts

In the following spectra, the absorbance is plotted as function of wavenumber. Spectrum 1 corresponds to the freshly prepared catalysts, and spectrum 2 to the same material after calcination at 1073K.

All catalyst present bands in the 3400-3800 cm -1 range of wavenumbers, specific of n(OH) vibration, and revealing the occurrence of surface hydroxyls groups.

The bands in the vicinity of 3660 cm -1 have been assigned in ref [START_REF] Ramis | Surface Acidity of the Layered Pyrophosphates of Quadrivalent Ti, Zr, Ge, and Sn and Their Activity in Some Acid-catalysed Reactions[END_REF] of main text to the OH stretching of free surface POH groups. The feed rate of reaction mixture is set between 3 and15 Nl.hr -1 , while the recycle flow-rate is maintained at 150 Nl.hr -1 so as to ensure the condition of equivalence with a perfectly stirred differential reactor.

(2) Gomonai, V.I., Doctoral (Chem.) Dissertation, Uzhgorod State Univ., 1990.

of the rates of oxidation 𝑤 NM P . This first function is minimized with respect to the eight first parameters ∆𝑆 L > , ∆𝑆 A > , 𝐸 L > , 𝐸 A > , ∆𝑆 ± , 𝐸 z ± and 𝛼 z for l = 2, and 𝑊 NM P |} , using 8 experimental rates of oxidation corresponding to Mn 2 P 2 O 7 and the 7 orthophosphates catalysts (see Table SI -3-3) and constraining 𝐸 A ± to remain positive in the range of 𝐸 …: values spanned by these catalysts.

The four remaining pyrophosphates are excluded from this first training set since they appear obviously out of the regular volcano trend observed, when plotting experimental 𝑤 NM P versus 𝐸 …: . We obtained for this step a relative rmsd of 6.3 %.

In the second step, a second objective function is built as the root mean square deviation of predicted from experimental selectivities, and minimized with respect to the 8 remaining parameters 𝐸 z ± and 𝛼 z for l = 5, 6, 7, 8 using all 24 independent experimental selectivities available (see Table S-3-3), constraining 𝐸 z ± for l = 5, 6, 7, 8 to remain positive in the range of 𝐸 …: values spanned by these catalysts. We obtained for this step a relative rmsd of 1.8 %.

Rate constants given by equation SI-4_17 are converted from theoretical units of molecule.site -1 .s -1 into experimental units of 10 9 mol.m -2 .s -1 using an average site density of 1,25.10 19 m -2 estimated according to the volumic density of P atoms reported in Table 1 of main text.

Best fitted parameters are listed in the following Table SI-4-1:

SI-5 Molecular models of key surface intermediates along the catalytic pathways

The atomistic models presented in this section were built with the MAPS platform from SCIENOMICS (http://www.scienomics.com) provided to HT and LRS according to their participation to the Scienomics Group of Scientific Excellence (SGSE). The configurations shown were not optimized.