Phuoc Dien Pham
email: phuoc-dien.pham@univ-montp2.fr

Sophie Monge
email: sophie.monge-darcos@univ-montp2.fr

Vincent Lapinte
email: vincent.lapinte@univ-montp2.fr

Yann Raoul
email: y.raoul@prolea.com

Jean Jacques Robin

Various radical polymerizations of glycerol-based monomers -Review

Keywords: Radical polymerization, glycerol derivatives, acrylate monomers, photopolymerization, acceptor-donor

This review aims at making the state of the art of the radical polymerizations of glycerolbased monomers, especially (meth)acrylate and vinylic/allylic derivatives. Acrylate and methacrylate derivatives bearing numerous secondary groups like ester, cyclocarbonate, acetal and urethane were involved in free radical polymerization and radical controlled polymerizations as well as in the photopolymerization. Glycerol-based monomers bearing a vinylic or an allylic polymerizable group polymerize by acceptor-donor process and produce alternating copolymers. From all these monomers, the synthesis of functionalized polymers was possible, leading to various valuable applications.

Introduction

Glycerol is one of the most important feedstock in the modern oleochemical industry [START_REF] Zheng | Commodity Chemicals Derived from Glycerol, an Important Biorefinery Feedstock[END_REF] obtained as a by-product in saponification of fats, in fatty acid/ester production [START_REF] Morrison | Glycerol[END_REF], in microbial fermentation [START_REF] Wang | Glycerol production by microbial fermentation: A review[END_REF] and in biodiesel manufacture by trans-esterification of vegetable oils. Due to its competitive cost, worldwide availability, and built-in functionality, the glycerin and its derivatives become useful for numerous commercial applications. As production of glycerin is nowadays exceeding the demand, the valorization of the latter has emerged as a challenging trend [START_REF] Behr | Improved utilisation of renewable resources: New important derivatives of glycerol[END_REF]. Thus the glycerol is currently used in industry as an intermediate in the synthesis of numerous compounds (glycerin carbonate (GC), solketal (Sol), acrylic acid, etc.) [START_REF] Zhou | Chemoselective catalytic conversion of glycerol as a biorenewable source to valuable commodity chemicals[END_REF]. For instance, GC is employed as solvent [6], surfactant [START_REF] Mouloungui | Study of the acyl transfer reaction: Structure and properties of glycerol carbonate esters[END_REF] or in the synthesis of polyurethane [START_REF] Vora | A new approach to 3-miktoarm star polymers using a combination of reversible addition-fragmentation chain transfer (RAFT) and ring opening polymerization (ROP) via "Click" chemistry[END_REF], polyether [START_REF] Rokicki | Hyperbranched aliphatic polyethers obtained from environmentally benign monomer: glycerol carbonate[END_REF] and polycarbonate [10]. More specifically a lot of glycerol derivatives bear a polymerizable unit such as acrylate, methacrylate, allylic or vinylic group as summarized in Figure 1. Among (meth)acrylates, some of them have not a secondary function: glycerol acrylate (GA), glycerol methacrylate (GMA), dihydroxypropylacrylate (DHPA), dihydroxyisopropyl methacrylate (DHIMA) whereas others have a secondary function: i) cyclocarbonate such as glycerol carbonate acrylate (GCA), glycerol carbonate methacrylate (GCMA), carbonate glycerol carbonate acrylate (CGCA) and cyclocarbonate carbamate methacrylate (CCCM), ii) acetal such as solketal methacrylate (SolMA), solketal acrylate (SolA), 5-membered cyclic acetal acrylate (5-CAA) and 6-membered cyclic acetal acrylate (6-CAA), and iii) various pendant chains bearing ether (GM-5), ester (GM-1, GM-2, GM-3, GM-4) or phosphine oxide group (GM-6). Di(meth)acrylates are also described including bis-glycerol methacrylate (Bis-GMA), 1,3glycerol dimethacrylate (GDM), ethyl -chloromethylacrylate-glycerol diacrylate (ECMA-GD), methyl -chloromethylacrylate-glycerol diacrylate (MCMA-GD), urethane dimethacrylate (UR-1 and UR-2) as well as tri(meth)acrylates: glycerol trimethacrylate (GTM), ethyl -chloromethylacrylate-glycerol triacrylate (ECMA-GT), methyl chloromethylacrylate-glycerol triacrylate (MCMA-GT), and urethane trimethacrylate (UR-3).

Several vinylic derivatives based on glycerol are also reported (vinylic solketal (vinylSol), vinylic glycerol carbonate (vinylGC), 1-propenyl glycerol carbonate (1-propenylGC)) as well as allylic one: allylic solketal (allylSol). Among these monomers, those bearing a cyclocarbonate secondary group are interesting and employed in pharmaceutical or cosmetic fields [START_REF] Klug | in: Cosmetic or pharmaceutical compositions, emulsions comprinsing modified polysiloxanes with at least one carbonate group[END_REF] thanks to their low vapor pressure, flammability, toxicity, good biodegradability and moisturizing ability [START_REF] Aresta | Valorization of bio-glycerol: New catalytic materials for the synthesis of glycerol carbonate via glycerolysis of urea[END_REF]. The cyclocarbonate function also brings high polarity and coordinating properties [START_REF] Britz | Blends of poly(meth)acrylates with 2-Oxo-(1,3)dioxolane side chains and lithium salts as lithium ion conductors[END_REF][START_REF] Rossi | Silicon-containing liquid polymer electrolytes for application in lithium ion batteries[END_REF][START_REF] Zhang | Synthesis and ionic conductivity of mixed substituted polysiloxanes with oligoethyleneoxy and cyclic carbonate substituents[END_REF][START_REF] Zhu | Synthesis of polysiloxanes bearing cyclic carbonate side-chains -dielectric properties and ionic conductivities of lithium triflate coplexes[END_REF].

Another peculiar type of monomers concerns that bearing acetal group which is hydrolyzable in diol to notably afford polymers bearing hydroxyl functionality leading to hydrophilicity and water-soluble properties. Such derivatives can be used in the pharmaceutical field such as contact lenses, dental materials, optical lenses, drug delivery, encapsulated cells [START_REF] Kuenzler | Contact lens materials[END_REF][START_REF] Schmaljohann | Thermo-and pH-responsive polymers in drug delivery[END_REF][START_REF] Kurisawa | Gene expression control by temperature with thermo-responsive polymeric gene carriers[END_REF] for instance, as well as in material field including hydrogels [START_REF] Arima | Complement Activation by Polymers Carrying Hydroxyl Groups[END_REF][START_REF] Zhang | Preparation of temperature-sensitive poly(N-isopropylacrylamide)/β-cyclodextrin-grafted polyethylenimine hydrogels for drug delivery[END_REF], chelating agents [START_REF] Ma | Synthesis of Magnetic Chelator for High-Capacity Immobilized Metal Affinity Adsorption of Protein by Cerium Initiated Graft Polymerization[END_REF][START_REF] Johnson | Polymeric drug carriers functionalized with pairwise arranged hydroxyl and/or carboxyl groups for platinum chelation[END_REF] and functionalized polymers [START_REF] Giardi | Synthesis of polyoxazolines using glycerol carbonate derivative and end chains functionalization via carbonate and isocyanate routes[END_REF]. Among all possibilities, some amphiphilic block copolymers bearing hydrophilic block with hydroxyl groups allow applications like thermoresponsive polymers [START_REF] Eggenhuisen | Libraries of Statistical Hydroxypropyl Acrylate Containing Copolymers with LCST Properties Prepared by NMP[END_REF][START_REF] Weaver | Stimulus-Responsive Water-Soluble Polymers Based on 2-Hydroxyethyl Methacrylate[END_REF]. The free radical polymerization (FRP) and the photopolymerization [START_REF] Andrzejewska | Photopolymerization kinetics of multifunctional monomers[END_REF][START_REF] Oster | Photopolymerization of vinyl monomers[END_REF] are the simplest ways to synthesize polymers with or without radical initiator under thermal or photo conditions (Scheme 1). The second process is appropriate to specific reactive groups such as allylic and vinylic unlike the first one which is employed with acrylic and methacrylic monomers. In some cases, the control of the molecular weight of the polymer with narrow mass distribution, the control of the composition and the architectures (linear diblock, triblock, Y-shaped triblock, 4 arm star copolymers, grafted copolymers) were achieved. Thus, the control/living radical polymerization (CRP) [START_REF] Braunecker | Controlled/living radical polymerization: Features, developments, and perspectives[END_REF][START_REF] Matyjaszewski | Atom Transfer Radical Polymerization[END_REF] using ATRP, and NMP was carried out. Finally, the acceptor-donor polymerization is described for the vinylic, allylic monomers based on glycerol such as allylSol, vinylSol, vinylGC, and 1-propenylGC (Figure 1). This review aims at providing an overview about the different radical polymerization processes of various monomers based on glycerol such as free radical polymerization, controlled radical polymerization, photo-polymerization and the acceptor-donor polymerization. Applications of the resulting materials will also be discussed.

Scheme 1.

Radical polymerizations of glycerol-based unsaturated monomers.

Free radical polymerization

Free radical polymerization (FRP) is a very widespread method of polymerization which has been used for a long time in industry for many reasons: FRP is easy to carry out, can be performed in bulk, in solution in various solvents, and also in dispersed media (suspension, emulsion…). It can be achieved for a large range of reaction temperatures (-100 to +200 °C).

Additionally, this polymerization technique can be carried out with a lot of monomers, even functionalized. Monomers can bear "standard" functional groups or heteroatoms such as phosphorus [START_REF] Zoulalian | Self-Assembly of Poly(ethylene glycol)-Poly(alkyl phosphonate) Terpolymers on Titanium Oxide Surfaces: Synthesis, Interface Characterization, Investigation of Nonfouling Properties, and Long-Term Stability[END_REF] or silicon [START_REF] Gayet | Polymer nanocomposite ionogels, high-performance electrolyte membranes[END_REF] for instance. Glycerol-based monomers, often bearing specific functional group(s) (Figure 2) were thus considered using free radical polymerization.

In general, FRP is described by three main steps. In the initiation step, radicals are produced and then add a first monomer molecule. Initiators are typically diazoic compounds such as azobisisobutyronitrile (AIBN), or peroxides, with benzoyl peroxide for instance. Redox systems and energy sources can also produce initiating species. The second step (propagation) corresponds to the growth of the macromolecular chains. Finally, two active chains react together in the termination step to produce macromolecular dead chains. Glycerol-based mono(meth)acrylate monomers were polymerized by free radical polymerization. AIBN is used as initiator in all cases, whereas solvent and temperature vary.

Polymerization of acrylates and methacrylates with 2-oxo-1,3-dioxolane (glycerol carbonate) was carried out [START_REF] Britz | Blends of Poly(meth)acrylates with 2-Oxo-(1,3)dioxolane Side Chains and Lithium Salts as Lithium Ion Conductors[END_REF]. Resulting poly(GCA) and poly(GCMA) exhibit a good thermal stability and the glass transition temperature range from 11 to 93 °C, depending on the structure of the glycerol-based monomer. Concerning the application, these polymers showed significant ion conductivity when blended with lithium salts.

In a fundamental study, the solketal methacrylate was polymerized by free radical polymerization [START_REF] Miranda | Binary copolymer reactivity of tert-butyl methacrylate, 2-(N,N-dimethylamino)ethyl methacrylate, solketal methacrylate, and 2-bromoethyl methacrylate[END_REF], and its reactivity ratio was determined using tert-butyl methacrylate (tBMA), 2-bromoethyl methacrylate (BEMA), or 2-(N,N-dimethylamino) ethyl methacrylate (DMAEMA) as co-monomer (Figure 3). This research aims at producing statistical polyampholytes from terpolymers. PGMA was prepared via FRP from corresponding monomer synthesized by enzymatic way in the presence of lipase catalyst known as Candida antarctica lipase B by trans-acylation of methyl acrylate [START_REF] Popescu | Free radical and nitroxide mediated polymerization of hydroxy-functional acrylates prepared via lipase-catalyzed transacylation reactions[END_REF]. The authors compared the polymers obtained from free radical polymerization and nitroxide mediated polymerization techniques in term of conversions, molecular weights and polydispersity indices.

The glycerol diacrylate (GDA) has also been studied in polymer field and was used as crosslinking agent for solid-phase synthesis. This monomer is introduced into polystyrene [START_REF] Roice | Poly(styrene-co-glycerol dimethacrylate): Synthesis, characterization, and application as a resin for gel-phase peptide synthesis[END_REF] or into polyvinyl pyrrol network [START_REF] Vijitha | Synthesis and characterization of glycerol dimethacrylate-4vinyl pyrrole[END_REF] using free radical aqueous suspension polymerization. The effect of amount of cross-linking agent on the swelling, loading and the mechanical stability of the resin is reported. Interestingly, the presence of hydroxyl groups coming from GDA moiety leads to an optimum hydrophobic-hydrophilic equilibrium of the resin. This permits a high swelling in different polar and non-polar solvents. The structure of the polymer provided easy diffusion of reagents and solvents through the resin matrix.

Poly(styrene-co-GDMA) proved to be equally efficient as commercial resins in peptide synthesis.

To conclude, glycerol-based monomers were involved in free radical polymerization and the resulting polymers proved to find interesting applications. Nevertheless, the main drawback of conventional radical polymerization is the lack of control over the molecular weight, and the impossibility to achieve complex polymeric structures. This led to the development of controlled radical polymerization techniques.

Controlled/living radical polymerization

The synthesis of well-defined macromolecules with controlled compositions, architectures, and functionalities has emerged as an important aspect of polymer science. The development of controlled/living radical polymerization (CRP) permitted to achieve well-defined polymers using a radical process easy to carry out. In the last twenty years, techniques such as atom transfer radical polymerization (ATRP) [START_REF] Matyjaszewski | Atom transfer radical polymerization[END_REF], nitroxide mediated polymerization (NMP) [START_REF] Hawker | New polymer synthesis by nitroxide mediated living radical polymerizations[END_REF], and reversible addition-fragmentation chain-transfer polymerization (RAFT) [START_REF] Barner-Kowollik | Hanbook of RAFT Polymerization[END_REF] led to the synthesis of complex architectures, with control over the molecular weight, low polydispersity indexes and complex architectures notably due to the possible functionalization of terminal end groups [START_REF] Monge | A new route for the modification of halogen end groups to amino end-functionalized poly(tert-butyl acrylate)s[END_REF]. All the CRP is based on a fast equilibrium between active and dormant species. The concentration in active species has to remain low during all the polymerization, minimizing termination and transfer reactions. Kinetics and thus control of the polymerization is easily determined by taking samples throughout the polymerization or using appropriate probes [START_REF] Darcos | In situ Fourier transform near infrared spectroscopy monitoring of copper mediated living radical polymerization[END_REF]. CRP of many monomers has been carried out and glycerol-based monomers were also considered. Almost all examples described in the literature deal with atom transfer radical polymerization (ATRP) technique. These examples will be discussed as a function of the obtained architecture. Linear diblock and triblock copolymers were prepared using solketal or glycerol carbonate (meth)acrylates. Diblock copolymers combining a biodegradable core and a hydrophilic shell were synthesized [START_REF] Kipping | Synthesis and characterization of particles consisting of a biodegradable poly(L-lactide) core and a functional hydrophilic shell[END_REF]. In such case, two polymerization techniques are involved, with the ring opening polymerization (ROP) of the lactic acid (LLA) first followed by the ATRP of the solketal acrylate (SolA). PLLA-b-PSolA diblock copolymers are thus obtained (Figure 4). These copolymers showed self-assembly property as they form aggregates in aqueous solution. Besides, degradation behavior of latter was also investigated. It was shown that enzymatic degradation proceeds much faster than hydrolytic one, within 4000 minutes. Double hydrophilic block copolymers were also prepared [START_REF] Tang | Double hydrophilic block copolymers PEO-b-PGA: synthesis, application as potential drug carrier and drug release via pH-sensitive linkage[END_REF]. The SolA was polymerized by ATRP in the presence of a poly(ethylene oxide) (PEO) macroinitiator. Then the hydrolysis of the acetal protecting group of the SolA moiety was carried out in the presence of a chlorhydric acid solution to afford a hydrophilic block (PEO-b-PGA). Finally, a reaction between the glycerol unit of the PEO-b-PGA with a fluorescent probe (pyrene derivative) led to an amphiphilic structure (Scheme 2). The micelles, with a PGA-Pyrene core and a PEO shell, self organize in THF/water solution. In acidic medium, acetal groups were deprotected to deliver the pyrene in the core of the micelles. The released ratio of pyrene varied as a function of the pH values. It is equal to 100% at pH 1.0, 80% at pH 5.0 and 9% at pH 7.4.

Moreover, the micelles are destroyed after 80% of acetal linkages cleaved. As a result, the hydrophilic PEO-b-PGA block copolymer offers a potential application as pH-sensitive drug delivery material for carbonyl bearing hydrophobic drugs. Scheme 2. Amphiphilic diblock copolymers bearing fluorescent probe from solketal moiety.

Fluorescent probe can also be introduced via the initiator, as in the case of perylene bis-imide (PBI) [START_REF] Xu | Functionalization of reduced graphene oxide nanosheets via stacking interactions with the fluorescent and water-soluble perylene bisimide-containing polymers[END_REF] which was converted into ATRP initiator for the controlled polymerization of SolA. After hydrolysis of the acetal protecting group of the SolA moiety, resulting polymer (PBIPGA) (Figure 5) becomes water soluble. Interestingly, the PBIPGA self-assemble onto reduced grapheme oxide (RGO) nano-sheets via  stacking interactions, resulting in non-covalent functionalization of RGO. Furthermore, the RGO-PBIPGA composites display good dispersity in an aqueous medium and very low cytotoxicity toward mouse 3T3 fibroblasts. As a consequence, this composite is potentially interesting for biomedical applications such as cell imaging and intracellular drug delivery.

Another kind of diblock copolymers results from the copolymerization of SolA and tBA by ATRP using bromine-terminated PSolA to initiate polymerization of tBA and vise versa [START_REF] Zhen | Atom transfer radical polymerization of solketal acrylate using cyclohexanone as solvent[END_REF]. The SolA units were hydrolyzed using chlorhydric acid without cleaving tBA moiety, whereas SolA and tBA were completely hydrolyzed to give acid acrylic and glycerol unit, respectively, in the presence of trifluoroacetic acid (Scheme 3). Scheme 3. Hydrolyzis of SolA moiety with different experimental conditions, consequently leading to different diblock copolymers.

Glyceryl methacrylate was also directly involved in the synthesis of diblock copolymers, associated with 2-hydroxypropyl methacrylate (HPMA) [START_REF] Save | Atom Transfer Radical Polymerization of Hydroxy-Functional Methacrylates at Ambient Temperature: Comparison of Glycerol Monomethacrylate with 2-Hydroxypropyl Methacrylate[END_REF]. Polymerization was achieved at room temperature. In addition, diblock copolymers were successfully synthesized in the presence of different poly(alkylene oxide)s as macroinitiators. For instance, poly(propylene oxide)-b-poly(glycerol methacrylate) (PPO-b-GMA) diblock copolymers were thermoresponsive, leading to reversible aggregation. The obtained average size of these aggregates is about 150-240 nm, depending on the block composition.

Linear triblock copolymers were also prepared by controlled radical polymerization, involving glycerol-based monomers. GMA was polymerized from a poly(propylene oxide) difunctional macroinitiator to afford water-soluble PGMA-b-PPO-b-PGMA triblock copolymers (Figure 6) [START_REF] Amado | Amphiphilic Water Soluble Triblock Copolymers Based on Poly(2,3-dihydroxypropyl methacrylate) and Poly(propylene oxide): Synthesis by Atom Transfer Radical Polymerization and Micellization in Aqueous Solutions[END_REF]. SolA was also used to produce triblock copolymers. Synthesis is accomplished with copper mediated living radical polymerization using a difunctional poly(vinyl acetate) as macroinitiator preliminary synthesized (Figure 7) [START_REF] Batt-Coutrot | Synthesis and properties of amphiphilic vinyl acetate triblock copolymers prepared by copper mediated living radical polymerisation[END_REF]. Solubility, film forming property, wettability, and thermal stability were investigated.

(AC)B(AC) type terpolymers reveal an enhanced solubility, film-forming property, and thermal stability thanks to the presence of additional methyl methacrylate repeat units. The lithium ion conductivity of all polymers bearing cyclocarbonate was determined when mixed with lithium bis(trifluoromethane) sulfonamide. In addition, ring-opening reaction of cyclocarbonate by mono-and diamine was achieved. Treatment with monoamine afforded polymer bearing hydroxyurethane whereas the reaction with diamine gave an insoluble and crosslinked polymers. Y-shaped triblock copolymers (PSolA-b-PEG-b-PSolA) (Figure 8)

were designed from a dibromine-PEG macroinitiator on which ATRP of SolA was carried out [START_REF] Yang | Synthesis of Y-shaped poly(solketal acrylate)-containing block copolymers and study on the thermoresponsive behavior for micellar aggregates[END_REF]. The Y-shaped triblock polymers exhibit micelles composed of PEG and PSolA as hydrophilic and hydrophobic block, respectively. The amphiphilic polymers self assemble, forming a core-corona type spherical micellar aggregate. The particle sizes depend on hydrophobic block length of the amphiphilic polymers and on the temperature. Furthermore, the encapsulation and release of nile red were also investigated and proved to be temperature dependent. As the biocompatibility of these copolymers to cells was revealed by MTT assay, these thermoresponsive polymer micelles can be potentially useful to produce nano-vehicles for biomaterial applications in the near future.

Star copolymers also represent architecture achievable by controlled radical polymerization.

Well-defined linear and star homo-and diblock-copolymers of PMA, PSolA, and P(MA-b-GA) n (with n, number of arms, ranging from 1 to 4) (Figure 9) were prepared by the Single-Electron Transfer Living Radical Polymerization (SET-LRP) using Cu(0) powder [START_REF] Whittaker | Synthesis of linear and 4-arm star block copolymers of poly(methyl acrylate-b-solketal acrylate) by SET-LRP at 25°C[END_REF]. All SET-LRP reveal first order kinetic and reach high conversions in a short time. Besides, the copper residue is easily removed, resulting in the obtaining of pure polymer. The molecular weight distribution is well controlled using the Cu(II)Br 2 /Me 6 -TREN catalytic complex. The hydrolysis of linear and star block copolymers afforded biocompatible amphiphilic P(MA-b-GA) n polymers. In aqueous solution, they form core-shell micelles and vesicles for the linear and star blocks, respectively, in which the hydrophilic block comes from the glyceryl side-groups. Consequently, the vesicle structure can be potentially used for drug delivery applications.

ATRP method has been broadly used for glycerol-based monomers, in particular for solketal acrylate, in order to obtain well-defined compositions, architectures, and molecular weights.

Many applications have been envisaged thanks to the architecture diversity. Nitroxide mediated polymerization (NMP) has been reported in the literature only for dihydroxypropyl acrylate (DHPA), and GMA. The kinetic studies of these hydroxyl functional monomers were investigated [START_REF] Popescu | Free radical and nitroxide mediated polymerization of hydroxy-functional acrylates prepared via lipase-catalyzed transacylation reactions[END_REF]. On the reverse, RAFT polymerization has not been considered.

To conclude, polymers based on glycerol derivatives have been well investigated in the literature using ATRP and NMP controlled radical polymerization. Different complex architectures were synthesized and polymers were obtained with a predictable molecular weight, and relatively narrow mass distribution indexes. In term of applications, the polymers bearing solketal pendant functionality exhibited amphiphilic and thermoresponsive behaviors.

They could be potentially applied in the biomedical field. Polymers bearing cyclocarbonate rings were used as low molecular weight polymer solvent for lithium ion conductivity.

Photopolymerization

The UV radiation is a simple and convenient energy which does not require expensive devices and requires low temperature and energy efficient processes [START_REF] Andrzejewska | Photopolymerization kinetics of multifunctional monomers[END_REF][START_REF] Oster | Photopolymerization of vinyl monomers[END_REF]. Thanks to its high output this special polymer processing is enjoying a new expansion and is applied at industrial scale for coating industry, paints or printing inks, adhesives, composite materials and dental restorative formulations and many others. Liquid resins can be converted into solid resins in a few tenths of a second making this process very attractive for the scientist community since almost three decades. The photopolymerization mechanism can be achieved either by polyaddition of double bonds under radical or cationic initiation. Glycerol (meth)acrylatebased monomers were involved in photopolymerization (Scheme 5). According to the target applications, these (meth)acrylate monomers bear one, two or three unsaturated groups for developing thermoplastic or thermoset materials. The impact of the polymer structure on the material properties is investigated as well as the effect of monomer structure on the polymerization kinetic. In order to explain the difference of reactivity of these glycerol derivatives in photopolymerization, the dipole moment, the intramolecular and the intermolecular effects that would alter the electron density on the radical and the double bond, the contribution of partially ionic intermediates will be evaluated.

Scheme 5. Photopolymerization of glycerol (meth)acrylate based monomers.

Six hydroxyl-bearing methacrylate monomers with various pendant chains bearing ether (GM-5), ester (GM-1, GM-2, GM-3, GM-4) or phosphine oxide group (GM-6) were successfully synthesized and used as reactive diluents in dental composites. The structurereactivity of the monomer relationship was also investigated [START_REF] Karahan | Development of reactive methacrylates based on glycidyl methacrylate[END_REF]. The photo-, homo-, and copolymerization behaviors of these monomers with (Bis-GMA) were studied to afford network polymers with appropriate mechanical properties for dental material applications.

Especially, the crosslinked polymer from the monomer GM-6 was reached as potential flame retardant thanks to its phosphine oxide groups. The influence of acetal groups on acrylate reactivity was also investigated by Kilambi et al. [START_REF] Kilambi | Influence of the secondary functionality on the radical-vinyl chemistry of highly reactive monoacrylates[END_REF]. Indeed, the high reactivity of several acrylates bearing 5-membered ring cycloacetal (5-CAA) and 6-membered ring acetal (6-CAA) in radical photopolymerization compared to hexyl acrylate was demonstrated. Jansen et al. also studied the photopolymerization of SolA and demonstrated the relationship between its low dipole moment and its low reactivity [START_REF] Jansen | Effect of Dipole Moment on the Maximum Rate of Photoinitiated Acrylate Polymerizations[END_REF].

Several teams are interested by the peculiar (meth)acrylate monomers bearing cyclocarbonate group including GCA, GCMA, CCCM and CGCA whose the structures are reported in Figure 1. The originality of these monomers is that they are not able of hydrogen bonding. Decker [START_REF] Moussa | Light-induced polymerization of new highly reactive acrylic monomers[END_REF][START_REF] Decker | A new class of highly reactive acrylic monomers, 2. Lightinduced copolymerization with difunctional oligomers[END_REF][START_REF] Decker | A new class of highly reactive acrylic monomers, 1. Lightinduced polymerization[END_REF] and Podszun [59,60] demonstrated the high reactivity of the acrylates bearing cyclocarbonate group relative to acrylates bearing a linear secondary functionality.

Additionally, Jansen et al. showed the direct relationship between the dipole moment and the rate polymerization [START_REF] Jansen | Effect of Dipole Moment on the Maximum Rate of Photoinitiated Acrylate Polymerizations[END_REF]. For instance, GCA with a high dipole moment is correlated to a high photopolymerization rate unlike SolA with a low dipole moment correlated to a low photopolymerization rate. Additionally, Berchtold et al. investigated the influence of the temperature on the cyclocarbonate-based acrylates [START_REF] Berchtold | Reactivity of Monovinyl (Meth)acrylates Containing Cyclic Carbonates[END_REF]. The latter proved to have little influence on the photopolymerization rate even if an increase is noted with this parameter.

Otherwise Beckel et al. [START_REF] Beckel | Evaluation of a Potential Ionic Contribution to the Polymerization of Highly Reactive (Meth)acrylate Monomers[END_REF] determined the anionic contribution related to the photopolymerization mechanism of (meth)acrylates in the presence of a strong acid (methanesulfonic acid). It could stabilize the radical species of the propagating (meth)acrylic chain leading to reduce the hydrogen abstraction reactions that contribute to the anionic characteristic. Kilambi et al. investigated the same influence of acid concentration in the photopolymerization rate of cyclocarbonate acrylates in the presence of highly polar solvent such as propylene carbonate [START_REF] Kilambi | Factors affecting the sensitivity to acid inhibition in novel acrylates characterized by secondary functionalities[END_REF].

Recently, Podgorski et al. proposed 1,3-glycerol dimethacrylate (GDM) as an alternative monomer for triethyleneglycol dimethylacrylate in dental mixtures. [START_REF] Maciej | Synthesis and characterization of acetyloxypropylene dimethacrylate as a new dental monomer[END_REF] This dimethylacrylate is a promising photocurable dental diluent due to its low viscosity, good mechanical and water uptake properties. Park [START_REF] Park | The influence of chemical structure on the properties in methacrylate-based dentin adhesives[END_REF] and Podgorski [START_REF] Podgórski | Structure-property relationship in new photo-cured dimethacrylatebased dental resins[END_REF] also reported the study of GDM as well as the tetramethacrylate (UR-3) [START_REF] Moszner | Synthesis and polymerisation of new multifunctional urethane methacrylates[END_REF]. As expected, the dimethacrylates have higher degree of double bond conversion related to the tetramethacrylate as a result of the higher density of the resulting network.

In summary, numerous studies have been realized on the theoretical aspects of the photopolymerization of the (meth)acrylate monomers based on glycerol to understand the influence of secondary function on polymerization rate. The rapid radical photopolymerization of these monomers is especially interesting for the elaboration of thermoset dental restoration applications as well as in adhesives, due to the presence of free hydroxyl groups in the glycerol part.

Polymerization by Acceptor -Donor process

Alkyl and vinyl ethers do not radically homopolymerize, but copolymerize readily with electron-deficient olefins such as fumarates or maleates. These facile reactions proceed through an electron donor-acceptor complex, which generates zwitterion or diradical tetramethylenes as initiating species. Alternating copolymers is a peculiar class of polymers obtained by the reaction between two monomers which cannot homopolymerize. One synthetic method consists in the polymerization between a monomer rich in electrons named donor and a monomer poor in electron named acceptor [START_REF] Braun | Polymers from non-homopolymerizable monomers by free radical processes[END_REF][START_REF] Rzaev | Complex-radical alternating copolymerization[END_REF][START_REF] Singh | Alternating copolymerization of donoracceptor monomers[END_REF]. Generally, this kind of copolymerization is initiated by a free radical or sometimes by a spontaneous reaction due to hydrogen transfer from donor to acceptor monomer (Scheme 6). The alternating copolymers between vinylSol and maleic anhydride was carried out by Crivello et al. [START_REF] Crivello | Transition metal-catalyzed isomerization and polymerization of allylic and enolic ethers[END_REF] Interestingly, a further crosslinking was performed leading to the formation of ester cross linkers. The hydrolysis of some anhydride groups permitted the formation of carboxylic acid groups which induced the hydrolysis of the ketal moieties. The resulting diols reacted with the residual anhydrides to provoke the ester crosslinking (Scheme 7.) The presence of both anhydride and acid sensitive protected diol groups in the same polymer thus makes this polymer a potential candidate for moisture curable coatings.

Scheme 7. Crosslinking of alternating copolymers based on vinylSol.

Ha et al. [START_REF] Ha | Synthesis of poly{[(2-oxo-1,3-dioxolan-4-yl)-methyl vinyl ether]-co-acrylonitrile} and its miscibility with SAN[END_REF] reported the miscibility of blends based on PVC and an alternating poly(vinylGC-co-AN) copolymer obtained by reaction between acrylonitrile and vinylGC. As a result, a blend of vinylGC and acrylonitrile (AN) is entirely miscible as the optical clarity, the thermal analysis, and the light scattering studies demonstrate it. Moreover, Ha et al. [START_REF] Ha | Synthesis of poly{[(2-oxo-1,3-dioxolan-4-yl)-methyl vinyl ether]-co-acrylonitrile} and its miscibility with SAN[END_REF] also described a blend of poly(vinylGC-co-AN)/SAN in order to develop an integrated process for the catalytic conversion of CO 2 to useful polymer materials. Their miscibility over the entire range of blend compositions was confirmed. Otherwise Moon et al. [START_REF] Moon | Synthesis of poly[(2-oxo-1,3-dioxolan-4-yl)methyl vinyl ether-co-N-phenylmaleimide] and its miscibility in blends with styrene-acrylonitrile or poly(vinyl chloride)[END_REF] investigated the synthesis of a copolymer bearing cyclic carbonate (vinylGC) and its miscibility with SAN or PVC. Interestingly, the vinylGC was synthesized from oxirane derivative using a catalytic conversion of carbon dioxide. The blended copolymer exhibited a higher thermal stability than poly(vinylGC).

To conclude, various polymers based on glycerol were synthesized via donor-acceptor polymerization process. This reaction is thermally or photochemically induced and conducted in solution as well as in liquid film of bulk reactants. Furthermore, the blended copolymer like poly(vinylGC-co-AN)/SAN was achieved in order to improve poly(vinylGC) thermal stability. Moreover, the potential moisture curable coatings based vinylSol was also described.

Conclusions and outlook

A variety of polymers containing glycerol derivatives can be synthesized via radical polymerization processes such as for instance free radical polymerization (FRP), controlled/living radical polymerization (CRP), photopolymerization as well as acceptordonor alternating polymerization. The FRP proved to be the simplest method in comparison with CRP in term of reaction conditions. Nevertheless only CRP technique allowed predictable molecular weight, well-defined compositions, and controlled architectures.

Thanks to this possibility, the obtained polymers can be employed in many different applications. Typically, the polymers containing acetal protecting functionality lead to amphiphilic structures, or thermoresponsive polymers. Hence, these kinds of polymer can find applications in the biomedical field such as drug delivery, for instance. Only ATRP and one example of NMP have been described in the literature for glycerol-based monomers, it would be a tremendous opportunity to use the RAFT process for this kind of monomers.

Furthermore, the photopolymerization is also a stimulant method which permitted a wide application in industry due to its main advantages: solvent free, energy efficient… However, most of the publications in this field are limited in fundamental studies. Only potential reported application using monomers multifunctionality bearing cyclocarbonate dealt with restorative dental materials to improve mechanical properties.

Figure 1 .

 1 Figure 1. Summary of monomers based on glycerol.

Figure 2 .

 2 Figure 2. Glycerol derivatives based-monomers involved in free radical polymerization.

Figure 3 .

 3 Figure 3. Monomers involved in the synthesis of statistical polymer with SolMA.

Figure 4 .

 4 Figure 4. PLLA-b-PSolA diblock copolymers prepared by atom transfer radical

Figure 5 .

 5 Figure 5. PSolA containing fluorescent segment brought by the difunctional initiator.

Figure 6 .

 6 Figure 6. PGMA-b-PPO-b-PGMA triblock copolymers from the ATRP of GMA.

Figure 7 .Scheme 4 .

 74 Figure 7. PSolA-b-PVA-b-PSolA triblock copolymers prepared by ATRP.

Figure 8 .

 8 Figure 8. Y-shaped PSolA-b-PEG-b-PSolA triblock copolymers prepared by ATRP.

Figure 9 .

 9 Figure 9. PSolA-based star copolymers prepared by ATRP.

  corresponding glycerol trimethacrylate (GTM) in dentin adhesives. The experimental adhesives exhibit significantly higher storage and rubbery modulus related to the widespread diethyleneglycol diacrylate photocrosslinker. During the study of hydroxyl-containing di(meth)acrylate monomers such as MCMA-D and ECMA-D, Avci et al. demonstrated the effect of the hydroxyl groups on monomer on the bulk reactivity [67]. The nature of the dimethacrylate group has little influence on the reactivity. Corresponding trimethacrylates (ECMA-T and MCMA-T) are less reactive than dimethacrylates. Moszner et al. successfully synthesized two urethane dimethacrylates based on glycerol (UR-1, UR-2) and a

Scheme 6 .

 6 Scheme 6. Synthesis of alternating copolymers by acceptor-donor polymerization.

Figure 10 .

 10 Figure 10. Alternating copolymers synthesized from vinylSol, vinylGC, and allylSol.


A Gramfort

D Strohmeier

J Haueisen

M Hämäläinen

M Kowalski

Time-Frequency Mixed-Norm Estimates: Sparse M/EEG imaging with non-stationary source activations

Keywords: Inverse problem, Magnetoencephalography (MEG), Electroencephalography (EEG), sparse structured priors, convex optimization, time-frequency, algorithms

Magnetoencephalography (MEG) and electroencephalography (EEG) allow functional brain imaging with high temporal resolution. While solving the inverse problem independently at every time point can give an image of the active brain at every millisecond, such a procedure does not capitalize on the temporal dynamics of the signal. Linear inverse methods (Minimumnorm, dSPM, sLORETA, beamformers) typically assume that the signal is stationary: regularization parameter and data covariance are independent of time and the time varying signal-to-noise ratio (SNR). Other recently proposed non-linear inverse solvers promoting focal activations estimate the sources in both space and time while also assuming stationary sources during a time interval. However such an hypothesis only holds for short time intervals. To overcome this limitation, we propose time-frequency mixed-norm estimates (TF-MxNE), which use time-frequency analysis to regularize the ill-posed inverse problem. This method makes use of structured sparse priors

Introduction

Distributed source models in magnetoencephalography and electroencephalography (collectively M/EEG) use thousands of current dipoles that are used as candidate sources to explain the M/EEG measurements. Those dipoles can be located on a dense three-dimensional grid within the brain volume, typically every 5 mm, or over a surface of the segmented cortical mantle [START_REF] Dale | Improved localization of cortical activity by combining EEG and MEG with MRI cortical surface reconstruction: a linear approach[END_REF], both of which can be automatically segmented from high-resolution anatomical Magnetic-Resonance Images (MRIs). Following Maxwell's equations, each dipole adds its contribution linearly to the measured signal. Note that this linearity of the forward problem is not a modeling assumption but a fact based on the fundamental physics of the problem.

The task in the inverse problem is to map the M/EEG measurements to the brain, i.e., to estimate the distribution of dipolar currents that can explain the measured data. Inverse methods that estimate distributed sources are commonly referred to as imaging methods. This is motivated by the fact that the current estimate explains the data and can be visualized as an image, at least at a given point in time. The orientations of the dipoles can be either considered to be known, e.g., by aligning them with the estimated cortical surface normals [START_REF] Dale | Improved localization of cortical activity by combining EEG and MEG with MRI cortical surface reconstruction: a linear approach[END_REF], in which case only the dipole amplitudes need to be estimated. Alternatively, the orientations can be considered as unknown in which case both amplitudes and orientations need to be estimated at each spatial location.

One of the challenges for distributed inverse methods is that the number of dipoles by far exceeds the number of M/EEG sensors: the problem is ill-posed. Therefore, constraints using a priori knowledge based on the characteristics of the actual source distributions are necessary. Common priors are based on the Frobenius norm and lead to a family of methods generally referred to as mininum norm estimators (MNE) [START_REF] Wang | Magnetic source images determined by a lead-field analysis: the unique minimum-norm least-squares estimation[END_REF][START_REF] Hämäläinen | Interpreting magnetic fields of the brain: minimum norm estimates[END_REF]. Minimum norm estimates can be converted into statistical parameter maps, which take into account the noise level, leading to noise-normalized methods such as dSPM [START_REF] Dale | Dynamic statistical parametric neurotechnique mapping: combining fMRI and MEG for high-resolution imaging of cortical activity[END_REF] or sLORETA [START_REF] Pascual-Marqui | Standardized low resolution brain electromagnetic tomography (sLORETA): technical details[END_REF]. While these methods have some benefits like simple implementation and a good robustness to noise, they do not take into account the natural assumption that only a few brain regions are typically active during a cognitive task. Interestingly, this latter assumption is what justifies a parametric method known as "dipole fitting" [START_REF] Scherg | Two bilateral sources of the late AEP as identified by a spatio-temporal dipole model[END_REF] routinely used in clinical practice. In order to promote such focal or sparse solutions within the distributed source model framework, one uses sparsity-inducing priors such as a p norm with p ≤ 1 [START_REF] Matsuura | Selective minimum-norm solution of the biomagnetic inverse problem[END_REF][START_REF] Gorodnitsky | Neuromagnetic source imaging with FOCUSS: a recursive weighted minimum norm algorithm[END_REF]. However, with such priors it is challenging to obtain consistent estimates of the source orientations [START_REF] Uutela | Visualization of magnetoencephalographic data using minimum current estimates[END_REF] as well as temporally coherent source estimates [START_REF] Ou | A distributed spatiotemporal EEG/MEG inverse solver[END_REF].

In order to promote spatio-temporally coherent focal estimates, several publications have proposed to constrain the active sources to remain the same over the time interval of interest [START_REF] Ou | A distributed spatiotemporal EEG/MEG inverse solver[END_REF][START_REF] Friston | Multiple sparse priors for the M/EEG inverse problem[END_REF][START_REF] Wipf | A unified Bayesian framework for MEG/EEG source imaging[END_REF][START_REF] Gramfort | Mixed-norm estimates for the M/EEG inverse problem using accelerated gradient methods[END_REF]. The implicit assumption is then that the sources are stationary. While this conjecture is reasonable for short time intervals, it is not a good model for realistic sources configurations where multiple transient sources activate sequentially during the analysis period, or simultaneously, before returning to baseline at different time instants.

When working with time series with transient and non-stationary effects, relevant signal processing tools are short time Fourier transforms (STFT) and wavelet decompositions. Contrary to a simple Fast Fourier Transform (FFT), they provide information localized in time and frequency (or scale). In particular, time-frequency decompositions, e.g., Morlet wavelet transforms, are routinely used in MEG and EEG analysis to study transient oscillatory signals. Such decompositions have been employed to analyze both sensorlevel data and source estimates, but no attempt has been made to use their output in constructing a regularizer for the inverse problem.

In this contribution, we address the problem of localizing non-stationary focal sources from M/EEG data using appropriate sparsity inducing norms. Extending the work from [START_REF] Gramfort | Mixed-norm estimates for the M/EEG inverse problem using accelerated gradient methods[END_REF] in which we coined the term Mixed-Norm Estimates (MxNE), we propose to use mixed-norms defined in terms of the time-frequency decompositions of the sources. We call this approach the Time-Frequency Mixed-Norms Estimates (TF-MxNE). The benefit is that the estimates can be obtained over longer time intervals while making stan-dard preprocessing such as filtering or time-frequency analysis on the sensors optional. The inverse problem is formulated as in [START_REF] Gramfort | Mixed-norm estimates for the M/EEG inverse problem using accelerated gradient methods[END_REF] as a convex optimization problem whose solutions are computed with an efficient solver based on proximal iterations.

We start with a detailed presentation of the problem and the algorithm. Next, we compare the characteristics and performance of various priors with help of realistic simulated data. Finally, we analyze publicly available MEG datasets (auditory and visual stimulations) demonstrating the benefit of TF-MxNE in terms of source localization and estimation of the time courses of the sources.

A preliminary version of this work was presented at the international conference on Information Processing in Medical Imaging (IPMI) [START_REF] Gramfort | Functional brain imaging with m/eeg using structured sparsity in time-frequency dictionaries[END_REF]. In this paper we improve the solver to support loose orientation constraints, depth compensation as well as a debiasing step to better estimate source amplitudes. We also analyze new experimental data.

Notation:

We indicate vectors with bold letters, a ∈ R N (resp. C N ) and matrices with capital bold letters, A ∈ R N ×N (resp. C N ×N ). a[i] stands for the i th entry in the vector, while A[i, •] and A[•, i] denote the i th row and i th column of a matrix, respectively. We denote A Fro the Frobenius norm,

A 2 Fro = N i,j=1 |A[i, j]| 2 , A 1 = N i,j=1 |A[i, j]|

General model and method

After a short introduction to Gabor time-frequency dictionaries for M/EEG signals, we present the details of our TF-MxNE inverse problem approach. We then detail the proposed optimization strategy, which uses proximal iterations.

Gabor dictionaries

Here we briefly present some important properties of Gabor dictionaries, see [START_REF] Daubechies | Ten lectures on Wavelets[END_REF] for more details. Given a signal observed over a time interval, its conventional Fourier transform estimates the frequency content but loses the time information. To analyze the evolution of the spectrum with time and hence the non-stationarity of the signal, Gabor introduced windowed Fourier atoms which correspond to a short-time Fourier transform (STFT) with a Gaussian window. In practice, for numerical computation, a challenge is to properly discretize the continuous STFT. The discrete STFT with a Gaussian window is also known as the discrete Gabor Transform [START_REF] Gabor | Theory of communication[END_REF].

The setting we are considering is the finite-dimensional one. Let g ∈ R T be a "mother" analysis window. Let f 0 ∈ N and k 0 ∈ N be the frequency and the time sampling rate in the time-frequency plane generated by the STFT, respectively. The family of the translations and modulations of the mother window generates a family of Gabor atoms (φ mf ) mf forming the dictionary Φ ∈ C T ×K , where K denotes the number of atoms. The atoms can be written as

φ mf [n] = g[n-mk 0 ]e i2πf 0 f n T , m ∈ {0, . . . , T k 0 -1}, f ∈ {0, . . . , T f 0 -1} . ( 1
)
If the product f 0 k 0 is small enough, i.e., the time-frequency plane is sufficiently sampled, the family (φ mf ) mf is a frame of R T , i.e., one can recover any signal x ∈ R T from its Gabor coefficients ( x, φ mf ) = Φ H x. More precisely, there exists two constants A, B > 0 such that [START_REF] Balian | Un principe d'incertitude en théorie du signal ou en mécanique quantique[END_REF]:

A x 2 2 ≤ m,f | x, φ mf | 2 ≤ B x 2 2 . (2)
When A = B, the frame is tight. When the vectors φ mf are normalized the frame is an orthogonal basis if and only if A = B = 1. The Balian-Low theorem says that it is impossible to construct a Gabor frame which is a basis. Consequently, a Gabor transform is redundant or overcomplete and there exists an infinitely number of ways to reconstruct x from a given family of Gabor atoms. In the following, the considered Φ dictionaries are tight frames. The canonical reconstruction of x from its Gabor coefficients requires a canonical dual window, denoted by g. Following (1) to define ( φmf ) mf we have:

x = m,f x, φ mf φmf = m,f x, φmf φ mf = Φ H x Φ = ΦH xΦ ,
where Φ is the Gabor dictionary formed with the dual windows. When the frame is tight, then we have g = g, and more particularly we have ΦΦ H = ΦΦ H Id 8 . The representation being redundant, for any x ∈ R T one can find a set of coefficients z mf such that x = m,f z mf φ mf , while the z mf verify some suitable properties dictated by the application. For example, it is particularly interesting for M/EEG to find a sparse representation of the signal. Indeed, a scalogram, sometimes simply called TF transform of the data in the MEG literature, generally exhibits a few peaks localized in the time-frequency domain. In other words, an M/EEG signal can be expressed as a linear combinations of a few oscillatory atoms. In order to demonstrate this, Fig. 1 shows the STFT of a single planar gradiometer channel MEG signal from a somatosensory experiment, the same STFT restricted to the 50 largest coefficients (approximately only 10% of the coefficients), and the signal reconstructed with only these coefficients compared to the original signal. We observe that the true signal can be well approximated by only a few coefficients, i.e., a few Gabor atoms. In the presence of white Gaussian noise, restricting the time-frequency representation of a signal to the largest coefficients denoises the data. This stems from the fact, that Gaussian white noise in not sparse in the time-frequency domain, but rather spreads energy uniformly over all time-frequency coefficients [START_REF] Trujillo-Barreto | Bayesian M/EEG source reconstruction with spatio-temporal priors[END_REF]. Thresholding or shrinking the coefficients therefore reduces noise and smoothes the data. This is further explained in the context of wavelet transforms in [START_REF] Donoho | De-noising by soft-thresholding[END_REF]. In practice, the Gabor coefficients are computed using the Fast Fourier Transform (FFT) and not by a multiplication by a Φ matrix as suggested above. Such operations can be efficiently implemented as in the LTFAT toolbox9 [START_REF] Soendergard | The linear time frequency toolbox[END_REF]. Another practical concern to keep in mind is the tradeoff between the size of the window g and the time shift k 0 . A long window will have a good frequency resolution and a limited time resolution. The time resolution can be improved with a small time shift leading however to a larger computational cost, both in time and memory. Finally, as any computation done with an FFT, the STFT implementations assume circular boundary conditions for the signal. To take this into account and avoid edge artifacts, the signal has to be windowed, e.g., using a Hann window.

The inverse problem with time-frequency dictionaries

The linearity of Maxwell's equations implies that the signals measured by M/EEG sensors are linear combinations of the electromagnetic fields produced by all current sources. The linear forward operator, called gain matrix, predicts the M/EEG measurements due to a configuration of sources based on a given volume conductor model [START_REF] Mosher | EEG and MEG: Forward solutions for inverse methods[END_REF]. Given such a linear forward operator G ∈ R N ×P , where N is the number of sensors and P the number of sources, the measurements M ∈ R N ×T (T number of time instants) are related to the source amplitudes X ∈ R P ×T by M = GX .

The computation of the gain matrix G, e.g., with a Boundary Element Method (BEM) [START_REF] Kybic | A common formalism for the integral formulations of the forward EEG problem[END_REF][START_REF] Gramfort | OpenMEEG: opensource software for quasistatic bioelectromagnetics[END_REF], requires modeling of the electromagnetic properties of the head [START_REF] Hämäläinen | Interpreting magnetic fields of the brain: minimum norm estimates[END_REF] such as the specification of the tissue conductivities. The matrix is then numerically computed. In the inverse problem one computes a best estimate of the neural currents, X , based on the measurements M. However, since P N , the problem is ill-posed and priors need to be imposed on X. Historically, the sources amplitudes were computed time instant by time instant using priors based on p norms. The 2 (Frobenius) norm leads to MNE, LORETA, dSPM, or sLORETA while several alternative solvers based on p norms with p ≤ 1 have also been proposed to promote sparse solutions [START_REF] Matsuura | Selective minimum-norm solution of the biomagnetic inverse problem[END_REF][START_REF] Gorodnitsky | Neuromagnetic source imaging with FOCUSS: a recursive weighted minimum norm algorithm[END_REF]. However, since such solvers work on an instant by instant basis they do not model the oscillatory nature of electromagnetic brain signals. Note that even if the 2 norm based methods work time instant by time instant, the estimates reflect the temporal characteristics of the data, since they are obtained by linear combinations of sensor data. This, however, implies that the parameters of the inverse solver are independent of time, which corresponds to assuming that the SNR is independent of time. Although MNE type approaches have been used with success, the assumption of constant SNR is clearly wrong since the signal amplitudes vary in time while the noise stays constant, or may be even smaller during an evoked response. The noise is usually estimated from baseline periods such as prestimulus intervals or periods when the brain is not yet responding to the stimulus.

Beyond single instant solvers, various sparsity-promoting approaches have been proposed [START_REF] Ou | A distributed spatiotemporal EEG/MEG inverse solver[END_REF][START_REF] Friston | Multiple sparse priors for the M/EEG inverse problem[END_REF][START_REF] Wipf | A unified Bayesian framework for MEG/EEG source imaging[END_REF]. Although they manage to capture the time courses of the activations more accuratly than the instantaneous sparse solvers, they implicitly assume that the active sources are the same over the entire time interval of interest. This also implies that if a source is detected as active at one time point, its activation will be non-zero during the entire time interval of interest. To go beyond this approach, we propose a solver which promotes on the one hand that the source configuration is spatially sparse, and on the other hand that the time course of each active dipole is a linear combination of a limited number of Gabor atoms, as suggested by Fig. 1. Since a Gabor oscillatory atom is localized in time, sources can be marked as active only during a short time period. The model reads:

M = GX + E = GZΦ H + E , (3)
where Φ H ∈ C K×T is a dictionary of K Gabor atoms, Z ∈ C P ×K are the coefficients of the decomposition, and E is additive white noise, E ∼ N (0, λI). Given a prior on Z, P(Z) ∼ exp(-Ω(Z)), the maximum a posteriori estimate (MAP) is obtained by solving:

Z = arg min Z 1 2 M -GZΦ H 2 Fro + λΩ(Z) , λ > 0 . (4)
If we consider Ω(Z) = Z 1 , (4) corresponds to a LASSO problem [START_REF] Tibshirani | Regression shrinkage and selection via the lasso[END_REF], a.k.a. Basis Pursuit Denoising (BPDN) [START_REF] Chen | Atomic decomposition by basis pursuit[END_REF], where features (or regressors) are spatio-temporal atoms. Similarly to the original formulation of MCE, i.e., 1 regularization without applying Φ, such a prior is likely to suffer from inconsistencies over time [START_REF] Ou | A distributed spatiotemporal EEG/MEG inverse solver[END_REF]. Indeed such a norm does not impose a structure for the non-zero coefficients: they are likely to be scattered all over Z (see Fig. 2). Therefore, simple 1 priors do not guarantee that only a few sources are active during the time window of interest. To promote this, one needs to employ mixed-norms such as the 21 norm [START_REF] Ou | A distributed spatiotemporal EEG/MEG inverse solver[END_REF][START_REF] Gramfort | Mixed-norm estimates for the M/EEG inverse problem using accelerated gradient methods[END_REF]. By doing so, the estimates have a sparse row structure (see Fig. 2). However the sum of 21 and 1 norms. The prior then reads:

λΩ(Z) = λ space Z 21 + λ time Z 1 , λ space > 0, λ time > 0 . (5)
A large regularization parameter λ space will lead to a spatially very sparse solution, while a large regularization parameter λ time will promote sources with smooth times series. This is due to the uniform spectrum of the noise (see Section 2.1) and the fact that a large λ time will promote source activations made up of few TF atoms, each of which has a smooth waveform.

Optimization strategy

The optimization strategy, which we propose for minimizing the cost function in [START_REF] Chen | Atomic decomposition by basis pursuit[END_REF], is based on the Fast Iterative Shrinkage Thresholding Algorithm (FISTA) [START_REF] Beck | A fast iterative shrinkage-thresholding algorithm for linear inverse problems[END_REF], a first-order schemes that handles the minimization of any cost function F that can be written as a sum of two terms: a smooth convex term f 1 with Lipschitz gradient and a convex term f 2 , potentially non-differentiable: F(Z) = f 1 (Z) + f 2 (Z). In order to apply FISTA, we need to be able to compute the so-called proximity operator associated with f 2 , i.e., the proximity operator associated with the composite 21 + 1 prior [START_REF] Gramfort | Functional brain imaging with m/eeg using structured sparsity in time-frequency dictionaries[END_REF].

Definition 1 (Proximity operator). Let ϕ : R M → R be a proper convex function. The proximity operator associated to ϕ, denoted by prox ϕ : R M → R M reads:

prox ϕ (Z) = arg min V∈R M 1 2 Z -V 2 2 + ϕ(V) .
While the proximity operators of mixed-norms relevant for M/EEG can be found in [START_REF] Gramfort | Mixed-norm estimates for the M/EEG inverse problem using accelerated gradient methods[END_REF], in the case of the composite prior in [START_REF] Cottereau | Phase delays within visual cortex shape the response to steady-state visual stimulation[END_REF], the proximity operator is given by the following lemma.

Lemma 1 (Proximity operator for

21 + 1 ). Let Y ∈ C P ×K be indexed by a double index (p, k). Z = prox λ . 1 +µ . 21 ) (Y) ∈ C P ×K is given for each coordinates (p, k) by Z[p, k] = Y [p, k] |Y [p, k]| (|Y [p, k]| -λ) + 1 - µ k (|Y [p, k]| -λ) 2 + +
.

where for x ∈ R, (x) + = max(x, 0) , and by convention 0 0 = 0 . This result is a corollary of the proximity operator derived for hierarchical group penalties recently proposed in [START_REF] Jenatton | Proximal methods for hierarchical sparse coding[END_REF]. The penalty described here can indeed be seen as a 2-level hierarchical structure, and the resulting proximity operator reduces to successively applying the 1 and 21 proximity operator. Both of these proximity operators are discussed in detail in [START_REF] Gramfort | Mixed-norm estimates for the M/EEG inverse problem using accelerated gradient methods[END_REF].

The pseudo code is provided in Algorithm 1. The Lipschitz constant L of the gradient of the smooth term in ( 4) is given by the square of the spectral norm of the linear operator Z → GZΦ H . We estimate it with the power iteration method.

Z o = Z 6: Z = prox µλΩ Y + µG T (M -GYΦ H )Φ 7: τ o = τ 8: τ = 1+ √ 1+4τ 2 2
9:

Y = Z + τo-1 τ (Z -Z o ) 10: end for

Specific modeling for M/EEG inverse problem

The M/EEG literature has shown that general solvers of the statistics literature need to be adapted to the specificities of the M/EEG inverse problem.

Crucial steps in the computation of the source estimates are noise whitening, depth compensation, handling of source orientations, and amplitude bias correction.

Spatial whitening

The model in (3) assumes that the additive noise is Gaussian white with E ∼ N (0, λI). This strong modeling assumption is made realistic by a whitening step that relies on estimating the noise covariance matrix. For this purpose, baseline data is employed, which is recorded while the subject is at rest e.g. during pre-stimulus periods. If only MEG is recorded, the noise covariance can be estimated from data recorded without subject, often called empty room data. This approach provides good estimates of the measurement noise level. Although the noise level depends on the signal frequency, one usually uses a single frequency-unspecific noise covariance matrix. An alternative approach for frequency-dependent spatial whitening is presented in [START_REF] Ramírez | Spectral signal space projection algorithm for frequency domain MEG and EEG denoising, whitening, and source imaging[END_REF].

The whitening step is particularly fundamental when different sensor types are used: EEG and MEG with gradiometers and magnetometers record signals with different units of measure and with different noise levels. The whitening step makes data recorded by different sensors comparable and adapted for joint estimation.

Source models with unconstrained orientations

When the source orientations given by the normals of the cortical mesh cannot be trusted, it can be interesting to relax this constraint by placing three orthogonal sources at each spatial location. When all three orientations are allowed to explain the data equivalently the model is called free orientation. Moreover, it can be of interest to have intermediate models using loose orientation constraints [START_REF] Lin | Distributed current estimates using cortical orientation constraints[END_REF]. However for such loose and free orientation models, the TF composite prior needs to be adapted. Let each source be indexed by a spatial location i and an orientation o ∈ {1, 2, 3}. Let o = 1 correspond to the orientation normal to the cortex, and o = 2 and o = 3 the two tangential orientations. We call 0 < ρ ≤ 1, the parameter controlling how loose is the orientation constraint. The 1 and 21 norms read:

Z 1 = i,k |Z[(i, o = 1), k]| 2 + 1 ρ |Z[(i, o = 2), k]| 2 + 1 ρ |Z[(i, o = 3), k]| 2 Z 21 = i k |Z[(i, o = 1), k]| 2 + 1 ρ |Z[(i, o = 2), k]| 2 + 1 ρ |Z[(i, o = 3), k]| 2 ,
where k indexes the TF coefficients. When ρ = 1 the orientation is free and it amounts to grouping the orientations in a common 2 norm such as in [START_REF] Ou | A distributed spatiotemporal EEG/MEG inverse solver[END_REF][START_REF] Haufe | Combining sparsity and rotational invariance in EEG/MEG source reconstruction[END_REF]. Such priors are a principled way of supporting loose orientation constraints in the context of non-2 priors.

Observe here that Z 1 is not an 1 norm per se. Indeed, it is a 21 norm, but we have chosen to keep the same notations as in the constrained orientation case for the sake of readability.

In practice, using free orientation models means that at a given location, the current dipoles selected to explain the data can have an orientation that varies in time similarly to the rotating dipole model employed in dipole fitting.

Depth compensation

The principal contribution to M/EEG data comes from superficial cortical gray matter: deep sources are attenuated due to their larger distance from the sensors. While it is common in statistics to scale the columns of the gain matrix such that G i 2 = 1, practice with M/EEG data shows that it is often not a good idea. The rationale in statistics is to avoid favoring regressors, here sources, just due to the amplitude of the corresponding column in the gain matrix. When doing this for M/EEG, it tends to favor too much very deep sources which are less likely to be visible with M/EEG. For this reason, a common practice with MNE type approaches is to use a softer depth bias compensation. Given a parameter 0

≤ γ ≤ 1, the three columns (G[•, (i, o = 1)], G[•, (i, o = 2)], G[•, (i, o = 3)]) of G for the three orientations at the same lo- cation are normalized by ( G[•, (i, o = 1)] 2 2 + G[•, (i, o = 2)] 2 2 + G[•, (i, o = 3)] 2
2 ) γ . If γ = 0 it corresponds to no depth bias compensation and γ = 1 leads to full scaling which may lead to spurious deep sources appearing in the results.

Source weighting: fMRI priors?

Mixed-norm regularizations [START_REF] Gramfort | Mixed-norm estimates for the M/EEG inverse problem using accelerated gradient methods[END_REF] can be written with spatially dependent scalar weights. It can be used to promote some sources by reducing their regularization. For example, given a weight vector w ∈ R P + , one weight per physical location, the TF-MxNE prior can be modified as:

λΩ w (Z) = P i=1 w[i] (λ space Z[i, •] 2 + λ time Z[i, •] 1 ) ,
where Z[i, •] stands for the ith row of Z. If w[i] is small the regularization for the source at location i will be small and the source i is likely to be selected to explain the data. Assuming that additional location information of the sources is available, e.g., from fMRI, information is known about the sources, such as fMRI localizations, it would be possible to inject this knowledge in the prior in order to have fMRI informed sparse estimates. Note that sparsity promoting priors do not lead to source estimates where every dipole in the source space has a non-zero activation. It means that, although some regions are promoted by the weights, they may not contain any estimated source. It indeed may happen that MEG misses sources, for example if they are radially oriented. In this sense the proposed weighted scheme does not act as a strong prior on the MEG source localization.

For computational reasons, one can also exploit fast solvers such as dSPM or sLORETA to derive scalar weights that can help reduce the number of candidate sources. Typically, one can threshold dSPM/sLORETA estimates and restrict the TF-MxNE solver to a small portion of the cortex, further improving the computational efficiency of the optimization algorithm. It corresponds to setting w[i] to infinity (or very large) if the i th spatial location yields very low dSPM values at all points in time.

Amplitude bias compensation

Methods based on 1 priors, such as TF-MxNE, are known to impose an amplitude bias on the solution. This is due to the general bias-variance tradeoff in statistical estimation. With 1 based priors, the high sparsity of the solution comes at the price of a strong amplitude bias. Given the waveforms for the selected sources it is possible to post-process them and correct the amplitude bias leading to meaningful amplitudes of the source activations. See [START_REF] Grosenick | A family of interpretable multivariate models for regression and classification of whole-brain fMRI data[END_REF] for an example of amplitude bias correction in the context of fMRI decoding.

A first natural approach to correct for the amplitude bias is to compute the least squares solution restricted to the active set of sources provided by the TF-MxNE solution. It amounts to computing a dipole fit with a known set of dipoles, which is no longer an ill-posed problem. However, this procedure affects the source time courses, and the signal smoothness promoted by the TF-MxNE is lost. Hence, rather than re-estimating the source time courses using least squares, we correct the amplitude bias by scaling the TF-MxNE results. For this purpose, we introduce a diagonal scaling matrix D, whose diagonal elements are scaling factors for all sources in the active set. These scaling factors are constrained to be above 1 to actually remove the bias, and are constant over time. Furthermore, in the case of free orientation, they are identical for all orientations at a given location in order to preserve the source characteristics and orientations estimated using TF-MxNE. The bias corrected source estimate X is computed using D as X = DX = DZΦ H . We estimated the scaling matrix D based on the following convex optimization problem:

D = arg min D M -GDX 2
Fro s.t.

D ij ≥ 1, i = j D ij = 0, i = j .
The optimization problem can also be solved efficiently with FISTA after writing the constraint on D as an indicator function over a convex set C = {D s.t. D ii ≥ 1, and

D ij = 0, if i = j}: f 2 (D) = ι C (D) = 0 if D ∈ C ∞ otherwise .

Practical details

This section presents the details in the efficient implementation of Algorithm 1. We also discuss the choice of the hyperparameters (regularization parameters).

Implementation

Algorithm 1 requires to compute Gabor transforms at each iteration which can be computationally demanding. However, due to the 21 sparsity inducing prior, only a few rows of Z have non-zero coefficients. The Gabor transform is therefore computed for only a limited number of rows, equivalently a small number of active sources. This makes the computation of YΦ H (cf. Algorithm 1 line 6) much faster.

Also when a tight frame is used, the 21 norm of a signal does not change when Φ is applied. This means that the 21 proximity operator can be applied to temporal data to discard some sources from the active set without computing the STFT. This comes from the fact that if prox 21 (x) = 0 for a time series x then prox 1 + 21 (Φx) = 0.

Since the proposed optimization problem is convex, the solution does not depend on initial conditions. Hence, in order to further reduce the computation time, it is beneficial to initialize the TF-MxNE solver with the 21 MxNE solution obtained with the same spatial regularization, since MxNE can be computed efficiently using active set strategies [START_REF] Gramfort | Mixed-norm estimates for the M/EEG inverse problem using accelerated gradient methods[END_REF]. Note again that the 21 MxNE solution is used as an initialization and not for restricting the source space.

Selection of the regularization parameters

Model selection in the present case amounts to setting the regularization parameters λ space and λ time , as well as the parameter of the Gabor transform, namely the time resolution with k 0 and the frequency resolution, function of the window length T . The parameter k 0 and T will depend on the length of a time interval during which signals can be considered stationary. A too high sampling of the time-frequency plane will also lead to high computational costs. The regularization parameters have an effect on the spatial sparsity, the number of active dipoles, and the temporal smoothness of the source time series. Different strategies exist to set such model parameters (crossvalidation, discrepancy principle etc.).

In the case of 21 priors, one can prove that there exists a value of λ max space for λ space such that if λ space ≥ λ max space , then Z * is filled with zeros, i.e., no source is active. This provides a convenient way to specify the regularization parameter as the ratio of λ space and λ max space , between 0 and 1. In the next section if λ space is given as a percentage it corresponds to this ratio, rescaled to percents. For convenience, the parameter λ time can then be also scaled by λ max space . The benefit of the reparametrization of the regularization parameter is that they become much less sensitive to the dataset. Assuming Φ is a tight frame then X 21 = XΦ 21 = Z 21 , then one can show based on the optimality conditions for the 21 mixed-norm [START_REF] Gramfort | Mixed-norm estimates for the M/EEG inverse problem using accelerated gradient methods[END_REF] that:

λ max space = max i (G T M)[i, •] 2

Results

In the following, we first evaluate the accuracy of our solver with simulations. We then apply our solver to two MEG/EEG datasets.

Simulation study

In order to have a reproducible and reasonably fast comparison of different priors, we generated a small simulation dataset with 20 EEG electrodes and 200 sources. Four of these sources were randomly selected to be active. The ECD waveforms (Fig. 4(a)) represent 1 high and 3 low frequency components. The time course of the oscillatory high frequency component is modeled by a Gabor atom, whereas the time courses of the low frequency components were obtained from a somatosensory evoked potential study [START_REF] Jaros | Nonlinear interactions of high-frequency oscillations in the human somatosensory system[END_REF] by fitting manually ECDs to the P15, N20 and P23 components. To make the comparison of the priors independent of the forward model and the sources spatial configuration, the linear forward operator was a random matrix, whose columns were normalized to 1. The scaling of the columns is not mandatory here but simplifies the parameter setting by making it independent of the number of simulated sensors. White Gaussian noise or realistic 1/f noise, simulated with an auto-regressive (AR) process of order five calibrated on real MEG data, was added to the signals to achieve a desired signal-to-noise ratio (SNR). Following the notation of (3), we define SNR as 20 log 10 ( M Fro / E Fro ).

Figure 3 presents the RMSE on the estimation for different solvers as a function of λ space (RMSE = X sim -X Ω 2 Fro ) for two different values of λ time , with and without correcting for the amplitude bias. High and low noise conditions were investigated, with both white and AR noise. The parameter λ space was chosen on a logarithmic grid from 10 -2 to 10 2.5 . The Gabor dictionary is tight, constructed with a 64 sample window g with k 0 = 4 samples time shift.

Many observations can be made from these results. First, the spatiotemporal solvers, namely MxNE and TF-MxNE, outperform instantaneous MNE and MCE. Second, we observe that a small value of λ time in TF-MxNE leads to results similar to MxNE which is fully expected since we use a tight frame that does not change the 2 norm of a signal. Finally, we observe that correcting for the amplitude bias clearly improves the results in all conditions and particularly with high values of λ space . This suggests to run the solver with a high value of λ space (high percentage of λ max space ) which also leads in practice to a faster convergence of the solver. Note also that, as expected, violation of the modeling assumption about the noise, namely Gaussian white noise, produces a degradation in the performance of the solver in the AR condition. Figure 4 shows the reconstructions for the best λ space according to Figure 3 for the 1 , 21 MxNE and the TF-MxNE (with and without bias correction). It can be observed, that the TF-MxNE method with the composite TF prior is able to reconstruct the smooth time course of the simulated sources contrary to 1 and 21 priors.

The TF composite prior was then challenged on a realistic EEG configuration with a 4-shell spherical head model (radii 0.90, 0.92, 0.96 and 1) and 60 electrodes placed according to the international 10-5 electrode system. The source waveforms were the same as before. The source space in Fig. 5 consisted of 152 sources in a regular grid (0.2 spacing) inside the innermost sphere. Source orientations were randomly selected. Figure 5 shows the head model and reconstructions obtained with 21 MxNE and the TF-MxNE. Even if the performance drops down due to the limited spatial resolution of EEG, the TF composite prior gives the best RMSE and is able to reconstruct and separate the high frequency component.

Experimental results with MEG/EEG data

We also applied the TF-MxNE solver to publicly available data: the auditory and visual conditions in the data shipped as sample data with the MNE software (http://martinos.org/mne/). In this MEG/EEG experiment, checkerboard patterns were presented into the left and right visual field, interspersed by tones to the left or right ear. The interval between the stimuli was 750 ms. Occasionally a smiley face was presented at the center of the visual field. The subject was asked to press a key with the right index finger as soon as possible after the appearance of the face. We will report results obtained on the responses evoked by the auditory stimuli presented to the left ear and visual stimuli shown in the left hemifield.

Data were collected in a magnetically shielded room using the wholehead Elekta Neuromag Vector View 306 MEG system (Neuromag Elekta LTD, Helsinki) equipped with 102 triple-sensor elements (two orthogonal planar gradiometers and one magnetometer per location). EEG data from a 60-channel electrode cap was acquired simultaneously with the MEG.

The signals were recorded with a bandpass of 0.1 -172 Hz, digitized at 600 samples/s and averaged offline triggered by the stimulus onset. All epochs containing EOG signals higher than 150 µV peak-to-peak amplitude were discarded from the averages. Peak-to-peak rejection parameters for EEG was set to 80 µV, for magnetometers 4000 fT and for gradiometers 2000 fT/cm. This resulted in respectively 55 and 67 averaged epochs for the left auditory and left visual conditions. Evoked data were baseline corrected using 200 ms of pre-stimulus data. The same data segment was used to estimate the noise covariance matrix for spatial whitening. The M/EEG recordings contained two bad channels (1 EEG and 1 MEG), leading to a total of 364 combined M/EEG channels. Signal-space projection (SSP) correction was applied to MEG magnetometers data to suppress environmental noise and biological artifacts [START_REF] Uusitalo | Signal-space projection method for separating meg or eeg into components[END_REF].

The anatomical MRI data were collected with a Siemens Trio 3T scanner with a T1-weighted sagittal MPRAGE protocol, which were employed for cortical surface reconstruction using FreeSurfer. Two multi-echo 3D Flash acquisitions were also performed to extract the inner skull surface for the 3-layers boundary-element model used for forward model computation. The dipolar sources were sampled over the cortical mantle with 5 mm average distance between sources leading to a total of 7498 dipoles.

Auditory data

The sources were estimated assuming loose orientation constraint (ρ = 0.2) and depth bias compensation with γ = 0.9. The Gabor dictionary was tight, constructed with a 16 samples ( 27 ms) long window g with k 0 = 4 samples time shift and f 0 = 1 sample frequency shift. A scalar weighting as described in Section 3.4, was performed with a dSPM solution obtained with the same depth weighting and loose orientation. The λ time parameter was set to 1%. Results are presented in Figure 6.

In order to compare the estimated sources of the TF-MxNE, we computed the solution using an 21 MxNE with the same spatial regularization. It amounts to setting λ time in the TF-MxNE to zero and hence not promoting any temporal regularization. This can be observed in the Figure 6 where both MxNE and TF-MxNE lead to bilateral auditory sources given an λ space = 50%. However the TF-MxNE, leads to smooth time courses and true zero activations during baseline. Note that any sparse solver that would only promote spatial sparsity without modeling the dynamics of the signal would fail to estimate true null activations during baseline.

Further looking at the time series of the sources, it can be observed that the contralateral auditory cortex activates before the ipsilateral. This is consistent with the literature. The peak to peak latency difference between right and left cortices for the main activity, known as the N100, can be clearly quantified on the smooth waveforms provided by the TF-MxNE. It is for this subject equal to 16 ms. These experimental results show how the proposed method can provide both fine spatial localization and temporal smoothness, which e.g. allows to fully exploit the temporal resolution of M/EEG for chronometry.

In order to illustrate, what fraction of the data has actually been explained by the bilateral auditory sources obtained with our choice of hyperparameters λ space and λ time , we present in Figure 6-a the original data restricted to the gradiometers as well as in Figure 6-b the data predicted by the sources. Note that all MEG and EEG sensors were used for the localization. However, we show only the gradiometer signals due to visualization purposes. We observe that with the regularization parameters used, the evoked components predicted by the sources are mostly before 150 ms, including the P50 and N100.

For comparison with alternative solvers, we show in Figure 7 the TF-MxNE source locations superposed to dSPM reconstructions and LCMV beamformer [START_REF] Veen | Localization of brain electrical activity via linearly constrained minimum variance spatial filtering[END_REF] outputs 100 ms after stimulation. Since LCMV does not exploit the cortical orientation information, the map peaks are located at the gyri, below the auditory cortex in the temporal lobe. The LCMV solution was here computed using only the left temporal channels in order to avoid running a beamformer on correlated sources, here coming from bilateral auditory cortices.

Visual data

We used the same model parameters (loose orientation, depth bias, scalar weighting, λ time and Gabor dictionary) as for the auditory condition. The spatial regularization was however changed to λ space = 30% of λ max space . Results are presented in Figure 8.

Figure 8(a) presents the raw evoked response, restricted to the gradiometers. Sources reconstructions lead to three dipoles. According the automatic parcellation of the cortex provided by FreeSurfer, two sources are localized in the early visual cortex V1, while a third one is positioned on the dorsal part of the secondary visual cortex (V2). One can then quantify the temporal latencies. Between the largest peaks for the early V1 sources the latency is zero which suggests that it could be the same source not properly modeled with a single cortically constraint dipole. Between the V1 sources and the later V2 source the latency is equal to 7 ms, which is very reasonable according to the literature [START_REF] Cottereau | Phase delays within visual cortex shape the response to steady-state visual stimulation[END_REF]. Here again one can observe the source time series are truely set to zero during baseline.

Using the implementation details provided in Section 4.1, the computation on each dataset presented in this section takes about 30 seconds on a standard laptop computer.

For comparison with alternative solvers, we show in Figure 9 the TF-MxNE source locations superposed to dSPM reconstructions and LCMV beamformer [START_REF] Veen | Localization of brain electrical activity via linearly constrained minimum variance spatial filtering[END_REF] outputs at the time instants corresponding to the two peaks in the evoked response. For LCMV, the data covariance was computed from 40 ms to 150 ms. One can observe that for the early peak all three methods agree to locate the main source in V1 along the calcarine fissure, although dSPM and LCMV lead to more smeared maps. For the second peak, it is pleasing to see how well LCMV and TF-MxNE agree on the site of the blue source, located in the dorsal part of V2 according to the FreeSurfer parcellation.

Discussion

While time-frequency analysis is commonly used in the context of M/EEG both in the sensor and source space, it has rarely been used to better model the source dynamics in the context of the inverse problem. Some earlier contributions, such as [START_REF] Durka | Multichannel matching pursuit and EEG inverse solutions[END_REF][START_REF] Geva | Spatio-temporal matching pursuit (SToMP) for multiple source estimation of evoked potentials[END_REF][START_REF] Lelic | Inverse modeling on decomposed electroencephalographic data: A way forward[END_REF][START_REF] Lina | Wavelet-based localization of oscillatory sources from magnetoencephalography data[END_REF], apply a two-step approach. First TF atoms are estimated from sensor data, typically with greedy approaches like Matching Pursuit. Subsequently, the inverse problem is solved on the selected components using parametric [START_REF] Geva | Spatio-temporal matching pursuit (SToMP) for multiple source estimation of evoked potentials[END_REF], scanning [START_REF] Lelic | Inverse modeling on decomposed electroencephalographic data: A way forward[END_REF] or distributed methods [START_REF] Durka | Multichannel matching pursuit and EEG inverse solutions[END_REF][START_REF] Lina | Wavelet-based localization of oscillatory sources from magnetoencephalography data[END_REF]. Such methods suffer from several limitations. They implicitly assume that the source waveforms correspond to single TF atoms, while real brain signals are better modeled by a combination of atoms as proposed in this contribution. In addition, estimation errors made in the first step have a direct impact on the accuracy of the sources estimates. This is a particularly critical issue since the first step does not take into account the biophysics of the problem, i.e., the solution of the forward problem.

Certainly motivated by the ability of dipole fitting methods to explain M/EEG data, spatial sparsity of source configurations has also been a recurrent assumption to improve the performance of the M/EEG distributed source models. However, variational formulations based on mixed-norms [START_REF] Ou | A distributed spatiotemporal EEG/MEG inverse solver[END_REF][START_REF] Gramfort | Mixed-norm estimates for the M/EEG inverse problem using accelerated gradient methods[END_REF], or Bayesian methods with sparsity inducing mechanisms [START_REF] Friston | Multiple sparse priors for the M/EEG inverse problem[END_REF][START_REF] Wipf | A unified Bayesian framework for MEG/EEG source imaging[END_REF] that have been proposed so far in the literature do not model the transient oscillatory nature of brain signals. These approaches make the strong assumption that the source time courses are stationary. For example, results are strongly dependent on the time interval considered. Also the solutions obtained by these solvers are invariant with respect to the permutation of the columns of M, i.e., the temporal sequence of the data is immaterial. This contribution does not suffer from this limitation by explicitly modeling the temporal dynamics of the signals.

To improve over MNE type methods, by removing the assumption of a constant SNR over time and to correlate the time instants together, one should mention recently proposed state-space models based on Kalman filters and smoothers [START_REF] Long | State-space solutions to the dynamic magnetoencephalography inverse problem using high performance computing[END_REF][START_REF] Lamus | A spatiotemporal dynamic distributed solution to the meg inverse problem[END_REF]. Such methods, although promising, are still computationally demanding and could suffer from more practical issues like rather long time series to estimate the latent parameters reliably and the necessity to work with fixed orientations.

In [START_REF] Trujillo-Barreto | Bayesian M/EEG source reconstruction with spatio-temporal priors[END_REF], an inverse solver that models the transient and non-stationary responses in M/EEG is proposed. A probabilistic model with wavelet shrinkage is employed to promote spatially smooth time courses. The estimation however relies on model approximations with no guarantee on the solution obtained. The most related work to ours, beyond the field of M/EEG, is probably [START_REF] Model | Signal reconstruction in sensor arrays using sparse representations[END_REF] where sparsity is also promoted on the TF decompositions. The related optimization problem, is however solved with a truncated Newton method which can only be applied to differentiable problems. The nondifferentiability of the cost function is tackled by using smooth approximation in the minimization. Moreover, Newton methods are known to be fast in the neighborhood of the solution, but little is known about the global convergence rate. In [START_REF] Nesterov | Cubic regularization of newton's method and its global performance[END_REF], it is proved that a suitable Newton technique has the same rate of convergence as the accelerated first order schemes like the one we are employing in TF-MxNE. In this contribution, we do not address the problem of learning spatial basis functions as in [START_REF] Valdés-Sosa | EEG source imaging with spatio-temporal tomographic nonnegative independent component analysis[END_REF] as doing so would make the cost function non-convex, which would deteriorate the speed of convergence and would also make the solver dependent on the initialization. However, using a pre-defined dictionary of spatial basis functions in line with [START_REF] Bolstad | Space-time event sparse penalization for magneto-/electroencephalography[END_REF][START_REF] Haufe | Large-scale EEG/MEG source localization with spatial flexibility[END_REF] and multiplying the gain matrix with this dictionary, our prior could be used to estimate spatially extended sources with temporally smooth waveforms. This would, however, be significantly more computationally expensive.

In this work, we demonstrated how physiologically motivated priors for brain activations can be accounted for in a mathematically principled framework in M/EEG source analysis. Using a composite prior, the sparsity of spatial patterns, the temporal smoothness, and the non-stationarity of the source signals were well recovered. Thanks to the structure of the cost function considered, mainly its convexity, an efficient optimization strategy was proposed. The problem being convex, the solver is not affected by improper initialization and cannot be trapped in local minima. Simulations indicated benefits of the approach over alternative solvers, while results with well understood MEG data confirm the accuracy of the reconstruction with real signals. Both results show that our solver is a promising new approach for mining M/EEG data.

TF-MxNE

  the 1 norm, and A 21 = N i=1 N j=1 |A[i, j]| 2 the 21 mixed norm. A T and A H stand for the matrix transpose and a Hermitian transpose, respectively.

Figure 1 :

 1 Figure 1: a) Short-time Fourier transform (STFT) of a single channel MEG signal sampled at 1000 Hz showing the sparse nature of the transformation (window size 64 time points and time shift k 0 = 16 samples). b) STFT restricted to the 50 largest coefficients c) Data and data reconstructed using only the 50 largest coefficients.

Figure 2 :

 2 Figure 2: Sparsity patterns promoted by the different priors: 2 all non-zero, 1 scattered and unstructured non-zero, 21 block row structure, and 21 + 1 block row structure with intra-row sparsity. Red color indicates non-zero coefficients.

Algorithm 1 1 :

 11 FISTA with TF dictionaries to minimize 4 Input: Measurements M, gain matrix G, regularization parameter λ > 0 and I the number of iterations. Output: Z Auxiliary variables : Y and Z o ∈ R P ×K , and τ and τ o ∈ R.2: Estimate the Lipschitz constant L with the power iteration method.3: Y = Z = Z, τ = 1, 0 < µ < L -1 4: for i = 1 to I do 5:

Figure 3 :

 3 Figure3: Comparison of RMSE in the source space as a function of λ space (SNR=6bB and SNR=0bB with Gaussian white noise and AR noise). TF priors improve the reconstruction and the best accuracy is obtained with the TF 21 + 1 prior for both noise conditions. The bias correction improves the performance of the reconstruction in both cases, particularly for high values of λ space .

21 +

 21 1 non-zeros (NB)

Figure 4 :

 4 Figure 4: Simulations results with SN R = 6 dB. (a) A simulated source activations. (b) Noiseless simulated measurements. (c) Simulated measurements corrupted by noise. (d-ef) Estimation with 1 prior. (g-h-i) Estimation with 21 prior [34]. (j-k-l) Estimation with composite TF prior. (m-n-o) Estimation with composite TF prior and debiasing.(f-i-l-o) show the sparsity patterns obtained by the 3 different priors as explained in Fig.2. Result (j) shows how the composite TF prior improves over (d) and (g). (l) presents also a higher level of sparsity compared to (f) and (i).

28

 28

Figure 5 :

 5 Figure 5: Results with real EEG lead field (SNR=3dB). The 4 dipoles are color coded. Magenta dots show the 3D grid of sources. Dark dots show the EEG sensors locations. Contrary to TF 21 + 1 , 21 fails to recover the deep green dipole time course.

Figure 7 :

 7 Figure 7: Comparison of TF-MxNE, dSPM, and LCMV for an auditory stimulation. The images show the TF-MxNE solutions overlayed with dSPM estimates (a) and LCMV output map (b) 100 ms after stimulation. Since LCMV does not exploit the cortical orientation information, the map peaks are located at the gyri, below the auditory cortex in the temporal lobe.

  GX TF-MxNE (explained data) (d) Source locations in the right occipital cortex (visual cortex) (e) Source locations in right hemisphere with V1 (red) and V2 (yellow) labels.

Figure 8 :

 8 Figure 8: Results obtained with TF-MxNE for a visual stimulation (checkerboard stimulus in the left visual field) with unfiltered MEG data. Estimation was performed with a loose orientation (parameter 0.2), with a depth compensation of 0.9 on a set of 7498 cortical locations (G ∈ R 364×22494 ). Estimation with λ space = 30% of λ max space leads to 3 active brain locations in the contralateral visual cortex (V1 and V2). TF-MxNE leads to smooth time courses and zeros during baseline.

We can however say nothing about Φ H Φ in general.

http://ltfat.sourceforge.net/

prior on Z does not produce denoised time series as it does not promote source estimates that are formed by a sum of a few Gabor atoms. In order to recover the sparse row structure, while simultaneously promoting sparsity of the decompositions, we propose to use a composite prior formed by the

Acknowledgements. This work was supported by National Center for Research Resources P41 RR014075-11, the National Institute of Biomedical Imaging and Bioengineering grants 5R01EB009048 and R01 EB006385, the German Research Foundation (Ha 2899/8-2), and the French ANR ViMAGINEANR-08-BLAN-0250-02.

 8(a). For the early peak, there is a very good agreement between all three methods although dSPM and LCMV maps are smeared in space. For the later peak, LCMV and TF-MxNE agree on the location of the blue source.