Marie-Hélène Morel
email: marie-helene.morel@inrae.fr

Justine Pincemaille

Laure Lecacheux

Paul Menut

Laurence Ramos

Amélie Banc

Thermodynamic insights on the liquid-liquid fractionation of gluten proteins in aqueous ethanol

Keywords: Liquid-liquid phase separation, Intrinsically disordered protein (IDP), gluten, glutenin polymers, Flory-Huggins interaction parameter

Highlights (85 characters including spaces)

 An ethanol/water extract containing 30% glutenin and 60% gliadin was obtained from gluten  The extract (45 g/L) underwent liquid-liquid phase separation (LLPS) below 14°C  Glutenin polymers were fractionated by LLPS according to their molecular weight  LLPS delivered fractions with large range in glutenin/gliadin weight ratios (up to 2.5)  Flory-Huggins lattice model properly reflected gluten proteins partitioning by LLPS Abstract (250 words) Wheat gluten includes two major proteins classes, gliadin (25-65 kg/moles) and glutenin polymers (100 to >2,000 kg/mol) each comprising several polypeptides routinely identified by size-exclusion chromatography and electrophoresis. Gluten proteins are rich in glutamine (30%) and contain several repeated sequences, linking them to the wide class of intrinsically disordered protein (IDP). Here we showed that an ethanol/water (EtOH/W, 50/50, v/v) extract of an industrial gluten comprising 1/3 of glutenin polymers and 2/3 of gliadin, underwent liquid-liquid phase separation (LLPS) below 14°C, leading to two coexisting phases, respectively rich and poor in protein. As the quenching depth increased, proteins of lower and lower molecular weight joined the rich phase, akin to what would have been obtained for a polydisperse polymer sample. Within the rich phase the mass ratio of glutenin over gliadin shifted from 2.5 to 0.5 when the temperature decreased from 14°C to -0.8°C. Concomitantly the concentration in glutenin polymers increased until 9°C and then stopped to evolve, suggesting that the LLPS line intersected the gelation line below 9°C. We showed that the partitioning of the glutenin polymers in the coexisting phases obeyed the Flory-Huggins (FH) lattice model, with different partial specific volumes in the coexisting phases and that the enthalpy of mixing of all gluten proteins with ethanol/water was size-independent. However, between 14°C and 5°C, some gliadins joined the rich phase, whereas on a strict thermodynamic basis and owing to their small-size they should not. The feature exemplified their specific interactions with glutenin polymers through weak forces.

compositions, being rich in glutamine (35-40%) and proline (12-27%) and poor in charged residues (5-6.5%), the majority of which being basic (basic/acidic ratio, 1.2-2.5). As another specific feature, gluten polypeptides comprise several repeated sequences of 6 to 8 amino-acids, made of glutamine, proline and phenylalanine (LMW-GS and gliadins) or glutamine, proline, glycine and tyrosine (HMW-GS) [START_REF] Wieser | Chemistry of gluten proteins[END_REF]. On the grounds of these features gluten proteins belongs to the class of intrinsically disordered proteins and should display extensive structural plasticity [START_REF] Rauscher | Proline and glycine control protein self-organization into elastomeric or amyloid fibrils[END_REF].

The way the gliadin and glutenin interact in the gluten viscoelastic mass is still poorly understood. It is likely that subtle balances of hydrogen, hydrophobic and ionic forces are at play, but this remains experimentally difficult to establish. At the origin of the problem is first the complexity of gluten composition, with several dozens of different components. Secondly, the lack of solubility of gluten protein in water has limited the use of classical scattering methods, especially because establishing the 3D structure of a protein dispersed in a denaturing solvent (to maintain solubility) is unlikely to provide relevant information about functionality. Assimilating gluten proteins to intrinsically disordered proteins (IDPs) opened up a new world of thought by challenging the classical structurefunction paradigm. It allows one to choose the solvent solely on the basis of thermodynamic grounds. Thus, while 50% ethanol in water induces denaturation and aggregation of most globular proteins, it can be assimilated to a good to theta solvent of gluten proteins, been known for a long time to be well suited to maintain them into solution [START_REF] Beckwith | Isolation and characterization of a high-molecular-weight protein from wheat gliadin[END_REF][START_REF] Jackson | Characterisation of high molecular weight gliadin and low-molecular-weight glutenin subunits of wheat endosperm by two-dimensional electrophoresis and the chromosomal localisation of their controlling genes[END_REF][START_REF] Huebner | Improved chromatographic separation and characterization of ethanol-soluble wheat proteins[END_REF]. In addition, a century ago, [START_REF] Dill | Preparation, solubility, and specific rotation of wheat gliadin[END_REF] showed that aqueous ethanol (55% to 70%, v/v) extracts of wheat gluten obtained at ambient temperature underwent liquid-liquid phase separation (LLPS) when temperature was lowered; the phenomena being fully reversible and repeatable at will (ibid; [START_REF] Dill | The behavior of the prolamins in mixed solvents[END_REF]. The solubility of polymers in a solvent is driven by thermodynamic considerations and relies on null to negative change of Gibbs free energy of mixing (G mix = H mix -TS mix ). In the case of moderately bad or theta solvent, the polymer transfer into a solvent phase is endothermic, and the solubility entirely relies on entropic gain. Mixing is always accompanied by a gain in combinatorial entropy (S mix >0), whose amplitude is inversely proportional to the polymer to solvent molar volume ratio. Thus, irrespective of chemical composition, large polymers are always less soluble than their small counterparts. They are the first to undergo phase separation when the temperature decreases, causing the Gibbs free energy of mixing to become positive except for the two coexisting phases where solvent-polymer or polymer-polymer interactions are privileged.

Recently, some of us [START_REF] Boire | Phase behaviour of a wheat protein isolate[END_REF] took advantage of the LLPS of a wheat flour protein extract (in 0.5 mM NaCl in 55% ethanol) to recover a soluble fraction including mainly gliadin (less than 10% of glutenin polymers) following a quench at 2°C,. They studied the LLPS behavior of this fraction and found that the glutenin polymers partitioning between the two coexisting phases was consistent with the Flory-Huggins lattice model developed for polymer solutions. In addition, some of us investigated the potential of LLPS to deliver protein samples of given glutenin/gliadin ratio, starting from an ethanol/water (50/50, v/v) wheat gluten extract. Using the same binary solvent and in the dilute regime, multi-angle light scattering analysis revealed that a balanced blend of gliadin/glutenin polymers displayed the typical features of Gaussian chains in a good solvent [START_REF] Dahesh | Polymeric assembly of gluten proteins in an aqueous ethanol solvent[END_REF]. In the semi-diluted to the concentrated regime (100-400 g/L) and by using a combination of various scattering techniques, large protein assemblies of radius of gyration ∼120 nm were identified (ibid). More recently we showed that these protein assemblies comprised glutenin polymers and specific ω-gliadins, interacting together through weak bonds [START_REF] Morel | Insight into gluten structure in a mild chaotropic solvent by asymmetrical flow fieldflow fractionation (AF4) and evidence of non-covalent assemblies between glutenin and ωgliadin[END_REF]. The dynamic of the phase separation, the rheological and gelling properties of these gliadin/glutenin samples have been recently investigated [START_REF] Dahesh | Spontaneous gelation of wheat gluten proteins in a food grade solvent[END_REF][START_REF] Banc | Small angle neutron scattering contrast variation reveals heterogeneities of interactions in protein gels[END_REF][START_REF] Banc | Phase separation dynamics of gluten protein mixtures[END_REF][START_REF] Ramos | Impact of the protein composition on the structure and viscoelasticity of polymer-like gluten gels[END_REF]). Here we present in details the main experimental parameters governing the extraction yield of gluten protein in ethanol/water and driving their fractionation by LLPS, highlighting the underlying thermodynamic causes.

2 Materials and methods 2.1 Chemicals Deionized water and ethanol 99.9% of analytical grade (CarloErba) were used to prepare the water/ethanol solvent solutions. Other chemicals were of ACS grade. Industrial gluten powder batches were courtesy of TEREOS-SYRAL (Aalst, Belgium). Gluten protein content was estimated by Kjeldahl analysis (AACC-Method 46-10). Most experiments, unless otherwise specified, were performed using a gluten batch including 5.89% water (wb) and 74.5% protein (wb) and comprising 42.5%, 43.8% and 13.6% glutenin polymers, gliadin and albumin/globulin, respectively.

2.2 Extraction and liquid-liquid fractionation of gluten protein in 50/50 (v/v) ethanol/water Typically, gluten (20 g) was placed in a 250 mL propylene Nalgene bottle, in which 200 mL of water/ethanol (50/50, v/v) solvent was added. The bottle was closed and vigorously shaken by hand for 60 s to disperse gluten. The bottle was then agitated for 19 h at 20°C, using a rotary shaker (Reax 2, Heidolph) rotating at 60 rpm. After 30 min of centrifugation (15 000 g, 20°C) the supernatant (parent extract, Par) was recovered, weighted and placed in a clean bottle, while the pellet (gluten gel, Gel) was sampled for high performance size-exclusion chromatography (HPSEC) analysis and then discarded. The bottle containing the parent extract was immersed in a stirred bath regulated in temperature (from for 12.5°C to -0.8°C) to trigger phase separation. The temperature quenching rate was about -0.5 °C/min. After one hour in the thermostated bath maintained at a quench temperature Tq, the bottle was centrifuged (3 000 g) for one hour at Tq. Then the upper (poor phase, P phase) and bottom (rich phase, R phase) phases were collected and weighted. Routinely 6 replicate bottles were run in parallel.

High Performance Size-Exclusion Chromatography (HPSEC)

HPSEC analysis of the proteins present in the gluten powder and in the collected Gel, Par, P et R phases was performed as described by Morel et al. (2000). Gluten powder (20 mg at 1g/mL) was first extracted at 60°C in sodium phosphate 100 mM, pH 6.9; SDS 1% (sample buffer below) by rotary shaking (Reax 2, Heidolph at 60 rpm for 80 min). The supernatant was recovered by centrifugation (30 min, 18 000rpm) for analysis of SDS-soluble proteins. The pellet was further extracted using 5 mL of the sample buffer applying a sonication step (Vibra Cell sonificator at 20 kHz for 90s) allowing one to recover the SDS-insoluble glutenin polymers into solution through mechanical disruption. The gluten gel (Gel) and the ethanol/water parent extract (Par), poor (P) and rich (R) phases (cf. 2.2) were simply diluted 201, 31, 41 and 201 times in the sample buffer supplemented with 6M urea before analysis. Routinely, 20µL was injected onto a size exclusion column TSK-G 4000-SW (Merck, France) (7.5 mm x 30 cm) coupled to a TSK 3000-SW (Merck, France) guard column (7.5 mm, 7.5 cm). Elution was performed with 0.1 M sodium phosphate buffer, 0.1% SDS and protein detection was monitored at 214 nm. The HPSEC UV trace was converted into protein mass concentration using a specific extinction coefficient of 18.51 L/cm/g (Morel & Bar-L'Helgouac'h, 2000). From the protein concentration, the volume of each samples (Gel, Par, P and R) was calculated considering the weighted mass and the densities of gluten protein and solvent (respectively, 1.32 and 0.923). Then the total mass of protein in each extract was calculated from the protein concentration deduced from HPSEC and the calculated sample volume. The apparent molecular weight of the eluted species was assessed after column calibration with known protein standards saturated with SDS, as described in [START_REF] Redl | Extrusion of wheat gluten plasticized with glycerol: Influence of process conditions on flow behavior, rheological properties, and molecular size distribution[END_REF].

SDS-PAGE

Freeze dried samples were suspended at 3 mg/mL in a loading buffer including 0.5 M Tris-glycine buffer (pH 6.8), 2% SDS (w/v), 5% beta-mercaptoethanol (v/v), 20% glycerol (v/v) and traces of pyronin Y. After 2 h shaking at ambient temperature, sample tubes were heated for 2.5 min in boiling water. The discontinuous SDS-PAGE buffer system of Laemli was used with a running gel (T=15%, C=0.5%) and a stacking gel (T=5%, C= 10%). Samples loading was 10 µL per well (comb of 10 wells). A homemade vertical gel (160 x 180 x 1.5 mm) was run at 20 mA/gel at 18°C and stopped 1 h after the tracking dye (pyronin Y) had reached the bottom of the gel. The gels were stained in 12.5% (w/v) trichloroacetic acid, 0.01% (w/v) Coomassie Brilliant Blue R250, and destained with distilled water before being scanned using a Biorad GS-700 gel scanner. Densitometry profiles of each sample track were extracted from inverted scanned images, using FiJi software and a macro computing the sum of pixels' grey level from a box of one pixel high and of the track width.

Results & discussion

3.1 Impact of solid/solvent ratio and ethanol content on gluten protein extraction A gluten powder including 84.5% protein (wb) of which 89% were soluble in SDS-phosphate buffer, one of the best extracting solvent for gluten protein, was used for studying the impact of solid/solvent ratio on the protein extraction yield. Extraction was performed for 19 h at 20°C, using aqueous ethanol (50/50, v/v) for solid/solvent ratios ranging from 10 g to 20 g/100 mL, in order to recover the soluble proteins fraction (thereafter parent fraction, Par) according to figure 1a. The percentage of gluten proteins recovered in Par averaged around 44±3% (Table 1). The protein concentration in Par increased proportionally to the solid/solvent ratio, while the volume ratio of Par over the remaining insoluble gel fraction (thereafter gel fraction, Gel) tended to decrease (Table 1). Attempts to further increase the protein concentration in Par by rising the solid/solvent ratio beyond 20 g/ 100 mL was not feasible because more and larger clumps formed, which contrarily to what was obtained at lower ratios, were impossible to disperse, even after prolonged shaking. HPSEC analysis of the proteins present in the Par and Gel fractions revealed that even if Gel included a majority of glutenin polymers it also comprised one third of gliadins, whereas Par included one third of glutenin polymers (Table I, Figure 1b). As already noted by several authors [START_REF] Beckwith | Isolation and characterization of a high-molecular-weight protein from wheat gliadin[END_REF][START_REF] Jackson | Characterisation of high molecular weight gliadin and low-molecular-weight glutenin subunits of wheat endosperm by two-dimensional electrophoresis and the chromosomal localisation of their controlling genes[END_REF] a clear cut separation of gliadin and glutenin polymers could thus not be obtained by using aqueous ethanol for gluten protein extraction.

Whatever the solid/solvent ratio, the Gel was enriched in protein (concentration 236± 16 g/L), while its volume ranged from 14 to 33% of the total sample volume, as the solid/solvent ratio doubled. Typically, the Gel consisted in a pale greyish cohesive mass of elastic consistency (Figure 1c). This elastic mass was highly solvated (solvent volume fraction φ = 0.8±0.1). Interestingly, it was significantly more solvated than typical gluten gel extracted from flour doughs (φ = 0.66 to 0.74, depending on the flour baking quality) [START_REF] Mccaig | Changes in the physical properties of gluten with aging of flour[END_REF]. This higher swelling ratio could be due to the better solvent quality of aqueous ethanol for gluten proteins as compared to pure water. The Gel displayed stable protein composition and concentration whatever the solid/solvent ratio. This may not be surprising, considering that hydrated gluten is a critical gel [START_REF] Ng | Power law gels at finite strains: The nonlinear rheology of gluten gels[END_REF]. Replacement of water by aqueous ethanol, a better solvent for gluten proteins simply shifted the sol-gel equilibrium of gluten, increasing the gel swelling capacity, while promoting gliadin and glutenin solubility. Because our ultimate aim was to fractionate the gluten proteins according to their molecular weight by liquid-liquid phase separation (LLPS), in the following we selected the lowest solid/solvent ratio (10 g/100 mL), as it gave the lowest supernatant concentration. Indeed, the partitioning efficiency between small and large polymers is better when starting from very diluted polymer solutions [START_REF] Flory | Thermodynamics of heterogeneous polymers and their solutions[END_REF]. According to the seminal work of [START_REF] Dill | Preparation, solubility, and specific rotation of wheat gliadin[END_REF] the cloud point temperature of gliadin suspended in ethanol/water mixtures is highly dependent on ethanol content. Isobaric cloud point temperature curves of gliadin follow a U-profile, nearly flat between 55 and 65% (v/v) ethanol where the solvent quality is the best possible one. This prompted us to examine the impact of the ethanol content in the solvent. As shown in figure 2, an increase of the ethanol/water ratio from 48/52 to 66/34 (v/v) slightly decreased the protein extraction yield. It also shifted the protein composition at the benefit of gliadin leading to a drop in the glutenin/gliadin ratio from 0.5 to less than 0.4 when the ethanol content increased. Since our objective was to obtain by LLPS samples especially enriched in glutenin polymers, the soluble protein extract obtained using 50/50 ethanol/water (v/v) was adopted in the following.

Fig. 2: Impact of solvent composition on gluten protein extraction yield and supernatant composition.

Extraction was performed using 10% of gluten in solvents of various ethanol contents (v/v), according to 2.2.

Impact of the quench depth on gluten fractionation by LLPS

The experiments reported hereafter were performed using the parent fraction obtained as described in 2.2. On average, 49 ± 3 % of the gluten proteins were recovered in the Parent (96 extractions), at a concentration of 46 ± 2 mg/ml (with 57 ± 6% gliadin, 31 ± 3% glutenin polymers and 11 ± 1.4 % albumin/globulin).

The Parent remained stable at ambient temperature for days but became cloudy when cooled below the cloud temperature, T cloud = 14°C, while subsequent sample warming reversed the phenomena. Phase Below T cloud , Parent demixing into well delineated bottom and upper phases (thereafter called rich and poor phases) could take hours, especially when the quench temperature (Tq) was closed to T cloud . To recover the demixed phases, samples, once tempered for one hour at a given quench temperature (Tq with Tq < T cloud ), were submitted to mild centrifugation (1 hour at 3 000g) at Tq. We verified that extended tempering (from 1 to 4 hours) or centrifugation did not change the respective volumes or compositions of the recovered phases (see Fig. 5ab, below). We plot in figure 3a the protein concentration of the rich and poor phases as a function of the quench depth ΔT=T cloud -Tq. Figures 3b and3c report the partitioning of the parent volume and protein content between the demixed rich and poor phases as a function of ΔT. Below a quench depth ΔT=5°C the protein concentration in the poor phase remained almost constant (41±2 mg/ml), close to that of the Parent one (46±2 mg/ml). Concomitantly, the protein concentration in the rich phase increased substantially from 84±7 to 173±27 mg/mL because of the large imbalance in the coexisting phases volumes (Fig. 3b). At a quench depth of 1.5°C, the protein concentration in the rich phase was already nearly twice as much as in the Parent phase (84 vs 46 mg/mL, Fig. 3a). Ultimately, at the highest quench depth tested (ΔT=14.8°C), about 80% of the Parent proteins joined the rich phase (Fig. 3c) whose protein concentration reached 376±32mg/mL (Fig. 3a). Standard errors from 6 replicate bottles.

3.3 Quantitative partitioning of gluten protein types depending on the quench depth Fig. 4 shows the HPSEC profiles of proteins from parent (yellow) and demixed phases (blue and orange for poor and rich phases, respectively) obtained at different quench temperatures. Profiles were expressed on a weight basis to reflect the partitioning of the parent proteins into rich and poor phases. Typically, looking at the parent phase (yellow profile), glutenin polymers appeared first for elution times between 8.0 to 13.2 min. The first fraction between 8.0 and 8.9 min, demarcated by vertical bars on Fig. 4, coincides with the column void volume and includes the glutenin polymers of the largest size (LS-glu polymers). The molecular weight of the LS-glu polymers could not be estimated precisely but should be comprised between 700to more than 2,000 kg/mol. Next, glutenin polymers ranging from 700 to 80 kg/mol were recovered (medium size, MS-glu polymers). The following small shoulder (13.2-14.2 min) corresponded to ω-gliadin with molecular mass of about 60 kg/mol. The high peak between 14.2 and 16.1 min included other types of gliadin (γ, β and α; from 45 to 18 kg/mol). Lastly albumin and globulin were eluted (16.1-17.80 min). Just below T cloud , following a quench depth ΔT=1.5 °C, the rich phase included almost exclusively LS-glu polymers. At ΔT=5°C, all LS-glu polymers from Parent joined the rich phase, as exemplified by the overlapped HPSEC profiles of Parent and rich phases, at the location of the exclusion peak. For deeper quenches, in addition to LS-glu polymers, medium-sized glutenin polymers joined the rich phase. Gliadins started to join the rich phase beyond a quench depth of 5°C, as revealed by the down shift of the main peak of the blue profile (poor phase) compared to the Parent one. At ΔT=14.8°C the rich phase included almost all proteins, the poor phase being diluted (protein concentration 7.4 mg/mL) and almost exclusively composed of γ-, α/β-gliadins (92%) supplemented by albumin/globulin. HPSEC analysis revealed that as the quench depth increased, gluten proteins tended to join the rich phase in sequence according to their molecular weights, from the largest to the smallest species. As a result, the composition of the rich phases progressively shifted towards decreasing proportion in glutenin (Fig. 5a andb). Thus, the liquid-liquid phase separation of an ethanol/water extract of gluten proteins allowed one to obtain protein extracts with a tunable proportion of glutenin from 0.6 to 65%. Note that ultimately extracts richer in glutenin polymers than the Gel fraction could be obtained (Table 1).

Fig. 5ab. Proportion in glutenin polymers (in % of total protein mass) and glutenin/gliadin mass ratio in the rich and poor phases recovered at different quench depths.

Filled and empty black symbols are from the present study, while red squares are from a previous study of our group [START_REF] Dahesh | Polymeric assembly of gluten proteins in an aqueous ethanol solvent[END_REF][START_REF] Dahesh | Spontaneous gelation of wheat gluten proteins in a food grade solvent[END_REF] obtained using another batch of gluten and letting the phase separation to proceed naturally at ΔT=10°C during 24H without centrifugation.

Polypeptides composition of parent and demixed phases

Figure 6 shows a picture of the SDS-PAGE obtained for gluten, Parent and some of the poor and rich phases recovered after applying quench depths, ΔT = 5°, 8° and 12°C. On the right hand side of the picture, the densitometry profile of the gluten protein sample is presented. High-molecular weight subunits (HMW-GS) and -gliadins were easily identified in the upper part of the gel while the patterns of the other gliadins and of the low-molecular weight glutenin subunits (LMW-GS) overlapped, preventing a specific band assignment. At the bottom of the gel, a protein doublet was identified as belonging to albumin/globulin proteins. Comparison of the rich and poor patterns showed that the 4 bands corresponding to the HMW-GS behaved as one, joining together the rich phase. The situation was more contrasted when considering the ω-gliadins pattern. The two upper bands (of higher molecular weight) disappeared from the poor phase at ΔT=8°C, while the two other seemed rather unaffected. Despite, their close molecular weight, ω-gliadins demonstrated different ability of to undergo phase separation.

The contrast between the different patterns were quantified by gel densitometry (Table 2). In the coexisting phases, the higher the relative proportion in LMW-GS/-gliadin, the lower was the proportion in HMW-GS. The -gliadin did not follow the same trend as the LMW-GS/-gliadin since it was preferentially enriched in the rich phases, like the HMW-GS.

As anticipated, a linear relationship (R 2 =0.996) was obtained between the content in glutenin polymers of the samples, determined from HPSEC analysis, and the proportion of the HMW-GS measured by densitometry (glutenin polymer a and HMW-GS*, Table2). Considering the total glutenin polymers content determined from the HPSEC profiles of the samples (Fig. 4 and5a) and assuming that it coincided to the sum of the HMW-and LMW-glutenin subunits

LMW-GS -gliadins

contents determined from gel densitometry, one could estimate the LMW-GS content. This assumption neglected the possible involvement of atypical gliadins within the glutenin polymer chain architecture. Comprising an odd number of cysteine residues, the atypical gliadins were found covalently linked to glutenin polymers through disulfide bonds [START_REF] Vensel | Protein composition of wheat gluten polymer fractions determined by quantitative two-dimensional gel electrophoresis and tandem mass spectrometry[END_REF]. They may act as terminator of the glutenin polymer chain propagation and therefore might be overrepresented in the small glutenin polymers, as compared to larger ones. The HMW-GS/LMW-GS glutenin subunits ratio for gluten, parent and rich phases were of the same order of magnitude, around 0.5, well within the typical range of wheat flours (HMW-GS/LMW-GS ratio about 0.46±0.05) as rigorously assessed by reverse-phase HPLC analysis [START_REF] Wieser | Correlations of the amount of gluten protein types to the technological properties of wheat flours determined on a micro-scale[END_REF]. Hence, even if the rich phases differed in their glutenin polymers contents (Fig. 5a) and size distributions (Fig. 4), they would include polymers of relatively identical subunits composition. By contrast the poor phases obtained at ΔT=8°C and in particular at ΔT =12°C, showed low to very low HMW-GS/LMW-GS ratios (0.34 and 0.05, respectively). The small glutenin polymers remaining in these two poor phases (see blue profiles in Fig. 4) would be enriched in gliadin terminators.

Coexisting concentration curves

Large supra-molecular protein assemblies (~ 120 nm) displaying an apparent molecular size of 26,000 kg/mol were previously identified in the rich phase recovered at ΔT=10°C [START_REF] Dahesh | Polymeric assembly of gluten proteins in an aqueous ethanol solvent[END_REF]. These assemblies were successfully isolated from gliadins and free glutenin polymers by asymmetrical flow field-flow fractionation in ethanol/water (50/50, v/v) [START_REF] Morel | Insight into gluten structure in a mild chaotropic solvent by asymmetrical flow fieldflow fractionation (AF4) and evidence of non-covalent assemblies between glutenin and ωgliadin[END_REF]. After HPSEC analysis in a SDS-buffer, which disrupted all the weak interactions, they were found to be composed of glutenin polymers and of fast eluting ω-gliadins (likely the two ones of higher molecular weight, identified by SDS-PAGE in Fig. 6).

The gluten protein assemblies, owing to their huge molecular size should be the first to phase separate as temperature is lowered, followed by free glutenin polymers and then free gliadins, in decreasing order of size and so ascending order of combinatorial entropy of mixing. Fig. 7ab presents the coexistence concentration curves for glutenin polymers (black), α/β and γgliadins (red) and ω-gliadin (green). Continuous lines are drawn as guide for the eyes to delineate the rich (full lines and symbols) and poor (dotted lines and empty symbols) branches arising from the liquid-liquid phase separation of the parent phase. For the glutenin polymers, a large concentration gap was observed between the rich and poor phases, even at a shallow quench depth (ΔT =1.5°C) (Fig. 7a), a feature theoretically consistent with the large size of the gluten protein assemblies and their low concentration in the Par extract [START_REF] Flory | Thermodynamics of heterogeneous polymers and their solutions[END_REF]. (Empty and filled symbols) poor and rich phases, dotted and continuous lines are guides for the eye to delineate the poor and rich branches, respectively; arrows (width in proportion with standard error) indicated proteins concentration in the parent phase. The concentration in the different proteins were calculated from the HPSEC profiles taking into account the following time limits: glutenin polymers (8.0-13.2 min), ω-gliadin (13.2-14.2min) and α/β and γ-gliadins (14.2-16.1).

Compared to the glutenin, the increases with ΔT of α/β and γ-gliadins concentration within the rich phases were first modest before intensifying beyond quench depths ΔT = 8-9°C (Fig. 4). Concomitantly, the gliadin poor branch showed an inflection point at ΔT=8-9°C (Tq=5°C) which may signal the onset of an additional demixing event. In support to this assumption, the demixing temperature of an almost pure sample of α/β and γ-gliadin, used in the same concentration range, was found to be about Tq=5.1±1.0 °C [START_REF] Morel | Insight into gluten structure in a mild chaotropic solvent by asymmetrical flow fieldflow fractionation (AF4) and evidence of non-covalent assemblies between glutenin and ωgliadin[END_REF]. For shallow quenches (0°C<ΔT<8°C ) the accumulation of α/β and γ-gliadins in the rich phase would t arise from their specific interaction with the glutenin polymers rather than from a decrease in solvent quality.

Poor branches of glutenin polymers and of ω-gliadin also showed an inflection point in between ΔT = 5-6°C (Fig. 7b), a feature that could signal the exhaustion of the gluten protein assemblies from the Par phase. Thereafter, for deeper quenches, free glutenin polymers and the two remaining ω-gliadins of smaller molecular weight identified above by SDS-PAGE, would begin to join the rich phase. Whether they will form new dense protein droplets or merge with the former ones is undecidable. However, multiphases coexistence is a hallmark of complex systems comprising several non-interacting components [START_REF] Feric | Coexisting liquid phases underlie nucleolar subcompartments[END_REF]). In the case of blends of multiple species (A and B) undergoing spinodal decomposition, as the quench depth increases, creation of inclusions of species B within an already formed droplet enriched in A, or alternatively successive trains of spinodal decomposition forming rich droplets of variable composition in A and B can be speculated [START_REF] Dignon | Biomolecular Phase Separation: From Molecular Driving Forces to Macroscopic Properties[END_REF]. At ΔT= 5°C, the concentration in glutenin polymers within the rich phase, stopped growing (vertical black branch) suggesting that the LLPS line could intersect the gelation line. Gelation of gluten protein samples, deriving from rich phases obtained after variable ΔT (4°C to 12°C) and suspended at 237 g/L in ethanol/water (50/50, v/v) was found to vary with the protein composition [START_REF] Banc | Phase separation dynamics of gluten protein mixtures[END_REF]. The samples comprising 120 and 145 g/L of glutenin polymers (ibid, Glu 57% and 66%) showed the typical mechanical spectra of gels, while the sample including only 99 g/L glutenin remained liquid (ibid, Glu ). In another study, the gluten proteins from the rich phase recovered at ΔT=10°C (red dot in Fig. 5ab) demonstrated spontaneous gelation at room temperature; in several hours to days depending on concentration [START_REF] Dahesh | Spontaneous gelation of wheat gluten proteins in a food grade solvent[END_REF]. The critical concentration for gelation, extrapolated in the limit of vanishing gelling time, was found to be 250±10 g/L (i.e., 125±5 g/L in glutenin, ibid). Based on these investigations, we propose that beyond a quench depth of 5°C and as the glutenin concentration within the rich phase attained 143 ± 6 g/L of glutenin, a sol-to-gel transition would occur. Whether gelation can precede, accompany or follow spinodal decomposition is well documented for synthetic polymers and for proteins in a binary solvent [START_REF] Bansil | Effects of gelation on spinodal decomposition kinetics in gelatin[END_REF][START_REF] Tanaka | Phase separation and gelation in gelatin gels[END_REF][START_REF] Kawanishi | Thermodynamic consideration of the sol-gel transition in polymer solutions[END_REF].

Thermodynamic insights of phase separation

For polymers carrying multiple sites for solvent interaction the balance between the enthalpy cost of mixing and the entropic gain is in favor of the former as the molecular size of the polymer grows and the temperature decreases. This property is at the basis of the fractionation of polydisperse polymers according to their size by liquid-liquid phase separation (LLPS). The lattice model for binary blends provides a simple formalism for fractionation of polymers by LLPS according to their molecular weight [START_REF] Flory | Thermodynamics of heterogeneous polymers and their solutions[END_REF]:

ln(𝜑 (𝑖) ′ 𝜑 (𝑖) ′′ ⁄ ) = -𝛼𝑥 (𝑖) (1)
where, prime and double prime refer to poor and rich phase, respectively and index (i) refers to individual molecular species, is the polymer volume fraction, ais the partitioning coefficient and x( i)

is the molar volume ratio of polymer (i) to the solvent:

𝑥 (𝑖) = 𝑀 𝑤(𝑖) ×𝑉 ̅ 𝑠𝑒𝑔 𝑀 𝑠𝑒𝑔 × 𝑉 ̅ 𝑠𝑜𝑙 (2)
In the frame of the lattice model, x( i) stands for the number of adjacent lattice sites a polymer molecule of size M w(i) occupies while a solvent molecule occupies one lattice site. Here M w(i) is the polymer (i) molecular weight, 𝑉 ̅ 𝑠𝑜𝑙 and 𝑉 ̅ 𝑠𝑒𝑔 are the molar volumes of the solvent and the polymer segment while M seg is the molecular weight of the polymer segment unit.

Using 𝜑 (𝑖) ′ = 𝐶 (𝑖) ′ ν ′ and 𝜑 (𝑖) ′′ = 𝐶 (𝑖) ′′ ν ′ ′ with C (i) the mass concentration of the molecular species (i) (in g/L) and ν' and ν'' the partial specific volume of glutenin polymers in the poor and rich phases, leads to :

ln ( 𝐶 (𝑖) ′ 𝐶 (𝑖) ′′ ) = -𝛼𝑥 (𝑖) -ln ( υ ′ υ ′′ ) (3)
Equation ( 1) applies for binary systems, composed of one solvent and one type of polydisperse polymer. For more complex systems involving several types of polymers or binary solvent, the Flory interaction parameter will not capture the multiple interactions at play. In our case, we found that the parent, poor and rich phases had the same ethanol content (checked by infrared spectroscopy, results not shown). Thus with respect to gluten proteins, the ethanol/water solvent would behave like a monomolecular solvent. Accordingly, its molar volume was calculated considering ideal mixing rules on a molar fraction basis. According to Eqs. 1 and 3 the partitioning of polymers between a rich phase and a poor phase follows an exponential law with their molecular weight; i.e. larger polymers are first to join the rich phase as temperature is lowered below the cloud temperature (T cloud ). Figure 8ab shows that Flory assumption well applies here since the logarithm of the phases concentration ratio linearly scales for glutenin polymers of molecular weights ranging from 80 kg/mol to 700 kg/mol. Beyond this size limit, elution of the glutenin polymers at the void volume of the HPSEC column, precluded accurate molecular weight determination.

As expected from lattice theory, the fitted slopes (α) increased with the quench depth (Fig. 8c) [START_REF] Flory | Thermodynamics of heterogeneous polymers and their solutions[END_REF].

According to equations 1 and 3, the origin of the linear fit should be zero, as long as the partial specific volume of the glutenin polymers is the same in the poor and rich phases (v '= v"). This was not the case in our experiments, especially when the quench depth increased. Fig. 8c shows that the ratio of the partial specific volumes of the glutenin polymers, as estimated from Eq. 3, increased up to 8 times when ΔT increased from 1.5 to 9 °C. Similar deviation from a zero intercept has already been reported for the fractionation of polydisperse polystyrene polymers in cyclohexane and was attributed to a coil to globule collapse of the polymer chains as they settled in the rich phases [START_REF] Breitenbach | Untersuchungen von Phasengleichgewichten an polymerlösungen mittels Säulenfraktionierung[END_REF]. Gluten protein assemblies detected in diluted suspensions of gluten proteins in ethanol/water (50/50, v/v) were identified as highly swollen with an apparent density of 5 to 7 mg/mL [START_REF] Dahesh | Polymeric assembly of gluten proteins in an aqueous ethanol solvent[END_REF][START_REF] Morel | Insight into gluten structure in a mild chaotropic solvent by asymmetrical flow fieldflow fractionation (AF4) and evidence of non-covalent assemblies between glutenin and ωgliadin[END_REF]. Their swelling factor is likely to depend on the osmotic pressure and transfer within the rich phase would lead to their shrinkage. (8a, 8b) Continuous lines are experimental data and dotted lines are the best fit with Eq.3. The typical HPSEC profile of the parent phase (upper insert) is given to help identify the considered respective ranges of α/β and γgliadins (light red area), -gliadin (light green area) and glutenin polymers (light blue area). The ln(C'/C'') parameter is the logarithm of the average protein concentration ratio obtained from 6 replicates. Only the glutenin polymers with Mw ranging from 98.7 to 700 kg/mol (2,519 < 𝑥 (𝑖) <17,760) were considered for data fitting according to Eq.3 (dotted lines). (8c) average partitioning coefficient (α, black symbols and polymer partial specific volume ratio (empty symbols) were calculated from Eq. 3, taking 𝑉 ̅ 𝑠𝑜𝑙 and 𝑉 ̅ 𝑠𝑒𝑔 equal to 27.6 and 79.2 mL/mol and M seg to 96.4 g/mol. Only regressions giving R 2 value > 0.99 are considered (ΔT=1.5 to 9 °C).

As shown in figure 8b, the behavior of ω-gliadins (Mw ~ 60 kg/mol) is very different from that of the other gliadins. Their partitioning toward the rich phase was more pronounced than it would have been for a glutenin polymer of similar molecular weight, since the experimental data laid below the linear fit. However, the feature was consistent with the privileged interaction of the fast eluting ω-gliadins with the glutenin polymer assemblies [START_REF] Morel | Insight into gluten structure in a mild chaotropic solvent by asymmetrical flow fieldflow fractionation (AF4) and evidence of non-covalent assemblies between glutenin and ωgliadin[END_REF]. For α/β and γ-gliadins (18 < Mw < 49.8 kg/mol) the partitioning did not follow the trend found for glutenin polymers. In other words, their tendency to join the rich phase appeared far lower than it would had been for a theoretical glutenin polymer of same mass (see the offset from the linear tendency, Fig. 8b).

To account for the multiple possible interactions between segments of polymer chains and solvent molecules through van der Waals forces, a single polymer-solvent interaction parameter (was introduced independently by Flory and Huggins [START_REF] Flory | Thermodynamics of high polymer solutions[END_REF][START_REF] Huggins | Solutions of long chain compounds[END_REF]. This phenomenological parameter, was calculated for the glutenin polymers, considering the fitted partitioning coefficients α (Eq. 3, Fig. 8c) and equation 4 where:

a 𝛼 = 𝜑 ′ (1 - 1 𝑥̅ 𝑛 ′ ) -𝜑 ′′ (1 - 1 𝑥̅ 𝑛 ′′ ) + 𝜒{(1 -𝜑 ′ ) 2 -(1 -𝜑 ′′ ) 2 } (4)
Where 𝑥̅ 𝑛 is the number average molecular size of the polymers within the range of interest in the poor (') and rich phase ('') and is the Flory-Huggins (FH) interaction parameter.

In order to estimate the protein volume fractions, the partial specific volume of the gluten proteins within the rich phases (ν") was taken as the inverse of the gluten protein density (1/ρ, 0.758 ml/g). To account for the shift in the partial specific volume ratios between rich and poor phases (Fig. 8c), the partial specific volume of the polymers in the poor phases (ν') were taken as 0.758 ml/g multiplied by the ratio ν'/ ν".

The numerical values of χ (half-filled red dots) are plotted against the inverse of the quenching temperature in Fig. 9. They exceeded 0.5 and tended to grow linearly when the quenching temperature decreased, at least until ΔT =7°C. This finding confirmed that glutenin polymers demixing was enthalpy driven, linked to a decreasing contribution of the mixing combinatorial entropy to the free energy of mixing as the solvent quality decreased.

The partitioning behavior of the α/β and γ-gliadins could also be analyzed on the basis of the Flory-Huggins lattice model, considering that the proteins eluted between 14.2 and 15.8 min (Fig. 4) behaved like a single polymer species, because of an overall low polydispersity index (PI=1.04). In that case, the Helmholtz free energy of mixing the pure protein in the solvent (Δf mix ), expressed per lattice site, is given by [START_REF] Flory | Thermodynamics of heterogeneous polymers and their solutions[END_REF]:

Δ𝑓 𝑚𝑖𝑥 (𝜑) = 𝑘𝑇[ 𝜑 𝑥 ln(𝜑) + (1 -𝜑) ln(1 -𝜑) + 𝜒𝜑(1 -𝜑)] (5)
with φ=φ' or φ'' and 𝑥 is the molar volume ratio of the polymer to the solvent (Eq. 2). Phase separation occurs when the free energy of mixing, instead of being a convex function of the polymer volume fraction, shows two local minima at φ = φ' and φ'', lying on a common straight line, whose slope is given by:

(Δ𝑓 𝑚𝑖𝑥 𝛿𝜑 (𝜑=𝜑 ′ ) ) ⁄ = (Δ𝑓 𝑚𝑖𝑥 𝛿𝜑 (𝜑=𝜑 ′′ ) ) ⁄ (6)
where the first derivative of the Helmholtz free energy of mixing is given by:

(𝛿Δ𝑓 𝑚𝑖𝑥 𝛿𝜑 ⁄ ) = 1 𝑥 ln(𝜑) + 1 𝑥 -ln(1 -𝜑) -1 + 𝜒(1 -2𝜑) with φ=φ' or φ'' (7)
Yielding:

𝜒 = 1 𝑥 ln( 𝜑 ′′ 𝜑 ′ )+ ln( 1-𝜑′ 1-𝜑′′ ) 2(𝜑 ′′ -𝜑 ′ ) (8)
The volume fractions for both phases were calculated from the known mass concentrations taking 0.758 mL/g as the gluten protein partial specific volume. The numerical values of χ as computed with Eq. 8 are given in Fig. 9 (full and half-filled black symbols). They displayed a monotonic evolution with 1/Tq roughly followed along a convex trend; first decreasing as the quench depth increased before increasing as ΔT increased further. In the limit of deep quenches (ΔT>9 °C), the χ values for the gliadin appeared in line with the trend observed for the glutenin polymers. The feature supports that interestingly the partitioning of the α/β and γ-gliadins and of glutenin polymers follows the same trend, albeit their different molecular weights.

From a theoretical standpoint, χ was originally envisioned as purely enthalpic in nature, with negative value indicating favorable interactions between polymer and solvent while when χ> 0.5, polymerpolymer interactions are favored at the expense of polymer-solvent mixing. In practice it is considered as a free-energy parameter including, both enthalpic ( H ) and entropic ( S ) contributions; χ S accounting for polymer conformational or solvent translational entropies. The enthalpic and entropic contributions of χ are estimated according to [START_REF] Schuld | Solvent quality as reflected in concentration-and temperaturedependent Flory-Huggins interaction parameters[END_REF]:

𝜒 = 𝜒 𝐻 + 𝜒 𝑠 with 𝜒 𝐻 = 𝑇 -1 (𝑑𝜒 𝑑𝑇 -1 ) ⁄ and 𝜒 𝑆 ≅ -Δ𝑆 (9)
Thus, the decrease of the χ values for the α/β and γ-gliadins, from T cloud up to ΔT = 7°C, derived from a negative χ H value counterbalanced by a positive χ s . In this range of quench depth, the partitioning of the α/β and γ-gliadins into the rich phase would be driven by their specific interactions with glutenin polymers (negative χ H ) in parallel to an increase of ordering (positive χ s ). Only at deeper quenches the partitioning of the α/β and γ-gliadins within the rich phases would become enthalpy driven (positive χ H ) as the solvent quality decreased.

The χ H andχ s contributions to the Flory-Huggins (FH) interaction parameter were determined by fitting the χ values obtained for the glutenin polymers and the ones obtained for the gliadin beyond ΔT > 7 °C (half-filled symbols, Fig. 9). according to:

𝜒 = 𝐴 𝑇 + 𝐵 ⁄ (10)
Best linear fit of the experimental data yielded χ H = 0.72 ±0.07 (at T cloud = 14°C) and χ S = -0.21±0.08.

Considering the negative value of χ S , it can be concluded that the solubility of gluten proteins in ethanol/water relies on a significant entropic gain, acting against phase separation. The feature is common in the case of co-solvency brought by binary solvents, composed of highly miscible solvent molecular species (like water and ethanol) able to interact with each other but also with the polymer chains. In this way, the addition of polymers increases the degree of disorder of the solvent compared with its pure state, thereby providing an appreciable entropic gain [START_REF] Zhang | Polymers with upper critical solution temperature behavior in alcohol/water solvent mixtures[END_REF].

The critical value of χ (χ c ) below which glutenin polymers and solvent will be miscible over the entire concentration range and the corresponding critical temperature (T C ) below which phase separation will occur, can be calculated according to [START_REF] Flory | Principles of polymer chemistry[END_REF]:

𝜒 𝑐 = 0.5 × (1 + 1 𝑥̅ 𝑛 0.5 ⁄ ) 2 (11)
and

𝑇 𝑐 = (𝜒 𝐻 ) (𝜒 𝑐 -𝜒 𝑠 ⁄ ) (in K) (12)
For the glutenin polymers, a theoretical critical temperature (T c , Eq. 11) of 13.2 °C (red star, Fig. 9) was calculated, closed to the cloud point temperature of 14°C, visually detected. The feature indicated that the concentration the glutenin polymers in the parent extract was closed to their critical concentration. This would also be the case for the α/β and γ-gliadins, since the calculated critical temperature (5°C, black star, Fig. 9) was also close to the 5.1°C obtained for a gluten protein sample including almost pure gliadins [START_REF] Morel | Insight into gluten structure in a mild chaotropic solvent by asymmetrical flow fieldflow fractionation (AF4) and evidence of non-covalent assemblies between glutenin and ωgliadin[END_REF]. The enthalpy of mixing gluten proteins with the ethanol/water solvent appeared size-independent as long as we admitted that the partial specific volume of the gluten protein assemblies could vary during LLPS. Thus, within the poor phase it would increase from 1.6 to 6.2 ml/g, as the protein concentration decreased from 41 to 26 g/L, while the gluten protein assemblies reaching the rich phase would undergo full deswelling.

The formalism followed to determine the range of χ values for the α/β and γ-gliadins undergoing LLPS at different quench depths, cannot be applied for the ω-gliadins, since HPSEC did not allow to discriminate them according to their specific involvement in the gluten protein assemblies or not.

Fig. 9: Flory-Higgins interaction parameters for gliadins and glutenin polymers

Conclusion

Gluten, an industrial commodity derived from wheat flour, is mainly composed of gliadin and glutenin polymers that were found in the present study to be partly soluble (by 45-50%) in an ethanol/water solvent (50/50, v/v). Even if gliadin was the main class of wheat storage proteins solubilized, the extract also included 30% of glutenin polymers. Their molecular weight distribution range did not differ from the ones of the SDS-soluble glutenin polymers from wheat flours. The gluten/solvent ratio (10-20 g/100 mL) and the ethanol content (48 to 66 %, v/v) of the solvent did not significantly impacted the extraction yield and composition of the extract. Due to the low solid/solvent ratio chosen (10 g/100mL) the obtained soluble extract was not very concentrated (~45 g/L), a situation very favorable for polymer fractionation by liquid-liquid phase separation (LLPS) induced by a temperature quench. We showed that dense phases substantially enriched in glutenin polymers could be recovered by applying shallow temperature quenches, a few degrees below the cloud temperature of the extract (T cloud = 14°C). Because of the very unbalanced volumes of the coexisting phases, fractions of high protein content (up to 384 g/L) and of tunable glutenin over gliadin ratio (from 0.5 to 2.5) could be recovered in a few hours. The fractionation method paves the way for a rigorous assessment of the impact of glutenin polymers content on the rheological properties of gluten proteins. Thermodynamic insights on gluten protein LLPS were examined at the level of each gluten protein classes. We found that the partitioning of the glutenin polymers between the rich and poor phases varied exponentially with their molecular weight, consistently with the lattice theory model introduced by Flory (Flory, 1994). Contrarily to what was observed for the α/β and γ-gliadins undergoing LLPS, the partial specific volume of the glutenin polymers varied whether they settled in the poor or in the rich phase and the difference widened when the quench depth increased. The change in the partial specific volume of the glutenin polymers during LLPS would be related to their involvement in highly swollen protein assemblies. It could even play a role in relation to the rising idea that de-solvation entropy drives liquid-liquid phase separation of proteins [START_REF] Harmon | Differential solvation of intrinsically disordered linkers drives the formation of spatially organized droplets in ternary systems of linear multivalent proteins[END_REF][START_REF] Park | Dehydration entropy drives liquid-liquid phase separation by molecular crowding[END_REF].

Thanks to the analysis of the temperature dependency of the enthalpic and entropic contributions of the Flory-Huggins interaction parameters, we were able to predict the critical temperatures of the different gluten protein classes. We showed that interactions between α/β and γ-gliadins and glutenin polymers caused the accumulation of the former in the rich phase, while the quench temperature was above the critical temperature of gliadins. At lower temperatures, for deeper quench depths (ΔT>9°C), we demonstrated that the α/β-and γ-gliadins remaining in the parent phase reached a metastable state and in turn phase separated, as individual sepcies. Thus within the rich phase droplets already formed at shallow quench depths, further decomposition might occur leading to the internal clustering of the gliadins. The glutenin concentration in the protein rich phase was found to plateau at about 143 g/L for quench temperature below 9°C. In our opinion, this would coincide to the critical glutenin polymer concentration leading to gelation, irrespective of the glutenin/gliadin ratio. Whether the transition from critical to viscoelastic gel occurred in the course of the present LLPS experiments, performed at a relatively slow quenching rate would deserve further exploration. Overall, fine-tuning of gluten protein LLPS in aqueous ethanol solutions might be a way to prepare micro-sized droplets with heterogeneous internal structure. This could open the gate to a better understanding of formation of cereal storage protein bodies, which may show (maize, sorghum) or not (wheat) a spatial gradient in the storage proteins distribution from the core to the periphery or even inner droplets of specific composition [START_REF] Lending | Structure of maize protein bodies and immunocytochemical localization of zeins[END_REF]. The interest in LLPS of storage proteins is a worthwhile cause, inasmuch that internal structure of protein bodies impact cereal protein digestibility [START_REF] Duodu | Effect of grain structure and cooking on sorghum and maize in vitro protein digestibility[END_REF][START_REF] Oria | A highly digestible sorghum mutant cultivar exhibits a unique folded structure of endosperm protein bodies[END_REF].

  Fig. 1ab: Sketche of the extraction procedure (a) and protein composition of the gluten powder, parent and gel phases (b), physical aspect of the Gel (c).

  (a) Gluten powder extracted with 50/50 ethanol/water (v/v) at 10 g/100 mL (cf.2.2). (b) Protein composition of gluten powder and of the Parent and Gel fractions as determined by HPSEC allowing to quantify the content in SDS-soluble (glu) glutenin, gliadin (gli), albumin/globulin (al/glo) and SDS-insoluble glutenin (glu *) (cf. 2.3). (c) Pictures of the Gel, deformed with a spatula (top) and relaxed after removing the spatula (bottom).

  Fig. 3abc: Protein concentration (a), percentage of the Parent volume (b) and percentage of the Parent protein (c), recovered in the rich and poor phases following a quench of depth ΔT.

Fig. 4 .

 4 Fig. 4. Protein size distribution profiles of Parent and demixed phases obtained at different quench depths.Quench depth (ΔT=T cloud -Tq) is indicated in each graph. Averaged HPSEC-UV trace from 6 replicates was converted into protein concentration and then into protein mass taking the volume of the phases into consideration. Results are given for a total of 4.57 g of proteins in 100 mL of parent phase. (Yellow) Parent; (Orange) rich phase; (blue) poor phase.

  elution time (min) elution time (min) elution time (min) elution time (min) elution time (min) elution time (min) elution time (min) elution time (min) elution time (min) elution time (min) elution time (min)

Fig. 6 .

 6 Fig. 6. SDS-PAGE analysis of gluten, parent and demixed phases.Electrophoresis in reducing conditions, 30 µg of protein per well. (Glu) gluten powder, (Par) parent phase, (P) poor and (R) rich phases. The quenching depth (ΔT) is indicated on the top of the picture of the gel. (HMW-GS, LMW-GS) high-and low-molecular-weight glutenin subunits, respectively. On the right of the picture the densitometry profile of the gluten sample is given as a typical example.

a)

  results from HPSEC analysis (Fig.4 and 5a). LMW-GS = [glutenin polymers a -HMW-GS*].

Fig. 7ab .

 7ab Fig. 7ab. Coexistence curves for glutenin polymers (black, diamond), α/β and γ-gliadins (red, square) and ω-gliadin (green, circle).

Fig. 8abc :

 8abc Fig. 8abc: Partitioning of proteins in rich and poor phases according to quenching temperatures (a,b) and parameters fitted according to Eq. 3 for glutenin polymers (c).

Table 1 :

 1 Gluten extraction with ethanol/water (50/50, v/v) at variable solid/solvent ratios.

	b c in mass percent of the total sample protein content. (glu) stands for glutenin polymers, (a) Gram of gluten powder (84.5% protein, wb) per 100 mL of added solvent (50/50, ethanol/water, v/v). concentration in g/L.
	(gli) for gliadins.

Table 2 :

 2 Proportions in gluten protein types estimated from SDS-PAGE densitometry (Fig.6) as compared to HPSEC analysis.

	Samples	HMW-GS	*	ω-gli*	LMW-GS γ-gli*	LMW-GS α/β-gli*	al/glo*	Glutenin polymers a	HMW-GS*/LMW-GS ratio
	Gluten	13.67		8.20	40.07	29.94	8.12	42.50	0.47
	Par	9.55		9.69	46.63	26.36	7.78	31.25	0.44
	p-ΔT 12°C	0.17		3.22	60.78	25.88	9.95	3.37	0.05
	p-ΔT 8°C	3.26		4.86	58.27	25.86	7.75	12.86	0.34
	r-ΔT 12°C	14.15		10.22	41.16	27.08	7.39	42.21	0.50
	r-ΔT 8°C	17.80		11.53	36.30	26.60	7.78	52.60	0.51
	r-ΔT 5°C	20.33		10.36	33.79	27.14	8.37	62.87	0.48
	(*) Relative amounts deduced from gel densitometry. (

Acknowledgments

Joelle Bonicel is warmly thanks for her technical assistance in the achievement of HPSEC analysis. We acknowledge the Doctoral School GAIA and the Labex Numev (ANR-10-LABX-20) for the funding of the PhD Grant of J. Pincemaille.


Trinh Van Chien

Hien Quoc Ngo

Symeon Chatzinotas

Marco Di Renzo

Björn Ottersten

RIS and Cell-Free Massive MIMO: A Marriage For Harsh Propagation Environments

This paper considers Cell-Free Massive Multiple Input Multiple Output (MIMO) systems with the assistance of an RIS for enhancing the system performance. Distributed maximum-ratio combining (MRC) is considered at the access points (APs). We introduce an aggregated channel estimation method that provides sufficient information for data processing. The considered system is studied by using asymptotic analysis which lets the number of APs and/or the number of RIS elements grow large. A lower bound for the channel capacity is obtained for a finite number of APs and engineered scattering elements of the RIS, and closed-form expression for the uplink ergodic net throughput is formulated. In addition, a simple scheme for controlling the configuration of the RIS scattering elements is proposed. Numerical results verify the effectiveness of the proposed system design and the benefits of using RISs in Cell-Free Massive MIMO systems are quantified.

I. INTRODUCTION

Cell-Free Massive Multiple Input Multiple Output (MIMO) has recently been introduced to reduce the intercell interference of colocated Massive MIMO architectures. This is a network deployment where a large number of access points (APs) are located in a given coverage area to serve a small number of users [START_REF] Ngo | Cell-free massive MIMO versus small cells[END_REF]. All the APs collaborate with each other via a backhaul network and serve all the users in the absence of cell boundaries. The system performance is enhanced in Cell-Free Massive MIMO systems because they inherit the benefits of the distributed MIMO and network MIMO architectures, but the users are also close to the APs. When each AP is equipped with a single antenna, maximumratio combining (MRC) results in a good net throughput for every user, while ensuring a low computational complexity and offering a distributed implementation that is convenient for scalability purposes. However, this network deployment cannot guarantee a good service under harsh propagation environments.

RIS is an emerging technology that is capable of shaping the radio waves at the electromagnetic level without applying digital signal processing methods and requiring power amplifiers [START_REF] Le | Robust probabilisticconstrained optimization for IRS-aided MISO communication systems[END_REF]. Each element of the RIS scatters (e.g., reflects) the incident signal without using radio frequency chains and This work of T. V. Chien, S. Chatzinotas, and B. Ottersten was supported by RISOTTI -Reconfigurable Intelligent Surfaces for Smart Cities under project FNR/C20/IS/14773976/RISOTTI. The work of H. Q. Ngo was supported by the UK Research and Innovation Future Leaders Fellowships under Grant MR/S017666/1. The work of M. Di Renzo was supported in part by the European Commission through the H2020 ARIADNE project under grant agreement number 871464 and through the H2020 RISE-6G project under grant agreement number 101017011. The long version of this paper was submitted to the IEEE Trans. Wireless Comm. [START_REF] Chien | Reconfigurable intelligent surface-assisted Cell-Free Massive MIMO systems over spatially-correlated channels[END_REF]. power amplification. Integrating an RIS into wireless networks introduces digitally controllable links that scale up with the number of engineered scattering elements of the RIS, whose estimation is, however, challenged by the lack of digital signal processing units at the RIS. For simplicity, the main attention has so far been concentrated on designing the phase shifts with perfect channel state information (CSI) [START_REF] Zhao | Intelligent reflecting surface enhanced wireless network: Two-timescale beamforming optimization[END_REF], [START_REF] Perović | Achievable rate optimization for MIMO systems with reconfigurable intelligent surfaces[END_REF] and the references therein. As far as the integration of Cell-Free Massive MIMO and RIS is concerned, recent works have formulated and solved optimization problems with different communication objectives under the assumption of perfect (and instantaneous) CSI [START_REF] Zhou | Achievable rate optimization for aerial intelligent reflecting surface-aided Cell-Free Massive MIMO system[END_REF], [START_REF] Zhang | Capacity improvement in wideband reconfigurable intelligent surface-aided cell-free network[END_REF]. Recent results have shown that designs for the phase shifts of the RIS elements based on statistical CSI may be of practical interest and provide good performance [START_REF] Van Chien | Outage probability analysis of IRSassisted systems under spatially correlated channels[END_REF]. In the depicted context, no prior work has analyzed an RIS-assisted Cell-Free Massive MIMO system in the presence of spatially-correlated channels.

In this work, we consider an RIS-assisted Cell-Free Massive MIMO under spatially correlated channels. We exploit a channel estimation scheme that estimates the aggregated channels including both the direct and indirect links. We analytically show that, even by using a low complexity MRC technique, the non-coherent interference, small-scale fading effects, and additive noise are averaged out when the number of APs and RIS elements increases. The received signal includes, hence, only the desired signal and the coherent interference. We derive a closed-form expression of the net throughput for the uplink data transmission. The impact of the array gain, coherent joint transmission, channel estimation errors, pilot contamination, spatial correlation, and phase shifts of the RIS, which determine the system performance, are explicitly observable in the obtained analytical expressions. With the aid of numerical simulations, we verify the effectiveness of the proposed channel estimation scheme and the accuracy of the closed-form expression of the net throughput. The obtained numerical results show that the use of RISs enhance the net throughput per user significantly, especially when the direct links are blocked with high probability.

Notation: Upper and lower bold letters denote matrices and vectors. The identity matrix of size N × N is denoted by I N . (•) * , (•) T , and (•) H are the complex conjugate, transpose, and Hermitian transpose. E{•} and Var{•} denote the expectation and variance of a random variable. The circularly symmetric Gaussian distribution is denoted by CN (•, •) and diag(x) is the diagonal matrix whose main diagonal is given by x. tr(•) is the trace operator. The Euclidean norm of vector x is x , and X is the spectral norm of matrix X. Finally, mod(•, •) is the modulus operation and • denotes the truncated argument.

II. SYSTEM MODEL, CHANNEL ESTIMATION, AND RIS

PHASE-SHIFT CONTROL We consider an RIS-assisted Cell-Free Massive MIMO system, where M APs connected to a central processing unit (CPU) serve K users on the same time and frequency resource. All APs and users are equipped with a single antenna and they are randomly located in the coverage area. The communication is assisted by an RIS that comprises N engineered scattering elements that can modify the phases of the incident signals. The matrix of phase shifts of the RIS is denoted by Φ Φ Φ = diag [e jθ1 , . . . , e jθ N ] T , where θ n ∈ [-π, π] is the phase shift applied by the n-th element of the RIS.

A. Channel Model

We assume a quasi-static block fading model where each coherence interval comprises τ c symbols. The APs have knowledge of only the channel statistics instead of the instantaneous channel realizations. Also, τ p symbols (τ p < τ c ) in each coherence interval are dedicated to the channel estimation and the remaining (τ c -τ p ) symbols are the data transmission.

The following notation is used: g mk is the channel between the user k and the AP m, which is the direct link [START_REF] Wu | Intelligent reflecting surface enhanced wireless network via joint active and passive beamforming[END_REF]; h m ∈ C N is the channel between the AP m and the RIS; and z k ∈ C N is the channel between the RIS and the user k. In this paper, we consider a realistic channel model by taking into account the spatial correlation among the scattering elements of the RIS, which is due to their sub-wavelength size, sub-wavelength inter-distance, and geometric layout. In an isotropic propagation environment, in particular, g mk , h m , and z k can be modeled as follows

g mk ∼ CN (0, β mk ), h m ∼ CN (0, R m ), z k ∼ CN (0, R k ), (1)
where β mk is the large-scale fading coefficient; R m ∈ C N ×N and R k ∈ C N ×N are the covariance matrices. The covariance matrices in [START_REF] Chien | Reconfigurable intelligent surface-assisted Cell-Free Massive MIMO systems over spatially-correlated channels[END_REF] correspond to a general model, which can be further particularized for application to typical RIS designs and propagation environments. A correlation model that is applicable to isotropic scattering with uniformly distributed multipath components in the half-space in front of the RIS was recently reported in [START_REF] Björnson | Rayleigh fading modeling and channel hardening for reconfigurable intelligent surfaces[END_REF], whose covariance matrices are

R m = α m d H d V R and R k = αk d H d V R, (2)
where α m , αmk ∈ C are the large-scale channel coefficients. The matrices in (2) assume that the size of each RIS element is d H × d V , with d H being the horizontal width and d V being the vertical height of each RIS element. In particular, the (l, t)-th element of the spatial correlation matrix

R ∈ C N ×N in (2) is [R] lt = sinc(2 u l -u t /λ)
, where λ is the wavelength and sinc(x) = sin(πx)/(πx) is the sinc function. The vector u x , x ∈ {l, t} is given by u

x = [0, mod (x -1, N H )d H , (x -1)/N H d V ] T ,
where N H and N V denote the total number of RIS elements in each row and column, respectively.

B. Uplink Pilot Training Phase

The channels are independently estimated from the τ p pilot sequences transmitted by the K users. All the users share the same τ p pilot sequences. In particular, φ φ φ k ∈ C τp with φ φ φ k 2 = 1 is defined as the pilot sequence allocated to the user k. We denote by P k the set of indices of the users (including the user k) that share the same pilot sequence as the user k. The pilot sequences are assumed to be mutually orthogonal such that the pilot reuse pattern is

φ φ φ H k φ φ φ k = 1, k ∈ P k . Otherwise, φ φ φ H k φ φ φ k = 0.
During the pilot training phase, all the K users transmit the pilot sequences to the M APs simultaneously. In particular, the user k transmits the pilot sequence √ τ p φ φ φ k . The received training signal at the AP m can be written as

y pm = K k=1 √ pτ p g mk φ φ φ k + K k=1 √ pτ p h H m Φ Φ Φz k φ φ φ k + w pm , (3)
where p is the normalized signal-to-noise ratio (SNR) of each pilot symbol, and w pm ∈ C τp is the additive noise at the AP m, which is distributed as w pm ∼ CN (0, I τp ). In order for the AP m to estimate the desired channels from the user k, the received training signal in (3) is projected on φ φ φ H k as

y pmk = φ φ φ H k y pm = √ pτ p g mk + h H m Φ Φ Φz k + k ∈P k \{k} √ pτ p g mk + h H m Φ Φ Φz k + w pmk , (4)
where w pmk = φ φ φ H k w pm ∼ CN (0, 1). We emphasize that the co-existence of the direct and indirect channels due to the presence of the RIS results in a complicated channel estimation process. In particular, the cascaded channel in (4) results in a nontrivial procedure to apply the minimum mean-square error (MMSE) estimation method, as reported in previous works, for processing the projected signals [START_REF] Ngo | Cell-free massive MIMO versus small cells[END_REF], [START_REF] Van Chien | Joint power allocation and load balancing optimization for energy-efficient cell-free Massive MIMO networks[END_REF]. Based on the specific signal structure in (4), we denote the channel between the AP m and the user k through the RIS as

u mk = g mk + h H m Φ Φ Φz k , (5)
which is referred to as the aggregated channel that comprises the direct and indirect link between the user k and the AP m. By capitalizing on the definition of the aggregated channel in ( 5), the required channels can be estimated in an effective manner even in the presence of the RIS. In particular, the aggregated channel in ( 5) is given by the product of weighted complex Gaussian and spatially correlated random variables, as given in [START_REF] Chien | Reconfigurable intelligent surface-assisted Cell-Free Massive MIMO systems over spatially-correlated channels[END_REF]. Conditioned on the phase shifts, we employ the linear MMSE method for estimating u mk at the AP. Despite the complex structure of the RIS-assisted channels, Lemma 1 provides analytical expressions of the estimated channels.

Lemma 1. By assuming that the AP m employs the linear MMSE estimation method based on the observation in (4), the estimate of the aggregate channel u mk is formulated as ûmk = E{y * pmk u mk }y pmk /E{|y pmk | 2 } = c mk y pmk , (6) where c mk = E{y * pmk u mk }/E{|y pmk | 2 } has the following closed-form expression

c mk = √ pτ p β mk + tr Φ Φ Φ H R m Φ Φ Φ R k pτ p k ∈P k β mk + tr Φ Φ Φ H R m Φ Φ Φ R k + 1 . (7)
The estimated channel in (6) has zero mean and variance γ mk equal to

γ mk = E{|û mk | 2 } = √ pτ p β mk + tr Φ Φ Φ H R m Φ Φ Φ R k c mk .
(8) Also, the channel estimation error e mk = u mk -ûmk and the channel estimate ûmk are uncorrelated. The channel estimation error has zero mean and variance equal to

E |e mk | 2 = β mk + tr Φ Φ Φ H R m Φ Φ Φ R k -γ mk . (9)
Proof. It is similar to the proof in [START_REF] Kay | Fundamentals of Statistical Signal Processing: Estimation Theory[END_REF], and is obtained by applying similar analytical steps to the received signal in [START_REF] Zhao | Intelligent reflecting surface enhanced wireless network: Two-timescale beamforming optimization[END_REF] and by taking into account the structure of the RIS-assisted channel and the spatial correlation matrices in [START_REF] Chien | Reconfigurable intelligent surface-assisted Cell-Free Massive MIMO systems over spatially-correlated channels[END_REF].

Lemma 1 shows that, by assuming Φ Φ Φ fixed, the aggregated channel in (5) can be estimated without increasing the pilot training overhead, as compared to a conventional Cell-Free Massive MIMO system. The obtained channel estimate in (6) unveils the relation ûmk = c mk c mk ûmk if the user k uses the same pilot sequence as the user k. Because of pilot contamination, it may be difficult to distinguish the signals of these two users. In the following, the analytical expression of the channel estimates in Lemma 1 are employed for signal detection in the uplink data transmission. They are used also to optimize the phase shifts of the RIS in order to minimize the channel estimation error and to evaluate the corresponding ergodic net throughput.

C. RIS Phase-Shift Control

Channel estimation is a critical aspect in Cell-Free Massive MIMO. As discussed in previous text, in many scenarios, non-orthogonal pilots have to be used. This causes pilot contamination, which may reduce the system performance significantly. In this section, we design an RIS-assisted phase shift control scheme that is aimed to improve the quality of channel estimation. To this end, we introduce the normalized mean square error (NMSE) of the channel estimate of the user k at the AP m as follows

NMSE mk = E{|e mk | 2 }/E{|u mk | 2 } = 1 - pτ p β mk + tr(Φ Φ Φ H R m Φ Φ Φ R k ) pτ p k ∈P k β mk + tr(Φ Φ Φ H R m Φ Φ Φ R k ) + 1 . ( 10
)
where the last equality is obtained from [START_REF] Wu | Intelligent reflecting surface enhanced wireless network via joint active and passive beamforming[END_REF]. We optimize the phase shift matrix Φ Φ Φ of the RIS so as to minimize the total NMSE obtained from all the users and all the APs as follows

minimize {θn} M m=1 K k=1 NMSE mk subject to -π ≤ θ n ≤ π, ∀n. (11)
The optimal phase shifts solution to problem [START_REF] Van Chien | Joint power allocation and load balancing optimization for energy-efficient cell-free Massive MIMO networks[END_REF] is obtained by exploiting the statistical CSI that include the large-scale fading coefficients and the covariance matrices. Problem ( 11) is a fractional program, whose globally-optimal solution is not simple to be obtained for an RIS with a large number of independently tunable elements. Nonetheless, in the special network setup where the direct links from the APs to the users are weak enough to be negligible with respect to the RISassisted links, the optimal solution to problem [START_REF] Van Chien | Joint power allocation and load balancing optimization for energy-efficient cell-free Massive MIMO networks[END_REF] is available in a closed-form expression as summarized in Corollary 1.

Corollary 1. If the direct links are weak enough to be negligible and the RIS-assisted channels are spatially correlated as formulated in (2), the optimal maximizer of the optimization problem in (11) is θ 1 = . . . = θ N , i.e., the equal phase shift design is optimal.

Proof. The proof follows by analyzing the objective function of problem [START_REF] Van Chien | Joint power allocation and load balancing optimization for energy-efficient cell-free Massive MIMO networks[END_REF] with respect to the phase-shift elements. The detailed proof is available in the journal version [START_REF] Chien | Reconfigurable intelligent surface-assisted Cell-Free Massive MIMO systems over spatially-correlated channels[END_REF].

Corollary 1 provides a simple but effective option to design the phase shifts of the RIS while ensuring the optimal estimation of the aggregated channels according to the sum-NMSE minimization criterion, provided that the direct link are completely blocked and the spatial correlation model in (2) holds true. Therefore, an efficient channel estimation protocol can be designed even in the presence of an RIS with a large number of engineered scattering elements. The numerical results in Section IV show that the phase shift design in Corollary 1 offers good gains in terms of net throughput even if the direct links are not negligible.

III. UPLINK DATA TRANSMISSION AND PERFORMANCE ANALYSIS WITH MR COMBINING In this section, we introduce a procedure to detect the uplink transmitted signals and derive an asymptotic closedform expression of the ergodic net throughput.

A. Uplink Data Transmission Phase

In the uplink, all the K users transmit their data to the M APs simultaneously. Specifically, the user k transmits a modulated symbol s k with E{|s k | 2 } = 1. This symbol is weighted by a power control factor √ η k , 0 ≤ η k ≤ 1. Then, the received baseband signal, y um ∈ C, at the AP m is

y m = √ ρ K k=1 √ η k u mk s k + w m , ( 12
)
where ρ is the normalized uplink SNR of each data symbol and w m is the normalized additive noise with w m ∼ CN (0, 1). For data detection, the MRC method is used at the CPU, i.e., ûmk , ∀m, k, in ( 6) is employed to detect the data transmitted by the user k. In mathematical terms, the corresponding decision statistic is

r k = √ ρ M m=1 K k =1 √ η k û * mk u mk s k + M m=1 û * mk w m . ( 13
)
Based on the observation r k , the uplink ergodic net throughput of the user k is analyzed in the next subsection.

B. Asymptotic Analysis

Since the number of APs, M , and the number of tunable elements of the RIS, N , can be large, we analyze the performance of two case studies: (i) N is fixed and M is large; and (ii) both N and M are large. The asymptotic analysis is conditioned upon a given setup of the CSI. To this end, the uplink weighted signal in ( 13) is split into three terms based on the pilot reuse set P k , as follows

r k = √ ρ k ∈P k M m=1 √ η k û * mk u mk s k T k1 + √ ρ k / ∈P k M m=1 √ η k û * mk u mk s k T k2 + M m=1 û * mk w m T k3 , (14)
where T k1 accounts for the signals received from all the users in P k , and T k2 accounts for the mutual interference from the users that are assigned orthogonal pilot sequences. The impact of the additive noise obtained after applying MR combining is given by T k3 . From ( 4)-( 6), we obtain the following identity

M m=1 √ η k û * mk u mk = k ∈P k \{k } M m=1 √ η k pτ p c mk u mk u * mk + M m=1 √ η k pτ p c mk |u mk | 2 + M m=1 √ η k c mk u mk w * pmk , (15)
1) Case I: N is fixed and M is large, i.e., M → ∞. In this case, we divide both sides of (15) by M and exploits Tchebyshev's theorem [START_REF] Cramér | Random variables and probability distributions[END_REF] 1 to obtain

1 M M m=1 √ η k û * mk u mk P ----→ M →∞ 1 M M m=1 √ η k pτ p c mk β mk + tr Φ Φ Φ H R m Φ Φ Φ R k , ( 16
)
where P -→ denotes the convergence in probability. 2 Note that the second and third terms in (15) converge to zero. By 1 Let X 1 , . . . , Xn be independent random variables such that E{X i } = xi and Var{X i } ≤ c < ∞. Then, Tchebyshev's theorem states

1 n n n =1 X n P ----→ n→∞ 1 n n xn .
2 A sequence {Xn} of random variables converges in probability to the random variable X if, for all > 0, it holds that limn→∞ Pr(|Xn -X| > ) = 0, where Pr(•) denotes the probability of an event.

inserting ( 16) into the decision variable in ( 14), we obtain the following deterministic value

1 M r k P ----→ M →∞ 1 M k ∈P k M m=1 √ η k pτ p ρ u c mk β mk +tr Φ Φ Φ H R m Φ Φ Φ R k s k , (17)
because T k2 /M → 0 and T k3 /M → 0 as M → ∞.

The result in (17) unveils that, for a fixed N , the channels become asymptotically orthogonal. In particular, the small scale fading, the non-coherent interference, and the additive noise vanish. The only residual impairment is the pilot contamination caused by the users that employ the same pilot sequence. Due to pilot contamination, the system performance cannot be improved by adding more APs if MRC is used. The contributions of both the direct and RIS-assisted indirect channels appear explicitly in (17) through β mk and

tr(Φ Φ Φ H R m Φ Φ Φ R k ), respectively .
2) Case II: Both N and M are large, i.e., N → ∞ and M → ∞. We first need some assumptions on the covariance matrices R m and R k , as summarized as follows.

Assumption 1. For m = 1, . . . , M and k = 1, . . . , K, the covariance matrices R m and R k are assumed to fulfill the following properties

lim sup N R m 2 < ∞, lim inf N 1 N tr(R m ) > 0, (18)
lim sup N R k 2 < ∞, lim inf N 1 N tr( R k ) > 0. (19)
The assumptions in (18) and (19) imply that the largest singular value and the sum of the eigenvalues (counted with their mutiplicity) of the N × N covariance matrices that characterize the spatial correlation among the channels of the RIS elements are finite and positive. Dividing both sides of (15) by M N and applying Tchebyshev's theorem, we obtain

1 M N M m=1 √ η k û * mk u mk P ----→ M →∞ N →∞ 1 M N M m=1 √ η k pτ p ρ u c mk tr Φ Φ Φ H R m Φ Φ Φ R k . (20) We observe that Φ Φ Φ H R m Φ Φ Φ R k is similar to R 1/2 k Φ Φ ΦR m Φ Φ Φ H R 1/2
k , which is a positive semi-definite matrix. 3 Because similar matrices have the same eigenvalues, it follows that tr Φ Φ Φ H R m Φ Φ Φ R k > 0. Based on Assumption 1, we obtain the following inequalities

1 N tr Φ Φ Φ H R m Φ Φ Φ R k (a) ≤ 1 N Φ Φ Φ 2 tr R m Φ Φ Φ R k (b) = 1 N tr Φ Φ Φ R k R m (c) ≤ 1 N R k 2 tr(R m ), (21)
where (a) is obtained by an inequality on the trace of the product of matrices; (b) follows because Φ Φ Φ 2 = 1; and (c) is obtained from Assumption 1. Based on Assumption 1, in addition, the last inequality in ( 14) is bounded by a positive constant. From (20) and ( 21), therefore, the decision variable in (14) can be formulated as

1 M N r k P ----→ M →∞ N →∞ 1 M N k ∈P k M m=1 √ η k pτ p ρ u c mk tr Φ Φ Φ H R m Φ Φ Φ R k s k . ( 22
)
The expression obtained in (22) reveals that, as M, N → ∞, the post-processed signal at the CPU consists of the desired signal of the intended user k and the interference from the other users in P k . Compared with (17), we observe that ( 22) is independent of the direct links and depends only on the RIS-assisted indirect links. This highlights the potentially promising contribution of an RIS, in the limiting regime M, N → ∞, for enhancing the system performance.

C. Uplink Ergodic Net Throughput Analysis with a Finite Number of APs and Phase Shifts

We now focus our attention on the practical setup in which M and N are both finite. By utilizing the user-and-then forget channel capacity bounding method [START_REF] Chien | Reconfigurable intelligent surface-assisted Cell-Free Massive MIMO systems over spatially-correlated channels[END_REF], the uplink ergodic net throughput of the user k can be computed in a closed-form expression for (23) as given in Theorem 1.

Theorem 1. If the CPU utilizes the MRC method, a lower bound closed-form expression for the uplink net throughput of the user k is given as follows

R k = Bν (1 -τ p /τ c ) log 2 (1 + SINR k ) , [Mbps], ( 23
)
where B is the system bandwidth measured in MHz and 0 ≤ ν ≤ 1 is the portion of each coherence interval that is dedicated to the uplink data transmission. The effective uplink signal-to-noise-plus-interference ratio (SINR) is

SINR k = ρη k M m=1 γ mk 2 (MI k + NO k ), (24)
where MI k is the mutual interference and the noise denoted by NO k are, respectively, given by

MI k = ρ u K k =1 M m=1 η k γ mk δ mk + pτ p ρ u k ∈P k M m=1 η k c 2 mk × tr(Θ Θ Θ 2 mk ) + pτ p ρ u K k =1 k ∈P k M m=1 M m =1 η k c mk c m k × tr(Θ Θ Θ mk Θ Θ Θ m k ) + pτ p ρ u k ∈P k \{k} η k M m=1 c mk δ mk 2 , (25)
NO k = M m=1 γ mk , (26)
with δ mk = β mk + tr Φ Φ Φ H R m Φ Φ Φ R k , c mk given in [START_REF] Zhang | Capacity improvement in wideband reconfigurable intelligent surface-aided cell-free network[END_REF], and γ mk given in [START_REF] Van Chien | Outage probability analysis of IRSassisted systems under spatially correlated channels[END_REF].

Proof. The main idea of proof is to average out the randomness by using the use-and-then-forget capacity bounding technique and fundamental properties of Massive MIMO. The detailed proof is available in the journal version [START_REF] Chien | Reconfigurable intelligent surface-assisted Cell-Free Massive MIMO systems over spatially-correlated channels[END_REF].

By direct inspection of the SINR in (24), the numerator increases with the square of the sum of the variances of the channel estimates, γ mk , ∀m, thanks to the joint coherent transmission. On the other hand, the first term in the denominator represents the power of the interference. Due to the limited and finite number of orthogonal pilot sequences being used, it represents the impact of pilot contamination. The last term is the additive noise. The SINR in (24) is a multivariate function of the matrix of phase shifts of the RIS and of the channel statistics, i.e., the channel covariance matrices. Compared with conventional Cell-Free Massive MIMO systems, the strength of the desired signal increases thanks to the assistance of an RIS. However, the coherent and non-coherent interference become more severe as well, due to the need of estimating both the direct and indirect links in the presence of an RIS.

IV. NUMERICAL RESULTS We consider a geographic area of size 1 km 2 that is wrapped around at the edges. The locations of 100 APs and 10 users are given in terms of (x, y) coordinates. To simulate a harsh communication environment, the APs are uniformly distributed in the sub-region x, y ∈ [-0.75, -0.5] km, while the users are uniformly distributed in the sub-region x, y ∈ [0.375, 0.75] km. An RIS with N = 900 is located at the origin, i.e., (x, y) = (0, 0). Each coherence interval comprises τ c = 200 symbols and τ p = 5 orthonormal pilot sequences. The large-scale fading coefficients α m and αmk are generated according to the three-slope propagation model in [START_REF] Ngo | Cell-free massive MIMO versus small cells[END_REF]. The large-scale fading coefficient β mk is formulated as β mk = βmk a mk , where βmk is generated by the three-slope propagation model in [START_REF] Ngo | Cell-free massive MIMO versus small cells[END_REF]. The binary variables a mk accounts for the probability that the direct links are unblocked, and it is defined as a mk = 1 with probability p. Otherwise, a mk = 0 with probability 1 -p, where p ∈ [0, 1] is the probability that the direct link is not blocked.The covariance matrices are generated according to the spatial correlation model in [START_REF] Ngo | Cell-free massive MIMO versus small cells[END_REF] with d H = d V = λ/4. The pilot power is 100 mW and ν = 1. The power control coefficients are η k = 1, ∀k. Without loss of generality, in particular, the N phase shifts in Φ Φ Φ are all set equal to π/4, except in Fig. 2(c) with different phase shifts. Three system configurations are considered for comparison: i) RIS-Assisted Cell-Free Massive MIMO: This is the proposed system model where the direct links are unblocked with probability p. It is denoted by "RIS-CellFree". ii) Conventional Cell-Free Massive MIMO: This is the same as the previous model with the only exception that the RIS is not deployed. It is denoted by "CellFree". iii) Cell-Free Massive MIMO without the direct links: This is the worst case study in which the direct links are blocked with unit probability and the uplink transmission is ensured only through the RIS. This setup is denoted by "RIS-CellFree-NoLOS". In Fig. 2(a), we illustrate the sum net throughput as a function of the probability p. In particular, the average sum net throughput is defined as K k=1 E{R k }. Cell-Free Massive MIMO provides the worst performance if the blocking probability is large (p is small). If the direct links are unreliable, as expected, the net throughput offered by Cell-Free Massive MIMO tends to zero if p → 0. In addition, the proposed RIS-assisted Cell-Free Massive MIMO setup offers the best net throughput, since it can overcome the unreliability of the direct links. An RIS is particularly useful if p is small since in this case the direct links are not able to support a high throughput. Fig. 2(b) compares the three considered systems in terms of sum net throughput (defined as K k=1 R k ) when p = 0.2. We observe the net advantage of the proposed RIS-assisted Cell-Free Massive MIMO system. The worstcase RIS-assisted Cell-Free Massive MIMO system setup (i.e., p = 0) outperforms the Cell-Free Massive MIMO setup in the absence of an RIS. Fig. 2(c) focuses on the RISasissted Cell-Free Massive MIMO setup, since it provides the best performance. We compare the sum net throughput as a function of the phase shifts of the RISs (random and uniform phase shifts according to Corollary 1 in the presence of spatially-correlated and spatially-independent fading channels according to [START_REF] Ngo | Cell-free massive MIMO versus small cells[END_REF]. In the presence of spatial correlation, we consider R m = α m d H d V I N and R k = αk d H d V I N , ∀m, k. If the spatial correlation is not considered, there is no significant difference between the random and uniform phase shifts setup. In the presence of spatial correlation, on the other hand, the proposed uniform phase shift design, which is obtained from Corollary 1, provides a much better throughput. This result highlights the relevance of using even simple optimization designs for RIS-assisted communications in the presence of spatial correlation.

V. CONCLUSION

We have considered an RIS-assisted Cell-Free Massive MIMO system and have introduced an efficient channel estimation scheme to overcome the high channel estimation overhead. An optimal design for the phase shifts of the RIS that minimizes the channel estimation error has been introduced and has been used for system analysis. Also, a closed-form expression of the ergodic net throughput for the uplink data transmission phase has been proposed. Based on them, the performance of RIS-assisted Cell-Free Massive MIMO has been analyzed as a function of the fading spatial correlation and the blocking probability of the direct AP-user links. The numerical results have shown that the presence of an RIS is very useful if the AP-user links are mostly unreliable with high probability.

Fig. 1 .

 1 Fig. 1. An RIS-assisted Cell-Free Massive MIMO system where M APs collaborate with each other to serve K distant users.

Fig. 2 .

 2 Fig. 2. The sum net throughput [Mbps]: (a) Average sum net throughput versus the unblocked probability of the direct links; (b) CDF of the sum net throughput with the unblocked probability of the direct links p = 0.2; (c) Sum net throughput as a function of the phase shift setup (equal or random) and the unblocked probability of the direct links p = 0.2.

Two matrices A and B of size N × N are similar if there exists an invertible N × N matrix U such that B = U -1 AU.