H Benhabiles-Mezhoud

J Kiener

V Tatischeff

A W Strong

Deexcitation nuclear gamma-ray line emission from low-energy cosmic rays in the inner Galaxy

Recent observations of high ionization rates of molecular hydrogen in diffuse interstellar clouds point to a distinct low-energy cosmic-ray component. Supposing that this component is made of nuclei, two models for the origin of such particles are explored and low-energy cosmic-ray spectra are calculated which, added to the standard cosmic ray spectra, produce the observed ionization rates. The clearest evidence of the presence of such low-energy nuclei between a few MeV per nucleon and several hundred MeV per nucleon in the interstellar medium would be a detection of nuclear γ-ray line emission in the range E γ ∼ 0.1 -10 MeV, which is strongly produced in their collisions with the

interstellar gas and dust. Using a recent γ-ray cross section compilation for nuclear collisions, γ-ray line emission spectra are calculated alongside with the high-energy γ-ray emission due to π 0 decay, the latter providing normalization of the absolute fluxes by comparison with Fermi-LAT observations of the diffuse emission above E γ = 0.1 GeV. Our predicted fluxes of strong nuclear γ-ray lines from the inner Galaxy are well below the detection sensitivies of INTEGRAL, but a detection, especially of the 4.4-MeV line, seems possible with new-generation γ-ray telescopes based on available technology. We predict also strong γ-ray continuum emission in the 1-8 MeV range, which in a large part of our model space for low-energy cosmic rays exceeds considerably estimated instrument sensitivities of future telescopes.

Subject headings: gamma rays:ISM -cosmic rays

Introduction

The component below about one GeV per nucleon of the interstellar cosmic-ray spectrum is probably the least well known. Direct observations are not conclusive because cosmic rays with these energies are effectively deflected by the action of the outstreaming solar wind such that the energy spectra measured inside the heliosphere are strongly altered with respect to the interstellar one. It results in an uncertainty of the particle flux which is increasing towards lower energies and below a few hundred MeV per nucleon the cosmic-ray composition and spectrum are essentially unknown. Indirect evidence for an important low-energy component comes from observations of the interstellar ionization rate. In particular, observations of the molecule H + 3 in diffuse interstellar clouds point to a cosmic-ray induced ionization rate which exceeds by more than an order of magnitude the one which is calculated for the standard cosmic-ray flux thought to be produced by diffuse shock acceleration in supernova remnants (Indriolo, Fields, & McCall (2009); Indriolo & McCall (2012)). However, neither the nature -electrons or nuclei -nor the spectrum of the low-energy component can be deduced from these observations. An independent argument for the existence of a significant component of low-energy cosmic-ray nuclei comes from the observed quasi-linear increase of Be with metallicity (Tatischeff & Kiener (2011)).

The most compelling indirect observation of this component would be a detection of nuclear γ-ray lines from the galactic disk. In fact, the strongest lines from the most abundant nuclei are preferentially produced in the cosmic-ray energy range below a few hundred MeV in collisions of protons and α-particles with interstellar matter. They are expected to be the same as frequently observed from strong solar flares, i.e. lines from the deexcitation of the first few levels in 12 C, 16 O, 20 Ne, 24 Mg, 28 Si and 56 Fe (Ramaty, Kozlovsky, & Lingenfelter (1979)). Indriolo, Fields, & McCall (2009) constructed cosmic-ray spectra which produce the observed ionization rate in diffuse clouds by adding a low-energy component to the standard cosmic-ray spectrum. With that they estimated fluxes of the 4.4-MeV line of 12 C and the 6.1-MeV line of 16 O from the central radian of the Galaxy close to the sensitivity limits of INTEGRAL. Besides strong narrow lines the nuclear line emission is composed of broad lines produced by the cosmic-ray heavy-ion component and of thousands of weaker lines which together form a quasi-continuum in the E γ ∼ 0.1 -10 MeV range.

In this paper, we calculate this total nuclear γ-ray line emission from the inner Galaxy by making use of the work of Indriolo, Fields, & McCall (2009), new developments in nuclear reaction cross sections (Murphy et al. (2009), Benhabiles et al. (2011)) and recent observations of cosmicray induced high-energy γ-ray emission with Fermi-LAT (Abdo et al. (2009); [START_REF] Strong | Proc. of the 12th ICATPP Conference on Cosmic Rays for Particle and Astroparticle Physics[END_REF]; Ackermann et al. (2012)). Finally, the possibilities for a detection of this emission with existing and eventual future space-borne γ-ray instruments are discussed.

Galactic cosmic-ray spectra and composition

Following the method of Indriolo, Fields, & McCall (2009), the interstellar cosmic-ray (CR) spectra for the calculation of the nuclear γ-ray line emission have been divided into two distinct components:

(1) Standard CRs whose origin and acceleration are widely attributed to massive stars through the action of their stellar winds followed by their explosions as supernova. Their composition and energy spectra are relatively well known above several GeV per nucleon from direct measurements while the extrapolation to lower energies can be done with propagation and demodulation calculations. An important indirect constraint on the CR spectra above several hundred MeV per nucleon is the diffuse interstellar high-energy γ-ray emission due to π 0 production and decay. This emission has been observed from several Galactic sources by Fermi-LAT, as for example from the local gas and dust (Abdo et al. (2009)) and from the inner Galaxy [START_REF] Strong | Proc. of the 12th ICATPP Conference on Cosmic Rays for Particle and Astroparticle Physics[END_REF]) which are of particular interest for the present study. These standard CRs produce an ionization rate well below the mean ionization rate deduced from H + 3 column densities observed in diffuse molecular clouds.

(2) A distinct low-energy component added to the standard CRs to reach the observed ionization rate. For this we adopt the values deduced from recent observations of the molecule H + 3 in diffuse interstellar clouds pointing to a mean ionization rate of molecular hydrogen of ζ 2 = 3-4×10 -16 s -1 [START_REF] Mccall | [END_REF]; Indriolo, Fields, & McCall (2009); Indriolo & McCall (2012)). Another requirement to this component is that the total interstellar CR spectra stay compatible with the observed CR spectra above 1 GeV per nucleon ensuring also that the high-energy γ-ray emissions stay compatible with the above cited Fermi-LAT observations.

Standard CRs

In the most important energy band for comparison with Fermi-LAT data, between about 1 GeV and a few hundred GeV per nucleon, data from some ten balloon-borne and satellite-borne instruments exist. We adopted the most recent observations of CR fluxes from the PAMELA detector onboard the Russian Resurs-DK1 spacecraft which provides high-precision data for protons and helium nuclei in the rigidity range R ∼ 1-1000 GV (Adriani et al. (2011)). Their proton fluxes agree in the kinetic energy range above E ∼ 30 GeV, where solar modulation effects should be negligible, with other recent data typically at the 20% level or better. The agreement is slightly worse for the helium fluxes, where differences up to 40% can be observed (see e.g. Fig. 1 in Adriani et al. (2011)).

To obtain the local interstellar (LIS) proton and helium fluxes the observed fluxes have been demodulated with the force-field model of Gleeson & Axford (1968) with a solar modulation parameter of Φ = 450 MV (Adriani et al. (2011)). This yields LIS proton fluxes in the range E = 0.89-1010 GeV and LIS helium fluxes in the range E = 0.36-505 GeV nucleon -1 . For the LIS helium fluxes we supposed a pure α-particle composition, ignoring an eventual small contribution of 3 He nuclei.

The extrapolation of the proton fluxes to lower energies was done with a standard leaky-box model for Galactic CR propagation. An unbroken power law in particle momentum as it results from diffuse shock acceleration (Jones & Ellison (1991); Dermer (2012)) has been used for the source spectrum: F(p) ∝ p -s . Using dF/dE = dF/dp dp/dE = dF/dp 1/v, the corresponding energy spectrum is F(E) = F 0 p -s /β, where v, p mean the particle velocity and momentum and β = v/c. F 0 and the spectral index s were then adjusted such that the propagated spectrum reproduces the LIS proton fluxes obtained from the PAMELA measurements. In the propagation calculation the escape length distribution was taken from the disk-halo diffusion model of Jones et al. (2001) and the proton inelastic cross sections from Moskalenko et al. (2002). The best adjustment was achieved with a source index s = 2.35. Figure 1 (left panel) shows the propagated proton spectrum together with the Pamela data and some other recent measurements (Menn et al. (2000) [IMAX]; [START_REF] Shikaze | [END_REF] [Bess]; Alcaraz et al. (2000a,b) [AMS]; Boezio et al. (1999) [Caprice]).

For the LIS helium spectrum we took the demodulated helium data of PAMELA from E = 0.36 -300 GeV nucleon -1 only, the two data points above that energy suffering a large uncertainty. Extrapolation of the LIS helium spectra to lower energies was done as in the proton case with the leaky-box model with the same parameters for the escape lengths. The inelastic cross section for α + 4 He was taken from the model of Tripathi, Cucinotta, & Wilson (1999) except below E α = 15 MeV where a transition to the experimental values around 10 MeV was operated. Here the best adjustment of the Pamela He data in the range E = 1-100 GeV nucleon -1 was found with a source index s = 2.32. As shown in Fig. 1, the propagated helium spectrum misses the relatively steep fall of the LIS spectrum below about 0.4 GeV per nucleon. This part is certainly subject to large demodulation uncertainties, the force-field model being only an approximation and may especially fail at these low energies. In the case of electrons, Strong et al. (2011) showed recently through constraints on the interstellar electron spectrum from synchrotron radiation that modulation may be significantly overestimated for them and that their true interstellar spectrum shows a sharp turn over below ∼1 GeV. However, as this part has only a minor influence on high-energy γ-ray emission calculations for standard CRs we simply opted for the values of the propagated spectrum. Above 300 GeV per nucleon we used the data of the ATIC-2 experiment (Panov et al. (2009)) which have smaller uncertainties than the Pamela data there and agree with them at lower energies. The spectra of heavier nuclei were taken identical to the helium spectrum with relative CR abundances given in table 9.1 of [START_REF] Longair | High Energy Astrophysics: Volume 1, Particles, photons and their detection[END_REF], averaged over the three energy bands.

With these CR spectra the γ-ray emissivity due to π 0 decay has been calculated with the model of Kamae et al. (2006) for the differential γ-ray emission cross section in the π 0 → 2γ decay. Cross section data for π 0 production in p + p reactions have been taken from the parametrization of Dermer (1986). For the p + 4 He reaction the inclusive cross section displayed in Fig. 4 of Murphy, Dermer, & Ramaty (1987) up to 2 GeV has been used. Above 6.3 GeV the cross section has been obtained by multiplying the p + p cross section with the multiplication factor m Hep of Mori (2009) and a linear interpolation was applied in the energy gap. All other reaction channels (α + 4 He and p,α + HI) were calculated with Mori's multiplication factors relative to the p + 4 He cross section. Here, HI stands for elements C, N, O, Ne, Mg, Al, Si, S, Ar, Ca, Fe, and Ni. For the ambient local interstellar composition, solar photospheric abundances have been used. In the CR composition we added 11 B and 52 Cr and 55 Mn as representatives of Sc-Mn, which are relatively abundant in cosmic rays. The result of the calculation is shown in Fig. 2 together with the Fermi-LAT data for the local diffuse γ-ray emissivity of Abdo et al. (2009) and their estimate of the bremsstrahlung component. The calculated curve reproduces remarkably well the observed emissivity in shape and absolute normalization, very close to the calculation of Abdo et al. (2009) albeit with slightly different CR spectra and composition.

The Galactic CR electron spectrum is constrained by its electromagnetic emissions induced during propagation in the Galaxy and by direct data from many balloon and satellite experiments (see e.g. Strong, Moskalenko & Reimer (2009)). The emissions span a wide range from the Radio band for the synchroton radiation to MeV and TeV energies for bremsstrahlung and inverse Compton scattering. For the standard CR electron spectrum we rely on results of the GALPROP model, which has been developed to calculate numerically the Galactic propagation of hadrons and leptons including secondary-particle production and propagation (Strong & Moskalenko (1998)). The model makes also detailed calculations of the electromagnetic emissions which are induced by primary and secondary particles. From the different CR electron spectra and propagation models studied in [START_REF] Strong | [END_REF], we adopt the conventional model. It fits the direct data above several GeV and agrees simultaneously with the Galactic synchrotron emission. The latter constraint requires a spectral break of the CR electron spectrum below several GeV, which was confirmed in a recent study by Strong et al. (2011).

Low-energy part

The above described standard hadronic CRs yield an ionization rate for molecular hydrogen of ζ 2 = 4.3×10 -17 s -1 integrated in the CR energy range E = 0.002 -10 3 GeV nucleon -1 , about a factor of ten below the most recent observed mean rates in diffuse molecular clouds ζ 2 = 4×10 -16 s -1 (Indriolo, Fields, & McCall (2009) ) andζ 2 = 3.5 +5.3 -3.0 ×10 -16 s -1 (Indriolo & McCall (2012)). For a direct comparison with Indriolo, Fields, & McCall (2009) who first proposed a distinct lowenergy CR (LECR) component above 2 MeV per nucleon to account for a high ζ 2 , we adopted in the following for calculations and figures the first value.

Using the conventional model electron spectrum of [START_REF] Strong | [END_REF] extrapolated down into the eV range, Padovani, Galli & Glassgold (2011) found a negligible contribution of electrons to the ionization rate in molecular clouds (ζ 2 ∼10 -19 s -1 ). There could still be an important CR electron component below 200 MeV to produce the observed ionization rate. Such a LECR electron component has for example been studied by [START_REF] Strong | [END_REF] and Porter et al. (2008) in order to reproduce the Comptel MeV-excess. But this is out of the scope of this paper and we will concentrate in the following on the hadronic component. Supposing that the extra ionization in these clouds is entirely due to low-energy nuclei we explored a range of spectra and compositions for such a LECR component. Scherer et al. (2008a) proposed that a CR population of anomalous cosmic rays (ACRs) accelerated in astrospheres of stars similar to the Sun makes up a significant part of the interstellar CR spectrum below E = 0.3 GeV, and when adding the contributions of other stars those ACRs may be the dominant component of LECRs. Such a scenario was discussed by Indriolo, Fields, & McCall (2009) for their low-energy power-law spectrum (dubbed "carrot" spectrum) with powerlaw index and normalization adjusted to produce the observed ionization rate. We adopted this model in the first scenario for our low-energy part of the interstellar spectrum.

The spectra and composition were derived in the following way: LECR source spectra were constructed from the ACR proton spectra at the helio-and astropause of Scherer et al. (2008a), which were constrained by recent data of the Voyager spacecraft. These spectra are well described by a power-law function with exponential cutoff at E c : F(E) = F 0 E -s e -E/Ec with a spectral index of s = 2.4. We explored ACR spectra for this index and an upper and lower limit of s = 2.7 and s = 2.0, respectively. For α particles and heavier nuclei the same index was used and the cutoff energy was taken following the species scaling of [START_REF] Steenberg | AIP Conf. Proc. 528, Acceleration and Transport of Energetic Particles Observed in the Heliopshere[END_REF]:

E c (A,Q)/E c (p) = (A/Q) -2γ/(γ+1)
with energies expressed in energy per nucleon and A, Q being the mass number and the charge state, respectively. For sake of simplicity we used γ = 1 and considered only singly-charged ions, which represent the majority of ACR ions in the heliosphere, reducing the scaling to E c (A,Z)/E c (p) = A -1 . E c was varied from 5 MeV to an upper limit, where the particle fluxes are still compatible with the observed CR helium and proton fluxes above 0.5 and 1 GeV per nucleon, respectively. The upper limits were found to be E c = 600, 1200 and 4800 MeV for spectral indices s = 2.0, 2.4 and 2.7, respectively. The ACR composition was taken from Cummings, Stone, & Steenberg (2002).

For the second scenario we used source spectra resulting from shock acceleration (SA-LECRs) as for the standard CRs, but with an energy cutoff E c : F(E) = F 0 R -s /β e -E/Ec . R is the particle rigidity; spectral index s and E c , expressed in energy per nucleon, were taken the same for all species. For the source spectral index we used s = 2.35, which is the index found for the standard CRs as our central value, and s = 2.0 and s = 2.7 falling in the typical range of models of diffuse shock-acceleration (see e.g. Jones & Ellison (1991); [START_REF] Schlickeiser | Cosmic Ray Astrophysics[END_REF]). The energy cutoff E c , which in this scenario may be considered arising from a tuncation in space or time in the acceleration process was varied from 5 MeV to an upper limit compatible with the observed CR helium and proton fluxes above 0.5 and 1 GeV per nucleon, respectively. It was found to be E c = 150 MeV for s = 2.7 and E c = 120 MeV for s = 2.35 and s = 2.0. For the composition of the source spectra we took the CR source abundances of [START_REF] Lund | Cosmic Abundances of Matter, A.I.P. Conf. Proc[END_REF] from table 9.1 in Reames (1999).

These spectra were then propagated for both scenarios with the leaky-box model to obtain the interstellar LECR spectra. All propagated LECR spectra were normalized to produce an ionization rate ζ 2 = 3.6×10 -16 s -1 in diffuse clouds, yielding ζ 2 = 4×10 -16 s -1 when added to the standard CRs. We calculated the ionization rate as Indriolo, Fields, & McCall (2009), assuming in particular that CRs of kinetic energy below a threshold energy of 2 MeV per nucleon cannot penetrate diffuse clouds and thus do not affect the ionization rate there.

Examples of LECR proton and α-particle spectra are shown in Fig. 3 together with the standard CRs alone. The curves with s = 2.35, E c = 120 MeV and s = 2.4, E c = 1.2 GeV for SA-LECRs and ACR-LECRs, respectively, represent upper limits still compatible with the demodulated proton and helium data. The ionization rates of LECRs added to the standard CRs in dense clouds, taking a low-energy threshold of 10 MeV as proposed in Indriolo, Fields, & McCall (2009), are displayed in Fig. 4. Nearly all spectra yield rates above the recommended value of van der Tak & van Dishoeck (2000) for molecular cloud cores and with the exception of ACR-LECRs with s = 2.7 cross their upper limit for the higher cutoff energies. We note, however, that the energy threshold of 10 MeV is an approximation based on stopping range calculations and certainly subject to considerable uncertainties. Furthermore, magnetic mirroring can decrease the CR ionization rate in molecular cloud cores by a factor of ∼3-4 (Padovani & Galli (2011)), which would bring all rates if not into the recommended range but at least below the upper limit discussed in van der Tak & van Dishoeck (2000).

Nuclear γ-ray line emission of the inner Galaxy

Nuclear line emission was calculated for the same reaction channels as for the calculation of π 0 production except the p + p channel which produces no γ-ray lines. Cross sections have been taken from the recent compilation of Murphy et al. (2009), which contains cross section excitation functions of about 140 relatively strongly produced γ-ray lines by energetic proton and α-particle reactions in astrophysical phenomena like solar flares. The excitation functions extend from reaction threshold to typically 1 GeV per nucleon which in most cases is amply sufficient for the present study where the nuclear γ-ray line emission is dominated by the low-energy part. Above that range, constant cross sections were supposed.

For all other γ-ray lines we relied entirely on calculations with the reaction code TALYS (Koning, Hilaire, & Duivestijn (2008)) as it was done in Murphy et al. (2009), with some slight parameter changes for 14 N, 20,22 Ne and 28 Si which gave an improved reproduction of new experimental γ-ray data for proton and α-particle induced reactions with these nuclei (Benhabiles et al. (2011)). This code works in the range up to 250 MeV for proton reactions and 62.5 MeV per nucleon for α-particle reactions. Above these energies a constant cross section was supposed. Lines emitted during the decay of radioactive nuclei were included in these calculations but not the 511-keV line and continuum emissions resulting from β + decay and subsequent positron annihilation. Also emission induced by π + -decay positrons was not included in the calculated spectra.

As described above, the cosmic-ray spectra were normalized such that an ionization rate ζ 2 = 4×10 -16 s -1 was obtained and that the calculated high-energy emission due to π 0 decay reproduced the Fermi observations of the diffuse γ-ray emission of the inner Galaxy [START_REF] Strong | Proc. of the 12th ICATPP Conference on Cosmic Rays for Particle and Astroparticle Physics[END_REF]). For the ambient medium composition solar abundances with two times solar metallicity were taken. The Fermi-LAT data and calculations for ACR-LECR spectra with s = 2.4 and two extreme cases of energy cutoffs E c = 5 MeV and E c = 1.2 GeV are shown in Fig. 5. When adding the inverse Compton component from high-energy leptons and the estimated emission due to unresolved point sources as presented in [START_REF] Strong | Proc. of the 12th ICATPP Conference on Cosmic Rays for Particle and Astroparticle Physics[END_REF] the reproduction is quite good for both cases. The π 0 production is apparently completely dominated by the standard CRs, even for LECRs with E c = 1.2 GeV which add only a small component below E ∼ 1 GeV.

Examples of nuclear γ-ray line emission for three cases of ACR-LECRs with s = 2.4 and the case with strongest emission s = 2.0, E c = 120 MeV of SA-LECRs are shown in Fig. 6. Emission due to the standard CRs alone, as described in Sect. 2.1, is also shown for comparison. The total γ-ray line intensity for the latter is dominated by interactions of the heavy-ion content due to their enrichment in standard CRs with respect to the interstellar medium. It accounts for about 70% of the total nuclear line emission. Some outstanding narrow lines, produced in proton-and α-particle induced reactions are nevertheless present. The emissions of ACR-LECRs with cutoff energies E c = 5 MeV and 25 MeV add practically only a narrow-line component, the quasi-continuum being completely dominated by the standard CR component of the CR spectra. In fact, this is mainly because heavy ions play a completely negligible role in low-energy cutoff ACR-LECRs due to their very low cutoff energies (E c (A,Z)/E c (p) = A -1 ; see Sect. 2.2). Their contribution is less than 10 -3 at E c ≤ 8 MeV, ∼6% at E c = 80 MeV and raises to ∼50% at E c = 1.2 GeV.

The overall emission is stronger for SA-LECRs because more particles above the reaction thresholds are present with respect to ACR-LECRs of the same cutoff energy. ACR-LECR spectra with their steep rise towards low energies contain many low-energy particles that contribute strongly to ionization but are very inefficient for γ-ray line production. This is true for the narrow lines induced by light ions and slightly more pronounced in the quasi-continuum, which reflects furthermore the importance of energetic heavy ions in SA-LECRs. They contribute between ∼50% at E c = 5 MeV and ∼60% at E c = 120 MeV to the nuclear γ-ray emission.

The strongest narrow lines in decreasing order are the 4.4-MeV line (essentially from 12 C and 11 B, produced in reactions with 12 C and 16 O), the 6.1-MeV line (mainly inelastic scattering off 16 O and spallation reaction products 15 O and 15 N) and the 1.63-MeV line (mainly from inelastic scattering reactions with 20 Ne). Integrated fluxes of these narrow lines and the integrated fluxes in energy ranges 1-3 MeV and 3-8 MeV are shown on Fig. 7. Although the continuum fluxes increase steadily with cutoff energy, the narrow-line fluxes tend to level or even slightly decline above E c ∼50 MeV. The latter can be explained by the cross section excitation functions of the dominant line-producing reactions which have their maxima generally well below 50 MeV. An exception are 6.13-MeV fluxes for SA-LECRs which include an important contribution of 16 O spallation reactions with higher thresholds.

Narrow-line and continuum fluxes have also been calculated with two spectra and compositions as proposed by Indriolo, Fields, & McCall (2009). The flux values for their "carrot" spectrum are within 20% of the values for ACR-LECRs with s = 2.7, E c = 8 MeV while the values for their power-law spectrum agree to better than 15% with the ones for SA-LECRs with s = 2.0, E c = 14 MeV. The 4.4-MeV and 6.1-MeV narrow line fluxes are smaller by about a factor of ten than the fluxes predicted by Indriolo, Fields, & McCall (2009) for the central radian of the Galaxy (330 • < l < 30 • , |b| < 10 • ) probably mainly because of the method of absolute normalization. In Indriolo, Fields, & McCall (2009) it was based on the column density of the target material in the central radian assumed to be 10 23 cm -2 . By normalizing the nuclear line emission to the π 0 -decay flux from the inner Galaxy as determined by Fermi-LAT [START_REF] Strong | Proc. of the 12th ICATPP Conference on Cosmic Rays for Particle and Astroparticle Physics[END_REF]) our 4.4-MeV and 6.1-MeV line fluxes are independent of the column density. They depend only on the spectral shape and composition of CRs and the metallicity of the inner Galaxy. Supposing that the intensity of the CRs which are responsible for the bulk of the π 0 production (E ∼0.5 -20 GeV per nucleon) is the same throughout the Galaxy, the Fermi-LAT data imply a mean column density of about 2×10 22 cm -2 in the inner Galaxy (300 • < l < 60 • , |b| < 10 • ). This explains at least a large part of the difference.

Discussion

Our estimations of nuclear cross section uncertainties can roughly be divided into two categories. (1) Strong lines whose cross sections excitation functions are explicitely listed in the compilation of Murphy et al. (2009). They are at least partly based on experimental data whose uncertainties are of the order of 15% or less. For the 1.63, 4.4 and 6.1-MeV lines, the bulk of the emission is induced by particles in the energy range where a wealth of experimental data exist, and the nuclear uncertainties on the narrow-line fluxes should be less than 20%. (2) The weak-line quasi-continuum which is taken from TALYS (Koning, Hilaire, & Duivestijn (2008)) calculations. The predictions of this code have been extensively compared to experimental γ-ray cross sections in several studies ( Benhabiles et al. (2011) Murphy et al. (2009) Tatischeff et al. (2006)) and a typical agreement of a factor of two or better for individual lines was found. As this quasi-continuum consists of thousands of individual lines, the overall uncertainty on the integrated cross section in the two energy bands 1-3 and 3-8 MeV should be smaller than that. A conservative estimate would be less than 50%. This weak-line emission component is maximal for the high-energy CR component (∼70%), while for the LECRs it represents ∼50% at the highest-energy cutoffs and at E c = 5 MeV drops to ∼20% and ∼25% for ACR-LECRs and SA-LECRs, respectively. It results in nuclear uncertainties of the 1-3 and 3-8 MeV fluxes of less than ∼40%.

Another uncertainty of the narrow-line fluxes comes from the metallicity of the ambient medium. Measurements generally agree for the average metallicity at the Galactocentric distance of the Sun to be equal or slightly below the solar system value (Daflon & Cunha (2004), [START_REF] Lemasle | [END_REF], Pedicelli et al. (2009), Rolleston et al. (2000)). The radial gradients however, especially towards the Galactic center (R = 0), appear less well defined. Using the available extracted radial gradients, and assuming they stay constant towards the center, extrapolation of element abundances to the center yield values ranging from [O/H] = 0 (Daflon & Cunha (2004)) to [Fe/H] = 1 (Pedicelli et al. (2009)), where [O/H]=log10((O/H) R=0 /(O/H) ⊙ ), and (O/H) ⊙ , the solar system value. The fact that the inner Galaxy largely dominates the γ-ray emission and positive metallicity gradients, however, make a metallicity higher than solar reasonable in this case. Our choice of taking [X/H] = 0.3, X = C-Ni for the average metallicity towards the inner Galaxy is only a guess with a conservative estimated uncertainty of a factor of two. This uncertainty affects, however, less the quasi-continuum fluxes which are to a great part due to the CR heavy-ion component. Adding that to the nuclear uncertainties a typical overall uncertainty of a factor of two is expected for the calculated narrow-line fluxes and between 50% and a factor of two for the continuum fluxes.

The explored parameter space for LECR spectra and composition spans a wide range in spectral forms and composition, from extremely metal-poor for ACRs with low E c to the metal-rich propagated CR source abundances for SA-LECRs. All spectra are compatible with the high-energy diffuse emission of the inner Galaxy observed by Fermi-LAT, the ionization rate of diffuse clouds, and recent satellite and balloon data of the CR proton and helium spectra. Estimated ionization rates of dense clouds, however, are for all CR spectra, except for the low cutoff-energy ACR-LECRs somewhat higher than the recommended value. These estimations are very uncertain, however, as they neglect the effects of magnetic fields on the CR propagation. Padovani & Galli (2011) recently found that magnetic mirroring always dominates over magnetic focusing, implying a reduction of the CR ionization rate in molecular cloud cores by a factor of ∼3-4.

It would of course be very surprising if the LECR component could be described by a single distribution with one definite spectral shape, index, cutoff energy and composition, and be exclusively composed of nuclei. However, whatever the nature of LECRs, they can probably be approximated by a combination of the spectra and compositions studied here and we believe that the nuclear γ-ray line emission falls into the range which we have explored and which is presented in Figs. 6,7. It is therefore interesting to compare our predictions with γ-ray instrument sensitivities.

Integrating over the solid angle of the inner Galaxy (300 • < l < 60 • , |b| < 10 • ) yields flux values (in cm -2 s -1 ) 0.73 times the values presented in Fig. 7. These values are unfortunately far below the narrow-line sensitivities of SPI/INTEGRAL for diffuse emissions: ∼5 × 10 -4 cm -2 s -1 for the 4.4-and 6.1-MeV lines with FWHM widths ∆E ∼100 keV and ∼2 × 10 -4 cm -2 s -1 for the 1.63-MeV line with ∆E ∼20 keV (Teegarden & Watanabe (2006)). An important fraction of the 6.1-MeV line may be very narrow, if much of the interstellar oxygen is locked up in dust grains (Tatischeff & Kiener (2004)). It would result in a sensitivity gain of a factor of 5 at maximum, but even that is still far from the highest estimated fluxes. In summary, no constraints on LECRs can be inferred from SPI/INTEGRAL observations of the inner Galaxy.

A definite detection of CR-induced nuclear γ-ray line emission or stringent constraints on the nature of LECRs in the case of non-detection has probably to wait for the next-generation telescopes with an expected increase in sensitivity of typically a factor of 10-100. We therefore compare the calculated fluxes to instrument proposals to the last call (2010) of ESA's 2015-2025 Cosmic Vision program. Of the three instruments with similar estimated sensitivities in the 0.1-10 MeV range which have been presented we focus on the CAPSiTT proposal [START_REF] Lebrun | Proceedings of the 8th INTEGRAL Workshop "The Restless Gamma-ray Universe[END_REF]) for which some of us participated substantially in the performance estimates.

The most promising narrow line comparing CAPSiTT sensitivities with the predicted fluxes is the 4.4-MeV line, the sensitivity gain for an eventual very-narrow-line component of the 6.1-MeV line is less significant for this instrument with an estimated energy resolution at 6 MeV of ∼30 keV. Assuming a uniform emission in the defined region, the CAPSiTT 3σ 5-y survey line sensitivity for the 4.4-MeV line (∆E = 100 keV FWHM) is ∼7 × 10 -6 cm -2 s -1 for such an extended emission, using only Compton events in the CAPSiTT instrument. This is at the limit for ACR-LECR spectra and would allow detection for SA-LECRs with cutoff energies above 15 MeV. Pair-creation events, which have not been considered for the proposal to ESA below E γ = 10 MeV can improve significantly the sensitivity for such extended emissions, increasing e.g. the effective detection area threefold at E = 5 MeV. This would allow a detection for nearly all considered LECR scenarios.

The broadband 1-3 and 3-8 MeV sensitivities in 5 years amount to ∼2 × 10 -5 cm -2 s -1 and ∼1 × 10 -5 cm -2 s -1 , respectively. The flux expected from the standard CRs alone is already slightly higher than the sensitivity in the 1-3 MeV band and well above the sensitivity in the 3-8 MeV band. A detection looks therefore very promising for all LECR scenarios. It may, however, be difficult to disentangle from other emissions like inverse Compton scattering of CR electrons and unidentified point sources, which probably dominate the diffuse flux of the inner Galaxy in the 1-10 MeV range. Analysis and discussion of the diffuse interstellar emission as measured by SPI/INTEGRAL and by COMPTEL/CGRO in this energy range can be found in Bouchet et al. (2011).

The nuclear γ-ray line emission has probably a latitude and longitude profile similar to that of the high-energy diffuse emission observed by Fermi-LAT [START_REF] Strong | Proc. of the 12th ICATPP Conference on Cosmic Rays for Particle and Astroparticle Physics[END_REF]). This emission is approximately constant along 320 • < l < 40 • and concentrated on the Galactic plane, where about 60% of the emission is in a band with b = ±1.5 • . Differences may arise because of the metallicity dependence of narrow-line fluxes, which is negligible for π 0 -decay emission, essentially produced in reactions with ambient H and He. Differences of ionization rates in diffuse interstellar clouds, sometimes observed for sight lines with small angular separation indicate a non-uniform localized distribution of LECRs (Indriolo & McCall (2012)). This is not surprising since low-energy particles have shorter ranges than the standard CRs resulting certainly also in an inhomogeneous distribution of LECRs in the inner Galaxy, concentrated around the acceleration sites. This would probably affect only slightly the integrated fluxes since they are tied to the mean ionization rate, but the spatial distribution of the nuclear line emission would then be different from the high-energy emission from π 0 decay. These effects would be favorable for an eventual observation by a Compton-imaging telescope like CAPSiTT, whose estimated angular resolution at 5 MeV (0.8 • ) is about a factor of three better than that of COMPTEL/CGRO and SPI/INTEGRAL. The point-source sensitivities of CAPSiTT in 5 years amount to ∼2 × 10 -7 cm -2 s -1 and ∼3 × 10 -7 cm -2 s -1 for the 4.4-MeV line and the 3-8 MeV band, respectively. Our calculated fluxes of the 4.4-MeV line from the inner Galaxy are typically 10-10 2 times higher than that. If e.g. the spatial extension of the inhomogeneities is of the order of the angular resolution, a detection gets very probable. Assuming that the large scatter of observed ionization rates in diffuse clouds (see e.g. Indriolo & McCall (2012)) reflects the LECR density, the 4.4-MeV flux variations may reach a factor of 10. This would be very favorable for a detection by a futur instrument with the characteristics of CAPSiTT. Our calculated fluxes in the 3-8 MeV band are a factor of 10 2 -10 3 higher than the point-source CAPSiTT sensitivity, and a concentration of a few percent of this flux in a point source would be detectable. Following the scatter of ionization rates, a factor-of-few variations in the fluxes are probable although a non-negligible part of the 3-8 MeV band is produced by standard CRs. That inevitably would be detected by a futur instrument with good angular resolution.

Conclusion

The presented nuclear line and continuum fluxes from the inner Galaxy induced by interactions of CRs in the interstellar medium were calculated for a large variety of LECR spectra and composition and it is probable that the real values fall somewhere inside the obtained ranges. The total nuclear γ-ray emission from the inner Galaxy is predicted to be in the range (0.1 -2.0) × 10 -5 cm -2 s -1 for the 4.4-MeV line, (0.1 -1.0) × 10 -5 cm -2 s -1 for the 6.1-MeV line and (0.3 -3.7) × 10 -6 cm -2 s -1 for the 1.63-MeV line. For the broad-band emissions we find (0.2 -1.3) and (0.3 -2.1) × 10 -4 cm -2 s -1 for the 1-3 and 3-8 MeV bands, respectively. An eventual observation of narrow nuclear γ-ray lines with future γ-ray space telescopes based on actually available technology appears to be in reach for a large range of LECR scenarios, although challenging for scenarios with low energy cutoffs, where fluxes do not exceed a few 10 -6 cm -2 s -1 from the inner Galaxy. It may also be promising for the continuum, especially the 3-8 MeV continuum flux which

Fig. 2 .

 2 Fig. 1.-Left: LIS proton fluxes from demodulated data of five recent experiments and propagated proton spectrum from the leaky-box model. The modulation potential for the force-field approximation used to demodulate the different data sets are indicated in the left figure. Right: LIS helium fluxes from the same experiments and from ATIC-2. No demodulation was applied to the ATIC-2 data because they start above 10 GeV per nucleon, where modulation should be negligible.

Fig. 5 .

 5 Fig. 5.-Fermi-LAT data of the diffuse γ-ray emission of the inner Galaxy (300 • < l < 60 • , | b | < 10 • ) (Strong (2011)) are shown by the black triangles. Solid lines: Calculated γ-ray emissions from π 0 decay with CR spectra including a LECR component following the model of Scherer et al. (2008a) for ACRs with s = 2.4, E c = 1200 MeV (blue) and s = 2.4, E c = 8 MeV (red). The curve for standard CRs is indistinguishable from the red curve, since for E c = 8 MeV the π 0 production is entirely due to the standard CRs. The dashed lines show the total emission with contributions of inverse Compton scattering and unresolved point sources as estimated in Strong (2011) added to the calculated π 0 -decay emission.

exceeds the estimated future instrument sensitivities by up to a factor of 20. This is especially true for an instrument with good angular resolution which allows the separation of this emission from other diffuse emissions and a very inhomogeneous spatial distribution of LECRs in the inner Galaxy could not be missed by such an instrument. ). The dotted line represents their upper limit (∼1.3×10 -16 s -1 ).


Parameter estimation from experimental data is a crucial problem in quantitative modeling of biochemical reaction networks. An especially important issue, raised by the complexity of the models and the challenging nature of the experimental data, is parameter identifiability. Despite several approaches proposed in the systems biology literature, no agreement exists on the analysis of structural and practical identifiability, and the relations among the two. In this paper we propose a mathematical framework for the analysis of identifiability of metabolic network models, establish basic results and methods for the structural and practical identifiability analysis of the class of so-called linlog models, and discuss the results on the basis of an artificial example.

Structural and practical identifiability of approximate metabolic network models

Sara Berthoumieux * Daniel Kahn * * Hidde de Jong * Eugenio Cinquemani * ,⋆

INTRODUCTION

Kinetic models of biochemical reaction systems usually contain a large number of parameters, many of which are difficult to measure in a direct way. This makes parameter estimation from experimental data a crucial problem for quantitative systems biology [START_REF] Ashyraliyev | Systems biology: Parameter estimation for biochemical models[END_REF][START_REF] Crampin | System identification challenges from systems biology[END_REF]. For all but the simplest systems, parameter estimation is a difficult problem. Models contain a large number of variables, whose dynamics evolve on different time-scales and are described by complex, nonlinear rate equations. Moreover, the available experimental data usually consist of noisy and incomplete measurements, obtained under different conditions and by means of heterogeneous experimental methods.

These difficulties have stimulated the use of approximate kinetic models, with rate equations following, for example, linear, piecewise-linear, linlog, or power-law kinetics [ [START_REF] Jong | Modeling and simulation of genetic regulatory systems: A literature review[END_REF][START_REF] Heijnen | Approximative kinetic formats used in metabolic network modeling[END_REF]. In general, these models have less parameters to estimate and, in many cases, parameter estimation can be reduced to linear or orthogonal regression [START_REF] Nikerel | Model reduction and a priori kinetic parameter identifiability analysis using metabolome time series for metabolic reaction networks with linlog kinetics[END_REF]. Therefore, approximate kinetic models have been chosen and successfully used for the quantitative modeling of both metabolic networks and gene regulatory networks [START_REF] Heijnen | Approximative kinetic formats used in metabolic network modeling[END_REF].

We are interested in how approximate kinetic models help addressing the key problem of the identifiability of biochemical networks in a principled and scalable way. We focus on the case where the structure of the model is fixed by a priori knowledge on the network (i.e. the chemical species considered and the reactions among them), and discuss the identifiability of the model parameters. That is, we are interested in the problem of unambiguously reconstructing the unknown parameter values from the observed network behavior.

⋆ Corresponding author, e-mail: eugenio.cinquemani@inria.fr

A distinction is usually made between structural (or a priori) and practical (or a posteriori) identifiability [Ljung, 1999, Walter and[START_REF] Walter | Identification of parametric models[END_REF]. Structural identifiability is an intrinsic property of the model family, guaranteeing that unique parameter reconstruction would be possible from perfect observations of the system response to an arbitrarily rich set of inputs. Practical identifiability refers to the ability of estimating unknown parameter values from the available experimental data within a prespecified degree of accuracy. In classical control theory, this concept is essentially related to the notion of persistence of excitation (see e.g. [START_REF] Ljung | System identification, theory for the user[END_REF] for the case of linear dynamical systems). However, in a biological context, this notion requires at least some nontrivial adaptation, and no agreement on the concepts of structural and practical identifiability exists.

The aim of this paper is to discuss identifiability of kinetic models of metabolism from a rigorous theoretical standpoint. While most of our definitions are of general applicability, identifiability results will be developed primarily for approximate kinetic models known as linlog models [START_REF] Visser | Dynamic simulation and metabolic re-design of a branched pathway using linlog kinetics[END_REF], whose pseudo-linear form enables us to apply tools from linear algebra and estimation theory in a straightforward manner. Similar results can be derived easily for many other approximate kinetic model classes in pseudo-linear form, such as the linear, loglin and generalized mass-action kinetic models [START_REF] Delgado | Metabolic control analysis using transient metabolite concentrations[END_REF][START_REF] Hatzimanikatis | Effects of spatiotemporal variations on metabolic control: Approximate analysis using (log)linear kinetic models[END_REF][START_REF] Savageau | Biochemical Systems Analysis: A Study of Function and Design in Molecular Biology[END_REF].

The contributions of the paper are as follows. First of all, we provide precise definitions and links between structural and practical identifiability of kinetic network models.

Then, based on this solid mathematical foundation, our analysis provides tools that are applicable to the currently available data sets, e.g., those obtained by means of recent high-throughput methods in biology [START_REF] Ishii | Multiple high-throughput analyses monitor the response of E. coli to perturbations[END_REF].

This sets the stage for rigorous kinetic model reduction, which will be addressed separately in another paper.

The paper is organized as follows. In Section 2 we discuss kinetic models as well as linlog and related approximations. In Section 3 we discuss the notions of structural and practical identifiability, and provide readily applicable methods for the identifiability analysis of linlog models.

Results are illustrated throughout by means of a small example. Conclusions and perspectives of the work are reported in Section 4. Proofs of theoretical results are rather straightforward and are omitted in the interest of space.

PARAMETER ESTIMATION IN LINLOG AND RELATED MODELING FRAMEWORKS

The dynamics of biochemical reaction networks are described by kinetic models having the form of systems of ordinary differential equations (ODEs) [START_REF] Heinrich | The Regulation of Cellular Systems[END_REF]. In this paper we focus on kinetic models of metabolism, where the rate functions describe enzymecatalyzed reactions. This leads to models of the general form:

ẋ = N • v(x, u, e), (1) with x(0) = x 0 ∈ R nx ≥0 , where x ∈ X ⊆ R nx ≥0 denotes the vector of (nonnegative) internal metabolite concentrations, u ∈ U ⊆ R nu ≥0 the vector of external metabolite concentra- tions, e ∈ E ⊆ R m
>0 the vector of enzyme concentrations, and v : R nx+nu+m ≥0 → V , with V ⊆ R m , the vector of reaction rate functions. N ∈ Z nx×m is a stoichiometry matrix.

Kinetic modeling formalisms differ in the choice of rate functions. Examples of classical functions in enzyme kinetics are the Michaelis-Menten, reversible Michaelis-Menten, and Monod-Wyman-Changeux rate laws [START_REF] Heinrich | The Regulation of Cellular Systems[END_REF]. Approximate formalisms simplify the mathematical form of the rate laws, in particular the nonlinear dependency of the reaction rates on the metabolite concentrations. Moreover, they usually assume all reactions to follow the same simplified kinetic format, thus giving a uniform structure to the models. For one type of approximate kinetic formalism, so-called linearlogarithmic (linlog) models, below we show that the parameter estimation problem can be formulated as multiple linear regression. A brief discussion of how this can be equally done for a number of other well-known approximate formalisms is reported at the end of the section.

The linear-logarithmic (linlog) approximation [Heijnen, 2005, Visser and[START_REF] Visser | Dynamic simulation and metabolic re-design of a branched pathway using linlog kinetics[END_REF]] expresses the reaction rates as proportional to the enzyme concentrations and to a linear function of the logarithms of internal and external metabolite concentrations, v(x, u, e) = diag(e) • a + B x • ln(x) + B u • ln(u) , (2) where diag(e) is the square diagonal matrix with the elements of e on the diagonal, and the logarithm of a vector means the vector of logarithms of its elements. Concavity of the function with respect to metabolite concentrations, and the fact that B x and B u are in direct relationship with the so-called system elasticities [Heijnen, 2005, Heinrich and[START_REF] Heinrich | The Regulation of Cellular Systems[END_REF]] make linlog modelling a very convenient (local) approximation of metabolic kinetics. For con-ciseness, in the sequel we shall often drop the dependence of v on (x, u, e) from the notation.

The identification of metabolic networks in the linlog formalism amounts to estimating the (generally unknown) parameters a ∈ R m , B x ∈ R m×nx and B u ∈ R m×nu from experimental data. In most experiments, concentrations of enzymes and external metabolites are under (partial) control of the experimentalist, and the concentrations of internal metabolites and metabolic fluxes are measured after the system has relaxed to the steady state N • v(x, u, e) = 0.

(3) In accordance with this, we shall assume that, from each of q ∈ N experiments, the data are (noisy) measurements (ṽ k , xk , ũk , ẽk ) of (v k , x k , u k , e k ), where the latter satisfy v k = v(x k , u k , e k ) and (3), with k = 1, . . . , q. Clearly the restriction to steady-state measurements limits the informativity of the data and may affect the identifiability of the models, as will be apparent in later sections.

For the purpose of parameter estimation, it is convenient to rewrite (2) in the form of a regression model:

(v/e) T = 1 ln(x) T ln(u) T • [a B x B u ] T (4)
where the ratio of two vectors (here v/e) denotes elementwise division. Let us use an upperbar to denote the mean of a quantity over its q experimental outcomes, for instance: v/e = (1/q) q k=1 v k /e k . By the linearity of (4), it holds that

(v/e) T = 1 ln(x) T ln(u) T • [a B x B u ] T . (5)
This allows one to rewrite (4) as a mean-removed model

v e - v e T = ln(x) -ln(x) ln(u) -ln(u) T • (B x ) T (B u ) T . (6)
We can now define the following estimation problem. Problem 1. Given the data matrices

     ṽ1 ẽ1 -ṽ ẽ T . . . ṽq ẽq -ṽ ẽ T      W ,     ln(x 1 ) -ln(x) T ln(ũ 1 ) -ln(ũ) T . . . . . . ln(x q ) -ln(x) T ln(ũ q ) -ln(ũ) T     Ỹ find parameters B [B x B u ] T minimizing ||| W -Ỹ • B|||, for some convenient (matrix) norm ||| • |||.
We fix this to be the Frobenious norm.

Notice that the parameter vector a no longer appears in the regression problem, but that an estimate of it can be recovered from estimates of

B = [B x B u ]
T via Eq. ( 5).

In practice, the experimental measurement error on the normalized fluxes v/e is often significantly larger than that on metabolite concentrations. In light of this [START_REF] Berthoumieux | Identification of linlog models of metabolic networks from incomplete high-throughput datasets[END_REF] we consider the linear model

W = W + ε = Y • B + ε,
where W and Y are the noiseless versions of W and Ỹ , respectively. ε ∈ R q×m is assumed to be a zero-mean Gaussian matrix with independent columns, each with q-dimensional positive definite covariance matrix Var(ε

•i ) = Σ ε•i , i = 1, . . . , m.
Thus, the problem is equivalent to the m least-squares problems of minimizing

||W •i -Y • B •i || Σε •i over B •i , with i = 1, . . . , m, where || • || Σε •i is the Σ ε•i -weighted L 2 -norm (i.e., for w ∈ R q , ||w|| Σε •i = w T Σ -1 ε•i w)
, and W •i and B •i are the ith columns of W and B, respectively. That is, the parameter estimation problem splits up into m smaller estimation problems, one for each reaction i, with i = 1, . . . , m.

In turn, each reaction i depends only on a (known) subset of metabolites. In this case, the elements of B •i corresponding to the metabolites not involved in the reaction can be set to zero, and the least-squares problem can be reduced accordingly. If C(i) ⊆ {1, . . . , n x + n u } denotes the indices of the relevant values in B .i , then the regression problem becomes min

B C(i)i ||W •i -Y •C(i) • B C(i)i || Σε •i , (7)
where the meaning of B C(i)i follows from the above and

Y •C(i)
is the submatrix of Y formed by the columns indexed by C(i).

Similar parameter estimation problems can be formulated for other approximate modeling formalisms, such as models linear in metabolite concentrations [START_REF] Delgado | Metabolic control analysis using transient metabolite concentrations[END_REF], log-lin models [START_REF] Hatzimanikatis | Effects of spatiotemporal variations on metabolic control: Approximate analysis using (log)linear kinetic models[END_REF], and generalized mass-action models [START_REF] Savageau | Biochemical Systems Analysis: A Study of Function and Design in Molecular Biology[END_REF]. Thus, the results that will be discussed below for linlog models can be adapted straightforwardly to the approximate kinetic modeling formalisms mentioned above.

On the contrary, most results are not applicable to formalisms where parameter estimation cannot be reduced to linear regression, such as classical Michaelis-Menten kinetics and convenience kinetics [START_REF] Liebermeister | Bringing metabolic networks to life: Convenience rate law and thermodynamic constraints[END_REF].

IDENTIFIABILITY OF LINLOG MODELS

The problem of identifiability refers to the ability to unambiguously extract parameter values of a model structure from experimental data. Here we focus on linlog models. We shall first discuss the problem from the perspective of structural identifiability. For practical purposes, this is equivalent to answering the question whether each parameter can be uniquely reconstructed from arbitrarily rich and errorless data sets. Structural identifiability forms the basis for studying practical identifiability, i.e. the ability to estimate parameter values from real data sets, which will be discussed further below.

The system, described by Equations ( 1)-( 3), is parametrized by the parameter vector p = [a nx+nu+1) . For simplicity, let us make the standing assumption that N diag(e)B x is invertible for all p ∈ P and e ∈ E (generalizations are possible). For varying values of p, this defines a class of models M p mapping inputs (e, u) ∈ E × U to outputs (J p (e, u), x p (e, u)) ∈ V × X, the latter being the steady-state values of v, x (in accordance with metabolic control analysis standards [START_REF] Heinrich | The Regulation of Cellular Systems[END_REF], symbol J is used in place of v for steady state fluxes). More precisely, we have M p : E × U → V × X : (e, u) → J p (e, u), x p (e, u) (8) where, plugging (2) into (3) to express ln x p (e, u), and writing e in place of diag(e) for shortness,

B x B u ] ∈ P ⊆ R m×(
J p (e, u) = e • (a + B x • ln x p (e, u) + B u • ln u), ln x p (e, u) = -(N eB x ) -1 • N e • (a + B u • ln u). (9)
Since we are observing the system in steady state, inputs and outputs are not time-varying signals but just fixed vectors.

In agreement with Section 2, where the identification problem is split into the identification of each reaction separately, we look at the identifiability of the parameters of the generic ith reaction, and say that a model is identifiable if all its reactions are.

Identifiability from a structural perspective

We adapt the definition from [START_REF] Ljung | System identification, theory for the user[END_REF] to our context as follows.

Definition 1. A reaction i of model M p is identifiable at p * if there exists I ⊆ E × U such that, for all p ∈ P , (J p ) i , x p | I = (J p * ) i , x p * | I ⇒ p i = p * i . (10)
Here J p and x p are seen as vector functions of e and u, subscript "i" indicates the ith element (row) of the vector, and "| I " indicates the restriction of the functions on I.

Turned another way, a reaction i is considered identifiable for a particular model parametrization p * if no p ∈ P with p i = p * i exists such that the velocity of the reaction predicted by M p and M p * from identical steady state concentrations is identical for all possible system inputs. Note that this definition is applicable to any form of the reaction rates (2) (provided suitable definition of the parameters p).

How can we apply Definition 1 to the analysis of identifiability of linlog models? The following proposition establishes the link between this definition and the uniqueness of the solution to Problem 1 (in the form (7)). Given the input set I = {(e 1 , u 1 ) , • • • , (e q , u q )} and a "true" parameter vector p * , let J k * and x k * denote the outputs J p * (e k , u k ) and x p * (e k , u k ), respectively, with k = 1, . . . , q. Proposition 1. A reaction i of M p is identifiable at p * if and only if there exists I ⊆ E × U such that the solution of the equation

W * •i = Y * B •i , with W * •i = J 1 * e 1 -J * e i • • • J q * e q -J * e i T , Y * = ln x 1 * -ln x * • • • ln x q * -ln x * ln u 1 -ln u • • • ln u q -ln u T , is unique in the parameters B .i = ([B x B u ] i. ) T . In turn, this happens if and only if Y * •C(i) (the submatrix of Y * formed by the columns indexed by C(i)) is full column- rank.
Clearly this condition is equivalent to the uniqueness of the solution of (7), the difference being that the identifiable parameter vector p i also contains the element a i .

A standard approach for studying the rank of Y *

•C(i) is the Singular Value Decomposition (SVD)

Y * •C(i) = U • diag(s 1 , s 2 , ..., s n b ) • V T (11)
with

n b = |C(i)|, U ∈ R q×n b and V ∈ R n b ×n b orthonormal matrices and s 1 ≥ • • • ≥ s n b ≥ 0 the singular values of Y * •C(i) .
In the presence of dependencies between the columns, there exists an index r with 1 ≤ r < n b such that s r+1 = • • • = s n b = 0, and Y * •C(i) is of rank r. In order to illustrate the identifiability properties of a metabolic reaction, consider the following example. Example 1. Consider the negative feedback network structure shown in Figure 1(a). The network includes n = 2 internal metabolites, no external metabolites, and m = 3 reactions (enzymes). The network structure and linlog model parameters are SVD analysis reveals that, in the latter case, s 2 ≪ s 1 , which hints at an ill-conditioned estimation problem. This point will be further developed in the next example.

N = 1 -1 0 0 1 -1 , a =
Second, some pathological parameterizations may give rise to a nonidentifiable model. Indeed, if a is in the span of B x , then the unique solution of N diag(e)(a + B x ln x) = 0 corresponds to the value of ln x satisfying B x ln x = -a, i.e. it is independent of e. Thus, no matter the number of experiments q, the rank of Y * is at most 1 (for a = 0, the solution is ln x = 0, i.e. Y is the zero matrix which has rank 0). Thus, in the light of Proposition 1, the model is not identifiable. In our example, this is the case for a 1 = -7.5491 and B 21 = -7.2961.

From the example above, it is clear that, for some specific values of the parameters, a model may be non-identifiable even if it is identifiable for non-pathological parametrizations. In the light of this, a convenient generalization of Definition 1 to the identifiability of a network structure can be obtained following [START_REF] Walter | Identification of parametric models[END_REF], who state that a model in M p is identifiable if it is identifiable at almost every p * ∈ P (in the sense of the Lebesgue measure on P ). Whence, the negative feedback network structure of Example 1 is identifiable in the sense of Walter and Pronzato. The identifiability criterion of Definition 1 does not hold for the "rare" parameter combinations p * such that a ∈ span(B x ). A second observation, following from the example above, is that our theoretical definition of identifiability is actually too restricted to be practically useful. If we look at Figure 1(c), we see that strong collinearities exist between the metabolite concentrations x 1 and x 2 . This makes the informativity of the experiments limited. Again, we will analyze later the implications of this on the achievable parameter estimation performance. Moreover, the definition assumes that the measurements are not corrupted by noise, which is even less realistic.

We therefore need to weaken the definition of identifiability in order to make it more suitable for applications to actual data on metabolism. While taking into account realistic assumptions on the experimental datasets, i.e., measurements available in a limited amount and affected by experimental error, this notion of identifiability should draw upon the theoretical notion of model identifiability discussed above.

Identifiability from a practical perspective

Let I be a fixed set of q inputs (external metabolites and enzyme concentrations), and let O be the set of the correponding system outputs (fluxes and steady-state concentrations of internal metabolites determined by M p * ).

Consider the problem of estimating the parameters B •i of reaction i given observations of I and O affected by measurement error. An estimator B•i of B •i is a function of the observations of I and O, well-defined for every possible (a priori unknown) value of p * ∈ R n b (compare [Ljung, 1999, §7.4]). Since, due to noise, the observations are stochastic variables, B•i is itself a stochastic variable. Therefore, one cannot hope to estimate B •i exactly, but only within a certain degree of approximation, except possibly for a few "adverse" outcomes of the measurement error. In this spirit, we define identifiability in terms of the existence of an estimator satisfying prespecified statistical requirements. In doing this, we restrict attention to the nonzero entries of B •i , i.e., B C(i)i . Let B i ⊂ R n b be a bounded neighbourhood of the origin, and let α ∈ (0, 1). Definition 2. For a given I ⊆ E × U , a reaction i of M p is identifiable at p * with uncertainty B i and confidence level 1 -α if there exists an estimator BC(i)i such that

P p * [ BC(i)i -B C(i)i ∈ B i ] ≥ 1 -α, (12)
where P p * is the probability measure induced by M p * .1

In this view, the experimentalist, or the modeler, sets the requirements (estimation accuracy and confidence level) that the estimates must fulfill in order to be useful, via the a priori specification of B i and α. Then, the possibility of fulfilling (12), i.e. the practical identifiability of the model, depends on the system itself and on the richness of the input set I. In general, the larger the I, the tighter the requirements that one can fulfill (i.e. the smaller the values of B i and α for which practical identifiability in the sense of Definition 2 holds). Note that this is conceptually different from what suggested by [START_REF] Raue | Structural and practical identifiability analysis of partially observed dynamical models by exploiting the profile likelihood[END_REF], where the definition of practical identifiability only requires that the uncertainty on the parameter estimates (as defined via the profile likelihood) is bounded (though arbitrarily large) for a specific (not necessarily optimal) choice of the estimator. In an alternative view, one may start from a given input set I, and look for the choices of α and B i that ensure satisfaction of ( 12). Here in turn, one may fix α and look for the B i that makes (12) achievable, or fix the acceptable estimation uncertainty B i and establish at what confidence level α this performance can be attained.

In all cases, the natural questions that arise are how Definition 2 can be verified in practice, how this notion of identifiability depends on the structural system identifiability discussed in the previous section, and what B i may look like. To answer these questions one needs to further specify the properties of the data, i.e., the (stochastic) "measurement model". Recall from Section 2 the following assumptions. The metabolite concentrations forming Y are measured without error. For i = 1, . . . , q, the noisy versions W.i of the

W •i obey the model W•i = W •i + ε •i = Y •C(i) B C(i)i + ε •i , (13)
with ε •i ∼ N (0, Σ εi ) and Σ εi > 0. The following proposition answers the questions above. Proposition 2. If a reaction i of M p is identifiable at p * in the sense of Definition 1 then, for every α ∈ (0, 1), it is identifiable in the sense of Definition 2 with confidence level at least 1 -α for any uncertainty set B i containing the (1 -α)-confidence ellipsoid of a zero-mean Gaussian distribution with variance (Y T

•C(i) Σ -1 εi Y •C(i) ) -1 .
The proof relies on the use of minimum variance estimators, as dictated by standard results in linear estimation theory [Ljung, 1999, Appendix II]. From now on, estimation will be performed based on this type of estimator.

Observe that B i implicitly depends on the choice of inputs I via the structure of W •i and Y •C(i) . Typically, the larger q, the smaller B i for a fixed α. We argue that similar identifiability results can be derived even in cases where the noise is not Gaussian and metabolite measurements are affected by stochastic error, at the price of a much more complicated characterization of B i . Finally, one may speak about identifiability of the whole model, e.g. by requiring that each reaction is individually identifiable with a given confidence level α and uncertainty set B i ∈ R m×n b . Alternatively, one may require that all reactions be simulsuch as BC(i)i B C(i)i (which makes the reaction identifiable for any α and B i but cannot be built without the knowledge of B C(i)i itself). Unfortunately, this is not a good choice in general, in that the uncertainty set B i may severely depend on p * , as we shall see later on in Example 2. Hence we stick to Definition 2 with the understanding that any such triviality is avoided.

taneously identifiable with confidence level 1 -α and a suitably defined joint uncertainty set.

There is a clear link between Definition 2 and singular values of the regression matrix Y •C(i) , i.e., the rank analysis of the previous section. Assume for simplicity that all error terms have the same variance Σ εi = σ 2 I, with σ > 0. From Proposition 2, if B C(i)i is identifiable with uncertainty B i and confidence 1 -α, then there exists an estimator (the minimum variance estimator) such that, with probability 1-α, the estimates of B C(i)i will lie in an ellipsoid centered at B C(i)i with axes length proportional to the square roots of the singular values of σ 2 (Y T

•C(i) Y •C(i) ) -1 , i.e., to σ[s -1 1 , . . . , s -1 n b ] (as before, s l denotes the lth singular value of Y •C(i) , l ∈ {1, . . . , n b }). If these eigenvalues are vastly different, the ellipsoid is skewed in the direction of the singular vectors of Y T

•C(i) Y •C(i) associated to the smallest singular values of Y •C(i) . In this sense, if a data matrix Y •C(i) is ill-conditioned (s 1 ≫ s n b ), i.e.
, some data vectors are nearly collinear, parameter B C(i)i appears as nonidentifiable even if rank(Y •C(i) ) = n b (that is, even if s n b > 0). In other words, this makes it possible to detect that, for some r < n b , discrepancies between values of s r+1 , • • • , s n b and s 1 can make the estimates of B C(i)i solving regression poorly determined. Example 2. To illustrate the implications for parameter identifiability of a poorly conditioned data matrix, consider the parameter estimation results from noisy and finite datasets for the two different identifiable parametrizations of the model of Example 1. As in Example 1, datapoints were simulated from random values of enzyme concentrations. Noise was added to W by drawing values from normal distributions with standard deviations proportional to the corresponding elements of W . Two different dataset sizes (q = 20 and q = 100) and two different noise levels (20% and 50%, meaning that 99% of the noise samples fall within 20% and 50% of the corresponding values in W ) were tested, for a total of 4 experimental scenarios for each model parameterization. For each scenario, 100 datasets were simulated and the corresponding estimates drawn for reaction 1 are reported in the scatter plots of Figure 2(a) (a 1 = 0.0297 and B 21 = -0.0073) and 2(b) (a 1 = 0.0297 and B 21 = -7.2961).

An immediate observation is that the shape of the 95%confidence ellipse of the parameter estimates is different for the two different model parameterizations. While estimation accuracy for B 11 and B 21 is comparable in the case of smaller feedback (B 21 = -0.0073), the shape of the uncertainty ellipse becomes very skewed in the case of stronger feedback (B 21 = -7.2961). In particular, in this case, estimation accuracy is much higher for B 11 than for B 21 regardless of the features of the dataset (note the change in the scale of the vertical axes of the plots from Figure 2(a) to Figure 2(b)). As apparent from the plots, larger and/or less noisy datasets improve estimation performance. However, ameliorating the estimation of B 21 requires extremely large/high-quality datasets. In other words, estimating the model (i.e. all its parameters) accurately requires a significant increase in the experimental effort, even if most parameters are easy to estimate (see also [START_REF] Gutenkunst | Universally sloppy parameter sensitivities in systems biology models[END_REF]). Following upon Example 1, it is also apparent that the skewed estimation uncertainty is related to the poor conditioning of the data matrix Y in the case of stronger feedback (the shape of the ellipsoid is determined by the ratio of the singular values of Y ). In terms of practical identifiability, assuming a modeler has set a maximum allowable uncertainty B 1 for some confidence level α, it is clear that in this case the system will not be practically identifiable (even if the model is structurally identifiable at the given p * ), unless B 1 is large enough, i.e. rather sloppy estimates are deemed acceptable.

CONCLUSIONS

In this paper we have discussed the notion of structural and practical identifiability for quantitative dynamical metabolic network models. Assuming a fixed model structure, we have introduced definitions of general applicability for the identifiability of the model parameters a priori and a posteriori from a given dataset. In the relevant case of approximate linlog modeling, using standard tools from linear algebra and estimation theory, we have linked the two concepts, provided conditions for structural and practical identifiability, and methods for testing these conditions in practice. Concepts and results were illustrated by a simple example. Current research directions include model reduction, as well as methods for handling missing data entries and noise on metabolite concentrations.

  x = B T ) where different values of a 1 ∈ R ≥0 and B 21 ∈ R ≥0 (the coefficient that determines the strength of the feedback regulation) will be considered. For all values of the enzyme concentrations e i > 0, with i = 1, 2, 3, and all a 1 , B 21 ∈ R ≥0 , the equation N v(x, e) = N diag(e)(a + B x ln x) = 0 yields a unique steady-state solution ln x = -(N diag(e)B x ) -1 N diag(e)a. One first observation is that different values of a 1 and B 21 may lead to very different properties of the matrix Y * even when this remains full rank, i.e. the system is structurally identifiable. For a 1 = 0.0297 and two different values of B 21 , scatter plots of the steady state solutions ln x from 1000 randomly generated samples of e are reported in Figures 1(b) (B 21 = -0.0073, weeker feedback action) and 1(c) (B 21 = -7.2961, stronger feedback action). In Figure 1(b), steady-state metabolite concentrations are spread over a two-dimensional region, while in Figure 1(c) they are essentially aligned along a one-dimensional line.

Fig. 1 .

 1 Fig. 1. An example of a small metabolic network with negative feedback. (a) Structure of the network. (b)-(c) Steady-state metabolite concentrations from 1000 randomly generated enzyme concentrations, for (b) a 1 = 0.0297, B 21 = -0.0073 and (c) a 1 = 0.0297, B 21 = -7.2961 (see Example 1). Main figures are crops of the scatter plots of the data (blue markers) with equally scaled axes, emphasizing the shape of the data clouds. All generated data points are shown in the inset scatter plots (red markers), in arbitrary axes.

Fig. 2 .

 2 Fig. 2. Estimates of the parameters (B x ) 1• = [ B 11 B 21 ], mediating the effects of x 1 and x 2 in reaction 1, from simulated steady-state data, for a linlog model of the network of Figure 1(a) with two different parametrizations: (a) a 1 = 0.0297, B 21 = -0.0073; (b) a 1 = 0.0297, B 21 = -7.2961. In each panel, scatter plots are reported for four different experimental scenarios specified by the couple (% noise level, number of datapoints q): (20%, 20) (blue), (20%, 100) (green), (50%, 20) (red), (50%, 100) (magenta). 95%confidence ellipses are drawn for each scenario (solid lines). Reference parameter values are indicated by horizontal and vertical dotted lines, and are given in Example 1 (see Example 2 for other details).

Strictly speaking, a better version of Definition

would require that condition (12) holds for all p * within a sufficiently large subset of P . This would automatically rule out trivial definitions of BC(i)i