Annamaria Quaranta

Anja Krieger-Liszkay

Andrew A Pascal

François Perreau

Bruno Robert

Mikas Vengris

Manuel J Llansola-Portoles

Bruno Rob- Ert †

Singlet fission in naturally-organized carotenoid molecules

Keywords: singlet exciton fission, daffodil chromoplast, aggregated carotenoids, triplet state

à la diffusion de documents scientifiques de niveau recherche, publiés ou non, émanant des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

INTRODUCTION

The bright colours of many flowers and fruits originate from the carotenoids accumulated in specialized plastids in the plant cell, termed chromoplasts. Chromoplasts are generally produced in flower petals and ripe fruit through a remodeling of the photosynthetic chloroplast 1 , involving chlorophyll degradation, down-regulation of genes encoding for photosynthetic proteins, and up-regulation of stress-response and carotenoid-synthesis genes (although differentiation may also occur from other plastid types 2 ). In mature chromoplasts, the carotenoid pigments are sequestered in the form of crystal-like structures (crystalloids), granules, or stromal plastoglobules 3 . The main function of the chromoplast is thought to be secondary signaling, attracting animals to flowers and fruits to aid pollination and seed dispersal, respectively. They are also the site of protein, lipid and carbohydrate biosynthesis 4 , and play an important role in fruit and flower scent (for instance, one of the main components of the tomato scent is a degradation product of chromoplastic lycopene) 5 . It is interesting to note that the chloroplast/chromoplast transition is in some cases a reversible process, such that chromoplasts can re-acquire photosynthetic activity and differentiate back into chloroplasts 1,2 .

Carotenoid molecules are tetraterpenoid derivatives formed by eight isoprene units. They display a wide variety of functions in photosynthetic tissues, some of which are also observed throughout nature. They play an important role in harvesting photons in the blue-green range in photosynthesis [START_REF] Britton | Natural Functions[END_REF] , and ensure efficient photoprotection of the photosynthetic apparatus -through the quenching of chlorophyll triplet [START_REF] Jahns | [END_REF] and singlet 8 excited states, and by short-circuiting lipid peroxidation 9 . Carotenoids exhibit strong electronic absorption corresponding to a transition from the ground state to the second excited state (S0S2), while transition to the first excited state (which has the same symmetry as the ground state) is one-photon forbidden 10 . After photoexcitation, isolated carotenoids decay rapidly to the ground state in a few picoseconds 11,12 . Although chromoplasts are not expected to display any specific functional photoactivity, we recently discovered that the lycopene crystalloids in tomato chromoplasts do indeed exhibit photophysics that are not observed in isolated carotenoids 13 . This self-associated lycopene dissipates excitation energy via singlet fission, a mechanism we tentatively proposed to be photoprotective. In the singlet exciton fission process, one singlet excited state splits into two excited triplets (T1 + T1), fulfilling the energetic condition E(S1)  2E(T1) [14][15][16] . Carotenoids in tightly-interacting assemblies in vitro have been shown to generate triplets by singlet fission -such as aggregated astaxantin 17,18 , β-carotene 19,20 , and zeaxanthin [21][START_REF] Wang | Proceedings of SPIE[END_REF][START_REF] Billsten | [END_REF] -but up to now, tomato chromoplasts represent the only known example of this process in a biological structure 13 .

We investigated the lutein/violaxanthin photochemistry in daffodil <Narcissus pseudonarcissus L.> chromoplasts [START_REF] Camara | International Review of Cytology[END_REF][START_REF] Liedvogel | [END_REF] , to determine whether singlet fission is unique to the lycopene in tomato crystalloids, or is a more general chromoplast phenomenon. Daffodil chromoplasts accumulate these carotenoids in a concentric membrane system, forming condensed structures involving carotenoids and acyl lipid molecules (they belong to the class of so-called membranous chromoplasts) [START_REF] Camara | International Review of Cytology[END_REF][START_REF] Liedvogel | [END_REF] . These structures allow close interactions between carotenoids, which perturb the absorption properties of lutein and violaxanthin. We show that singlet fission indeed occurs in these plastids, suggesting that this process may be common to all carotenoidcontaining chromoplasts. Our experiments reveal the singlet fission mechanism in unprecedented detail. The information on specific energy levels of intermediate species allows us to identify the entangled triplet pair, and to pinpoint the spectral signatures of the disentangled triplet excited state(s) of the excitonic lutein/violaxanthin system.

MATERIALS AND METHODS

Chromoplasts of the inner coronae of fully-opened daffodil flowers (Narcissus pseudonarcissus L. cv Dutch Master), at developmental stage IV 26 , were isolated according to the method of Liedvogel et. al. [START_REF] Liedvogel | [END_REF] using a step gradient in which chromoplasts accumulate at the 15 % / 30 % sucrose interphase 27 . The isolated chromoplasts were suspended in 100 mM Tris/HCl pH 7.4, 10 mM MgCl2, 1 mM 1,4-dithioerythritol. The preparations were stored at -80 °C.

For pigment determination, the preparation was centrifuged, the supernatant discarded and the pelleted chromoplasts were resuspended in 2 mL acetone. 5 µL was then subjected to liquid chromatography at a flow rate of 0.5 mL min -1 , using an HPLC column (Uptisphere Strategy C18-HQ 250 × 3 mm 3 µm, Interchim, Montluçon, France) with a liquid chromatography chain HPLC/UV (Shimadzu, Tokyo, Japan) which included two pumps (LC-20AD), a sample manager (SIL-20AC HT), a column oven (CTO-20A) and a UV diode-array detector (UVSPD-M20A) monitoring the absorbance at 450 nm. Mobile phase A was composed of 10 % water and 90 % acetonitrile with 0.5 % acetic acid, mobile phase B of ethyl acetate containing 0.5 % acetic acid, and the elution gradient was as follows: 1 min 90 % A : 10 % B, then a linear gradient to 5 % A : 95 % B after 25 min. Concentrations of each pigment were calculated based on their peak areas, calibrated using injected standards or (for some minor pigments) with known absorption factors for each compound at 450 nm.

Absorption spectra were measured using a Varian Cary E5 scanning spectrophotometer. Low temperature measurements were performed in a helium bath cryostat (Maico Metriks, Tartu, Estonia); 60 % glycerol (v/v) was added to the sample prior to freezing.

Circular dichroism spectra were recorded on an Aviv 215 spectropolarimeter (Aviv Biomedical). Near-UV CD spectra were recorded between 300 and 600 nm using a square cell with a 1 cm path length, and an averaging time of 1 s per step. Three consecutive scans from each sample were merged to produce an averaged spectrum, which was then corrected against buffer baselines measured under the same conditions. Dichroism value is measured in mdeg.

Resonance Raman spectra were recorded at 77 K with laser excitations obtained from Coherent Ar+ (Sabre) and Kr + (Innova 90) lasers. Output laser powers of 10-100 mW were attenuated to < 5 mW at the sample. Scattered light was collected at 90֯ to the incident light, and focused into a Jobin-Yvon U1000 double-grating spectrometer (1800 grooves/mm) equipped with a red-sensitive, back-illuminated, LN2-cooled CCD camera. The sample stability and integrity were assessed based on the stability of the Raman signal.

Nano-to-millisecond transient absorption experiments were performed on an Edinburgh Instruments LP920 Flash Photolysis Spectrometer system, incorporating a Continuum Surelite OPO. The OPO was pumped by a Q-switched Nd:YAG laser operating at 355 nm. The LP920 system is equipped with a 450 W Xenon arc lamp as the probe for the transient absorption measurements; this lamp was pulsed in the time range 10 ns to 100 µs. Detection in the LP920 system was performed either via a 500-nm-blazed Czerny-Turner monochromator (bandwidth 1-5 nm) coupled with a Hamamatsu R928 photomultiplier tube (kinetics mode), or via a 500-nm-blazed spectrograph (bandwidth 5 nm) coupled with a water-cooled ICCD nanosecond Andor DH720 camera (spectral mode). Excitation energy was approximately 15 mJ / pulse. Femto-to-nanosecond time-resolved transient absorption spectra were recorded on a setup described in detail elsewhere 28 . Pump and probe pulses were derived from an amplified Ti:Sapphire laser (Libra, Coherent). Tunable pump pulses with center wavelength at 490 nm were obtained from an optical parametric amplifier (Topas, Light Conversion). White light continuum generated in CaF2 from the fundamental output of the laser was used to probe the absorption changes in the 350-750 nm range. The time resolution of the instrument is ca. 120 fs. Excitation energies were set to approximately 200 nJ / pulse.

RESULTS AND DISCUSSION

The carotenoid content of daffodil chromoplasts was determined by HPLC (Supporting information, Figure S1). The carotenoids with significant concentrations are lutein and violaxanthin, with a ratio circa 5:4. Lutein and violaxanthin in cyclohexane both present the typical absorption spectrum for carotenoids, with four bands corresponding to the vibronic sublevels 0-3, 0-2, 0-1, and 0-0 of the S0S2 electronic transition -at 395, 416, 440 & 474 nm and 393, 422, 452 & 471 nm, respectively (Fig. 1a) 29 . Figure 1b shows the absorption spectrum of broken daffodil chromoplasts in 99 % ethanol, which disrupts the chromoplast and solubilizes the carotenoids. This spectrum, with maxima at 394, 418, 442 & 471 nm, is essentially identical to a linear combination of isolated lutein and violaxanthin using the proportion obtained from the HPLC analysis (supporting information, Figure S2). The absorption spectrum of intact chromoplasts at room temperature exhibits broad vibronic transitions peaking at 393, 422, 447 and 481 nm (Fig. 1c). The 0-0 band at 481 nm is redshifted relative to the spectrum of solubilized chromoplasts (note that the absorption of lutein may be observed as high as 485 nm in high polarizability solvents 30 ). The intensity of the 0-3 and 0-2 electronic transitions is significantly higher in chromoplasts than in monomeric carotenoids (cf. Fig. 1a,b,c), which may indicate some perturbation of the electronic transition caused by packing interactions 31 . The absorption spectrum of daffodil chromoplasts recorded at cryogenic temperature (4.2 K) reveals a complex picture, indicating the presence of several carotenoid populations (Fig. 1d) which can be resolved by second derivative analysis (Fig. 1e). The red part of the absorption spectrum shows the presence of two species with close electronic transitions: a dominant species (species 1) absorbing at 422, 451 & 487 nm, and a minor one (species 2) with absorption at 430, 463 & 500 nm (0-2, 0-1 & 0-0, respectively). An additional band is present at 402 nm whose second derivative shows no clear structure. This suggests that this 402-nm band is not a 0-3 vibronic transition of species 1 & 2, but could rather arise from a third, blue-shifted carotenoid population (species 3). Overall, the absorption at 4.2 K presents a characteristic bathochromic shift relative to the spectrum at room temperature 32,33 (note red-shift of the major 0-0 band from 481 to 487 nm). The circular dichroism of daffodil chromoplasts at room temperature (Fig 1f) exhibits chirality across the whole spectral range, with a dominant S-shaped feature on the blue side of the envelope (maximum 397 nm, minimum 375 nm). The CD spectrum of solubilized chromoplasts displays no chirality (data not shown). The observed chirality for intact chromoplasts indicates that all three species are constituted from aggregated (interacting) carotenoid molecules. Species 3 has a blue shifted 0-0 transition, around 402 nm at 4.2 K, and displays a large CD signature. This is representative of H-aggregates, where exciton coupling is strong and extends over multiple chromophores. The absorption of such strongly-coupled H-aggregates collapses into a single electronic transition, with the large hypsochromic shift observed here (circa 5000 cm -1 relative to monomeric violaxanthin and lutein) 34 . Species 1 & 2 conserve the vibronic peaks of monomeric carotenoids, and exhibit only a small bathochromic shift relative to them. Carotenoid aggregates that have a predominant bathochromic shift have been traditionally interpreted as J-aggregates [START_REF] Köhn | Carotenoids: Volume 4: Natural Functions[END_REF][START_REF] Wang | [END_REF] . However, species 1 & 2 do not resemble the typical narrow exciton-coupled bands expected for strongly-coupled J-aggregates 17,18,37 . Their absorption and CD spectra more closely resemble model spectra of weaklycoupled H-aggregates of lutein 31 . Further studies are necessary to determine the precise molecular interaction(s) involved in species 1 & 2. We used resonance Raman spectroscopy 38 at three different excitation wavelengths to characterize the three carotenoid species in daffodil chromoplasts. The high selectivity of this technique provides information for pigment populations individually, by exciting close to their respective absorption maxima. Figure 2 shows the resonance Raman spectra at 77 K for isolated all-trans lutein and violaxanthin in cyclohexane, measured at 488.0 nm, and for intact chromoplasts in buffer excited at 488.0, 501.7 and 406.7 nm. The position of the 1 vibrational mode depends on the length of the π-electron conjugated chain in isolated carotenoids 30,39 . This band occurs at 1530.2 and 1529.1 cm -1 for violaxanthin and lutein in cyclohexane, respectively. The 1 frequency in daffodil chromoplasts exhibits a clear dependence on the excitation wavelength, peaking at 1530.4, 1526.7 & 1528.0 cm -1 for excitation at 488.0, 501.7 & 406.7 nm, respectively (see figure 2d). These different frequencies support the presence of three non-identical carotenoid species -note that the three excitation wavelengths correspond to the absorption peaks of species 1-3 (487, 500 & 402 nm, respectively). Due to the excitonic nature of aggregates, it is not possible to relate the 1 frequency to that of monomeric violaxanthin or lutein. However, the position of the 1 mode permits us to confirm the presence of more than one aggregate species with very close absorption.

Around 1160 cm -1 , a group of bands termed 2 constitutes a fingerprint for the isomerization state (cis/trans) of the carotenoid carbon backbone [40][41][42] . Chromoplasts display the same modes at all wavelengths used, with a dominant band at 1058 cm -1 and two satellites circa 1188 & 1214 cm -1 . The absence of a vibrational mode at 1130 cm -1 in daffodil chromoplasts, as for the isolated standards, confirms that both lutein and violaxanthin are in the all-trans configuration 40- 42 . Observed around 1000 cm -1 , the 3 band arises from inplane rocking vibrations of the methyl groups attached to the conjugated chain, coupled with in-plane bending modes of the adjacent C-H's 43 . The 3 mode is a fingerprint of the conformation of carotenoid conjugated end-cycles 44 . This band can thus be used to track the presence of lutein and/or violaxanthin in the chromoplast spectra, as a result of differences in their end-rings. Lutein in cyclohexane shows a sharp asymmetric band peaking at 1003.1 cm -1 with a small shoulder ca. 1007.9 cm -1 , while violaxanthin presents a broad band peaking at 1006.5 cm -1 with a clear shoulder ca. 1002.9 cm -1 . Chromoplasts excited at 501.7 nm exhibit a sharp vibrational mode at 1004.4 cm -1 , resembling the spectrum of monomeric lutein, thus indicating that this carotenoid is dominant in species 2. At 488.0 & 406.7 nm, the 3 band is wider with a prominent shoulder extending up to 1007.2 cm -1 -indicating the presence of both violaxanthin and lutein in species 1 & 3 (see figure 2c). The 4 mode around 960 cm -1 arises from C-H out-of-plane wagging motions coupled with C=C torsional modes, 39 and this region is an indicator of out-of-the-plane distortions of the carotenoid conjugated chain 39,45 . In chromoplasts, this region contains small but clear modes at 951.8 and 963.8 cm -1 , which are also present for isolated lutein and violaxanthin. Thus, all three species appear to experience no significant distortions from their relaxed conformation. Figure 3 shows the time-resolved absorption data obtained for chromoplasts in the 100 fs -5 ns time window, measured with an excitation pulse centered at 488 nm to investigate the photochemistry of the major population of carotenoid aggregates (species 1). The dataset overview (Figure 3a) shows that different parts of the ground-state bleaching (GSB) region in the 350-475 nm range decay with markedly different kinetics, and the excited state absorption (ESA) band in the 500-600 nm region exhibits two prominent peaks that reach their maxima at different times. Both these trends can be observed in more detail in the kinetic traces, presented in Figure 3b. The two ESA bands evolve at different rates, with the 552-nm peak appearing and then decaying slightly before the one at 528 nm. The complex evolution of these ESA bands is clearly observed in time-gated spectra (Figure 3c). The initial S2 excited state, observed by its negative stimulated emission (SE) band visible above 505 nm in the 120-fs spectrum (black trace), relaxes very quickly -evolving into a broad, almost featureless band peaking around 560 nm (red trace, 300 fs). This broad feature can be attributed to a vibrationally-hot S1 state. The evolution of this state is different from that observed for isolated carotenoids -instead of Wavenumber (cm -1 )

Wavenumber (cm -1 )

 1  3  4  3  2  1 c) d) b) Lutein Violaxanthin Species 2
Species 1

Wavenumber (cm -1 ) Species 3

K

narrowing to a well-defined S1 ESA peak, a nearly-symmetrical doublet appears with maxima at 528 and 552 nm (green trace, 2 ps). These peaks decay at different rates, so that the 528-nm maximum dominates the 15-ps spectrum (blue trace). At later times the ESA band shifts further to the blue, where it partly decays within 100 ps (cyan trace).

However, a small fraction of this signal remains beyond the 5-ns time window of the experiment (eg. pink trace, 2-5 ns). The evolution at longer timescales is again different from that seen for monomeric carotenoids in solution -including violaxanthin and lutein -which return to the ground state in tens of ps 12 . Global analysis using an evolutionary model yields the Evolution-Associated Difference Spectra (EADS) shown in Figure 3d, reflecting a cascade of events with the following time constants: 150 fs  500 fs  8.1 ps  24 ps  250 ps  non-decaying. These first two EADS at 150 and 500 fs (black, red spectra) are similar to those observed for monomeric carotenoids, reflecting analogous evolution over the same timescales 12 . The first step corresponds to the decay of the S2 state (150 fs, black spectrum), somewhat obscured by cross-phase modulation between pump and probe. The second step corresponds to cooling of the resulting hot S1 state (500 fs, red spectrum), whereby ESA loses its red shoulder at 550-650 nm. The third EADS (8.1 ps, green spectrum) presents two peaks of similar intensity at 528 and 552 nm, while the GSB in the blue part of the spectrum remains unchanged relative to the 500-fs spectrum. In the fourth EADS (24 ps, blue spectrum) a substantial decrease of relative intensity of the 552-nm peak is observed. The fifth EADS (250 ps, cyan spectrum) presents a distorted shape of the GSB on its red wing, being attenuated by about 90 %, possibly indicating the presence of a transient species with ESA overlapping in this region and pushing the signal up. The peaks at 528 and 552 nm disappear completely, while a broad ESA band peaks at ca. 507 nm. This sixth EADS (ꝏ, pink spectrum) is similar to that at 250 ps and persists beyond 5 ns (time window of the experiment). It should be noted that the two excited species with ESA peaking at 528 and 552 nm decay in parallel, and they cannot be separated in this sequential model. This model is therefore unable to describe the dynamics after photoexcitation completely, while it is adequate to perform an initial assessment of the excited species that will lead to the design of a representative model (see target analysis, below and Fig 5). In monomeric carotenoids in solution, S1 decays to the ground state in a few picoseconds, with no long-living species. The ESA position and decay rate for S1 in monomeric violaxanthin and lutein are 532 nm & 24 ps, and 558 nm & 14 ps, respectively 12 . In addition, the GSB intensity for isolated carotenoids decays in proportion to the S1 decay. In daffodil chromoplasts, the unchanged intensity of GSB between 500 fs and 8.1 ps indicates that the recombination of excited states is somewhat compensated by the formation of additional excited species depleting the ground state. It suggests that the two ESA cannot be assigned to independent S1 states decaying to the ground state. In the literature, similar ESA doublets have been reported for aggregates of -carotene 20 and zeaxanthin [START_REF] Wang | Proceedings of SPIE[END_REF][START_REF] Billsten | [END_REF] . Hence, the red ESA may be assigned to S1 (generated by excitonic coupling) decaying in a few ps, while the blue ESA can be assigned to a triplet state decaying over a wide range of timescales, from 10s to 100s of ps. Hence, the ESA observed for daffodil chromoplasts at 528 nm may indicate the presence of triplet states.

In order to analyze the fate of the excitation at longer timescales, transient absorption spectra of daffodil chromoplasts in the ns-to-s range were measured in buffer solutions saturated with oxygen, air and argon (Fig. 4). All transient absorption kinetics can be fitted to a mono-exponential decay, and the gated spectra for chromoplasts at three different oxygen concentrations are nearly identical. They all show a wide negative region around 350-440 nm, attributable to GSB. In the 440-500 nm region, the ground state bleaching is mixed up with ESA, with negative contributions at 480 and 444 nm corresponding to the 0-0 and 0-1 vibronic peaks of the S0S2 electronic transitions, respectively. The ESA in the region 500-530 nm exhibits a peak at 505 nm, accompanied by a signal at 520 nm appearing as a peak in argon or a shoulder in air and oxygen. These signals decay with lifetimes of 180 ns, 700 ns, and 3.5 s in oxygen, air and argon, respectively. The dependence of this lifetime on oxygen concentration (Fig. 4), a molecule known to enhance inter-system crossing in carotenoid triplet states 46,47 , is consistent with the attribution of these ESA to T1Tn transitions. In monomers of the plant antenna protein LHCII, no significant difference is observed between lutein and violaxanthin triplets, in the native protein and the npq1lut2 mutant, respectively -T1Tn peaks at 508 nm in both cases 48 . The weak excitonic coupling of species 1 and 2 may introduce changes in their T1Tn ESA position. It therefore seems reasonable to assign the higher-energy T1Tn transition at 505 nm to associated lutein/violaxanthin triplets in species 1 (0-0 at 481 nm a room temperature), while associated lutein (species 2; 500 nm) contributes to the lower-energy triplet signal at 520 nm. The only known process capable of forming triplet states in carotenoids in the absence of a photo-sensitizer is singlet exciton fission. This process has been observed in several cases where intermolecular interactions occur between carotenoids 13,[17][18][19][20][21][START_REF] Wang | Proceedings of SPIE[END_REF][START_REF] Billsten | [END_REF] . In acenes, the mechanism of triplet formation is fairly well understood. A singlet-excited state (S1) on one pigment shares its energy with a neighboring pigment in the ground state (S0), resulting in a correlated pair of triplets 1 (TT) on these two molecules, entangled with spin coherence and overall zero spin. This strongly-entangled state has a singlet character, so that its decay to the singlet ground state is allowed with rates reminiscent of singlet excited states (ps-ns) 14,15 . Conversely, the intermediates involved in the evolution of this 1 (TT) state remain controversial. A weakly-entangled triplet (T1…T1) is necessary as an intermediate step, which may either yield two independent triplet states (T1 + T1), or recombine by triplet-triplet annihilation. The spectral position of the 1 (TT) signature is not well established -for acenes it appears blue-shifted with respect to an independent triplet, while it has been proposed as red-shifted for terrylenediimide dimers 49 and β-carotene aggregates 14,15 .

In carotenoids, photoexcitation leads to the formation of the second excited state (S2), with an energy of the lowest vibrational level at 2.57 eV. The (S1) singlet state consists of two coupled triplets residing on the same isolated molecule 10 . Consequently, the mechanism of singlet fission in daffodil chromoplasts must take into account the contribution of the second excited state, while also allowing for the excitonic character of S2, S1 and T1. Excitation energy is transferred from S2 to the higher vibrational levels of the first excited state (hot S1) in 150 femtoseconds, whereupon two pathways compete in ca. 500 fs: i) decay to the lowest vibrational states of S1, and ii) formation of an entangled triplet state 1 (TT) between two neighboring carotenoid molecules, with both triplets residing on the same excitonic manifold. The relaxed S1 state decays in 11 ps to the ground state (S0), whereas the entangled triplets 1 (TT) either decay in 36 ps to the ground state (transition allowed due to singlet character), or dissociate into weakly-entangled and then independent triplets (T1…T1 & T1 + T1, respectively). Target analysis based on this model, presented in Figure 5a, produces a remarkably good fit to the femtosecond data (see solid lines in Fig. 5b), and the resulting spectra exhibit the exact features expected from the model. This is clear from the hot S1, the relaxed S1 and the 1 (TT) states, which exhibit very well-separated ESA. Hot S1 has a broad and featureless induced absorption with a pronounced red shoulder, the relaxed S1 state exhibits a narrow and well-defined peak at 553 nm, and the 1 (TT) state peaks at 528 nm. This 1 (TT) decays following two competitive channels, with a 36-ps decay to the ground state, or towards (T1…T1) in 240 ps. The weakly-entangled triplets (T1…T1) and independent triplets (T1+T1) are modeled as two exponentially-decaying states with identical spectra. They show a broad induced absorption band peaking around 510 nm, with no ESA above 550 nm at long timescales. This set of experiments additionally a fairly precise estimation of the energy of the T1 state of lutein/violaxanthin aggregates, a parameter which has never been measured directly. The energy of the S2 state is 2.57 eV, as obtained from the absorption spectrum, and that of the relaxed S1 state has been estimated at 1.55 eV 50 . The energy level of 1 (TT) can be assumed to be equal to S1 within a difference depending on the entanglement energy (0.15eV) 14,15 , while the independent triplets must obey the energetic condition E(S1)  2E(T1). Hence the 1 (TT) energy is 1.550.20 eV, while the T1 energy level for aggregated carotenoids in daffodils is 0.780.10 eV (Scheme 1).

CONCLUSIONS

In daffodil chromoplasts, the interactions of lutein and violaxanthin lead to the existence of three different species of carotenoids aggregates, the most abundant being a mixture of lutein/violaxanthin with 0-0 transition at 481 nm. The photodynamics of this species induces the emergence of µs-lived triplet states produced via singlet fission, as already observed for self-associated lycopene in ripe tomato fruit. This first report of singlet fission occurring in mature flower plastids suggests that triplet fission is a phenomenon general to carotenoid aggregates in chromoplasts. While in tomatoes, lycopene interactions induce both a large redshift and strong distortions of the lycopene molecule 13 , in daffodil chromoplasts singlet fission is observed from lutein/violaxanthin aggregates which exhibit only a moderate red-shift in absorption transitions, and no molecular distortion of the interacting carotenoids. The fact that we observe singlet fission in both fruits and, for the first time here, flowers, occurring in carotenoid assemblies displaying different organizations and levels of interaction (strong and medium-to-weak, respectively), suggests a specific role of singlet fission in plastids. Finally, the structural organization of lutein/violaxanthin in daffodil chromoplasts produces well-defined electronic levels, allowing the energetic pathways following carotenoid singlet fission to be disentangled unequivocally for the first time.

ASSOCIATED CONTENT

Supporting Information. HPLC determination of carotenoid content in daffodil chromoplasts (Fig. S1); linear combination of in vitro absorption spectra of isolated carotenoids, to reproduce the spectrum of solubilized chromoplasts (Fig. S2). This material is available free of charge via the Internet at http://pubs.acs.org.

Figure 1 :

 1 Figure 1: Room temperature absorption spectra of (a) monomeric lutein and violaxanthin in n-hexane; (b) solubilized chromoplasts; and (c) intact daffodil chromoplasts in buffer solution. (d) 4.2 K absorption spectrum of daffodil chromoplasts, and (e) its second derivative. (f) CD spectrum of daffodil chromoplasts in buffer solution, obtained with 0.3 OD at 478 nm.

Figure 2 :

 2 Figure 2: 77 K resonance Raman spectra of violaxanthin (black) and lutein (grey) in cyclohexane, excited at 488.0 nm, and of intact chromoplasts excited at 406.7, 488.0 & 501.7 nm (blue, green, red respectively), showing the 930-1050 and 1100-1600 cm -1 regions (a, b respectively). The 3 and 1 regions are expanded in panels c and d, respectively (normalized to peak height).

Figure 3 :

 3 Figure 3: Femtosecond-to-nanosecond transient absorption data of daffodil chromoplasts obtained with excitation at 488 nm: (a) dataset overview; (b) selected absorption-difference kinetic traces (symbols -experimental data; solid lines -results of global fit, as per d); (c) time-gated spectra at selected delay times, as indicated; (d) EADS estimated by global fitting, using 6-component evolutionary model (decay times indicated in the legend). Greyed regions indicate spectral range corrupted by scattering of the pump pulses.

Figure 4 :

 4 Figure 4: Transient absorption in the ns-to-s window for daffodil chromoplasts upon 488-nm excitation: gated spectra and time evolution for 405 and 500 nm in (a) oxygen, (b) air, and (c) argon.

Figure 5 :

 5 Figure 5: Target model (a) and species-associated difference spectra (b) resulting from target analysis of femtosecond pump-probe data on daffodil chromoplasts. The estimated time constants are shown next to the corresponding arrows in panel a.

Scheme 1 :

 1 Scheme 1: Energetic pathways in aggregated lutein & violaxanthin after absorption of one photon (lifetimes as obtained from target analysis).

ACKNOWLEDGMENTS

This work was supported by the French Infrastructure for Integrated Structural Biology (FRISBI) ANR-10-INBS-05, and the CNRS as part of its 80PRIME interdisciplinary program. M.V. was supported by the European Regional Development Fund according to the supported activity 'Research Projects Implemented by World-class Researcher Groups' under Measure No. 01.2.2-LMT-K-718, grant number 01.2.2-LMT-K-718-01-0014. A.K.L. and the IJPB are supported by Saclay Plant Sciences-SPS (ANR-17-EUR-0007). This work benefited from the support of IJPB's Plant Observatory technological platforms.

AUTHOR INFORMATION

Corresponding Author * Mikas Vengris, Laser Research Center, Faculty of Physics, Vilnius University, Saulėtekio Ave. 10, LT-10223, Vilnius, Lithuania. e-mail: mikas.vengris@ff.vu.lt

Author Contributions

The manuscript was written through contributions of all authors.


Zhiwei Wang

Keywords: 2010 Mathematics Subject Classification. Primary 11N25 ; Secondary 11K65 largest prime factor, consecutive integers, positive proportion, sieve. L'auteur est partiellement soutenu

Let P + y (n) denote the largest prime factor p of n with p y. We prove that there exists a positive proportion of integers n such that P + y (n) < P + y (n + 1) for y = x α when α is small. Especially, the proportion is larger than 1/4 when α tends to 0, which improves our previous result.

Introduction

Ce travail est la continuation de notre article précédent [START_REF] Wang | Sur les plus grands facteurs premiers d'entiers consécutifs[END_REF] sur la fonction P + y (n), qui est définie comme le plus grand facteur premier d'entier n inférieur à y : P + y (n) := max{p | n : p y} avec la convention que P + y (n) = 1 si le plus petit facteur premier de n est strictement supérieur à y. Cette fonction est étudiée initialement par Rivat [START_REF] Rivat | On pseudo-random properties of P (n) and P (n + 1)[END_REF] pour approcher la conjecture d'Erdős et Pomerance [START_REF] Erdős | On the largest prime factors of n and n+1[END_REF] sur les plus grands facteurs premiers de deux entiers consécutifs :

(1. [START_REF] Balog | On triplets with descending largest prime factors[END_REF] |{n x : P + (n) < P + (n + 1)}| ∼ 1 2 x, où P + (n) désigne le plus grand facteur premier d'entier n avec la convention que P + (1) = 1. Rivat [START_REF] Rivat | On pseudo-random properties of P (n) and P (n + 1)[END_REF] En particulier, on a, dans le même domaine en (x, y),

(1.3) n x : P + y (n) < P + y (n + 1) ∼ 1 2 x Pour y = x, la formule (1.
3) est équivalente à la conjecture d'Erdős et Pomerance (1.1). Ainsi, il est intéressant d'étendre la domaine de y. Dans [START_REF] Wang | Sur les plus grands facteurs premiers d'entiers consécutifs[END_REF], nous démontrons le résultat suivant.

Théorème (B). Soit α ∈]0, 1]. Il existe C(α) > 0 tel que n x : P + y (n) < P + y (n + 1) {C(α) + o(1)}x (1.4) pour x → ∞ et y = x α .
Le Théorème (B) fournit des expressions explicites de C(α) admissibles. Cependant, on a lim α→0+ C(α) = 0, ce qui est contraire à l'intuition. Cela est dû au système de poids que nous utilisons qui est perfectible pour des petites valeurs de α. Dans cette note, nous corrigeons ce défaut par une méthode différente.

Théorème 1. Soit α ∈ ]0, 1 6 ]. Il existe une constante C 0 (α) > 0, définie en (3. [2,[START_REF] Erdős | On the largest prime factors of n and n+1[END_REF][START_REF] Teräväinen | On binary correlations of multiplicative functions[END_REF][START_REF] Wang | On the largest prime factors of consecutive integers in short intervals[END_REF][START_REF] Wang | Sur les plus grands facteurs premiers d'entiers consécutifs[END_REF]. De plus, pour les plus grands facteurs premier de trois entiers consécutifs, on peut voir [START_REF] Balog | On triplets with descending largest prime factors[END_REF][START_REF] Erdős | On the largest prime factors of n and n+1[END_REF][START_REF] Teräväinen | On binary correlations of multiplicative functions[END_REF][START_REF] Wang | Sur les plus grands facteurs premiers d'entiers consécutifs[END_REF].

Remerciements. Ce travail a été réalisé sous la direction de mes directeurs de thèse Cécile Dartyge et Jie Wu. Je les remercie vivement pour les nombreuses suggestions cruciales qu'ils ont proposées dans l'élaboration de ce travail.

Lemmes Fondamentaux

Nous rappelons tout d'abord un résultat d'Iwaniec sur le crible linéaire. Nous énonçons ici seulement la minoration car seule celle-ci sera utilisée dans cet article.

Soient A une suite finie d'entiers, P un ensemble de nombres premiers, z 2 un nombre réel, d un entier sans facteur carré dont les facteurs premiers appartiennent à P. Notons On a ainsi [START_REF] Iwaniec | Rosser's sieve[END_REF] Lemme 2.1. On suppose qu'il existe une constante L 2 telle que

u p<v 1 - w(p) p -1 log v log u 1 + L log u pour tout v > u 2. Alors pour D z 2, s = log D/ log z, on a S(A; P, z) > XV (z) f (s) + O 1 (log D) 1/3 - d<D, d|P (z)
|r(A, d)|, où f (s) est la fonction monotone croissante définie comme la solution continue du système de l'équation différentielle aux différences

F (s) = 2e γ /s, f (s) = 0, 0 < s 2, (sF (s)) = f (s -1), (sf (s)) = F (s -1), s 2. Remarque. f (s) vérifie f (s) = 0 pour 0 < s 2 et f (s) > 0 pour s > 2. En particulier, on a f (s) = 1 + O(e -s ) s → ∞.
(2.3) ϑ 0 (λ, u) := ρ(u/λ) + u 0 ρ(t/λ)ω(u -t)dt
avec la convention ϑ 0 (0, u) = 0. La fonction de Dickman ρ(u) est définie comme l'unique solution continue de l'équation différentielle aux différences

(2.4) ρ(u) = 1 si 0 u 1, uρ (u) = -ρ(u -1) si u > 1.
La fonction de Buchstab ω(u) est définie comme la solution continue du système

uω(u) = 1 si 1 u 2, (uω(u)) = ω(u -1) si u > 2.
De plus, nous prolongeons ω(u) par 0 pour u < 1.

Le troisième lemme [9, Proposition 1] est le théorème de type Bombieri-Vinogradov pour S(x; y, z). Lemme 2.3. Soit ε > 0. Pour tout A > 0, il existe une constante B = B(A) > 0 telle que

q x 1/2 /(log x) B max t x max (a, q)=1 n∈S(t; y,z) n≡a(mod q) 1 - 1 ϕ(q) n∈S(t; y,z) (n, q)=1 1 A,ε x (log x) A a lieu uniformement pour (2.5) 2 z y x et exp{(log x) 2/5+ε } y x,
où ϕ(q) est la fonction d'Euler.

Démonstration du Théorème 1

Pour α ∈ ]0, 1 6 ], soient K ∈ N * et β, δ ∈ R + * trois paramètres vérifiant les conditions : (3.1)

K (2α) -1 -2, 0 < β < α, ( 2α
) -1 -K δ (2α) -1 -α -1 βK. Posons y = x α , z = x β , ω(n; y, z) := p|n, z<p y 1.
Différente de [START_REF] Wang | On the largest prime factors of consecutive integers in short intervals[END_REF] et [START_REF] Wang | Sur les plus grands facteurs premiers d'entiers consécutifs[END_REF], notre point de départ est l'inégalité suivante :

(3.2)

n x P + y (n)<P + y (n+1)
1 k K n∈S(x; y,z), ω(n+1; y,z)=k

1 k K z<p 1 <•••<p k y n∈S(x; y,z) p 1 •••p k |(n+1) tel que (n+1)/(p 1 •••p k ) ∈S(x; y,z) 1 k K z<p 1 <•••<p k y S(A (b k ); P, y), où P := {p : p > z}, b k := p 1 • • • p k , A (b k ) := {(n + 1)/b k : n ∈ S(x; y, z) et n + 1 ≡ 0 (mod b k )}.
Remarque. Dans (3.2), on aurait pu utiliser la minoration

n x P + y (n)<P + y (n+1) 1 k K n∈S(x; y,z) ω(n+1; y,z)=k 1 + n∈S(x; y,z) ω(n+1; y,z) K+1
1 et minorer la deuxième somme du membre de droite de la même façon que dans [START_REF] Wang | Sur les plus grands facteurs premiers d'entiers consécutifs[END_REF] à l'aide des Lemme 2.2 et Lemme 2.3. Nous avons choisi de ne pas tenir compte de ce deuxième terme car sa contribution est très petite par rapport au premier terme lorsque α est petit.

Pour d | P (y, z), on a |A d (b k )| = |{(n + 1)/b k : n ∈ S(x; y, z), n + 1 ≡ 0 (mod b k d)}| = 1 ϕ(b k d) n∈S(x;y,z) (n,b k d)=1 1 + E(S(x; y, z); -1, b k d), où (3.3) E(S(x; y, z); -1, b k d) := n∈S(x;y,z) n≡-1(mod b k d) 1 - 1 ϕ(b k d) n∈S(x;y,z) (n,b k d)=1 1.
En remarquant que pour n ∈ S(x; y, z) on a automatiquement (n, b k d) = 1, on a donc

(3.4) |A d (b k )| = w(d) d X + r(A (b k ), d) avec X = Ψ 0 (x; y, z) ϕ(b k ) , w(d) = ϕ(b k )d ϕ(b k d) • 1 d , r(A (b k ), d) = E(S(x; y, z); -1, b k d).
En appliquant le Lemme 2.

n x P + y (n)<P + y (n+1) 1 S 1 -S 2 , où S 1 := {f (δ) + o(1)} k K z<p 1 <•••<p k y Ψ 0 (x; y, z) ϕ(b k ) p|P (y,z) 1 - w(p) p , S 2 := k K z<p 1 <•••<p k y d<D, d|P (y,z) |E(S(x; y, z); -1, b k d)|.
En appliquant le Lemme 2.3 et l'inégalité de Cauchy-Schwarz, on peut montrer

(3.6) S 2 K q x 1/2 (log x) -B τ K+1 (q)|E(S(x; y, z); -1, q)| K x (log x) 2 , où τ k (n) := n 1 •••n k =n 1.
Il reste à évaluer le terme principal S 1 . D'après la formule de Mertens, on a p|P (y,z)

1 - w(p) p = p|p 1 •••p k z<p y 1 - 1 p p p 1 •••p k z<p y 1 - 1 p -1 = 1 j k (p j -1) 2 p j (p j -2) z<p y 1 - 1 p 1 - 1 (p -1) 2 = ϕ(p 1 • • • p k ) p 1 • • • p k • log z log y 1 + O K 1 log z .
En reportant dans la définition de S 1 , le Lemme 2.2 et le théorème des nombres premiers nous permettent d'en déduire que (3.7)

S 1 = {f (δ) + o K (1)} k K z<p 1 <•••<p k y Ψ 0 (x; y, z) p 1 • • • p k • log z log y = {1 + o(1)}xf (δ)ϑ 0 β α , 1 α β α k K 1 k! log α β

Références

  montre que Théorème (A). Soient (a, b(b + 1)) = 1 et f y (n) := 1 si P + y (n + 1) > P + y (n), -1 si P + y (n + 1) < P + y (n). Alors l'inégalité 1 an+b x f y (an + b) a,b x exp -

A

  d := a ∈ A : d | a , P P (z) := p<z, p∈P p. On souhaite évaluer S(A; P, z) := |{a ∈ A : (a, P P (z)) = 1}|. On suppose que |A d | vérifie une formule de la forme |A d | = w(d) d X + r(A, d) pour d | P P (z), où X est une approximation de |A| indépendante de d, w une fonction multiplicative vérifiant 0 < w(p) < p pour p ∈ P, w(d)d -1 X un terme principal et r(A, d) un terme d'erreur que l'on espère petit en moyenne sur d. De plus, on définit V (z) := p<z, p∈P 1 -w(p) p .

  Pour z < y x, on définit P (y, z) := z<p y p, S(x; y, z) := {n x : (n, P (y, z)) = 1} (2.1) l'ensemble des entiers sans facteur premier dans l'intervalle ]z, y] et n'excédant pas x. On note son cardinal (2.2) Ψ 0 (x; y, z) := |S(x; y, z)|. Alors Ψ 0 (x; y, z) est évaluée par le lemme suivant (voir [6, Exercice 299]). Lemme 2.2. On a Ψ 0 (x; y, z) = ϑ 0 (λ, u)x{1 + O(1/ log z)} uniformément pour y z 2 et x yz, où u := (log x)/ log y, λ := (log z)/ log y et

  1 avec D = (x 1/2 (log x) -B )/b k , on obtient S(A (b k ); P, y) X p|P (y,z) |r(A (b k ), d)|.En reportant dans (3.2), il suit(3.5)

k.

  En reportant (3.6) et (3.7) dans (3.5), on obtient (1.5) avec C 0 (α) := max dv = e γ , et en prenant β = α/2, K = [1/(2α)] -2, δ = (2α) -1 -3 4 K, les formules (3.9), (3.10) et (3.8) nous permettent de déduire (1.6).

√log x -


Florian Richoux
email: florian.richoux@aist.go.jp

Jean-Franc ¸ois Baffier

Learning Interpretable Error Functions for Combinatorial Optimization Problem Modeling

In Constraint Programming, constraints are usually represented as predicates allowing or forbidding combinations of values. However, some algorithms exploit a finer representation: error functions. Their usage comes with a price though: it makes problem modeling significantly harder. Here, we propose a method to automatically learn an error function corresponding to a constraint, given a function deciding if assignments are valid or not. This is, to the best of our knowledge, the first attempt to automatically learn error functions for hard constraints. Our method uses a variant of neural networks we named Interpretable Compositional Networks, allowing us to get interpretable results, unlike regular artificial neural networks. Experiments on 5 different constraints show that our system can learn functions that scale to high dimensions, and can learn fairly good functions over incomplete spaces.

Introduction

Twenty years separate Freuder's papers [1997] and [2018], both about the grand challenges Constraint Programming (CP) must tackle "to be pioneer of a new usability science and to go on to engineering usability" [Freuder, 2007]. To respond to the lack of a "Model and Run" approach in CP [Puget, 2004;Wallace, 2003], several languages have been developed since the late 2000's, such as ESSENCE [Frisch et al., 2008], XCSP [Boussemart et al., 2016] or MiniZinc [Nethercote et al., 2007]. However, they require users to have deep expertise on global constraints and to know how well these constraints, and their associated mechanisms such as propagators, are suiting the solver. We are still far from the original Holy Grail of CP: "the user states the problem, the computer solves it" [Freuder, 1997].

This paper makes a contribution in automatic CP problem modeling. We focus on Error Function Satisfaction and Optimization Problems we defined in the next section. Compare to classical Constraint Satisfaction and Constrained Optimization Problems, they rely on a finer structure about the * Contact Author problem: the cost functions network, which is, for this work, an ordered structure over invalid assignments that constraint solver can exploit efficiently to improve the search.

In this paper, we propose a method to learn error functions automatically; a direction that, to the best of our knowledge, had not been explored in Constraint Programming. We focus here on "easy-to-use" aspects of Constraint Programming.

To sum up, the main contributions of this paper are: 1. to give the first formal definition of Error Function Satisfaction Problems and Error Function Optimization Problems, 2. to introduce Interpretable Compositional Networks, a variant of neural networks to get interpretable results, 3. to propose an architecture of Interpretable Compositional Network to learn error functions, and 4. to provide a proof of concept by learning Interpretable Compositional Network models of error functions, using a genetic algorithm, and to show that most of models give scalable functions, and remain fairly effective using incomplete training sets.

Error Function and Optimization Problems

Constraint Satisfaction Problem (CSP) and Constrained Optimization Problem (COP) are hard constraint-based problems defined upon a classical constraint network, where constraints can be seen as predicates allowing or forbidding some combinations of variable assignments.

Likewise, Error Function Satisfaction Problem (EFSP) and Error Function Optimization Problem (EFOP) are hard constraint-based problems defined upon a specific constraint network named cost function network [START_REF] Cooper | [END_REF]. Constraints are then represented by cost functions f :

D 1 × D 2 × . . . × D n → E,
where D i is the domain of i-th variable in the constraint scope, n the number of variables (i.e., the size of this scope) and E the set of possible costs.

A cost function network is a quadruplet V, D, F, S where V is a set of variables, D the set of domains for each variable, i.e., the sets of values each variable can take, F the set of cost functions and S a cost structure. A cost structure is also a quadruplet S = E, ⊕, ⊥, with E the totally ordered set of possible costs, ⊕ a commutative, associative, and monotone aggregation operator and ⊥ and the neutral and absorbing elements of ⊕, respectively.

In Constraint Programming, cost functions are often associated to soft constraints: they can be interpreted as preferences over valid or acceptable assignments. However, this is not necessarily the case: it depends on the cost structure. For instance, the classical cost structure S t/f = {true, f alse}, ∧, true, f alse makes the cost function network equivalent to a classical constraint network, so dealing with hard constraints.

Here, we consider particular cost functions that represent hard constraints only, by considering the additive cost structure S + = R, +, 0, ∞ . The additive cost structure produces useful cost function networks capturing problems such as Maximum Probability Explanation (MPE) in Bayesian networks and Maximum A Posteriori (MAP) problems in Markov random fields [Hurley et al., 2016].

In this paper, an error function is a cost function defined in a cost function network with the additive cost structure S + . Intuitively, error functions are preferences over invalid assignments. Let f c be an error function representing a constraint c and x c be an assignment of variables in the scope of c. Then f c ( x c ) = 0 iff x c satisfies the constraint c. For all invalid assignments i c , f c ( i c ) > 0 such that the closer f c ( i c ) is to 0, the closer i c is to satisfy c.

The goal of this paper is not to study the advantages of such cost function networks over regular constraint networks. Without formally defining EFSP and EFOP problems, some studies illustrate that solvers (in particular, metaheuristics) can exploit this structure efficiently leading to state-of-the-art experimental results, both in sequential [START_REF] Codognet | Yet another local search method for constraint solving[END_REF] and parallel solving [START_REF] Caniou | [END_REF]. In addition, our Experiment 3 shows that error functions representing the classic AllDifferent constraint gives models that clearly outperformed a model based on a regular constraint networks in terms of runtimes, for models with either hand-crafted or learned error functions.

Let x be a variable assignment, and denote by x c the projection of x over variables in the scope of a constraint c. We can now define the EFSP and EFOP problems. With the system we propose in this paper, users provide the usual constraint network V, D, C , and it computes the equivalent cost function networks V, D, F, S + . Learned error functions composing the set F are independent of the number of variables in constraints scope, and are expressed in an interpretable way: users can understand these functions and easily modify them at will. This way, users can have the power of EFSP and EFOP with the same modeling effort as for CSP and COP.

Related works

This work belongs to one of the three directions identified by Freuder [2007]: Automation, i.e., "automating efficient and effective modeling and solving". To the best of our knowledge, few efforts have been done on the modeling side.

Another of these three directions which is slightly related is Acquisition described by Freuder to be "acquiring a complete and correct representation of real problems". Remarkable efforts on this topic have been done by Bessiere's research team, for instance with constraints learning by induction from positive and negative examples [Bessiere et al., 2005] and with interactive queries asked to users [START_REF] Bessiere | [END_REF], and with constraint network learning also through with interactive queries [START_REF] Bessiere | [END_REF].

Model Seeker [START_REF] Beldiceanu | A model seeker: Extracting global constraint models from positive examples[END_REF]] is a passive learning system taking positive examples only, which are certainly easier for users to provide. It transforms examples into data adapted to the Global Constraint Catalog [START_REF] Beldiceanu | Nicolas Beldiceanu, Mats Carlsson, Sophie Demassey, and Thierry Petit. Global constraint catalog: Past, present and future[END_REF], then generate and simplify candidates by eliminating dominated ones. Model Seeker is particularly efficient to find a good inner structure of the target constraint network. Teso [2019] gives a good survey on Constraint Learning with this interesting remark: "A major bottleneck of [constraint-based problem modeling] is that obtaining a formal constraint theory is non-obvious: designing an appropriate, working constraint satisfaction or optimization problem requires both domain and modeling expertise. For this reason, in many cases a modeling expert is hired and has to interact with domain expert to acquire informal requirements and turn them into a valid constraint theory. This process can be expensive and time-consuming."

We can consider that Constraint Acquisition, or Constraint Learning, focuses on modeling expertise and puts domain expertise on background: users would not be able to understand and modify a learned model without the help of a modeling expert. The goal of these systems is mainly to simplify the interaction between the domain and the modeling experts.

Our work is taking the opposite direction: we focus on domain expertise and put modeling expertise on background. With our system, users always have the control over constraints' representation, which can be modified at will to fit needs related to their domain expertise. Constraint Implementation Learning is what best describes this research topic.

Method design

The main result of this paper is to propose a method to automatically learn an error function representing a constraint, to make easier the modeling of EFSP/EFOP. We are tackling a regression problem since the goal is to find a function that outputs a target value. Before diving into the description of our method, we need to introduce some essential notions.

Definitions

We propose a method to automatically learn an error function from the concept of a constraint. As described in [START_REF] Bessiere | Christian Bessiere, Frederic Koriche, Nadjib Lazaar, and Barry O'Sullivan[END_REF], the concept of a constraint is a Boolean function that, given an assignment x, outputs true if x satisfies the constraint, and false otherwise. Concepts are the predicate representation of constraints referred at the beginning of Section 2.

Our method learns error functions in a supervised fashion, searching for a function computing the Hamming cost of each assignment. The Hamming cost of an assignment x is the minimum number of variables in x to reassign to get a solution, i.e., a variable assignment satisfying the considered constraint. If x is a solution, then its Hamming cost is 0.

Given the number of variables of a constraint and their domain, the constraint assignment space is the set of couples ( x, b) where x is an assignment and b the Boolean output of the concept applied on x. Such constraint assignment spaces can be generated from concepts. These spaces are said to be complete if and only if they contain all possible assignments, i.e., all combinations of possible values of variables in the scope of the constraint. Otherwise, spaces are said to be incomplete.

In this work, we consider an error function to be a (nonlinear) combination of elementary operations. Complete spaces are intuitively good training sets since it is easy to compute the exact Hamming cost of their elements. We also consider assignments from incomplete spaces where their Hamming cost has been approximated regarding a subset of solutions in the constraint assignment space, in case the exact Hamming cost function is unknown.

Main result

To learn an error function as a non-linear combination of elementary operations, we propose a network inspired by Compositional Pattern-Producing Networks (CPPN). CPPNs [Stanley, 2007] are themselves a variant of artificial neural networks. While neurons in regular neural networks usually contain sigmoid-like functions only (such as ReLU, i.e. Rectified Linear Unit), CPPN's neurons can contain many other kinds of function: sigmoids, Gaussians, trigonometric functions, and linear functions among others. CPPNs are often used to generate 2D or 3D images by applying the function modeled by a CPPN giving each pixel individually as input, instead of considering all pixels at once. This simple trick allows the learned CPPN model to produce images of any resolution.

We propose our variant by taking these two principles from CPPN: having neurons containing one operation among many possible ones, and handling inputs in a sizeindependent fashion. Due to their interpretable nature, we named our variant Interpretable Compositional Networks (ICN). In this paper, our ICNs are composed of four layers, each of them having a specific purpose and themselves composed of neurons applying a unique operation each. All neurons from a layer are linked to all neurons from the next layer. The weight on each link is purely binary: its value is either 0 or 1. This restriction is crucial to obtain interpretable functions. A weight between neurons n 1 and n 2 with the value 1 means that the neuron n 2 from layer l + 1 takes as input the output of the neuron n 1 from layer l. Weight with the value 0 means that n 2 discards the output of n 1 .

Here is our method workflow in 4 points: 1. Users provide a regular constraint network V, D, C where C is a set of concepts representing constraints.

2. For each constraint concept c, we generate its ICN input space X, which is either a complete or incomplete constraint assignment space. Those input spaces are our training sets. If the space is complete, then the Hamming cost of each assignment can be pre-computed before learning our ICN model. Otherwise, the incomplete space is composed of randomly drawn assignments and only an approximation of their Hamming cost can be pre-computed.

3. We learn the weights of our ICN model in a supervised fashion, with the following loss function:

loss = x∈X (|ICN ( x) -Hamming( x)|) + R(ICN ) (1)
where X is the constraint assignment space, ICN ( x) the output of the ICN model giving x ∈ X as an input, Hamming( x) the pre-computed Hamming cost of x (only approximated if X is incomplete), and R(ICN) is a regularization between 0 and 0.9 to favor short ICNs, i.e., with as few elementary operations as possible, such that R(ICN ) = 0.9 × Number of selected elementary operations Maximal number of elementary operations . 4. We have hold-out test sets of assignments from larger dimensions to evaluate the quality of our learned error functions.

Notice we also have a hold-out validation set to fix the values of our hyperparameters, as described in Section 4.3.

Figure 1 is a schematic representation of our network. It takes as input an assignment of n variables, i.e., a vector of n integers. The first layer, called transformation layer, is composed of 18 transformation operations, each of them applied element-wise on each element of the input vector. Such operations are for instance the maximum between the i-th and i + 1-th elements of the input vector, or the number of j-th elements of the vector smaller than the i-th element such that j > i holds. This layer is composed of both linear and nonlinear operations. If an operation is selected (i.e., it has an outgoing weight equals to 1), it outputs a vector of n integers.

If k transformation operations are selected, then the next layer gets k vectors of n integers as input. This layer is the arithmetic layer. Its goal is to apply a simple arithmetic operation in a component-wise fashion on all i-th element of our k vectors to get one vector of n integers at the end, combining previous transformations into a unique vector. We have considered only 2 arithmetic operations so far: the addition and the multiplication.

The output of the arithmetic layer is given to the aggregation layer. This layer crunches the whole vector into a unique integer. At the moment, the aggregation layer is composed of 2 operations: Sum computing the sum of input values and Count >0 counting the number of input values strictly greater than 0.

Finally, the computed scalar is transmitted to the comparison layer with 9 operations. Examples of these operations are the identity, or the absolute value of the input minus a given parameter. This layer compares its input with an external parameter value, or the number of variables of the problem, or the domain size, among others.

All elementary operations in our model are generic: we do not choose them to fit one or several particular constraints. Due to the page limit, we cannot give a comprehensive list of the 18 transformation and 9 comparison operations here. Although an in-depth study of the elementary operations properties would be interesting, this is out of the scope of this paper: its goal is to show that learning interpretable error functions via a generic ICN is possible, and in the same way results with neural networks do not always use ReLU as an activation function, there is no reason to reduce ICN to its current 31 elementary operations or even a 4-layer architecture. Such elements can be changed by users to best fit their needs.

To have simple models of error functions, operations of the arithmetic, the aggregation, and the comparison layers are mutually exclusive, meaning that precisely one operation is selected for each of these layers. However, many operations from the transformation layer can be selected to compose the error function. Combined with the choice of having binary weights, it allows us to have a very comprehensible combination of elementary operations to model an error function, making it readable and intelligible by a human being. For instance, the most frequenly learned error function is Count >0 |{j | x[j] = x[i] and j > i}| for AllDifferent, and Euclidian p n i=1 x[i] for LinearSum, i.e., the Euclidian division of n i=1 x[i] -p by the maximal domain size, with a parameter p equals to the right hand side constant of the linear equation. Thus, once the model of an error function is learned, users have the choice to run the network in a feed-forward fashion to compute the error function, or to reimplement it directly in a programming language. Users can use our system to find error functions automatically, but they can also use it as a decision support system to find promising error functions that they may modify and adapt by hand.

Learning with Genetic Algorithms

Like any neural network, learning an error function through an ICN boils down to learning the value of its weights. Many of our elementary operations are discrete, therefore are not differentiable. Then, we cannot use a back-propagation algorithm to learn the ICN's weights. This is why we use a genetic algorithm for this task.

Since our weights are binary, we represent individuals of our genetic algorithm by a binary vector, each bit correspond-ing to one operation in the four layers indicating if the operation is selected to be part of the error function.

We randomly generate an initial population of 160 individuals, check and fix them if they do not satisfy the mutually exclusive constraint of the comparison layer. Then, we run the genetic algorithm to produce at most 800 generations before outputting its best individual according to our fitness function.

Our genetic algorithm is rather simple: The fitness function is the loss function of our supervised learning depicted by Equation 1. Selection is made by a tournament selection between 2 individuals. Variation is done by a onepoint crossover operation and a one-flip mutation operation, both crafted to always produce new individuals verifying the mutually exclusive constraint of the comparison layer. The crossover rate is fixed at 0.4, and exactly one bit is mutated for each selected individual with a mutation rate of 1. Replacement is done by an elitist merge, keeping 17% of the best individuals from the old generation into the new one, and a deterministic tournament truncates the new population to 160 individuals. The algorithm stops before reaching 800 generations if no improvements have been done in the last 50 generations. We use the framework EVOLVING OBJECTS [Keijzer et al., 2002] to code our genetic algorithm.

Our hyperparameters, i.e., the population size, the maximal number of generations, the number of steady generations before early stop, the crossover, mutation and replacement rates, and the size of tournaments have been chosen using ParamILS [Hutter et al., 2009], trained one week on one CPU over a large range of values for each hyperparameter. We use the same training instance used for Experiment 1, and new, larger instances as a hold-out validation set. These instances have been chosen because they are larger than our training instances and each of them contains about 4∼5% of solutions, which is significantly less than the 10∼20% of solutions in training instances.

Experiments

To show the versatility of our method, we tested it on five very different constraints: AllDifferent, Ordered, LinearSum, NoOverlap1D, and Minimum. According to XCSP specifications [Boussemart et al., 2016] 1 , those global constraints belong to four different families: Comparison (AllDifferent and Ordered), Counting/Summing (LinearSum), Packing/Scheduling (NoOverlap1D) and Connection (Minimum). Again according to XCSP specifications, these five constraints are among the twenty most popular and common constraints. We give a brief description of those five constraints below:

• AllDifferent ensures that variables must all be assigned to different values.

• LinearSum ensures that the equation x 1 + x 2 + . . . + x n = p holds, with the parameter p a given integer.

• Minimum ensures that the minimum value of an assignment verifies a given numerical condition. In this paper, we choose to consider that the minimum value must be greater than or equals to a given parameter p.

• NoOverlap1D is considering variables as tasks, starting from a certain time (their value) and each with a given length p (their parameter). The constraint ensures that no tasks are overlapping, i.e., for all indexes i, j ∈ {1, n} with n the number of variables, we have x i + p i ≤ x j or x j + p j ≤ x i . To have a simpler code, we have considered in our system that all tasks have the same length p.

• Ordered ensures that an assignment of n variables (x 1 , . . . , x n ) must be ordered, given a total order. In this paper, we choose the total order ≤. Thus, for all indexes i, j ∈ {1, n}, i < j implies x i ≤ x j .

Experimental protocols

We conducted three experiments, with two of them requiring samplings. These samplings have been done using Latin hypercube sampling to have a good diversity among drawn assignments. When we need to sample the same number k solutions and non-solutions, we draw assignments until we get k of solutions and k non-solutions. Due to stochastic learning, all learning and testing have been done 100 times. We did not re-run batches of experiments to keep the ones with the best results, as it should always be the case with such experimental protocols.

All experiments have been done on a computer with a Core i9 9900 CPU and 32 GB of RAM, running on Ubuntu 20.04. Programs have been compiled with GCC with the 03 optimization option. Our entire system, its C++ source code, experimental setups, and the results files are accessible on GitHub 2 .

Experiment 1: scaling The goal of this experiment is to show that learned error functions scale to high-dimensional constraints, indicating that learned error functions are independent of the size of the constraint scope.

For this experiment, error functions are learned upon a small, complete constraint assignment space, composed of about 500∼600 assignments and containing about 10∼20% of solutions. For each constraint, we run 100 error function learning over pre-computed complete constraint assignment space. Then, we compute the test error of these learned error functions over a sampled test set. Sampled test sets contain 10,000 solutions and 10,000 non-solutions, with 100 variables on domains of size 100, belonging to a constraint assignment space of size 100 100 = 10 200 , thus greatly larger than training spaces containing 500∼600 assignments.

We show normalized mean training and test errors: first, we compute the mean error among all assignments composing the training or the test set. Then, we divide it by the number of variables composing the assignments. Indeed, having a mean error of 5 on assignments with 100 variables and 10 variables is significantly different: the first one indicates a mean error every 20 variables, the second a mean error one in two variables.

2 https://github.com/richoux/LearningErrorFunctions/tree/1.1

Experiment 2: learning over incomplete spaces If, for any reasons, it is not possible to build a complete constraint assignment space, a robust system must be able to learn effective error functions upon large, incomplete spaces where the exact Hamming cost of their assignments is unknown.

In this experiment, we built pre-sampled training spaces by sampling 10,000 solutions and 10,000 non-solutions on large constraint assignment spaces of size between 10 12 and 10 13 , and with solution rates from 0.15% to 2.10 -7 %. Then, we approximate the Hamming cost of each non-solution by computing their Hamming distance with the closest solution among the 10,000 ones, and learn error functions on these 20,000 assignments and their estimated Hamming cost. Like for Experiment 1, we run 100 error functions learning of these pre-sampled incomplete spaces, so that each learning relies on the same training set. Finally, we evaluate the learned error functions over the same test sets than Experiment 1.

Experiment 3: learned error functions to solve problems The goal of this experiment is to assess that learned error function can effectively be used to solve toy problems. Here, we use a local search solver to solve Sudoku.

Sudoku is a puzzle game where all numbers in the same row, the same column and the same sub-square must be different. There, it can be modeled as a satisfaction problem using the AllDifferent constraint only. We run 100 resolutions of random 9 × 9 and 16 × 16 Sudoku grids, with a timeout of 10 seconds. If no solutions have been found within 10 seconds, we consider the run to be unsolved.

We consider the mean and median run-time to compare different representations of the AllDifferent constraint. We have two baselines: 1., a pure CSP model where constraints are predicates, and 2., an EFSP model with an efficient handcrafted error function representing AllDifferent. We compare those with two models using error functions learned with our system to represent AllDifferent: a., our EFSP model using the most frequently learned error function from the previous experiments and run through our neural network in a fastforward fashion, and b., our EFSP model with the same error function but directly hard-coded in C++. The solver and its parameters remain the same: the only thing that is modified in these four different models is the expression of the AllDifferent constraint.

Results

Experiments 1 & 2 Table 1 shows the training errors of Experiment 1, where error functions have been learned 100 times for each constraint. The first column contains the normalized mean training error of the most frequently learned error function among the 100 runs, with its frequency in parenthesis. Next columns concern the median, the mean and the standard deviation. Comparing Table 1 and the first half of Table 2 lead us to conclude that our system is able to learn error functions that scale for most constraint, namely AllDifferent, Linear-Sum and Minimum. Their median training errors in Table 1 are perfect of almost perfect, so as their median test errors on greatly larger constraint assignment spaces.

Results are not as good for NoOverlap1D and Ordered, which are clearly the most intrinsically combinatorial constraints among our five ones. One could think that our system is overfitting on its training set, but results from Experiment 2 lead us to another conclusion.

To see this, let's observe these Experiment 2's results by comparing the first and the second half of Table 2. Error functions learned over incomplete training spaces are as good as the ones learned over small complete spaces for LinearSum and Minimum. We observe significant improvements of the median and the mean for NoOverlap1D (36.07% and 55.76%) and Ordered (52.83% and 49.05%). This is due not because error functions from Experiment 1 were overfitting, but because spaces from Experiment 1 were too small for these highly combinatorial constraints, containing too few different combinations and Hamming cost patterns.

Experiment 3

The goal of this experiment is not to be state-of-the-art in terms of run-times for solving Sudoku, but to compare the average run-times of the same solver on four nearly identical Sudoku models presented in Section 5.1. For the model with a hand-crafted error function, we implemented the primal graph based violation error of AllDifferent from [START_REF] Petit | [END_REF]. This function simply outputs the number of couples with identical values within a given assignment. To run this experiment, we used the framework GHOST from Richoux et al. [2016], which includes a local search algorithm able to handle both CSP and EFSP models. Rows in gray means that some runs hit the 10-second timeout.

Table 3 shows that models with error functions clearly outperformed the model with the constraint represented as a predicate. Over 100 runs, no error function-based models hit the 10s timeout, but 4 runs of the regular constraint network model timed out on the 9 × 9 grid, and all of them on the 16 × 16 grid. Moreover, the learned error function hardcoded in C++ is nearly as efficient as the hand-crafted one (also coded in C++). The difference of runtimes between the learned error function hard-coded and computed through the ICN gives us an idea of the overhead of computing such a function through the ICN.

Discussions

Like Freuder [2007] wrote: "This research program is not easy because 'ease of use' is not a science." However, we believe our result is a step toward the 'ease of use' of Constraint Programming, and in particular about EFSP and EFOP. Our system is excellent for learning error functions of simple constraints over complete spaces. For intrinsically combinatorial constraints, learning over large, incomplete spaces should be favored. One of the most significant results in this paper is that our system outputs interpretable results, unlike regular artificial neural networks. Error functions output by our system are intelligible. This allows our system to have two operating modes: 1. a fully automatic system, where error functions are learned and called within our system, being completely transparent to users who only need to furnish a concept function for each constraint, in addition to the regular sets of variables V and domains D, and 2. a decision support system, where users can look at a set of proposed error functions, pick up and modify the one they prefer.

The current limitation of our system is that it struggles to learn high-quality error function for very combinatorial constraints, such as Ordered and, in particular, NoOverlap1D. By combining results from Experiments 1 and 2, we can conclude that: 1. our system is not overfitting but need more diverse and expressive operations to learn a high-quality error function for such constraints, and 2. the Hamming cost is certainly not the better choice to represent their assignment error.

An extension of our work would be to do reinforcement learning rather than supervision learning based on the Hamming cost. Learning via reinforcement learning would allow finding error functions that are more adapted to the chosen solver.

  Problem: ERROR FUNCTION SATISFACTION PROBLEM Input: A cost function network V, D, F, S + . Question: Does a variable assignment x exist such that ∀f c ∈ F, f c ( x c ) = 0 holds? Problem: ERROR FUNCTION OPTIMIZATION PROBLEM Input: A cost function network V, D, F, S + and an objective function o. Question: Find a variable assignment x maximizing or minimizing the value of o( x) such that ∀f c ∈ F, f c ( x c ) = 0 holds.

Figure 1 :

 1 Figure 1: Our 4-layer network. Layers with blue neurons have mutually exclusive operations.

Table 2 :

 2 Table 2 contains the normalized mean test errors of error functions learned with Experiments 1 and 2, with their median, mean and standard deviation. The normalized mean test error of the most frequently learned error function for each constraint in each experiment has been isolated in the first column of number, for comparison. Test errors in high dimensions of error functions learned with Experiments 1 and 2.

	Constraints most (freq.) median mean std dev.
	AllDifferent 0	(98) 0		5.001 35.185
	LinearSum 0.004 (48) 0.059	0.032	0.027
	Minimum	0	(71) 0		0.026	0.044
	NoOverlap 0.039 (32) 0.074	0.074	0.030
	Ordered	0.020 (100) 0.020	0.020	0
	Table 1: Training errors (100 runs) of Experiment 1, over small and
	complete constraint assignment spaces.
	Exp.	Constraints	most freq median	mean std dev
		AllDifferent 0	0	0.017 0.119
		LinearSum	3×10 -4	0.019	0.179 0.341
	1	Minimum	0	0	1.435 4.866
		NoOverlap	0.268	0.316	0.486 0.682
		Ordered	0.106	0.106	0.106 0
		AllDifferent 0.052	0.052	0.052 0
		LinearSum	3×10 -4	3×10 -4 0.200 0.629
	2	Minimum	0	0	0.193 0.978
		NoOverlap	0.202	0.202	0.215 0.020
		Ordered	0.050	0.050	0.054 0.008

Table 3 :

 3 Mean run-times in milliseconds over 100 runs to solve Sudoku with 4 different representations of the AllDifferent constraint.

	Sudoku	Error Function	mean median	std dev
		nothing (CSP) 624.41	217.21 1,1196.80
	9 × 9	hand-crafted fast-forward	33.86 55.57	32.00 49.27	10.34 41.17
		hard-coded	33.48	31.76	9.48
		nothing (CSP)	-	-	-
	16 × 16	hand-crafted fast-forward	432.34 825.61	393.95 774.36	164.06 271.09
		hard-coded	537.48	539.86	162.03

see also http://xcsp.org/specifications