Valentin Pannetier
email: valentin.pannetier@math.u-bordeaux.fr

Michael Leguèbe

Yves Coudière

Richard D Walton

Philippe Dhiver

Delphine Feuerstein

Diego Amaro

Modeling Cardiac Stimulation by a Pacemaker, with Accurate Tissue-Electrode Interface

de niveau recherche, publiés ou non, émanant des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

Introduction

An implantable pacemaker aims to restore a cardiac beat when the intrinsic conduction system fails. It sends energy to the heart in the form of a voltage pulse for a certain duration via pacing leads implanted inside the heart. The device is programmed to deliver enough energy to trigger a cardiac depolarization (which is called capture). For battery saving reason, the energy must be just above the cut-off threshold between capturing and non-capturing regions of the voltage -duration plane, known as the Lapicque curve [START_REF] Blair | On the intensity-time relations for stimulation by electric currents[END_REF]. We aim to reproduce by computer modeling and simulation the experimental tests run to identify this curve. This would ultimately enable manufacturers to try out and optimize the delivered energy of several leads with different electrode shapes in-silico, instead of prototyping, manufacturing and testing leads experimentally.

We propose to couple a computational model of cardiac excitable cells (tissue) in an electrolyte (blood) to a model of the pacemaker's circuitry, through interfaces located at the metal electrodes of the leads. Modeling these interfaces by themselves is an extensive field of research, because many complex biophysical phenomena occur when an electrode delivers current into an electrolyte, as compiled in [START_REF] Grimnes | Bioimpedance and Bioelectricity Basics, 3rd edn[END_REF], and because a large part of the energy is dissipated by the contacts. It is then crucial to model correctly these bio-electrode interfaces, in order to reproduce quantitatively the threshold voltages. Cardiac simulations generally consider only the myocardium, and apply an artificial, irrealistic stimulation. Instead, we introduce boundary conditions that include impedance of the contacts, and act as current source of the cardiac model, following the approach of Somersalo et al. [START_REF] Somersalo | Existence and uniqueness for electrode models for electric current computed tomography[END_REF].

In this paper, Sects. 2, 3, and 4, describe the complete modeling approach, model calibration with respect to experimental data, and comparison of computed Lapicque curves to experimental ones.

Pacemaker Model Coupled to Cardiac Ionic Equations

Pacemaker Circuit. We model the electronics of the pacemaker with the circuit shown in Fig. 1 (left), which appeared to be standard among manufacturers. We focus on the short time interval dedicated to pacing, ignoring the sensing ability of the device, as it is not linked directly to the triggering of a cardiac action potential. A tank capacitor C pulse (∼ 10 µF) is charged (not modelled here), and then the actual pacing follows the three steps below (see Fig. 2).

Pulse. When the device triggers a stimulation, the tank capacitor acts as source of current. In the meantime, a secondary capacitance C ocd is charged. The duration of this step is programmed by the clinician, and can range from 0.25 to 2 ms. The amplitude of the stimulation, proportional to the charge of the tank capacitor, is also programmable, typically between 0.25 V and 2 V. Switch. There is a short transition phase (in the order of hundreds of µs) during which the two switches are open (see Fig. 1, left, black part). Output Capacitor Discharge (OCD). The capacitor C pulse is isolated from the circuit, and C ocd discharges. This step lasts for 13 ms and is intended to discharge the equivalent capacity due to the polarization of the electrodes that occurs during the pulse (see below).

After the OCD, the stimulating part of the device is switched off, and is replaced by a sensing circuit until the next stimulation. In the meantime, the tank capacitor C pulse is recharged. During the pulse and OCD steps, the circuit can be simplified into an equivalent RC series circuit. In the following, we simply denote by R and C the equivalent resistance and capacitance of this simplified circuit, specifically, 1 C = 1 C pulse + 1 C ocd during the pulse and 1 C = 1 C ocd during OCD, while R = R pulse + R gnd during the pulse and R = R gnd during OCD.

Modeling the Bio-Electrode Contacts. Multiple electrochemical interactions occur in the vicinity of the surface of an electrode placed in an electrolyte when a current is applied, as described extensively in [START_REF] Grimnes | Bioimpedance and Bioelectricity Basics, 3rd edn[END_REF]. These interactions are modeled by equivalent electric circuits, called contact models. There exists a large collection of contact models, as reviewed in [START_REF] Grimnes | Bioimpedance and Bioelectricity Basics, 3rd edn[END_REF]. We choose to model the tip and ring contacts by a resistance in parallel with a capacitance (parallel R-C, see Fig. 1, green part). This model is a good compromise that can reproduce experimental measures while relying on a limited number of parameters that can, in consequence, be identified. For instance, a simple resistive-only model of contact cannot reproduce the charging behavior of the bio-electrode contact, either during the pacing step nor during the switch step. On the other hand, the parameters of a complex Cole-Cole model cannot be identified from our data.

We state the equations of the contacts and pacemaker model, during pulse and OCD, as follows:

dI dt + 1 τ I = 1 R tot V 1 τ 1 + V 2 τ 2 , ( 1a
)
dV i dt + V i τ i = I C i , i = 1, 2 (1b)
where the current I, and voltages V 2 and V 1 are such as represented on Fig. 1.

The other constants are

τ 1 = R 1 C 1 , τ 2 = R 2 C 2 , τ = R tot C tot , with R tot = R + R w1 + R w2 + R bath and 1 Ctot = 1 C + 1 C1 + 1 C2 .
Experimentally, the voltage

V = V 1 + V 2 + (R bath + R w1 + R w2 )I is measured
, where R bath is the equivalent resistances of the saline solution.

In absence of applied current, the two contacts naturally discharge. Hence, we assume that V 1 (t = 0) = V 2 (t = 0) = 0. Meanwhile, the main capacitor is charged to its nominal capacity, corresponding to a voltage V stim . In consequence, a current I(t = 0) = -V stim /R tot initially flows out of it. This voltage V stim is referred to as the amplitude of stimulation.

Pacemaker and Cardiac Tissue -Surrogate Model. To account for the excitability of the cardiac tissue, the resistance R bath is not sufficient. It is replaced by a cardiac cell ionic model. The ionic model is adapted to account for the propagation of the electric current in the extracellular medium, of equivalent conductance g e . The membrane is modeled by a capacitance C m in parallel with the total ionic current I ion . The intracellular path of conduction through gap junctions between cells is not taken into account in this simplification. We use the Beeler-Reuter ionic model [START_REF] Beeler | Reconstruction of the action potential of ventricular myocardial fibres[END_REF], and we denote the transmembrane voltage by V m .

As the Beeler-Reuter model is written in units per cm 2 of cell membrane, we introduce a scaling parameter S, representing a given surface of cell membrane. We assume that the total current I spreads uniformly on this surface. We also introduce the total conductance g e of the extracardiac medium between the electrodes. Assuming that a volume of 1 mm 3 is relevant to this approximation, we obtain orders of magnitude of 100 cm 2 for S and 10 -2 mS for g e . However, the excited medium consists of blood and cardiac tissue, which is highly heterogeneous. Therefore, these estimates can only be seen as a starting point for any calibration procedure, and the current model as a surrogate to a complete 3D bidomain model. The coupled 0D model reads as follow

dI dt + 1 τ I = 1 R tot V 1 τ 1 + V 2 τ 2 + I ion (h, V m ) C m + V m τ m , ( 2a
)
dV i dt + V i τ i = I C i , i = 1, 2 (2b)
dV m dt + I ion (h, V m ) C m + V m τ m = I SC m , ( 2c
) dh dt + g(h, V m ) = 0, (2d)
with τ , τ 1 , τ 2 as above, but with

R tot = R + R w1 + R w2 and 1 Ctot = 1 C + 1 C1 + 1 C2 + 1
SCm , and τ m := SCm ge . The ionic current I ion , state variables h and evolution function g are given by the Beeler-Reuter model. The voltage measured in experiments is now

V = V 1 + V 2 + (R w1 + R w2 )I + V m .
In addition to the parameters S and g e , we need to determine the parameters R 1 , τ 1 and R 2 , τ 2 that characterize each electrode.

Calibration from Bench Experiments

In order to evaluate the accuracy of the parallel R-C model for the electrode polarization, we measure the current I and the voltage V delivered by pacemaker leads directly in a saline solution, without cardiac tissue (Fig. 3, left). We used a Microport CRM Borea DR pacemaker, with a Vega lead in a solution with measured conductivity 2.85 mS cm -1 , to obtain 9 datasets of current and voltage (I meas and V meas ). The datasets differ exclusively by the amplitude (1, 2 or 4 V) and duration (0.25, 0.5 or 1 ms) of the delivered pulse. The experimental setup can be modelled by the circuit shown in Figs. 1 (left andmiddle), where we used the equivalent circuit from the previous section for the pacemaker, and R bath denotes the resistance of the saline medium. The corresponding differential equations system (1) has 5 parameters to be identified: R bath and the properties of the contacts R 1 , R 2 , τ 1 and τ 2 . The parameter R bath was deduced from the measured initial conditions, as the electrodes are considered depolarized: R bath = Vmeas(t=0) Imeas(t=0) -R w2 -R w1 . The four remaining parameters θ = (R 1 , R 2 , τ 1 , τ 2 ) were calibrated using the following cost function:

J(θ) = s ω s 1 s (t) V (θ) -V meas 2 V meas 2 + I(θ) -I meas 2 I meas 2 ,
where ω is a vector of three weights of sum 1 that allows to give more importance on a specific step, s ∈ {pulse, switch, OCD}. The function 1 s restricts the computation of the relative difference of the 2 norms of voltage and current on a specific step, V meas and I meas are the experimental data of the voltage and current, and V (θ) and I(θ) are the numerical solutions of the voltage and current, respectively. We used the SciPy implementation of the L-BFGS-B algorithm to minimize the cost function J within a subset of R 4 for θ. Each evaluation of J requires to solve system (1) that is linear in I, V 1 and V 2 . We solve it semi-analytically using eigenvalue decomposition at each time step in a dedicated Python code. The convexity of J was not established. Hence, we first optimize the logarithm of each parameter (R ∈ [1, 10 5 ] Ω for the resistances, and τ ∈ [10 -6 , 1] s for the characteristic times), in order to fix a correct order of magnitude, and then fit the values of the parameters. First, the model was calibrated using the 9 data files separately. For each file, and for a given choice of Ω, the L-BFGS-B algorithm was started at a value of the parameters θ for which the cost function J(θ) is the lowest among 10 000 random samples of θ. Then, we minimized J using the 9 files all at once, with the average value of the previously found parameters as starting point. We set ω = (0.45, 0.45, 0.1), so as to emphasize the steps in which the information can be more easily recovered.

The results of the fits are illustrated on Fig. 4, for both voltage and current. We found that one of the couples of parameters is well determined: R 2 = 26.0 ± 1.77 Ω, τ 2 = 0.0548 ± 6.93 10 -3 ms (C 2 = 2.12 ± 0.34 µF). However, we found values of R 1 with much larger variability across the files: between 2 and 6 kΩ, with average ∼ 4 kΩ, which is one order of magnitude larger than the total impedance measured by the device (940 Ω). We believe that it is due to the oversimplification of our contact model, leading to this equivalent resistance having no real physical meaning. We tried to use more complex models to reproduce the data during the OCD step, to overcome this limitation of the simple R-C parallel model. However, we faced problems of identifiability of the parameters during the calibration (ie, finding different local minima for repeated computations).

Anyway, the currents produced by the parallel R-C model are close to the data, and consequently, we consider this model to be accurate enough as a current source for the cardiac model. An example of simulated transmembrane voltage is given on Fig. 5 (left), for 5 stimulations of 428.5 mV and 1 ms, with S = 50 cm 2 and g e = 0.01 mS at 90 bpm. This set of parameters results in several captures and non-captures, due to variables of the Beeler-Reuter model that do not return to their original state between each pulse. This shows that a single stimulation with these parameters is not sufficient to trigger a complete depolarization.

From transmembrane voltages, we extract the ratio of successful captures (60% in the illustrated case). We can compute this ratio for several amplitudes V stim and durations d of stimulation, following a 2D dichotomy algorithm in the Lapicque space (Fig. 5 right). The algorithm detects when the capture ratio exceeds 50%, and returns a collection of intervals that should contain the Lapicque curve. We then compute the rheobase V rh (V) and chronaxie T ch (S) of Lapcique's law V stim (d) = V rh (1 + T ch /d) that best match the midpoints of these intervals. Our model seems in agreement with this law, for the set of parameters used to generate Fig. 5.

We then computed several Lapicque curves, for other values of S and g e , and compared the results with experimental data described in the next paragraph. Ex-vivo Measurements. Isolated coronary-perfused ventricular wedge preparations were carried out on two sheep hearts, as described in [START_REF] Walton | Compartmentalized structure of the moderator band provides a unique substrate for macroreentrant ventricular tachycardia[END_REF]. Perfusion leaks along cut surfaces were carefully occluded for homogeneous perfusion resistance across the preparation. The wedge preparation was then stretched on to a frame to immobilize the tissue into a bath of saline solution (Fig. 3, right).

To study the MicroPort lead (contacts characteristics and voltages applied to the tissue), three lead implantation regions were identified (RV apex, RV septum and RV basal sites) and the lead was implanted in the heart by fully deploying the screw electrode perpendicularly to the tissue's surface. Electrocardiography (ECG) and electrograms (EGM) measured with the lead were simultaneously recorded through an independent acquisition system (PowerLab, ADInstruments). The EGMs were also measured with the pacemaker.

We tested a range of pulse durations and stimulus amplitudes to localize the stimulation threshold. Due to the limited number of steps in stimulation amplitude and duration of the pacemaker, the measured data points cover only sparsely the Lapicque plane. Nevertheless, these data can later be used for fitting the computational model responses to lead stimuli.

Comparison Between Experiment and Simulation. From experimental data, we identified, for each pacing pulse duration, the minimum pulse amplitude that captured. Below this pulse amplitude, any stimulation fails to capture the tissue. We could therefore define a region, delimited on the top by the lowest capturing amplitude, and on the bottom by the highest non-capturing amplitude, in which the Lapicque curve should be found (Fig. 6).

We show from the 0D model that the parameter S has great influence, since it determines the amount of current seen by square centimeter of membrane. On the contrary, we found very little influence of the parameter g e on the output (not shown). We also studied the influence of parameters of the bio-electrode contact, by varying the resistance R 1 . This is justified by the fact that the tip electrode is not any more in contact with a saline solution only, but with the cardiac tissue too (Fig. 6 right). The fits of Lapicque law V stim = V rh (1 + T ch /d) do not fall within these regions for all durations. This is expected because the parameters S ad g e of the 0D model cannot be completely characterized, as surrogates of the geometry and electric conductivity coefficients of the real tissue. Additionally, these parameters may depend on the amplitude of the stimulation, as the excited volume of tissue may change. As a consequence, further calibration work is required, especially to determine the parameter S which plays a crucial role in this 0D model.

Conclusion

In this paper, an innovative model is proposed, which couples a pacemaker electrical circuit to equations modeling cardiac action potential, accounting for the complex physics at the bio-electrode interfaces. Animal experiments were conducted which provide information on the electrical function of a cardiac tissue sample stimulated by a commercial pacemaker. It allows us to compare Lapicque curves obtained with the numerical model with experimental capture data.

It is noticeable that the literature provides no model for coupling a pacemaker to a cardiac tissue, because all studies focus either on characterizing in-depth the physics of the bio-electrode contact, or on characterizing the onset of cardiac activation in the tissue sample (virtual electrodes, make and break activations, etc.). To our knowledge, this work provides a first possible model that links the energy delivered by the pacemaker to the onset of activation, trying to evaluate with a model the real distribution of energy between all the concerned elements, pacemaker circuit, bio-electrical contact, and cardiac tissue. The model proposed in this article includes many simplifications, for example it assumes that the contact properties calibrated from bench experiments do not change when the contact is made with cardiac tissue. However, we prefered not to use a more descriptive model that could not be validated against the experimental data at hand.

We are currently deriving and implementing a more reliable 3D model, with bidomain equations for the cardiac tissue in blood. This model will characterize more accurately the delivery of current into cardiac cells, accounting for the spatial distribution of current. In any case, in order to derive such a 3D model and computational solver, we need to calibrate the bio-physical parameters with respect to experimental conditions. The calibration and validation process will benefit from data that were recorded, alongside the EGMs described in the previous section. In addition, optical mapping data were also recorded. They may allow us to better monitor the propagation of the action potential in the tissue. High resolution 9.4T MR images of the anatomy were also acquire, so as to better characterize the tissue structure, prior to simulations.

We also plan to use better techniques than a 2D dichotomy to localize the simulated Lapicque curves, which will reduce the computational burden that will arise when using 3D bidomain equations.

Fig. 1 .

 1 Fig. 1. Equivalent circuit of a pacemaker (black, left) coupled to either a resistive medium (blue, middle), or a single cell model (red, right) via connecting dots. On the left circuit, S pulse is closed and S ocd is open during the pulse step, S pulse and S ocd are open during switch step, S pulse is open and S ocd is closed during the OCD step. The resistances R pulse and R gnd are located within the canister of the pacemaker, whereas the resistances Rw1 and Rw2 are equivalent resistances for the wires between the canister and the leads.

Fig. 2 .

 2 Fig. 2. Current (left) and voltage (right) delivered by a pacemaker lead in a saline solution (black solid line), between the tip and ring electrodes. Red dashed lines are the result of the calibration (Sect. 3). The voltage and duration set on the device were 2 V and 1 ms. (Color figure online)

Fig. 3 .

 3 Fig. 3. Left: sketch of the bench test used to calibrate the electrode contact model. Right: right ventricle endocardial view of the tissue wedge preparation with the pacemaker electrode implanted in the septum (pointed by arrow). Other implantations sites were located in the apex (×) and at the base of the ventricle (+).

Fig. 4 .

 4 Fig. 4. Fitted currents and voltages when calibrating parameters θ = (R1, R2, τ1, τ2). Gray solid lines are one experimental dataset, blue dotted lines are computed with θ calibrated with this dataset, and red solid lines from the calibration of the 9 datasets together. (Color figure online)

Fig. 5 .

 5 Fig. 5. Transmembrane voltage (left) computed for a single set of parameters (Vstim, d, S, ge); and Lapicque curve (right) obtained after detection of a 50% capture ratio. Blue points are for 0% capture, and red points for 100% capture. (Color figure online)

Fig. 6 .

 6 Fig. 6. Comparison of Lapicque curves (solid lines) computed with different values of S (left) and R1 (right). The shaded regions indicate where the experimental capture threshold is located, for different stimulation sites on the tissue wedge.

Acknowledgements. This work was supported by the H2020 EU SimCardioTest project (Digital transformation in Health and Care SC1-DTH-06-2020; grant agreement number 101016496). This study received financial support from the French Government as part of the "Investments of the Future" program managed by the National Research Agency (ANR), Grant reference ANR-10-IAHU-04. Experiments presented in this paper were partially carried out using the PlaFRIM experimental testbed, supported by Inria, CNRS (LABRI and IMB), Université de Bordeaux, Bordeaux INP and Conseil Régional d'Aquitaine.


Saeed Shoaie

Pouyan Ghaffari

Petia Kovatcheva-Datchary

Adil Mardinoglu

Partho Sen

Estelle Pujos-Guillot

Tomas De Wouters

Catherine Juste

Salwa Rizkalla

Julien Chilloux

Lesley Hoyles

Jeremy K Nicholson

ANR Microobese Consortium

Joel Dore

Marc E Dumas

Karine Clement

Fredrik Bäckhed

Jens Nielsen
email: nielsenj@chalmers.se

Quantifying diet-induced metabolic changes of the human gut microbiome

The human gut microbiome is known to be associated with various human disorders, but a major challenge is to go beyond association studies and elucidate causalities. Mathematical modeling of the human gut microbiome at a genome-scale is a useful tool to decipher microbe-microbe, diet-microbe and microbe-host interactions. Here, we describe the CASINO (Community and Systems-level Interactive Optimization) toolbox, a comprehensive computational platform for analysis of microbial communities through metabolic modeling. We first validated the toolbox by simulating and testing the performance of single bacteria and whole communities in in vitro. Focusing on metabolic interactions between the diet, gut microbiota and host metabolism, we demonstrated the predictive power of the toolbox in a diet-intervention study of 45 obese and overweight individuals, and validated our predictions by fecal and blood metabolomics data. Thus, modeling could quantitatively describe altered fecal and serum amino acid levels in response to diet intervention.

Introduction

Increasing evidence indicates that changes in the composition of the human gut microbiota affect host metabolism and are associated with a variety of diseases [START_REF] Backhed | Host-bacterial mutualism in the human intestine[END_REF][START_REF] Qin | Alterations of the human gut microbiome in liver cirrhosis[END_REF]. Changes in diet have been shown to rapidly affect the composition of the gut microbiota [START_REF] David | Diet rapidly and reproducibly alters the human gut microbiome[END_REF][START_REF] Wu | Linking long-term dietary patterns with gut microbial enterotypes[END_REF]. Furthermore, microbiota-diet interactions impact host physiology through the generation of a number of bioactive metabolites [START_REF] Cotillard | Dietary intervention impact on gut microbial gene richness[END_REF][START_REF] Chatelier | Richness of human gut microbiome correlates with metabolic markers[END_REF][START_REF] Nicholson | Host-gut microbiota metabolic interactions[END_REF][START_REF] Wu | Comparative metabolomics in vegans and omnivores reveal constraints on diet-dependent gut microbiota metabolite production[END_REF]. For example, short-chain fatty acids (SCFAs), which are generated by microbial fermentation of dietary polysaccharides in the gut, are an important energy source for colonocytes and also function as signaling molecules, modulating intestinal inflammation and metabolism [START_REF] Donohoe | The microbiome and butyrate regulate energy metabolism and autophagy in the mammalian colon[END_REF][START_REF] Fernandes | Adiposity, gut microbiota and faecal short chain fatty acids are linked in adult humans[END_REF][START_REF] Samuel | Effects of the gut microbiota on host adiposity are modulated by the short-chain fatty-acid binding G protein-coupled receptor, Gpr41[END_REF]Smith et al., 2013b;[START_REF] Tolhurst | Short-chain fatty acids stimulate glucagon-like peptide-1 secretion via the G-protein-coupled receptor FFAR2[END_REF]. In addition, dietary proteins and amino acids are important substrates for microbial fermentation in the colon [START_REF] Cummings | Colonic microflora: nutrition and health[END_REF] where they also serve as an important nitrogen source for the microbiota and support the growth of the microbiota and the host [START_REF] Wu | Amino acids: metabolism, functions, and nutrition[END_REF].

By quantifying the release and consumption of metabolites by the gut microbiota, it may be possible to elucidate interactions between the gut microbiota and host metabolism [START_REF] Tremaroli | Functional interactions between the gut microbiota and host metabolism[END_REF]. This information would allow identification of diagnostic biomarkers and may provide insight into the role of the gut microbiota in disease progression [START_REF] Karlsson | Gut metagenome in European women with normal, impaired and diabetic glucose control[END_REF][START_REF] Qin | Alterations of the human gut microbiome in liver cirrhosis[END_REF][START_REF] Zeller | Potential of fecal microbiota for early-stage detection of colorectal cancer[END_REF]. A predictive systems-level model of the human gut microbiome is required to elucidate causalities and quantify the interactions between microbes, host and diet [START_REF] Greenblum | Towards a predictive systems-level model of the human microbiome: progress, challenges, and opportunities[END_REF][START_REF] Manor | Mapping the Inner Workings of the Microbiome: Genomic-and Metagenomic-Based Study of Metabolism and Metabolic Interactions in the Human Microbiome[END_REF][START_REF] Shoaie | Elucidating the interactions between the human gut microbiota and its host through metabolic modeling[END_REF].

A genome-scale metabolic model (GEM) is an integrative platform for exploring genotypephenotype relationships and metabolic differences between different clinical conditions [START_REF] Ghaffari | Identifying anti-growth factors for human cancer cell lines through genome-scale metabolic modeling[END_REF][START_REF] Henry | High-throughput generation, optimization and analysis of genome-scale metabolic models[END_REF][START_REF] Mardinoglu | Genomescale metabolic modelling of hepatocytes reveals serine deficiency in patients with non-alcoholic fatty liver disease[END_REF][START_REF] Mardinoglu | New paradigms for metabolic modeling of human cells[END_REF][START_REF] Monk | Optimizing genome-scale network reconstructions[END_REF][START_REF] Shoaie | Elucidating the interactions between the human gut microbiota and its host through metabolic modeling[END_REF]. We previously reconstructed GEMs to study the interactions between Bacteroides thetaiotaomicron and Eubacterium rectale [START_REF] Shoaie | Understanding the interactions between bacteria in the human gut through metabolic modeling[END_REF], representatives of Bacteroidetes and Firmicutes, the two dominant phyla in the human gut [START_REF] Huttenhower | Structure, function and diversity of the healthy human microbiome[END_REF], and between Bifidobacterium adolescentis and Faecalibacterium prausnitzii [START_REF] El-Semman | Genome-scale metabolic reconstructions of Bifidobacterium adolescentis L2-32 and Faecalibacterium prausnitzii A2-165 and their interaction[END_REF], also dominant and dietary-responsive gut microorganisms [START_REF] Walker | Dominant and diet-responsive groups of bacteria within the human colonic microbiota[END_REF]. In both studies, we manually identified the interactions between the bacteria and quantified the consumption and production rates of the defined interacting metabolites for each bacterial species. Although other studies have been conducted for communities of two and three species [START_REF] Harcombe | Metabolic resource allocation in individual microbes determines ecosystem interactions and spatial dynamics[END_REF][START_REF] Zomorrodi | OptCom: a multi-level optimization framework for the metabolic modeling and analysis of microbial communities[END_REF], these approaches cannot be expanded to simulate the interactions of a large number of species representing the complex gut ecosystem. We therefore developed the CASINO (Community And Systems-level INteractive Optimization) toolbox, which comprises an optimization algorithm integrated with diet analysis to predict the phenotypes and related dietary intake within the human gut microbiota. The toolbox was tested using both data from in vitro experiments and results from a nutritional intervention study of subjects with varying gut microbial gene richness.

Results

CASINO toolbox

We first developed an optimization algorithm in the CASINO toolbox, which is based on a collaborative and a multi-dimensional distributed approach [START_REF] Grimm | Pattern-oriented modeling of agent-based complex systems: lessons from ecology[END_REF]. It takes into account both collaboration between the multiple species as well as the fact that each individual species seeks to optimize its growth individually. Although GEMs are linear models, the presence of several GEMs in the overall community model means that the optimization of community biomass production becomes a non-linear problem. We therefore separated the community model into systems-level (representing the community) and organism-level (representing each species), which allowed us to linearize the optimization problem.

Simulations using CASINO start with an initialization stage that defines a primary profile of the systems-level topology (i.e. which species are present and how do they interact). This leads to the construction of a community matrix that defines effectors and receptors, with effectors being species that produce metabolites and receptors being species that consume metabolites. Following definition of the topology, the initialization step calculates metabolite production by each species using organism-level optimization. Thereafter, CASINO performs iterative multilevel optimization to calculate the relative uptake of carbohydrates by each species, until the total community biomass production is optimized. In this study, this calculation is constrained by the relative abundance of each species (Experimental Procedures & Figure S1).

To evaluate CASINO, we used the RAVEN toolbox [START_REF] Agren | The RAVEN toolbox and its use for generating a genome-scale metabolic model for Penicillium chrysogenum[END_REF] to update and significantly expand the content of our previously published GEMs for B. thetaiotaomicron, E. rectale, B. adolescentis and F. prausnitzii and to generate a GEM for Ruminococcus bromii, a representative of Clostridiales and a key gut symbiont [START_REF] Ze | Ruminococcus bromii is a keystone species for the degradation of resistant starch in the human colon[END_REF]. All GEMs were manually curated for functionality based on literature information. We defined a set of metabolic tasks, e.g. generation of biomass precursors (Table S1), to further investigate the functionality of the GEMs and checked that the resulting models could perform the defined tasks (Experimental Procedures).

The GEMs were functionally validated using experimental data for each of the five bacteria. We quantified the abundance of the bacteria by 16S ribosomal RNA (rRNA) quantitative PCR (qPCR) at baseline and after 24 h of growth in selected media (Table S2; Extended Experimental Procedures). We performed targeted metabolomics to quantify products of the fermentative activity of the studied bacteria, specifically the SCFAs butyrate, acetate and propionate and 15 different amino acids, and consumption of carbohydrates (glucose, maltose, cellobiose and starch). Flux constraints were imposed using the metabolomics profiles of the growth media, and maximum growth of each bacterial species was set as an objective function to simulate the predictive power of the corresponding model (Experimental Procedures and Figure 1A). The experimental data confirmed that the GEMs predicted the metabolism and biomass growth for each bacterial species almost perfectly (Figures 1B-1D). The GEMs correctly predicted that acetate can be produced by B. adolescentis, B. thetaiotaomicron and R. bromii, butyrate can be produced by E. rectale and F. prausnitzii and propionate can only be produced by B. thetaiotaomicron (Figure 1B). Our simulations also predicted that these five bacteria synthesize significantly higher level of essential amino acids (valine, leucine, methionine, lysine and phenylalanine) compared to non-essential amino acids (serine, tyrosine and threonine) (Figure 1D).

Revealing the interactions between constituents of in-vitro microbial communities

To test the performance of CASINO, we simulated the interactions between the microbes in two microbial communities that differed only in one bacterial species: EBBR (E. rectale, B. adolescentis, B. thetaiotaomicron and R. bromii) and FBBR (F. prausnitzii, B. adolescentis, B. thetaiotaomicron and R. bromii) (Figure 2A). The simulated values were validated by culturing EBBR and FBBR communities in selected media. We quantified the abundance of individual bacterial species in each community by 16S rRNA-based qPCR (Table S3). We also performed targeted metabolomics to quantify the production of SCFAs and amino acids and the consumption of carbohydrates (starch and cellobiose) for each community (Extended Experimental Procedures).

Our model simulations correctly predicted the net production of the metabolites produced by each community and showed that the communities synthesized more essential amino acids than nonessential amino acids (Figure 2B). More importantly, the simulations enabled quantification of the contribution of each individual bacterial species to the overall microbial conversion in the communities and showed that two of the species in each community dominated. Specifically, we predicted that B. thetaiotaomicron and E. rectale synthesized 41% and 36% of the amino acids in the EBBR community, respectively, and B. thetaiotaomicron and F. prausnitzii synthesized 39% and 47% of the amino acids in the FBBR community, respectively (Figure S2). We also predicted that E. rectale mainly contributed to the synthesis of valine, leucine, phenylalanine and methionine in the EBBR community, while F. prausnitzii was the major contributor to the production of valine and leucine in the FBBR community. Furthermore, the experimental data showed that substitution of E. rectale with F. prausnitzii decreased the level of butyrate in the media, due to the higher capacity of E. rectale for butyrate production [START_REF] Louis | Diversity of human colonic butyrateproducing bacteria revealed by analysis of the butyryl-CoA:acetate CoA-transferase gene[END_REF], and also the model simulations showed a slightly lower butyrate production by the FBBR community compared with the EBBR community.

Next we calculated the centrality scores for each bacterial species in order to identify which species have a dominant role in the overall metabolic conversion in each community (Extended Experimental Procedures). We observed that E. rectale and B. thetaiotaomicron were the main receptor and effector, respectively and thus represent key species (Figure 3A). We then evaluated the sensitivity of the optimization algorithm in CASINO by adding bacteria in three steps to each of these two bacteria, culminating in reconstruction of the two in vitro communities. We calculated the SCFA levels for each step (Figure 3B). Addition of B. adolescentis to E. rectale in the EBBR community resulted in reduced production of butyrate and increased production of propionate and acetate. Addition of F. prausnitzii to B. thetaiotaomicron in the FBBR community resulted in reduced production of propionate and acetate and increased butyrate production. The levels of the SCFAs changed further when the other species were added (Figure 3B).

Analyzing the effect of gene richness and diet on gut microbiota composition

To further evaluate CASINO, we examined data from a clinical study where 45 overweight and obese individuals were subjected to an energy-restricted high-protein diet with low glycemic index for 6 weeks (Figure 4A; clinical data in Table S4). These patients had previously been stratified based on their gut microbial gene richness into 'low gene count' (LGC; n = 18) and 'high gene count' (HGC; n = 27) based on a cut-off threshold of 480,000 genes [START_REF] Cotillard | Dietary intervention impact on gut microbial gene richness[END_REF]. LGC demonstrated a worse metabolic profile compared with HGC individuals [START_REF] Cotillard | Dietary intervention impact on gut microbial gene richness[END_REF].

Analysis of metagenomics data before and after the diet intervention showed that six species dominated in all subjects: Escherichia coli and F. prausnitzii and four species associated with Clostridia, Bacteroides, Bifidobacteria, and Lactobacillus [START_REF] Cotillard | Dietary intervention impact on gut microbial gene richness[END_REF]. To obtain quantitative data of these species, we analyzed fecal samples by 16S rRNA qPCR before and after the dietary intervention (Table S5 and Extended Experimental Procedures) and used these results to calculate the distribution of biomass between the species (Figure 4B). We observed significant differences in abundance for B. adolescentis, F. prausnitzi and E. rectale at baseline and for B. adolescentis and L. reuteri after 6 weeks between LGC and HGC individuals. In HGC individuals, levels of B. thetaiotaomicron significantly increased and L. reuteri and F. prausnitzii significantly decreased after 6 weeks of dietary intervention compared with baseline, whereas a significant decrease in LGC individuals was only seen for L. reuteri.

Diet interventions alter amino acid and SCFA levels in HGC and LGC individuals

In order to simulate the effect of the diet on the overall gut microbiota metabolism, we used representatives of the most abundant microbial groups that we had also modeled in vitro, i.e. B. thetaiotaomicron, B. adolescentis, F. prausnitzi, E. rectale as described above and Lactobacillus reuteri, for which we reconstructed a GEM. We also performed simulations with inclusion of E. coli, but as this species had no major impact on the production of SCFAs and amino acids (data not shown), we did not include this species in our further analysis. Using CASINO, we simulated the effect of diet on the human gut microbiota composition at baseline and after the dietary intervention for 44 of the subjects (registered diet information in Table S6). In order to translate the diets into metabolites that can be utilized by the five gut bacterial species, we computed the dietary macronutrients of 24 different food items (Table S7) and used this information in a diet allocation algorithm in CASINO. Using this algorithm, CASINO predicted that there was a decrease in carbohydrate consumption and an increase in amino acid consumption for all individuals after 6 weeks of dietary intervention (Figure 4C). The intake of fiber from bread and potatoes, and to lesser extent rice, cereals and snacks, was decreased but fiber from fruits and vegetables was increased in agreement with the dietary recommendation given to the patient during the intervention (Figure 4C).

For the simulations, we assumed that carbohydrates and fibers were hydrolyzed to glucose to the same degree in all subjects, allowing us to calculate the relative amount of glucose available to the gut microbiome in each subject. We further assumed that glucose was the limiting substrate for the gut microbiome. We first used CASINO to quantify the community interactions and the relative glucose uptake by the individual species. We used the calculated values of species abundance in this process. Thereafter, we used CASINO to repeat the simulations, but now allowing the individual species to consume amino acids in the same ratio as their glucose uptake. The amount of available amino acids was calculated from the diet composition using CASINO.

With this approach we could simulate the profile of three SCFAs and 14 amino acids produced by the gut ecosystem as well as the contribution of each microbial species to the overall metabolite production of the ecosystem at baseline and after 6 weeks of dietary intervention for each individual. By plotting average profiles for all the subjects, we found that the levels of the SCFAs and amino acids produced by the gut microbiota were significantly decreased after dietary intervention when both LGC and HGC groups were combined (Figure 5A, decrease in the y-axis direction), but the greatest reductions were observed in LGC individuals (Figure 5A, increase in the x-axis direction).

In order to experimentally evaluate our predictions on altered metabolite production by the gut ecosystem, we performed metabolomics analysis of fecal samples obtained from the HGC and LGC individuals at baseline and after 6 weeks of dietary intervention. These data confirmed the many of the predicted simulations by CASINO, i.e. alanine, proline, glycine, serine, phenylalanine, and tyrosine all showed a decreased level in response to the diet intervention in both HGC and LGC subjects but, with a larger decrease in the LGC subjects (Figure 5B, shift down-wards right). To test the significance of these changes for each group of subjects and between the two time points we calculated p-values using a student t-test, and except for alanine these changes were significant for the different groups (Figure 5C). Measured serine levels were significantly higher in LGC than HGC individuals at baseline but not after 6 weeks of dietary intervention (Figure 5B), in agreement with the predicted results (Figure 5A). Also, measured phenylalanine levels were significantly higher in LGC individuals than in HGC individuals at baseline but lower in LGC individuals compared to HGC subjects after 6 weeks of dietary intervention, in agreement with predicted results.

In addition to predicting changes in some of the metabolites in response to dietary intervention, the model could also be used to predict the relative contribution of each bacterial species to production of specific metabolites, allowing us to quantitatively access how a variation in the gut microbiome correlates with metabolite production. Thus, we predicted the contribution of each bacterial species to phenylalanine levels in the gut ecosystem and showed that 23% of the total phenylalanine is produced by B. adolescentis and 26% by E. rectale in HGC individuals at baseline while this contribution increased for B. adolescentis to 26% and decreased for E. rectale by 21% after 6 weeks of dietary intervention (Figure S3). For LGC individuals, the contribution of E. rectale to phenylalanine production was 29% at baseline and decreased to 15% after 6 weeks of dietary intervention (Figure S3).

Serum metabolomics confirm model predictions and associate with clinical parameters

Although the model simulations could correctly predict changes in several of the metabolites in the feces, we noted that the model simulations did not accurately predict changes in all the measured metabolites, which may be a result of differential absorption by the host. We therefore evaluated whether the model could predict changes in the serum. We used metabolomics to analyze serum of the 45 subjects at baseline and after 6 weeks of dietary intervention and found an excellent correspondence between the model predictions (Figure 5A) and the measured changes (Figure 6A). The serum levels of 10 detected amino acids decreased in response to the dietary intervention in the LGC subjects (Figure 6A). Furthermore, in agreement with the model predictions, there was a decrease in acetate in response to the diet intervention in all subjects (Figure 6A). In addition, we observed that phenylalanine levels were higher in LGC subjects compared to HGC subjects at baseline, but the level of phenylalanine decreased in LGC subjects after 6 weeks of dietary intervention (Figure 6Aand 6B). We also observed that levels of valine, leucine and alanine were higher in LGC subjects at baseline (Figure 6A and6B).

To evaluate if these changes in serum metabolite levels may have any clinical relevance, we analyzed the correlations between the levels of the 10 amino acids in the serum and bioclinical parameters of the subjects at baseline (Figure 6C). Here we found that the serum phenylalanine levels were positively correlated with clinical variables related to body corpulence (BMI, DXAmeasure fat mass, waist circumference, leptin), insulin resistance, blood lipid homeostasis (serum triglycerides and cholesterol), and low-grade inflammation (hsCRP). The serum levels of valine, leucine were also positively correlated with BMI.

Model-based diet design to improve metabolism of LGC individuals

Finally, assuming that LGC subjects have a non-optimal gut microbiome metabolism (associated with a clinically altered metabolism) we wanted to identify which dietary change would improve the metabolism of their gut microbiome. We thus made the assumption that an adapted dietary recommendation in LCG subjects provided at baseline would enable them to reach the "optimal" gut microbiome metabolism of HGC subjects after 6 weeks of dietary intervention, which is indeed associated with an improved metabolic phenotype.

.We then used CASINO and the abundance of the five different species, i.e. B. thetaiotaomicron, B. adolescentis, F. prausnitzi, E. rectale and L.reuteri, to predict the relative consumption of 8 essential amino acids by the gut microbiome in the LGC subjects at baseline (base phenotype in Figure 7A) and in the HGC subjects at week 6 (improved phenotype in Figure 7A).

From this model analysis, we found that the gut microbiome of HGC individuals had a higher consumption of these 8 essential amino acids at week 6 compared to that of the LGC subjects at baseline. An incremental augmentation of these amino acids would permit to acquire a similar metabolism of the gut microbiome in LGC and HGC subjects (Figure 7A). Many different combinations of food sources could fulfill such a requirement for essential amino acids. However, in an attempt to identify some overall guidelines, we correlated the difference between the two different requirements of amino acids with the composition of these amino acids in different food types (Table S8). This showed that LGC individuals should significantly increase consumption of dairy products, vegetables, white meat, fish pulses, eggs, oils and butter. In the meantime, they should considerably reduce intake of pastries, bread and rice to improve and slightly reduce intake of cereals and nuts (Figure 7B).

Discussion

The overall metabolism of the gut microbiome can be modeled in one of two ways: 1) by using a lumped model of all the metabolic reactions active in the different gut microorganisms or 2) by compartmentalizing the metabolism according to the individual microorganisms. The latter is clearly a better reflection of the true biological system and it also ensures that redox and energy balances are constrained within each organism considered. We therefore used this approach to model the metabolism of the human microbiome and reconstructed GEMs for individual species from the predominant phyla in the human gut. We identified which species to include in our analysis based on their abundance in the gut ecosystem. Thus, we reconstructed GEMs for five species that are representative bacteria of the dominant phyla in the human gut and we hypothesized that the reactions included in our models cover most of the metabolic functions that are present in the human gut. Compared with earlier attempts to model the human gut metabolism using GEMs, i.e. the COMETS algorithm [START_REF] Harcombe | Metabolic resource allocation in individual microbes determines ecosystem interactions and spatial dynamics[END_REF], CASINO allows inclusion of several species in the simulations. Furthermore, it is scalable and enables expansion to include even more than the five species that we considered in this study.

In order to evaluate whether we are covering the main metabolic functions, we simulated the effect of different diets, studied the interactions between the microbes and host in response to the diet, and quantified the contribution of each bacterial species to the fecal metabolite profiling. The model simulations matched fecal metabolomics data, but more importantly it correctly correlated with changes in serum levels of 10 amino acids and one SCFA (acetate). Thus, the model captures some major metabolic functions of the human gut microbiome. In the future, the selection of species to be considered should be expanded, in particular to ensure representation of more specific metabolic functions such as vitamin biosynthesis and bile acid metabolism.

The consistency between model predictions of metabolite productions and measurement of changes in metabolite levels in feces and serum suggests that our modeling correctly predicts the overall carbon fluxes in the gut ecosystem. Furthermore, our simulations enabled quantification of how the individual species compete for nutrients and produce different metabolites that may serve as nutrients for other species or be absorbed by the host. Studying the gut metabolism with our holistic approach also enabled understanding of metabolic shifts under different clinical conditions, and hereby could provide a direct link between the gut microbiome metabolism and serum chemistry. Thus, our simulations suggest that the gut microbiome may contribute to altered levels of several amino acids in the serum, including phenylalanine and branched chain amino acids. This is in line with an early report, using germ free mice, showing that the microbiota of the large intestine increased the free amino acid level in the gastrointestinal tract [START_REF] Macfarlane | Contribution of the microflora to proteolysis in the human large intestine[END_REF]. A later study showed that bacteria in the human large intestine take up peptides and amino acids, and convert these to different amino acids and SCFAs [START_REF] Smith | Enumeration of amino acid fermenting bacteria in the human large intestine: effects of pH and starch on peptide metabolism and dissimilation of amino acids[END_REF]. This study also showed that the production of amino acids was dependent on the composition of starch, proteins and peptides, and hence will be dependent on the dietary composition. Further confirmation of our findings is a documented in a recent review on the role of microbial amino acid metabolism in host metabolism, that provide a summary on a number of findings related to the role of the microbiota in the large intestine on production of not only SCFAs but also amino acids that are subsequently taken up by the host [START_REF] Neis | The Role of Microbial Amino Acid Metabolism in Host Metabolism[END_REF]. However, we are aware that more experiments are required to validate if the gut microbiome actually contribute to host amino acid metabolism.

In our study, obese individuals with a LGC microbiome, associated with more impaired metabolic phenotype compared with those with HGC, had elevated levels of these amino acids. A previous study has shown that phenylalanine is associated with type 2 diabetes (T2D), and serum levels of this essential amino acid are five to seven fold higher in individuals at risk of T2D compared to control subjects [START_REF] Wang | Metabolite profiles and the risk of developing diabetes[END_REF]. Furthermore, the serum level of branched chain amino acids (valine, leucine and isoleucine) has been found to correlate with insulin resistance [START_REF] Newgard | A branched-chain amino acid-related metabolic signature that differentiates obese and lean humans and contributes to insulin resistance[END_REF]. Circulating levels of leucine, arginine, valine, proline, phenylalanine, isoleucine and lysine are also significantly associated with an increased risk of hypertriglyceridemia in diabetic subjects (Mook-Kanamori et al., 2014), a phenotype we noted in the subjects with the LGC microbiome. Thus, our simulations point to two important findings. First it suggests that the gut microbiota in LGC individuals may contribute to increased serum levels of many amino acids that have been found correlated with metabolic diseases such as T2D, and we believe this is consistent with the deteriorated glucose homeostasis related to insulin resistance observed for the LGC subjects at baseline both in French and Danish subjects [START_REF] Cotillard | Dietary intervention impact on gut microbial gene richness[END_REF][START_REF] Chatelier | Richness of human gut microbiome correlates with metabolic markers[END_REF]. The HGC subjects on the other hand, with their more gene rich gut microbiome may have a better conversion of amino acids, resulting in lower levels of these in the plasma. Furthermore, the HGC microbiome also have higher capacity to produce SCFAs that are important energy sources for the colonocytes as well as they function as signaling molecules, modulating intestinal inflammation and metabolism [START_REF] Donohoe | The microbiome and butyrate regulate energy metabolism and autophagy in the mammalian colon[END_REF][START_REF] Fernandes | Adiposity, gut microbiota and faecal short chain fatty acids are linked in adult humans[END_REF][START_REF] Samuel | Effects of the gut microbiota on host adiposity are modulated by the short-chain fatty-acid binding G protein-coupled receptor, Gpr41[END_REF]Smith et al., 2013b;[START_REF] Tolhurst | Short-chain fatty acids stimulate glucagon-like peptide-1 secretion via the G-protein-coupled receptor FFAR2[END_REF]. Second, in agreement with the original study on the LGC in comparison to the HGC individuals, it highlights that LGC subjects can benefit from a dietary intervention that improve their gut microbiome metabolism, paving the way for a personalized approach in these subjects. Indeed, despite higher levels of some amino acids at baseline, LGC subjects had a larger decrease in the levels of a range of metabolites that are positively correlated with insulin resistance markers and cardio-metabolic risk factors.

Using our approach we also predicted the relative contribution of each bacterial species to production of specific metabolites, and studied how this variation in the gut microbiome is correlated with specific metabolite production. Information generated from CASINO may therefore be extended for rational design of prebiotics as well as for identifying novel beneficial bacteria that can be used to fortify the microbiota to improve the gut microbiome metabolism. Importantly, rational design of microbiota interventions requires knowledge of diet as demonstrated in a study of children with kwashiorkor which showed that a disrupted microbiome that can be reversed by dietary interventions (Smith et al., 2013a). Interestingly, transferring the microbiota from children with kwashiorkor to germ-free mice in combination with Malawian diet resulted in marked weight loss in recipient mice associated with perturbations in amino acids. As we demonstrated CASINO can also be used to predict dietary changes required in order to ensure a certain profile of the gut metabolism, here represented as a specific consumption of 8 essential amino acids. The gut microbiome may change in response to dietary modulation, something that our simulations are not capturing. This study also emphasizes the importance of developing accurate tools to properly record dietary intakes in different populations. With more data it will probably be possible to also predict how the diet influences gut microbiome changes, and hereby CASINO may assist in the development of a precision medicine approach to treat metabolic diseases associated with dysfunction of the gut microbiota.

In conclusion, we demonstrate how we can use model simulations to predict metabolic interactions within the gut microbiome and hereby assist in generating mechanistic insight into the contribution of individual species of the gut microbiome to the overall metabolism of the ecosystem and the host. Furthermore, focused on the diet and on host and gut microbiota metabolic interactions we show how the gut ecosystem and the individual members of the gut microbiota contribute to the host metabolism. CASINO may thus constitutes a valuable tool for enriching the information content provided by gut metagenome analysis, and hereby advance our understanding on how this important metabolic organ contributes to disease development and thus facilitate personalized interventions based on the microbiome.

Experimental procedures

Reconstruction/updating of GEMs The B. adolescentis, B. thetaiotaomicron, E. rectale and F. Prausnitzii metabolic models were already published [START_REF] El-Semman | Genome-scale metabolic reconstructions of Bifidobacterium adolescentis L2-32 and Faecalibacterium prausnitzii A2-165 and their interaction[END_REF][START_REF] Shoaie | Understanding the interactions between bacteria in the human gut through metabolic modeling[END_REF].These models were validated individually and updated based on extensive bibliometric survey of the literature and databases [START_REF] Henry | High-throughput generation, optimization and analysis of genome-scale metabolic models[END_REF]. The R. bromii and L. reuteri metabolic networks were reconstructed based on automatic and manual curation, considering information available in the literature and databases. For model reconstructions, updating and quality checks, the RAVEN toolbox was used [START_REF] Agren | The RAVEN toolbox and its use for generating a genome-scale metabolic model for Penicillium chrysogenum[END_REF]. The defined metabolic task file was implemented to perform the gap filling process. This task file includes synthesis of amino acids, nucleotides and carbohydrates. A published metabolic model for E. coli was used for simulation [START_REF] Monk | Genome-scale metabolic reconstructions of multiple Escherichia coli strains highlight strain-specific adaptations to nutritional environments[END_REF].

Validation of GEMs

The metabolic models for B. adolescentis, B. thetaiotaomicron, E. rectale and F. prausnitzii and R. bromii were used to predict the experimental phenotypes based on the individual in-vitro data. The growth for each model was predicted as an objective function based on the available media components from the experiments. Measured amino acids and SCFAs secretion were used to constrain the models for individual simulations. The lower and upper bounds for the uptake and secretion reactions were assigned based on standard deviation of the measured metabolites in experimental data. The results were tested to see the consistency of the model through fixation of non-growth associated maintenance ATP.

CASINO toolbox

CASINO defines the primary topology of the community and identifies optimum solution by applying a multi-dimensional optimization procedure in two successive and connected stages: Initialization and Community optimization. The convexity of the solution space and linearity of the parameters are maintained by separating the problem into two independent optimizing dimensions, systems-level and organism-level, with three classes of variables: inputs, outputs and connecting parameters. Maximum biomass production was considered as objective function at both levels. At the organism level, each species seeks to maximize its biomass yield while at the system level the community seeks to stay in optimum balanced condition by synchronizing the competition between species (Figure S1).

Initialization process assumes a structure of the community as a complex network and starts with an activation cascade supporting a specific threshold. Species are categorized into primary (grows independently by taking up system-level inputs) and non-primary (growth is dependent on connecting inputs, i.e. metabolites produced by other species) classes. The activation begins with identifying primary species within the community and activates them by providing required resources to grow. The compounds produced by activated species are added to the resource pool and the community is screened again to find non-primary species that now can grow based on the updated resource pool. This cascade of activation is repeated until the whole network has been activated. Now a feasible profile of community topology is constructed and this results in a definition of the community constraints matrix and the community objective function. At the end of the initialization process the system is locally optimum (species grow on maximum biomass yield rate) but globally non-optimum (resources distribution between species do not satisfy community optimum conditions).

Community optimization, a multi-level iterative process, starts based on the initial feasible space and a community objective function defined by the initialization procedure. In each step, the community level biomass production is optimized to find the optimum distribution of resources between species. Relative carbohydrate uptake rates are calculated in each intermediate systems-level optimum condition and these are used by the species to reach organism-level optimality. Metabolites secretion rates obtained after organism-level optimization is used to expand the boundaries of the intermediate feasible space and find a new system-level optimum. This iterative process continues until the solution converges to local and global optima.

Objective function in this algorithm is considered as maximum biomass yield (local force) at the organism-level and maximum community biomass (community force) at the system-level. Summation of these two forces defines the direction of optimization through expansion of the feasible space. The community force is adjusted by centrality scores assigned to each species based on the network topology of the community established in the initialization procedure. Two centrality degrees, power centrality and degree centrality, are used to calculate controlling power of species as effectors and receptors. 𝛿 is a binary vector that activates certain parameter in the objective equation. This binary vector is corrected by imposing the relative centrality scores of species. μ(i) is the biomass yield of individuals and μ PL is biomass yield of species belong to PL.

Statistical Analysis

Peaks obtained from 144 samples (HGC and LGC subjects) at two different time points (0 hours and 6 weeks) was aligned and subjected to retention time correction, XCMS [START_REF] Smith | XCMS: Processing mass spectrometry data for metabolite profiling using Nonlinear peak alignment, matching, and identification[END_REF]. About 295 features (mz) with aligned peak were normalized by Quantile normalization method [START_REF] Amaratunga | Analysis of data from viral DNA microchips[END_REF]. Mass fragmentation spectra of each featured peak were matched against the reference spectra obtained from Metlin [START_REF] Tautenhahn | An accelerated workflow for untargeted metabolomics using the METLIN database[END_REF] and HMDB [START_REF] Wishart | HMDB 3.0-The Human Metabolome Database in 2013[END_REF] using exact template matching approach [START_REF] Pavlidis | Using ANOVA for gene selection from microarray studies of the nervous system[END_REF]. Spectra with 90% match probability were considered for further analysis. Again, the retention index of the matched peak(s) (mz) was verified with metabolites in HMDB [START_REF] Wishart | HMDB 3.0-The Human Metabolome Database in 2013[END_REF].Thus, spectral match and retention index ensure identification and annotation of metabolites. Intensity (or expression) of 15 metabolites of interest was extracted for HGC and LGC subjects at two time points. Two sample t-tests were performed for detecting the significant metabolites and P value < 0.05 was considered to be significant. All data are shown as mean ± 1 SD. The results were compared with data from triplicate in-vitro experiments for EBBR and FBBR communities. In CASINO the interactions of the bacteria as well as the phenotype of the community were identified using an optimization algorithm. Growth of each bacterium had local optimum whereas the community had global optimum. The community optimum was detected by the intersection point of the fixed constraints for the community and the calculated dynamic constraints, which was obtained by summation of the local and community forces. See Figure S1 for development of the toolbox. See Table S3 We found that synthesis of essential AAs (histidine, isoleucine, leucine, lysine, methionine, phenylalanine, threonine and valine) produced by the communities is higher than the production of non-essential AAs (alanine, glutamate, glycine, proline, serine and tyrosine). Data are presented as mean value ± SD. See Figure S2 for details of AAs prediction. profile upon adding different species to the community. First the most important receptor and effector in the communities were identified using the result of Figure 3A. 1 mmol/l of glucose was used for all the simulations and the SCFAs profile were predicted. Following identification of the dominant receptor and effector, the key species, the other species were added to the community one by one until the EBBR and FBBR community were reconstructed. Comparison between the simulations showed that the SCFA profile is very sensitivity of absence and existence of species with respect to their abundance and interactions. S5 for more details. (c) A diet algorithm was developed and implemented for prediction of the macromolecules present in different food sources and this allowed further conversion of diets to three main categories of macronutrients carbohydrates, fiber and AAs. See table S6 and S7 for more details.

Figure 5 The effect of the 6 weeks diet interventions in predictions and fecal metabolomics of HGC and LGC individuals (a) Summary of average phenotypic predictions for baseline and after 6 weeks. The group of metabolites in the top right hand of the figure denotes the predictions at baseline and at the bottom left represents predictions after 6 week of dietary intervention. Subtracting the log 10 average metabolite fluxes for HGC from LGC is represented on the x-axis and the summation on the y-axis. The x-axis shows the ratio of predicted metabolite levels between HGC and LGC and the y-axis the sum of predicted metabolite levels in the two groups. The colors show the metabolites distance from zero in y-axis (from dark blue at the top to dark red at the bottom). (b) Metabolomics analysis of fecal samples obtained from HGC and LGC subjects. The differences are shown for 14 detected AAs as well as for butyrate. (c) p-values based on Students t-test for specification of significantly changed metabolites for four different comparisons. B stands for baseline (before the dietary intervention), w6 stands for week 6 (after the dietary intervention). S8 for more details.

  Maximize 𝑆 = [𝜑 × 𝜃 ′ ]. 𝐶 𝐸 . 𝜔 𝜇 . 𝛿 Where 𝑠 = 𝛼. 𝑋 + 𝛽. 𝑍 𝑖𝑛 + 𝛾. 𝑌 + 𝛿. 𝑍 𝑜𝑢𝑡 ω μ = 𝜇(i) min(μ PL ) 𝜑 = [𝛼, 𝛽, 𝛾, 𝛿] ∶ 𝐶𝑜𝑒𝑓𝑓𝑖𝑐𝑖𝑒𝑛𝑡𝑠 𝑚𝑎𝑡𝑟𝑖𝑥 𝜃 = [𝑋, 𝑌, 𝑍] ∶ 𝐼𝑛𝑝𝑢𝑡𝑠, 𝑜𝑢𝑡𝑝𝑢𝑡𝑠 𝑎𝑛𝑑 𝑐𝑜𝑛𝑛𝑒𝑐𝑡𝑖𝑛𝑔 𝑝𝑎𝑟𝑎𝑚𝑒𝑡𝑒𝑟𝑠 𝑃𝐿 ∶ 𝑃𝑟𝑖𝑚𝑎𝑟𝑦 𝑠𝑝𝑒𝑐𝑖𝑒𝑠 𝑙𝑖𝑠𝑡

Figure 1

 1 Figure 1 Validation of the GEMs (a) Each GEM was validated based on the rRNA and metabolomics data generated by in vitro experiments. The byproducts and the substrate usage were constrained in the models and the growth rate was compared with the experimental data. (b) Predicted and measured SCFA levels by the individual bacteria. Propionate was produced only by B. thetaiotaomicron, while acetate was produced by Actinobacteria and Bacteroidetes phyla. Butyrate production was mainly produced by the bacteria from the class clostridia. (c) Predicted and measured biomass at the end of the fermentations. Growth was set as an objective function for each model and the predicted growth was compared with the experimental data. (d) Predicted and measured levels of AAs by the individual bacteria. Each model also predicted AA levels and the details of 15 significant amino acids produced are shown for each bacterium. The predicted and experimental values showed that all AAs could be produced in the range of experimental data with specific optimum solution. Data are presented as mean value ± SD. See also TableS1and S2 for more details.

Figure 2

 2 Figure 2 Validation of the CASINO Toolbox (a) Two in-silico microbial communities EBBR (E. rectale + B. adolescentis + B. thetaiotaomicron + R. bromii) and FBBR (F. prausnitzii + B. adolescentis + B. thetaiotaomicron + R. bromii) were designed and simulated using the CASINO Toolbox.The results were compared with data from triplicate in-vitro experiments for EBBR and FBBR communities. In CASINO the interactions of the bacteria as well as the phenotype of the community were identified using an optimization algorithm. Growth of each bacterium had local optimum whereas the community had global optimum. The community optimum was detected by the intersection point of the fixed constraints for the community and the calculated dynamic constraints, which was obtained by summation of the local and community forces. See FigureS1for development of the toolbox. See TableS3for more details. (b) Predicted and measured levels of SCFA and AAs by the two in-silico microbial communities including EBBR (E. rectale, B. adolescentis, B. thetaiotaomicron and R. bromii) and FBBR (F. prausnitzii, B. adolescentis, B. thetaiotaomicron and R. bromii).We found that synthesis of essential AAs (histidine, isoleucine, leucine, lysine, methionine, phenylalanine, threonine and valine) produced by the communities is higher than the production of non-essential AAs (alanine, glutamate, glycine, proline, serine and tyrosine). Data are presented as mean value ± SD. See FigureS2for details of AAs prediction.

Figure 3

 3 Figure 3 The network structure influence and sensitivity analysis on CASINO.(a) The community of B. thetaiotaomicron, B. adolescentis, F. prausnitzii, E.rectale and R.bromii were tested based on being receptor (receiving metabolites from the other microbes) or effector (producing metabolites where consumed by receptors). Two methods of centrality were tested on these networks (power centrality and degree

Figure 4

 4 Figure 4 The effect of the 6 weeks diet interventions in HGC and LGC individuals. (a) For each food source the major macronutrients were quantified and this enabled using the CASINO Toolbox to study the effect of the diet on the gut microbiota composition of subjects classified based on their microbial gene richness. (b) Abundance of species before and after diet interventions in HGC and LGC subjects. Data are shown as box-plots with E. coli (red), B. thetaiotaomicron (yellow), B. adolescentis (jade), L. reuteri (brown), F. prausnitzi (blue) and E. rectale (green). The heat map shows the p-values for four different comparison of the species levels (each row associated with the corresponding species indicated in the left part of the figure). See TableS5for more details. (c) A diet algorithm was developed and implemented for prediction of the macromolecules present in different food sources and this allowed further conversion of diets to three main categories of macronutrients carbohydrates, fiber and AAs. See tableS6and S7 for more details.

Figure 6

 6 Figure 6 The serum metabolomics validated the predictions and fecal metabolomics (a) Metabolomics of the serum samples obtained from HGC and LGC individuals. The 10 AAs and acetate were quantified at the baseline (b) and week 6 (w6).

Figure 7

 7 Figure 7 Modeling the dietary composition to transfer LGC individuals towards improved phenotype (a) The yellow circles specify the simulated consumption of the 8 essential AAs by the gut microbiome of the LGC individuals at baseline (base phenotype) and the green circles specify the simulated consumption of the 8 essential AAs for the HGC individuals at week 6 (improved phenotype). (b) After calculating the required amount of 8 essential AAs at base and improved phenotype, both patterns were correlated with composition of AAs in different food categories. The direction of the Corr Improved -Corr Base indicates the positive/negative effect of different food sources to improve the phenotype of LGC subjects. (Corr Base : correlation between patternAAs in base phenotype and food, Corr Improved : correlation between pattern AAs in improved phenotype and food). See TableS8for more details.

  for more details. (b) Predicted and measured levels of SCFA and AAs by the two in-silico microbial communities including EBBR (E. rectale, B. adolescentis, B. thetaiotaomicron and R.

bromii) and FBBR (F. prausnitzii, B. adolescentis, B. thetaiotaomicron and R. bromii).

Acknowledgements

We thank Dr. Rosie Perkins (University of Gothenburg) for editing the manuscript, Dr Ana Luisa Neves for NMR sample preparation and Suwanee Jansa-Ard for metabolomics sample preparation. This project was supported by Knut and Alice Wallenberg Foundation, the Bill & Melinda Gates Foundation, Torsten Söderbergs Stiftelse, European Commission FP7 project METACARDIS with the grant agreement HEALTH-F4-2012-305312, Fondation Coeur et Artères and French National Agency of Research (ANR-MicroObese, and "Investissements d'avenir" ANR-10-IAHU-05).

ANR MicroObese consortium members

Verger E, Cotillard A, Kennedy SP, Prifti E, Pons N, Le Chatelier E, Almeida M, Quinquis B, Levenez F, Galleron N, Gougis S, Batto JM, Ehrlich D, Blottière H, Lepage P, Maguin E, Van de Guchte M,. Zucker JD.

Author contribution

JN conceived and designed the study. SS and PG developed the toolbox. SS and AM reconstructed and updated the GEMs. FB and PKD designed the in vitro experiments and PKD performed the in-vitro experiments, and analyzed the samples. KC designed the clinical study and coordinated sample analysis. JPF performed qPCR on fecal samples. KC and SR performed the clinical investigation. CJ and TW performed the fecal water extraction and EPG performed untargeted metabolomics on fecal samples. SS, PG, PKD and AM analyzed the data. PS performed in-silico analysis on untargeted metabolomics. JC, LH, ALN performed NMR, JC and MED assigned the NMR data, JC, LH JKN and MED analyzed the data. JN, FB, KC, AM and SS wrote the manuscript.

Competing financial interests

The authors declare no competing financial interests.


Paul Lerner

Olivier Ferret

Multimodal inverse cloze task

 1in answering questions about named entities grounded in a visual context using a textual Knowledge Base. Therefore, the interaction between the modalities in this task is paramount and must be captured with complex fusion models. As these models require a lot of training data, we design this pre-training task from existing work in textual Question Answering [START_REF] Lee | Latent Retrieval for Weakly Supervised Open Domain Question Answering[END_REF]. It consists in considering a sentence as a pseudo-question and its context as a pseudo-relevant passage and is extended by considering images near texts in multimodal documents. Our method is applicable to different neural network architectures and leads to a 9% relative-MRR and 15% relative-F1 gain for retrieval and reading comprehension, respectively, over a no-pre-training baseline.

"Fig. 1 :

 1 Fig. 1: Example of a question about a named entity with a visual context along with relevant passages from a Knowledge Base.

  Multimodal Inverse Cloze Task for Knowledge-based Visual Question Answering Paul Lerner 1 , Olivier Ferret 2 , and Camille Guinaudeau 1 1 Université Paris-Saclay, CNRS, LISN, 91400, Orsay, France 2 Université Paris-Saclay, CEA, List, F-91120, Palaiseau, France paul.lerner@lisn.upsaclay.fr We present a new pre-training method, Multimodal Inverse Cloze Task, for Knowledge-based Visual Question Answering about named Entities (KVQAE, [2]). KVQAE is a recently introduced task that consists, as illustrated by Figure