Guillaume ; Gillet

Olivier Vitrac
email: olivier.vitrac@agroparistech.fr

Delphine Tissier

Philippe Saillard

Stephane Desobry

Development of decision tools to

Keywords: Packaging -food simulants food contact materials, migration modeling, compliance testing, sanitary survey. 31 32

Introduction 33

For the last three decades, consumer exposure to chemicals from food packaging materials 34 has attracted much public attention and interest of European regulatory authorities. Thus,35 although the EU regulation is still in a consolidation phase, 502 substances (including 230 36 monomers and 272 additives) among the 937 which are positively listed in EU directives on 37 plastics in contact with food (European Commission, 2008), are subjected to specific 38 migration limits (SML) due to toxicological concern. The recent authorization of substances 39 not belonging to the positive list, when they are located behind one or more layers, the so-40 called functional barrier, which prevents their migration into foods or food simulants above a 41 detectable level (European Commission, 2007), complicated the assessment of compliance 42 testing and requirements of traceability (European Commission, 2004). The introduction of 43 modeling in directive 2002/72/EC (European Commission, 2002) as an alternative to time 44 consuming migration testing was an appropriate response to the increasing sophistication of 45 plastics materials on the market including multilayer materials, active or intelligent 46 packaging, recycled polymers, nanocomposites… Despite significant collaborative work on 47 on specific experiments has been applied on finished food products on the market (Vitrac et 73 al., 2007) and on consumer exposure [START_REF] Vitrac | Design of Safe Packaging Materials Under Uncertainty[END_REF] only in rare occasions where a 74 significant expertise on physicochemical properties of packaging materials was available. 75 76 This paper proposes to revise the paradigm of migration modeling by introducing modeling 77 along with deformulation experiments and refined prediction scenarios, recently made 78 available. The main objective was to combine several rapid experimental and numerical 79 techniques to improve the expected utility value up to 1:3 (Table I), in order to justify the 80 additional costs induced by modeling (software licenses, training…). The final approach will 81 be integrated in an expert system developed by INRA including databases and simulation 82 tools. In this new framework, the issue is not to minimize the deviation between the predicted 83 value and the real one, but determining the fastest decision path among possible alternatives 84 to demonstrate the compliance in the first place. The demonstration relies either on a 85 predicted value lower than SML or to a maximum initial concentration value in the polymer 86 (MIC), acknowledged by the supplier to be higher than the true amount (Figure 1). The 87 proposed four steps algorithm is particularly comprehensive for monolayer materials as few 88 input parameters, denoted x i , and including initial concentrations, partition coefficients and 89 diffusion coefficients, are involved: 90 -Step 1: finding the most influencing set of inputs, X = [x 1 x 2 …x n ] for the considered 91 conditions of migration (food packaging geometry, contact time and temperature…). 92 -Step 2: finding the optimal set x°, which leads to migration = SML (idem with MIC). 93 -Step 3: for all inputs in the set, finding all conditions "x i <x i °" or "x i >x i °" to guarantee 94 migration<SML. 95 -Step 4: check whether all conditions are likely according to values available in 96 databases, additional modeling or to experiments. 97 The first step is the most critical because the number of influencing inputs, n, and their 98 identities are not a priori known. Our objective was to analyze all four steps for polyolefins in 99 contact with 50% ethanol as simulant of dairy products (European Commission, 2007). This 100 situation was of significant concern as polyolefin materials represent nearly one third of all 101 food packaging applications and currently no specific guidance exists to predict the migration 102 in such a food simulant. Additional interest in such semi-crystalline materials was triggered 103 by the recent validation of several alternatives to conventional migration testing. 104

Overestimates (Begley et al., 2005) and likely estimates (Vitrac et al., 2006) of diffusion 105 coefficients have been published for all substances, which do not cause blooming effects. The 106 French national project ACTIA RA 0522 "decision tools for compliance testing" (ACTIA, 107 2005) contributed to generate missing steps in the predictive approach of the migration from 108 polyolefins starting from the diagram depicted in Figure 1: i) rapid deformulation of materials 109 from Fourier Transform Infrared spectroscopy (FTIR) absorption spectra of polymer extracts 110 (Gillet et al., 2009a), ii) prediction of partitioning between polyolefins and most simulants 111 (alcohols, water and water-ethanol mixtures) from molecular simulations of molecular pair 112 interactions [START_REF] Vitrac | Prediction of partition coefficients between food simulants and 585 packaging materials[END_REF]Gillet et al., 2009b and2009c). The importance of 113 deformulation to identify possible migrants and their extent has been already discussed in the 114 conclusions of the EU AIR research program CT94-1025 [START_REF] Feigenbaum | Safety and quality of food contact 552 materials. Part 1: Evaluation of analytical strategies to introduce migration testing into good 553 manufacturing practice[END_REF]. For 115 polymer additives comprising several repeated chemical patterns, e.g. di-tert-butylphenol in 116 hindered phenolic antioxidants, FTIR offers additional worst-case scenarios based on an 117 equivalent amount patterns, considered as high diffusing individual molecules. The direct 118 calculation of activity coefficients at atomistic scale for both flexible migrants and hindered 119 migrants in various food simulants and in polyolefins widens dramatically the range of 120 applicability of modeling for monolayer and multilayer materials (Gillet et al., 2009b;Gillet 121 et al., 2009c;[START_REF] Vitrac | Prediction of partition coefficients between food simulants and 585 packaging materials[END_REF][START_REF] Vitrac | An off-lattice Flory-Huggins approach of the partitioning of bulky 588 solutes between polymers and interacting liquids[END_REF]. As the approach was validated 122 in more than 51 different conditions and did not require any fitting procedure, it would 123 replace advantageously both rough estimations suggested in the "practical guide" to European 124 Directives (European Commission, 2003) and the semi-empirical group contribution devised 125 for the partitioning of volatile substances [START_REF] Baner | Prediction of solute partition coefficients between polyolefins 520 and alcohols using the regular solution theory and group contribution methods[END_REF]. Besides, direct 126 calculations provided two important insights at molecular scale, which were underestimated 127 in previous works (Gillet et al., 2009b). On one hand, additives comprising a repeated pattern 128 have a higher chemical affinity than their constitutive patterns for simulants. This entropic 129 effect, which does not depend on chemical composition effects, was related to the increased 130 number of micro-configurations to recreate simulant-simulant interactions (e.g. hydrogen 131 bonding) when a large guest molecule was inserted among smaller host ones. The practical 132 consequence is that the partitioning of large additives cannot be extrapolated from volatile 133 substances. On the other hand, strong deviations from ideality in water-ethanol mixtures were 134 found to be responsible for highly non-linear variations of partitioning with the volume 135 fraction in ethanol, and where in contradiction from prior extrapolations in pure water. The 136 practical consequence is that the lowest chemical affinity for the simulant is expected to be 137 obtained with 20%-30% ethanol rather than in pure water Gillet et al. (2009c). 138 139The paper is organized as follows. The first section details the methodologies, 140 assumptions and scenarios to assess the compliance according to the decision tree depicted in 141 Figure 1 and the proposed four steps algorithm. As a quantification of the initial amount may 142 be required, a fast extraction technique under pressure was compared with the conventional 143 reflux method on a typical high density polyethylene (HDPE). Eight typical migrants of 144 polyolefins were particularly considered. The photo-initiator used in printing inks, 2-145 isopropylthioxanthone (2-ITX), was included because it was recently found in several food 146 products including dairy ones in spite of not being authorized in EU [START_REF] Bagnati | Direct analysis of 517 isopropylthioxanthone (ITX) in milk by high-performance liquid chromatography/tandem 518 mass spectrometry[END_REF] et al., 2007;[START_REF] Rothenbacher | 2-isopropylthioxanthone (2-ITX) in food and 577 food packaging materials on the German market[END_REF]. The main experimental and simulated results are 148 presented in the second section. Biases due to possible incomplete extractions of semi-149 crystalline polymers, which are insoluble in solvents, were tabulated according to the number 150 of extraction steps as suggested by [START_REF] Choller | Safety and quality 530 of plastic food contact materials. Optimization of extraction time and extraction yield, based 531 on arithmetic rules derived from mathematical description of diffusion. Application to control 532 strategies[END_REF]. The effect of combined sources of 151 uncertainty was analyzed by ranking all prediction scenarios according to their ability to 152 demonstrate the compliance and according to the number of uncertain inputs. The main 153 interest of the approach is that the typical generated diagrams can be used to extrapolate the 154 effect of uncertainty for migrants different from those studied by simple geometric 155 constructions. The findings and provides recommendations for its generalization in good 156 manufacturing practices and quality auditing purposes are summarized in the last section. 157 158

Materials and methods 159

Tested substances 160

Tested substances are listed in Table II. They were chosen to be representative of most 161 technological functions in polyolefins. Irganox 1330 and Irganox 3114 were proposed as 162 typical phenolic hindered antioxidants including a repeated chemical pattern. As they include 163 three BHT patterns, they can be easily mistaken for dispersed BHT molecules in FTIR 164 spectroscopy, with an equivalent concentration assessed three times higher than the true one. 165

Erucamide is an unsaturated long chain carboxylic acid amide used as slip agent or lubricant 166 in plastic films. At concentrations above 10 3 mg•kg -1 in polyolefins, it is only partly 167 compatible with polyethylene (Gillet et al., 2009c), its migration does not obey to Fickian 168 diffusion and migration modeling cannot be applied. As a result, it was only used in extraction experiments as a prototype of substance poorly dispersed in polyolefins. Two 170 typical low molecular weights surrogates, which are not intentionally added substances into 171 plastics, were also included: 2-ITX and diphenylmethane. The latter was introduced as 172 partitioning data have already been published by [START_REF] Baner | Prediction of solute partition coefficients between polyolefins 520 and alcohols using the regular solution theory and group contribution methods[END_REF]

Migration modeling 199

Migration modeling was performed under the same assumptions described by Vitrac and 200 Hayert (2006) and by Vitrac et al. (2007) for similar applications. Migration in food or in a 201 food simulant was assumed to be one-dimensional without any reaction and mass losses with 202 surroundings. The material was considered to be homogeneous and the relaxation of polymer 203 chains was assumed unaffected by the desorption process. The concentration profile in 204 substance i, ( )

, i
C Fo x , along the packaging thickness, with SI units in kg•m -3 , evolved 205 therefore with the dimensionless diffusion time, denoted Fo , as: 206 and x is the dimensionless position between 0 and 1 within the film. 208 D i,P is the diffusion coefficient with SI units in m 2 •s -1 ; P l is the packaging thickness and t is 209 the contact time. Only a pervious contact is assumed on food side, at x=1: 210

( ) ( ) 2 2 , , i i C Fo x C Fo x Fo x ∂ ∂ = ∂ ∂ (
( ) ( ) , / 1, 1, i F i F P i F Fo x Fo C C Bi K C x Fo C Fo Fo x = ∂ ∂   = - = ⋅ ⋅ = -   ∂ ∂ (2) 211
where Bi is the mass Biot number assigned to 10 3 to prevent a significant mass transfer 212 resistance at the interface, as the liquid simulant was stirred. 213

( ) ( ) , / , , 1, i F P i F i P K C Fo C x Fo = → ∞ = → ∞
is the partition coefficient between the food product or 214 food simulant, denoted F, and the packaging material, denoted P. The concentration in food 215 or in food simulant, ( ) F C Fo , was obtained from the mass balance t=0 and t by assuming a 216 uniform initial concentration profile equal to 0

, i P C : 217 ( ) ( ) 1 0 / , 0 , F P F i P i L C Fo C C x Fo dx ⋅ = - ⋅ ∫ (3) 218 where ( ) / F P F p L V A l = ⋅
is the volume dilution ratio. V F and A are respectively the volume of 219 III for each substance and in Table IV for each scenario. A similar , i P D 229 value was assigned to Irganox 1330 and Irganox 3114 as they had a very similar chemical 230 structure and as the expected difference is lower than the common experimental error as 231 assessed in (Vitrac et al. 2006).

, / i F P K values were calculated in the framework of a 232 generalized Flory-Huggins theory at atomistic scale as discussed in Vitrac and Gillet (2008) 233 and in Gillet et al. (2009b and2009c). Desorption was assumed to occur in water-ethanol 234 (50:50 v/v) to simulate a contact with a dairy product (European Commission 2007). As a 235 possible worst case of packaging materials in contact with a dairy product, additives were 236 assumed to be used in a rigid container of 200 µm thickness and in contact with a low volume 237 of food of 100 mL. The conditions of test corresponded to a contact of 10 days at 40°C. By 238 contrast, low molecular weight contaminants were assumed to be located in thin 50 µm liner 239 in contact with a product of 1 L. To simulate more realistic conditions in the latter case, a 240 contact with a dairy product at 4°C during 90 days was assumed. 241

The scenarios numbered from 1 to 8 were located at the vertices of a cube (Figure 2). The 242 face 1,4,7,3 corresponded to uncertainties on transport properties. The opposite face 2,6,8,5 243 combined in addition uncertainty on composition, 0 , i P C . In absence of a transport property, a 244 robust overestimate was used instead. As no recommendation has been specifically proposed 245 for 50% ethanol used as simulant, an upper bound of 10 3 was chosen for , / i F P K + instead of the 246 unitary value recommended for simulant D (olive oil) in the "Practical Guide" (European 247

Commission 2003) as extreme worst case. This choice was consistent with experimental , / i F P K 248 values found higher than 1 between polyolefins and ethanol (Vitrac et al. 2007, Gillet et al. 249 2009c). As detailed in table II, in absence of data on initial concentrations for migrants used 250 Scenarios 9 and 10 were special cases where the identity of the migrant i was unknown. In 255 this work, they were applied to phenolic antioxidants similar to those previously studied 256 (Irganox 1330, Irganox 3114 and BHT) but for which only an estimate of the concentration in 257 di-tert-butylphenol (BHT pattern) was assumed to be known, as detailed in Table IV. This 258 situation was consistent with a concentration determination by FTIR (Gillet et al., 2009a). 259 260

Results and discussion 261

Extraction yields 262 A list of substances, which may migrate from the contact material to the food, is 263 conventionally required before starting any assessment of specific migration limits. This 264 approach is well suited for substances, which are intentionally added to materials in contact 265 with food. This information is obtained either by partial disclosure of confidential information 266 or by a deformulation step. For practical cases taken in consideration by the Laboratoire 267 National d'Essais in France, the formulation is transmitted for less than 5 % of tested samples 268 (Table I). As a result, deformulation appears currently as the most robust approach to start the 269 decision tree depicted in Figure 1. To be cost-efficient, this additional step must be rapid and 270 predictive. The difficulty arises with semi-crystalline materials such as polyolefins, which 271 cannot be dissolved in most of solvents below their melting point. Figure 3 C →∞ , given by the overall mass balance defined in Equation 4, were also 297 highly different: 298

( ) 0 , , / , / 1 1 Fo i F i P F P i F P C C L K →∞ = ⋅ + (4) 299
As a same dilution ratio, L F/P , and a same initial concentration, 0 , i P C , were applied to all 300 migrants, the differences in calculated migration rates were only related to partitioning 301 between HDPE and 50% ethanol, K i,F/P used as input. Despite all migrants consisted in the 302 repetition of same BHT pattern, their chemical affinity for 50% ethanol varied significantly 303 with their size. This effect has been extensively studied in Gillet et al. (2009b and2009c) and 304 has been proposed to explain the relative high chemical affinity of large additives for mixtures 305 consisting in much smaller molecules. Similar considerations are detailed by Fill and 306 Bromberg (2003). Accordingly, the effect of solute size on K i,F/P was related to a much higher 307 dispersion of chemical interactions between the migrant and the simulant when the migrant 308 was larger. Indeed, a large molecule as Irganox 3114 consisting in 3 BHT patterns exposed a 309 much lower specific surface to simulant than a single BHT molecule and consequently 310 perturbed less the cooperative network of hydrogen bonds in water-ethanol mixtures. It is 311 underlined that the current approach neglects possible reactions, which could reduce the size 312 In detail, the migration of large additives, Irganox 1330 and Irganox 3114, was controlled 316 by both the migrant chemical affinity for 50% ethanol (which favors the migration 317 comparatively to BHT) and by the limiting diffusion in the film. . The uncertainty effects 318 were however not symmetric. As illustrated in Figure 4b, concentrations of Irganox 1330 and 319 3114 in simulant were greater than 40 % of their equilibrium values, e.g. maximum possible 320 concentration in simulant. Overestimating diffusion coefficients would therefore increase the 321 concentration in food simulation at a maximum by a factor 1/0.4=2.5. In contrast, 322 overestimating the partitioning by several orders of magnitude might lead to a complete 323 extraction, that is 1/50 the initial concentration, that is overestimating the expected migration 324 by one order of magnitude (Figure 4a). 325

It is worth noting that such a very low partitioning (Table III) with the simulant requires a 326 particular care in the writing of the boundary condition (Equation 2) and in the resolution of 327 the complete set of transport equations. It is expected that the requested time to reach the 328 equilibrium is lower when the chemical affinity is very low for F. Simply explained, the 329 equilibrium was reached all the more rapidly when the amount of migrant to be desorbed to 330 reach the equilibrium was lower. In the absence of chemical affinity of the simulant, the 331 macroscopic equilibrium would be thus obtained instantaneously as soon as a single molecule 332 would have left the polymer. These features, which may appear counter-intuitive, are 333 discussed by [START_REF] Vitrac | Identification of diffusion transport properties from 592 desorption/sorption kinetics: an analysis based on a new approximation of Fick equation 593 during solid-liquid contact[END_REF]. 334

The practical consequence was that the mis-identification in FTIR of Irganox 1330 as 335 BHT, with a three times higher molar concentration (note that the mass concentration would 336 remain unchanged), would not be necessarily a worst case as the affinity of Irganox 1330 for 337 the simulant would be also higher. Figure 5, the values corresponding to scenarios numbered from 1 to 8 (see Table IV) are 348 plotted at the vertices of a rectangle. The extents of variations represent the effects of 349 overestimations when one or two input values were missing and when strong overestimates 350 were used instead (see Table III). 351

For a small migrant as BHT, equilibrium was reached in less than two days so that only 352 L F/P introduced only some additional effects. They were related to two complementary 362 contributions: i) lowering the effect of a poor chemical affinity for simulant at high L F/P 363 values and ii) increasing the dilution of the amount of migrated substances. Depending on the 364 units used to express the specific migration limit (SML), in mg•dm -2 of packaging or mg•kg -1 365 of packaging, the one or the other led to a more conservative assumption. As we expressed 366 SML in mg•kg -1 , L F/P =50 was chosen in the remainder of the work as a worst case 367 assumption. The conventional assumption 1 kg in contact with 6 dm 2 would lead to a 368 concentration in food simulant down to twice lower. 369 370

The effect of cumulating up to three sources of uncertainty made the utility of some 371 scenarios questionable according to the expected migration limit of concern (SML or 372 detection limit). To generate a practical sensitivity diagram depicting all combined sources of 373 uncertainty, a simplified descriptor, derived from Equation 2, 10 , / i F P Fo K ⋅ was used instead of 374 isocontours. For low molecular weight contaminants with a contact time of 90 days ( values (3,5,7,8), as generally recommended for migration 383 modeling, were able to demonstrate the compliance only when the migration limit was as high 384 as the overall migration limit in EU regulation, 60 mg•kg -1 , that is for Irganox 1330. By 385 contrast, refined scenarios 1,2,4,6 involving at least one likely transport properties were 386 always applicable to demonstrate compliance. Irganox 3114 and DPM were special cases, 387 where compliance could not be determined when more than one source of uncertainty was 388 included. 389 390

Decision diagrams in MIC 391

The usability of previous sensitivity diagrams by end-users, having not access to the initial 392 formulation of their materials, was improved by replacing the concentration in food simulant 393 by the equivalent maximum concentration in the packaging material, MIC, which matches the 394 limit of concern. It was calculated as: 395

( ) 0 , , limit of concern i P Fo i F C MIC C = ⋅
(5) 396 MIC values are more practical as they can be calculated a priori for a given food 397 packaging application and they can be compared subsequently with formulation data derived 398 from general manufacturing practices or from a partial formulation disclosure (Figure 1). The 399 interesting feature is that the provider needs only to acknowledge that it uses the substance at 400 a lower concentration than MIC. To make robust conclusions, the uncertainty associated to 401 each MIC estimates must be accounted. MIC values corresponding to Figure 6 are plotted in 402 Figure 7 in increasing order. When implemented as decision diagram, each bin represents a 403 condition node: is the initial concentration higher than the value above each bin? If the answer 404 is yes, the question must be repeated with the next bin on the right. If the answer is no, the 405 compliance is demonstrated. For large additives, the MIC values were progressively 406 increasing as the migration model was sensitive to all sources of uncertainty. For low 407 molecular weight contaminants and additives, changes between values were most abrupt. 408

Ordering scenarios a priori guaranteed the decision-maker that acquiring new information 409 will improve the decision. A strict improvement is obtained by merging scenarios with similar 410 values. As it is physically expected, , / i F P

K

was the most valuable information for small 411 migrants while it was , i P D for large additives including a repeated pattern. 412 413

Decision diagrams for phenolic antioxidants mistaken with BHT 414

When the identity of the substance was not known as it might occur in FTIR 415 deformulation, a large molecule (in this work a large hindered phenolic antioxidant) was 416 replaced by much smaller one (BHT). This approximation suited particularly for Irganox 417 1330 and Irganox 3114. Nevertheless two additional difficulties appeared: the initial 418 concentration was expressed in equivalent concentration in BHT patterns and the limit of 419 concern was assumed to be the SML of BHT. The reliability of such assumptions to 420 demonstrate the compliance is depicted in Figure 8. As the SML of BHT was much lower 421 than the one of real additives, the approximation guaranteed a conservative decision limit. By 422 introducing the expected transport properties of BHT, From the point of view of decision-maker, the dependence diagram depicted in Figure 2 434 should be started in the reverse order with the lowest knowledge level. For antioxidants, it 435 could be started with scenarios 9 and 10 as a homologous molecule was available. For other 436 molecules, it should be started with the scenario 8. The best scenario was the one, which 437 made it possible to demonstrate the compliance and was connected with shortest path length 438 (decision path length) to the initial scenario. In other words, the optimal scenario should 439 minimize the extent of additional information necessary to demonstrate the compliance. 440 441 Previous figures 5-6 cannot be used directly to determine the shortest path length as they 442 focus on the values of physicochemical properties instead of on the amount of available 443 information. In practice, the lengths of all possible paths (8) were read in Figure 2, they did 444 not depend on the considered substance or on the considered problem. The efficiency of each 445 corresponding scenario (ability to demonstrate the compliance) depends both on the 446 considered case, e.g. substance, and on the decision limit value. In order to derive almost self-447 similar diagrams, i.e. reusable for other applications, where the shortest distance can be easily 448 read, all scenarios were ranked according to their equivalent decision path length and plotted 449 in Figure 9 as concentrations normalized according to the considered SML. Starting from the 450 left, the shortest path length was consequently defined as the first abscissa with a normalized 451 concentration lower than 1. Except when the SML was as high as the global migration limit 452 (Irganox 1330), a path of length 4 was required, which indicated that no more than two 453 sources of uncertainty could be combined. On one hand, it could be argued, as a rule of 454 thumb, that when partial formulation disclosure is not accepted or available, at least one 455 transport property must be known with enough accuracy. On the other hand, when the initial 456 concentration is available, conventional overestimates, although inaccurate, are expected to be 457 able to demonstrate the compliance not in all cases but with a significant expectation, between 458 50% and 75%. 459

For polyolefins (one third of total food contact applications), the current state of the art 460 suggests that scenario 2 and possibly scenario 1 should be nearly available for all practical 461 cases, while minimizing the risk of false negative (considered as "not able to demonstrate the 462 compliance" when it is compliant). In addition, the risk of false positives would be very low 463 while appropriate safety margin on likely transport properties are included in the decision. 464 465

Conclusions and future prospects 466

This paper addresses a general framework to combine migration modeling and experiments 467 into the same decision tool to demonstrate the compliance of thermoplastics materials in 468 contact with food with regulatory limits. The objective was to select the necessary 469 experiments and simulation scenarios among several alternatives to make a decision: "it is 470 compliant" or "the compliance cannot be demonstrated with this approach". In comparison to 471 conventional approaches used in migration modeling, the ambition was neither to assess the 472 real migration nor to overestimate all uncontrolled factors but to dramatically increase the 473 designing an end user-oriented approach, which was: 475 -more efficient (more robust, it takes less time to accomplish the demonstration, it requires 476 fewer assumptions, formulation disclosure is not systematically required…); 477 -easier to learn (based on simple diagrams, include a robust physical reasoning based in 478 particular on minimal number of dimensionless quantities); 479 -more satisfying to use (the method toggles when required from empirical models of 480 transport properties to more sophisticated ones based on a molecular simulation). 481

The method was tested in complicated cases involving the migration of low molecular 482 weight contaminants and hindered phenolic antioxidants into the simulant proposed for dairy 483 products, 50% ethanol, for which the possibilities of modeling had not yet been examined. 484

Although the current results might find an application in the previous 2-ITX crisis, the finality 485 was mainly methodological. All ten scenarios incorporating up to three sources of uncertainty 486 has been incorporated into an on line semi-supervised software on the SAFE FOOD 487 PACKAGING PORTAL (INRA 2008). It generates all proposed diagrams in less than 3 488 minutes so that the overall decision tree can be easily integrated along with the process of 489 identification and quantification of substances and along with the audit of the provider. 490 491 Apart from the decision tool concept, this work brought two important insights for both 492 compliance testing and future sanitary survey. The migration of large additives such 493 antioxidants is mainly controlled by the resistance of material to diffusion. This low effect of 494 partitioning is exacerbated by their higher chemical affinity for food simulant consisting in 495 much smaller molecules. As an example of consequence, long time contacts and heating 496 should lead to higher contamination of food products by large additives than by smaller ones 497 even if they had a similar chemical composition. As this description neglects a possible 498 -2 ‡ † experimental partition coefficient between HDPE and ethanol is available [START_REF] Baner | Prediction of solute partition coefficients between polyolefins 520 and alcohols using the regular solution theory and group contribution methods[END_REF]. IV. Combination of input values associated to each scenario. The index i in scenarios 1..8 corresponds to an identified substance, whereas 630 the index h in scenarios 9..10 corresponds to an homologous molecule (BHT for antioxidants) with a similar mid-infrared vibrational spectrum 631 (Gillet et al., 2009a). 632 with one. 666

Scenario 0 , i P C value , / i F P K value , i P D value 1 0 , i P C , / i F P K , i P D 2 ( ) 0 , i P C + , / i F P K , i P D 3 0 , i P C , / i F P K + , i P D 4 0 , i P C , / i F P K , i P D + 5 ( ) 0 , i P C + , / i F P K + , i P D 6 ( ) 0 , i P C + , / i F P K , i P D + 7 0 , i P C , / i F P K + , i P D + 8 ( ) 0 , i P C + , / i F P K + , i P D +
The corresponding limit for the real additives is also depicted. 667

F

  and the contact surface area between P and F. Equations 1-3 were integrated in time using a 220 finite volume technique documented and implemented online in QSPR-MS tools of the Safe 221 Food Packaging Portal (INRA 2008). All simulations results were exported as datasheets to Ten scenarios with different sources of uncertainty were considered for five additives as 226 plotted in Figure 2. For each quantity X required for modeling, two values were considered: a 227 likely value X and a robust overestimate X + such that X X + < . The corresponding values are 228 gathered in Table

  Irganox 1330, Irganox 3114, BHT), maximum values recommended by the 251 providers were used. For low molecular weight contaminants (2-ITX and diphenylmethane), 252 concentrations in agreement with the extensive study of[START_REF] Rothenbacher | 2-isopropylthioxanthone (2-ITX) in food and 577 food packaging materials on the German market[END_REF] were 253 used. 254

  Figure 4 plots the dimensionless kinetics corresponding to the most complete scenario (i.e. 289

  of migrants and which could modify accordingly the chemical affinity of fragments for the

  Introducing simultaneously a significant amount of uncertainty for the diffusion coefficient 341 and partition coefficient, as it could be in practice, complicates arriving at clear conclusions 342 regarding the possible exceedance of a given threshold. The effects of combined variations in two dilution ratios: L F/P =50 and L F/P =83. These ratios matched the 345 conventional assumption used in conservative consumer exposure assessment in EU 346 (European Commission 2002): 1 kg on food in contact with 6 dm 2 of packaging material. In 347

  uncertainty on partitioning was able to change the concentration in food simulant after 10 353 days. By choosing , / i F P K + (scenarios 3,5,7,8), the extraction was almost complete and the 354 concentration in simulant was solely controlled by the maximum permitted concentration in 355 HDPE. The overestimation factor reached up two decades. The behavior was opposite for 356 large antioxidants and their migration rates were almost independent of partitioning. The 357 isocontours appeared therefore mainly parallel to , / i F P K axis with a maximum overestimation 358 factor of 5. According to Equation4, it is worth to notice that decreasing , / worst case recommended in "practical guide" (2003)) did not change previous 360 conclusions as L F/P values higher than 10 led to almost complete extraction for any , /

  equivalent 423 concentration in BHT pattern, the predictions of scenarios 9 and 10 were close to those of 424 2 respectively (Figure5). As a result, the compliance could be demonstrated 425 easily with scenario 9 and 10. 426The proposed approximation seemed valid when BHT replaces molecules based on the 427 repetition of several BHT patterns. It should not be generalized to any series of additives. Figure8. As it led to equivalent MIC values lower than 200, it was not suitable for 430 antioxidants, which are conventionally used at concentrations much higher. 431 432 Shortest paths to demonstrate the compliance 433

  -based compliance testing. The improvements consisted mainly in 474

  , a dedicated experimental verification is encouraged. Finally, the 499 results highlighted that the crude assumptions recommended in the current "PracticalGuide" 500 (2003) for modeling, based mainly on simulant D and strong overestimates of diffusion 501 coefficients, failed to deserve the demonstration of compliance for intentionally added 502 substances in most of tested cased. As modeling combined with fast deformulation techniques 503 appear to be the mainly independent approach to audit the compliance, further extension of 504 this work to other polymers and multilayer materials is highly desirable. 505 506 Acknowledgment 507We would like to thank the Association de Coordination Technique pour l'Industrie 508 Agroalimentaire and the Association Nationale pour la Recherche Technique for their 509 financial support. We also want to thank the Institut National des Sciences Appliquées from 510 Lyon, and particularly Hervé Perrier-Camby, for processing the studied polymer materials.

607

  ‡The values were derived from blind samples tested during the project RA05

  614 ‡ threshold value corresponding to the migration of an authorized below the detection level through a functional barrier as described in article 10 of directive 2007/19/EC 615 (EC, 2007).

Figure 1 :Figure 2 .Figure 3 .Figure 4 .FigureFigure 6 .PageFigure 7 .Figure 8 .

 1234678 Figure 1: Decision diagram to demonstrate the compliance of plastic materials according to 635

668Figure 9 .

 9 Figure 9. Decision diagrams expressed as minimum path lengths in the influence diagram 669

  Fig.2

  Fig.3

  Fig.6

Table 375 III

 375

	376	), a similar descriptor was proposed: 90 , / i F P Fo K ⋅	. Values of ( , Fo Fo Fo 10 90 , i F C =	)	against the previous
	377	criteria are plotted in Figure 6 for 0 , i P C	and ( ) 0 , i P C	+ . Confidence intervals caused by
	378	uncertainties on transport properties,	K	, / i F P	and , i P D , are also presented to highlight that the

Table I .

 I Transition and overall probabilities of the decision tree plotted in Figure 1. The final probability to demonstrate the compliance is 603 denoted p compliance . The values obtained with the current methodologies used at LNE (Laboratoire National de métrologie et d'Essais) and the 604 expected values with the new proposed paradigm are compared. 605 606 were in particular inferred from the number and type of compliance requests received by LNE during the last five years.

			Current state	Expected value with the new
		Notations	at LNE † † † †	methodology ‡ ‡ ‡ ‡
		a	5%	50%
	Transition probabilities	b c d SML	50% 10% 90%	75% 80% 90%
		d MIC	50%	50%
		p 1 = a	5%	50%
		p 2,SML = a•b	2.5%	37.5%
		p 3,SML = a•b•c	0.25%	30%
	Overall probabilities	p 4,SML = a•b•c•d SML p 2, MIC = a•(1-b) p 3,MIC = a•(1-b)•c	0.20% 2.5% 0.25%	27% 12.5% 10%
		p 4,MIC = a•(1-b)•c•d MIC	0.13%	5%
		p compliance = p 4,SML + p 4,MIC	0.33%	32%
	†The values

Table 2

 2

	). Scenarios 9 and

http://mc.manuscriptcentral.com/tfac Email: fac@tandf.co.ukFood Additives and Contaminants


Bjørn Schmidt

Jens Højslev Petersen

Christian Bender Koch

David Plackett

Nini R Johansen

Vimal Katiyar

Erik Huusfeldt Larsen

Combining Asymmetrical Flow Field-Flow Fractionation with Light Scattering and Inductively Coupled Plasma Mass Spectrometric Detection for Characterization of Nanoclay used in Biopolymer Nanocomposites

Keywords: Food simulants, Migration, Packaging Biopolymer nanocomposites, asymmetrical flow field-flow fractionation, transmission electron microscopy, inductively coupled plasma mass spectrometry, aspect ratio, migration .manuscriptcentral

Introduction

The packaging industry accounts for about 41 % of total global plastic production and out of this 47 % are presently used in food packaging [START_REF] Fomin | Biodegradable polymers, their present state and future prospects[END_REF]. The plastics used for food packaging are at present largely fossil fuel derived. However, polymers obtained from renewable plant resources or biopolymers are now available and are useful in some applications. The sources include fermentation products from agricultural crops and wastes, starch, plant oils, cellulose, gelatine and chitosan. Unfortunately, undesirable biopolymer properties such as brittleness, low thermal stability, and high gas permeability can restrict their practical use. A solution to these problems is to improve the chemical and physical properties of the biopolymer by incorporating nanofillers [START_REF] Fomin | Biodegradable polymers, their present state and future prospects[END_REF]Ray and Bousmina 2005). Layered aluminium silicates (clays) can be used as nanofillers, but other minerals, such as layered double hydroxides (LDH) (Taviot-Gueho and Leroux 2006) may also be used. Evidently, there is a potential for combining green materials with nanotechnology for the development of a new generation of advanced packaging materials. The overall goal of the on-going NanoPack project (www.nanopack.dk) is to develop a technological basis for a cost-efficient production and use of biopolymer nanocomposites produced from renewable resources for use in the food packaging industry. These novel nanocomposites must meet the demands made by consumers and various authorities in terms of functionality, sustainability and safety.

Many successful attempts have been made to develop nanocomposite films, and mostly x-ray diffraction and electron microscopy have been used for characterization of the composite structure (Pluta et al. 2002;Pospisil et al. 2004;Rhim et al. 2006;Zhang et al. 2007). Other methods including dielectric spectroscopy (Pluta et al. 2007) and solid state NMR spectroscopy [START_REF] Bourbigot | Processing and nanodispersion: A quantitative approach for polylactide nanocomposite[END_REF] have also been explored.

To our knowledge no attempt has been made to characterize materials migrating from nanocomposites into liquid food simulants such as ethanol-water mixtures. The most widely investigated clays for nanocomposites are montmorillonites that have been chemically modified by long hydrocarbon side-chain quaternary ammonium compounds (e.g. Cloisite®, Southern Clay Products, Gonzales, Texas). Depending on the modifier chain length, these clays exhibit more or less hydrophobic surface characteristics. Unfortunately, none of the organomodifiers of this type are yet to be found on the positive list in the current EU food contact plastics legislation (EU 2002).

From a food safety point of view it is important to characterize migrates from nanocomposites containing clays as fillers. As a general rule, nanocomposites must comply with the total migration limit of 10 mg dm -2 (EU 2002). In order to further characterize nanoparticles contained in such migrates, asymmetrical flow field-flow fractionation (AF 4 ) is well-suited for determining the size distribution of particles suspended in the food stimulant used for the migration study. In AF 4 an external field force or gradient, rather than partitioning between phases, causes differential retention of nanoparticles according to their hydrodynamic radius [START_REF] Giddings | Field-flow fractionation: analysis of macromolecular, colloidal, and particulate materials[END_REF]. Following the time and size-resolved elution, the analyte particles are carried to one or several detectors. Detection techniques commonly used in conjunction with AF 4 include multi angle light scattering (MALS), or ultraviolet (UV) and fluorescence (FL) spectrometry. In order to acquire data on the sizes (radii) of the separated nanoparticles the MALS detector is particularly useful. This non-destructive detection technique employs a laser beam, which is passed through a flow cell in which the suspended analyte particles cause the incident laser light to be scattered. The intensity of the scattered light is measured simultaneously by several photodiodes placed at different scattering angles in the plane around the cell. By measuring the angular dependence of the scattered light, it is possible to deduce the radius of the particles (Wyatt 1998). The radius of the particles measured by this technique is the so-called root mean square (rms) radius, which is a mass-averaged squared distance of each scattering point in the particle from the particle's center of mass. The theory and application of field-flow fractionation was originally developed by Giddings [START_REF] Giddings | A New Separation Concept Based on a Coupling of Concentration and Flow Nonuniformities[END_REF] and is today widely used with MALS detection for particle characterization in environmental analysis. Furthermore, coupling the non destructive technique AF 4 -MALS to inductively coupled plasma mass spectrometry (ICP-MS) gives an additional detection and quantisation possibility of a wide range of elements contained in the separated nanoparticles such as nanoclay platelets [START_REF] Baalousha | Size fractionation and characterization of natural colloids by flow-field flow fractionation coupled to multi-angle laser light scattering[END_REF][START_REF] Beckett | Use of field-flow fractionation to study pollutant-colloid interactions[END_REF][START_REF] Hassellöv | Determination of continuous size and trace element distribution of colloidal material in natural water by on-line coupling of flow field-flow fractionation with ICPMS[END_REF]Kammer 2005).

The aim of this study was to characterize migrates obtained from PLA nanocomposites with respect to total migration as well as particle size and elemental composition of nanosized matter. For this purpose an analytical platform consisting of AF 4 coupled with MALS and ICP-MS was established and applied. Some of the challenges regarding the correctness of the analytical characterization of particle size were associated with the high aspect ratio of the nanoclay platelets and the wide size range of these platelets in suspension.

Materials and methods

Chemicals

Sodium dodecyl sulphate (SDS, ReagentPlus ≥98.5 %) and sodium azide (SigmaUltra) were purchased from Sigma-Aldrich (Milwaukee, Wisconsin, USA). Water was obtained from Milli-Q Element water purification system (Millipore, Massachusetts, USA). Ethanol (gradient grade for liquid chromatography ≥99.9 %) and 25 % ammonia solution were obtained from Merck (Darmstadt, Germany). Paraffin oil was obtained from Sigma-Aldrich (Milwaukee, Wisconsin, USA). Glass distilled water was used for the migration study. Nitric acid (PlasmaPURE 67-69 %) and hydrochloric acid (PlasmaPURE 34-37 %) from SCP Science (Champlain, New York, USA) as well as hydrogen peroxide (Suprapur 30 %) and hydrofluoric acid (Suprapur 40 %) from Merck (Darmstadt, Germany) were used for dissolution of the clay minerals.

Clays and polylactide

Two types of Cloisite® clays were used in the experiments, Cloisite®Na + and Cloisite®30B.

Cloisite®Na + is natural montmorillonite clay and Cloisite®30B is montmorillonite, which has been organically modified with methyl tallow bis-2-hydroxyethyl quarternary ammonium cations. The PLA used in this work was L-9000, an L-polylactide obtained from Biomer (Krailling, Germany).

Polylactide and polylactide/Cloisite®30B extrusion

In order to produce films for the migration studies, PLA granulate was first dried under vacuum at 40 o C for three hours. Similarly, the Cloisite®30B clay was dried under vacuum at 100 o C for three hours. For practical reasons, PLA and PLA-Cloisite®30B granulates were initially prepared from the starting polymer using a Haake twin-screw extruder with counter-rotating screws. The Cloisite®30B clay were added to the PLA granules at a 5 % w/w loading and physically mixed by mechanical tumbling in a plastic container for five min. In order to ensure adequate adhesion of the clays to the polymer granules, the latter were first mixed with ~0.4 % w/w paraffin oil before clay addition. The twin-screw extruder has four zones and these were typically set to temperatures in the 150-190 o C range with increasing temperature from the in-feed to the out-feed. In each case the extrudate was immediately cooled by passing rapidly through a water bath before pelletising. After drying the pellets at 40 o C under vacuum for a minimum of three hours, the resulting dried pellets were fed into a Haake single screw extruder with five temperature zones set in the range 150-195 o C with increasing temperature from in-feed to out-feed. The extrudate from the single-screw extruder was fed through a film die and then through chill rolls set at 50 o C to generate film with thickness in the range of 100-200 µm and width up to 12 cm.

Characterization of polylactide/Cloisite®30B

The nanocomposite used in this study was examined for intercalation and exfoliation of nanoclay platelets in the PLA matrix by x-ray diffraction (XRD) and transmission electron microscopy (TEM).

X-ray diffraction

As shown in equation 1, Bragg's law can be used to calculate the spacing between planes in the lattice of the clays. In equation 1, n is an integer , λ is the wavelength of the x-rays, d is the spacing between planes and θ is the angle between the incident x-rays and the scattering planes. Nanocomposite films were studied by x-ray diffraction using a Siemens D500 diffractometer (Bruker, Karlsruhe, Germany) equipped with a Co tube and a diffracted beam monochromator.

Transmission electron microscopy

The morphology of PLA/Cloisite®30B nanocomposites was examined by performing bright field imaging using a Tecnai G 2 transmission electron microscope (FEI Company, Hillsboro, Oregon, USA) optimized for this purpose. The microscope was operated at an accelerating voltage of 200 kV. Ultrathin sections of PLA/Cloisite®30B nanocomposites (~ 80 nm thickness) were prepared at room temperature using a Leica EM UC6 ultramicrotome (Leica Microsystems, Wetzlar ,Germany) equipped with a glass knife. The microtomed slices were transferred directly from a dry glass knife to copper grids for examination in the TEM.

Acid digestion-and inductively coupled plasma mass spectrometry of clays

The two clays were digested according to method 3052 of the United States Environmental Protection Agency (United States Environmental Protection Agency 1996) using a mixture of hydrogen peroxide, nitric-, hydrochloric-and hydrofluoric acids. For determination of elements in the acid-digested clays or in the eluting clay fractions from the AF 4 , an Agilent 7500ce ICP-MS system (Agilent Technologies, California, USA), equipped with a Micromist nebulizer, was operated in the non-cell mode. The on-line and real-time coupling of the AF 4 -MALS system with the ICP-MS detector was achieved by a short length of capillary tubing from the outlet of the MALS detector to the inlet of the nebuliser of the ICP-MS. The ICP-MS operating conditions used are given in Table 1.

Instrumentation for separation and size determination of clays

Centrifugation

The system consisted of a centrifuge model 3-18K from Sigma Laborzentrifugen (Osterode, Germany) with an 11180 swing bucket rotor. Sample tubes were 15 ml Sarstedt (Nümbrecht, Germany) tubes.

Asymmetrical flow field-flow fractionation with multi angle light scattering detection

The AF 4 instrument used was an Eclipse 3 system equipped with a Dawn HELEOS MALS detector (Wyatt Technology Europe GmbH, Dernbach, Germany), which was used in combination with an Agilent 1200 HPLC system (Agilent Technologies, Santa Clara, California, USA). The AF 4 and MALS operating conditions are given in Table 1.

Migration study

Studies of migration from the biopolymer were carried out in triplicate using a mixture 95 % ethanol and 5 % glass distilled water mixture as a simulant for fatty foods according to the current EU legislation (EU 2002) and European Standard (CEN 2002;2004). One dm 2 of each of the biopolymer nanocomposite and the pure PLA films were totally immersed into the simulant at 40°C for 10 days. Furthermore, a food simulant blank and food simulant spiked with 2 mg of Cloisite®30B were included in the experiments. One full set of samples were used for gravimetric determination of total migration according to the standardised procedure. Using a second set of samples the migrates were evaporated almost to dryness, re-suspended in 20 ml food simulant and analysed by AF 4 -MALS-ICP-MS without any further sample preparation. Prior to analysis the spiked sample was diluted 1:9 with simulant and the analyses were performed using 500 µl injections. The ICP-MS signal intensity of 26 Mg was used for on-line detection of clay, due to its high abundance in the clay and at the same time low content in the blank food simulant.

Determination of aspect ratio of clay by way of centrifugation

Cloisite®Na + (50 mg) was weighed into to each of nine centrifuge tubes and five ml of water and one drop of 25% ammonia solution were added. The ammonia solution was added to enhance the stability of the suspension by electrostatic repulsion of the negatively charged nanoclay platelets at alkaline pH. The suspensions were ultrasonicated for 15 min and mixed overnight attached to a vertical rotating wheel. The next day, the nine solutions were centrifuged at 4700 rpm for two, five, 10, 15, 20, 25, 50, 100 and 200 minutes respectively. From the top one cm, 750 µl supernatant was sampled and 20 µl were injected into the AF 4 -MALS for analysis.

Results and discussion

Characterization of biopolymer nanocomposite

PLA/Cloisite®30B structure by wide angle x-ray diffraction

The X-ray diffractogram in Figure 1 (curve B) indicates that Cloisite®30B derivatised with the quaternary ammonium compound had a fairly well-ordered intercalated system with a basal spacing between platelets of approximately 1.8 nm (as indicated at 5.7°2θ) and very little low-angle scattering. In contrast, curve A demonstrates that the PLA/clay nanocomposite was highly disordered with a basal spacing of 1.6 nm (as indicated at 6.3°2θ) and an extensive low-angle scatter from exfoliated plates and PLA. In order to acquire representative TEM images of the nanocomposite, nine images of two nanocomposites were recorded and a representative image is shown in Figure 2. Single exfoliated clay platelets can be observed as needle-shaped bodies in the figure, which suggests that the clay was largely exfoliated in the polymer matrix. Some occurring darker areas also indicated intercalation of clay platelets. The information obtained from the TEM-images therefore confirmed the XRD data, namely that the desired exfoliation of clay indeed occurred in the PLA matrix. This information together with the features of the TEM image also indicated that the individual platelets had large aspect ratios.

Clay analysis using inductively coupled plasma mass spectrometry

Mass spectrometric characterization of elemental composition of clays

The ICP-MS spectra of acid-digested Cloisite®Na + shown in Figure 3 A-B were obtained following subtraction of the acid blank from the mass spectrum recorded for the sample. The mass spectrum for Cloisite®30B (data not shown) was similar to that obtained for Cloisite®Na + . The mass spectrum in Figure 3 showed that Cloisite®Na + contained a range of elements that were useful as a multi-variable elemental fingerprint of the clay. In addition to the presence of isotopes of major

Selective on-line detection of clay by inductively coupled mass spectrometry

A fractogram of a nanoclay suspension using the hyphenated AF 4 -MALS-ICP-MS system is presented in Figure 4 showing overlaid intensity profiles corresponding to the 90° light scattering

Migration study

The results obtained from the migration study in Table 2 report values for total migration, particle sizes of migrated matter, and indicate whether or not clay was detected by ICP-MS in the migrates.

Particle sizes of migrated matter were deduced from the cumulative number fraction as depicted in 2) that there was a large difference between the pure PLA and the nanocomposite with respect to the particle size distributions and also with respect to the upper particle size limit in the corresponding migrates. Finally, the results of this migration study showed that the amount of total migration differed between the two samples. Both sample materials however, were in compliance with the total migration limit of 10 mg dm -2 . Analysis of the migrate obtained from the nanocomposite material using the full hyphenated analytical system is shown in Figure 6. The fractogram shows that the migrate indeed contained nanosized particulates as indicated by the light scattering signal and that their sizes ranged from 50-800 nm. The 20 th to 80 th percentiles of the size distribution were 190 nm and 230 nm, respectively, in rms radii as shown in Table 2. The fractogram however, also demonstrates absence of clay as no increase in the time resolved ICP-MS signal was seen for 26 Mg within the elution time window. In order to confirm this finding the migrates were separately submitted to acid digestion and the isotope pattern following ICP-MS analysis indeed showed that clay was absent.

Šimon et al. ( 2008) discussed in theory the possible migration of particles from nanocomposites and concluded that only very small particles with a diameter of about 1 nm could migrate. A prerequisite for this conclusion was that the polymer matrix had a relative low dynamic viscosity and that it did not interact with the nanoparticles. To our knowledge the migration study presented here is the first attempt to characterize migrates from nanocomposites by using size separation and determination by AF 4 -MALS-ICP-MS and our results confirm the theoretical predictions by Šimon et al.

Clay aspect ratio

Information on clay aspect ratio, α, was of importance for the correctness of the determination of nanoparticle radius by light scattering techniques. This was mainly because the algorithms applied to the collected raw light scattering data did not immediately apply to estimation of radius of particles with large aspect ratios. The aspect ratio is given by α = d pl /h where h is the thickness and d pl the diameter of the platelet. To prevent large clay platelets from precipitating in the AF 4 channel, removal of this fraction of the suspensions was necessary prior to running the analyses and was

ρ π η ∆         = (2)
This equation allows calculation of the time t s it takes for a spherical particle with a diameter d to travel from a given surface in a tube (fill-up point r in ) to a chosen distance below that surface (sampling point r out ) at constant force (revolutions per min, rpm). The difference in density of the particles relative to the surrounding medium is given by ∆ρ and η is the viscosity of the medium.

The constant values used in the equation are r out = 12.5 cm, r in = 11.5 cm, ∆ρ = 1.86 g ml -1 , rpm = 4700 and η = 0.000894 Pa s. For non-spherical particles, like nanoclay platelets, equation ( 2) can not be used. The effective sedimentation diameter d s for platelets has been described as [START_REF] Jennings | Particle-Size Measurement -the Equivalent Spherical Diameter[END_REF]:

α α arctan 5 . 1 pl s d d = (3)
where d pl is the diameter and α is the aspect ratio of a platelet. Combination of equations ( 2) and ( 3

        ∆         = α α ρ π η (4)
which allows for the estimation of the equivalent travel time for a platelet with diameter d pl and an aspect ratio α.

The two clays included in this study originated from the same montmorillonite source. The determination of α, which was carried out for Cloisite®Na + because of its immediate compatibility with aqueous solvents, was based on data obtained from centrifugation tests. The apparent rms radii of the nanoclay platelets in the supernatants obtained following variation of centrifugation time were determined by AF 4 -MALS analysis. The 98 th percentile of the cumulative number fraction was used to characterize the upper cut-off radius of the size distribution (Figure 5 A). The use of the 98 th percentile rather than the maximum value of determined radii improved the robustness of this metric by excluding a few large agglomerated particles from the measured size range.

The experimental fit between centrifugation time and 98 th percentile diameter is shown in Figure 7 curve A using 2 times the rms radius from AF 4 -MALS as the equivalent diameter, d rms . This value is used as the best estimate of d pl in equation 4. Using d rms as d pl of the platelets we iterate an α value by estimating t s values which most closely fitted the experimental data. The best fit was obtained by an α value of 320 and can be seen in Figure 7 curve B. The fit obtained using an aspect ratio of 320 however, is valid only for particle sizes (d rms ) above 450 nm. This may be due to a different aspect ratio of the fraction of smaller nanoclay platelets but further iterations of the αvalue did not lead to any successful match between calculated t s values and experimental set t s values for the smaller nanoparticles.

Conclusions

This study demonstrated that AF 4 with MALS detection was useful for characterizing the size of nanoparticles contained in migrates from nanocomposites of PLA and organomodified montmorillonite clay as filler. However, this coupled instrumentation alone did not provide any information on the identity of the nanoparticulates occurring in the food simulant. This limitation was overcome by coupling the AF 4 -MALS to the element-selective ICP-MS detector, which provided additional information on trace elements known to be naturally present in the clay. The hyphenated analytical system was used for selective detection and size characterization of suspensions of the unmodified montmorillonite Cloisite®Na + . Furthermore, the analytical system was applied for characterization of migrates from nanocomposite films of PLA and the organomodified Cloisite®30B montmorrilonite clay used as filler. The results demonstrated that nanoparticulates at 50-800 nm in radius indeed migrated from the nanocomposite but ICP-MS signals corresponding to clay minerals were absent.

A possible drawback to an accurate size determination of nanoclay platelets by the MALS detector is their large aspect ratio. The algorithms of the manufacturer's software for radius calculations of nanoparticles were in theory valid only for specific particle shapes and not for platelets. The experimental determination of the aspect ratio by centrifugation experiments led to a value of about 320, demonstrating a highly anisotropic particle morphology. Applying this aspect ratio made possible a reliable determination of platelet size by AF 4 -MALS.

  millilitres of the same migrates were acid digested and a range of elements characteristic to clay were analysed by ICP-MS for additional characterization of the elemental composition of clay.

  ICP-MS signals. Isotopes of a minor and a major element, 90 Zr and 26 Mg, were selected for on-line detection of Cloisite®Na + by ICP-MS coupled with AF 4 -MALS, because of their low background level or high abundance in the clay matrix, respectively. The fact that the light scattering and ICP-MS signals coincide in time show that particles detected by the MALS detector are in fact clay platelets. Finally, the rms radii derived from the MALS detection show that the platelet radii of the large unresolved AF 4 peak eluting between 10 and 45 min (Figure 4) are in the range 20-250 nm.

Figure 5 B

 5 Figure 5 B. The results show that nano-sized particles were detected in the blank food simulant,

  Figure 1. X-ray diffraction patterns of, A, PLA containing 5 % Cloisite®30B and B, pure

Figure 2 .

 2 Figure 2. TEM micrograph of a biopolymer nanocomposite consisting of PLA with 5 %

Figure 3 .

 3 Figure 3. Mass spectra of Cloisite®Na + . The spectra have been divided into two. A, signal intensity

Figure 4 .

 4 Figure 4. AF 4 -MALS-ICP-MS detection of nanoclay in the supernatant of a suspension of 50 mg Cloisite® Na + in 5 ml water following centrifugation for 50 minutes at 4700 rpm. The four curves represent a main constituent ( 26 Mg, left y-axis), a trace constituent ( 90 Zr, right y-axis), the 90° light

Figure 5 .

 5 Figure 5. Evaluation of the cumulative number fraction curves from suspensions of nanoclays. (A) Supernatant following centrifugation for 10 minutes of Cloisite®Na + and (B) Food simulant spiked

Figure 6 .Figure 7 .

 67 Figure 6. Combined AF 4 -MALS and ICP-MS data for a migrate from PLA containing 5 %

  Figure 4.

  Figure 5.

  Figure 6.

  elements like28 Si,27 Al,58 Fe and 26 Mg in the clays, the lower signal intensities of the mass spectra demonstrate that Cloisite®Na + also contained a range of trace or ultra-trace isotopes like 7 Li, 45 Sc, 48 Ti, 51 V, 55 Mn, 63 Cu, 64 Zn, 69 Ga, 88 Sr, 89 Y, 90 Zr, 93 Nb, 120 Sn, 138 Ba, 139 La, 140 Ce, 141 Pr, 144 Nd, 156 Gd,180 Hf, 181 Ta, 208 Pb, 232 Th and 238 U.

Table 1 .

 1 Operating conditions for the AF 4 -MALS and the ICP-MS AF 4 -MALS operating conditions Channel Short channel with a "wide" 350 µm thick spacer and a 10 kDa pore size regenerated cellulose membrane

			Food Additives and Contaminants		Page 30 of 30
	Carrier		0.05% SDS and 200 ppm sodium azide in water
	Fitting algorithm		Berry 2 nd order polynomial fit
	Collection interval		1 second
	Flow Field-flow fractionation program
	Time (min)	Detector Flow	Cross Flow	Focus Flow	Injection Flow	Mode
		(ml min -1 )	(ml min -1 )	(ml min -1 )	(ml min -1 )
	0-2	1	0.5	0	0	Elution
	2-3	1	-	2	0	Focus
	3-9	1	-	2	0.2	Focus+Injection
	9-10	1	-	2	0	Focus
	10-30	1	0.5-0.1 linear gradient	0	0	Elution
	30-52	1	0.1	0	0	Elution
	52-61	1	0	0	0	Elution
	61-70	1	0	0	0.2	Elution+Injection
	ICP-MS operating conditions
			Mass scan m/z 6-254	Time resolved detection of m/z 26 and 90
	RF power		1500 W	1550 W
	Carrier gas		0.83 L min -1	0.88 L min -1
	Makeup gas		0.33 L min -1	0.31 L min -1
	Number of points per mass	3	1
	Acquisition time		111 s	4198 s
	Number of repetitions	1	1

http://mc.manuscriptcentral.com/tfac Email: fac@tandf.co.uk

Table 2 .

 2 Results from the migration studies Values are corrected for the blank simulant b Percentiles are given for the size distribution of measured radius of nanoparticulates c Based on signal intensities of 26 Mg being significantly higher blank

			Parameter measured
	Sample	Total migration	rms radius range	rms radius	Cloisite®30B
		(mg dm -2 )	corresponding to	corresponding to	detected in
			20 th -80 th	90 th percentile b	migrates by
			percentile b	(nm)	ICP-MS after
			(nm)		acid digestion c
	Blank simulant	0.4	20-31	40	No
	Pure PLA	1.7 a ± 0.6	20-40	100	No
	PLA/5%	6.7 a ± 0.5	190-230	235	No
	Cloisite®30B
	Spiked simulant	2.1 a ± 0.5	60-110	130	Yes
	a

Page 31 of 30 http://mc.manuscriptcentral.com/tfac Email: fac@tandf.co.uk Food Additives and Contaminants

http://mc.manuscriptcentral.com/tfac Email: fac@tandf.co.ukFood Additives and Contaminants

Acknowledgements


Lei Yu

Jean-Pierre Barbot

Gang Zheng

Hong Sun

J-P Barbot

Compressive Sensing for Cluster Structured Sparse Signals: Variational Bayes Approach

Keywords: Compressive Sensing, Cluster Structure, Variational Bayes. I. INTRODUCTION

Compressive Sensing (CS) provides a new paradigm of sub-Nyquist sampling which can be considered as an alternative to Nyquist sampling theorem. In particular, providing that signals are with sparse representations in some known space (or domain), information can be perfectly preserved even with small amount of measurements captured by random projections. Besides sparsity prior of signals, the inherent structure property underline some specific signals is often exploited to enhance the reconstruction accuracy and promote the ability of recovery. In this paper, we are aiming to take into account the cluster structure property of sparse signals, of which the nonzero coefficients appear in clustered blocks. By modeling simultaneously both sparsity and cluster prior within a hierarchical statistical Bayesian framework, a nonparametric algorithm can be obtained through variational Bayes approach to recover original sparse signals. The proposed algorithm could be slightly considered as a generalization of Bayesian CS, but with a consideration on cluster property. Consequently, the performance of the proposed algorithm is at least as good as Bayesian CS, which is verified by the experimental results.

yuleiwhu@gmail.com, barbot@ensea.fr, gang.zheng@inria.fr, hongsun@whu.edu.cn. Block-CoSaMP [START_REF] Baraniuk | Model-based compressive sensing[END_REF], [START_REF] Eldar | Block-sparse signals: Uncertainty relations and efficient recovery[END_REF] Yes

Dynamical Programming [START_REF] Cevher | Recovery of clustered sparse signals from compressive measurements[END_REF] Yes X X LaMP [START_REF] Cevher | Sparse signal recovery using markov random fields[END_REF] Yes X X X

Struct OMP [START_REF] Huang | Learning with structured sparsity[END_REF] Yes X X

CluSS MCMC [START_REF] Yu | Bayesian compressive sensing for cluster structured sparse signals[END_REF] Yes X X X X

1 " " denotes necessary for the corresponding algorithm.

2 "X" denotes unnecessary for the corresponding algorithm.

3 "-" means no consideration for this algorithm.

[23], Bayesian CS [START_REF] Ji | Bayesian compressive sensing[END_REF], Iterative Hard Threshold (IHT) algorithm [START_REF] Blumensath | Iterative hard thresholding and l0 regularisation[END_REF], etc. While besides sparse prior, inherent structures underlying the sparse patterns have been widely employed to improve the recovery accuracy and promote the efficiency, [START_REF] Baraniuk | Model-based compressive sensing[END_REF], [START_REF] Blumensath | Sampling theorems for signals from the union of finite-dimensional linear subspaces[END_REF], [START_REF] Cevher | Sparse signal recovery using markov random fields[END_REF], [START_REF] Cevher | Recovery of clustered sparse signals from compressive measurements[END_REF], [START_REF] Eldar | Robust recovery of signals from a structured union of subspaces[END_REF], [START_REF] He | Exploiting structure in wavelet-based bayesian compressive sensing[END_REF], [START_REF] Eldar | Block-sparse signals: Uncertainty relations and efficient recovery[END_REF].

In this paper, we focus on the cluster structured sparse signals, of which significant coefficients appear in clustered blocks. This kind of sparse pattern is often exploited in many concrete applications, such as multi-band signals, gene expression levels, source localization in sensor networks, MIMO channel equalization, magnetoencephalography [START_REF] Baraniuk | Model-based compressive sensing[END_REF], [START_REF] Blumensath | Sampling theorems for signals from the union of finite-dimensional linear subspaces[END_REF], [START_REF] Eldar | Block-sparse signals: Uncertainty relations and efficient recovery[END_REF]. Existing algorithms designed for cluster structured sparse signals always require lots of pre-defined information (Tab. I), such as (a) number of clusters; (b) size of each clusters; (c) positions where clusters are; (d) number of significant coefficients (Sparsity). However, these priors can never be known in real applications, and thus a nonparametric recovery algorithm for cluster structured sparse signal is appealed in practical problems.

A. Motivation 1) From graphical Bayesian model to CS: Considering the process of CS measurement as a hierarchical Bayesian model, namely, a graphical model [START_REF] Cevher | Sparse signal recovery and acquisition with graphical models[END_REF], it provides a new framework for CS [START_REF] Tipping | Sparse bayesian learning and the relevance vector machine[END_REF], [START_REF] Ji | Bayesian compressive sensing[END_REF], [START_REF] Babacan | Bayesian compressive sensing using laplace priors[END_REF] and leads to a nonparametric recovery algorithm. In this framework, the sparse constraint is injected through some sparse priors: a Gaussian distribution together with an Inverse Gamma on the invariance, a Laplace distribution, etc. The interpretation of CS with Bayesian model provides a systematic framework, where one could conveniently consider other priors, such as structures on sparse pattern [START_REF] He | Exploiting structure in wavelet-based bayesian compressive sensing[END_REF], dependencies between multiple signals [START_REF] Ji | Multitask compressive sensing[END_REF], and so on. Moreover, rather than providing a point estimation for sparse signals, a full posterior density function is provided, which yields "error bars" on the estimated sparse signals. These error bars can be used to give a sense of confidence of the recovered sparse signals.

2) From latent variable model to structures: In a probabilistic, Bayesian approach, through Graphical Models (GMs) [START_REF] Cevher | Sparse signal recovery and acquisition with graphical models[END_REF], [START_REF] Bishop | Pattern Recognition and Machine Learning (Information Science and Statistics)[END_REF], latent variables are often exploited to describe the dependencies (or joint probability distributions)

between observations and parameters. This method is usually called latent variable model [START_REF] Blei | Latent dirichlet allocation[END_REF], and possibly, results in some non-parametric approaches to Bayesian estimators. By imposing geometrical relations underlying the sparse pattern, structures of the sparse coefficients can be conveniently described by latent variable model [START_REF] Cevher | Sparse signal recovery using markov random fields[END_REF], [START_REF] He | Exploiting structure in wavelet-based bayesian compressive sensing[END_REF].

3) From MCMC to VB: In the last work [START_REF] Yu | Bayesian compressive sensing for cluster structured sparse signals[END_REF], a Bayesian approach to reconstruct the cluster structured sparse signals from compressed measurements has been proposed, where an MCMC re-sampling procedure is exploited for Bayesian inference. It is well known that even though MCMC is capable to find the global solution via infinite MCMC iterations, it cannot guarantee the convergence in finite iterations, which is not applicable all the time.

Consequently, we turn to modifying the Bayesian model to make it conjugate and thus can lead to a deterministic solution, through a variational Bayesian method [START_REF] Beal | Variational algorithms for approximate bayesian inference[END_REF], of which the main idea is to optimize the lower bound of the log marginal likelihood function, and simultaneously give a maximum of the posterior distribution.

B. Contribution of this paper

The main contribution of this paper is to exploit the statistical graphical model to describe the cluster structured sparse signals and hence lead to a deterministic algorithm through the variational Bayesian method. The idea of this work is largely inspired by [START_REF] He | Tree-structured compressive sensing with variational bayesian analysis[END_REF], where the resulted algorithm is dedicated to solve the tree structured sparse inverse problems. Even though, it is different from the work of [START_REF] He | Tree-structured compressive sensing with variational bayesian analysis[END_REF]: the considered structures are different, where the tree structure considered in [START_REF] He | Tree-structured compressive sensing with variational bayesian analysis[END_REF] is a directional graphical model but cluster structure in this paper is more likely an undirectional graphical model, and thus it results different models.

C. Outline

The following sections will introduce the proposed algorithm in detail. In section II, the Bayesian clustered sparsity model is addressed, where both the sparse prior and cluster prior are considered. Then using variational Bayes method, the inference of the introduced Bayesian model is implemented in section III. After that, in section IV, some simulations are presented to show the performance of the proposed algorithm. The paper ends with a conclusion.

II. BAYESIAN CLUSTERED SPARSITY MODEL FOR CS

In the framework of CS, the sampling process can be modeled as a vector y ∈ R m , captured by a multiplication of sensing matrix Φ ∈ R m×n and the original signals θ ∈ R n , then plus an error, ǫ, as follows:

y = Φθ + ǫ (1)
Suppose that the perturbation ǫ is white, i.e. ǫ ∼ N (0, σ 0 I) and thus y ∼ N (Φθ, σ 0 I), with I an all one vector with appropriate dimension. To infer the posterior of the noise variance σ 0 , a Gamma prior is assigned on the inverse of noise variance α 0 = σ -1 0 , i.e. α 0 ∼ Gamma(c, d), conjugate to Gaussian distribution.

In order to introduce both the sparse and cluster prior inside the Bayesian model, we exploit a latent variable z to indicate whether the corresponding element of θ is nonzero, i.e. θ = w • z, where • is point wise multiplication and w ∼ N (0, σ). Meanwhile, Gamma prior is assigned on the inverse of weight variance α = σ -1 , i.e. α i ∼ Gamma(a, b), with α i the i-th element of α, i ∈ {1, ..., n}. The overall prior on w with respect to a, b can be evaluated analytically through the integration over α, and it corresponds to the Student-t distribution [START_REF] Tipping | Sparse bayesian learning and the relevance vector machine[END_REF]. With appropriate choice of a, b, the Student-t distribution strongly peaked about w = 0, and thus the overall prior on w favors sparseness.

Meanwhile, we suppose that the latent variable z is drawn from a Bernoulli distribution with probability π, i.e. for each element of z, z i ∼ Bernoulli(π i ). To model the cluster prior over the coefficients θ, relations between the current component θ i and its neighbors should be considered, called the cluster pattern of θ i . Define the neighborhood of location i over the coefficients θ, U i = {j|D(i, j) = 1, j = i} with D(i, j) the Euclidean distance between i and j. Denote V N the set of all locations over the coefficients θ, V n = {1, ..., n}, then define

J i U i ∩ V n .
Hence we can denote θ Ji the set of components located in the neighbor of the ith coefficient θ i . And denote ν i = j∈Ji z j the number of nonzero neighbors for the ith element θ i . Then we can use this counter ν to categorize the relations into 3 different cluster patterns, i.e. Pattern 0: "strongly eliminate", when ν i = 0, namely, the ith element should have large probability of being eliminated, see Fig. 1a; Pattern 1: "weakly eliminate", when ν i = 1, namely, it takes small probability of being eliminated (because of sparsity), see Fig. 1b; Pattern 2: "strongly plump", when ν i = 2, namely, it should be plumped up with large probability, see Fig. 1c. Then according to the cluster patterns, the mixing weight π i is chosen by the following pattern selection procedure:

π i =          π 0 i , if Pattern 0 π 1 i , if Pattern 1 π 2 i , if Pattern 2 (2)
where π 0 i , π 1 i and π 2 i are drawn from different Beta distribution 1 .

π 0 i ∼ Beta(e 0 , f 0 ) π 1 i ∼ Beta(e 1 , f 1 )
π 2 i ∼ Beta(e 2 , f 2 ) (3)
1 Since Beta distribution is a conjugate prior to Bernoulli likelihood with p the model parameters.

In order to clarify the dependance between the random variables, the distributions for π could be rewritten as follows:

π i | e,f ,zJ i ∼ p(π i |e, f , z Ji ) (4)
where e {e 0 , e 1 , e 2 }, f {f 0 , f 1 , f 2 }. Above all, the overall Bayesian clustered sparsity model is as follows:

θ = w • z w ∼ N (0, α -1 ) α i ∼ Gamma(a, b) z i ∼ Bernoulli(π i ) y ∼ N (Φθ, α -1 0 I) α 0 ∼ Gamma(c, d) (5)
where π i is described by ( 2) and ( 3). Moreover, it could be depicted by a graphical model in Fig. 2.

III. VARIATIONAL BAYESIAN INFERENCE

In this section, we derive a variational Bayesian (VB) algorithm [START_REF] Beal | Variational algorithms for approximate bayesian inference[END_REF] to implement the inference for the Bayesian clustered sparsity model. The main idea of VB is to optimize the lower bound of the log marginal likelihood function, and simultaneously give a maximum for the posterior distribution. Based on the Bayesian clustered sparsity model ( 5), the complete marginal likelihood can be written as follows:

p(y|M) = dθdΘp(θ, y, Θ|M) (6)
with M = {a, b, c, d, e, f } the hyperparameters and Θ = {α, α 0 , π} the unknown random variables.

In the following, we define • x the expectation with respect to random variables2 x, x the updated estimation for random variable x, and y -k = yi =k z i w i φ i the contribution of the k-th element of sparse signal on the measurement y.

A. The VB-E Step 1) Update for w:

q(w) ∼ exp n i=1 ln p(w i |α i ) αi exp ( ln p(y|w, z, α 0 ) z,α0 )
The posterior of w can be shown to be multivariate Gaussian distribution, with mean µ and variation Σ, i.e., w ∼ N (µ, Σ), with

Σ = ( Ã + α0 ZΦ T ΦZ ) -1 µ = α0 Σ ZΦ T y (7)
where A = diag(α) and Z = diag(z). And given the updated value z, we can derive

ZΦ T ΦZ = (Φ T Φ) • zz T + diag(z • (1 -z)
Consequently, we can derive the update for w as follows:

w = µ (8)
2) Update for z: For each element of z, the posterior could be given as

q(z i ) ∼ exp( ln p(z i |π i ) πi ) exp( ln p(y -i |z i , w i , α 0 ) wi,α0 )
Thus the probability that z i = 1 is proportional to

p(z i = 1) ∼ exp ( ln π i ) exp - α0 2 ( w 2 i φ T i φ i -2 wi φ T i y -i ) (9)
where w 2 i = w2 i +σ ii with σ ii the i-th element of diagonal entries of Σ. The probability that z i = 0 is proportional to

p(z i = 0) ∼ exp( ln(1 -π i ) ) (10)
where the update for π i could be referred in [START_REF] Eldar | Robust recovery of signals from a structured union of subspaces[END_REF] in following subsections.

Thus, the update for z can be easily obtained

zi = p(z i = 1) p(z i = 0) + p(z i = 1) (11)
B. The VB-M Step 1) Update for α: For each element of α i , the posterior could be given as

q(α i ) ∼ p(α i |a, b) exp( ln p(w i |α i ) wi ) ∼ Gamma (α i |a ′ , b ′ ) (12)
with

a ′ = a + 1 2 , b ′ = b + w 2 i 2
where the conjugate property between Gamma prior and Gaussian distribution is used. Thus, the update for α i is obtained:

αi = a ′ /b ′ (13)
2) Update for α 0 :

q(α 0 ) ∼ p(α 0 |c, d) exp ( ln p(y|w, z, α 0 ) w,z ) ∼ Gamma(α 0 |c ′ , d ′ ) (14)
with

c ′ = c + m 2 , d ′ = d + y -Φ(w • z) 2 w,z 2
where the conjugate property between Gamma prior and Gaussian distribution is used. And given the updated value of w, z, the expectation could be derived:

y -Φ(w • z) 2 w,z = y T y -2( w • z) T Φ T y + I T [ zz T • ww T • (Φ T Φ)]I (15)
where zz T = zz T + diag(z • (1z)) and ww T = w wT + Σ.

Thus the update for α 0 could be obtained:

α0 = c ′ /d ′ (16)
3) Update for π: For each of the π i , given the updated value for z, then we can easily calculate portion of sparsity pattern effecting on π i :

           p(Pattern 0) = (1 -zi-1 )(1 -zi+1 ) p(Pattern 1) = zi-1 (1 -zi+1 ) + (1 -zi-1 )z i+1 p(Pattern 2) = zi-1 zi+1
Thus the posterior of π j i , with j denoting the sparsity pattern, could be written as:

q(π j i ) ∼ p(π j i |e j , f j ) exp( ln p(z i |π j i ) zi ) ∼ Beta(π j i |e ′j i , f ′j i )
where e ′j i = e j + p(Pattern j)z i

f ′j i = f j + p(Pattern j)(1 -zi ) thus ln π j i = ψ(e ′j i ) -ψ(e ′j i + f ′j i ) ln(1 -π j i ) = ψ(f ′j i ) -ψ(e ′j i + f ′j i )
where ψ(x) = d dx ln Γ(x) is a digamma function, and then the update for π i could be obtained:

ln π i = 2 j=0 p(Pattern j) ln π j i ln(1 -π i ) = 2 j=0 p(Pattern j) ln(1 -π j i ) (17)

C. Summary of the algorithm and acceleration

Given observation y and sensing matrix Φ, the algorithm could be summarized as Algorithm 1. Update unknown parameters via equations ( 13), ( 16) and ( 17):

Θ = { α, α0 , π} 3:
Update latent parameters via equations ( 8) and ( 11):

x = { w, z} 4: until C
At the iteration t, define the residual as

Res(t) = y -Φ θ(t) thus C Res(t) √ mσ 0
where σ 0 is the invariance of noise.

1) Comparison to BCS:

In the framework of BCS, the inverse problem (reconstruction) is solved by a Sparse Bayesian Learning (SBL) [START_REF] Wipf | Sparse bayesian learning for basis selection[END_REF], which has been proven to be capable to find the unique sparse solution (Theorem 2 and 3 in [START_REF] Wipf | Sparse bayesian learning for basis selection[END_REF]).

On the other hand, by setting z ≡ 1, the proposed model would degenerate into BCS model, i.e. without consideration on clusters. Roughly, BCS model could be considered as a special case of Bayesian clustered sparsity model, where the prior on clusters provides a guidance of how to assign the priority of choosing bases: "considering first your neighbors". While, in BCS model, this priority does not exist. Consequently, the performance of the proposed algorithm is at least as good as that of BCS and the worst case is when all significant entries of sparse signal are distributed isolated.

2) Acceleration: In the current form, the algorithm requires an inverse problem to update the value of Σ, which is an n × n matrix, and thus requires an O(n 3 ) operations. This can be problematic since in cases of large length of signals, n might be quit large. To alleviate this problem, we compute the Σ as follows:

Σ = ( Ã + α0 ZΦ T ΦZ ) -1 = ( Ã + α0 ZΦ T Φ Z + diag(z • (1 -z)) • (Φ T Φ)) -1 = (B + α0 ΦT Φ) -1 (18)
where Φ = Φ Z and

B = Ã + diag(z • (1 -z)) • (Φ T Φ)
which is a diagonal matrix, thus its inverse can be easily computed by directly inverse the elements located at the diagonal, D = B -1 . Then using the inverse identity property, it has

Σ = D -D ΦT (α -1 0 I + ΦD ΦT ) -1 ΦD (19)
where matrix α -1 0 I + ΦD ΦT is with dimension of m × m, which reduces the operations to O(m 3 ), with m ≪ n.

IV. EXPERIMENTS

To distinguish the algorithm proposed in this paper from the former CluSS algorithm in [START_REF] Yu | Bayesian compressive sensing for cluster structured sparse signals[END_REF], for simplicity, we denote the proposed algorithm as CluSS-VB. All the hyperparameters are fixed for every experiments as follows:

a = b = c = d = 1e -6
, (e 0 , f 0 ) = (1/3, 2/3), (e 1 , f 1 ) = (1/3, 1/3) and (e 2 , f 2 ) = (2/3, 1/3). The following experiments are organized as follows. A first glance on the performance of CluSS-VB on synthetic cluster structured sparse signals is given. Afterwards, with respect to the oversampling rate, defined as m/s, we compare the recovery accuracy between CluSS-VB and other state-of-the-art CS algorithms, respectively, BP [START_REF] Chen | Atomic decomposition by basis pursuit[END_REF], CoSaMP [START_REF] Needell | Cosamp: Iterative signal recovery from incomplete and inaccurate samples[END_REF], and Baysian Compressive Sensing (BCS) [START_REF] Ji | Bayesian compressive sensing[END_REF]. Then in order to verify the performance for mismatched models, we compare the reconstructions by each of the algorithms with augmenting number of clusters. Meanwhile, the robust to measurement noises is considered and an application on real musical signal is given.

In the following experiments, all sensing matrix are constructed through Gaussian ensemble with normalized row vectors. Moreover, if without clarifying, the measurements are corrupted by a white noise with variance σ 0 = 0.01. Note that the accuracy of reconstruction is evaluated through SNR between original sparse signal θ and its reconstruction θ:

SNR = 20 log 10 θ θ -θ (20)

A. General view

First, we shall give a glance on the performance of the proposed algorithm and show the basis selection procedure during iterations. The original sparse signal is generated with length n = 256, sparsity s = 30 and clusters k = 2,

where cluster information on sizes, locations and numbers is chosen totally randomly and thus is unknown 3 .

Moreover, the spikes are randomly generated by Gaussian distribution which coincides with the sparsity model and only m = 2s = 60 measurements are obtained by random projections.

In Fig. 3, the inference of z and w are shown along with the iterations. Apparently, the result shows that basis choosing priority is considering the neighbors first. The final result is shown in Fig. 4, where the convergence evolution for Res(t) is also shown. Moreover, the inference for measurement noise invariance σ0 = 0.0077, which coincides with the original set σ 0 = 0.01.

Meanwhile, we also exploit the other state-of-the-art algorithms (BP, CoSaMP, BCS) respectively to reconstruct the same sparse signal, shown in Fig. 5. From the comparison, except for CluSS-VB, all other algorithms cannot be able to reconstruct this sparse signal with only m = 60 measurements. On the other hand, like BCS, CluSS-VB can also provide "error bars" for the final reconstruction, which can be used to evaluate the confidence for the estimation. 3 For Block-CoSaMP proposed in [START_REF] Baraniuk | Model-based compressive sensing[END_REF], it is impossible to recover this kind of cluster structured signals successfully, and thus we did not make comparison with Block-CoSaMP.

B. Performance with respect to oversampling rate m/s

In order to deeply verify the performance of CluSS-VB, we shall take into account the oversampling rate m/s, which determines how many measurements are enough to reach the successful reconstruction. The cluster structured sparse signals are randomly generated with length n = 256 and sparsity s = 30, and their nonzero entries drawn from a Gaussian distribution with zero mean and unit variance. Respectively, sparse signals with 1 cluster and 2 clusters are both considered. In each case, for fixed cluster structured sparse signal, the number of measurements is ranging from 50 to 150 with step 10, and for each step, we run the program 100 times with different sensing matrix.

Fig. 6 shows the corresponding results, where (a) is for 1 cluster case and (b) is for 2 clusters case. The results

show that CluSS-VB improves both the accuracy and the ability of reconstruction: even with very low oversampling rate, it can obtain desirable reconstructions. And more interestingly, when the number of clusters k goes to the sparsity s, the performance of CluSS-VB tends to converge to the performance of BCS, which coincides with the discussion in section III.

E. Experiments on real musical signals

In the last experiment, we apply the proposed algorithm on real musical signals, which have the property of cluster structured sparsity if considering signals in the frequency domain 4 , as shown in Fig. 9, where the significant spectrums are almost clustered together. We choose a clip of music of Mozart played by flute as the test example.

The CS procedure is carried out piecewise with length n = 256 for each of the pieces. Then varying the number of measurements obtained by random projections, from m = ⌊n/2⌋ to ⌊n/5⌋, corrupted by white noises with variance σ 0 = 0.01, respectively, we use the different algorithms to recover the original musical signals piecewise from the compressed measurements and then concatenate the recovered pieces. As shown in Fig. 10, with only m = ⌊n/5⌋ measurements, the spectrograms of recovered signals are given for each algorithm. It is shown that CluSS-VB can desirably preserve the clusters of the original sparsity while suppress the isolated spikes, and hence can give better

Fig. 1 .

 1 Fig. 1. Three different cluster pattern for 1D signals. (a) Pattern 0; (b) Pattern 1; (c) Pattern 2.

Fig. 2 .

 2 Fig. 2. Graphic depiction for Bayesian clustered sparsity model.

Algorithm 1

 1 CS for Cluster Structured Sparse Signals via VB Initialization The hyperparameters M = {a, b, c, d, e, f }, the pre-estimated w(0) = Φ T y, z(0) = w max w , a stop criterion C.

Fig. 3 .

 3 Fig. 3. Inference of z (left column) and w(right column) along with iterations: (from top to bottom) 5, 10, 20, 40, 60 iterations.

Fig. 4 .

 4 Fig. 4. Performance of CluSS-VB. (a) Inference of latent variable z; (b) Recovery of sparse signal; (c) Convergence of Res(t).

Fig. 5 .

 5 Fig. 5. General view of performance of CluSS-VB. (a) Original Signal, with length n = 256, sparsity s = 30 and clusters k = 2. (b) Recovery via BP and its SN R = 7.0449dB; (c) Recovery via CoSaMP and its SN R = 2.8955dB; (d) Recovery via BCS and its SN R = 7.2022dB; (e) Recovery via CluSS-VB and its SN R = 30.5702dB.

  Sparse signal with length n = 256, sparsity s = 30, and k = 1 clusters.

  Sparse signal with length n = 256, sparsity s = 30, and k = 2 clusters.

Fig. 6 .

 6 Fig. 6. Performance comparison with respect to oversampling rate.

  : 20 log( y / ǫ ) (dB) SNR for recovery: 20 log( θ / θ -θ ) (dB)

Fig. 7 .

 7 Fig. 7. Robustness to different level of measurement noise.

  with 80 measurements.

  with 100 measurements.

Fig. 8 .

 8 Fig. 8. Effects of clusters on the performance.

Fig. 9 .

 9 Fig. 9. Spectrogram of a musical signal.

Fig. 10 .

 10 Fig. 10. Spectrogram of reconstructions of musical signals via CluSS-VB (left-top), BCS (right-top), BP (left-bottom) and CoSaMP (rightbottom).

TABLE I COMPARISON

 I BETWEEN DIFFERENT CS RECOVERY ALGORITHMS.

	Algorithms	Is for Cluster?	Num. of Clusters	Size of Clusters	Fixed Cluster Positions	Sparsity
	Greedy algorithms (CoSaMP, OMP, etc)	No	-	-	-
	Linear Programming (BP, etc)	No	-	-	-
	Iterative Thresholding (IHT, IST)	No	-	-	-	X
	Bayesian CS	No	-	-	-	X

TABLE II RECOVERY

 II SNR FOR DIFFERENT NUMBER OF MEASUREMENTS. (DB)

	Measurements	⌊n/2⌋	⌊n/3⌋	⌊n/4⌋	⌊n/5⌋
	BP	12.6384	9.8339	7.3692	5.8631
	CoSaMP	10.5388	7.9345	4.8592	4.4899
	BCS	13.1023	9.688	7.0389	5.4191
	CluSS-VB	12.7616	11.1035 9.2251 7.5384

Denote • if there is only one random variable inside the triangle bracket.

Other transforms are not considered since it is not the important issue in this experiment.

This work was supported by "Bourses Doctorales en Alternance", PEPS-A2SDC of INSIS-CNRS, NSFC (No. 60872131) and FEDER through CPER 2007-2013 with INRIA Lille-Nord Europe.