Gaëlle Recher

Pascal Coumailleau

Denis Rouède

François Tiaho
email: francois.tiaho@univ-rennes1.fr

Structural origin of the drastic modification of second harmonic generation intensity pattern occurring in tail muscles of climax stages xenopus tadpoles

Keywords: skeletal muscle, second harmonic generation microscopy, two photon excitation fluorescence, electron microscopy, xenopus, apoptosis

Second harmonic generation (SHG) microscopy is a powerful tool for studying submicron architecture of muscles tissues. Using this technique, we show that the canonical single frequency sarcomeric SHG intensity pattern (SHG-IP) of premetamorphic xenopus tadpole tail muscles is converted to double frequency (2f) sarcomeric SHG-IP in metamorphic climax stages due to massive physiological muscle proteolysis. This conversion was found to rise from 7% in premetamorphic muscles to about 97% in fragmented muscular apoptotic bodies. Moreover a 66% conversion was also found in non-fragmented metamorphic tail muscles. Also, a strong correlation between predominant 2f sarcomeric SHG-IPs and myofibrillar misalignment is established with electron microscopy. Experimental and theoretical results demonstrate the higher sensitivity and the supra resolution power of SHG microscopy over TPEF to reveal 3D myofibrillar misalignment. From this study, we suggest that 2f sarcomeric SHG-IP could be used as signature of triad defect and disruption of excitation-contraction coupling. As the mechanism of muscle proteolysis is similar to that found in mdx mouse muscles, we further suggest that xenopus tadpole tail resorption at climax stages could be used as an alternative or complementary model of Duchene muscular dystrophy.

Introduction

Excitation-contraction coupling at sarcomeric level required triads composed of a central transverse tubule (T-tubule), flanked by two closely apposed sarcoplasmic reticulum (SR) terminal cisternae. Triads and myofibrils are regularly ordered in adult skeletal muscles and this organization is developmentally regulated [START_REF] Ezerman | Differentiation of the sarcoplasmic reticulum and T system in developing chick skeletal muscle in vitro[END_REF][START_REF] Veratti | Investigations on the fine structure of striated muscle fiber read before the Reale Istituto Lombardo, 13 March 1902[END_REF]. Structural intimate relation between the two has been first observed at Z-bands of amphibian larvae tail muscles [START_REF] Porter | Studies on the endoplasmic reticulum. III. Its form and distribution in striated muscle cells[END_REF], suggesting that the functional efficiency of triads could be coupled to perfect structural myofibrillar alignment at their A-, M-, I-and Zbands. This perfect myofibrillar arrangement is altered during muscle development, damage and regeneration [START_REF] Takekura | Eccentric exercise-induced morphological changes in the membrane systems involved in excitation-contraction coupling in rat skeletal muscle[END_REF][START_REF] Yeung | Development of T-tubular vacuoles in eccentrically damaged mouse muscle fibres[END_REF]. Since myofibrillar misalignment could damage triads and disrupt excitation-contraction coupling, a non-invasive label-free method that can probe in vivo or in vitro the 3D muscular micrometric architecture in thick muscle tissues is therefore expected to provide paramount structural information concerning local proteolysis and alteration occurring at triads.

Second harmonic generation (SHG) microscopy is a low-invasive technique that relies on a nonlinear optical interaction with hyperpolarizable endogeneous molecules like collagen and myosin causing scattered coherent radiation at twice the fundamental frequency [START_REF] Campagnola | Three-dimensional high-resolution second-harmonic generation imaging of endogenous structural proteins in biological tissues[END_REF]. Hence, SHG highlights the myosin filaments of sarcomeric Abands. It has proved to be an extremely beneficial contrast mechanism for label-free imaging of skeletal muscle in vivo, in disease state as well as during dynamic contraction [START_REF] Buttgereit | From chaos to split-ups--SHG microscopy reveals a specific remodelling mechanism in ageing dystrophic muscle[END_REF][START_REF] Llewellyn | Minimally invasive highspeed imaging of sarcomere contractile dynamics in mice and humans[END_REF][START_REF] Nucciotti | Probing myosin structural conformation in vivo by second-harmonic generation microscopy[END_REF][START_REF] Plotnikov | Measurement of muscle disease by quantitative second-harmonic generation imaging[END_REF][START_REF] Prent | Intermyofilament dynamics of myocytes revealed by second harmonic generation microscopy[END_REF][START_REF] Ralston | Detection and imaging of non-contractile inclusions and sarcomeric anomalies in skeletal muscle by second harmonic generation combined with two-photon excited fluorescence[END_REF]Rouede et al., 2013b;[START_REF] Rouede | Myofibrillar misalignment correlated to triad disappearance of mdx mouse gastrocnemius muscle probed by SHG microscopy[END_REF]. Our research group has previously shown that the predominant feature of sarcomeric SHG intensity pattern (SHG-IP) is a single peak per sarcomere i.e. single frequency (1f) in normal physiological condition [START_REF] Recher | Three distinct sarcomeric patterns of skeletal muscle revealed by SHG and TPEF microscopy[END_REF] and that this pattern is converted to two peaks per sarcomere i.e. double frequency (2f) upon proteolysis, mechanical and chemical stresses, laser induced photodamage and muscular dystrophy (Recher et al., 2011a;Recher et al., 2011b;Rouede et al., 2013b). More recently we have shown that these patterns originate from myofibrillar misalignment (Rouede et al., 2013a) and are characterized by triads alteration in mdx mouse model of Duchene muscular dystrophy [START_REF] Rouede | Myofibrillar misalignment correlated to triad disappearance of mdx mouse gastrocnemius muscle probed by SHG microscopy[END_REF] suggesting that SHG microscopy could detect structural disruption of excitation-contraction coupling. However the mechanism underlying 2f sarcomeric SHG-IPs is still not completely elucidated. Several mechanisms have been proposed to explain the ultrastructural origin of this 2f sarcomeric SHG-IP encompassing centrosymmetry of antiparallel myosin tails arrangement at the M-band in normal cardiac muscles [START_REF] Boulesteix | Second-harmonic microscopy of unstained living cardiac myocytes: measurements of sarcomere length with 20-nm accuracy[END_REF][START_REF] Garcia-Canadilla | Automated cardiac sarcomere analysis from second harmonic generation images[END_REF][START_REF] Plotnikov | Characterization of the myosin-based source for second-harmonic generation from muscle sarcomeres[END_REF], disease dependent myofibrillar deformation or change of myosin orientation [START_REF] Friedrich | Microarchitecture is severely compromised but motor protein function is preserved in dystrophic mdx skeletal muscle[END_REF], proteolysis induced intra-sarcomeric misalignment of myosin thick filaments [START_REF] Rouede | Modeling of supramolecular centrosymmetry effect on sarcomeric SHG intensity pattern of skeletal muscles[END_REF] and myofibrillar misalignment (Rouede et al., 2013a). Difficulties in correlating SHG and electron microscopy (EM) images represent the major technological limitation. If EM is the gold standard imaging technique when nanometric resolution is needed, its multiple-step processing of biological samples frequently causes alteration of tissue structure [START_REF] Endo | Entry of a Dye into the Sarcotubular System of Muscle[END_REF][START_REF] Jayasinghe | Three-dimensional reconstruction and analysis of the tubular system of vertebrate skeletal muscle[END_REF][START_REF] Launikonis | Tubular system volume changes in twitch fibres from toad and rat skeletal muscle assessed by confocal microscopy[END_REF][START_REF] Soeller | Examination of the transverse tubular system in living cardiac rat myocytes by 2-photon microscopy and digital image-processing techniques[END_REF] and observations are often limited to ultra-thin slices (about 100 nm) of muscle tissue. In these conditions, correlation between 2f sarcomeric SHG-IP and myofibrillar misalignment at EM level is only possible when 2f sarcomeric SHG-IP is abundant. Based on our previous result showing that muscle proteolysis increases the percentage of 2f sarcomeric SHG-IP (Recher et al., 2011b;[START_REF] Recher | Three distinct sarcomeric patterns of skeletal muscle revealed by SHG and TPEF microscopy[END_REF]Rouede et al., 2013a;[START_REF] Rouede | Modeling of supramolecular centrosymmetry effect on sarcomeric SHG intensity pattern of skeletal muscles[END_REF], we hypothesize that amphibian tadpole tail which undergoes physiological muscle proteolysis with multiple apoptotic pathways during spontaneous metamorphosis [START_REF] Kerr | An electron-microscope study of cell deletion in the anuran tadpole tail during spontaneous metamorphosis with special reference to apoptosis of striated muscle fibers[END_REF], might be a good tissue source of 2f sarcomeric SHG-IPs enabling precise determination of their ultra-structural origin.

Amphibian metamorphosis is characterized by several physiological transformations occurring in a predetermined order 1) proliferation and differentiation of adult tissues triggered by thyroid hormone (TH) signaling and 2) delayed degeneration of many larval tissues [START_REF] Ishizuya-Oka | Amphibian organ remodeling during metamorphosis: insight into thyroid hormone-induced apoptosis[END_REF][START_REF] Nakajima | Regulation of thyroid hormone sensitivity by differential expression of the thyroid hormone receptor during Xenopus metamorphosis[END_REF]. Hence low levels of TH at prometamorphic stages induced hind limb morphogenesis while tail muscle cell death is delayed until metamorphic climax stages when high levels of TH are achieved. High TH level induced spontaneous cell death during tail involution of the anuran xenopus tadpole, is considered as a model for describing and defining cellular events of physiological program cell death (apoptosis) [START_REF] Berry | The expression pattern of thyroid hormone response genes in the tadpole tail identifies multiple resorption programs[END_REF][START_REF] Du Pasquier | Developmental cell death during Xenopus metamorphosis involves BID cleavage and caspase 2 and 8 activation[END_REF][START_REF] Kerr | Apoptosis: a basic biological phenomenon with wide-ranging implications in tissue kinetics[END_REF][START_REF] Kerr | An electron-microscope study of cell deletion in the anuran tadpole tail during spontaneous metamorphosis with special reference to apoptosis of striated muscle fibers[END_REF][START_REF] Rowe | Apoptosis of tail muscle during amphibian metamorphosis involves a caspase 9dependent mechanism[END_REF][START_REF] Thornberry | The caspase family of cysteine proteases[END_REF]. More recently it has been suggested that TH induced reactive oxygen species (ROS) dependent proteolysis is another pathway that plays an important role in amphibian tadpoles tail regression during metamorphosis [START_REF] Inoue | Free radical theory of apoptosis and metamorphosis[END_REF][START_REF] Johnson | Reactive oxygen species and antioxidant defenses in tail of tadpoles, Xenopus laevis[END_REF][START_REF] Kashiwagi | Thyroxine enhancement and the role of reactive oxygen species in tadpole tail apoptosis[END_REF][START_REF] Menon | Oxidative stress, tissue remodeling and regression during amphibian metamorphosis[END_REF]. Despite the general agreement that TH initiates the molecular program leading to apoptosis and tail resorption, the downstream molecular mechanisms and their spatio-temporal characteristic are still not elucidated. Interestingly, we have previously shown in premetamorphic xenopus tadpole that oxidative stress induces tail muscle proteolysis that is characterized by 2f sarcomeric SHG-IPs (Recher et al., 2011a) suggesting that the appearance of this pattern could be used to determine the contribution of ROS dependent proteolysis in climax stage metamorphosis. Altogether both caspase-dependent and ROS-dependent proteolysis processes that occur during xenopus tail muscle involution are supportive of our hypothesis that climax stage xenopus tail muscles might be a rich source of 2f sarcomeric SHG-IPs.

In order to gain insight into the ultrastructural mechanism responsible for 2f sarcomeric SHG-IPs occurring in tail muscles proteolysis during metamorphosis, we undertook SHG, TPEF and electron microscopy (EM) studies from premetamorphic, prometamorphic and metamorphic climax stage of xenopus tadpoles.

Experimental methods

Tissue preparation

Nieuwkoop and Faber (NF) premetamorphic stage 50-53, prometamorphic stage 55-57 and metamorphic climax stage 60-63 Xenopus laevis tadpoles (national breeding facility of xenopus animals in Rennes, France) were anesthetized in MS222 (0.5 mg.mL -1 ). At least 10-15 tadpoles for each group were used. For SHG imaging, fresh tails were cut and mounted in an imaging chamber stabilized and sealed between two coverslips (POC System, Pecon), in Mark's Modified Ringer (MMR). For metamorphic climax stage tadpoles, only tissues from the distal third of their tails, characterized by advanced regression, were used for imaging.

For α-actinin and myosin immunohistochemistry experiments, MS222-euthanized whole tadpoles were fixed at 4•C overnight in MMR containing 4% paraformaldehyd (PFA) then rinsed at least three times with MMR. 10-15 µm cryostat sections of tail tissues, mounted on gelatin-coated coverslips, were permeabilized with buffer containing 0.5% triton and 5% fetal calf serum for 10 minutes. They were then incubated at 4°C over night in MMR containing 0.1% triton, 5% fetal calf serum and one of the following primary antibody. Sarcomeric Zlines were targeted with α-actinin antibody (1:100, mouse monoclonal IgM, ab9465, Abcam, Cambridge, MA, USA) and A-bands were targeted with myosin antibody (1:60, MF20 mouse monoclonal IgG developed by D. A. Fischman or 1:50, F59 mouse monoclonal IgG developed by F.E. Stockdale, both obtained from the Developmental Studies Hybridoma Bank, University of Iowa). Slices were rinsed three times in MMR and incubated at 4°C over night in Ringer containing 0.1% triton, 5% fetal calf serum and secondary antibody (1:100, Alexa Fluor 594 goat antimouse IgG, Molecular Probes, Eugene, OR, USA). After being rinsed three times, the tissue was covered with a drop of mounting medium (Vectashield, Vector Burlingame, CA) and a second coverslip was laid on. After being sealed, preparations were imaged. For transmission electron microscopy (TEM), blocs of approximately 1 mm 3 from PFAprefixed 3 premetamorphic stage 53 and 3 metamorphic climax stage 63 Xenopus laevis tadpoles tail muscles were cut and PBS-washed. The blocs were further fixed in glutaraldehyde (2.5%, 4°C, in 0.2 M cacodylate buffer, pH 7.2) for one hour and washed in cacodylate buffer. The next steps of sample preparation leading to TEM imaging were previously reported [START_REF] Recher | Skeletal muscle sarcomeric SHG patterns photo-conversion by femtosecond infrared laser[END_REF]. Briefly, blocs were post-fixed in osmium tetroxide (2%) and dehydrated in a graded acetone series. They were impregnated firstly into a mix of Epon, Araldite and acetone and secondly into a mix of Epon, Araldite and DMP30 then laid down into molds and polymerized. After bloc desiccation, ultra fine slices were made with an ultra-microtome (Ultracut E, Reichert), mounted on a copper grid and stained with uranyl acetate (2.5%). Slices were then dipped in lead citrate and finally rinsed. Imaging was performed using a JEOL 100CXII microscope at 80 kV accelerating voltage, at the TEM mRIC platform (http://microscopie.univ-rennes1.fr/).

Imaging system

Images were acquired on PIXEL (http://pixel.univ-rennes1.fr/) facility of GIS EUROPIA, (University of Rennes1, France) based on a Leica TCS SP2 confocal scanning head coupled to a DMIRE2 inverted microscope and equipped with an IR 80 MHz femtosecond laser (MAITAI, Spectra Physics). A laser beam tuned at 940 nm was focalized through a 60X water immersion objective (LUMFl 60XW, NA = 1.1, Olympus). The pulse width at the sample is typically 250 fs. A motorized rotation stage (PR50CC, Newport) equipped with an achromatic zero-order Quartz-MgF 2 half-wave plate was used in order to adjust the mean power of the laser beam at the sample. To avoid photodomage (Recher et al 2011) the incident power never exceed 20 mW in living tail muscle tissue. However, for fixed and thin muscle tissue (used for imunolabeling), 5-10 mW was used to obtain well-contrasted and stable SHG images. This power always induces photo bleaching of the TPEF signal. SHG signal was collected in a forward direction using 0.90 NA multi-immersion condenser (model 501000, Leica). BG39 bandpass and a 470 ± 5 nm interference filters were placed in front of the PMT. For dual SHG and TPEF imaging a laser beam tuned at 860 nm was used to enable efficient excitation of Alexa Fluor 594 goat antimouse IgG antibody. BG39 bandpass and 430 ± 5 nm interference filters were placed in front of the transmission PMT to record the forward SHG signal while BG39 bandpass and 600-700 bandpass filters were placed in front of the backward PMT to record TPEF signal.

Image analysis and theoretical simulations

TPEF and SHG images were processed by a fast Fourier transform filter and were analyzed using open source ImageJ software (http://rsb.info.nih.gov/ij/). As previously reported, the percentage of 2f sarcomeric SHG-IPs was quantified for each frame using either pixel intensity profile analysis along myofibrils main axis or a FFT analysis followed by a two-band pass Gabor filtering method [START_REF] Recher | Three distinct sarcomeric patterns of skeletal muscle revealed by SHG and TPEF microscopy[END_REF].

Quantification of myofibrillar misalignment was undertaken by measuring the average misregistered distance separating adjacent sarcomeric A-bands from random chosen TEM images of both premetamorphic and metamorphic Xenopus laevis tadpoles tail muscles.

Theoretical calculation of SHG and TPEF intensities emitted by a bundle of myofibrils has been done elsewhere (Rouede et al., 2013a;Rouede et al., 2013b;[START_REF] Rouede | Modeling of supramolecular centrosymmetry effect on sarcomeric SHG intensity pattern of skeletal muscles[END_REF][START_REF] Rouede | Myofibrillar misalignment correlated to triad disappearance of mdx mouse gastrocnemius muscle probed by SHG microscopy[END_REF]. It is obtained by summation of the amplitude (SHG) or the intensity (TPEF) emitted by all nonlinear sources within the excitation volume (PSF) for each position of the laser beam along the sarcomere. Eq. 3 of [START_REF] Rouede | Myofibrillar misalignment correlated to triad disappearance of mdx mouse gastrocnemius muscle probed by SHG microscopy[END_REF] [START_REF] Rouede | Myofibrillar misalignment correlated to triad disappearance of mdx mouse gastrocnemius muscle probed by SHG microscopy[END_REF] is used for the calculation of theoretical SHG, Z-band and A-band TPEF intensity patterns.

Results

Fragmentation of muscle fibers into a number of so-called apoptotic bodies with wellpreserved cross-striations is an histological characteristic change occurring in metamorphosing tail muscles at climax stages 60-64 [START_REF] Kerr | An electron-microscope study of cell deletion in the anuran tadpole tail during spontaneous metamorphosis with special reference to apoptosis of striated muscle fibers[END_REF]. Typical Z-stack sequence of SHG images taken from live 63 climax stage tadpole tail muscles, within these apoptotic bodies is illustrating in panels A and B of Fig. 1.

Several disorganized fragments of myofibers (dotted line) and well-preserved myofiber (continuous line) can be seen respectively in the ventral and dorsal regions (Fig. 1B). Corresponding sarcomeric SHG-IP are 2f (ventral) and 1f (dorsal) as indicated in Fig. 1C. In apoptotic bodies the percentage of 2f sarcomeric SHG-IP was found to be 97,7 ± 0,5 % (n= 64 images). SHG intensity was very high in apoptotic bodies probably due to an increased condensation of sarcomeric A-bands (see Fig. 2). They are also regions within apoptotic bodies where striation pattern is lost indicating advanced proteolysis (data not shown).

It has been reported at EM level that, following fiber fragmentations, apoptotic bodies of amphibian tail muscle were invaded by macrophages, engulfed into their vacuoles called phagosomes where proteolysis leads to progressive loss of structural organization [START_REF] Kerr | An electron-microscope study of cell deletion in the anuran tadpole tail during spontaneous metamorphosis with special reference to apoptosis of striated muscle fibers[END_REF]. Since apoptotic bodies are characterized by predominant 2f sarcomeric SHG-IP due to drastic proteolysis, we undertook EM studies in order to determine the structural origin of such pattern. Typical EM image of apoptotic body from metamorphic climax stage 63 xenopus tail is illustrated in Fig. 2. A membrane-delimited vesicle or vacuole containing several fragments of myofibrils and a cluster of dead mitochondria can be observed inside the cytoplasm of a macrophage. It becomes evident, at higher magnification, and after a 50degree counter-clock rotation of the whole image that fragments of myofibrils are not registered at their Z-bands (Fig. 2B). We found that adjacent sarcomeres are misaligned by more than 150 nm (Fig. 2C). We found that the average myofibrillar misalignment was 277 ± 20 nm (n = 28 images). As a consequence, triads composed of the association between a central T-tubule and two SR cisternae are completely absent. SR is also absent between condensed myofibrils. Altogether EM images obtained from phagosomes show a good correlation between predominant 2f sarcomeric SHG-IP, myofibrillar misalignment and triads disappearance in agreement with results reported for mdx muscle [START_REF] Rouede | Myofibrillar misalignment correlated to triad disappearance of mdx mouse gastrocnemius muscle probed by SHG microscopy[END_REF]. It has also been shown that TH induced spatial dilatation of SR is the mechanism that precedes and generates muscle fragmentation and apoptotic bodies [START_REF] Kerr | An electron-microscope study of cell deletion in the anuran tadpole tail during spontaneous metamorphosis with special reference to apoptosis of striated muscle fibers[END_REF]. To study the effect of local SR dilatation prior muscle fragmentation on the appearance of 2f sarcomeric SHG-IP, we compare the percentage of this pattern on 3 characteristic different developmental stages. To this aim, we choose stages 50-53 premetamorphic (no TH), stages 55-57 prometamorphic (low level of TH induced hindlimb growth) and metamorphic climax stages 60-63 (high level of TH induced tail resorption) xenopus tadpoles. The results are shown in Fig. 3A. We found that stages 50-53 premetamorphic and stages 55-57 prometamorphic tadpoles non-fragmented muscles were characterized by scarce 2f sarcomeric SHG-IP (image not shown) in agreement with our previous report (Recher et al., 2011a). Quantifications of the percentage of 2f sarcomeric SHG-IP in non-fragmented muscles (Fig. 3A) were found to be respectively 7,1 ± 0,9 % (n = 214 images) and 7,3 ± 1,4 % (n = 82 images) and were not significantly different. However, drastic and significant increase in the proportion of 2f sarcomeric SHG-IP (66,1 ± 3,0 %; n = 124 images) was observed in nonfragmented muscles at metamorphic climax stage 60-63 at time where TH has been shown to induce drastic increase in TH receptors (TR) expression in tail muscles [START_REF] Fujimoto | Expression of matrix metalloproteinase genes in regressing or remodeling organs during amphibian metamorphosis[END_REF]. This latter result is illustrated in Fig. 3B as a 3D orthogonal views from 6 µm thick Zstack. The XY image (left panel) shows a mixed 1f and 2f sarcomeric SHG-IP located respectively on top and bottom part of the image. Two orthogonal XZ slices (b1) and (b2) at horizontal lines 1 and 2 in (b0) show that 2f sarcomeric SHG-IP is partially or fully observed through the 6 µm thickness of the stack. Individual myofibrils could not be resolved in the YZ slices (b3) and (b4) corresponding to transverse section through adjacent sarcomeres. Therefore no information concerning SHG-IP modification was available in YZ views. Based on the drastic increase in 2f sarcomeric SHG-IP in XY and XZ views of SHG images of non-fragmented muscle at climax stage, we were confident that its structural origin could be determined at EM level. To this aim, we compared EM images of non-fragmented muscles from premetamorphic and metamorphic climax stage tadpoles as illustrated in Fig. 3C. Whereas the former exhibit well registered myofibrils with clear visible triads, the latter were characterized by misaligned adjacent sarcomeres with disrupted and severed triads. The frequency of misaligned adjacent sarcomeres was significantly greater (p < 0,001) in metamorphic climax compared to premetamorphic stage tadpoles (Fig. 3D). The average misaligned values were found to be respectively 108 ± 10 nm (n = 144 images) and 21 ± 3 nm (n = 194 images). Interestingly at these metamorphic climax stages the SR exhibits swelling suggesting that this dilatation might be the cue of triad disruptions leading to myofibrillar misalignment and appearance of 2f sarcomeric SHG-IP.

Altogether these results indicate that 2f sarcomeric SHG-IP is well correlated to myofibrillar misalignment and triads damage observed at EM level in non-fragmented tail muscles of metamorphic climax stage tadpoles. Moreover 2f sarcomeric SHG-IP appeared to be a signature of SR swelling that induced triad disruptions and myofibrillar misalignment.

Since disruption of triads resulting in myofibrillar misalignment may have detrimental consequence to excitation-contraction coupling with either muscle weakness or disease, mapping and quantification of label-free 2f sarcomeric SHG-IPs may emerge as powerful biomedical probe. We next compare sensitivity of SHG over TPEF microscopy to detect myofibrillar misalignments. To this aim we use sarcomeric Z-band α-actinin and A-band myosin fluorescent probes. For normal tail muscles of stage 57 xenopus tadpoles, we found that TPEF signals from α-actinin (Z-band marker) generate banding patterns alternating with that of SHG (originating from A-bands) as illustrated in Fig. 4A-C in agreement with our previous result in normal mouse [START_REF] Rouede | Myofibrillar misalignment correlated to triad disappearance of mdx mouse gastrocnemius muscle probed by SHG microscopy[END_REF] and xenopus [START_REF] Recher | Three distinct sarcomeric patterns of skeletal muscle revealed by SHG and TPEF microscopy[END_REF] gastrocnemius muscles. SHG signals (Fig. 4D-I) were found to fully colocalize with myosin-based TPEF from MF20 antibody targeting myosin tails [START_REF] Shimizu | Axial arrangement of the myosin rod in vertebrate thick filaments: immunoelectron microscopy with a monoclonal antibody to light meromyosin[END_REF]. A discrepancy was observed between SHG and myosin-TPEF for F59 antibody targeting myosin cross bridges [START_REF] Miller | Evolutionarily conserved sequences of striated muscle myosin heavy chain isoforms. Epitope mapping by cDNA expression[END_REF]. In this latter case TPEF signal was reduced at the M-bands (centers of myosin thick filaments voided of myosin cross bridges) as expected from the optical resolution of the microscope (0.4 µm). Altogether the localizations of myosin-based TPEF signals are in agreement with their structural origin. For SHG signal a paradoxical bright intensity was observed at the Mband despite the centrosymmetry resulting from antiparallel organization of myosin tail. This over intensity of the SHG signal at the M-band was suggested to originate from particular constructive interferences between harmonic waves (Rouede et al., 2013a;[START_REF] Rouede | Modeling of supramolecular centrosymmetry effect on sarcomeric SHG intensity pattern of skeletal muscles[END_REF]. For proteolyzed non-fragmented muscle fibers of metamorphic climax stages xenopus tails, several regions of 2f sarcomeric SHG patterns are correlated to misaligned Z-bands revealed by α-actinin-based TPEF signal (see large box ROIs in Fig. 5A-C). We also noticed that α-actinin TPEF signals partially reveal regions of 2f sarcomeric SHG in agreement with our previous report on mdx mouse muscle [START_REF] Rouede | Myofibrillar misalignment correlated to triad disappearance of mdx mouse gastrocnemius muscle probed by SHG microscopy[END_REF]. However, myosin-based TPEF signal originating from myosin tails labeled with MF20 antibody could not reveal regions characterized by 2f sarcomeric SHG-IP (indicated by arrows in Fig. 5D-F). In apoptotic bodies from late climax metamorphic tail muscles (Fig. 6A-C), abundant Zband misalignments revealed by α-actinin TPEF signal were often correlated with 2f sarcomeric SHG-IP (Fig. 6D-F). Moreover, we also found regions with 2f sarcomeric SHG-IP without α-actinin TPEF signal indicating complete proteolysis of sarcomeric Z-band αactinin preceding proteolysis of A-band (see rectangle ROIs in Fig. 6D-E).

Altogether these results demonstrate that SHG-IP is more sensitive than α-actinin TPEF signal to report myofibrillar misalignment whereas myosin based TPEF signal is not suitable to report this misalignment.

Discussion

In this report, we show that the canonical 1f sarcomeric SHG-IP of pretamorphic xenopus tadpole tail muscles is converted to 2f sarcomeric SHG-IP at metamorphic climax stages due to massive physiological muscle proteolysis resulting in myofibrillar misalignment observed with EM. We have also compared for the first time, myosin based TPEF and SHG intensity patterns in normal and metarmophosing tail muscles leading to the conclusion that SHG microscopy is more sensitive than TPEF microscopy to detect and to map sarcomere pattern disruption.

The drastic increase in 2f sarcomeric SHG-IP occurring at metamorphic climax stages 60-63 in tail muscle of xenopus tadpoles coincides with the increase in both TH level and TR expression [START_REF] Fujimoto | Expression of matrix metalloproteinase genes in regressing or remodeling organs during amphibian metamorphosis[END_REF][START_REF] Nakajima | Regulation of thyroid hormone sensitivity by differential expression of the thyroid hormone receptor during Xenopus metamorphosis[END_REF] that trigger muscle proteolysis and program cell death. This result is in agreement with our previous report showing that post mortem proteolysis of adult xenopus and rat gastrocnemius muscles results in preponderant 2f sarcomeric SHG-IP [START_REF] Recher | Three distinct sarcomeric patterns of skeletal muscle revealed by SHG and TPEF microscopy[END_REF]. Predominance of 2f sarcomeric SHG-IP in apoptotic bodies of live xenopus tails muscle reported herein was crucial for the ultrastructural determination of their origin at EM level. This has enabled the successful correlation of 2f sarcomeric SHG-IPs to EM myofibrillar misalignment in agreement with our previous report in mdx muscle [START_REF] Rouede | Myofibrillar misalignment correlated to triad disappearance of mdx mouse gastrocnemius muscle probed by SHG microscopy[END_REF].

Whereas in human Duchenne muscular dystrophy (DMD) muscles, hyper contraction and stretching of the muscle at different locations were responsible for the breakdown of the intermediate filament protein desmin that connect adjacent sarcomeric Z-band [START_REF] Cullen | The distribution of desmin and titin in normal and dystrophic human muscle[END_REF], in metamorphosing xenopus tadpole, a different mechanism of myofibrillar misalignment is involved. It has been suggested that muscle fragmentation into apoptotic bodies is due to dilatation and confluence of elements of SR triggered by high levels of TH signaling [START_REF] Kerr | An electron-microscope study of cell deletion in the anuran tadpole tail during spontaneous metamorphosis with special reference to apoptosis of striated muscle fibers[END_REF]. SR swelling was clearly observed in non-fragmented muscles of xenopus climax tadpoles (Fig. 3) and is sufficient to drive T-tubules disruption and myofibrillar misalignment that we observed both in fragmented and non-fragmented muscles. Altogether the result suggests that the increase in 2f sarcomeric SHG-IP in non-fragmented muscles (66%) might precede their appearance in apoptotic bodies (97%) where further amplification of myofibrillar misalignment is generated by macrophage induced phagocytosis and efficient proteolysis. Since the cascade of molecular events linking TH signaling to SR dilatation is not elucidated, 2f sarcomeric SHG-IP offers a new method to decipher the underlying mechanism. Interestingly, we have previously shown in premetamorphic xenopus tail muscles that 1f sarcomeric SHG-IP could be converted into 2f sarcomeric SHG-IP by ROS (Recher et al., 2011a). Also, a good correlation between a chronic increase in concentration of ROS [START_REF] Tidball | The role of free radicals in the pathophysiology of muscular dystrophy[END_REF] and an increase in 2f sarcomeric SHG-IP [START_REF] Buttgereit | From chaos to split-ups--SHG microscopy reveals a specific remodelling mechanism in ageing dystrophic muscle[END_REF][START_REF] Friedrich | Microarchitecture is severely compromised but motor protein function is preserved in dystrophic mdx skeletal muscle[END_REF][START_REF] Garbe | Automated multiscale morphometry of muscle disease from second harmonic generation microscopy using tensor-based image processing[END_REF]Rouede et al., 2013a;[START_REF] Rouede | Myofibrillar misalignment correlated to triad disappearance of mdx mouse gastrocnemius muscle probed by SHG microscopy[END_REF] can be established in mdx muscular dystrophy. For metamorphic climax stage of xenopus tadpoles, the mechanism of TH induced ROS production is not elucidated. TH can exert direct action on mitochondria uncoupling and generate ROS, but it can also up regulate the expression of NO synthase enzymes leading to NO induce mitochondria ROS production [START_REF] Inoue | Free radical theory of apoptosis and metamorphosis[END_REF][START_REF] Johnson | Reactive oxygen species and antioxidant defenses in tail of tadpoles, Xenopus laevis[END_REF][START_REF] Kashiwagi | Thyroxine enhancement and the role of reactive oxygen species in tadpole tail apoptosis[END_REF][START_REF] Menon | Oxidative stress, tissue remodeling and regression during amphibian metamorphosis[END_REF].

In this study, we have also compared for the first time, myosin based TPEF and SHG intensity patterns in normal and metarmophosing tail muscles. For normal muscle, we found a good correlation between TPEF-IP signals based on primary MF20 monoclonal antibody targeting myosin tails [START_REF] Shimizu | Axial arrangement of the myosin rod in vertebrate thick filaments: immunoelectron microscopy with a monoclonal antibody to light meromyosin[END_REF] and 1f SHG-IP (Fig. 4). However, for F59 monoclonal antibody targeting only myosin cross-bridge [START_REF] Miller | Evolutionarily conserved sequences of striated muscle myosin heavy chain isoforms. Epitope mapping by cDNA expression[END_REF] with no signal emitted at the M-band, the expected TPEF intensity patterns were 2f due to the good optical resolution (400 nm) of the imaging system. 1f sarcomeric SHG-IP in thick muscle tissue has been previously explained to be based on particular coherent property of the diffracted SHG signal taking into account the centrosymmetric properties of the M-band at the middle of each sarcomere [START_REF] Rouede | Modeling of supramolecular centrosymmetry effect on sarcomeric SHG intensity pattern of skeletal muscles[END_REF]. For metarmophic climax stages characterized by muscle proteolysis, MF 20 based TPEF-IP remains 1f (Fig. 5) whereas sarcomeric SHG-IP is converted to 2f. In order to understand this discrepancy and to clearly correlate incoherent TPEF and coherent SHG optical signals to myofibrillar misalignment, we undertook theoretical simulation. Theoretical intensity patterns of SHG, MF20-and alpha α-actinin based TPEF signals are obtained by calculation [START_REF] Rouede | Myofibrillar misalignment correlated to triad disappearance of mdx mouse gastrocnemius muscle probed by SHG microscopy[END_REF] of the optical intensity emitted from a bundle of myofibrils centered at x=y=z=0 as indicated in the first panel (I) of (Fig. 7). Laser beam is propagating in the z direction and is focused at the center (z=0) of the bundle with myofibrillar main axis along the x direction. For simplicity of the calculation, myofibrillar section is chosen rectangular (1×1µm) [START_REF] Rouede | Myofibrillar misalignment correlated to triad disappearance of mdx mouse gastrocnemius muscle probed by SHG microscopy[END_REF]. We have considered the case of well-aligned myofibrils (Fig. 7a) and of a single misaligned myofibril (Fig. 7b). We found a good correlation between experimental (Fig. 4 and Fig. 5) and theoretical simulation (Fig. 7). For well-aligned myofibrils (as indicated in panel (a)) intensity patterns are 1f for both SHG (III) and TPEF (IV, V) in agreement with our previous report concerning normal gastrocnemius muscle [START_REF] Rouede | Modeling of supramolecular centrosymmetry effect on sarcomeric SHG intensity pattern of skeletal muscles[END_REF]. For a single myofibril misaligned by ∆ = 0.3 µm as indicated in panel (b), SHG-IP is 2f (III) and TPEF is 1f (IV, V). Such discrepancy originates from respectively the incoherent and coherent properties of the two signals. Because optical resolution (0.4 µm) of our imaging system is greater than the myofibrillar misalignment (∆ = 0.3 µm), it is expected that this misalignment cannot be detected with TPEF and therefore the pattern is 1f. Due to its coherent nature, the SHG process is more sensitive to such small displacement (Rouede et al., 2013a;[START_REF] Rouede | Myofibrillar misalignment correlated to triad disappearance of mdx mouse gastrocnemius muscle probed by SHG microscopy[END_REF]. Theoretical simulation shows that 2f sarcomeric SHG-IP can be obtained for ∆ ≥ 0.2 µm. This result clearly demonstrates that myofibrillar displacement below the optical resolution observed with EM is responsible for SHG-and TPEF-IPs of figure 5 highlighting the supra resolution power of SHG microscopy. We can also notice a good correlation between predicted misalignment (∆ ≥ 0.2 µm) that induces 2f sarcomeric SHG-IP and observed average misalignment at EM level (∆ EM ≥ 0.28 µm, Fig. 2) in apoptotic bodies. For non-fragmented muscles, the observed average misalignment at EM level (∆ EM = 0.11 µm, Fig. 3) is probably under estimated and such discrepancy is probably due to the differential effect of EM process on tissues with different myofibrillar density. Theoretical simulation shows that for ∆ ≥ 0.5 µm sarcomeric α-actinin based intensity pattern becomes 2f while MF20 remains 1f due respectively to the discontinuity of adjacent Z-bands and to the lack of discontinuity of adjacent A-bands. This study therefore highlights the power of SHG over TPEF signal to map myofibrillar misalignment. Hence, being label-free, SHG is an intravital high contrast imaging modality that is more sensitive than TPEF in reporting myofibrillar misalignment.

Myofibrillar misalignment described herein in metamorphic climax stage xenopus tail muscles at EM level has been previously reported in human dystrophic and DMD muscles [START_REF] Cullen | The distribution of desmin and titin in normal and dystrophic human muscle[END_REF] and has been suggested to underlie the enhanced force loss observed after eccentric contractions [START_REF] Blaauw | Eccentric contractions lead to myofibrillar dysfunction in muscular dystrophy[END_REF]. In addition to the correlation between 2f sarcomeric SHG-IP and A-band misalignment, we found a second correlation between 2f sarcomeric SHG-IP and triad reduction in agreement with our previous study in mdx muscle [START_REF] Rouede | Myofibrillar misalignment correlated to triad disappearance of mdx mouse gastrocnemius muscle probed by SHG microscopy[END_REF]. These latters are composed of a central T-tubule, flanked by two closely apposed sarcoplasmic reticulum terminal cisternae [START_REF] Porter | Studies on the endoplasmic reticulum. III. Its form and distribution in striated muscle cells[END_REF] and are essential for physiological excitation-contraction coupling. In normal muscle, triads and myofibrils are regularly ordered in adult skeletal muscles [START_REF] Franzini-Armstrong | Discrimination between fast-and slow-twitch fibres of guinea pig skeletal muscle using the relative surface density of junctional transverse tubule membrane[END_REF][START_REF] Page | A comparison of the fine structures of frog slow and twitch muscle fibers[END_REF][START_REF] Porter | Studies on the endoplasmic reticulum. III. Its form and distribution in striated muscle cells[END_REF][START_REF] Veratti | Investigations on the fine structure of striated muscle fiber read before the Reale Istituto Lombardo, 13 March 1902[END_REF] and are developmentally regulated [START_REF] Ezerman | Differentiation of the sarcoplasmic reticulum and T system in developing chick skeletal muscle in vitro[END_REF][START_REF] Veratti | Investigations on the fine structure of striated muscle fiber read before the Reale Istituto Lombardo, 13 March 1902[END_REF]. Structural intimate relation between the two was first observed at the Z-bands level in amphibian tadpole tail muscles [START_REF] Porter | Studies on the endoplasmic reticulum. III. Its form and distribution in striated muscle cells[END_REF], suggesting that the functional efficiency of triads could be coupled to structural perfect myofibrillar alignment at their A-, M-, I-and Z-bands. We suggest that 2f sarcomeric SHG-IP generated by A-bands misalignment should reflect structural alteration of triads in adult skeletal muscle fibers. This is confirmed herein by the strong correlation between observed triad alteration at EM level in apoptotic bodies and non-fragmented tail muscles of climax stage xenopus tadpole and the corresponding preponderant 2f sarcomeric SHG-IPs.

This study confirms and extends our previous report demonstrating that 2f sarcomeric SHG-IP are ultra-structurally correlated to myofibrillar misalignment in mdx mouse model of human DMD [START_REF] Rouede | Myofibrillar misalignment correlated to triad disappearance of mdx mouse gastrocnemius muscle probed by SHG microscopy[END_REF]. Moreover this study suggests similarities between physiological xenopus tadpole tail muscle resorption and DMD. Like in mdx mouse, 2f sarcomeric SHG-IP in xenopus tadpoles tail resorption is probably a ROS dependent signature of proteolysis induced by cytosolic TH signaling. Unlike the milder proteolysis and muscle damage of mdx mouse model compared to DMD boys [START_REF] Partridge | The mdx mouse model as a surrogate for Duchenne muscular dystrophy[END_REF], xenopus tail resorption is characterized by massive muscle proteolysis triggered by hormonal TH signaling. Therefore xenopus resorption tail could be an alternative organism model aimed at exploring new therapy for DMD. The advantages of xenopus tadpoles over mouse encompasses 1) thousand of embryos per mating 2) external embryonic development 3) optical transparency for SHG imaging 4) easy access to bath apply pharmacological agents 5) study of gene function in developmental processes.

Conclusions

In this report, we show that the canonical single frequency sarcomeric SHG-IP of pretamorphic xenopus tadpole tail muscles is converted to double frequency sarcomeric SHG-IP in metamorphic climax stages due to massive muscle proteolysis induced myofibrillar misalignment observed at EM level. Using dual myosin-based SHG and TPEF we demonstrate experimentally, for the first time to our knowledge, the better sensitivity of the coherent SHG over the incoherent TPEF signal to reveal molecular organization in muscle tissue. Theoretical simulation shows that a single myofibril misalignment as small as 200 nm produces 2f sarcomeric SHG-IP.

Whereas 2f sarcomeric SHG-IP is predominantly observed in either pathological or experimental induced muscle proteolysis, we report for the first time that physiological metamorphosis induces predominance of 2f sarcomeric SHG-IP in xenopus tadpole tail muscle. We also suggest that 2f sarcomeric SHG-IP, which is correlated to myofibrillar misalignment, can be used as a physiological marker of muscle oxidation, proteolysis and disruption of excitation-contraction coupling.

From this study, we suggest that 2f sarcomeric SHG-IP could be used as signature of triad defect and disruption of excitation-contraction coupling and that xenopus tadpole tail resorption could be use as an alternative or complementary model of Duchenne muscular dystrophy.

Figure 1 :

 1 Figure 1: Typical SHG images illustrating fragmentation of muscle fibers occurring at metamorphic climax stage 63 in live xenopus tail muscle. (A) Sequential SHG images from a Z-stack and spaced by 4 µm as indicated. The indicated scale bar in the first panel is 2 µm and is for all panels. (B) Enlarge view of the image taken at z = 12 µm. Continuous and dotted lines indicates portions of respectively well-preserved and fragmented fibers. Scale bar is 2 µm. (C) Continuous and dotted line curves are SHG intensity profiles in arbitrary units (a.u.) along indicated lines in (B). Note that continuous and dotted line curves correspond respectively to continuous and dotted lines of (B). Note also that sarcomeres of well-preserved and fragmented portions of fibers have respectively 1f and 2f sarcomeric SHG-IP. Note the indicated orientation of the tail in B (D=dorsal; V=ventral; R=rostral; C=caudal).

Figure 2 :

 2 Figure 2: TEM images of metamorphic climax stage 63 tail muscle fragment and quantification of myofibrillar misalignment. (A) Typical TEM images showing a vesicular packed muscle fragments (apoptotic body) inside the cytoplasm of a macrophage. Several fragments of myofibrils (mf) and a cluster of dead mitochondria (mt) can be observed. A Nucleus (Nu) with condensed chromatin can also be seen (lower left corner). (B) Higher magnification of the image region indicated by the rectangular ROI in A. Note the counter clock 50 degrees rotation of the ROI. Note also the longitudinal displacement of myofibrils relative to each other, their condensation (lost of sarcoplasmic reticulum), the presence of half sarcomeres with lost of A-band thick filaments (asterisks), absence of triads (association between two sarcoplasmic reticulum cisternae and a central T-tubule). Note also sarcomere fractured at either I-(arrowhead) or A-bands (arrows). Scale bar are 2 µm. (C) Histogram representing frequency distribution of misaligned sarcomeres (sub-class increment is 50 nm).

Figure 3 :

 3 Figure 3: Relationship between SHG and TEM images of non-fragmented tail muscles from different developmental stages xenopus tadpoles. (A) Percentage of 2f sarcomeric SHG-IP of premetamophic (stages 50-53), prometamorphic (stages 55-57) and metamorphic climax (stages 60-63) stages of xenopus tadpoles. (B) 3D orthogonal XY, XZ and YZ views of SHG images from live non-fragmented tail muscle region of metamorphic climax stage 63 tadpole. Panel (b0) represents XY view (25 µm × 25 µm) of the last slice of a Z-stack. Panels (b1) and (b2) are XZ reslice views (25 µm × 6 µm) obtained respectively along horizontal lines 1 and 2 indicated in the XY view (b0). Panels (b3) and (b4) are the YZ reslice views (25 µm × 6 µm) obtained respectively along vertical lines 3 and 4 indicated in the XY view (b0). (C) TEM images from tail muscles, at premetamorphic stage 50 (panels (c1), (c2)) and metamorphic climax stage 63 (panels (c3), (c4)). (c2) and (c4) are enlarged images of ROIs indicated in panels (c1) and (c3). Note that sarcomeres are well registered and triads are preserved in panels (c1) and (c2) from premetamorphic tadpole tail muscle. In contrast myofibril are misaligned and severed triads can be observed in panels (c3) and (c4) from metamorphic climax stage 63. Scale bars = 2 µm. (D) Histogram representing frequency distribution of misaligned sarcomeres from premetamorphic (black bars) and climax stage (grey bars) xenopus tadpoles (sub-class increment is 50 nm).

Figure 4 :

 4 Figure 4: Typical duplex TPEF (A, D, G), SHG (B, E, H) and merge (C, F, I) images of non-metamorphosing tail muscles of xenopus tadpoles. Note that for merge images (C, F, J) TPEF is green and SHG is red. (A, B, C) are respectively multiplex α-actinin TPEF, SHG and merge images from the same field of view of a stage 57 prometamorphic xenopus tadpole. The dotted segments are ROIs delimiting the same Z-bands of contiguous sarcomeres belonging to different myofibrils that cannot be resolved. (D, E, F) are respectively multiplex myosin TPEF (MF-20 antibody), SHG and merge images from the same field of view of a stage 57 prometamorphic xenopus tadpole. The small rectangle boxes are ROIs indicating the same contiguous A-bands of adjacent myofobrils. Note that MF20 antibody labels the myosin molecules tails and therefore produces a good colocalization of TPEF and SHG signals. (G, H, I) are respectively multiplex myosin TPEF (F59 antibody), SHG and merge images from the same field of view

Figure 5 :

 5 Figure 5: Typical duplex TPEF (A, D), SHG (B, E) and merge (C, F) images from α-actinin and myosin of proteolyzed non-fragmented stage 60 xenopus tadpole tail muscles. Note that for merge images (C, F) TPEF is in green and SHG is in red color. (A, B, C) are respectively multiplex α-actinin TPEF, SHG and merge images from the same field of view. Large rectangle ROIs represent the same regions where Z-band misalignments and 2f sarcomeric SHG-IP are correlated. Small rectangle ROIs represent the same regions with registered Z-bands and good alternation of α-actinin TPEF and SHG signals. (D, E, F) are respectively multiplex myosin TPEF (MF-20 antibody), SHG and merge images from the same field of view. Note that arrows indicate examples of region with 2f sarcomeric SHG-IP not observed in corresponding TPEF images. Scale bars are 2 µm.

Figure 6 :

 6 Figure 6: Typical duplex α-actinin TPEF (A, D), SHG (B, E) and merge (C, F) images of tail muscle fragments forming apoptotic bodies from climax stage 64 xenopus tadpole. Note that for merge images (C, F) TPEF is in green and SHG is in red color. Note also the indicated orientation of the tail in C (D=dorsal; V=ventral; R=rostral; C=caudal) and that proteolyzed dorsal dermal collagen is observed only in SHG image. Note also the presence of dorsal row of slow myofibers that are preserved. (D, E, F) are respectively multiplex α-actinin TPEF, SHG and merge images from the same field of view. Note the global correlation between misaligned α-actinin TPEF signal and 2f sarcomeric SHG-IP. Rectangle ROIs represent the same region characterized by disappearance of Z-band α-actinin signal suggesting its complete proteolysis whereas sarcomeric A-band SHG signal is still preserved. Scale bars = 2 µm.

Figure 7 :

 7 Figure 7: Macroscopic diagrams illustrating the better sensitivity of SHG over TPEF to detect myofibrillar misalignment. (aI) and (bI) are macroscopic xyz view representing respectively well aligned myofibrils and single central myofibril misaligned by ∆=0.3 µm. (aII) and (bII) are corresponding xy views at the center of the bundles (z=0). For SHG and TPEF simulations, laser beam is propagating along z direction and is focused at the center of each bundle (z=0). (aIII) and (bIII) are corresponding SHG patterns. (aIV, aV) and (bIV, bV) are corresponding TPEF patterns from respectively α-actinin and myosin tails. Note that macroscopic diagrams (aI, bI) represent bundles of 3 sarcomeres (L=2 µm) from 5×5 myofibrils of rectangular size 1×1 µm in y and z directions. Note also that misalignment (bI, bII) only concerns the central myofibril as indicated. For schematic views (aI, bI, aII, bII), A-bands (1.6 µm) consist of well-ordered thick filaments that are represented by a double color code (blue and red) to account for polarity inversion at the M-band (0.15 µm, grey color) and α-actinin Z-bands (0.1 µm) are in green color. Note that for each pattern (III, IV, V), intensity profile (white line) calculated at y=z=0 is shown in inset.

Acknowledgements

We thank Matthieu Fabrega for his technical help. This work was supported by Région Bretagne, Rennes Métropole, Conseil Général d'Ille-et-Villaine, CRITT Santé Bretagne, CBB Developpement Bretagne, Ministère de l'Enseignement Supérieur et de la Recherche and the European Union Federal Funds FEDER.


Romain Guillaume
email: romain.guillaume@irit.fr

Caroline Thierry
email: caroline.thierry@irit.fr

Paweł Zieliński
email: pawel.zielinski@pwr.wroc.pl

Robust Production Plan with Periodic Order Quantity under Uncertain Cumulative Demands

Keywords: Supply Chain, Production Planning, Uncertainty, Scenario Optimization

In this paper, we are interested in a production planning process in collaborative supply chains. More precisely, we consider supply chains, where actors use Manufacturing Resource Planning process (MRPII). Moreover, these actors collaborate by sharing procurement plans. We focus on a supplier, who applies the Periodic Order Quantity (POQ) rule to plan a production integrating the uncertain procurement plan sent by her/his customer. The uncertainty of the procurement plan is expressed by closed intervals on the cumulative demands. In order to choose a robust production plan, under the interval uncertainty representation, the min-max criterion is applied.

We propose algorithms for determining the set of possible costs of a given production plan -due to the uncertainty on the cumulative demands. We then construct algorithms for computing a robust production plan with respect to the min-max criterion: the algorithm based on iterative adding constraints and the polynomial algorithms under certain realistic assumptions.

I. INTRODUCTION

Companies today evolve in high competitive context that obliges the companies to collaborate with their suppliers and customers and creates uncertainty on the demand. Due to the well-known bullwhip effect [START_REF] Lee | Information Distortion in a Supply Chain: The Bullwhip Effect[END_REF], this uncertainty induces supply chain risks as backordering, obsolete inventory. Sharing information on the demand and the collaboration with the suppliers are ways to reduce this risk.

Most companies use Manufacturing Resource Planning (MRPII) to plan and control all resources of a manufacturing company. MRPII is composed of three processes (the production process, the procurement process and the distribution process) and three levels [START_REF] Arnold | Introduction to Materials Management[END_REF]: the strategic level (Sales and Operation Plan-S&OP), the tactical level (Master Production Scheduling (MPS) and Material Requirement Planning (MRP)) and the operational level (detailed scheduling and shop floor control). Within MRP process, different lot sizing rules exist for purchased or produced items, as Fixed Order Quantity (FOQ), Lot-for-Lot (L4L), Minimal Order Quantity (MOQ), Periodic Order Quantity (POQ), etc. In collaborative supply chains, collaboration is usually characterized by a set of pointto-point (customer/supplier) relationships with partial information sharing. More precisely, the collaboration process in supply chains, where actors use MRPII, is realized by sharing procurement plans through the supply chain. The procurement plan can take into account uncertainty [START_REF] Zouggar | Supply Chain Reactivity Assessment Regarding Two Negotiated Commitments: Frozen Horizon and Flexibility Rate[END_REF]. Thus, the problem is how to integrate this information in a production planning process. In this paper, we focus on cases where no probability distribution is available to model the uncertainty. In this context, the uncertainty is modeled by specifying a set of all possible realizations of the demand, called scenarios.

In the literature, the planning processes of MRPII have been extended to take into account the imprecision on quantities of period demands (MPS and MRP) [START_REF] Fargier | The Use of Possibilistic Decision Theory in Manufacturing Planning and Control: Recent Results in Fuzzy Master Production Scheduling[END_REF], [START_REF] Mula | Material Requirement Planning with fuzzy constraints and fuzzy coefficients[END_REF], [START_REF] Tavakkoli-Moghaddam | A Fuzzy Aggregate Production Planning Model for Make-to-Stock Environments[END_REF], [START_REF] Guillaume | A robust lot sizing problem with ill-known demands[END_REF], quantities of period demands and uncertain orders (MRP) [START_REF] Grabot | Integration of uncertain and imprecise orders in the MRP method[END_REF] and the imprecision on order quantities and dates with uncertain order (MRP) [START_REF] Guillaume | Modelling of ill-known requirements and integration in production planning[END_REF].

To deal with the uncertainty in the production planning, three approaches can be distinguished: computing the possible inventory and backordering levels over all scenarios to help the decision maker to choose a production plan ( [START_REF] Fargier | The Use of Possibilistic Decision Theory in Manufacturing Planning and Control: Recent Results in Fuzzy Master Production Scheduling[END_REF], [START_REF] Grabot | Integration of uncertain and imprecise orders in the MRP method[END_REF], [START_REF] Guillaume | Modelling of ill-known requirements and integration in production planning[END_REF]), computing an optimal solution for one of possible demand scenarios [START_REF] Mula | Material Requirement Planning with fuzzy constraints and fuzzy coefficients[END_REF], [START_REF] Tavakkoli-Moghaddam | A Fuzzy Aggregate Production Planning Model for Make-to-Stock Environments[END_REF] and the robust optimization [START_REF] Kouvelis | Robust Discrete Optimization and its applications[END_REF] under the scenario uncertainty representation, more precisely using the min-max criterion [START_REF] Guillaume | A robust lot sizing problem with ill-known demands[END_REF]. Under this criterion, we seek a solution that minimizes the largest cost over all scenarios. The cost function in production planning is the sum of inventory and backordering costs over the planning horizon.

The aim of this paper is to investigate the MRP process with the POQ rule under imprecision on cumulative demands. The POQ rule consists in producing a quantity equal to the gross requirements for periods minus any items in onhand inventory plus any additional items needed to replenish safety stock if it has fallen below its desired level. To adapt the MRP with the POQ rule to the uncertain context, we have to consider the problem with backordering. Indeed, the problem without backordering is not satisfactory due to the fact that a solution method (without backordering) consists in applying the rule to the maximal cumulative demands. Thus, this solution method induces too much inventory and does not consider preferences of the decision maker between possible inventory and backordering levels. In the model proposed in this paper, the imprecision on the demand is represented by cumulative demand intervals. Such modeling allows us to describe the imprecision on order quantities and dates.

The paper is organized as follows. Section II presents the problem under consideration with the precise demands (parameters). In Section III, we formulate the problem under the scenario uncertainty model in the robust optimization setting. We adopt min-max criterion to choose a robust production plan. We then propose algorithms for evaluating a given production plan (for determining optimal interval containing all possible values of costs of the production plan) and for computing an optimal robust production plan.

II. THE DETERMINISTIC PROBLEM

In this section, a deterministic version of the problem under consideration, i.e. the problem in which all parameters are precisely known in advance.

Given + 1 periods. For period , = 0, . . . , , let be the demand in period , ≥ 0, the production amount in period . Furthermore, we are given a periodicity , ∈ ℕ, such that:

≥ 0 if = ⋅ ; otherwise (if ∕ = ⋅ )
= 0 for = 0, . . . , and = 0, . . . , where = / (we assume without loss of generality that is divisible by ). Now, the set feasible production amounts ⊆ ℝ +1 ≥0 can be defined as follows:

= { = ( 0 , . . . , ) : ≥ 0 for = ⋅ , = 0 for ∕ = ⋅ , = 0, . . . , , = 0, . . . , }. Set D = ∑ =0 and X = ∑ =0
, D and X stand for the cumulative demand up to period and the production level up to period , respectively. Obviously, X -1 ≤ X and D -1 ≤ D , = 1, . . . , . The nonnegative costs of carrying one unit of inventory from period to period + 1 are given, denoted by , and all the inventory costs are equal for every period. The nonnegative costs of backordering one unit from period + 1 to period are given, denoted by , and all the backorder costs are equal for every period. Furthermore, we assume that ≤ . The nonnegative real function ( , ) represents either the cost of storing inventory from period to period + 1 or the cost of backordering quantity from period + 1 to period , namely

(X , D ) = (X -D ) if X ≥ D ; (D -X ) otherwise. The function has the form (X , D ) = max{ (X -D ), (D -X )}.
The optimization problem with the precise parameters consists in finding a feasible production plan = ( 0 , . . . , ), ∈ , that minimizes the total cost of storage and backordering subject to the conditions of satisfying each demand, that is

min ∈ ( ) = min ∈ ∑ =0 ( ∑ =0 , ∑ =0 ) = min ∈ ∑ =0 (X , D ). (1)
It is easily seen that when = 1, the problem ( 1) is equivalent to the classical lot sizing with backordering with the Lot-For-Lot (L4L) rule (see, e.g., [START_REF] Johnson | Operations Research in Production Planning, Scheduling and Inventory Control[END_REF], [START_REF] Wagner | Dynamic Version of the Economic Lot Size Model[END_REF], [START_REF] Florian | Deterministic Production Planning: Algorithms and Complexity[END_REF]). The problem (1) can be formulated as the minimum cost flow problem (see, e.g., [START_REF] Ahuja | Network Flows: theory, algorithms, and applications[END_REF]):

min ∑ =0 ( + ) s.t. -= ∑ =0 ( -), = 0, . . . , , = 0, ∕ = ⋅ , = 0, . . . , , = 0, . . . , , , , ≥ 0, = 0, . . . , . (2)
The problem (2) can solved in ( ). For each = 0, . . . , -1, we determine periods:

ℎ = max{ : ∈ { ⋅ , . . . , ( + 1) ⋅ -1}, ( -⋅ ) ≤ (( + 1) ⋅ -)}.
An optimal production plan to ( 2) is computed by the following formula for = 0, . . . , and = 0, . . . , :

= ⎧    ⎨    ⎩ 0 if ∕ = ⋅ , D ℎ if = ⋅ and = 0, D ℎ -D ℎ -1 if = ⋅ and 0 < < , D ⋅ -D ℎ -1 if = ⋅ and = .
We have assumed that an initial inventory and an initial backorder are equal to zero. Otherwise, one can easily modify the above method to cope with > 0 or > 0.

III. ROBUST VERSION OF THE PROBLEM

In Section II, we have assumed the all input parameters in problem (1) are precisely known. However, in real life this is rarely the case. Here, we admit uncertainty on the demands.

A. Model of uncertainty

One of the simplest form of the uncertainty representations is modeling the imprecise demands ˜ , = 0, . . . , , as closed intervals [ , ], ≥ 0, where and are a minimal and a maximal possible values of demand ˜ in period , respectively. So, assigning some interval [ , ] to demand ˜ means that it will take some value within the interval, but it is not possible to predict at present which one, i.e.

∈ [ , ]. From the above model of uncertainty, it follows that the imprecision of cumulative demand

D = ∑ =0 increases in subsequent periods, i.e. D ∈ [ ∑ =0 , ∑ =0
]. In fact, practitioners often express the knowledge on demand uncertainty by the range. The demand can be interpreted in two different ways: a demand in the period or a cumulative demand. For instance, a practitioner expresses the uncertainty on the demand by range ±Δ. If this uncertainty is interpreted as the one on demands in periods, then it leads to the cumulative demands with increasing uncertainty (see Fig. 1a)). This case is unrealistic compared to the case when the uncertainty is interpreted as the one on cumulative demands (see Fig. 1b)). So in this paper, the uncertainty of the demands ˜ is described by the uncertainty on the cumulative demands, modeled by intervals [D ; D ], instead of the uncertainty on demands in periods, Among the scenarios of Γ, we distinguish the ones called extreme scenarios. Each extreme scenario = (D ) =0 belongs to the set of scenarios defined by the following recurrence formula:

D ∈ { {D 0 , D 0 } if = 0, {max{D -1 , D }, D } if = 1, . . . , . (3)
We will denoted by Γ ext , the set of extreme scenarios. Clearly, In order to choose a robust production plan, one of robust criteria, called the min-max can be adopted (see, e.g. [START_REF] Kouvelis | Robust Discrete Optimization and its applications[END_REF]).

Γ
In the min-max version of problem (1), we seek a feasible production plan with the minimum the worst total cost over all scenarios, that is

ROB : min ∈ ( ) = min ∈ max ∈Γ ( , ) = min ∈ max ∈Γ ∑ =0 (X , D ( )).
In other words, we wish to find among all production plans the one that minimizes the maximum production plan cost over all scenarios, that minimizes ( ), ( ) is the maximal cost of production plan . An optimal solution to the problem ROB is called optimal robust production plan.

Here and subsequently (as in Section II), we assume that an initial inventory and an initial backorder are equal to zero. Otherwise, one can modify the algorithms presented in this section to cope with the case > 0 or > 0.

B. Evaluating a Given Production Plan

In this section, we will be concerned with evaluating a given production plan. We will propose methods for computing the optimal interval containing all possible values of costs of the production plan.

Let * ∈ be a given production plan. A scenario ∈ Γ that minimizes the total cost ( * , ) of the production plan * is called optimistic scenario. A scenario ∈ Γ that maximizes the total cost ( * , ) of the production plan * is called the worst case scenario. Thus, the optimal interval containing all possible values of costs of the production plan * is of form: [ ( * , ), ( * , )]. We begin with a result on function ( * , ) = ∑ =0 (X , D ( )): Proposition 1: Function ( * , ) is convex on Γ for any fixed production plan * ∈ . Proposition 1 follows by similar arguments as in [START_REF] Guillaume | A robust lot sizing problem with ill-known demands[END_REF], i.e. function (X * -D ( )) and (D ( ) -X * ) are convex on Γ and so max{ (X * -D ( )), (D ( ) -X * )} and ∑ =0 max{ (X * -D ( )), (D ( ) -X * )} are convex.

1) Computing an Optimistic Scenario:

The problem of determining an optimistic scenario = (D ) =0 for a given production plan * ∈ , i.e. the problem

( * , ) = min ∈Γ ( * , ), (4)
can be formulated by a linear programming problem:

min ∑ =0 ( + ) s.t. -= D -X * , = 0, . . . , , D -1 ≤ D , = 1, . . . , , D ≤ D ≤ D , = 0, . . . , . , ≥ 0, = 0, . . . , (5)
If D , and is an optimal solution to (5), then = (D ) =0 is an optimistic scenario for * . Furthermore, since , ≥ 0, either > 0 or > 0, which means that for storing from period to +1 and backordering from period +1 to are not performed simultaneously. However, determining an optimistic scenario for * can be improved to ( ), since one can give the explicit form of an optimistic scenario = (D ) =0 and thus an optimal solution to (5):

D = ⎧  ⎨  ⎩ D if X * < D , X * if D ≤ X * ≤ D , = 0, . . . , . D if X * > D , (6)
The form [START_REF] Tavakkoli-Moghaddam | A Fuzzy Aggregate Production Planning Model for Make-to-Stock Environments[END_REF] follows from inequalities:

D -1 ≤ D , D -1 ≤ D , X * -1 ≤ X * , = 1
, . . . , , and Proposition 1.

2) Computing a Worst Case Scenario:

We now pass on to the problem of computing a worst case scenario = (D ) =0 for a given production plan * ∈ , i.e.

( * , ) = max ∈Γ ( * , ). (7)
Problem ( 7) is more difficult than the one of computing an optimistic scenario and thus a solution algorithm is much more involved. The following proposition allows us to construct an efficient algorithm for the problem under consideration: Proposition 2: A worst case scenario is an extreme one, i.e.

∈ Γ ext . Proof: Proposition follows by the same method as [START_REF] Guillaume | A robust lot sizing problem with ill-known demands[END_REF]. Function ( * , ) attains its maximum in convex set Γ. An easy computation shows that scenarios ∈ Γ ext are the vertices of Γ. From the above and Proposition 1, it follows that ( * , ) attains the maximum value at a vertex of Γ (see, e.g., [START_REF] Martos | Nonlinear programming theory and methods[END_REF]).

We now construct an algorithm for the problem [START_REF] Guillaume | A robust lot sizing problem with ill-known demands[END_REF], based on a dynamic programming technique. Let be the set of feasible cumulative demand levels in period , = 0, . . . , . Namely D ∈ if and only if D is the th component of an extreme scenario that belongs to the set generated by formula (3) (Γ ext ). Let C -1 (D -1 ) be the maximal cost of a given production plan * over periods , . . . , , when the cumulative demand level up to period -1 is equal to

D -1 , D -1 ∈ -1 , C -1 : -1 → ℝ ≥0 . Set -1 = {0}.
We see at once that:

C (D ) = 0, D ∈ , (8)
C -1 (D -1 ) = max ⎧ ⎨ ⎩ (X * , D ) + C (D ) (X * , max{D -1 , D })+ +C (max{D -1 , D }) ⎫ ⎬ ⎭ D -1 ∈ -1 , = , . . . , 0. (9)
The maximal cost of production plan * over period 0, . . . , is equal to C -1 (0), C -1 (0) = ( * , ), which is computed by the backward recursion ( 8) and ( 9). Worst case scenario = (D ) =0 for * can be determined by a forward recursion technique. It is sufficient to store for each D -1 ∈ -1 the value for which the maximum in ( 9) is attained, that is D is either max{D -1 , D } or D .

Let us analyze the running time of the dynamic programming based algorithm. Building the sets of feasible cumulative demand levels 0 , . . . , according to formula (3) and computing ( * , ) by the backward recursion ( 8) and ( 9) can be done in ( ⋅ max =0,..., | |). Determining by by a forward recursion takes ( ). It is easy to check that max =0,..., | | ≤ + 2. Hence, the total running time of the algorithm is ( 2 ).

C. Computing an Optimal Robust Production Plan

In this section, we will propose algorithms for computing an optimal robust production plan to problem ROB. We will first construct an algorithm for problem ROB without any additional restrictions on ROB. We will then put some restrictions on ROB and obtain more efficient algorithms for computing an optimal robust production plan.

1) Algorithm for the General Problem ROB: We now give an iterative algorithm for the problem ROB based on iterative adding constraints for min-max problems proposed in [START_REF] Shimizu | Necessary Conditions for Min-Max Problems and Algorithms by a Relaxation Procedure[END_REF]. Similar methods were developed for min-max regret linear programming problems with an interval objective function [START_REF] Inuiguchi | Minimax regret solution to linear programming problems with an interval objective function[END_REF], [START_REF] Mausser | A heuristic to minimax absolute regret for linear programs with interval objective function coefficients[END_REF] and a production planning problem with interval demands [START_REF] Guillaume | A robust lot sizing problem with ill-known demands[END_REF]. Let us perform a relaxation of the problem ROB that consists in replacing a given cumulative demand scenario set Γ with a discrete scenario set

Γ dis = { 1 , . . . , }, Γ dis ⊆ Γ: RX-ROB: ˆ = min s.t. ≥ ( , ) ∀ ∈ Γ dis , ∈ , (10)
where = (D ) =0 , = ( ) =0 . Since Γ dis ⊆ Γ, ˆ is a lower bound on the maximal cost of an optimal robust production plan to problem ROB, ˆ ≤ ( ) = ( , ). Note that the constraint ≥ ( , ) called scenario cut, associated with is not a linear constraint. Each scenario cut associated with can be linearized by replacing it with the following +2 constraints and 2 +2 new decision variables ( and ):

≥ ∑ =0 ( + ), -= D - ∑ =0 , = 0, . . . , , , ≥ 0, = 0, . . . , . (11)
Replacing each scenario cut ≥ ( , ) by ( 11) in [START_REF] Kouvelis | Robust Discrete Optimization and its applications[END_REF] leads to the following linear program:

ˆ = min s.t. ≥ ∑ =0 ( + ), ∀ ∈ Γ dis , -= D - ∑ =0 , = 0, . . . , , ∀ ∈ Γ dis , , ≥ 0, = 0, . . . , , ∀ ∈ Γ dis , = 0, ∕ = ⋅ , = 0, . . . , , = 0, . . . , , ≥ 0, = 0, . . . , . (12)
The iterative algorithm for the problem ROB (Algorithm 1) starts with zero lower bound = 0, a candidate ˆ ∈ for an optimal solution for ROB (any solution in ) and empty discrete scenario set, Γ dis = ∅. At each iteration, a worst case scenario for ˆ is determined by the dynamic programming algorithm presented in Section III-B2. The value of (ˆ ) = (ˆ , ) is an upper bound on ( ), ( ) ≤ (ˆ ). If a termination criterion is fulfilled (Step 3), for a given precision > 0, then the algorithm stops with production plan ˆ being an approximation of an optimal robust production plan . Otherwise the worst case scenario is added to Γ dis , the scenario cut corresponding to is appended to problem RX-ROB or equivalently to linear programming problem [START_REF] Wagner | Dynamic Version of the Economic Lot Size Model[END_REF] . Next the updated linear programming problem ( 12) is solved to obtain a better candidate ˆ for an optimal robust production plan to problem ROB and new lower bound = ˆ . Since set Γ dis is updated during the course of the algorithm, the computed values of lower bounds {ˆ } form a nondecreasing sequence of their values. Then new iteration is started. Algorithm 1 terminates in a finite number of iterations for any given > 0. In order to show this, the same reasoning applies as those given in [START_REF] Geoffrion | Generalized Benders Decomposition[END_REF]Theorem 2.5], [START_REF] Shimizu | Necessary Conditions for Min-Max Problems and Algorithms by a Relaxation Procedure[END_REF]Theorem 3] and [START_REF] Guillaume | A robust lot sizing problem with ill-known demands[END_REF]Theorem 3]. It is worth pointing out that the running time of each iteration highly depends on Step 2 and Step 5, where a worst case scenario is computed and linear programming problem ( 12) is solved. The running time of Step 2 is

( 2 ), linear program (12) (Step 5) can be solved by using a specially-tuned method for this kind of problems [START_REF] Ahuja | Minimax linear programming problem[END_REF] or by some standard off-the-shelf LP solvers. Hence, each iteration of Algorithm 1 can be done efficiently (in a polynomial time).

2) Algorithms for Special Cases of ROB: Consider the problem ROB, when periodicity = 1. In this case, we can apply a method proposed in [START_REF] Guillaume | A robust lot sizing problem with ill-known demands[END_REF] for computing an optimal robust production plan = ( ) =0 :

X = D + D + , = 0, . . . , . (13)
Set 0 = X 0 and = X -X -1 for = 1, . . . , . Consider the problem ROB, when > 1 (if = 1 one can use [START_REF] Florian | Deterministic Production Planning: Algorithms and Complexity[END_REF]) and the bounds of cumulative demand intervals are such that D -1 ≤ D for every = 1, . . . , . According to the above assumption, the set of cumulative demand scenarios Γ

has the form Γ = [D 0 , D 0 ] × ⋅ ⋅ ⋅ × [D , D ]. Now each = (D ) =0 ∈ Γ has components such that D -1 ≤ D , = 1
, . . . , . Hence, and the periodicity in the problem ROB, it follows that one can decompose the problem into + 1 separate subproblems, that is:

= arg min ∈ max ∈Γ ∑ =0 (X , D ( )) = = -1 ∑ =0 min X ⋅ ∈ max ∈Γ ( +1)⋅ -1 ∑ = ⋅ (X ⋅ , D ( )) + min X ⋅ ∈ max ∈Γ ⋅ (X ⋅ , D ⋅ ( )), (14)
where

= [D ⋅ , D ⋅ ] ∪ ⋅ ⋅ ⋅ ∪ [D ( +1)⋅ -1 , D ( +1)⋅ -1 ], Γ = [D ⋅ , D ⋅ ] × ⋅ ⋅ ⋅ × [D ( +1)⋅ -1 , D ( +1)⋅ -1 ], = 0, . . . , -1, = [D ⋅ , D ⋅ ] and Γ = [D ⋅ , D ⋅ ]; Γ = Γ 0 × ⋅ ⋅ ⋅ × Γ .
Obviously, the above possible cumulative production levels are nondecreasing sequence of their values, X 0 ≤ X ≤ ⋅ ⋅ ⋅ ≤ X ⋅ . Therefore, we need only to solve the + 1 separate subproblems. The last subproblem is trivial, i.e. an optimal cumulative production level X ⋅ can be computed by formula:

X ⋅ = D ⋅ + D ⋅ + . (15)
It remains to solve the subproblems for = 0, . . . , -1, i.e. For each ℎ = ⋅ , . . . , ( + 1) ⋅ -1, we compute a possible cumulative production level

min X ⋅ ∈ max ∈Γ ( +1)⋅ -1 ∑ = ⋅ (X ⋅ , D ( )). (16)
X ℎ ⋅ = D ℎ + D ℎ + ∈ [D ℎ , D ℎ ].
Note that for X ℎ ⋅ equality ℎ (X ℎ ⋅ , D ℎ ) = ℎ (X ℎ ⋅ , D ℎ ) holds and its worst case scenario has form (D ⋅ , . . . , D ℎ-1 , D ℎ , . . . , D ( +1)⋅ -1 ). Thus, the worst case cost for X ℎ ⋅ over the set of scenarios Γ , denoted by ℎ , is as follows:

ℎ (X ℎ ⋅ ) = max ∈Γ ( +1)⋅ -1 ∑ = ⋅ (X ℎ ⋅ , D ( )) = ℎ-1 ∑ = ⋅ (X ℎ ⋅ , D ) + ( +1)⋅ -1 ∑ =ℎ (X ℎ ⋅ , D ).
We then determine an optimal cumulative production level X ⋅ for the th subproblem ( 16) by formula:

X ⋅ = arg min ℎ∈{ ⋅ ,...,( +1)⋅ -1} ℎ (X ℎ ⋅ ). (17)
Using ( 15) and ( 17), one can easily compute an optimal robust production plan . Thus, ( 14) can be solved in ( ).

IV. CONCLUSION In this paper, we have studied the possibility of integrating the uncertainty into planning process. We have proposed efficient algorithms for solving the production planning problem using MRP process with the POQ rule under uncertain cumulative demands, modeled by closed intervals, with the minmax criterion. Namely, the algorithm based on iterative adding constraints. It is worth pointing out that each iteration can done in a polynomial time. Moreover, we have provided algorithms for two special cases of the problem under consideration. The first case corresponds to the Lot-For-Lot (L4L) rule. In the second one the range of uncertainty is small compared to the demand.

An interesting topic for further research is investigating the production planning problem using MRP process under uncertain cumulative demands, considered in this paper, with other rules such as: Fixed Order Quantity (FOQ) and Minimal Order Quantity (MOQ), and with production constraints.

Fig. 1 .

 1 Fig.1.An example: (a) the case with the uncertainty on demands in periods and the resulting uncertain cumulative demands, (b) the case with the uncertainty on cumulative demands.

  ext ⊆ Γ. The cumulative demand and the demand in period under scenario are denoted by D ( ), D ( ) ∈ [D ; D ], and ( ), respectively, ( ) = D ( ) -D -1 ( ). Clearly, for every ∈ Γ it holds D -1 ( ) ≤ D ( ), = 1, . . . , . The function (X , D ( )) = max{ (X -D ( )), (D ( )-X )}, represents either the cost of storing inventory from period to period + 1 or the cost of backordering quantity from period + 1 to period under scenario . Now ( , ) denotes the total cost of a production plan ∈ under scenario , i.e. ( , ) = ∑ =0 (X , D ( )). The set feasible production amounts is the same as in Section II.

Algorithm 1 :

 1 Solving problem ROB.Input: [D , D , ], = 0, . . . , , , initial production plan ˆ ∈ , a tolerance > 0. Output: A production plan, an approximation of an optimal robust production plan, and its worst case scenario. Step 0. := 0, := 0, Γdis := ∅. Step 1. := ˆ . Step 2. Compute a worst case scenario for by the dynamic programming algorithm presented in Section III-B2. Step 3. Δ := ( , ) -. If > 1 then Δ := Δ/ . If Δ ≤ then output , and STOP. Step 4. := + 1, := , Γdis := Γdis ∪ { } and append scenario cut ≥ ( , ) to problem RX-ROB. Step 5. Compute an optimal solution (ˆ , ˆ ) to RX-ROB (linear programming problem (12)), := ˆ , and go to Step 1.

  Consider the th subproblem. Since X = X ⋅ in period , ∈ { ⋅ + 1, . . . , ( + 1) ⋅ -1}, worst case scenarios ∈ Γ for X ⋅ belong to the set[D ⋅ , D ⋅ ] × ⋅ ⋅ ⋅ × [D ℎ-1 , D ℎ-1 ] × [D ℎ , D ℎ ] × [D ℎ+1 , D ℎ+1 ] × ⋅ ⋅ ⋅ × [D ( +1)⋅ -1 , D ( +1)⋅ -1 ] ⊆ Γ , assuming that X ⋅ ∈ [D ℎ , D ℎ ].