M. D Kerim

L ; Evanno

L Ferrié

Stereodivergent Total Syntheses of (+)-Mycaperoxides C, D, G Methyl Ester and (-)-Mycaperoxide B

Mycaperoxides are natural endoperoxides isolated from different Mycale genus sponges, showing significant antiviral or antibacterial activities. We report herein the first total syntheses of representative congeners of this family from sclareol using a stereodivergent approach. Thus, an innovative oxidative ring expansion of cyclobutanol was used to bring the 1,2-dioxane subunit, and a Mukaiyama aldol reaction on peroxycarbenium species was utilized to install the propionic acid subunit. During the study toward (+)-mycaperoxide D methyl ester (2), the isolation of the eight possible diastereomers under their ethyl thioester form allowed to build a pertinent database for further NMR assignment studies. Thus, we completed the total syntheses of (+)-mycaperoxides D, C, G methyl ester, and (-)-mycaperoxide B in 11 to 15 steps, confirming their original assignment.

INTRODUCTION

The endoperoxide motif is an uncommon chemical function in natural products. However, a relatively large group of marine norsesterterpene metabolites represent a significant class of natural endoperoxides. [START_REF] Casteel | Peroxy Natural Products[END_REF][START_REF] Casteel | Peroxy Natural Products[END_REF][START_REF] Liu | Peroxy Natural Products[END_REF] Indeed, more than 70 members have been isolated from various sponges of the class Demospongiae, including the genera Prianos, Sigmosceptrella, Lacuntrulia, Diacarnus or Mycale. Structurally, these endoperoxides have a common backbone consisting of a 1,2-dioxane ring substituted by a propionic acid subunit on one side and a sesquiterpene subunit connected vie a methylene bridge on the other side. The main difference between all members lies in the nature of the sesquiterpene side chain and, in particular, the degree of cyclization of the chain, the intervention of methyl transposition, and the enantiomeric series of this chain. The absolute configuration of the stereogenic centers of this subclass of natural products is not identical from one product to another, depending mainly on their origin. Despite many members, this subclass has been barely studied at the synthetic level. Indeed, Harwood reported several synthetic studies towards the total synthesis of mycaperoxide B, leading, with difficulty, to the obtention of four inseparable diastereomers, in about 16 steps in total, sharing the gross structure of the target. However, none of them corresponded to the natural product. [START_REF] Harwood | Singlet Oxygen Cycloaddition to Dienes: A Biomimetic Approach to Marine Endoperoxides?[END_REF][START_REF] Silva | Probing a Biomimetic Approach to Mycaperoxide B: Hydroperoxidation Studies[END_REF][START_REF] Drew | Convergent Synthesis of the Putative Biogenetic Precursor of Mycaperoxide B and a Norsesterterpene Triene Isolated from an Australian Sponge[END_REF][START_REF] Silva | Towards the Total Synthesis of Mycaperoxide B: Probing Biosynthetic Rationale[END_REF] Later, Seifert reported the total synthesis of diacarnoxide C in about 20 steps but obtained as an inseparable mixture of cis and trans-1,2-dioxanes with epimerization at C6 (Figure 1). [START_REF] Schneider | Towards the Total Synthesis of the Norsesterterpene Diacarnoxide C[END_REF] Recently, we have developed a synthetic method to construct 1,2-dioxanes from the oxidative expansion of cyclobutanols, [START_REF] López | Oxidative Ring Expansion of Cyclobutanols: Access to Functionalized 1,2-Dioxanes[END_REF] followed by the alkylation of the resulting peroxyacetals. [START_REF] López | Oxidative Ring Expansion of Cyclobutanols: Access to Functionalized 1,2-Dioxanes[END_REF][START_REF] Pinet | Access to Functionalized 3,5-Disubstituted 1,2-Dioxolanes under Mild Conditions through Indium(III) Chloride/Trimethylsilyl Chloride or Scandium(III) Triflate Catalysis[END_REF][START_REF] Pinet | Synthesis of 3,5-Disubstituted 1,2-Dioxolanes through the Use of Acetoxy Peroxyacetals[END_REF] Thus, we felt it would be particularly well suited to the total synthesis of several members of this family of norsesterterpene endoperoxides. Mycaperoxides were targeted as some members have been reported to have significant antiviral activities (IC50 = 0.25-1.0 µg/mL against vesicular stomatitis virus and herpes simplex virus type-1) and antibacterial properties (inhibition of the growth of gram-positive bacteria Bacillus subtilis and Staphylococcus aureus).

RESULTS AND DISCUSSION

In terms of synthetic access, these molecules, mainly mycaperoxides B, C, D, and G, would be strategically accessible from sclareol (Salvia sclarea) used as starting synthon. Concerning the chirality control strategy, mycaperoxide B (4) would be directly accessible from sclareol as the latter corresponds to the proper enantiomeric series of the decalin moiety. For the mycaperoxides C, D, and G, their antipodes (ent-mycaperoxides C, D, and G) would be accessible by total synthesis from sclareol. For these three compounds, methyl esters (1), (2), and (3) were targeted as the natural products were isolated and characterized in this form after derivatization with diazomethane. [START_REF] Capon | Mycaperoxides F and G and a Related Norterpene Ketone from Southern Australian Marine Sponges, Mycale Species[END_REF][START_REF] Capon | Cyclic Peroxides and Related Norterpenes from a Southern Australian Marine Sponge, Mycale Sp[END_REF] Thus, a stereodivergent total synthesis was collectively designed for all these products. They could be obtained by a Mukaiyama aldol reaction performed on an acetoxy-endoperoxyketal. [START_REF] López | Oxidative Ring Expansion of Cyclobutanols: Access to Functionalized 1,2-Dioxanes[END_REF][START_REF] Pinet | Access to Functionalized 3,5-Disubstituted 1,2-Dioxolanes under Mild Conditions through Indium(III) Chloride/Trimethylsilyl Chloride or Scandium(III) Triflate Catalysis[END_REF][START_REF] Pinet | Synthesis of 3,5-Disubstituted 1,2-Dioxolanes through the Use of Acetoxy Peroxyacetals[END_REF] The 1,2dioxane ring of the latter could be constructed from a cyclobutanol by ring expansion with molecular oxygen mediated by Co(acac)2. [START_REF] López | Oxidative Ring Expansion of Cyclobutanols: Access to Functionalized 1,2-Dioxanes[END_REF] The cyclobutanol could be easily obtained from a methyl ketone derived from sclareol. The four required decalins would thus be prepared from the same starting synthon via divergent functional arrangements.

Scheme 1. Stereodivergent strategy from sclareol to access ent-mycaperoxides C, D, and G, and mycaperoxide B

The synthesis of the different targeted consisted first of the functionalization of sclareol (5). The latter underwent an oxidative cleavage with KMnO4 leading to methyl ketone 6. [START_REF] Ruzicka | Mitteilung). Oxydation Des Sclareols Mit Kaliumpermanganat[END_REF] A Wittig olefination then led to the formation of cyclopropylidene 7, [START_REF] Camacho | Palladium-Catalyzed Addition of Alcohol Pronucleophiles to Alkylidenecyclopropanes[END_REF] the precursor of mycaperoxide D. To obtain the other congeners, a hydroxyl elimination gave rise to a mixture of isomers 8 and 8' in an 85:15 ratio and thus gave access to the other required decalins. Epoxidation of the 8/8' mixture followed by Wittig olefination provided the three separable isomers 9a, 9b, and 9c in a 61:30:9 ratio. Compounds 9a and 9b were subjected to reduction with LiAlH4. From 9a, hydride attack led exclusively to 10a, the precursor of ent-mycaperoxide C. The same reaction performed on 9b gave only 10b, the precursor of mycaperoxide B and ent-mycaperoxide G. The regioselectivity of hydride attack is in agreement with the Fürst-Plattner rule. [START_REF] Fürst | Über Steroide und Sexualhormone. 160. Mitteilung. 2α, 3α-und 2β, 3β-Oxidochlolestane; Konfiguration der 2-Oxy-cholestane[END_REF] The synthesis of ent-mycaperoxide D was firstly studied for its better accessibility. Thus, after the protection of the hydroxyl in the form of a silyl ether (TES), the cyclopropylidene was epoxidized. In situ acid-catalyzed rearrangement led then directly to cyclobutanone 11. After reduction with NaBH4, the resulting cyclobutanol underwent a ring enlargement reaction in the presence of Co(acac)2, affording 1,2-dioxanol 12 in good yield. [START_REF] López | Oxidative Ring Expansion of Cyclobutanols: Access to Functionalized 1,2-Dioxanes[END_REF] We observed that the decalin slightly induced stereoselectivity during the cyclobutanol oxidative ring opening. The hydroxyl was acetylated in the presence of pyridine, [START_REF] Jamey | From Ring-Expansion to Ring-Contraction: Synthesis of γ-Lactones from Cyclobutanols and Relative Stability of Five-and Six-Membered Endoperoxides toward Organic Bases[END_REF] and the resulting acetoxy-endoperoxyacetal 13 reacted with silylketenethioketal 14 [START_REF] Deroy | Total Synthesis of (+)-Cystothiazole A[END_REF] as well as a catalytic amount of Sc(OTf)3. [START_REF] López | Oxidative Ring Expansion of Cyclobutanols: Access to Functionalized 1,2-Dioxanes[END_REF] At this step, the first resolution of diastereomers [2,3-syn/2,3-anti and cis/trans-1,2-dioxanes] was thus performed by silica gel chromatography. Therefore, three fractions (F1, F2, and F3) were obtained in 8, 32, and 31% yield, respectively; each containing a mixture of defined diastereomers among the eight possible (Scheme 3, A). The stereoselectivity of the reaction was estimated to be 4:1 in favor of cis-1,2-dioxanes, but no selectivity was observed for the syn/anti relationship at C2-C3. After cleavage of the silyl ether by acid methanolysis performed independently on each fraction (F1, F2 and F3), a further separation of diastereomers was again done to deliver thioesters 16ah. All cis-1,2-dioxanes could be isolated as single stereoisomers but, trans-1,2-dioxanes, do not allow a separation of their 2,3,6epimer and were thus obtained as a mixture of two diastereoisomers. For the trans-1,2-dioxanes mixtures, a selectivity about 55:45 in favor of one particular 2,3,6-isomer, allowed the differentiation of carbon for attribution.

The structures of the different diastereomers were attributed as follow by using Capon's empirical rules [START_REF] Capon | Structural and Stereochemical Studies on Marine Norterpene Cyclic Peroxides[END_REF] based on the [START_REF] Capon | Mycaperoxides F and G and a Related Norterpene Ketone from Southern Australian Marine Sponges, Mycale Species[END_REF] C NMR chemical shift. Firstly, cis/trans-1,2-dioxanes were assigned thanks to [START_REF] Capon | Mycaperoxides F and G and a Related Norterpene Ketone from Southern Australian Marine Sponges, Mycale Species[END_REF] C NMR and the characteristic chemical shifts of the methyl groups at position C6 (20.0 -21.0 ppm for transdioxanes vs. 23.5 -23.7 ppm for cis-dioxanes) and of position C7 (43.0 -43.9 ppm for trans-dioxanes vs. 37.7 -39.1 ppm for cis-dioxanes); Secondly, the 2,3-syn or 2,3-anti configurations were then attributed based on the characteristic chemical shifts of the methyl groups at position C2 (13.2 -13.3 ppm for 2,3-anti vs. 14.4 -14.9 ppm for 2,3-syn products) and of position C4 (23.4 -23.7 ppm for 2,3-anti vs. 22.6 -22.8 ppm for 2,3-syn products). Thirdly, 2,3,6-epimer pairs [2,3-syn, trans-dioxane 15a vs. 15b, 2,3-syn, cis-dioxane 15c vs. 15d, etc.] were the most difficult to discriminate due to the slight chemical shift differences due to the spacer between decalin and 1,2-dioxane units. Structure of cis-dioxanes 16g and 16h were deduced from the similarities observed on [START_REF] Capon | Mycaperoxides F and G and a Related Norterpene Ketone from Southern Australian Marine Sponges, Mycale Species[END_REF] C NMR with natural (-)mycaperoxide D. Although tight chemical shift differences at C7 and C8, these positions are characteristic (± 0.2 -0.8 ppm difference for each 2,3,6-epimer pairs) and cis-dioxanes 16c and 16d were next assigned. A more significant difference can be observed at C11 for cis-dioxanes 16c and 16g versus transdioxanes 16d and 16h (+1.6 ppm), which is indicative of a hydrogen bond of the hydroxyl at C10. Indeed, DFT calculation at the B3LYP/6-31G(d,p) showed that among the lowest energy conformers, an evident hydrogen bond exists between this hydroxyl group and carbonyl group at C1 with a distance of 2.18 Å. (Scheme 2) This particularity also explains the ability of cisdioxanes 6c-d and 6g-h to be isolatable as lone diastereomers compared to protected compounds 15a-h and trans-dioxanes where this hydrogen bond is not possible. Thus, collectively taking into account the 13 C-NMR chemical shifts at positions C7, C8, and C11, 2,3,6-epimers 16a, 16b, 16e, and 16f were assigned. Afterward, the synthesis of (+)-mycaperoxide C methyl ester (1) was investigated. Thus, endoperoxyacetal 19 was obtained in five steps from 10b, using a similar sequence to those described for the total synthesis of 2. A notable difference between the two synthetic routes lies in using TMS as the protecting group for the hydroxy group at C10. The conformation of 10b did not allow the introduction of TES; consequently, only the use of TMS-triflate with Et3N succeeded in the protection step. However, the use of TMS led to an unexpected consequence. Upon Mukaiyama-type aldol reaction, in the presence of Sc(OTf)3 on 19, the hydroxyl at C10 was eliminated. The use of less active Y(OTf)3 10 allowed the aldolization to be carried out without significant side-reaction (Scheme 3, B). The selectivities were comparable for the stereocenters at C6 and C2 to that was obtained in the total synthesis of 2.. However, during this synthesis, the formation of the cis-dioxane rings (C3) was predominant with a 2:1 cis:trans ratio. As described in the total synthesis of 2, we obtained a mixture of diastereomers distributed over three fractions. The fraction F3 also contained the isomer corresponding to the natural product according to [START_REF] Capon | Mycaperoxides F and G and a Related Norterpene Ketone from Southern Australian Marine Sponges, Mycale Species[END_REF] C NMR data. After deprotection of the alcohol function, the compound of interest 21h was separated from its isomer 21g. The final transesterification performed on 21h led to the synthesis of the (-)-mycaperoxide methyl ester C (1). Spectral data matched to the described natural product and, again, the optical rotation [+50 (c 0.198, CHCl3); litt: 12 -71 (c 1.1, CHCl3)] confirmed the synthesis of ent-mycaperoxide C methyl ester and the original assignment of mycaperoxide C (Scheme 3, B).

Regarding (+)-mycaperoxide methyl ester G ( 3) and (-)-mycaperoxide B (4), their total synthesis was planned from compound 10a. The same reaction sequence used to prepare intermediates 13 and 19 (Schemes 3A and 3B) led to the formation of 24 in five steps. However, here again, a particular reactivity of substrate 10a led us to some adaptations. Hydroxyl in position C9 was reluctant to protection and systematically underwent an elimination reaction during all attempts to form a silyl ether. It was thus chosen to leave it free for the rest of the synthesis, its tertiary position inducing a steric hindrance making it less accessible and therefore less reactive. Indeed, the acetylation of 23 was carried out chemoselectively on the hydroxyl of 1,2-dioxane. To avoid further elimination reactions of the hydroxyl at the C9 position during the Mukaiyama aldol reaction, Y(OTf)3 was used again, which led to the obtention of two fractions, F1 and F2 in 87% yield (1,2-dioxane cis:trans = 3:2; 2,3syn:2,3-anti = 1:1). The fraction of interest, F2, contained a mixture of three products, including isomers 25e and 25f, which corresponded to the relative stereochemistry of mycaperoxide B ( 4) and (+)-mycaperoxide G methyl ester (3) respectively. The utilization of preparative HPLC for fraction F2 led to the separation of isomer 25h [START_REF] Capon | The differentiation of 25h from 25g was performed with the help of optical rotation. Indeed, the absolute configuration of the 1,2-dioxane ring determines the optical rotation sign[END_REF] from the mixture of isomers 25e/25f. Thioesters 25e and 25f were saponified, providing mycaperoxide B (4) in an inseparable mixture with 26e isomer in about 53:47 ratio. Assignment of isomers 25e and 25f (or 26e and 4) was tricky due to tenuous differences in the [START_REF] Capon | Mycaperoxides F and G and a Related Norterpene Ketone from Southern Australian Marine Sponges, Mycale Species[END_REF] From the fraction F2 containing isomers 25e-g, elimination of hydroxyl group induced by the action of SOCl2 led selectively to the tetrasubstituted olefin. Compound 27g and a mixture of 27e and 27f were thus obtained. The combination of 27e/27f was transesterified, finalizing the synthesis of the (+)mycaperoxide methyl ester G (3) as a mixture with 28e. As for 26e and 4, the [START_REF] Capon | Mycaperoxides F and G and a Related Norterpene Ketone from Southern Australian Marine Sponges, Mycale Species[END_REF] C NMR data for products 3 and 28e were very close. The most significant differences were observed for the ethylenic carbons at C9 and C10 [ = 0.2 ppm]. The synthetic product with NMR signals matching the most with the natural product was attributed to structure 3 (Scheme 3, C).

CONCLUSIONS

In summary, we were able to construct for the first time several members of the mycaperoxides from sclareol, using a strategy employing an oxidative expansion of cyclobutanol to 1,2dioxane followed by a Mukaiyama aldol reaction on a peroxycarbenium species as key steps. Thus, the methyl esters of (+)mycaperoxides C (1) and D (2) were obtained in 15 and 11 steps, respectively. The methyl ester of (+)-mycaperoxide G (3) and (-)-mycaperoxide B (4) were obtained in 14 and 13 steps but combined with the 2,3,6-epimer. Original assignments were thus confirmed. We highlighted the utility of yttrium compared to scandium in the Lewis acid-catalyzed Mukaiyama aldol reaction, which allowed smoother conditions reducing elimination side reaction with tertiary alkoxides. Although our approach did not allow us to control all the stereogenic centers of the synthetic targets, the different mycaperoxides studied were obtained efficiently accordingly to literature precedents. [START_REF] Silva | Towards the Total Synthesis of Mycaperoxide B: Probing Biosynthetic Rationale[END_REF][START_REF] Schneider | Towards the Total Synthesis of the Norsesterterpene Diacarnoxide C[END_REF] In addition, obtaining several well-identified diastereomers provides a valuable NMR database for the assignment of new natural endoperoxides of the same family. 2,3-syn/2,3-anti and cis/trans-1,2-dioxanes relationships are easily determined by this method. However, the assignment of the 2,3,6-epimers is more challenging, leading to small differences in NMR spectral data. Moreover, the affected positions are also a function of the nature of the decalin unit, which makes a general rule impossible for these pairs of diastereomers. This project also underlines the lack of enantiomeric methods to build organic peroxides, with few reported examples. [START_REF] Lifchits | The Cinchona Primary Amine-Catalyzed Asymmetric Epoxidation and Hydroperoxidation of α,β-Unsaturated Carbonyl Compounds with Hydrogen Peroxide[END_REF][START_REF] Hu | Catalytic Enantioselective Peroxidation of α,β-Unsaturated Aldehydes for the Asymmetric Synthesis of Biologically Important Chiral Endoperoxides[END_REF][START_REF] Lu | Catalytic Enantioselective Peroxidation of α,β-Unsaturated Ketones[END_REF] Therefore, significant work must be undertaken to overcome the actual limitations in the field. [START_REF] Ferrié | Chapter Two -Advances in the Synthesis of 1,2-Dioxolanes and 1,2-Dioxanes[END_REF] Such effort will be reported from our laboratory in due course.

ASSOCIATED CONTENT

The Supporting Information is available free of charge on the ACS Publications website at DOI. General information, detailed experimental procedures, characterization data, cartesian coordinates of 16h, and copies of 1 H and [START_REF] Capon | Mycaperoxides F and G and a Related Norterpene Ketone from Southern Australian Marine Sponges, Mycale Species[END_REF] C NMR spectra of all compounds. (PDF).

12

 12

Figure 1 .

 1 Figure 1. Mycaperoxides A-D et G, and diacarnoxide C.

Scheme 2 .

 2 Scheme 2. Divergent synthesis of sesquiterpene moiety of ent-mycaperoxydes C, D, G et mycaperoxyde B from sclareol.

Figure 2 .

 2 Figure 2. Evidence of a hydrogen bond on diastereoisomer 16h between the carbonyl group of thioester at C1 and hydroxy group at C10. DFT Energy minimization at the theoretical level B3LYP/6-31G(d,p).

Scheme 3 .

 3 Scheme 3. Stereodivergent total of mycaperoxides: A (+)-mycaperoxide D methyl ester; B (+)-mycaperoxide C methyl ester; C (+)-mycaperoxide G methyl ester and (-)-mycaperoxide B.

ACKNOWLEDGMENT

We thank the CNRS for financing the project MYCOXY through a call of project on natural substances.We thank Karine Leblanc and Somia Rharrabti (BioCIS, Châtenay-Malabry)for HRMS analysis and preparative HPLC support, and Jean-Christophe Jullian and Rémi Franco (BioCIS, Châtenay-Malabry) for NMR services.

AUTHOR INFORMATION

Corresponding Author * E-mail: laurent.ferrie@universite-paris-saclay.fr ORCID Laurent Ferrié: 0000-0002-1171-205X Laurent Evanno: 0000-0001-5857-9383

Notes

The authors declare no competing financial interest.

ABBREVIATIONS

acac, acetylacetonate; mCPBA meta-chloroperbenzoic acid; DFT, density functional theory; HPLC, High-performance liquid chromatography; NMR, Nuclear magnetic resonance; TES, triethylsilyl; TMS, trimethylsilyl.


Léon Faure

Bastien Mollet

Wolfram Liebermeister

Jean-Loup Faulon
email: jean-loup.faulon@inrae.fr

Hybrid models enabling neural computations with metabolic networks

Keywords: Artificial Neural Network, Metabolic Network, Mechanistic Modeling, Flux Balance Analysis, Scientific Machine Learning, Hybrid Modeling AMN: Artificial Metabolic Network, ANN: Artificial Neural Network, FBA: Flux Balance Analysis, GD: Gradient Descent, LP(QP): Linear (Quadratic) Programming, ML: Machine Learning, MM: Mechanistic Modeling, PINN: Physics Informed Neural Network, RNN: Recurrent Neural Network. R 2

Constraint-based mechanistic models have largely been exploited to predict the phenotype of microorganisms in different environments. However, phenotype predictions are limited in quality unless labor intensive experiments including the measurement of media uptake fluxes, are performed.

We show how hybrid -mechanistic and neural -models provide ways to improve phenotype predictions. Our hybrid models named Artificial Metabolic Networks (AMNs) surrogate constraintbased modeling, make metabolic networks suitable for backpropagation and, consequently, can serve as an architecture for machine learning. We first show how learning principles brought by AMNs can replace the optimization principle of constraint-based modeling with excellent performances for various in silico training sets. We then illustrate how AMNs outperform mechanistic models with Escherichia coli growth rates measured in 110 different media compositions reaching regression coefficients > 0.76 on cross-validation data. We expect our hybrid AMN models to enhance constraintbased modeling and to prompt new biotechnological applications.

Introduction

The increasing amounts of data available for biological research bring the challenge of data integration with Machine Learning (ML) to accelerate the discovery process. The most compelling achievement within this grand challenge is protein folding recently cracked by AlphaFold 1 , which in the last CASP14 competition predicted structures with a precision similar to structures determined experimentally.

Following this foot step, one may wonder if in the future we will be able to use ML to accurately model whole cell behaviors.

The curse of dimensionality [START_REF] Bellman | Dynamic Programming[END_REF] , i.e., the fact that fitting many parameters may require prohibitively large data sets, is perhaps the biggest hurdle that prevents using ML to build cell models. Obviously, cells are far more complex than single proteins and since the amount of data needed for ML training grows exponentially with the dimensionality 2 , as of today, ML methods have not been used alone to model cellular dynamics.

For the past decades mechanistic (mathematical) models (MMs) have been developed to simulate whole-cell dynamics (cf. Thornburg et al. [START_REF] Thornburg | Fundamental behaviors emerge from simulations of a living minimal cell[END_REF] for one of the latest models). These models encompass metabolism, signal transduction, as well as gene and RNA regulation and expression. Cellular dynamics being tremendously complex, MMs are generally based on strong assumptions and oversimplifications and ultimately suffer from capacities of making predictions beyond the assumptions and the data used to build them.

The MM and ML approaches are based on two contrary paradigms. While the former is aimed at understanding biological phenomena with physical and biochemical detail, it has difficulties handling complex systems; the latter can accurately predict the outcomes of complex biological processes even without an understanding of the underlying mechanisms, but require large training sets. The pros of one are the cons of the other, suggesting that hybrid approaches should be developed. In particular, MMs may be used to tackle the dimensionality curse of ML methods. For instance, one can use MMs to produce in silico data which can be added to experimental data increasing the training set sizes for ML. However, with that strategy, if the model is inaccurate, ML will be trained on erroneous data. One can also embed MMs within the ML process, in this strategy ML and MMs are trained together and the model parameters can be estimated through training, increasing the model predictive capacities.

The issue with this strategy is the difficulty of making MMs amenable to training.

In the current paper we propose an MM-ML hybrid approach in which a whole-cell constraint-based model (CBM), Flux Balance Analysis (FBA), is embedded with ML. For the past three decades, FBA has been the main approach to study the relationship between nutrient uptake and the metabolic phenotype (i.e., the metabolic fluxes distribution) of a given organism, e.g., E. coli, with a model iteratively refined over the years [START_REF] Reed | Thirteen years of building constraint-based in silico models of Escherichia coli[END_REF] . FBA assumes that the metabolic phenotype is at steady state, i.e., a phenotype that is constant in time and in which all compounds are mass-balanced. Usually such a steady state is assumed to be reached in the mid-exponential growth phase. The search for a steady state happens in the space of possible solutions that satisfies the constraints of the metabolic model, i.e., the mass-balance constraints according to the stoichiometric matrix as well as upper and lower bounds for each flux in the distribution. In addition, FBA employs an optimality principle, with one principal objective (usually the 'biomass' production flux) and possibly secondary objectives (e.g., minimize the sum of fluxes in parsimonious Flux Balance Analysis (pFBA), or the flux of a metabolite of interest). For example, to predict a growth rate in a given environment, one would designate the 'biomass' flux (describing the growth rate) as the objective to maximize and set up non-zero upper bounds on some uptake fluxes. These bounds depend on the amounts of molecule transporters, which may change between conditions depending on the cell metabolic strategy. Therefore, finding realistic, condition-dependent bounds on the uptake fluxes requires labor-intensive measurements. In more sophisticated CBM approaches such as molecular crowding FBA [START_REF] Beg | Intracellular crowding defines the mode and sequence of substrate uptake by Escherichia coli and constrains its metabolic activity[END_REF] (mcFBA) or Resource Balance Analysis [START_REF] Goelzer | Quantitative prediction of genome-wide resource allocation in bacteria[END_REF] (RBA), constraints on the resource availability and allocation are added to obtain more biologically plausible metabolic phenotypes, but parameterizing such models requires additional data.

To set bounds and validate prediction results, fluxomic data can be obtained from experiments. With isotopic labeling (like 13-C), one can follow the path of a nutrient in the metabolic network [START_REF] Niedenführ | How to measure metabolic fluxes: a taxonomic guide for 13C fluxomics[END_REF] . Based on metabolomics, one can derive metabolic fluxes from metabolite concentrations and possibly in a timeresolved approach [START_REF] Willemsen | MetDFBA: incorporating time-resolved metabolomics measurements into dynamic flux balance analysis[END_REF] . With transcriptomics, RNA sequencing data is used as input for models to estimate fluxes indirectly, even at the single-cell level [START_REF] Alghamdi | A graph neural network model to estimate cell-wise metabolic flux using singlecell RNA-seq data[END_REF] . Data obtained from several -omics methods can be integrated, constituting a state-of-the-art multi-omics data integration for flux prediction [START_REF] Kim | Multi-omics integration accurately predicts cellular state in unexplored conditions for Escherichia coli[END_REF][START_REF] Lewis | Integration of machine learning and genome-scale metabolic modeling identifies multi-omics biomarkers for radiation resistance[END_REF][START_REF] Zampieri | Machine and deep learning meet genomescale metabolic modeling[END_REF] , which usually makes use of FBA alongside ML methods. Indeed, ML approaches were naturally developed to efficiently integrate such data and to enhance predictions by FBA. However, as described by Sahu et al. [START_REF] Sahu | Advances in flux balance analysis by integrating machine learning and mechanism-based models[END_REF] , the interplay between FBA and ML still shows a gap: some approaches use ML results as input for FBA, others use FBA results as input for ML, but none of them embed FBA into ML, as we do in this paper with the Artificial Metabolic Networks (AMNs) hybrid models.

Hybrid models have recently been developed under different names in biology for signaling pathways and gene-regulatory networks (Knowledge Primed Neural Network 14 , Biologically-Informed Neural Networks [START_REF] Lagergren | Biologically-informed neural networks guide mechanistic modeling from sparse experimental data[END_REF] ) and in physics where models solving partial differential equations (PDEs), such as Physics Informed Neural Network [START_REF] Raissi | Physics-informed neural networks: A deep learning framework for solving forward and inverse problems involving nonlinear partial differential equations[END_REF] (PINN), are available in open source repositories such as SciML.ai [START_REF]SciML: Open Source Software for Scientific Machine Learning[END_REF] . The goal of these emerging hybrid modeling methods is to generate models that comply well with observations or experimental results via ML, but that also use mechanistic insights from MM. The advantages of hybrid models are two-fold: they can be used to parametrize MM methods through direct training and therefore increasing MM predictability, and they enable ML methods to overcome the dimensionality curse by being trained on smaller datasets because of the constraints brought by MM.

The AMN models shown here fit in the emerging hybrid modeling field. AMNs bridge the gap between ML and FBA by computing steady-state metabolic phenotypes with different methods, that can be embedded with ML: a method based on recurrent neural network (RNN), and two trainable methods solving respectively linear programs (LP) and quadratic programs (QP). All these methods are relying on custom loss functions surrogating the FBA constraints. By doing so, our AMNs are mechanistic models, determined by the stoichiometry and other FBA constraints, and also ML models, as they are used as a learning architecture.

We showcase our AMNs with a critical limitation of classical FBA for experimentalists. As already discussed, realistic and condition-dependent bounds on medium uptake fluxes are critical for growth rate (or other fluxes) predictions, but there is no simple conversion from extracellular concentrations, i.e., the controlled experimental setting, to such bounds on uptake fluxes. In fact, these bounds depend on the internal allocation of protein resources to transporters, which itself depends on the cell's metabolic state -the state we are trying to predict. Some methods, such as satFBA [START_REF] Müller | Resource allocation in metabolic networks: kinetic optimization and approximations by FBA[END_REF] , assumes fixed transporter levels and converts medium concentrations to possible uptake fluxes by kinetic models, relying on a Michaelis-Menten value for each uptake reaction. Consequently, this MM approach relies on more assumptions and requires more data than a classical FBA, and does not integrate any large-scale regulation or resource allocation phenomena. More advanced and accurate CBM methods such as Resource Balance Analysis (RBA) are more difficult to set-up and harder to parametrize. As a result, many users rely on classical FBA, which cannot reliably predict growth rates from media compositions. AMNs are used in this study for tackling the same issue: predicting metabolites uptake fluxes from their external concentrations. To do so, where satFBA uses transporter kinetics with parameters that need to be acquired through additional experimental measurements, AMNs use a pre-processing neural layer that is accessible for learning. This neural pre-processing layer aims to capture, effectively, all effects of transporter kinetics and resource allocation in a particular experimental setting, predicting the adequate input for a metabolic model to give the most accurate steady-state phenotype prediction possible. Consequently, AMNs provide a new paradigm for the prediction of fluxes and growth rates: instead of relying on a constrained optimization principle performed for each condition (as in classical FBA), we use a learning procedure on a set of example flux distributions that attempts to generalize the best model for accurately predicting the metabolic phenotype of an organism in different conditions.

Results

As an introductory work, we show how the aforementioned idea of knowledge-primed neural networks inspired an in-house method based on branching ratios of metabolites going to different reactions (Figure 1). This preliminary work having limitations, we then showcase more versatile methods that surrogate FBA with gradient-backpropagation compatibility: a method for solving the linear program (LP) of FBA, and a quadratic program (QP) method solving FBA when measured fluxes (like growth rate) are provided. These two solvers have adequate gradient and loss functions satisfying the constraints of an FBA framework (Figure 2). Then, we used these methods inside hybrid models that can directly learn from sets of flux distributions (Figure 3). These flux distributions used as learning references (i.e., the training sets) are either acquired experimentally or produced through FBA simulations (Table 1 and Figure 4). Finally, we developed a non-trainable AMN reservoir to showcase how the predictive power of FBA can be improved. Indeed, once the AMN has been trained on adequate in silico FBA data, we can fix its parameters, resulting in a gradient backpropagation compatible reservoir that mimics FBA (Figure 5). This reservoir can then be used to tackle the abovementioned issue of unknown uptake fluxes: adding a pre-processing neural layer and training this layer with an experimental dataset, one can predict uptake fluxes from external metabolite concentrations. This neural layer can be reused by any FBA user to improve the predictive power of a metabolic model, with an adequate experimental set-up.

Alternative computation methods to surrogate FBA

The first method (Fig. 1), inspired by previous work on signaling networks [START_REF] Nilsson | Artificial neural networks enable genome-scale simulations of intracellular signaling[END_REF] , learns consensual flux branching ratios found in example flux distributions (i.e., the training set) and represents these ratios in the form of a weight matrix. Since the mass conservation law is the central rule when satisfying metabolic networks constraints, these ratios play a key role in the determination of the metabolic phenotype, i.e., the paths taken by metabolites in the organism. In this approach, we assume that the flux branching ratios remain similar between flux distributions with different bounds on the uptake fluxes. In other words, we postulate that most of the flux distribution examples of our training set resemble a consensual metabolic state, with consensual branching ratios. Consequently, the method learns from example flux distributions, a weight matrix representing these consensual branching ratios. A simple toy model network is shown to demonstrate the functioning of this first method in Fig. 1a. With the usual FBA method, we obtain the flux distribution maximizing the flux of the 𝑣 3 reaction (representing a classical 'biomass' reaction) with an uptake reaction 𝑣 1 (representing a classical uptake reaction) by solving a linear program (Fig. 1b). At steady state, metabolite production fluxes can be calculated from reaction fluxes via the transition matrix 𝑃 𝑣→𝑚 from fluxes to metabolites, and similarly, reaction fluxes can be calculated from metabolite production fluxes using the matrices 𝑃 𝑚→𝑣 and 𝑊 𝑟 (Fig. 1c). For a given metabolite j and a reaction i, the weight wji (in matrix 𝑊 𝑟 ) represents the fraction of metabolite j production flux going to reaction i.

To learn the weights w ji , one first needs to translate the model into a neural-network-like architecture (Fig. 1d). Precisely, in Fig. 1d we start from an initial set of given fluxes and then propagate knowledge about the fluxes through the entire network, each layer corresponding to one step in a discrete flux propagation. Mathematically, each layer is composed of two simple operations that update the M and V vectors, respectively representing metabolites production fluxes and reaction fluxes. Those operations are repeated until convergence, with 𝑣 1 being constant because it is known from the start (i.e., the input of the network). The network shown in Fig. 1d is the unrolled representation of the RNN (named AMN-Wt) depicted in Fig. 1e. As detailed in the Method section -AMN architectures and parameters, the matrix 𝑊 𝑟 can be learned through training, assuming the learned weights are a consensus of the different 𝑊 𝑟 matrices observed in the training set. The steady state fluxes of AMN-Wt iterated 30 times (or more, here we stopped at 30 because the reference data values were reached) equal to those obtained in the reference data (Fig. 1f). 𝑉 𝑖𝑛 is the medium represented by a vector of nutrient uptake fluxes (here 𝑣 1 = 0.1, symbol "-" indicates that no value is provided, in practice we use an 'infinity' value to represent an unbounded flux). The steady-state solution 𝑉 𝑜𝑢𝑡 is calculated by solving a linear program maximizing the objective 𝑐 𝑇 𝑉 = 𝑣 3 , here we used the Cobrapy package [START_REF] Ebrahim | COnstraints-Based Reconstruction and Analysis for Python[END_REF] to compute 𝑉 𝑜𝑢𝑡 (simplex-based algorithm). c. Stoichiometric model matrices. 𝑃 𝑣→𝑚 is the matrix to compute metabolite production fluxes from reaction fluxes, 𝑃 𝑚→𝑣 is a matrix to compute reaction fluxes from the production fluxes of the reaction substrates. When a metabolite is the substrate of several reactions, each reaction gets a fraction of the metabolite production flux, this is depicted in matrix 𝑊 𝑟 (r indicates this matrix is used in recurrence). 𝑀 𝑜𝑢𝑡 is the vector of metabolite production fluxes at steady state, the operator ⊙ stands for element-wise matrix product (Hadamard product), and 𝑅𝑒𝐿𝑈(𝑥) = 𝑚𝑎𝑥 (0, 𝑥). sj,i corresponds to the value of S at the j th row (metabolite) and i th column (flux) and zi is the number of negative elements in column i of S. d. Unrolled neural network built from the model. In this theoretical example an initial flux vector is set, representing only the uptake fluxes, and is mapped to a flux vector covering the entire network. In the initial layer ( 𝑙 0 ) only 𝑣 1 has a value (𝑣 1 0 ). In layer 1, 𝑣 1 value is passed to 𝑚 1 , the production flux for metabolite 1, 𝑚 While the AMN-Wt method is simple to implement it suffers from a drawback. As we shall see later in Table 1, a consensus set of weights leads to a solution when upper bounds (UB) for uptake fluxes are provided, but not when exact bounds (EB) for uptake fluxes are given. Supplementary Fig. S2 shows an example in which two flux distributions, with different EBs on uptake fluxes, lead to different weights for the same networks. Consequently, we cannot assume that AMN-Wt can handle all possible flux distributions in the EB case. To overcome this shortcoming, we next present two alternative methods for solving FBA optimization problems that can accommodate both EB and UB cases for uptake fluxes.

These methods are much closer to the LP or QP optimizations behind FBA, their input can be exact or upper bounds on fluxes. Let us note here that the alternative methods described next are not AMNs per se, but are mechanistic models (MM) with gradient backpropagation compatibility that 'replace' the Simplex-based LP or QP solving algorithms in our AMNs. Importantly, the two methods address two important tasks in flux modeling: optimizing a flux distribution for maximal biomass rate (in LP), as in classical FBA, and fitting a stationary flux distribution to partial flux data (in QP).

The first FBA solver (LP solver, Fig. 2b), derived from a method proposed by Yang et al. [START_REF] Yang | A new neural network for solving quadratic programming problems with equality and inequality constraints[END_REF] , handles linear problems using exact constraint bounds (EBs) or upper bounds (UBs) for uptake fluxes (Vin). That method make use of RNNs for optimization, which is a long-standing field of research [START_REF] Jin | A survey on projection neural networks and their applications[END_REF] inspired by the pioneering work of Hopfield and Tank [START_REF] Hopfield | Neural" computation of decisions in optimization problems[END_REF] . Later, these RNNs were showcased to perform well for solving linear programs [START_REF] Wang | Recurrent neural networks for linear programming: Analysis and design principles[END_REF] and simpler and more efficient solutions were developed over the years [START_REF] Yang | A new neural network for solving quadratic programming problems with equality and inequality constraints[END_REF][START_REF] Ghasabi-Oskoei | An efficient simplified neural network for solving linear and quadratic programming problems[END_REF] . It is important to point out at this stage that these RNNs are non-trainable networks and differ from the contemporary RNNs used in ML. The RNNs developed for optimization are instead recurrent procedures iteratively updating the solution of linear programs. Nonetheless, as we shall see thereafter, the LP solver of Fig. 2b can backpropagate a gradient and can therefore be connected to trainable layers.

As the AMN-Wt method, the LP solver, iteratively computes fluxes to come closer to the steady-state solution (Vout). However, calculations are more sophisticated than for the AMN-Wt, the method integrates the same objective function (c) than the classical FBA Simplex solver and iteratively update the flux vector (V), and the vector (M) representing the dual problem variables also named metabolites shadow prices [START_REF] Varma | Metabolic capabilities of Escherichia coli: I. synthesis of biosynthetic precursors and cofactors[END_REF] (cf. Method -LP-solver -for further details). To assess the validity of the LP solver, we calculated 100 different growth rates for E. coli core model [START_REF] Orth | Reconstruction and Use of Microbial Metabolic Networks: the Core Escherichia coli Metabolic Model as an Educational Guide[END_REF] taken from the BIGG database 21 using the Simplex solver available in the Cobrapy package [START_REF] Ebrahim | COnstraints-Based Reconstruction and Analysis for Python[END_REF] . The E. coli core model 29 contains 154 reactions, 72 metabolites (after duplicating bidirectional reactions), including 20 metabolites, which can be imported to the cell through uptake reactions. We generated 100 different uptake fluxes combinations, varying the 20 uptake fluxes of E. coli core model (cf. Method -Generating Training sets -for further details). Results presented in Fig. 2d exhibits almost identical results as the reference data, with the same initialization, both for known uptake fluxes with EB and for unknown uptake fluxes with UB. The second solver is loosely inspired by the work on Physics-Informed Neural Networks (PINNs), which has been developed to solve partial differential equations matching a small set of observations [START_REF] Cuomo | Scientific Machine Learning Through Physics-Informed Neural Networks: Where we are and What's Next[END_REF] . With PINNs, solutions are first approximated with a neural network and then refined to fulfill the constraints imposed by the differential equations and the boundary conditions. Refining the solutions necessitates the computation of three loss functions, the first is related to the observed data, the second to the boundary conditions and the third to the differential equations. As detailed in the Method -Loss functions derivation and QP solver -sections, we similarly compute losses for the measured fluxes, the EB or UB boundary constraints, and the metabolic network stoichiometry and flux positivity constraints. As in PINN we next compute the gradient on these losses to refine the solution vector V. Unlike with the LP solver, in the present case we do not provide an objective (c) to maximize, but the actual targeted measured data (Vout), consequently the method is named QP 10 because it is equivalent to solving an FBA problem with a quadratic program. Our QP solver reaches acceptable losses as the LP one (Fig. 2d) but requires more iterations.

AMNs: metabolic & neural hybrid models for predictive power with mechanistic insights

While the LP and QP solvers perform well, their main weakness is the number of iterations needed to reach satisfactory performances, more than 10,000 (Fig. 2e). Since our goal is to integrate such methods in a learning architecture, this drawback has to be tackled. As illustrated in Fig. 3a, our solution is to improve our initial guesses for fluxes, by training a prior network (a classical dense ANN architecture) to compute initial values for all fluxes (V 0 ) from exact or upper bounds on uptake fluxes (Vin). This prior network is trained along with the LP or QP solvers with the goal of finding V 0 values that are as close as possible to the optimum flux values (Vout in the reference data set). This is achieved by setting the iteration number of the LP or QP solvers to low values.

In the remainder of the paper, we name Artificial Metabolic Network (AMN), the hybrid model shown in Fig. 3a composed of a neural network followed by a LP or QP solver. AMN results are represented in Fig. 3b with prior networks of different hidden layer sizes. We note that the solver acts as a mechanistic layer as it contains the stoichiometry and bounds of a metabolic network, whereas the prior network acts as a neural layer that can be trained. The performances of all AMN architectures are given in Table 1 using FBA simulated data on two different E. coli metabolic models, E. coli core model 29 and iML1515 [START_REF] Monk | iML1515, a knowledgebase that computes Escherichia coli traits[END_REF] . These models are composed respectively of 154 reactions involving 72 metabolites, and 3682 reactions involving 1877 metabolites (after duplicating bidirectional reactions). In both cases the Simplex-based solver of Cobrapy [START_REF] Ebrahim | COnstraints-Based Reconstruction and Analysis for Python[END_REF] was run to optimize growth rates for different media. Each medium was composed of metabolites found in M9 (minimal media) and additional metabolites (carbon sources, amino acids) crossing the cell membrane.

Table 1. AMNs performances. (1) All SBML models describing different E. coli metabolic models were downloaded from the BiGG database [START_REF] Norsigian | BiGG Models 2020: multi-strain genome-scale models and expansion across the phylogenetic tree[END_REF] , 'Core' stands for the E. coli core model [START_REF] Orth | Reconstruction and Use of Microbial Metabolic Networks: the Core Escherichia coli Metabolic Model as an Educational Guide[END_REF] , EB (UB) stands for exact bounds (upper bounds) for medium uptake fluxes, the iML1515 20 model was reduced following the procedure described in Methods -'Making metabolic models suitable for neural computations'. (2) Training set size and range for the number of metabolites added to M9 growth medium, and method used when running Cobrapy 22 . ( 3) YES or NO if the model contains a neural layer or a mechanistic layer. ( 4) MM stands for Mechanistic Model, AMN stands for Artificial Metabolic Network. Architectures are described in Methods -'AMN architectures and parameters'. ( 5) Number of trainable parameters and epochs, in all cases dropout = 0.25, batch size = 5, the optimizer is Adam and the loss function is the mean squared error between predicted and provided fluxes to which loss constraints are added, see Methods -'Loss functions derivation' for additional details. ( 6) Regression coefficient and Loss values for training set (R²), and cross-validation sets (Q²) between provided growth rate and predicted growth rate. (7) Regression coefficient and Loss values for growth rates for independent test sets not seen during training. Test set sizes are 10% of training set sizes. For ( 6) and ( 7) the performance is displayed as the mean over 5 folds (or over a training set when no cross-validation scheme is performed, i.e., for the MM performances). All architectures presented in Table 1 exhibit excellent regression coefficients and losses for training sets, validation sets and test sets, and this for both models E. coli core 29 and iML1515 [START_REF] Monk | iML1515, a knowledgebase that computes Escherichia coli traits[END_REF] . The MM architectures (without neural layers) can compute losses for metabolic constraints (boundary and stoichiometry) but cannot be used to make predictions on cross-validation or test sets. The AMN architectures can be used to both make predictions (Q 2 ) and compute losses. It is interesting to observe the good performances of AMN-Wt when UB training sets are provided. Indeed, while counterexamples can be found for which AMN-Wt will not work with EB training sets (cf. Fig. S2), we argue in Supplementary Information -AMN-Wt architecture -that AMN-Wt is able to handle UB training sets because the initial inputs (UB values for uptake fluxes) are transformed into suitable exact bound values during training. Yet, the consensual weight matrix (Wr in Fig. 1 or Fig. S1) calculated during training does not have a direct physical meaning. Weights in Wr arise at branch point metabolites that are consumed by at least two reactions (like m1 and m2 in Fig. 1). Consequently, the weights should correspond to flux split ratios, for instance, taking the example of Fig. 1, the metabolite production flux m1 is spliced into w21m1 and w31m1. We show in Fig. S3, that the flux split ratios are conserved for nutrients leading to different metabolite production fluxes if the Michaelis-Menten kinetics parameters (Vmax and Km) of the enzymes catalyzing the reactions involved in the split remain constant. However, Chubukov et al. [START_REF] Chubukov | Transcriptional regulation is insufficient to explain substrate-induced flux changes in Bacillus subtilis[END_REF] have shown experimentally that it was not the case (for B. subtilis) and different nutrients do provide different ratios. This behavior is due to varying enzyme activities, which themselves depend on enzyme concentrations, post-translational modification, and gene regulations. Chubukov et al. [START_REF] Chubukov | Transcriptional regulation is insufficient to explain substrate-induced flux changes in Bacillus subtilis[END_REF] clearly show with experimental evidence that different nutrients give rise to different concentrations for many enzymes, implying that nutrients do have an effect on gene regulations.

In classical FBA, all regulations explaining a given flux distribution are completely ignored, and the flux distribution computation entirely relies on manually setting bounds on uptake fluxes. Therefore, when performing classical FBA, one needs to consider each condition independently from the other, to reach different metabolic phenotypes. Similarly, AMNs attempt to take regulations into account in the neural layer, while keeping the mechanistic layer for metabolic phenotype computations. However, unlike in FBA, AMNs attempt to learn the relationship between a set of flux bounds (on uptake reactions) and the whole steady-state metabolic phenotype. As a result, AMNs are generalizing this relationship for a set of conditions, and not just one as in FBA.

We note that even though weights do not have a physical meaning, AMN-Wt stills exhibits excellent performances, showing that the consensual Wr matrix and the initial V 0 vector (computed through a neural layer from the upper bound Vin) are performant enough. We also note that the weight issue does not arise with AMN-LP and AMN-QP as these architectures do not rely on flux split ratios.

AMNs can be directly trained on experimental datasets with good predictive power

To train AMNs on an experimental dataset, we grew E. coli DH5-alpha in 110 different media compositions, with M9 supplemented with 4 amino acids as a basis and 10 different carbon sources as possibly added nutrients. From 1 up to 4 carbon sources were simultaneously added in the medium at a concentration of 0.4 g.L -1 (more details in Methods -Culture conditions). We determined which compositions to test by choosing all the 1-carbon source media compositions and randomly picking one hundred of the 2-, 3-and 4-carbon sources media compositions (more details in Methods -Generation of an experimental training set). The growth of E. coli was monitored in 96-well plates, by measuring the optical density at 600nm (OD600) over 24 hours. The raw OD600 was then passed to a maximal growth rate determination method, based on a linear regression performed on log(OD600) data (more details in Methods -Growth rate determination).

The experimental dataset was used to train all AMN architectures (-LP, -QP, -Wt). For AMN -LP/-QP we used the architecture plotted in Fig. 3, for AMN-Wt we used the architecture of Fig. S1. In all cases the mechanistic layer was derived from the stoichiometric matrix of the iML1515 20 E. coli reduced model (cf. Method section -Making metabolic networks suitable for neural computations). For all AMNs upper bounds for nutrient uptake fluxes' upper bounds were provided, as exact values remain unknown. Results are provided in Fig. 4. S1) is composed of an input layer of size 38, an hidden layer of size 500, and an output layer of size 550 corresponding to all fluxes of the iML1515 20 reduced model, the size of the recurrent Wr matrix of Fig. S1 is 550x550.

Fig. 4: Benchmarking growth rate predictions by AMNs with experimental measurements. In all panels the experimental measurements were carried on E. coli grown in M9 mixed carbon sources (strain DH5-alpha, model iML1515 20 ). Training and 10-fold stratified cross-validation were performed 3 times with different initial random seeds each for 1000 epochs. All points plotted correspond to predicted values not present in the training set. This was compiled using all predicted values obtained for each cross-validation set. In all cases, means are plotted for both axes (measured and predicted), error bars are standard deviations. Supplementary File 'Data_Fig4.xlsx' compiles raw results used for this figure (see Data availability). a. Performance of the AMN-LP architecture. The neural layer of AMN-LP is composed of an input layer of size 38, an hidden layer of size 500, and an output layer of size 550 corresponding to all fluxes and 1083 metabolite shadow prices of a reduced iML1515 model 20 . The mechanistic layer takes as input the 550+1083 outputs of the neural layer and minimizes the loss between measured and predicted 'biomass' reaction fluxes and the losses of the metabolic network constraints. b. Performance of the AMN-QP architecture. The neural layer of AMN-QP is composed of an input layer of size 38, an hidden layer of size 500, and an output layer of size 550 corresponding to all fluxes of the iML1515 20 reduced model. c. Performance of the AMN-Wt architecture. The neural layer of AMN-Wt (matrix Wi in Figure

AMNs can be used in a reservoir computing framework to enhance the predictive power of traditional FBA solvers

Once an AMN has been trained on a large dataset of simulated data, we can fix its parameters and exploit it in subsequent further learning. As already mentioned in the introduction section, the uptake fluxes of E. coli nutrients, as well as their relation to external concentrations, remain largely unknown: the quantitative rate for each compound may vary between growth media, unlike in classical FBA calculations, where the same upper bound (or zero, if a compound is absent) is used in all cases. This is a common flaw of FBA, that is only partially remedied in satFBA [START_REF] Müller | Resource allocation in metabolic networks: kinetic optimization and approximations by FBA[END_REF] or mcFBA [START_REF] Beg | Intracellular crowding defines the mode and sequence of substrate uptake by Escherichia coli and constrains its metabolic activity[END_REF] , as those relies on supplementary knowledge (and additional experimental measurements). To use the architectures of Fig. 3, we thus need to first convert nutrient concentrations into uptake fluxes. This is achieved by adding a neural layer that precedes the AMN (Fig. 5a). In this architecture the AMN is no longer trainable, only the prior neural layer is, and the AMN acts as a reservoir as in reservoir computing [START_REF] Tanaka | Recent advances in physical reservoir computing: A review[END_REF] .

Among the various architectures we benchmarked, we took as a reservoir the AMN-QP of Table 1 trained on iML1515 [START_REF] Monk | iML1515, a knowledgebase that computes Escherichia coli traits[END_REF] UB dataset.

The predictive power of our AMN-reservoir was assessed by a regression coefficient R 2 = 0.97 (Fig. 5c) on training set, and Q 2 = 0.76 (Fig. 5b) on a test set built from aggregated validation sets during 10-fold cross-validation. To compare with an 'out-of-the-box' FBA approach, we first used an arbitrary value for the present compounds, and a zero value for the absent compounds. Applying these values on the corresponding uptake fluxes upper bounds, and optimizing the 'biomass' reaction flux produced results poorly correlated with experimental measurements. To enhance these results, we searched optimized upper bound values for uptake fluxes (see Methods -Optimization of uptake fluxes upper bounds in FBA). It resulted in a R² = 0.51 performance (Fig. 5f), showing a much weaker fit than our AMNreservoir. Overall, our results indicate that AMNs can be used to substantially increase the predictive capabilities of FBA without costly experimental work. [START_REF] Ebrahim | COnstraints-Based Reconstruction and Analysis for Python[END_REF] , using different scalers on Cmed to compute Vin. Here, Cobrapy [START_REF] Ebrahim | COnstraints-Based Reconstruction and Analysis for Python[END_REF] was run taking as input upper bound uptake fluxes for added metabolites in medium. Intake flux values were set to , and then scaled (0 when the metabolite was absent in the medium and 1*scaler when it was present -see section Method, Optimization of uptake fluxes upper bounds in FBA, for details on optimization procedure)

Fig. 5: Reservoir computing for improving the predictive power of FBA modeling (strain DH5-alpha, model iML1515 20 ). a. AMN used in reservoir computing. The non-trainable AMN-reservoir (dotted gray box) is the AMN shown in Figure 3.b (AMN-QP trained on iML1515 with UB bounds with performances given in the second to last row of Table 1) and is connected to a prior trainable network which purpose is to compute medium uptake fluxes (Vin) from the concentration (Cmed) of metabolites added to medium. Taking as input Vin, the non-trainable reservoir prints out all fluxes including the growth rate. b. Performance of the reservoir computing model shown in panel a, in terms of predicted growth rate vs. measured growth rate. As in Fig. 4, all points plotted correspond to predicted values not present in the training set. This was compiled using all predicted values obtained for each validation set in 10-fold stratified cross validation repeated 3 times with different initial random seeds. c. Performance of the AMN-reservoir with training points instead of predictions. All points plotted correspond to computed values in the same cross-validation scheme as in panel b, using the points of the training sets instead of validation sets. d. Performance of Cobrapy 22 when extracting Vin from panel b predictions to be used as inputs. e. Performance of Cobrapy 22 when extracting Vin from panel c computations to be used as inputs. f. R² between measured growth and computed growth by Cobrapy

Discussion and conclusion

In this study, we showed how a neural network approach, with metabolic networks as a learning architecture, can be used to address metabolic modeling problems. Previous work on RNNs and PINNs for solving constrained optimization problems was re-used and adapted to develop three models (AMN-Wt, -LP and -QP) enabling gradient backpropagation within metabolic networks. The models exhibited excellent performances on FBA generated training sets (Table 1). We also demonstrated the models could directly be trained on an experimental E. coli growth rate dataset with good predictive abilities (Fig. 4). For improved scalability and adaptability, we trained an AMN-reservoir on a large FBA generated training set, and used the reservoir to improve FBA predictions on the experimental dataset. Figure 5 shows that our AMNs far outperform the results obtained by the simplex-based Cobrapy [START_REF] Ebrahim | COnstraints-Based Reconstruction and Analysis for Python[END_REF] FBA solver.

Determining uptake fluxes is a core experimental work required for making FBA predictions realistic.

Here, we developed approaches that get rid of such needs for reaching plausible fluxes distributions.

We did so by backpropagating the error on the growth rate, to find complex relationships between the medium concentrations and the medium uptake fluxes. To this end, we demonstrated the high predictive power of AMNs, and their re-usability in classical FBA approaches. Indeed, FBA developers and users may now make use of our AMN-reservoir method for relating medium uptake fluxes to growth medium concentrations. In this regard, Supplementary File 'Data_Fig5.xlsx' (see Data availability) gives uptake fluxes for the metabolites used in our benchmarking work with E. coli (Fig. 5), these upper bounds for uptake fluxes that can directly be used by Cobrapy [START_REF] Ebrahim | COnstraints-Based Reconstruction and Analysis for Python[END_REF] to reproduce Fig. 5d and5e.

Making FBA suitable for training like we have done in this paper opens the door to improve FBA models.

For instance, in addition to uptake fluxes, AMNs could also be used to search for the coefficients of the biomass reaction appropriate to best fit measurements. So far, these coefficients are derived based on literature, but also using experimental data: growth rate, flux, and macromolecular fractions measures can help finding optimal coefficients 20 . However, these experiments are limited in number and performed once for the reaction parametrization, in a single experimental setup, meaning these coefficients are hardly extrapolated to all possible conditions. Some studies already underline this issue and attempt to efficiently integrate experimental data in the biomass reaction parametrization [START_REF] Lachance | Generating biomass objective functions for genome-scale metabolic models from experimental data[END_REF] . With AMNs, a trainable layer containing the reaction's coefficients could be added, adapting the biomass reaction to any experimental setup.

Another task that AMNs could handle is gene regulation which is not explicitly taken into account in classical FBA. Here a set of genes (including enzymes) corresponding to operons could be encoded via trainable layers either turned on or off. Such AMNs could be trained on a variety of experimental inputs (wider than the carbon source composition as shown in this study) to grasp the complexity of the regulation processes happening in the cell to better explain the end-point metabolic steady-state phenotype of the organism.

Returning back to the curse of dimensionality issue mentioned in the introduction, we systematically searched training set sizes for which 'black-box' ML methods would yield performances similar to our AMN hybrid models. To that end, we trained a simple dense ANN model on training sets of increasing sizes. Results obtained with E. coli core [START_REF] Orth | Reconstruction and Use of Microbial Metabolic Networks: the Core Escherichia coli Metabolic Model as an Educational Guide[END_REF] show that at least 500,000 entries are needed in the training sets to reach losses below 10 -1 (cf. in Supplementary information -AMN and ANN training set sizes and Fig. S7), which according to Table 1 are still 2 orders of magnitude higher than all AMNs losses trained on only 1000 entries. This clearly demonstrates the capacity of hybrid models to reduce training set sizes by constraining the search space through the mechanistic layer.

Beyond FBA and black-box ML improvements, AMNs can also be exploited for industrial applications. Indeed, since arbitrary objective functions can be designed and AMNs can be directly be trained on experimental measurements, AMNs can also be used to optimize media for the bioproduction of compounds of interest or to find the best gene deletion and insertion strategy in typical metabolic engineering projects. In this latter case, a trainable 0/1 layer turning on or off reactions can be added prior to the mechanistic layers of our AMNs. Another potential application is the engineering of microorganism-based decision-making devices for the multiplexed detection of metabolic biomarkers or environmental pollutants. Here, AMNs could be used to search for internal metabolite production fluxes enabling one to differentiate positive samples containing biomarkers or pollutants from negative ones. Such a device has already been engineered in cell-free systems [START_REF] Pandi | Metabolic perceptrons for neural computing in biological systems[END_REF] , and AMNs could be used to build a similar device in vivo by adding a trainable layer after the mechanistic layer which purpose would be to select metabolite production fluxes that best split positive from negative samples.

Methods

Making metabolic networks suitable for neural computations

The set-up of our AMNs requires all reactions to be unidirectional, that is, the solutions must show positive-only fluxes (which is not guaranteed by usual genome-scale metabolic models). To split reactions of a given metabolic network into separate forward and reverse reactions, we wrote a standardization script that loads an SBML model into Cobrapy [START_REF] Ebrahim | COnstraints-Based Reconstruction and Analysis for Python[END_REF] and screens for all two-sided reactions, then duplicating them into two separate reactions; and writes a new version of the model with bidirectional reactions split into separate forward and backward reactions. To avoid confusion, we add a suffix for these reactions, either "for" or "rev" respectively designating the original forward reaction and the reversed reaction. The uptake reactions were also duplicated, even if encoded as one-sided, and their suffix was set to "i" for inflow reactions (addition of matter to the system), and "o" for outflow reactions (removal of matter from the system).

As detailed in the next subsection, our unidirectional models are used to build training sets. The duplicated iML1515 [START_REF] Monk | iML1515, a knowledgebase that computes Escherichia coli traits[END_REF]

Generation of training sets with FBA

Our reference flux data are obtained from FBA simulations, using the GNU Linear Programming Kit ('GLPK', a simplex-based method) on Cobrapy [START_REF] Ebrahim | COnstraints-Based Reconstruction and Analysis for Python[END_REF] , with different models of different sizes. Throughout this paper, when 'reference data' is mentioned, it refers to data computed with this method.

Reference data for metabolic flux distributions were generated using models downloaded from the BiGG database [START_REF] Norsigian | BiGG Models 2020: multi-strain genome-scale models and expansion across the phylogenetic tree[END_REF] . The models were used to generate data using Cobrapy [START_REF] Ebrahim | COnstraints-Based Reconstruction and Analysis for Python[END_REF] following a precise set of rules. First, we identified essential uptake reactions for the models we used (E. coli core 29 and iML1515 [START_REF] Monk | iML1515, a knowledgebase that computes Escherichia coli traits[END_REF] ). Precisely, if one of these reactions has its flux upper bound set to zero, the 'biomass' reaction optimization is impossible, even if all other uptake fluxes bounds are set to a high value, e.g., 1000. In other words, we identified the minimal uptake fluxes enabling growth according to the models. During reference data generation, the upper bounds on these reactions were set to 10. For E. coli core [START_REF] Orth | Reconstruction and Use of Microbial Metabolic Networks: the Core Escherichia coli Metabolic Model as an Educational Guide[END_REF] we found 7 of such obligate reactions (for the uptake of CO2, H+, H20, NH4, O2, Phosphate, and Glycerol as the carbon source). For iML1515 [START_REF] Monk | iML1515, a knowledgebase that computes Escherichia coli traits[END_REF] we had the same 7 obligate reactions and additional salt and ions uptake reactions (for the uptake of Fe2+, Fe3+, Mn2+, Zinc, Mg, Calcium, Ni2+, Cu2+, Selenate, Co2+, Molybdate, Sulfate, K+, Sodium, Chloride, Tungstate, Selenite). With iML1515 [START_REF] Monk | iML1515, a knowledgebase that computes Escherichia coli traits[END_REF] , we also added as obligate reactions the uptake of four amino acids (Alanine, Proline, Threonine and Glycine) in order to be consistent with our experimental training set where the four amino acids were systematically added to M9 (cf. next subsection).

To generate different media compositions, we added to the obligate reactions, a set of variable uptake reactions. For the E. coli core model [START_REF] Orth | Reconstruction and Use of Microbial Metabolic Networks: the Core Escherichia coli Metabolic Model as an Educational Guide[END_REF] we added 13 variable uptake reactions (for the uptake of Acetate, Acetaldehyde, Oxoglutarate, Ethanol, Formate, Fructose, Fumarate, Glutamine, Glutamate, Lactate, Malate, Pyruvate, and Succinate). For each generated medium, variable uptake reactions were selected following a binomial distribution B(n, p) with n=13 and p = 0.5, p being a tunable parameter related to the ratio of selected reaction. Consequently, the mean number of selected variable uptake reactions was n x p = 6.5. Next for each selected reaction, the upper bound value of the reaction flux was drawn randomly between 2 and 10. For the iML1515 [START_REF] Monk | iML1515, a knowledgebase that computes Escherichia coli traits[END_REF] model, to limit the combinatorial search space, the selected variable uptake reactions were those of the experimental training set and consequently between 1 and 4 variable uptake reaction were added (cf. next subsection). The upper bound values for each selected variable reaction were chosen randomly between 0 and 2.2 (0 excluded). The 2.2 threshold was chosen to produce predicted growth rates that were in the range of those observed experimentally.

After generating the set of media for E. coli core [START_REF] Orth | Reconstruction and Use of Microbial Metabolic Networks: the Core Escherichia coli Metabolic Model as an Educational Guide[END_REF] and iML1515 [START_REF] Monk | iML1515, a knowledgebase that computes Escherichia coli traits[END_REF] we ran Cobrapy [START_REF] Ebrahim | COnstraints-Based Reconstruction and Analysis for Python[END_REF] for each medium with FBA and recorded the steady state value for all fluxes including the growth rate (flux of the 'biomass' reaction). These rates were used as a training set for all models presented in Table 1. For all UB training sets, the variable uptake flux values were those used by Cobrapy [START_REF] Ebrahim | COnstraints-Based Reconstruction and Analysis for Python[END_REF]

Generation of an experimental training set

Ten carbon sources (Ribose, Maltose, Melibiose, Trehalose, Fructose, Galactose, Acetate, Lactate, Succinate, Pyruvate) were picked for being the variables of our training sets. These could ensure observable growth as a sole carbon source with a concentration of 0.4 g.L -1 in our M9 preparations.

The selected carbon sources enter different parts of the metabolic network: 6 sugars enter the upper glycolysis pathway, and 4 acids enter the lower glycolysis pathway or the TCA cycle. With a binary (i.e., presence or absence of each carbon source) approach when generating the combinations to test for making the experimental training set, we generated all possible combinations of 1, 2, 3 or 4 carbon sources simultaneously present in the medium. Naturally, we picked all 1-carbon source media combinations for experimental determination (only 10 points). Then, we randomly selected 100 more combinations to experimentally determine, by randomly picking 20 points from the 2-, 40 points from number of preserved replicates range from 2 to 8, with an average of 4.6 (± 1.6) replicates per medium composition. The mean and standard deviations over replicates were computed to be used for training AMNs and making figures. The code for processing raw data and all the raw data itself is available in the Github package (see Code availability and Data availability sections).

Loss functions derivation

Loss functions are necessary to assess the performances of all MM (MM-LP, -QP) and AMN architectures (AMN-Wt, -QP, and -LP) and also to compute the gradients of the QP solvers.

To perform the derivation, we assume we have a metabolic model with 𝑛 reactions and 𝑚 metabolites. Let 𝑉 = (𝑣 1 , … , 𝑣 𝑛 ) 𝑇 be the reaction flux vector and 𝑆 the 𝑚 × 𝑛 stoichiometric matrix of the model. We assume some metabolites can be imported in the model, with each a corresponding uptake reaction. Let 𝑉 𝑖𝑛 be the vector of 𝑛 𝑖𝑛 upper bounds (or exact values) for these uptake reactions, and let 𝑃 𝑖𝑛 the 𝑛 𝑖𝑛 × 𝑛 projection matrix such that 𝑉 𝑖𝑛 = 𝑃 𝑖𝑛 𝑉. We also assume some reaction fluxes have been experimentally measured, let 𝑉 𝑜𝑢𝑡 be the vector of measured values. With 𝑃 𝑜𝑢𝑡 the 𝑛 𝑜𝑢𝑡 × 𝑛 projection matrix for measured fluxes, 𝑉 is calculated solving the following quadratic program (QP):

𝑚𝑖𝑛( ‖𝑃 𝑜𝑢𝑡 𝑉 -𝑉 𝑜𝑢𝑡 ‖ 2 ) (1)
subjected to:

𝑆 𝑉 = 0 𝑃 𝑖𝑛 𝑉 ≤ 𝑉 𝑖𝑛 𝑉 ≥ 0
For each solution (𝑉) of eq. ( 1), four losses apply. The first is related to the fit to the reference targeted values and the three additional losses are related to the boundary, stoichiometric and flux positivity constraints of the metabolic network.

The first loss is simply the Mean Square Error (MSE) between predictions (𝑉) and measurements(𝑉 𝑜𝑢𝑡 ):

𝐿 1 = 1 𝑛 𝑜𝑢𝑡 ‖𝑃 𝑜𝑢𝑡 𝑉 -𝑉 𝑜𝑢𝑡 ‖ 2 (2)
The second loss is linked to the network stoichiometric constraint (𝑆 𝑉 = 0), which in its normalized form (loss per constraint) is:

𝐿 2 = 1 𝑚 ‖𝑆𝑉‖ 2 (3)
The third loss evaluates how well boundary constraints are respected (𝑃 𝑖𝑛 𝑉 ≤ 𝑉 𝑖𝑛 ):

𝐿 3 = 1 𝑛 𝑖𝑛 ‖𝑅𝑒𝐿𝑈( 𝑃 𝑖𝑛 𝑉 -𝑉 𝑖𝑛 )‖ 2 (4)
The last loss enforces all fluxes to be positive:

𝐿 4 = 1 𝑛 ‖𝑅𝑒𝐿𝑈(-𝑉)‖ 2 (5)
We note that when exact bound are provided, 𝑃 𝑖𝑛 𝑉 = 𝑉 𝑖𝑛 , and 𝐿 3 becomes obsolete as the values of 𝑉 corresponding to 𝑉 𝑖𝑛 are not updated by the LP/QP solvers and AMN programs.

The total loss (𝐿) is calculated as the mean sum of the previous losses:

𝐿 = (𝐿 1 + 𝐿 2 + 𝐿 3 + 𝐿 4 ) / 4𝑁 (6)
where 𝑁 is the sample batch size.

𝐿 can also be computed as the MSE between the vectors ) and (𝑉 𝑜𝑢𝑡 , 0, 0, 0)

QP solver

The QP solver solves the quadratic program given by eq. ( 1) in the Loss functions derivation subsection.

While the QP system can be solved by a simplex algorithm, solutions can also be approximated calculating 𝑉 minimizing the loss of eq. ( 6):

𝑚𝑖𝑛 ( 1 𝑛 𝑜𝑢𝑡 ‖𝑃 𝑜𝑢𝑡 𝑉 -𝑉 𝑜𝑢𝑡 ‖ 2 + 1 𝑚 ‖𝑆𝑉‖ 2 + 1 𝑛 𝑖𝑛 ‖𝑅𝑒𝐿𝑈( 𝑃 𝑖𝑛 𝑉 -𝑉 𝑖𝑛 )‖ 2 + 1 𝑛 ‖𝑅𝑒𝐿𝑈(-𝑉)‖ 2 ) (8)
vector 𝑉 can thus be found solving:

𝜕 ( 1 𝑛 𝑜𝑢𝑡 ‖𝑃 𝑜𝑢𝑡 𝑉-𝑉 𝑜𝑢𝑡 ‖ 2 + 1 𝑚 ‖𝑆𝑉‖ 2 + 1 𝑛 𝑖𝑛 ‖𝑅𝑒𝐿𝑈( 𝑃 𝑖𝑛 𝑉-𝑉 𝑖𝑛 )‖ 2 + 1 𝑛 ‖𝑅𝑒𝐿𝑈(-𝑉)‖ 2 ) 𝜕𝑉 = 0 (9)
𝑉 satisfying (3) can be found iteratively:

𝑉 (𝑡+1) = 𝑉 (𝑡) + 𝑑𝑡 𝛻𝑉 (10) 𝑉 (0) = 𝑃 𝑖𝑛 𝑇 𝑉 𝑖𝑛
where 𝑡 is the iteration number, 𝑑𝑡 the learning rate, and:

𝛻𝑉 =
where When exact bounds are provided (EB case), the gradient becomes:

𝛻𝑉 = ( 2 𝑛 𝑜𝑢𝑡 𝑃 𝑜𝑢𝑡 𝑇 (𝑃 𝑜𝑢𝑡 𝑉 -𝑉 𝑜𝑢𝑡 ) + 2 𝑚 𝑆 𝑇 𝑆𝑉 - 2 𝑛 𝐷 𝑉 𝑅𝑒𝐿𝑈(-𝑉)) ⨀ (1 𝑛×1 -𝑃 𝑖𝑛 𝑇 1 𝑛 𝑖𝑛 ×1 ) (12)
where ⨀ stands for Hadamard product and 1 𝑛×1 (1 𝑛 𝑖𝑛 ×1 ) is a constant (=1) vector of dimension

𝑛 (𝑛 𝑖𝑛 )
Details on gradient calculations can be found in the Supplementary Information -section QP solver equations.

LP solver

When uptake fluxes are known (EB method), the linear optimization problem solved with FBA for a metabolic network with n fluxes and m metabolites can be written as:

𝑀𝑎𝑥: 𝑐 𝑇 𝑉 (13)
𝑠𝑡: 𝑆𝑉 = 𝑏 𝑉 ≥ 0 and its dual form:

𝑀𝑖𝑛: 𝑏 𝑇 𝑀 𝑠𝑡: 𝑆 𝑇 𝑀 ≤ 𝑐 𝑇 (14)
In the case where only uptake fluxes upper bounds are known (UB method) it can be written as:

𝑀𝑎𝑥: 𝑐 𝑇 𝑉 𝑠𝑡: 𝑆 𝑖𝑞 𝑉 ≤ 𝑏 𝑎𝑛𝑑 𝑆 𝑒𝑞 𝑉 = 0 (15) 𝑉 ≥ 0 its dual form being eq. ( 14).

In eqs. (13-15), V is the vector of fluxes, M its dual, referred to as metabolite shadow prices, S the stoichiometric matrix, b a vector of dimension m with bi corresponding to input fluxes of medium metabolites and c the objective vector of dimension n. In the UB case, Siq is the stoichiometric matrix involving only internal reactions and Seq is the stoichiometric matrix with uptake reactions. Both matrices that can be derived from S are given in Figs. S3 and S4 in Supplementary Information. Keeping in mind that b takes positive values only for metabolites transported by uptake reactions, consequently bi = 0 for all internal metabolites. Using the notation of Fig. 1 we note that 𝑏 = 𝑃 𝑚→𝑣 𝑉 0 , where V 0 is the vector of uptake fluxes. In this work c is a null vector with 1 for the 'biomass' reaction flux.

The usual method to solve this linear programming problem is the Simplex algorithm [START_REF] Karloff | The Simplex Algorithm[END_REF] , which is fast and scalable. But simplex has no gradient, nor backpropagation compatibility. Thus, it cannot be integrated in AMNs and in any training procedure.

The optimization-RNN proposed by Yang et al. [START_REF] Yang | A new neural network for solving quadratic programming problems with equality and inequality constraints[END_REF] is an alternative to solve LP problems, as this method is compatible with backpropagation. The general solving happens by iteratively updating V the vector of all reaction fluxes and M the dual vector of shadow prices using the following gradients: Further details on this derivation and in particular the form of the matrices 𝑆 𝑒𝑞 , 𝑆 𝑖𝑞 and vectors 𝑏, 𝑑 in the EB and UB cases can be found in the Supplementary Information -section LP-solver equations.

The method converges for both EB and UB cases to the global optimum independently of the initialization. The learning rate was arbitrarily fixed to 0.3 for both EB and UB.

For both Fig. 2d and Fig. 2e the results were obtained with 100 simulations and 10 6 iterations. The reference values were obtained using Cobrapy [START_REF] Ebrahim | COnstraints-Based Reconstruction and Analysis for Python[END_REF] and the E. coli core model [START_REF] Orth | Reconstruction and Use of Microbial Metabolic Networks: the Core Escherichia coli Metabolic Model as an Educational Guide[END_REF] . Each replicate was generated with a different medium computed with the training set generation method (see Method section 'Training Set Generation').

AMN architectures and parameters

As listed in 1e and detailed in Supplementary Information -AMN-Wt architecture. The architectures of AMN-LP/QP are shown in Fig. 3a. All ANMs take as input a vector of fluxes of size 𝑛 𝑖𝑛 for medium uptake fluxes (Vin), then transforms via a dense neural network the input vector into an initial vector of size n for all fluxes (V 0 ), which is refined through an iterative procedure computing 𝑉 (𝑡+1) from 𝑉 (𝑡) . With all AMNs a 𝑛 𝑖𝑛 x 𝑛 weight matrix transforming Vin to V 0 is learned during training, and we name this transforming layer the neural layer. With AMN-LP/QP, 𝑉 (𝑡) is iteratively updated in a mechanistic layer by the gradient (𝛻𝑉) of LP/QP solvers (cf. previous subsections in Method). With AMN-Wt, the mechanistic layer computes 𝑉 (𝑡+1) from 𝑉 (𝑡) using the transformations shown in Supplementary Information - For all architecture, we use the Mean Squared Error (MSE, eq. ( 2)) on all fluxes as the objective function to minimize while learning. In all AMN architectures we add to the MSE loss function, the terms corresponding to the 3 losses derived from the constraints of the metabolic model (eqs. (3-5)).

The parameters used when training AMNs can be decomposed into two categories:

1. Simulated data parameters. As described in the previous section (Generation of Training Sets with FBA), we can tune the size of the training set to be generated. We can also modify the mean number of selected variable intake medium fluxes, and the number of levels (i.e., the resolution) of the fluxes. We can also modify the actual variable intake medium flux list, but this modifies the architecture of the model (initial layer size), so we kept the same list for each model in the present work.

2.

Learning parameters. During learning on simulated data, an AMN has a small set of parameters to tune: the number and size of hidden layers, the number of epochs, the batch size, the dropout ratio and the number of folds in cross-validation.

Optimization of uptake fluxes upper bounds in FBA

The goal of this optimization was to find the best scaler for fluxes to best match experimentally determined growth rates, by using 'out-of-the-box' FBA, simply informing the presence or absence of the flux according to the experimental medium composition. The optimal scaler was found using the Cobrapy software package [START_REF] Ebrahim | COnstraints-Based Reconstruction and Analysis for Python[END_REF] by simply searching for the maximum R² between experimental and FBApredicted growth rates. The plot showing how the R² behaves as a function of the scaler is shown in Fig. 5f.

Fig. 1 .

 1 Fig. 1. Computing steady-state fluxes with stoichiometric and neural models. a. Simple toy stoichiometric model. The model is unidirectional, all flux values are positive. b. Steady-state solution fluxes maximizing 𝑣 3 . At steady state, the reaction fluxes (𝑣 𝑖 ) must satisfy stationarity conditions that guarantee mass balance of all metabolites, this is depicted by the equation 𝑆𝑉 = 0, where 𝑆 is the stoichiometric matrix representing the connectivity of the model and 𝑉 the vector of fluxes to be calculated.𝑉 𝑖𝑛 is the medium represented by a vector of nutrient uptake fluxes (here 𝑣 1 = 0.1, symbol "-" indicates that no value is provided, in practice we use an 'infinity' value to represent an unbounded flux). The steady-state solution 𝑉 𝑜𝑢𝑡 is calculated by solving a linear program maximizing the objective 𝑐 𝑇 𝑉 = 𝑣 3 , here we used the Cobrapy package[START_REF] Ebrahim | COnstraints-Based Reconstruction and Analysis for Python[END_REF] to compute 𝑉 𝑜𝑢𝑡 (simplex-based algorithm). c. Stoichiometric model matrices. 𝑃 𝑣→𝑚 is the matrix to compute metabolite production fluxes from reaction fluxes, 𝑃 𝑚→𝑣 is a matrix to compute reaction fluxes from the production fluxes of the reaction substrates. When a metabolite is the substrate of several reactions, each reaction gets a fraction of the metabolite production flux, this is depicted in matrix 𝑊 𝑟 (r indicates this matrix is used in recurrence). 𝑀 𝑜𝑢𝑡 is the vector of metabolite production fluxes at steady state, the operator ⊙ stands for element-wise matrix product (Hadamard product), and 𝑅𝑒𝐿𝑈(𝑥) = 𝑚𝑎𝑥 (0, 𝑥). sj,i corresponds to the value of S at the j th row (metabolite) and i th column (flux) and zi is the number of negative elements in column i of S. d. Unrolled neural network built from the model. In this theoretical example an initial flux vector is set, representing only the uptake fluxes, and is mapped to a flux vector covering the entire network. In the initial layer ( 𝑙 0 ) only 𝑣 1 has a value (𝑣 1 0 ). In layer 1, 𝑣 1 value is passed to 𝑚 1 , the production flux for metabolite 1, 𝑚 1 1 = 𝑣 1 0 . Subsequently a fraction (𝑤 21 ) of 𝑚 1 goes to 𝑣 2 and the other fraction (𝑤 31 ) to 𝑣 3 , therefore 𝑣 2 1 = 𝑤 21 𝑚 1 1 , 𝑣 3 1 = 𝑤 31 𝑚 1 1 with 𝑤 21 +𝑤 31 = 1. Additionally, 𝑣 1 remains as in 𝑙 0 : 𝑣 1 1 = 𝑣 1 0 . In layer 2, 𝑣 1 continues to feed 𝑚 1 , 𝑣 2 is passed on to 𝑚 2 and 𝑚 3 , therefore, 𝑚 2

Fig. 2 .

 2 Fig. 2. FBA solver architectures and performances. a. Schematic procedure for the Simplex solvers. From Vin, which is a vector describing the bounds of some uptake fluxes, the solvers reach a steady-state solution (Vout) optimizing the objective function c, satisfying the constraints and bounds of the network (cf. Fig. 1 panel b). One of the most popular solvers is a simplex-based method in Cobrapy 22 , taken as reference data for the results shown here. b. Schematic for LP-solver architecture. This in-house solver surrogates the simplex-based algorithm. Following Yang et al. 23 and as further detailed in the Method section (LP-solver Architecture), the full flux distribution V is updated by ∇V and the metabolites shadow prices M by ∇M through products of matrices derived from the stoichiometry of the network. c. Schematic for QP-solver architecture. Here target reference fluxes are given to the solver. The computed fluxes are fitted to the Reference Targets by the means of a custom loss function integrating also the input constraints along with the stoichiometric constraint of the metabolic network. The flux vector V is updated by ∇V which is the gradient minimizing the loss function (cf. Method section -QP Architecture -for further details). d. Matching of growth rate between the Simplex solver and the LP/QP solvers. EB refers to "Exact Bound", a set-up in which Vin is composed of the exact flux values for each uptake flux. UB refers to "Upper Bound", in which Vin is composed of upper threshold values for each uptake flux. e. Loss vs. Solver iterations. QP takes 1 million iterations to reach close to zero values, whereas LP takes 10,000 iterations. Plotted is the mean and standard error (95% confidence interval) across all elements of the set of 100 simulations.

Fig. 3 :

 3 Fig. 3: AMN architecture enabling training. a. Basic AMN architecture principle. Vin is the vector representing upper or exact bounds on uptake fluxes. Each white or blue circle represents one reaction. V 0 contains the initial fluxes predicted by a dense neural network. Vout is the final flux vector, computed by the LP or QP solvers. b. Loss vs. Learning Epochs number for E. coli core model 29 . The architecture of panel a (with 4 iterations of the QP solver as the mechanistic layer, with upper bounds) was trained for 50 epochs, measuring the loss on Vout after each epoch. Plotted is the loss mean and standard error (95% confidence interval) over a training set of 1000 examples generated with the E. coli core model 29 (see Methods -Generation of training sets with FBA). The different curves correspond to different hidden layer sizes used between Vin and V 0 .

  model is quite large, comprising 3682 reactions and 1877 metabolites. A substantial number of reactions in that model have a zero flux for many different media and it is unnecessary to keep these reactions during the training process. Prior training, we therefore generate a reduced model by removing zero flux reactions along with the metabolites no longer contributing to any reactions. Using that procedure, we were able to reduce iML1515 20 model to only 550 reactions and 1083 metabolites.

  to generate the training set. For EB training sets, the variable uptake flux values were those calculated by Cobrapy 22 at steady state.

  𝐷 𝑖𝑛 = 𝑅𝑒𝐿𝑈(𝑃 𝑖𝑛 𝑉-𝑉 𝑖𝑛 ) 𝑅𝑒𝐿𝑈(𝑃 𝑖𝑛 𝑉-𝑉 𝑖𝑛 ) and 𝐷 𝑉 = 𝑅𝑒𝐿𝑈(-𝑉) 𝑅𝑒𝐿𝑈(-𝑉) using an Hadamard division:

  𝛻𝑉 = ((𝐼 -𝑃)[ 𝑐 𝐹𝐵𝐴 -𝑆 𝑖𝑞 𝑇 𝑅 𝑖𝑞 ] -𝑄𝑅 𝑖𝑛𝑡 ) 𝑅 𝑖𝑞 = 𝑅𝑒𝐿𝑈(𝑀 + 𝑆 𝑖𝑞 𝑉 -𝑑) 𝑅 𝑒𝑞 = 𝑆 𝑒𝑞 𝑉 -𝑏 𝑄

  Fig. S1, which include a 𝑛 x 𝑚 weight matrix (𝑊 𝑟 ). That weight matrix is learned during training. With all AMN architectures, the values of V corresponding to Vin are not updated in the neural or mechanistic layers when training with exact values for medium uptake (EB training sets).

  𝑤 42 = 1, other fluxes remain the same as in 𝑙[START_REF] Jumper | Highly accurate protein structure prediction with AlphaFold[END_REF] . In layer 3, 𝑚 4 receives input from 𝑣 4 which in turn activates 𝑣 5 . In layer 4, 𝑚 1 receives input from both 𝑣 1 and 𝑣 5 : 𝑚1 4 = 𝑣 1 3 + 𝑣 5[START_REF] Thornburg | Fundamental behaviors emerge from simulations of a living minimal cell[END_REF] . e. Recurrent neural network representation. At each iteration step 𝑉 𝑙 and 𝑀 𝑙 are computed using matrices 𝑃 𝑣→𝑚 and 𝑃 𝑚→𝑣 of panel c, while the matrix 𝑊 𝑟 can be learned by training with experimental data or model simulations. For example, setting 𝑣 1 0 = 0.1 and searching weights for which 𝑣 3 𝑛 = 0.5 one finds 𝑤 21 = 0.74, 𝑤 31 = 0.26, 𝑤 32 = 0.20, and 𝑤 42 = 0.80. Taking these weights, the 𝑉 𝑛 values obtained for n=30 match those of panel B. f. Heatmap for flux values up to 30 iterations.

	3 2 = 𝑤 32 𝑚 2 2 , 𝑣 4 2 = 𝑤 42 𝑚 2 2 with 𝑤 32 +

1 1 = 𝑣 1 0 . Subsequently a fraction (𝑤 21 ) of 𝑚 1 goes to 𝑣 2 and the other fraction (𝑤 31 ) to 𝑣 3 , therefore 𝑣 2 1 = 𝑤 21 𝑚 1 1 , 𝑣 3 1 = 𝑤 31 𝑚 1 1 with 𝑤 21 +𝑤 31 = 1. Additionally, 𝑣 1 remains as in 𝑙 0 : 𝑣 1 1 = 𝑣 1 0 . In layer 2, 𝑣 1 continues to feed 𝑚 1 , 𝑣 2 is passed on to 𝑚 2 and 𝑚 3 , therefore, 𝑚 2 2 = 𝑣 2 1 and 𝑚 3 2 = 𝑣 2 1 , 𝑚 2 then goes to 𝑣 3 and 𝑣 4 , and as in the previous layer we have 𝑣

Table 1 ,

 1 we propose three AMNs architectures: AMN-Wt, AMN-LP and AMN-QP. The AMNs are run with training sets with exact medium values given at steady state (EB) or only upper bound values (UB). The architecture of RNN-Wt is shown in Fig.

Acknowledgements

JLF would like to acknowledge funding provided by the ANR funding agency grant numbers ANR-18-CE44-0015 (SynBioDiag project) and ANR-21-CE45-0021-01 (AMN project) and the BIOS UE HORIZON program (grant number 101070281). LF is supported by INRAE's MICA department and by INRAE's metaprogram DIGIT-BIO. BM is supported by an Ecole Normale Supérieure Scholarship. We thank Aymeric Gaudin (CentraleSupélec Engineering School) for early development in reservoir computing with AMN, Ivan Radkevich (University of Paris Saclay) for his work on custom RNN cells and Tom Lorthios and Hadi Jbara (AgroParisTech and University of Paris Saclay) for their help on collecting data for experimental training sets.

Data availability

All raw and processed experimental data can be found in our Github repository at https://github.com/brsynth/amn_release/tree/main/Dataset_experimental. Some supplementary files (including 'Data_Fig4.xlsx' and 'Data_Fig5.xlsx') and all other datasets used to produce figures of this study can be also found in our Github repository, at https://github.com/brsynth/amn_release/tree/main/Result.

Code availability

All codes are available within a documented GitHub repository at https://github.com/brsynth/amn_release. The GitHub repository includes tutorials in Google Colab notebooks. All AMN codes make use of Cobrapy 22 , numpy [START_REF] Harris | Array programming with NumPy[END_REF] , scipy [START_REF] Virtanen | SciPy 1.0: fundamental algorithms for scientific computing in Python[END_REF] , Pandas 38 , tensorflow [START_REF] Abadi | TensorFlow: Large-Scale Machine Learning on Heterogeneous Distributed Systems[END_REF] , sci-kit learn [START_REF] Pedregosa | Scikit-learn: Machine learning in Python[END_REF] and keras 41 libraries.

the 3-and 40 points from the 4-carbon source combinations sets. The python scripts to generate these combinations and pick the ones for making our experimental training set are available in our Github package (see Codes availability section). After picking the combinations to test, we experimentally determined the maximum specific growth rate of E. coli for each combination of carbon sources in M9 (see next two subsections). The mean over replicates for each media composition was computed as the corresponding growth rate value to make the final experimental training set (cf. Method -Growth rate determination).

Culture conditions

The base medium for culturing E. coli DH5-α (DH5a) was a M9 medium prepared with those final concentrations: 100µM CaCl2, 2mM MgSO4, 1X M9 salts (3 g.L -1 KH2PO4, 8.5 g.L -1 Na2HPO4 2H2O, 0.5 g.L -1 NaCl, 1g.L NH4Cl), 1X trace elements (15 mg.L -1 Na2EDTA 2H20, 4.5 mg.L -1 ZnSO4 7H2O, 0.3 mg.L - 1 CoCl2 6H2O, 1 mg.L -1 MnCl2 4H2O, 1 mg.L -1 H3BO3, 0.4mg.L -1 Na2MoO4 2H20, 3 mg.L -1 FeSO4 7H2O, 0.3 mg.L -1 CuSO4 5H2O; solution adjusted at pH=4 and stored at 4°C), 1 mg.L -1 Thiamine-HCl and 0.04g.L -1 amino acid mix so that each amino acid (L-Alanine, L-Proline, L-Threonine, Glycine) was at a final concentration of 5 mg.L -1 in the medium. The additional carbon sources that could be added were individually set to a final concentration of 0.4 g.L -1 . The pH was adjusted at 7.4 prior to a 0.22µm filter sterilization of the medium. Pre-cultures were recovered from glycerol -80°C stocks, grew in Luria-Bertani (LB) broth overday, then used as inoculate in 200µL M9 (supplemented with variable compounds) in 96 U-bottom wells plates overnight, then passed to a replicate of the plate on the next day. The temperature was set to 37°C in a plate reader (BioTek HTX Synergy), with continuous orbital shaking, allowing aerobic growth for 24 hours. A monitoring every 10 minutes of the optical density at 600 nm (OD600) was performed. A figure for summarizing the experimental workflow is available in Supplementary Fig. S6.

Growth rate determination

The maximal growth rate was determined by sliding a window of 1 hour-size, performing a linear regression on the log(OD600) data in each window. We then retrieve the maximum specific growth rate as the maximum regression coefficient over all windows. If several growth phases are visible, one can omit a part of the growth curve for the maximal growth rate determination (for this study we always retrieved the maximal growth rate on the first growth phase, so as we are certain that the media contains all added carbon sources). 8 replicates for each medium composition were performed (on a single column of a 96-well plate). Outliers were manually removed after visual inspection of the growth curves or clear statistical deviation of the computed growth rate from the remaining replicates. The

Author contributions

LF and JLF wrote the core of the text of the manuscript. JLF designed the study and wrote the QP solver of Fig. 2 and all the AMN and RC codes used to produce data for Figs. 1, 3, 4, 5 and Table 1. BM wrote the LP-solver producing Figs. 2, S4 and S5 and wrote the corresponding part in the method section. LF benchmarked all codes, wrote the codes transforming SBML models into unidirectional networks and processing experimental data, and handled the colab and git implementations. LF also performed all experimental work reported in Figs. 4 and5 and wrote the corresponding experimental method section. WL contributed to designing the project and was involved in the discussions and writing the manuscript. All authors read, edited, and approved the manuscript.