Malek Kammoun

Jerome Piquereau

Lydie Nadal-Desbarats

Sandra Même

Maud Beuvin

Gisèle Bonne

Vladimir Veksler

Yann Le Fur

Philippe Pouletaut

William Même

Frederic Szeremeta

Jean-Marc Constans

Elizabeth S Bruinsma

Molly H Nelson Holte

Zeynab Najafova

Steven A Johnsen

Malayannan Subramaniam

John R Hawse
email: hawse.john@mayo.edu

Sabine F Bensamoun
email: sabine.bensamoun@utc.fr

Novel Role of Tieg1 in Muscle Metabolism and Mitochondrial Oxidative Capacities

Keywords: Tieg1, Klf10, Skeletal Muscle, Mitochondria, Metabolism

Aim: Tieg1 is involved in multiple signaling pathways, human diseases, and is highly expressed in muscle where its functions are poorly understood.

Methods:

We have utilized Tieg1 KO mice to identify novel and important roles for this transcription factor in regulating muscle ultrastructure, metabolism and mitochondrial functions in the soleus and EDL muscles. RNA sequencing, immunoblotting, TEM, MRI, NMR, histochemical and mitochondrial function assays were performed.

Results: Loss of Tieg1 expression resulted in altered sarcomere organization and a significant decrease in mitochondrial number. Histochemical analyses demonstrated an absence of SDH staining and a decrease in COX enzyme activity in KO soleus with similar, but diminished, effects in the EDL. Decreased complex I, COX and citrate synthase activities were detected in the soleus muscle of KO mice indicating altered mitochondrial function. Complex I activity was also diminished in KO EDL. Significant decreases in citrate synthase and respiratory chain complex activities were identified in KO soleus. 1 H-NMR spectra revealed no significant metabolic difference between WT and KO muscles. However, 31 P spectra revealed a significant decrease in phosphocreatine and ATP. Altered expression of 279 genes, many of which play roles in mitochondrial and muscle function, were identified in KO soleus muscle. Ultimately, all of these changes resulted in an exercise intolerance phenotype in Tieg1 KO mice.

Conclusion:

Our findings have implicated novel roles for Tieg1 in muscle including regulation of gene expression, metabolic activity and organization of tissue ultrastructure. This muscle phenotype resembles diseases associated with exercise intolerance and myopathies of unknown consequence.

INTRODUCTION

TGFβ inducible early gene-1 (Tieg1), also known as Krüppel-like factor 10 (Klf10), is a zinc finger-containing transcription factor known to be involved in multiple signaling pathways and disease processes [START_REF] Subramaniam | Functional role of KLF10 in multiple disease processes[END_REF][START_REF] Hawse | Estrogen-TGFbeta cross-talk in bone and other cell types: role of TIEG, Runx2, and other transcription factors[END_REF] . Tieg1 was originally discovered in human osteoblast cells as an early response gene following TGFβ treatment [START_REF] Subramaniam | Identification of a novel TGF-beta-regulated gene encoding a putative zinc finger protein in human osteoblasts[END_REF] . Clinically, altered expression levels of Tieg1, as well as single nucleotide polymorphisms (SNPs) in the Tieg1 gene, are associated with osteoporosis in humans [START_REF] Yerges | Candidate gene analysis of femoral neck trabecular and cortical volumetric bone mineral density in older men[END_REF][START_REF] Hopwood | Gene expression profile of the bone microenvironment in human fragility fracture bone[END_REF] . Further, SNPs in Tieg1 have been identified in patients with hypertrophic cardiomyopathy [START_REF] Bos | TGFbeta-inducible early gene-1 (TIEG1) mutations in hypertrophic cardiomyopathy[END_REF] . Finally, decreased expression of Tieg1 is observed during breast cancer progression [START_REF] Subramaniam | Tissue, cell type, and breast cancer stagespecific expression of a TGF-beta inducible early transcription factor gene[END_REF][START_REF] Reinholz | Differential gene expression of TGF beta inducible early gene (TIEG), Smad7, Smad2 and Bard1 in normal and malignant breast tissue[END_REF] and SNPs in the Tieg1 gene have been shown to predict risk of ovarian cancer [START_REF] Jim | Common Genetic Variation in Circadian Rhythm Genes and Risk of Epithelial Ovarian Cancer (EOC)[END_REF] .

Many of these clinical findings have been corroborated through the use of a TIEG1 knockout mouse model system. Specifically, global deletion of Tieg1 in mice results in an osteopenic skeletal phenotype [START_REF] Bensamoun | TGFbeta inducible early gene-1 knockout mice display defects in bone strength and microarchitecture[END_REF][START_REF] Hawse | TIEG-null mice display an osteopenic gender-specific phenotype[END_REF] , attenuates estrogen signaling in bone tissues [START_REF] Hawse | TGFbeta inducible early gene-1 plays an important role in mediating estrogen signaling in the skeleton[END_REF][START_REF] Hawse | Estrogen receptor beta isoformspecific induction of transforming growth factor beta-inducible early gene-1 in human osteoblast cells: an essential role for the activation function 1 domain[END_REF] and leads to hypertrophy of the heart [START_REF] Rajamannan | TGFbeta inducible early gene-1 (TIEG1) and cardiac hypertrophy: Discovery and characterization of a novel signaling pathway[END_REF] . Further characterization of this gene has revealed a central role for Tieg1 in regulating canonical Wnt signaling [START_REF] Subramaniam | TIEG1 modulates beta-catenin sub-cellular localization and enhances Wnt signaling in bone[END_REF][START_REF] Subramaniam | TIEG and estrogen modulate SOST expression in the murine skeleton[END_REF] , TGFβ signaling [START_REF] Johnsen | TGFbeta inducible early gene enhances TGFbeta/Smad-dependent transcriptional responses[END_REF][START_REF] Johnsen | Transcriptional regulation of Smad2 is required for enhancement of TGFβ/Smad signaling by TGFβ inducible early gene[END_REF][START_REF] Johnsen | Modulation of transforming growth factor beta (TGFbeta)/Smad transcriptional responses through targeted degradation of TGFbeta-inducible early gene-1 by human seven in absentia homologue[END_REF] as well as BMP signaling [START_REF] Hefferan | Cytokine-specific induction of the TGF-beta inducible early gene (TIEG): regulation by specific members of the TGF-beta family[END_REF][START_REF] Subramaniam | TIEG1 enhances Osterix expression and mediates its induction by TGFbeta and BMP2 in osteoblasts[END_REF][START_REF] Hawse | TIEG1/KLF10 modulates Runx2 expression and activity in osteoblasts[END_REF][START_REF] Yadav | Characterization of BMP signaling dependent osteogenesis using a BMP depletable avianized bone marrow stromal cell line (TVA-BMSC)[END_REF][START_REF] Prashar | Microarray meta-analysis identifies evolutionarily conserved BMP signaling targets in developing long bones[END_REF] in multiple tissues and cell types.

As part of the original description and characterization of the Tieg1 gene, it was shown to be most highly expressed in skeletal muscle [START_REF] Subramaniam | Identification of a novel TGF-beta-regulated gene encoding a putative zinc finger protein in human osteoblasts[END_REF] . However, little is known about its involvement in mediating the development, repair or function of this tissue. Recent studies performed by our group have shown that deletion of Tieg1 in mice results in different functional responses of slow and fast twitch muscle fibers [START_REF] Kammoun | Impact of TIEG1 Deletion on the Passive Mechanical Properties of Fast and Slow Twitch Skeletal Muscles in Female Mice[END_REF][START_REF] Kammoun | Development of a novel multiphysical approach for the characterization of mechanical properties of musculotendinous tissues[END_REF] . In parallel to these mechanical tests, imaging analyses (MRI, 9.4T) revealed changes in the texture profile between wild-type (WT) and Tieg1 knockout (KO) mouse skeletal muscles [START_REF] Kammoun | Impact of TIEG1 on the structural properties of fast and slow twitch skeletal muscle[END_REF] . Histochemical analyses demonstrated hyperplasia and hypertrophy of all fiber types in Tieg1 KO muscles (soleus and extensor digitorum longus: EDL) with a significant increase in the percentage of glycolytic fibers and a decrease in the percentage of oxidative fibers [START_REF] Kammoun | Impact of TIEG1 Deletion on the Passive Mechanical Properties of Fast and Slow Twitch Skeletal Muscles in Female Mice[END_REF] . At the cellular level, others have implicated Tieg1 in myogenic differentiation given its role in downstream regulation of myostatin and TGFβ signaling in myoblasts [START_REF] Miyake | Possible role of TIEG1 as a feedback regulator of myostatin and TGF-beta in myoblasts[END_REF][START_REF] Miyake | TIEG1 negatively controls the myoblast pool indispensable for fusion during myogenic differentiation of C2C12 cells[END_REF] .

Here, we expand upon our previous studies to comprehensively characterize the skeletal muscle defects observed in Tieg1 KO mice and to begin to elucidate the molecular mechanisms by which deletion of a single nuclear transcription factor results in a severe skeletal muscle phenotype that resembles multiple forms of mitochondrial-related myopathies in humans. Specifically, we have utilized RNA sequencing (RNAseq) to characterize the Tieg1-regulated transcriptome in skeletal muscle and have identified mis-regulation of multiple genes implicated in mitochondrial function and muscle biology. Given these findings, we have further characterized the role of Tieg1 in skeletal muscle function and have identified significant defects with regard to mitochondrial oxidative capacities and metabolism in Tieg1 KO skeletal muscle. The present experiments shed light on the hierarchy of molecular changes occurring in Tieg1 KO muscles and enhance our understanding of muscle metabolism and mitochondrial oxidative capacities which are essential towards the goal of developing novel approaches for the treatment of myopathies and muscle-related diseases.

MATERIALS AND METHODS

Mice and Study design

The generation of Tieg1 KO mice has been previously described [START_REF] Subramaniam | TIEG1 null mouse-derived osteoblasts are defective in mineralization and in support of osteoclast differentiation in vitro[END_REF] and initial characterization of Tieg1 expression and function in skeletal muscle were previously reported [START_REF] Subramaniam | Identification of a novel TGF-beta-regulated gene encoding a putative zinc finger protein in human osteoblasts[END_REF][START_REF] Kammoun | Impact of TIEG1 Deletion on the Passive Mechanical Properties of Fast and Slow Twitch Skeletal Muscles in Female Mice[END_REF][START_REF] Kammoun | Impact of TIEG1 on the structural properties of fast and slow twitch skeletal muscle[END_REF] . For the studies presented here, littermate female animals from heterozygous breeding pairs were utilized and sacrificed for analysis at approximately 3 months of age. All mice were maintained in a temperature controlled room (22 ± 2°C) with light/dark cycle of 12 hours.

RNA isolation and sequencing

For the RNA sequencing studies, a total of 6 WT and 6 Tieg1 KO mice were utilized.

Immediately following sacrifice, the soleus muscle was dissected from mice, placed in trizol and sheared using an Ultra-Turrax homogenizer prior to extraction of total RNA using a miRNeasy Mini Kit (Qiagen, Germantown, MD). Total RNA was quantitated with a nanodrop and equal amounts of RNA from 2 individual WT or Tieg1 KO mice were pooled to generate a total of 3 replicates per genotype. All samples were tested for RNA integrity before submitting for sequencing. mRNA sequencing was performed by the Mayo Clinic Molecular Biology Sequencing Core Facility using standard procedures. RNA libraries were prepared using 200 ng of good quality total RNA according to the manufacturer's instructions for the TruSeq RNA Sample Prep Kit v2 (Illumina, San Diego, CA) employing poly-A mRNA enrichment using oligo dT magnetic beads. The final adapter-modified cDNA fragments were enriched by 12 cycles of PCR using Illumina TruSeq PCR primers. The concentration and size distribution of the completed libraries was determined using a Fragment Analyzer (AATI, Ankeny, IA) and Qubit fluorometry (Invitrogen, Carlsbad, CA). Libraries were sequenced at 40-50 million fragment reads per sample following Illumina's standard protocol using the Illumina cBot and HiSeq 3000/4000 PE Cluster Kit. The flow cells were sequenced as 100 X 2 paired end reads on an Illumina HiSeq 4000 using HiSeq 3000/4000 sequencing kit and HCS v3.3.52 collection software. Base-calling was performed using Illumina's RTA version 2.7.3.

RT-PCR

RNA was harvested from an additional 5 WT soleus and 8 EDL muscles as described above.

One µg of RNA was reverse transcribed using the iScript™ cDNA Synthesis Kit (Bio-Rad, Hercules, CA) and real-time PCR was performed in triplicate for Tieg1 and β-Actin using a Bio-Rad iCycler and a PerfeCTa™ SYBR Green Fast Mix™ for iQ real-time PCR kit (Quanta Biosciences, Gaithersburg, MD) as specified by the manufacturer. Cycling conditions were as follows: 95°C for 30 seconds followed by 40 cycles of 95°C for 3 seconds and 60°C for 30 seconds. Melt curves were generated to ensure amplification of a single PCR product.

Quantitation of the PCR results was calculated based on the threshold cycle (Ct) following normalization to β-Actin and averaged among muscle type. PCR primers were designed using Primer3 software (http://frodo.wi.mit.edu/primer3/) and were purchased from Integrated DNA Technologies (Coralville, IA). Primer sequences were as follows: Tieg1 F: GGGACAAAGCTGAGCAGAGT; Tieg1 R: GGAAGCAAGCTGTCATCCTC; β-Actin F: AATCGTGCGTGACATCAAAGAG; β-Actin R: GCCATCTCCTGCTCGAAGTC.

Histochemical staining analysis

Soleus and EDL muscles were harvested from Tieg1 KO (N=10) and WT (N=10) mice and immediately frozen in an isopentane bath that was cooled by liquid nitrogen and subsequently stored at -80°C until the time of analysis. Ten m serial sections were cut perpendicular to the muscle using a Microm HM 500M cryotome at -20°C. Analysis of succinate dehydrogenase (SDH), cytochrome c oxidase (COX) and menadione [START_REF] Wattenberg | Effects of coenzyme Q10 and menadione on succinic dehydrogenase activity as measured by tetrazolium salt reduction[END_REF] were performed with histochemical staining to identify potential alterations in mitochondrial metabolism. Menadione is a quinone that is active during mitochondrial electron transport and has been reported to serve as an electron acceptor for several metalloproteins. The study of Wattenberg et al. [START_REF] Wattenberg | Effects of coenzyme Q10 and menadione on succinic dehydrogenase activity as measured by tetrazolium salt reduction[END_REF] has demonstrated that menadione enhances the activity of SDH. Additionally, the endoplasmic reticulum (Dinitrophenylhydrazine), connective tissue (Gömöri trichrome), glycogen content (Periodic acid-Schiff) and lipid content (Oil Red O) were determined. Each staining procedure was performed using standard histochemical methods 32 on two independent tissue sections which were subsequently visualized under a Leica microscope using a 20 x objective.

High-resolution images were captured with an Olympus cooled digital camera (DP-72) and interpretation of staining was based on the color histogram of the image using ImageJ software (US National Institutes of Health, Bethesda, MD, United States) [START_REF] Schneider | NIH Image to ImageJ: 25 years of image analysis[END_REF] . This analysis separately measures the intensity in the red, green and blue channels of the image.

Transmission Electron Microscopy (TEM)

Prior to processing for TEM, the soleus (N=7) and the EDL (N=7) were dissected from WT and Tieg1 KO mice with both extremities pinned to avoid muscle contraction and to maintain original length. Muscles were then immediately placed in fixative [1% (vol/vol) glutaraldehyde and 4% (vol/vol) paraformaldehyde in 0.1 M phosphate buffer, pH 7.2] [START_REF] Mcdowell | Histologic fixatives suitable for diagnostic light and electron microscopy[END_REF] and incubated overnight at 4°C. Subsequently, two pieces (1 mm x 2 mm) of tissue were dissected from the middle of the tissue and rinsed for 30 min in three changes of 0.1 M phosphate buffer, pH 7.2, followed by a 1h secondary fixation in phosphate-buffered 1% OsO4 and 30 min in 1% uranyl acetate at room temperature. Following dehydration in a series of ethanol washes, tissue was embedded in EMbed 812 resin (EMS, Hatfield PA) and polymerized at 60°C for 18h [START_REF] Spurr | A low-viscosity epoxy resin embedding medium for electron microscopy[END_REF] . Ultrathin (90 nm) sections were cut using a Leica UC7 ultramicrotome (Buffalo Grove, IL), placed on 200 mesh copper grids, and stained with lead citrate. Five micrographs of each specimen were randomly captured across the muscle tissue using a JEOL 1400Plus TEM (Peabody, MA), operating at 80 kV with a magnification of 42000x. ImageJ 1.46/Java 8 software (National Institute of Health, Bethesda, MD, United States) [START_REF] Schneider | NIH Image to ImageJ: 25 years of image analysis[END_REF] was used to manually quantify the mitochondria area and the thick myofilament diameter.

Enzyme assays

Complete enzyme extractions from pieces of frozen soleus and EDL skeletal muscle tissue (N=7 WT, N=7 KO) were performed in an ice-cold buffer (50 mg ml -1 ; containing (in mM): HEPES 5 (pH 8.7), EGTA 1, dithiothreitol 1, and 0.1% Triton X-100) using a Bertin Precellys 24 homogenizer (Bertin, Montigny-le-Bretonneux, France). Protein concentrations were assessed using the bicinchoninic acid assay. Total enzyme activities of citrate synthase (CS), mitochondrial NADH: coenzyme Q oxidoreductase (Complex I), mitochondrial succinic dehydrogenase (SDH, Complex II) and cytochrome oxidase (COX, Complex IV) were determined (30°C, pH 7.5) using standard spectrophotometric protocols as previously described [START_REF] Sousa | Subcellular creatine kinase alterations. Implications in heart failure[END_REF][START_REF] Estornell | Assay conditions for the mitochondrial NADH:coenzyme Q oxidoreductase[END_REF] . Activities were reported as μmol min -1 / g protein (μmol min -1 , international units, IU). Mann Whitney tests were performed to compare enzyme activities as a function of genotype. A P-value < 0.05 was considered significant.

Mitochondrial respiration

Respiratory parameters of total mitochondrial populations were studied in situ in saponinpermeabilized fibers isolated from soleus and EDL muscles (N = 7 WT, N = 7 Tieg1 KO). We followed the protocol described previously [START_REF] Veksler | Mitochondrial respiratory parameters in cardiac tissue: a novel method of assessment by using saponin-skinned fibers[END_REF][START_REF] Kuznetsov | Analysis of mitochondrial function in situ in permeabilized muscle fibers, tissues and cells[END_REF] using a Clark electrode (Strathkelvin Instruments, Glasgow, Scotland) in a water-jacketed oxygraphic cell filled with 3 ml of respiration solution containing Complex I substrates (10 mM glutamate + 4 mM malate) at 22°C with continuous stirring. Respiration rates were determined in the absence of ADP (basal respiration, Vo) and in the presence of a saturating (2 mM) ADP concentration.

Functional activities of different complexes within the mitochondrial respiratory chain were estimated in the presence of succinate + amytal (or Amobarbital) to assess Complex IIsupported respiration, and tetramethyl-p-phenylenediamine (TMPD) + ascorbate to assess the Complex IV-supported respiration. Respiration rates were expressed as μmol O2 min -1 g -1 fiber dry weight. Mann Whitney tests were performed to compare mitochondrial respiration as a function of genotype. A P-value < 0.05 was considered significant.

Immunoblotting

Frozen soleus and EDL tissue samples from WT and Tieg1 KO mice (N=6) were homogenized (Bertin Precellys 24) in ice cold buffer containing (in mM): HEPES 50 (pH=7.1), KCl, ethylenediaminetetraacetic acid (EDTA), β-glycerophosphate, orthovanadate 1, dithithreitol 1, sodium fluoride 50, Na pyrophosphate 5, phenylmethylsulfonyl fluoride 0.2 and an anti-protease cocktail (Calbiochem 539134) in the presence of Triton X-100 0.1%.

Equal amounts of protein extracts were separated on SDS-polyacrylamide gels (10 to 12%) and transferred to polyvinylidene difluoride membranes for Western blotting. After 1h of blocking in PBS containing TWEEN20 (0.1%) and nonfat milk (5%), the membranes were incubated overnight at 4°C with primary antibodies (mitochondrial OXPHOS (complex I subunit NDUFB8, complex II subunit CII 30, complex III sunbunit core 2 and complex V subunit alpha): MS604 (mitosciences), Citrate synthase: ab96600 (Abcam), COX IV: ab14744 (Abcam), MFN2: H00009927-MO1 (Abnova), PGC-1: sc13067 (Santa Cruz), OPA1: 612606 (Becton Dickinson), DRP1: 611112 (BD Biosciences). Oxyblot™protein oxidation detection kit (Millipore S7150) was used for protein carbonylation assessment.

After washing, the membranes were incubated with secondary antibodies coupled with horseradish peroxidase for 1h at room temperature and visualized using chemiluminescent substrate (Luminata TM Western Chemiluminescent HRP Substrates, Millipore). Light emission was detected by autoradiography and quantified using an image-analysis system (Bio-Rad). Mann Whitney tests were performed to compare protein abundance as a function of genotype. A P-value < 0.05 was considered significant.

In vivo Magnetic Resonance Spectroscopy

Twenty four mice were used (N=12 WT and N=12 Tieg1 KO) for magnetic resonance (MR) acquisitions that were performed on a 9.4 Tesla Biospec 94/21 superconducting magnet (Bruker, Wissembourg, France) with a shielded gradient set (950 mT.m -1 maximum gradient amplitude). A home built homogeneous scroll double tuned coil (10 mm inner diameter) was developed. The resonance frequencies were 400 MHz for the proton and 162 MHz for the phosphorus. Animals were anesthetized by inhalation of 2.5% isoflurane then maintained during MR experiments at 1.5% (0.5 L.min -1 mixed in air and oxygen with 1:1 ratio). The physiological body temperature was maintained inside the magnet by circulating warm water.

A pressure sensor was used to monitor the respiration cycle and obtain the respiration frequency.

Localized 1 H spectroscopy was performed with a PRESS (point resolved spectroscopy) sequence (TE/TR, 16 ms/4000 ms, 256 accumulations, 17min) with VAPOR module (variable pulse power and optimized relaxation delays) for water suppression. A 22 mm 3 voxel was centered on the middle of the hindlimb.

31

P spectroscopy was then performed with a PRESS sequence (TR = 1500ms, 1 pulse BP of 100µs 200 accumulations, DS = 2, TA = 5min). Due to the lack of sensitivity of 31 P spectroscopy the entire hindlimb was used for the acquisition (1000mm 3 ). Spectra were analyzed with a homemade (CRMBM, Aix Marseille university) CSIAPO software used for quest interface [START_REF] Ratiney | Semi-parametric time-domain quantification of HR-MAS data from prostate tissue[END_REF] that works in time domain. n a 1 H spectrum the assessment of choline, creatine, taurine and extramyocellular lipids were made. On a 31 P spectrum, different metabolites are quantifiable: Inorganic phosphate, phosphocreatine, P, ATP and ATPQuantification of the area under peaks was performed in arbitrary units (AU). A Mann Whitney statistical test was used to compare those areas between both muscle genotypes.

H-NMR Metabolomics

The protocol for extraction of metabolites from soleus (N=18) and EDL (N=18) muscles was adapted from the Folch type two-step procedure described by Wu et al. [START_REF] Wu | High-throughput tissue extraction protocol for NMR-and MS-based metabolomics[END_REF] for sample mass < 50 mg. Prior to NMR analysis, dried samples were reconstituted in 220µL of deuterated buffer containing TSP [3-trimethylsilylpropionic acid] at 145µM and transferred to conventional 3mm NMR tubes (CortecNet, Paris, France). [START_REF] Subramaniam | Functional role of KLF10 in multiple disease processes[END_REF] H-NMR spectra were obtained with a Bruker DRX-600 AVANCE-III HD spectrometer (Bruker SADIS, Wissembourg, France), operating at 14 T, with a TCI cryoprobe. Standard water suppressed 1 H-NMR spectra were acquired at 298K using a "noesypr1d" pulse sequence with relaxation delay of 20s. Spectra were processed using Topspin version 3.2 software (Bruker Daltonik, Karlsruhe, Germany). [START_REF] Subramaniam | Functional role of KLF10 in multiple disease processes[END_REF] H-NMR spectra were automatically reduced to ASCII files using AMIX Software package (Analysis of MIXture, version 3.9.14, Bruker Biospin, Karlsruhe, Germany) as previously described [START_REF] Beauclercq | Serum and Muscle Metabolomics for the Prediction of Ultimate pH, a Key Factor for Chicken-Meat Quality[END_REF] . Spectral intensities were scaled to the total spectral intensity and the resulting data (supplementary data S1 A-B) was analyzed by multivariate and univariate statistical analyses.

Multivariate and Univariate analysis

Principal Component Analyses (PCA) and Partial Least Square Discriminant Analyses (PLS-DA) were performed using SIMCA-P + Software (version 13.0, Umetrics, Umeå, Sweden) on the NMR dataset containing soleus and EDL as a function of genotype. All data sets were scaled to unit of variance allowing all metabolites to become equally important. PCA describes the total variability within the dataset and can be used as an informative indicator while the PLS-DA is used to predict the spectral features (metabolites) that define separation between groups (phenotypes) [START_REF] Worley | PCA as a practical indicator of OPLS-DA model reliability[END_REF] . The overall quality of the models was judged by cumulative R 2 (goodness of fit) and cumulative Q 2 (goodness of prediction). Variable Importance in the Projection (VIP) > 1 were considered as the most contributive in the genotype separation.

Student's t-test were performed using MetaboAnalyst [START_REF] Xia | MetaboAnalyst 3.0--making metabolomics more meaningful[END_REF] for all VIP measures. A P-value < 0.05 was considered significant.

Following statistical analysis, heatmaps were generated to visually depict changes in metabolites (VIP > 1 and a P-value < 0.05) as a function of mouse genotype (Supplementary data S1 A-B).

Bioinformatic analysis

The FASTQ files were mapped to the mouse genome (assembly mm9) using TopHat Gappedread mapper [START_REF] Kim | TopHat2: accurate alignment of transcriptomes in the presence of insertions, deletions and gene fusions[END_REF] with very sensitive Bowtie2 settings on Galaxy Platform (version 2.1.0) [START_REF] Goecks | Galaxy: a comprehensive approach for supporting accessible, reproducible, and transparent computational research in the life sciences[END_REF] .

The read counting was performed via HTSeq 47 (version 0.6.1galaxy3) with the following parameters: -f bam -r name -s no -a 10 -t exon -m union. The count files were subsequently subjected for differential analysis using the DESeq2 Package [START_REF] Love | Moderated estimation of fold change and dispersion for RNA-seq data with DESeq2[END_REF] on R (Bioconductor version 3.4.1). The threshold values to determine the differentially expressed genes were as follows: abs (log2fold change) > 0.8 and padj < 0.05. The pathway analysis on the regulated genes was performed using "Enrichr" [START_REF] Chen | Enrichr: interactive and collaborative HTML5 gene list enrichment analysis tool[END_REF][START_REF] Kuleshov | Enrichr: a comprehensive gene set enrichment analysis web server 2016 update[END_REF] .

Treadmill exercise

The physical function of WT and Tieg1 KO mice was assessed using a motorized treadmill (Columbus Instruments, Columbus, OH) test of exercise ability. Briefly, 15 WT and 15 Tieg1 KO female mice are acclimated to the treadmill for 2 consecutive days immediately followed by a test day. Acclimation consisted of 5 minutes of running at a speed of 5 m per minute for the first two minutes, 7 m per minute for the next 2 minutes and 9 m per minute for the last minute at an incline grade of 5%. For the test, mice ran on the treadmill to exhaustion with a 5% incline and starting speed of 5 m per minute for two minutes. The speed was increased by 2 m per minute every two minutes until the mice were exhausted. Exhaustion was defined as the inability of a mouse to stay on the treadmill in spite of a mild electrical stimulus and mechanical prodding. Total distance run was recorded for each mouse and group averages were calculated. Mann Whitney test were performed to compare the total distance run between genotypes. A P-value < 0.05 was considered significant.

RESULTS

Gene expression profile changes in Tieg1 KO muscles

Given the morphological and functional changes observed in the skeletal muscles of Tieg1 KO mice [START_REF] Kammoun | Impact of TIEG1 Deletion on the Passive Mechanical Properties of Fast and Slow Twitch Skeletal Muscles in Female Mice[END_REF][START_REF] Kammoun | Impact of TIEG1 on the structural properties of fast and slow twitch skeletal muscle[END_REF] , we sought to elucidate Tieg1-dependent molecular mechanisms controlling skeletal muscle maintenance and function. We first assessed the expression levels of Tieg1 in soleus and EDL muscle isolated from WT female mice. Tieg1 expression was shown to be approximately 4-fold higher in soleus muscle compared to EDL muscle (Figure 1A). In light of these results, we next assessed the transcriptome-wide effect of Tieg1 loss in the soleus of 12 week old mice (Figure 1B -1C). In total, 279 genes were differentially regulated between the WT and Tieg1 KO muscle (Figure 1D). Although the number of downregulated genes were substantially fewer, the pathway analyses indicated that a number of these genes are known targets of the Myod1 transcription factor (Supplementary data, S2) which was shown to be upregulated in the muscle of Tieg1 KO mice (Log2FC=1.46, padj=3.06E-16).

Interestingly, some of the downregulated genes (Pln, Tfcp2l1, Ide, Bcl2l1) are also known components of mitochondria (Supplementary data, S2) whereas most of the upregulated genes were associated with collagen and extracellular matrix organization (Supplementary data, S3).

Mitochondrial metabolism deficiencies in

Disorganization of Tieg1 KO muscle ultrastructure

Given the observed differences in staining patterns for important mitochondrial enzymes, we next sought to analyze the ultrastructure of the soleus and EDL muscle fibers using TEM.

Longitudinal sections revealed disorganization of the muscle structure, and were suggestive of a contracted phenotype, in both Tieg1 KO soleus and EDL muscles (Figure 4) compared to WT littermates. Specifically, smaller sarcomere lengths (1.5m for KO and 2.5m for WT), disappearance of the I band, shortening of the H band and changes in the shape of mitochondria where observed in Tieg1 KO soleus (Figure 4D) and EDL muscles (Figure 4J).

The quantification of mitochondrial area (A) revealed a significant (P < 0.001) decrease for soleus and EDL Tieg1 KO muscles (AWT_soleus = 0.079 ± 0.003 m 2 vs ATieg1 KO_soleus = 0.054 ± 0.002 m 2 and AWT_EDL = 0.08 ± 0.003m 2 vs ATieg1 KO_EDL = 0.045 ± 0.002 m 2 ).

Mitochondrial area was classified within different ranges based on Picard's study [START_REF] Picard | Trans-mitochondrial coordination of cristae at regulated membrane junctions[END_REF] and a non-Gaussian curve was observed in both muscles with smaller mitochondria (Asoleus < 0.04 µm 2 : 37% vs 14% and AEDL < 0.04 µm 2 : 58% vs 23%), and a complete absence of mitochondria with an area greater than 0.15 m 2 , in Tieg1 KO mice.

Changes in mitochondrial activity in Tieg1 KO muscle

A significant decrease in citrate synthase (CS) and respiratory chain complex activities were observed in Tieg1 KO soleus muscle (Figure 5A -5C). This result is in agreement with the decreased area occupied by mitochondria and the weaker COX and SDH staining (Figure 2) observed in Tieg1 KO soleus fibers. No significant changes were observed for COX and CS activities in Tieg1 KO EDL muscle (Figure 5D and5F). However, NADH: ubiquinone oxidoreductase (Complex I) (Figure 5B and5E) activity was decreased in both Tieg1 KO soleus and EDL muscles.

Respiration rates of total mitochondrial populations using substrates for complexes I, II and IV were markedly lower in permeabilized Tieg1 KO soleus muscle (Figure 6A) compared with WT controls. In contrast, mitochondria in Tieg1 KO EDL muscle did not exhibit any differences in respiration (Figure 6B) as compared with WT. These results are in agreement with the enzyme activities showing a significant decrease only for the Tieg1 KO soleus.

In soleus, protein content (Figure 7) of all complexes measured (I, II, III, V) was slightly lower in Tieg1 KO muscle though these changes did not reach statistical significance except for complex IV (COX). No significant differences were observed for the protein of mitochondrial fusion (OPA1), mitofusin2 (Mfn2) as well as mitochondrial fission (DRP1) and biogenesis (PGC1-) in the soleus muscle. For the EDL muscle, no significant differences were observed for any of these proteins between WT and Tieg1 KO mice. Total protein carbonylation was determined to assess oxidative stress in skeletal muscles. No significant difference in the levels of oxidative stress were detected between WT and Tieg1 KO soleus or EDL muscles (Figure 8).

Lack of Tieg1 induces metabolic changes

1 H spectra showed no significant difference for choline, creatine, taurine and extramyocellular lipids between WT and Tieg1 KO muscle (Figure 9B). However, [START_REF] Wattenberg | Effects of coenzyme Q10 and menadione on succinic dehydrogenase activity as measured by tetrazolium salt reduction[END_REF] P spectra revealed a significant difference in the levels of phosphocreatine and ATP(Figure 10B). The PCA model obtained with the 1 H-NMR data set containing WT soleus, Tieg1 KO soleus, WT EDL and Tieg1 KO EDL revealed significant separation between soleus and EDL muscles (PC1) explaining 50% of the total variability within the dataset (Figure 11). Mouse genotype (PC2) explained 24% of the total variability within the data (Figure 11). These results highlight differences between WT and Tieg1 KO metabolic profiles. PLS-DA models obtained with 1 H-NMR datasets from soleus muscles (Supplementary data, S1 A, predictive ability Q²=0.72, goodness of fit R²Y(cum)=0.99) and EDL muscles (Supplementary data S1 B, predictive ability Q²=0.84, goodness of fit R²Y(cum)=0.99) revealed a clear separation between WT and KO genotypes regardless of muscle type suggesting that deletion of Tieg1 results in significant changes in the metabolic profiles of skeletal muscle. The PLS-DA models allow for the identification of very important metabolites in projection (VIPs) whose changes are likely to be biologically relevant (Supplementary data, S4). Figure 12 depicts heatmaps (12A for soleus muscle, 12B for EDL muscle) generated for metabolites with a VIP > 1 in the PLS-DA models. For both Tieg1 KO soleus and EDL muscles, there were more down regulated metabolites compared to WT muscles. Of note, there were differences in the metabolites that were identified to be differentially regulated by loss of Tieg1 expression between the two muscle types.

Loss of Tieg1 expression results in exercise intolerance

In light of the many defects observed in Tieg1 KO skeletal muscle including gene expression changes, muscle ultrastructure and mitochondrial activity and metabolism, we sought to determine if these defects impacted muscle function in vivo. Treadmill exercise tests were performed to determine the distance run prior to exhaustion. These tests revealed significant decreases in the distance run by Tieg1 KO mice compared to WT littermates (Figure 13). On average, WT mice ran 309 m prior to exhaustion while Tieg1 KO mice ran only 214 m equating to an approximately 25% decrease in exercise ability (Figure 13).

DISCUSSION

We have previously shown that deletion of the Tieg1 gene in mice results in muscle fiber hypertrophy, glycolytic fiber type switching, texture profile changes and defects in the functional properties of skeletal muscle [START_REF] Kammoun | Impact of TIEG1 Deletion on the Passive Mechanical Properties of Fast and Slow Twitch Skeletal Muscles in Female Mice[END_REF] . Here, we have expanded upon these original discoveries to further investigate the molecular and cellular basis for the observed skeletal muscle defects. Specifically, we have shown that Tieg1 plays a central role in regulating mitochondrial oxidative capacities and metabolism.

Through gene expression profiling, we observed up-regulation of the Myod1 "myogenic differentiation 1" gene in Tieg1 KO muscles. This gene is well known to be a regulator of myogenic differentiation [START_REF] Fuchtbauer | MyoD and myogenin are coexpressed in regenerating skeletal muscle of the mouse. Developmental dynamics : an official publication of the American Association of Anatomists[END_REF][START_REF] Kammoun | A network-based approach for predicting Hsp27 knock-out targets in mouse skeletal muscles[END_REF][START_REF] Blum | The role of MyoD1 and histone modifications in the activation of muscle enhancers[END_REF] and its mis-regulation may contribute to the observed muscle fiber hypertrophy in Tieg1 KO mice. However, this assumption must be confirmed by additional experiments to identify the exact connection between Tieg1 and Myod1. Moreover, a disorganization of the muscle ultrastructure and a contracted phenotype were observed in Tieg1 KO skeletal muscle as defined by smaller sarcomeres, the absence of I bands and the shortening of H bands. Unger's group showed that chaperone proteins were capable of aggregating and binding to titin filaments resulting in shorter sarcomeres in hereditary myopathies (Duchenne, Titinopathy, etc …) [START_REF] Unger | Translocation of molecular chaperones to the titin springs is common in skeletal myopathy patients and affects sarcomere function[END_REF] . However, RNAseq analysis did not identify significant differences in Hsp expression levels and it will be of interest to quantify the abundance of HSPs in Tieg1 KO muscle in future studies.

Pathway analyses of genes identified to be differentially expressed in Tieg1 KO skeletal muscle revealed alterations in the expression levels of 48 genes that are known to be related to muscle diseases (such as Duchenne muscular dystrophy, Becker muscular dystrophy, etc …).

Interestingly, the genes implicated in these muscle diseases were primarily collagens (Col1a1, Col1a2, Col5a1, Col6a2, etc …) which were all upregulated in Tieg1 KO mice. A study by Bönnemann et al. has previously shown that collagen VI-related myopathies encompass a spectrum of diseases ranging from severe Ullrich muscular dystrophy to mild Bethlem myopathy [START_REF] Bönnemann | The collagen VI-related myopathies: muscle meets its matrix[END_REF] . Moreover, it has been demonstrated that Col1a1 is the most commonly evaluated collagen gene for fibrosis [START_REF] Seet | Upregulation of distinct collagen transcripts in post-surgery scar tissue: a study of conjunctival fibrosis[END_REF] , and our results are also in agreement with DiMario's study showing an increase of Col1a1 gene expression, and protein content, in Tieg1 KO muscle (tibialis anterior and diaphragm) [START_REF] Dimario | KLF10 Gene Expression Modulates Fibrosis in Dystrophic Skeletal Muscle[END_REF] . In addition, this phenomenon is more accentuated in the Tieg1 KO/mdx (Duchenne muscular dystrophy) mouse model [START_REF] Dimario | KLF10 Gene Expression Modulates Fibrosis in Dystrophic Skeletal Muscle[END_REF] .

Tieg1 KO mice also exhibit increased expression of TGFβ within the skeletal muscle. We have specifically demonstrated that Tieg1 represses the expression of Smad 7 [START_REF] Johnsen | TGFbeta inducible early gene enhances TGFbeta/Smad-dependent transcriptional responses[END_REF] and induces the expression of Smad 2 [START_REF] Johnsen | Transcriptional regulation of Smad2 is required for enhancement of TGFβ/Smad signaling by TGFβ inducible early gene[END_REF] . We have also shown that over-expression of Tieg1 mimics the effects of TGF treatment in multiple cell types [START_REF] Hefferan | Overexpression of a nuclear protein, TIEG, mimics transforming growth factor-beta action in human osteoblast cells[END_REF] . Overall, we have shown that Tieg1 serves as a positive feedback loop for TGF signaling. Based on these previous observations, it is possible that increased expression of TGF in the soleus of Tieg1 KO mice is a compensatory mechanism for diminished TGF pathway activity in the absence of Tieg1 expression. Similar results were observed after muscle rotator cuff tears in a rat model which results in a significant increase of fibrosis and muscular dystrophy. To explain this phenomenon, Liu et al [START_REF] Liu | Upregulation of transforming growth factor-β signaling in a rat model of rotator cuff tears[END_REF] have linked TGF to the Akt/mTOR pathway and suggest that TGFβ may be the primary factor regulating muscle changes after muscle rotator cuff tears. In general, TGFβ pathways are highly conserved pathways that exert a potent level of control over muscle gene expression and are critical regulators of fibrosis in multiple organ systems.

In addition to muscular dystrophy, mitochondrial myopathy pathways were also identified from the gene expression data. For instance, Rmrp (RNA component of mitochondrial RNAase P) was identified in these pathways and was significantly down-regulated in Tieg1 KO soleus muscle. The decreased expression of mitochondrial complex proteins have confirmed the involvement of Tieg1 with regard to mitochondrial oxidative capacities and may explain some of the structural and functional changes of the mitochondria observed in Tieg1 KO muscle. These mitochondrial disease are also related to another disease called "mitochondrial complex I deficiency". In the present study, complex I activity and respiration were found to be decreased in Tieg1 KO muscle, further demonstrating an important role for Tieg1 in mitochondrial activity and function. Such a specific impairment of complex I function has been observed in several diseases including cardiomyopathy and heart failure [START_REF] Marin-Garcia | Regional distribution of mitochondrial dysfunction and apoptotic remodeling in pacing-induced heart failure[END_REF][START_REF] Scheubel | Dysfunction of mitochondrial respiratory chain complex I in human failing myocardium is not due to disturbed mitochondrial gene expression[END_REF] . Indeed, complex I specific and functional activities could be decreased, when other ETC (electron transport chain) complexes or total mitochondrial mass (usually estimated by citrate synthase activity) are normal. This is the case, for example, in generalized myopathy in advanced heart failure [START_REF] Barreiro | Therapeutic Approaches in Mitochondrial Dysfunction, Proteolysis, and Structural Alterations of Diaphragm and Gastrocnemius in Rats With Chronic Heart Failure[END_REF] . Similarly, complex I in skeletal muscle mitochondria is the primary target of the immunosuppressive molecule cyclosporine toxicity [START_REF] Sanchez | Effect of cyclosporin A and its vehicle on cardiac and skeletal muscle mitochondria: relationship to efficacy of the respiratory chain[END_REF] . A marked decrease in complex I activity was found in skeletal muscle in acute intermittent porphyria (a disease primarily affecting the nervous system), an autosomal dominant metabolic disease [START_REF] Homedan | Mitochondrial energetic defects in muscle and brain of a Hmbs-/-mouse model of acute intermittent porphyria[END_REF] . Such a specific inhibition of complex I has an important impact on mitochondrial function in striated muscle since this complex is the rate limiting step of the electron flux through the ETC.

Similarly, deficiencies in SDH (Succinate dehydrogenase) have been reported in human mitochondrial diseases [START_REF] Kollberg | Transient restoration of succinate dehydrogenase activity after rhabdomyolysis in iron-sulphur cluster deficiency myopathy[END_REF][START_REF] Haller | Deficiency of skeletal muscle succinate dehydrogenase and aconitase. Pathophysiology of exercise in a novel human muscle oxidative defect[END_REF][START_REF] Vladutiu | Succinate dehydrogenase deficiency[END_REF] [START_REF] Vladutiu | Succinate dehydrogenase deficiency[END_REF] exhibiting muscle pain, exercise intolerance and cardiomyopathy, etc. As is the case with Tieg1 KO mice as reported here, a decrease in all respiratory chain complexes was associated with a deficiency of SDH activity (23% of patients) showing an important reduction for 12 patients. The results reported in these clinical studies have described skeletal muscle disorders resulting from mitochondrial defects closely resembling the muscle phenotype of Tieg1 KO mice.

Oxygen consumption rate and the enzyme activity analyses reflect clear decreases in Tieg1 KO muscle, whereas the protein content of respiratory chain complexes do not (though all complexes do exhibit a slight but non-significant decrease). Indeed, the respiration rates do not depend only on absolute mitochondrial protein content, but rather are a function of numerous enzyme activities, lipid environment, inner membrane properties etc. Even specific enzyme activities do not exactly reflect the enzyme protein content. Moreover, mitochondrial complexes represent structures formed by various subunits. Generally, the antibodies against one subunit do not necessarily reveal protein content of the whole complex. It is for these reasons that investigation of mitochondrial function requires the complete analysis of oxygen consumption rates and specific enzyme activities in addition to protein content determination in order to fully understand the basis for a given mitochondrial defect. The same rationale concerns the SDH staining. This method determines SDH activity rather than SDH protein content. Moreover, the staining reveals the enzyme activity in situ rather than specific activity of purified SDH. The results obtained show that Tieg1 deficiency primarily alters the functional properties of mitochondrial proteins rather than their expression levels. Interestingly, we have observed that loss of Tieg1 expression has a strong muscle-type specific effect in that, the soleus is more affected than that of the EDL. The differences in the severity of this muscle phenotype between the soleus and EDL may be due to the observed higher expression of Tieg1 in WT soleus muscle compared to WT EDL muscle. It will be of interest to further characterize the glycolytic switch of oxidative muscle. These observations are in agreement with DiMario's study where loss of Tieg1 expression resulted in a pronounced fibrotic phenotype in slow twitch (diaphragm) muscle compared to fast twitch muscle (tibialis anterior) [START_REF] Dimario | KLF10 Gene Expression Modulates Fibrosis in Dystrophic Skeletal Muscle[END_REF] . However, through 1 H-NMR metabolic profiling, we have shown that both slow and fast twitch muscles are affected by deletion of Tieg1. Metabolites such as fumarate, lactate and pyruvate were up-regulated in Tieg1 KO EDL muscle. These metabolites are involved in the Krebs cycle and reflect alterations in mitochondria metabolism. Nicotinurate was up-regulated, and beta-hydroxy-butyrate was down-regulated, in both muscles of Tieg1 KO mice. Nicotinurate is used to diagnose disorders associated with mitochondrial fatty acid beta-oxidation and beta-hydroxy-butyrate, a ketone body, which is produced by beta-oxydation as an energy source. Furthermore, in vivo [START_REF] Wattenberg | Effects of coenzyme Q10 and menadione on succinic dehydrogenase activity as measured by tetrazolium salt reduction[END_REF] The present study has demonstrated a new role for Tieg1 in the regulation of skeletal muscle energy metabolism, specifically in oxidative muscle which strongly depends on mitochondrial function. The observed differences in mitochondrial activity could in part be explained by the known fiber type switching that occurs in Tieg1 KO mice [START_REF] Kammoun | Impact of TIEG1 Deletion on the Passive Mechanical Properties of Fast and Slow Twitch Skeletal Muscles in Female Mice[END_REF] . Thus, the Tieg1 KO mice mimic human diseases associated with muscle fatigue and exercise intolerance and could allow for the study of mitochondrial myopathies associated with deficiency of SDH in skeletal muscle. Indeed, we have demonstrated that loss of TIEG expression results in an exercise intolerance phenotype in mice that is likely explained by the many defects in mitochondrial function described here. For these reasons, Tieg1 KO mice represent a novel model for further understanding muscle pathologies and mitochondrial dysfunction. These findings have laid the groundwork for future efforts aimed at identifying specific genes and their associated pathways that can be pharmacologically manipulated to reverse/alleviate these defects in Tieg1 KO mice and to evaluate their efficacy for the treatment of muscle disorders in humans. Densitometry was used to quantitate protein levels and is depicted in the bar graphs for soleus (A-B) and EDL (C-D). Data are represented in arbitrary unit (A.U.) as the mean ± SEM and significant differences between genotypes were determined using Mann Whitney test ( † : P < 0.01).

  Animals had free access to water and were fed with standard laboratory chow ad libitum. This study was carried out in strict accordance with the recommendations in the Guide for the Care and Use of Laboratory Animals of the National Institutes of Health. The protocol was approved by the French ministry of higher education, research and innovation (Permit Number: DUO-4776) and the Mayo Clinic Institutional Animal Care and Use Committee (Permit Number: A9615).

  Tieg1 KO muscles Given the importance of mitochondrial function in skeletal muscle, we next histologically assessed the activity of some complexes of the mitochondrial electron transfer chain within muscle biopsies of WT and Tieg1 KO mice. Remarkably, as shown in Figure 2A -2B, SDH staining was completely absent in soleus muscle of Tieg1 KO mice while normal staining patterns were observed in WT littermates. In addition, COX (Figure 2C -2D) and menadione (Figure 2E -2F) staining was less intense in Tieg1 KO mice compared to WT controls. No significant differences in the staining for endoplasmic reticulum, connective tissue, glycogen content or lipids were observed (Figure 2G -2L). Decreased staining for SDH, COX and menadione were also observed in the EDL muscle of Tieg1 KO mice, albeit to a lesser degree (Figure 3).

  P-NMR spectroscopy revealed decreased ATP (adenosine triphosphate) and PCr levels that are involved in energy metabolism, again implicating an important role for Tieg1 in mitochondrial function. The amino-acids valine, leucine, lysine and glutamate were also shown to be down-regulated in both soleus and EDL muscles and are associated with protein metabolism. This metabolic imbalance in both Tieg1 KO muscle types can be related to mitochondrial pathways confirming the strong implication of the Tieg1 gene in the function and activity of mitochondria.

FIGURE 1

 1 FIGURE 1 Differential gene expression analysis in the skeletal muscle of 3-month old Tieg1

FIGURE 2

 2 FIGURE 2Histochemical analysis of Tieg1 KO (N=10) and WT (N=10) soleus muscles.

FIGURE 3

 3 FIGURE 3Histochemical analysis of Tieg1 KO (N=10) and WT (N=10) EDL muscles.

FIGURE 4

 4 FIGURE 4 Transmission electron microscopy analyses of soleus and EDL muscles isolated

FIGURE 5

 5 FIGURE 5 Mitochondrial enzyme activities in skeletal muscle. Citrate synthase activity (A

FIGURE 6

 6 FIGURE 6Mitochondrial respiratory activities in skeletal muscle. Mitochondrial complex

FIGURE 7

 7 FIGURE 7 Mitochondrial protein expression in soleus and EDL muscle. Relative expression

FIGURE 8

 8 FIGURE 8 Total protein carbonylation in skeletal muscle. Protein carbonylation was

FIGURE 9 1 H

 1 FIGURE 9 1 H spectral analysis of skeletal muscle. (A) In vivo hindlimb muscle 1 H spectra

FIGURE 10 31 P

 31 FIGURE 10 31 P spectral analysis of skeletal muscle. (A) In vivo hindlimb muscle 31 P spectra

FIGURE 11

 11 FIGURE11 Principal Component Analysis (PCA) of the 1 H-NMR spectral data. WT soleus

FIGURE 12

 12 FIGURE 12 Metabolic analysis of soleus and EDL muscles. Heatmaps are represented to

FIGURE 13

 13 FIGURE13 Treadmill exercise test of endurance. Three month old WT (N=15) and Tieg1

  . Haller et al. have reported on a 22 year old man with life-long exercise intolerance marked by premature fatigue 67 . Additionally, Kollberg et al. described a 77 year old woman with muscle fatigue and dyspnea resulting from SDH deficiency 66 . Moreover, Vladutiu et al. have assessed qualitative and quantitative SDH expression in muscle biopsies from 108 patients (age range: 1 to 68 years) with mitochondrial myopathies

ACKNOWLEDGMENTS

This work was carried out and funded in the framework of the Labex MS2T. It was supported by the French Government, through the program "Investments for the future" managed by the National Agency for Research (Reference ANR-11-IDEX-0004-02).

This work was also supported by the National Institutes of Health (R01 DE14036 to JRH and MS) and the Eisenberg Foundation (JRH and MS). We would also like to thank Dr. Thomas Spelsberg for his excellent mentorship.

We thank the Mayo Microscopy and Cell Analysis Core for experimental and technical support.

CONFLICT OF INTEREST

The authors declare no competing interests.


El-Houcine Bergou
email: elhoucine.bergou@inrae.fr

Youssef Diouane
email: youssef.diouane@isae-supaero.fr

Vyacheslav Kungurtsev
email: vyacheslav.kungurtsev@fel.cvut.cz

Complexity iteration analysis for strongly convex multi-objective optimization using a Newton path-following procedure

Keywords: Multi-objective optimization, strongly convex optimization, path-following, Newton method, complexity analysis Mathematics Subject Classification (2010) 49M05, 49M15, 90C06, 90C60

In this note, we consider the iteration complexity of solving strongly convex multi-objective optimization problems. We discuss the precise meaning of this problem, noting that its definition is ambiguous, and focus on the most natural notion of finding a set of Pareto optimal points across a grid

Introduction

Consider the following multi-objective optimization problem:

min x∈R n f (x) := {f i (x)} i=1,...,m (1)
where f : R n → R m is a strongly convex and twice continuously differentiable function. Our target is to find weak Pareto-optimality points for problem [START_REF] Allgower | Introduction to numerical continuation methods[END_REF].

We recall that weak Pareto-optimality holds at a point x ∈ R n if for all d ∈ R n , there exists an i ∈ {1, . . . , m} such that

∇f i ( x) d ≥ 0.
For single objective optimization, worst-case iteration complexity quantifies the number of iterations that could be necessary, in the worst-case (i.e., for the most ill-behaved problems), for an algorithm to achieve a certain level of satisfaction of an approximate measure of optimality, typically a small norm for the gradient [START_REF] Nesterov | Introductory lectures on convex optimization: A basic course[END_REF]. Classically, the multi-objective optimization community had not considered attempting to derive bounds on iteration complexity for problems in vector optimization. Recently, however, works have appeared, see [START_REF] Fliege | Complexity of gradient descent for multiobjective optimization[END_REF][START_REF] Calderòn | On high-order model regularization for multiobjective optimization[END_REF],

which consider the iteration complexity of gradient descent for multi-objective optimization. In both papers rates were derived for obtaining some point satisfying approximate weak Pareto-optimality.

To the best of our knowledge, only [START_REF] Fliege | Complexity of gradient descent for multiobjective optimization[END_REF] considers iteration complexity for specifically strongly convex objectives. However, in deriving their complexity result, convergence of the algorithm to some Pareto optimal point is assumed.

Moreover, we believe that a much stronger and more meaningful result can be shown by considering the precise meaning of the problem. In particular, consider the so-called scalarized problem, parametrized by {λ i } i=1,..,m

min x∈R n m i=1 λ i f i (x), (2)
for any {λ i } i=1,..,m ∈ D = {{λ i } i=1,..,m ∈ R m / 0 ≤ λ i ≤ 1, m i=1 λ i = 1}, the unit simplex of R m . A stationary point of the problem ( 2) is also Pareto optimal for [START_REF] Allgower | Introduction to numerical continuation methods[END_REF]. Thus, one can find a Pareto optimal point, at least for strongly convex multi-objective problems, by simply choosing any arbitrary convex combination {λ i } i=1,..,m ∈ D and solving the resulting mono-objective problem, thus the worst-case iteration complexity of finding some Pareto optimal point is already a known problem, it corresponds to the worst-case iteration complexity of solving a single objective strongly convex optimization problem.

In the multi-objective optimization literature, e.g., [START_REF] Marler | Survey of multi-objective optimization methods for engineering[END_REF], scalarization is typically, at most, one step in the process of finding the solution of a multi-objective problem, where the ultimate definition of a solution can vary. In particular, it can be that the goal of the optimization is (a) tracing the Pareto front itself, so in some sense finding all, or some adequate approximation to all, stationary points, or (b) finding an appropriately best point along the Pareto front through some secondary metrics, or (c) using an interactive environment with a human participant who grades potential solutions.

In this note, we shall concern ourselves with the first task: establish complexity bounds for some appropriate notion of finding the entire Pareto front.

To this end, we define the problem as, for all λ ∈ Λ, finding min

x∈R n m i=1 λ i f i (x) (3)
where Λ ⊂ D is some finite grid of elements. In the definition of D, given the constraint on the sum, we can consider m -1 dimensions as free which in turn entirely determine the remaining λ i . We thus divide each side of the hypercube [0, 1] by some desired width of the grid d, and thus there are

1 d m-1
total possible grid points, where a denotes the greatest integer less than or equal to a. Conversely, if there is some desired grid pre-defined by the user, we can define the quantity d denoting the maximum width between two neighbors on the grid, that is

d = min (λ,λ )∈Λ , λ =λ max 1≤i≤m |λ i -λ i | .
We would like to emphasize that the procedure of obtaining a Pareto front by solving a set of scalarized single-objective reformulations is suited strictly to strongly convex objectives, in the non-convex case this strict correspondence between the two may be lost.

We organize this note as follows. In Section 2, we describe our algorithm to solve the strongly convex multi-objective optimization problem. We explain how to use a Newton path-following procedure to find the entire Pareto front.

Section 3 addresses the convergence of our algorithm by characterizing its iteration complexity. A numerical illustration on a simple example demonstrating the efficiency of our approach is given in Section 4. Conclusions are given in Section 5.

2 Path-following for finding the entire Pareto front.

Recall that for a fixed λ 0 ∈ Λ, using only first order information one can solve a strongly convex optimization problem of the type (3) at best linearly (with a gradient descent based method, see for instance [2, Theorem 3.18]).

Namely, in order to obtain a point distance from the solution of (3) for a fixed In this note, we propose finding the entire Pareto front by performing pathfollowing, a method based on the implicit function theorem. Later we will show that the proposed strategy will reduce the overall iteration complexity drastically relative to naively solving every scalarized problem separately. To start with, for some initial grid point λ (0) ∈ Λ we obtain the solution x (0) ∈ R n of the problem (3), (for instance, by using a gradient descent method with the following stopping criterion m j=1 λ (0)

λ 0 ∈ Λ, O(log( 1
j ∇f i (x) ≤ ). Note that such point x (0)
gives the Pareto optimal point of the problem (1) associated with λ (0) . Now, let λ (1) be one of the closest neighbors to λ (0) in the finite grid Λ, our goal is to apply a predictor-corrector scheme to compute a new x (1) corresponding to an approximate solution to the scalarized problem associated with λ (1) .

Path-following, or tracing a set of solutions for a parametrized nonlinear system of equations across a range of parameters, is an important algorithmic tool, for which an introduction can be found in [START_REF] Allgower | Introduction to numerical continuation methods[END_REF]. Closest to our work, a predictor-corrector path-following procedure for strongly convex optimization problems (interpreted as strongly regular variational inequalities) is given in [START_REF] Dontchev | An Euler-Newton continuation method for tracking solution trajectories of parametric variational inequalities[END_REF]. In this work it is shown that, for this parametric problem, a property of uniform strong regularity holds and a procedure involving one tangential predictor (Euler) and one corrector (Newton) step result in a sequence of iterates whose distance to a set of solutions to the parametric variational inequality is of the order of d 4 , where recall that d is, in this context, the grid spacing.

Thus there exists C such that if d ≤ C( ) 1/4 , a set of solutions with approximate optimality across a set of parameters can be found. If applied to the multi-objective Pareto front context, the number of Euler-Newton continuation steps is the number of grid points, which corresponds to

d -1 ≥ C -1 -1/4 .
If the desired grid is already small enough, then it is clear that this pathfollowing procedure outperforms the naive method of solving the standalone Algorithm 1: A Newton path-following for finding the Pareto front.

Input: Let Λ = {λ (0) , . . . , λ (p-1) } ⊂ R m + be some finite grid of cardinality p satisfying: for all j = 0, . . . , p -1, m i λ

(j) i = 1, λ (j+1) is one of the closest neighbors of λ (j) not yet visited.
Output: The entire Pareto front by performing path-following associated with Λ:

x (0) , x (1) , . . . , x (p-1) .

Compute an initial Pareto optimal point x (0) , i.e.,

x (0) = arg min x f (0) (x), where f (0) (x) = m i=1 λ (0) i f i (x). ( 4
)
Set k = 0.

Step 1: Compute a predictor x(k+1) , i.e.,

x(k+1) = x (k) -   m j=1 λ (k) j ∇ 2 f j x (k)   -1 m i=1 λ (k+1) i -λ (k) i ∇f i x (k) (5)
Step 2: Apply a Newton correction to compute x (k+1) , i.e., starting from x(k+1) run the Newton method to find

x (k+1) = arg min x f (k+1) (x), where f (k+1) (x) = m i=1 λ (k+1) i f i (x). (6)
If k = p -1 then Stop, otherwise increment k by 1 and go to Step 1.

problem at every grid point. Otherwise, it depends on the magnitude of the desired number of additional grid points required to perform path-following.

We consider an alternative predictor-corrector scheme that is more aggressive in its use of potentially longer tangential steps and multiple Newton iterations. In particular, this is more suitable for obtaining the set of solutions across the Pareto front with the tightest iteration complexity bound. This predictor-corrector procedure will be repeated until we traverse the entire set Λ. A formal description of the algorithm is given as Algorithm 1.

The predictor step (Step 1 of Algorithm 1) is motivated by the implicit function theorem. Therefore, first, let's recall the implicit function theorem adapted to our context.

Theorem 1 Let g : R n+m → R n be a continuously differentiable function for a parametrized system of equations,

g(x, λ) = 0, where x ∈ R n and λ ∈ R m .
Consider that there exists a solution satisfying g(x 0 , λ 0 ) = 0. If the Jacobian matrix J g,x (x 0 , λ 0 ) of g with respect to x is invertible, then there exists an open neighborhood B ⊂ R m such that there exists a unique continuously differentiable path x(λ) defined on λ ∈ B with x(λ 0 ) = x 0 and g( x(λ), λ) = 0 for all λ ∈ B. Furthermore, it holds that the derivative of x(λ) over B is given by

∂ x ∂λ (λ) = -[J g,x ( x(λ), λ)] -1 ∂g ∂λ ( x(λ), λ). (7)
We consider applying Theorem 1 to the optimality conditions of (2) given by the following parametrized system of equations

g(x, λ) = m i=1 λ i ∇f i (x) = 0.
Precisely, for a given iteration index k, consider that we have a solution x (k) at some λ (k) ∈ Λ, i.e., m j=1 λ (k)

j ∇f i (x (k) ) = 0. ( 8
)
Since we consider strongly convex objectives, it holds that the matrix

m i=1 λ (k) i ∇ 2 f i (x)
is invertible for all x ∈ R n , in particular the inverse norm is bounded by the inverse of the weighted sum of the strong convexity constants of {f i }. Thus by Theorem 1 we have that there exists some ball B (k) around λ (k) and a unique path x(λ) such that x(λ

(k) ) = x (k) and m i=1 λ i ∇f i ( x(λ)) = 0, for all λ ∈ B (k) .
Furthermore, the derivative of the path given by ( 7) is defined for all λ ∈ B (k) to satisfy,

∂ x ∂λ (λ) = -   m j=1 λ j ∇ 2 f j ( x(λ))   -1 [∇f 1 ( x(λ)) , . . . , ∇f m ( x(λ))] .
Consider now a Taylor expansion of x(λ) along λ ∈ B (k) from the base point k) . This is given by

x(λ (k) ) = x (
x(λ) = x (k) -   m j=1 λ (k) j ∇ 2 f j (x (k)   -1 m i=1 λ i -λ (k) i ∇f i x (k) + O λ -λ (k) 2 .
Motivated by the discussion on Newton's method applied to path-following in [4, Chapter 5], we define a predictor x(k) (λ) by computing the first order Taylor approximation,

x(k) (λ) = x (k) -   m j=1 λ (k) j ∇ 2 f j (x (k)   -1 m i=1 λ i -λ (k) i ∇f i x (k) (9)
which is precisely the "tangent continuation method" with the order p = 2 as given in [START_REF] Deuflhard | Newton methods for nonlinear problems: affine invariance and adaptive algorithms[END_REF]Page 239]. Assuming that λ (k+1) is close enough to λ (k) (i.e., λ (k+1) ∈ B (k) ), the predictor step x(k+1) given by ( 5) in Algorithm 1 is defined as

x(k+1) = x(k) λ (k+1) .
Let η (k) be the norm of the residual x(λ (k+1) ) -x(k+1) , i.e.,

η (k) = x(λ (k+1) ) -x(k+1) .
There is a remaining algorithmic necessity before this becomes practical as the predictor x(k+1) does not necessarily satisfy the desired level of stationarity. In fact, to achieve a point closer to the actual solution, we consider a "corrector" step x (k+1) using a sequence of Newton iterations. Under a set of conditions, the ordinary Newton method is quadratically convergent towards the solution starting from the predicted point if this point is sufficiently close to the solution. Thus, we require that the predictor x(k+1) is sufficiently accurate and determine the size of the step λ (k+1) -λ (k) appropriately. Note that by definition of x we have

x (k+1) = x(λ (k+1) ).
3 Characterizing the Complexity of Algorithm 1

Based on the ideas above, we can consider iteration complexity in a new sense.

For a given iteration k, consider having an approximate solution to (3) for a particular λ (k) , up to a desired optimality tolerance . Then consider pathfollowing from λ (k) to some λ (k+1) where λ (k+1) -λ (k) is small enough (in terms of desired grid-spacing d) to be able to determine the associated solution on the Pareto front. The same procedure is repeated across all the grid Λ until all solutions of the Pareto front have been found.

Before developing our complexity analysis, we formally state our working assumptions on the objective function f .

Assumption 31

The objective function f : R n → R m is twice continuously differentiable and strongly convex. In particular, there exist two positive con-stants c > 0 and L > 0, such for all i ∈ {1, . . . , m}, x ∈ R n and y ∈ R n ,

c y 2 ≤ y ∇ 2 f i (x)y ≤ L y 2 . ( 10
)
In other words, for all i ∈ {1, . . . , m}, the eigenvalues of the Hessian of f i are uniformly bounded from below by c, and above by L, everywhere.

This implies the following condition regarding scaling invariance properties appropriate for Newton methods [START_REF] Deuflhard | Newton methods for nonlinear problems: affine invariance and adaptive algorithms[END_REF].

Lemma 1 Consider Assumption 31. For all λ ∈ Λ and x, y ∈ R n , one has

  m j=1 λ j ∇ 2 f j (x)   -1   m j=1 λ j ∇f j (y) - m j=1 λ j ∇f j (x)   ≤ L c x -y . ( 11
)
In this case, the mapping x → p i=1 λ i ∇ 2 f i (x) is said to be affine covariant Lipschitz.

Proof Let λ ∈ Λ, then from Assumption 31, we conclude that the mapping

x → m i=1 λ i f i (x) is L-smooth (i.e, its gradient is Lipschitz with constant L)
and that it is c-strongly convex. In fact, for all x and y ∈ R n , we have that this mapping satisfies

  m j=1 λ j ∇f j (y) - m j=1 λ j ∇f j (x)   ≤ L x -y ,
and that the inverse of its Hessian for all x is bounded by 1 c , i.e.,

  m j=1 λ j ∇ 2 f j (x)   -1 ≤ 1 c .
By combining these two inequalities with the Cauchy Schwarz inequality, we get (11).

In the next Lemma, we will show that the remainder η (k) is at most of the same order as the distance between λ (k+1) and λ (k) . This will be instrumental in giving sufficient conditions on the grid spacing necessary to ensure quadratic local convergence of the ordinary Newton method when it is applied to solve (6) starting from x(k+1) (see Lemma 3).

Lemma 2 Consider Assumption 31. Then there exists a constant η > 0 such that for all k, one has

η (k) ≤ η λ (k+1) -λ (k) .
Proof In fact, one has

η (k) = x (k+1) -x(k+1) ≤ x (k+1) -x (k) + ∂ x ∂λ (λ (k) ) (λ (k+1) -λ (k) ) ≤ x(λ (k+1) ) -x(λ (k) ) + max λ∈Λ ∂ x ∂λ (λ) λ (k+1) -λ (k) ) .
On the other hand, by Theorem 1, the function λ → x(λ) is continuously differentiable for all λ ∈ D, thus, it is bounded and Lipschitz continuous over the compact set D, i.e., there exists c 1 > 0 and c 2 > 0 such that

x(λ (k+1) ) -x(λ (k) ) ≤ c 1 λ (k+1) -λ (k) ) and max λ∈Λ ∂ x ∂λ (λ) ≤ c 2 .
Hence,

η (k) ≤ (c 1 + c 2 ) λ (k+1) -λ (k) ) ,
which completes the proof.

The next result will show that to ultimately get an -Pareto optimal solution, applying the correction step will require a number of iterations of the ordinary Newton method of order log log 1 . We will start by recalling some sufficient conditions concerning the local convergence of the Newton method adapted to our setting. These results are from [START_REF] Deuflhard | Newton methods for nonlinear problems: affine invariance and adaptive algorithms[END_REF]. In [4, Theorem 5.2], the authors gave sufficient conditions on the functions {f j (•)} and the distance between λ (k+1) and λ (k) to ensure the convergence of the ordinary Newton method applied to solve (6) starting from x(k+1) , and in [4, Theorem 2.3] the authors showed the classic local quadratic convergence property of the Newton method when it is used to solve [START_REF] Fliege | Complexity of gradient descent for multiobjective optimization[END_REF]. In fact under the strong convexity assumption and as long as the distance between λ (k+1) and λ (k) is sufficiently small, the ordinary Newton method converges quadratically to the minimizer of ( 6) from the starting point x(k+1) .

Lemma 3 Let Assumption 31 hold. For a given iteration index k, consider

λ (k) ∈ Λ and x (k) such that m j=1 λ (k) j ∇f i x (k) ≤ . Let λ (k+1) ∈ Λ such that λ (k+1) -λ (k) ≤ c ηL , ( 12
)
where η is as in Lemma 2. Then the ordinary Newton method with the starting point x(k+1) (as given by ( 5)) converges and the computational cost of achieving a solution point x (k+1) such that

m j=1 λ (k+1) j ∇f i x (k+1) ≤ , is of order log log 1 .
Proof Let x (k+1) j be the j th iterate produced by an ordinary Newton method starting from x(k+1) . By Lemma 1, the mapping x →

p i=1 λ (k+1) i ∇ 2 f i (x)
is affine covariant Lipschitz. Hence, using [4, Theorem 5.2] and with x(k+1) being the predictor step defined as [START_REF] Dontchev | An Euler-Newton continuation method for tracking solution trajectories of parametric variational inequalities[END_REF], one can deduce that the sequence x (k+1) j j∈N generated by the ordinary Newton method starting from x(k+1) converges towards the solution x (k+1) (i.e., x (k+1) = lim j→∞ x (k+1) j ). Using Lemma 2 and (12), one gets

x (k+1) 0 -x (k+1) = x(k+1) -x (k+1) = η (k) ≤ η λ (k+1) -λ (k) ≤ c L . ( 13
)
In this case, using [4, Theorem 2.3], the Newton method converges quadratically, i.e.,

x (k+1) j+1 -x (k+1) ≤ L 2c x (k+1) j -x (k+1) 2 .
Hence,

x (k+1) j -x (k+1) ≤ L 2c 2 j-1 x (k+1) 0 -x (k+1) 2 j .
Thus, using (13), one deduces that

x (k+1) j -x (k+1) ≤ c L 2 -2 j-1 .
Thus,

m j=1 λ (k) j ∇f i x (k+1) j = m j=1 λ (k) j ∇f i x (k+1) j - m j=1 λ (k) j ∇f i x (k+1) ≤ L x (k+1) j -x (k+1) ≤ c2 -2 j-1 .
This implies that the computational cost of achieving the desired level of stationarity is of order log log 1 .

Thus the iteration complexity of Algorithm 1 is just of the order of complexity for solving a standalone strongly convex problem (i.e., computing x (0) ) added with 1 d m-1 multiplied by the cost of a predictor and the iterated Newton step. We formalize this with the following theorem, Theorem 2 Let Assumption 31 hold. Define N to be the number of iterations required to obtain x (0) , a point that has distance at most from the optimal point corresponding to (2) at λ (0) . Assume that the maximum width between any two neighbors on the grid Λ is,

d ≤ min c ηL , d , (14)
with d the minimal desired distance between lattice points.

Then, the overall iteration complexity of Algorithm 1 is

N + O 1 d m-1 log log 1 .
Proof First, note that, for each iteration k of Algorithm 1, the complexity of the predictor step is constant as its computational cost does not depend on .

For the corrector Newton step, since one has

λ (k+1) -λ (k) ≤ d ≤ c ηL ,
Lemma 3 implies that the complexity of running the Newton iterations until approximate optimality is O log log 1 . The proof is thus completed since the total number of lattice points in the grid Λ is at most

1 d m-1 .
Note that by using a gradient solver, the first term N is of order log(1/ ).

Hence, one can see that the complexity is generally favorable compared to the naive method of solving the strongly convex problem at every grid point separately, as log log 1 log(1/ ) for small .

Remark 1 Note that both the naive method of solving every problem across the grid points and path-following are both about equally parallelizable with perfect speedup as long as the number of grid points is larger than the number of processors. We can split the grid into disjoint components, and each processor finds one point in its part of the convex hull of allowable {λ i } and proceeds to path-follow across section of the grid assigned to it.

Numerical Illustration

To show the numerical performance of our approach compared to the naive method (which corresponds to the Gradient Descent method applied sequentially to the set of problems (2) defined for varying λ), we consider a simple problem given in [START_REF] Huband | A review of multiobjective test problems and a scalable test problem toolkit[END_REF], defined by

f (x) = (x 1 -1) 2 + (x 1 -x 2 ) 2 , (x 2 -3) 2 + (x 1 -x 2 ) 2 .
Since we have two objective functions, the vector λ has two components λ 1 and λ 2 where λ 1 + λ 2 = 1. In our experiment, we discretize λ in a uniform grid with a grid step-size d (the desired distance between the lattice points).

In our Matlab illustration, we will call Multi-GD the naive method and GD+Newton Pathfollowing the implementation of our Algorithm 1 (where we used the standard Gradient Descent method to find the first Pareto optimal point and then apply the Newton path-following procedure). For the Gradient Descent method, we used a random initial point x 0 and a stepsize equal to 1/λ max where λ max is the maximum eigenvalue for the Hessians of f 1 and f 2 .

We stopped the methods when the norm of the gradient was less than 10 -7 .

The obtained results are shown in Figure 1. One can see that both methods are able to find a similar Pareto Front (independently of the value of the grid spacing d). In term of the elapsed CPU time to find the front, our proposed algorithm can be seen to be faster than the naive method. In particular, one can see that for some values of d, the method GD+Newton Pathfollowing can obtain the approximately optimal solution with up to 10 times faster total run-time than the Multi-GD method. We conducted other experiments

(not reported here) on many toy problems and in all of them our method was outperforming the naive method in run-time while finding essentially the same front.

Conclusion

In this note, we studied the iteration complexity of a class of strongly convex multi-objective optimization problems. We observed that the notion of iteration complexity is not uniquely defined as there can be varying possible criteria of what it means to solve a multi-objective optimization problem. By working with the most context-independent criterion (namely, finding the set of all Pareto optimal points on a front), we demonstrated that finding the solution of one scalarized problem and then path-following across the grid to obtain the others is superior to finding the solution of every problem independently.

  / )) iterations must be taken, with each iteration involving the computation of one gradient vector. As a result, naively, one can obtain the entire Pareto front by solving each of the 1 d m-1 scalarized problems defined across the grid points independently with a variant of gradient descent, to obtain an overall complexity of O log(1/ ) 1 d m-1 .

  Pareto Front with d = 10 -1 .

  Pareto Front with d = 10 -2 .

  CPU time in seconds.

Fig. 1

 1 Fig. 1 Pareto Front and CPU time comparison, using Multi-GD and GD+ Newton Path-following, for different values of d.

Acknowledgements

We would like to thank two anonymous referees for their careful readings and corrections that helped us to improve our manuscript significantly. E. Bergou received support from the AgreenSkills+ fellowship programme which has received funding from the EU's Seventh Framework Programme under grant agreement No FP7-609398 (AgreenSkills+ contract). V. Kungurtsev received support from the OP VVV project CZ.02.1.01/0.0/0.0/16 019/0000765 "Research Center for Informatics".