Valerie Girardin
email: valerie.girardin@unicaen.fr

Théo Grente
email: theo.grente@unicaen.fr

Nathalie Niquil
email: nathalie.niquil@unicaen.fr

Philippe Regnault
email: philippe.regnault@univ-reims.fr

Valérie Girardin

Comparing and Updating R Packages using MCMC Algorithms for Linear Inverse Modeling of Metabolic Networks

Keywords: Metabolic, trophic, biochemical and urban networks, Linear inverse modeling, High-dimensional polytope, Mirror and Billard walks, R package MSC Classification: 92-04 , 92-08 , 92-10 , 65C05 , 60J20

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

Introduction

According to Encyclopedia Britannica, metabolism is "the sum of the chemical reactions that take place within each cell of a living organism and that provide energy for vital processes and for synthesizing new organic material". However, metabolic systems cover many different concepts used in different domains.

The three main functions of metabolism in biology are: the conversion of energy in food into energy available to run cellular processes; the conversion of food into building blocks for proteins, lipids, nucleic acids, and some carbohydrates; the elimination of metabolic wastes. The development of molecular biology during the 1990's led to the characterization of first partial and then complete genomes in microbes, that opened a new area of research on microbial metabolic modeling; see [START_REF] Covert | Metabolic modeling of microbial strains in silico[END_REF]. The linear constraints on the flows within microbial cells are physico-chemical constraints -e.g. diffusion rates, topo-biological constraints -e.g. derived from the crowding of molecules inside cells, environmental constraints -e.g. nutrient availability, and regulatory constraints -e.g. enzyme regulation; see [START_REF] Price | Genome-scale models of microbial cells: evaluating the consequences of constraints[END_REF]. This type of model spread out of microbial biology with the development of overall sequencing of metazoan species and in particular of humans; see [START_REF] Duarte | Global reconstruction of the human metabolic network based on genomic and bibliomic data[END_REF]. The networks of metabolic flows within cells are called either metabolic networks in [START_REF] Bonarius | Flux analysis of underdetermined metabolic networks: the quest for the missing constraints[END_REF], metabolic reaction networks in [START_REF] Covert | Metabolic modeling of microbial strains in silico[END_REF] or biochemical networks in [START_REF] Orth | Reconstruction and use of microbial metabolic networks: the core Escherichia coli metabolic model as an educational guide[END_REF]; to avoid any confusion, we will here refer to them as the latter.

The Encyclopedia of Earth Science enlarges the definition: "Metabolism refers to the use of energy for the production and assimilation of food and for locomotion, maintenance, growth, reproduction, and other processes that characterize life". With this definition, the concept extends from biology to ecology, and allows the definition not only of cell or organism metabolism, but also of ecosystem metabolism. In ecosystems, the main metabolic processes are production and respiration. Because most of the production is ingested by consumers, the metabolic system can be represented as a network of flows from preys to predators, including waste production and respiration, that defines the food web.

The definition of metabolic networks also extends to economy, giving rise to the concept of urban metabolism, for studying material and energy exchanges in cities, with the metaphor of megapolis seen as bodies. The first discussions on urban metabolism go a long way back, to Karl Marx in 1881; see [START_REF] Zhang | Urban metabolism: A review of research methodologies[END_REF]. The concept came back in economic studies in the 1960's in the context of the growing pollution of air and water in USA cities. Since then, the literature on the subject developed and in 2008, a devoted congress highlighted urban metabolism as having a worldwide influence; see [START_REF] Havránek | ConAccount 2008: Urban Metabolism, Measuring the Ecological City[END_REF]. The study of urban metabolism facing strong environmental issues is particularly relevant in China, in a context of rapidly growing megapolis; see [START_REF] Zhang | Urban metabolism: A review of research methodologies[END_REF] and the references therein. A trend of studies merged economical approaches with concepts and methods derived from ecology; see, e.g. [START_REF] Chen | Information indices from ecological network analysis for urban metabolic system[END_REF]. Indeed, networks combining energy and matter flow circulation present common features with trophic networks. The other way round, ecological network analysis was strongly influenced by econometry with the application of Leontief input-output analysis, first in [START_REF] Hannon | The structure of ecosystems[END_REF]. These intertwined approaches of network analysis led in particular to the recent project UNCNET; see [START_REF] Guéret | Deliverable D2/3 of the UNCNET project: Using probability approaches to inform, revise, and improve contributions on the respective nitrogen flows[END_REF][START_REF] Kaltenegger | Urban nitrogen budgets: Evaluating and comparing the path of nitrogen through cities for improved management[END_REF].

In the three above types of metabolic networks, even if the goals may appear as different, flows can be represented by interactions, and the issuing mathematical models for the networks are graphs whose vertices are the compartments and edges the interactions. In order to specify the unknown flow values, a class of methods involving Linear Inverse Modeling (LIM) was developed; see [START_REF] Niquil | Inverse Modeling[END_REF] and the references therein. It relies on the principle of steady state mass balance, i.e., the sum of the inflows and outflows through the components of the system equals the rate of change in their standing stocks, most often considered negligible; the nature of the "mass" depends on the type of system. This yields a set of linear equations known as the mass balance equations (MBE). Through acquired knowledge, from dynamics field measurements, laboratory experiments, bibliographical knowledge, or simulation, equality constraints are added to the modeling. Moreover, bounds are imposed on the flows (biochemical networks) or linear combinations of flows are imposed to remain between bounds (urban and trophic). The set of these linear equations and inequations defines a convex bounded multidimensional polyhedron -a polytope, within which lie all realistic solutions to the problem.

Quite recently, the so-called sampling methods were developed to describe this kind of polytopes by calculating a representative sample of solutions through the Monte Carlo approach, and more efficiently the Monte Carlo Markov Chain approach (MCMC); see [START_REF] Kones | Are network indices robust indicators of food web functioning? A Monte Carlo approach[END_REF][START_REF] Leguerrier | Construction et étude d'un modèle de réseau trophique de la vasière de Brouage (bassin de Marennes Oléron[END_REF], and [START_REF] Van Oevelen | Quantifying Food Web Flows Using Linear Inverse Models[END_REF]. Thus, LIM-MCMC methods involve mass balanced models that consider the uncertainty in the data and link them pertinently to the variability of the living; see [START_REF] Grami | Functional Effects of Parasites on Food Web Properties during the Spring Diatom Bloom in Lake Pavin: A Linear Inverse Modeling Analysis[END_REF][START_REF] Niquil | Inverse Modeling[END_REF][START_REF] Van Oevelen | Quantifying Food Web Flows Using Linear Inverse Models[END_REF]. In all MCMC methods, an ergodic time-reversible Markov chain is designed, with the desired asymptotic distribution, leading to a uniform sampling of the polytope. Further, two reflective MCMC algorithms particularly suited to sampling highly complex LIM polytopes, the Mirror Walk (MiW) in [START_REF] Van Oevelen | Quantifying Food Web Flows Using Linear Inverse Models[END_REF] and of the Billard Walk (BiW) in [START_REF] Polyak | Billiard walk -a new sampling algorithm for control and optimization[END_REF], were developed independently. The former is designed in an R package called {limsolve} -including the function xsample()-especially intended for trophic modeling. The latter is implemented in [START_REF] Chalkis | {volesti}: Volume Approximation and Sampling for Convex Polytopes in R[END_REF] for sampling polytopes in order to estimate their volume, with an application to biochemical networks in [START_REF] Chalkis | Geometric algorithms for sampling the flux space of metabolic networks[END_REF]. Both include parameters, called jump lengths, to be fitted.

Recently, both [START_REF] Herrmann | Flux sampling is a powerful tool to study metabolism under changing environmental conditions[END_REF] and [START_REF] Fallahi | A comparison of Monte Carlo sampling methods for metabolic network models[END_REF] have undertaken rigorous comparisons of several sampling algorithms for biochemical networks, and concluded that the Coordinate Hit-and-Run with Rounding (CHRR) algorithm is the most efficient in terms of both computation time and multiple convergence diagnostics. Unfortunately, its only available implementation, the function chrrExpSampler() of the {COBRA} Matlab toolbox, works only for bounds on the flows, and not for more intricate inequality constraints, so cannot be applied to urban or trophic network models.

Regarding xsample(), [START_REF] Fallahi | A comparison of Monte Carlo sampling methods for metabolic network models[END_REF] was only able to successfully apply it to the Escherichia coli network, out of ten considered genome scale systems, and hence concluded that an efficient reflective MCMC algorithm implementation allowing for the framework of [START_REF] Van Oevelen | Quantifying Food Web Flows Using Linear Inverse Models[END_REF] to be applied at a large scale was lacking. One of the main goals of the present paper is to fill this gap. In this aim, we will describe and compare sampling methods, theoretically and with respect to computation time and statistical diagnostics of convergence. This will lead us to develop an updated version of the R package {limsolve} and function xsample(), called {samplelim} and rlim(), that uses the framework of the R package {volesti} designed in [START_REF] Chalkis | {volesti}: Volume Approximation and Sampling for Convex Polytopes in R[END_REF].

The paper is organized as follows. Section 2 is devoted to the construction of the polytopes of solutions to metabolic networks constraints, whichever be their type, with examples of real data trophic, urban and biochemical networks. Section 3 details different MCMC methods for sampling polytopes, from the classical Hit and Run and its variants, to the most recent reflective Hamiltonian MCMC methods, together with available implementations. The updated implementation {samplelim} is proposed in Section 3.4. For similar computation times, sampling performances of the functions of the R packages {samplelim}, {limsolve}, {volesti}, plus chrrExpSampler(), are thoroughly compared in Section 4, on the real data examples of Section 2. The methodology inspired by both [START_REF] Fallahi | A comparison of Monte Carlo sampling methods for metabolic network models[END_REF] and [START_REF] Herrmann | Flux sampling is a powerful tool to study metabolism under changing environmental conditions[END_REF], that also involves a new numeric index that we call range coverage, is presented in Section 4.1. In particular, Section 4.2 is devoted to show that {samplelim} is much faster than {limsolve}, while retaining its specific desirable features, which makes it a very competitive MCMC package for sampling LIM polytopes, especially for trophic networks.

From the constraints to the polytope of solutions

As highlighted above, from a mathematical point of view, a metabolic network is a valued oriented graph. After its theoretical construction from the physical or biological constraints, we will introduce realistic models for real data of the three different types of metabolic networks.

Let us call the graph (V, E, f ). Its vertices i ∈ V are the compartments, E ⊆ V ×V is the set of oriented edges with n = |E|, and f = (f ij ) ij∈E is the vector of flow values associated to edges. Precisely, ij = (i, j) ∈ E denotes the edge from compartment i to compartment j and f ij is the amounts of matter or energy transiting from i to j, the so-called flow -or flux. In many situations, all masses B i of the compartments can be measured with accuracy, while the flows f ij are much more difficult to evaluate. Therefore, we will adopt here the standard viewpoint that masses are given and flows are unknown.

The vector of flows f = (f ij ) ij∈E , has non negative components satisfying Kirchoff law on all compartments, yielding the MBE

j∈N + (i) f ij - j∈N -(i) f ji = 0, i ∈ V,
where the successors and predecessors of compartment i are N + (i) = {j ∈ V : ij ∈ E} and N -(i) = {j ∈ V : ji ∈ E}. Acquired knowledge on the network, coming from field measurements, laboratory experiments or any other source, may yield additional linear equations constraints on flows. All of them can be aggregated into

Af = b, (1)
where A is an m × n matrix and b an m-dimensional vector, with m the number of linear equations. Note that we are only interested in the under-determined problem, that is m < n: the knowledge on flows is not exhaustive, the solution is not unique.

Note that for all quantities to be well-defined, constraints for all the flows to be positive have to be added, say, f ij ≥ ε ij , for all i ∈ E, where the size of ε ij is to be fitted to the range of f ij . Other constraints imposed by acquired knowledge take the form of bounds on linear combinations of the flows. In mathematical words, this sums up to

Gf ≥ h, (2)
where G is a k × n matrix and h a k-dimensional vector, with k the number of such inequations. Apart from the MBE involving several flows in each equation, most constraints in trophic and urban networks involve only a small number of flows, and matrices A and G are mainly composed of 0. This is all the more true in the biochemical case where the inequality constraints are only bounds on the flows, so that G is of dimension 2n × n, with only one 1 per row corresponding to a lower bound and one -1 to an upper bound.

The set

S = {f = (f ij ) ij∈E ∈ R n : f ij ≥ ε ij , Af = b, Gf ≥ h}, (3)
of flows satisfying both constraints (1) and ( 2) is the intersection of hyperplanes and half-spaces so is a convex polyhedron. In addition, (1) and (2) always induce bounds on the flows, say

m ij ≤ f ij ≤ M ij , ij ∈ E, (4)
where

R ij = M ij -m ij is called the range of the flow f ij . Therefore S is a bounded polyhedron -called a polytope.
The polytope S is a subset of R n . Still, its dimension is n -ρ, where ρ is the rank of A. In other words, S = f * + A S, where f * ∈ S is any admissible flow, A is the n × (n -ρ)-matrix whose columns form a basis of the kernel of A, and

S = { f ∈ R n-ρ : G f ≥ h}, (5)
is the so-called reduced polytope, where G = G A and h = Gf * -h. Thus, up to an affine transformation, the polytope S with empty interior in R n is reduced to a full dimensional polytope S in R n-ρ . This proves to be a corner stone of most sampling methods.

We will conduct the comparison of packages on three different networks: a trophic network in the English Channel and its aggregated version, an urban network for nitrogen exchanges in Vienna and a biochemical network for Escherichia coli. Example 1. A trophic network of a marine ecosystem at Courseulles-sur-mer, France, in the English Channel, with 19 compartments and 144 flows from [START_REF] Noguès | Cumulative effects of marine renewable energy and climate change on ecosystem properties: Sensitivity of ecological network analysis[END_REF].

The ranges of the flows go from O(10 -5 ) to O( 102 ). The dimension of A is 25×144, corresponding to the 19 MBE plus 6 equations, each involving only one flow -so the values of these 6 flows are actually known. The dimension of G is 361 × 144. Here 144 rows of this matrix ensure positivity -each consists in only 0 but for one 1, while the other rows involve more than one flow and contain coefficients different from 1 or -1. Once reduced, the polytope is of dimension 119 and its ranges goes from O(10 -2 ) to O(10 2 ). Example 2. The trophic network of Example 1 aggregated according to a realistic ecological point of view into a network of 6 compartments and 28 flows; see [START_REF] Caputo | Analysis of trophic networks: an optimisation approach[END_REF], where all network constraints are detailed.

The ranges of the flows go from O(10 0 ) to O( 102 ). The dimension of A is 6 × 28, corresponding to 6 MBE, one for each compartment. The dimension of G is 72 × 28, of which 28 rows simply ensure positivity of the flows. The 44 other rows correspond to more complex constraints each involving more than one flow. Once reduced, the polytope is of dimension 22 with ranges going from O(10 0 ) to O(10 2 ). Example 3. An urban model of nitrogen exchanges in Vienna megapolis with 13 compartments and 69 flows; see [START_REF] Guéret | Deliverable D2/3 of the UNCNET project: Using probability approaches to inform, revise, and improve contributions on the respective nitrogen flows[END_REF] and [START_REF] Kaltenegger | Urban nitrogen budgets: Evaluating and comparing the path of nitrogen through cities for improved management[END_REF].

In [START_REF] Guéret | Deliverable D2/3 of the UNCNET project: Using probability approaches to inform, revise, and improve contributions on the respective nitrogen flows[END_REF], a LIM approach is considered to take into account the uncertainty on the flows measurements. In this network, the ranges of the flows go from O(10 2 ) to O(10 8 ), an artificial upper bound. The dimension of A is 16 × 69 with 13 MBE and 3 other equations, each involving several flows. The dimension of G is 302 × 69. Indeed, 69 inequations ensure that the flows are bounded by 10 8 , and 69 ensure that they are positive -note that 138 actually appear because of duplicates due to the way the model is declared. Then 91 inequations give more precise bounds on single flows, and 2 involve several flows. The reduced polytope is of dimension 53 with the ranges going from O(10 6 ) to O(10 9 ). Example 4. The core Escherichia coli metabolic model with 72 compartments (metabolites) and 95 flows; see [START_REF] Orth | Reconstruction and use of microbial metabolic networks: the core Escherichia coli metabolic model as an educational guide[END_REF]. The biochemical reactions are based on biochemical, genetic and genomic information, such as exact reaction stoichiometry, reaction reversibility, or relationships between genes and proteins.

The ranges of the flows go from O(10 0 ) to O(10 3 ), an artificial upper bound. The dimension of A is 72 × 95 with 72 MBE. All rows of G, of dimension 190 × 95, count only one non null coefficient, corresponding to bounds on the flows: 95 rows contain one -1 for upper bounds, while the other 95 contain one 1 for lower bounds. The rank of the matrix A is 71, once equations hidden in inequations are added to the constraints. These equations come from 8 flows with ranges 0. The reduced polytope is of dimension 24 with all ranges of order O(10 2 ).

MCMC for sampling polytopes

Despite the simplicity of its implicit definition in (3), the polytope of admissible flows is a complex, high dimensional, geometric object. Hence, its analytical description is an unreachable goal. A commonly used alternative consists in exhibiting sets of points taken at random in the polytope, obtained through the celebrated MCMC methods.

This section aims at presenting MCMC algorithms, together with some of their implementations, that are of common use in the context of metabolic network studies. An introductory discussion is proposed in Section 3.1 on the main points of attention to take care of when using MCMC methods. Then we focus on the Hit and Run algorithm and some of its variants in Section 3.2, and on two reflective Hamiltonian MCMC algorithms, the Mirror walk (MiW) and the Billard Walk (BiW), in Section 3.3. Section 3.4 will provide a unifying overview of three of these algorithms together with their most recent available implementations, the R packages {limsolve} and {volesti} and the Matlab toolbox {COBRA}, leading to the main innovation of this paper, an updated implementation of MiW and BiW, that we call {samplelim}. Its interface is in the R Statistical Software [START_REF] Core | R: A Language and Environment for Statistical Computing[END_REF], but computation is handled in the C++ programming language as in {volesti}, which significantly reduces computation time.

Note that all MCMC methods discussed here generate points inside a set with a non-empty interior. Hence, below, the "polytope" is, unless explicitly stated, the "reduced polytope" S in [START_REF] Chen | Information indices from ecological network analysis for urban metabolic system[END_REF] and not the original one of (3). For the sake of simplicity, it will still be denoted by S, with dimension n (instead of n -ρ).

MCMC uniform sampling of polytopes

First, basic Monte Carlo approaches have been considered in the ecological literature; see [START_REF] Kones | Gaining insight into food webs reconstructed by the inverse method[END_REF][START_REF] Leguerrier | Construction et étude d'un modèle de réseau trophique de la vasière de Brouage (bassin de Marennes Oléron[END_REF]. Unfortunately, they fail to generate efficiently uniform samples in high dimensional polytopes. Indeed, the ratio between the volume of S and the volume of the smallest hyper-rectangle ij∈E [m ij , M ij ] that contains S generally decreases exponentially fast with the number of dimensions, making the classical rejection method inefficient.

As an alternative, numerous MCMC methods have been introduced for generating samples drawn uniformly into a polytope. In all of them, an ergodic time-reversible Markov chain is designed, with the desired asymptotic distribution -here, the uniform distribution over the polytope. Thus a set of N 0 candidate points in the polytope is obtained by simulating the N 0 first realizations -draws-of the chain, called the pilot sample.

Such sets of points have two characteristics that are to be properly taken into account in order to get adequate results. Obviously, the first coordinates of the Markov chain generally do not follow its asymptotic distribution, so the first draws have to be discarded. The number M of such draws, called the burn-in period, intuitively represents the time to wait before convergence is quite achieved. Also, the successive coordinates of a Markov chain are correlated, so a certain number of draws have to be discarded after each drawn point before adding a next point to the returned sample. This number T is called the thinning parameter. Thus, the actual number of retained points from N 0 iterations of an MCMC algorithm is

N 1 ≃ (N 0 -M )/T. (6)
In other words, from any pilot sample

f (1:N0) = (f (1) , . . . , f (N0) ), (7)
with elements for each flow f

(1:N0) ij = (f (1) ij , . . . , f (N0) ij
), we obtain the actual burn-in free and thinned N 1 -sample in the polytope

(f (M +1) , f (M +T +1) , . . . , f (M +(N1-1)T +1) ). ( 8
)
All MCMC methods require a starting point for the chain. This first point is generally obtained by using a numerical solver. For instance, {limsolve} uses the LSEI (Least Squares with Equalities and Inequalities) algorithm of [START_REF] Haskell | An algorithm for linear least squares problems with equality and nonnegativity constraints[END_REF] to find a particular solution, whereas {volesti} computes the Chebyshev ball of the polytope and uses its center as the starting point.

Hit and Run and its variants

The symmetric mixing algorithm nowadays known as the Hit and Run algorithm (HR) was introduced by [START_REF] Smith | Efficient Monte-Carlo procedures for generating points uniformly distributed over bounded regions[END_REF]. The basic principle is to randomly choose both a direction to move forward (or backward) from a point of the polytope and a segment length to generate another one. The induced Markov chain formed by the iteration of this process is proven to converge to the uniform distribution in the polytope.

Precisely, in the HR procedure, from an initial point f (1) , points f (2) , . . . , f (N ) , are iteratively built by repeating the same successive actions, of which Algorithm 1 presents the pseudo-code:

1. Choose a random direction d i in the n dimensional space.

2. Determine the two intersection points I (i)-and I (i)+ between the line passing through f (i-1) and directed by d i , and the borders of the polytope.

Sample uniformly a point in the segment [I

(i)-, I (i)+ ].
Generating a random direction at Step 1 means simulating a realization of the uniform distribution on the unit sphere of R n . The main point here is that the distribution of the direction has to be symmetric (i.e., directions d and -d have the same probability), so that the transition kernel of the underlying Markov chain is symmetric.

When the polytope is very anisotropic, HR needs many iterations to achieve a uniform distribution of points. Indeed the convergence time or mixing time is of the order O(n 2 R 2 max /R 2 min ), where R max /R min , called the sandwiching ratio of the body, is the largest range of flows over the smallest. Examples 1 to 4 show that this ratio may be quite huge. Sampling directions utterly randomly leads to going into thin directions, thus inducing a long time to sample the large ones. Indeed, the algorithm spends a lot of time bouncing on closed borders of the polytope. This is a strong drawback of the method, somewhat addressed by variants.

The CHR, also called Gibbs Sampler, proposed in [START_REF] Turchin | On the computation of multidimensional integrals by the Monte-Carlo method[END_REF], is similar to HR, except for Step 1 to be replaced by:

1. Choose uniformly a random direction along one of the n dimensions.

The overhead per step for CHR with respect to HR is reduced from O(n 2 ) to O(n), since the update follows only coordinate directions. On the contrary, this leads to extremely correlated samples and induces a high thinning parameter.

Introduced in [START_REF] Haraldsdóttir | CHRR: coordinate hit-and-run with rounding for uniform sampling of constraint-based models[END_REF], CHRR begins by a pre-processing step of rounding of the polytope followed by a CHR exploration of the rounded polytope, before mapping samples back to the original polytope. The rounding step consists in computing the maximum volume ellipsoid inscribed in the polytope and applying the transformation that maps this ellipsoid to the unit ball. Note that a remarkable property of CHRR is its guaranteed distributional convergence to the uniform; see [START_REF] Laddha | Convergence of Gibbs sampling: coordinate Hitand-Run mixes fast[END_REF] for the mixing time.

More recent versions of HR are the artificially centered hit-and-run (ACHR), and its improved version optimized general parallel (OPTGP); see [START_REF] Herrmann | Flux sampling is a powerful tool to study metabolism under changing environmental conditions[END_REF]. Still, CHRR is found to perform the best for sampling LIM polytopes; see [START_REF] Fallahi | A comparison of Monte Carlo sampling methods for metabolic network models[END_REF] and [START_REF] Herrmann | Flux sampling is a powerful tool to study metabolism under changing environmental conditions[END_REF] and the references therein.

Reflective MCMC Algorithms

A well-documented pitfall of HR-like algorithms is the important auto-correlation of the samples; see [START_REF] Polyak | Billiard walk -a new sampling algorithm for control and optimization[END_REF] and the reference therein. This is particularly true for highly anisotropic LIM polytopes, with narrow angles inducing regions in which HR algorithms get stuck.

Reflective Hamiltonian MCMC algorithms rely on a reflection mechanism on the boundaries in order to avoid these traps and thus reduce auto-correlation; see [START_REF] Chevallier | Hamiltonian Monte Carlo with boundary reflections, and application to polytope volume calculations[END_REF] and the references therein. We will focus on two of them, the MiW of [START_REF] Van Oevelen | Quantifying Food Web Flows Using Linear Inverse Models[END_REF] in ecology, and the BiW in [START_REF] Polyak | Billiard walk -a new sampling algorithm for control and optimization[END_REF], independently designed but relying on a similar principle. The BiW Algorithm 1 Hit and Run algorithm pseudo-code Require: f (1) , an initial point in S. Ensure: f (1:N ) , a sequence of N points in S.

for i ∈ 2 : N do d i ← uniformly sampled in the unit sphere of R n .

λ (i)+ ← min{λ > 0 : f (i-1) + λd i ∈ Fr(S)}, I (i)+ = f (i-1) + λ (i)+ d i . λ (i)-← min{λ > 0 : f (i-1) -λd i ∈ Fr(S)}, I (i)-= f (i-1) -λ (i)-d i . f (i) ← uniformly sampled in [I (i)+ , I (i)-]
. end for Fig. 1: Reflection of the path in the MiW and BiW algorithms algorithm has been implemented in [START_REF] Chalkis | Geometric algorithms for sampling the flux space of metabolic networks[END_REF] for sampling polytopes in order to estimate their volumes.

As in HR, the general principle of MiW and BiW consists in drawing a direction and a path length. However, if the trajectory reaches a border of the polytope before the random path length is achieved, then it is reflected on this border, moving forward this reflected direction for the remaining distance. This is repeated until the trajectory finally reaches a point in the polytope; see Figure 1. Then, new direction and path length are randomly chosen, and so on. One major difference between HR and reflective algorithms is that the latter needs a parameter to be fitted, called the jump length of the algorithm.

Algorithm 2 Reflective Hamiltonian algorithms MiW and BiW pseudo-code

Require: f (1) , an initial point in S. Ensure: f (1:N ) , a sequence of N points in S for i ∈ 2 : N do d i ← sampled in unit sphere of R n . L i ← a random length for the trajectory. d ← d i (initial direction to follow) L ← L i (total length of the trajectory inside the polytope) f ← f (i-1)+dL while f ̸ ∈ S do λ * ← min{λ > 0 : f + dλ ∈ Fr(S)} z ← the normal vector to the hyperplan that reflects the trajectory

f * ← f + dλ * (intersection point in the hyperplan) d ← (d -2d T zz) (update direction to the reflected one) L ← (L -λ * ) (update the length of trajectory) f ← f * + dL (update candidate) end while f (i) ← f end for
The difference between MiW and BiW lies in the determination of both the direction d and the length L, also different from HR; see Algorithm 2 for details.

In the MiW of [START_REF] Van Oevelen | Quantifying Food Web Flows Using Linear Inverse Models[END_REF], a vector v is drawn from a centered non correlated Gaussian vector

V = (V 1 , . . . , V n ) with Var V i = σ 2 i . The vector σ 2 = (σ 2 1 , . . . , σ 2 n
) is called the jump length. Then d and L are obtained as the direction and norm of v. The distribution of the sampled points is stated in [START_REF] Van Oevelen | Quantifying Food Web Flows Using Linear Inverse Models[END_REF] to converge to the uniform distribution on the polytope.

It is worth noting for sound comparison with the BiW, that when σ i = σ for all i, then the squared norm ||V || 2 is σ 2 χ 2 (n) distributed, and hence its root L = ||V || has a Nakagami distribution with expectation and variance

EL = Γ((n + 1)/2) Γ(n/2) √ 2σ, VarL = σ 2 n -2 Γ((n + 1)/2) 2 Γ(n/2) 2 .
Thanks to the law of large numbers, L ∼ √ nσ for large n, and VarL converges to zero. In the BiW of [START_REF] Polyak | Billiard walk -a new sampling algorithm for control and optimization[END_REF], the direction is drawn from the uniform distribution on the unit sphere of R n and the path length is

L = τ log (1/U ), ( 9
)
where U is drawn from a uniform distribution on [0, 1], and the parameter τ > 0 is called the jump length too. Note that L has an exponential distribution with parameter 1/τ , with expectation τ and variance τ 2 , whatever be the dimension n of the polytope. The distribution of the sampled points is proven in [START_REF] Polyak | Billiard walk -a new sampling algorithm for control and optimization[END_REF] to converge to the uniform distribution on the polytope.

In both procedures, a major issue may occur depending on the choice of the jump length. Too small, a very large number of points would be needed to correctly sample the whole polytope. Too large, it would result into an important number of reflections and thus unnecessarily slow the process.

Implementation variants

According to the thorough studies [START_REF] Fallahi | A comparison of Monte Carlo sampling methods for metabolic network models[END_REF] and [START_REF] Herrmann | Flux sampling is a powerful tool to study metabolism under changing environmental conditions[END_REF], the CHRR is the best of HR variants for sampling LIM polytopes. So we have chosen to compare only CHRR to the reflective MiW and BiW. Available implementations of these algorithms, in various languages, generally addressed to a specific community, differ in their performances or outputs.

The implementation of CHRR from the MatLab {COBRA} toolbox [START_REF] Haraldsdóttir | CHRR: coordinate hit-and-run with rounding for uniform sampling of constraint-based models[END_REF] is, up to our knowledge, the only available one for sampling LIM polytopes. The function chrrExpSampler() takes a COBRA model structure as input and outputs the desired number of points. The number of iterations and the thinning parameter can be chosen. Unfortunately, this only works for upper and lower bounds on the flows, that is for biochemical models to which the implementation has been fitted. It does not work for inequations combining flows, which prevents its use for the more complex trophic or urban networks.

In contrast, [START_REF] Van Den Meersche | xsample(): An R Function for Sampling Linear Inverse Problems[END_REF] proposed, together with the MiW algorithm, an R package called {limsolve}, directed to sample LIM polytopes for trophic ecosystems. The function xsample() takes one point of the polytope as input, and yields a set of randomly sampled points in the polytope as output. Among other parameters, the number of iterations and the burn-in period can be chosen. Three different random walks are implemented: HR, CHR and MiW. The latter, that was new, is also the most commonly used. An associated jump length vector can also be chosen, with default value the vector with coordinates one tenth of each range of the reduced polytope, easily computed through the function xranges(). This R package also provides annex functions that are of special interest for researchers in ecology. In particular, the function setup() takes as input a declaration file describing the LIM constraints and outputs the matrix representation. The declaration files are written in a more natural language than matrices and vectors, which makes the package easy to use for practitioners, furthermore thanks to recent improvements in [START_REF] Gerber | Automated workflow for incorporation and evaluation of data uncertainty in ecological networks with autoLIMR[END_REF]. Unfortunately, {limsolve} is extremely slow, mainly because it is entirely coded in R. Note also that, according to our study, the function xranges() sometimes yields absurdly large ranges, that make the default jump length vector unreliable and lead to even longer computation times.

For BiW, we use the R package {volesti} of [START_REF] Chalkis | {volesti}: Volume Approximation and Sampling for Convex Polytopes in R[END_REF], whose prime purpose is to sample general polytopes in order to approximate their volume. Still it can be used for LIM polytopes for all types of metabolic models; see [START_REF] Chalkis | Geometric algorithms for sampling the flux space of metabolic networks[END_REF]. The language R is here only an easy-to-use interface, while the computation part is coded in C++. Since C++ is a low level programming language -a language whose functions are closer to the processor's instructions, this leads to way more efficient computation. The function sample points() takes as input the inequations defining the reduced polytope and the number of points to be sampled. Although polytopes may be defined through several representations, we will stick to the definition [START_REF] Chen | Information indices from ecological network analysis for urban metabolic system[END_REF], and hence to the function Hpolytope(). While CHR, HR, and BiW can all be chosen, the latter is the default when targeting a uniform distribution. The starting point of the chain, the thinning parameter -there called walk length, the burn-in period, and the jump length τ in [START_REF] Fallahi | A comparison of Monte Carlo sampling methods for metabolic network models[END_REF] are to be chosen. By default, the thinning parameter is set to 1, the burn-in period to 0, and the initial point f (0) is the center of the Chebyshev ball -the largest ball that fits inside the polytope-with radius r, obtained through the function inner ball(). The default value for τ is τ = 4 √ nr.

(10) Since it is of the order of the smallest range, this may lead to a poor exploration of the largest ranges, especially when the sandwiching ratio is large.

When analyzing the source code of the CRAN version of {volesti}1 , we uncovered two sticking points. First, the path length is said to be sampled as L = τ U , instead of equation [START_REF] Fallahi | A comparison of Monte Carlo sampling methods for metabolic network models[END_REF] as in [START_REF] Polyak | Billiard walk -a new sampling algorithm for control and optimization[END_REF]; this might be due to a typographical error since no justification is given for the change. Second, if the bound on the number of reflections typical of the BiW is reached, the point generated at the previous iteration is simply added to the returned sample before computing a new direction and length, leading to duplicates in the returned sample, thus clearly not uniform.

More recently, [START_REF] Chalkis | Geometric algorithms for sampling the flux space of metabolic networks[END_REF] proposed a multiphase Monte Carlo sampling algorithm addressed to the metabolic biochemical community. Its backbone is an accelerated variant of the BiW, performing both point and direction updates more efficiently by We have taken into account all the qualities of the above implementations in order to produce a more efficient one. Indeed the new package {samplelim} combines the performance of {volesti} and the convenient features of {limsolve}. It keeps their shared ability to handle all types of LIM polytopes, which the MatLab implementation chrrExpSampler() of CHRR fails to do. Technically, {samplelim} is built from a fork of {volesti} version 1.1.2-3, of which all functions non necessary for sampling LIM polytopes have been removed. On the one hand, {samplelim} uses the same C++ structure as {volesti} in order to perform efficient computations of the MiW inside the polytope. Moreover, it allows one to choose between BiW and MiW. On the other hand, it contains all the annex functions available in {limsolve} in order to facilitate its practical use. Also {samplelim} uses the {Rglpk} package to solve the linear systems leading to the computation of the theorical ranges of each flow; this solves the issue of xranges() of {limsolve} using {lpSolve} that sometimes returns absurd values. In the BiW version of {samplelim}, we have also rectified one of the two sticking points noted above in the CRAN version of {volesti}: the path length is sampled according to equation ( 9), with a logarithm as prescribed in [START_REF] Polyak | Billiard walk -a new sampling algorithm for control and optimization[END_REF].

For easy comparison, Table 1 summarizes for each of the four packages above: the random walk used for exploring the polytope; the programming language in which the package is developed; the addressed community, where Generic stands for all types of metabolic networks, and generic+t indicates functions specific to the trophic systems community on top of the generic framework; the polytope given as input to the sampling function, either the general polytope defined in (3) or the reduced one in [START_REF] Chen | Information indices from ecological network analysis for urban metabolic system[END_REF]; the paper associated to the package.

Comparison of implementations on real data

This section aims at comparing through simulation some implementations of the MCMC sampling algorithms described in Section 3. Precisely, we focus here on the comparison between the updated implementation of MiW and BiW provided in the R packages {limsolve} with the function xsample(), {samplelim} with rlim(), and {volesti} with sample points(), and with the CHRR implementation of the library {COBRA} of MatLab with the function chrrExpSampler().

Computation time and sampling performances are the two essential aspects of sampling functions to take into account for a sound comparison. Clearly, they do not evolve independently. Generically speaking, higher computation time yields better sampling performances, but their relationships take various forms, depending, in particular, on the polytope shape and on the parameters of the sampling functions -here the jump lengths of MiW and BiW. Exploring exhaustively these relationships for the three sampling functions would require a huge, prohibitive, number of simulations.

Instead, in Section 4.2, we first compare rlim() with xsample(). This comparison aims at showing that both these MiW implementations lead to the same quality of samples, but that rlim() is much faster. Since the very long computation time of {limsolve} hinders its use on large real data, the comparison is performed only on the aggregated trophic model of Example 2. The aim of Section 4.3 is to compare more thoroughly the implementations of MiW in {samplelim}, BiW in {volesti}, and CHRR in {COBRA}. Precisely, rlim() and sample_points() will be compared to chrrExpSampler() for sampling the polytope associated to the biochemical network of Example 4. Inequlity contrainsts of Examples 1 and 3 involve combinations of several flows, that are not supported by chrrExpSampler(). Hence, only rlim() and sample_points() will be applied in both examples, and their performances compared.

In Sections 4.2 and 4.3, the sampling performances will be investigated through a methodology inspired by [START_REF] Fallahi | A comparison of Monte Carlo sampling methods for metabolic network models[END_REF] and [START_REF] Herrmann | Flux sampling is a powerful tool to study metabolism under changing environmental conditions[END_REF] of convergence and uniformity statistical diagnostics, presented in Section 4.1: the Raftery-Lewis and Geweke diagnostics and the Effective Sample Size, and a new numeric index that we call Range Coverage.

When comparing computation time performances, the CPU-time is used. All samples have been generated using a computer with an Intel Core i5 processor (i5-8600K, 3.60GHz× 6) and 15.5Go RAM. The R session used version 3.10.3 of both BLAS and LAPACK.

Sampling performance diagnostics

The Raftery and Lewis (RL) and Geweke diagnostics, and the Effective Sample Size (ESS) are all diagnostics of convergence. They are applied to the sequences of draws, and performed separately for each flow, inducing some insight on the uniformity of the distribution. They are completed by the new index Range Coverage (RC) based on the ranges of the flows. Additional diagnostics can be found in [START_REF] Fallahi | A comparison of Monte Carlo sampling methods for metabolic network models[END_REF] and [START_REF] Herrmann | Flux sampling is a powerful tool to study metabolism under changing environmental conditions[END_REF], such as the Hellinger distance or the Interval Based Scale Reduction Factor. Several visual diagnostics are also used in [START_REF] Gerber | Automated workflow for incorporation and evaluation of data uncertainty in ecological networks with autoLIMR[END_REF], with comments on the ambiguity of their interpretation. The diagnostics that we selected will prove sufficient for discriminating the performances of the different implementations for the mere purpose of comparison.

The Raftery and Lewis diagnostics

The RL introduced in [START_REF] Raftery | The number of iterations, convergence diagnostics and generic Metropolis algorithms[END_REF] computes estimates of the thinning parameter, the burn-in period, and even the number of draws to be performed in an MCMC process. Its basis is that a correct sample of the objective distribution should give precise estimates of its quantiles.

For any pilot N 0 -sample f

(1:N0) ij
of a given flow f ij in [START_REF] Covert | Metabolic modeling of microbial strains in silico[END_REF] and T ∈ N * , the sequence S T defined by

S T,l = 1 ]-∞, q] f (1+T l) ij , l ≥ 0,
is a sequence of Bernoulli variables indicating whether each T -spaced sampled flow f

(1+T l) ij
is smaller than some quantile q or not. The RL diagnostics processes as follows: 1. Compute the empirical estimator q of the quantile q of f ij of order α ∈ (0, 1), where α is to be set by the user. 2. Let T be the smallest integer such that the amount of auto-correlation in S T is small enough to be overlooked; precisely, such that the penalized likelihood ratio between the order 2 and order 1 Markov models for (S T,l ) is negative. 3. Set M = T M ′ , where M ′ is the smallest integer such that the M ′ -th power P M ′ of the (empirical) transition matrix P of S T is close enough to its limit, whose both rows are (α, 1 -α). The closeness criteria is

∥e i P M ′ -(α, 1 -α)∥ < ϵ, for i = 1, 2
, where e 1 = (1, 0), e 2 = (0, 1), ∥ • ∥ is the euclidean norm on R 2 , and ε > 0 is a parameter to be set by the user. 4. Set N RL 0 = M + T N ′ , where N ′ is computed as follows. For any given N ′ , the central limit theorem for Markov chains applied to the vector (S T,l ) l=0,...,N ′ -1 yields that the empirical mean 1

N ′ N ′ -1 l=0 1 ]-∞, q] f (1+T l) ij
belongs to an interval [α -r, α + r] with probability p N ′ that increases to 1 as N ′ goes to infinity, where the value of r is to be set by the user. The value to be retained for N ′ is the smallest integer such that p N ′ is greater than some prescribed value p close to 1. The integers T and M returned by the RL diagnostics are estimates of the thinning parameter and burn-in period to be applied to the pilot sample. The integer N RL 0 is an estimate of the size of a pilot sample to draw for ensuring a good precision in estimating a prescribed quantile of the flows.

From a computational point of view, the parameters T , M , and N RL 0 are derived for all flows at once thanks to the function raftery.diag() of the R package {coda}. This requires five arguments: the pilot sample of points, the order α of the quantile to be estimated at first step, the tolerance level ε, the tolerance r required on the order of the quantile and the probability confidence level p. Note that the same function determines an estimate N min 0 of what would be the size of an idealized pilot sample, i.e., without correlation. It is determined as in Step 4 above, by assuming that (S T,l ) is an independent sequence instead of a Markov chain of order 1. We will not use this feature below, since it is not pertinent for comparison purposes.

Finally, the dependence factor

I = N RL 0 N min 0 ( 11
)
assesses the extent to which auto-correlation inflates the required sample size; values larger than 5 indicate a strong auto-correlation. In Sections 4.2 and 4.3 below, for the RL diagnostics, we have chosen the (default) parameters α = 0.025, ε = 0.001, r = 0.005, p = 0.95. Note that the RL diagnostics will here only be used to compare the quality of the different pilot samples through the values of M , T , N RL 0 and I in [START_REF] Gerber | Automated workflow for incorporation and evaluation of data uncertainty in ecological networks with autoLIMR[END_REF] and not to perform further sampling of size N 1 .

The Geweke diagnostics

The Geweke diagnostics introduced in [START_REF] Geweke | Evaluating the accuracy of sampling-based approaches to calculating posterior moments[END_REF] is a classical statistical test of comparison of means, applied to prescribed proportions of the first and last points generated by an MCMC algorithm, say 0 < p 1 < 1 and 0 < p 2 < 1, with p 1 + p 2 < 1. If the difference of means is too large, the null hypothesis of convergence is rejected.

For any pilot N 0 -sample f

(1:N0) ij
of a given flow f ij generated by an MCMC procedure,

Z = f (2) ij -f (1) ij se (1) ij + se (2) ij
is the standardized difference between the two empirical means, say f

(k) ij = 1 np k N p k l=1 f (l) ij for k = 1, 2,
of the first and last parts of the sample, where se

(1)

ij and se

(2) ij are the standard errors. Under the null hypothesis, Z is asymptotically normal. Thus, the convergence of the algorithm is rejected if |Z| exceeds a prescribed quantile of the normal distribution, usually 1.28 corresponding to a significance level of 0.2. In practice, several independent pilot samples are drawn and the number of times the null hypothesis that the draws are stationary is rejected is computed. The number of rejections over the number of pilot samples is theoretically close to 0.2, corresponding to the significance level of the test.

In Sections 4.2 and 4.3 below, parameters for the Geweke diagnostics have been set to p 1 = 0.1 and p 2 = 0.5, with quantile 1.28. A number of rejections of 2 out of 10 at most is expected for good sampling.

The Effective Sample Size

The effective sample size (ESS) of a correlated sample obtained through any procedure is the theoretical number of independent draws that would yield the same variance; for the MCMC procedure, see, e.g., [START_REF] Gelman | Bayesian Data Analysis[END_REF]. In mathematical words, for any N -sample f

(1:N ) ij of the flows, ESS(f (1:N ) ij ) = N 1 + 2 N -1 k=1 ρ k (f (1:N ) ij ) ,
where ρ k (f

(1:N ) ij ) = corr(f (1:(N -k)) ij , f ((1+k):N ) ij
) is the auto-correlation of order k of This score reflects the amount of auto-correlation of the MCMC draws: the closer the ESS is to N , the less the MCMC draws are auto-correlated. Up to our knowledge, no numeric bound is available in the literature to assess the quality of the score.

f (1:N ) ij . n rlim

The Range Coverage

The bounds m ij and M ij on the flows in ( 4) can be accurately computed from the matrices A et G by using a linear programming solver. This yields theoretical ranges

R ij = M ij -m ij for all the flows.
Here we propose to compute empirical estimates of these ranges from a sample f

(1:N ) ij
, and to compare them to the theoretical ones. Precisely, let us define the range coverages as

RC ij = M ij -m ij M ij -m ij , (12)
where

m ij = min l∈1:N f (l) ij and M ij = max l∈1:N f (l)
ij , for ij ∈ E. This produces a new index RC for assessing sampling performance. Quite obviously, the closer these rates are to 1, the wider is the percentage of area explored in the polytope.

Comparison of the functions rlim() and xsample()

We will show here that rlim() leads to the same quality of samples as xsample(), but much faster. Due to the very long computation time of xsample(), this comparison can only be conducted on relatively small models, here the aggregated trophic model of Example 2. The comparison is divided in two distinct parts, one to compare computation time and one to compare the quality of samples.

The first part consists in comparing the computation time2 needed to obtain pilot samples of size N 0 varying from 50 to 50 000, when using the default jump length, that is the same for both implementations. For each sample size, 10 replicates are computed. Figure 2 presents the mean computation time, expressed in seconds, associated to the 10 replicates for each sample size. Clearly, rlim() outperforms xsample() by a factor close to 10 for sample sizes greater than 1 000. This is essentially explained by the low-level optimized implementation of rlim() in C++ compared to the highlevel implementation of xsample() in direct R programming. The performance gap between C or C++ and R is a widely discussed topic on developers blogs, especially in the MCMC community; see, e.g. [START_REF] Kinlay | A comparison of Programming Languages[END_REF][START_REF] Wilkinson | Gibbs sampler in various languages[END_REF], where examples are detailed for which the outperformance factor goes from 50 to 300. Hybrid approaches using seamless C++ integration tools for R, as used in {samplelim}, offer substantial performance improvements while preserving the easy-to-use interface of R. Additional simulation, not included here, suggests that the speed difference between rlim() and xsample() increases with the size of the jump length and also with the complexity, dimension, and shape, of the polytope. They confirm comments in [START_REF] Fallahi | A comparison of Monte Carlo sampling methods for metabolic network models[END_REF], where the authors failed to use xsample() on most real-size metabolic biochemical systems, due to unreasonable computation time.

The second part of the comparison between xsample() and rlim() concerns sampling performances, based the diagnostics and indexes presented above in Section 4.1. 2, while the ESS and RC given in Figure 3 are the means of these indexes on the 10 samples.

These diagnostic outputs attest that both xsample() and rlim() yield MCMC samples of the polytope of Example 2 with similarly good convergence and autocorrelation properties. Small values for both the thinning parameter T and index I of RL diagnostics and the high values of ESS ratios for a large majority of flows indicate moderate auto-correlations. Small values of the burn-in period M , and the number of Geweke rejections close to the expected 2, that is 2.29 for xsample() and 2.39 for rlim(), indicate that convergence to the limit distribution is a statistically sound assumption even without discarding the first points of the pilot sample. The mean range coverages, with no RC lower than half the theorical range, indicate that an important part of the polytope has been explored. Further, these diagnostics highlight the very close behavior of xsample() and rlim(), since both lead to very similar performance results on all flows. The small residual differences may simply originate from the difference of the used pseudo-random number generating processes.

Comparison of the functions rlim(), sample points() and chrrExpSampler()

This section aims to compare the sampling performances of rlim() (MiW in {samplelim}) with the functions sample points() (BiW in {volesti}) and chrrExpSampler() (CHRR in {COBRA}). We will consider comparison on the basis of similar computation times: the parameters of the algorithms for the studied examples will be tuned in order for the simulation of a fixed number of points of the polytope to take the same time for the three functions. Then the quality of the samples will be compared. Since the CHRR implementation only allows for lower and upper bounds on the flows, its comparison will be possible only on the biochemical model of Example 4. First, the two others will be compared on the trophic model of Example 1 and the urban model of Example 3.

Both MiW and BiW depend on jump lengths, that have an important impact on the computation time and the quality of the sample. In sample points(), the default value of the jump length τ is given by [START_REF] Gelman | Bayesian Data Analysis[END_REF], while for rlim() the default jump length σ is one tenth of the vector of the reduced polytope ranges. Comparing results obtained for a vector parameter with respect to a scalar one turns out to be difficult to handle. Fortunately, for all σ = (σ 1 , . . . , σ n ) we have been able to determine numerically a scalar value σ such that (σ, . . . , σ) leads to a similar computation time for rlim(), with similar sampling performances; for the sake of shortness, these comparison results are not given here. This replacement questions the choice of a vector jump length in MiW with respect to a scalar one in BiW. Yet, since it is an open question to find a closed-form expression for σ as a function of σ, the easy-to-compute vector value is to be favored by practitioners. For the purpose of comparison, σ will below be replaced by the corresponding computed σ. In other words, we will consider dual values τ and σ leading to similar computation times for rlim() and sample points(), meaning that only the quality of the output samples will remain to be compared.

Further, MiW and BiW algorithms then differ mainly in the choice of the distribution of the path length L: Nakagami for MiW, and exponential for BiW; see Section 3.3. The path length value directly impacts the number of reflections of the trajectories, while the computation time is roughly a linear function of the number of reflections. Thus, we expected that, for a similar computation time, the path length mean would be essentially similar, and hence the samples produced by rlim() and sample points() would present similar qualities. Still, the sampling performances of the latter appear to be highly penalized by the included bound on the maximal number of reflections. Practically, for large values of τ , the probability that a trajectory exceeds the maximal number of reflections increases and may get close to 1; for these trajectories, the returned sample point is then the previous point. This often happens in the simulated sample, resulting in very poor sampling performances with respect to all statistical indexes.

The following three-step procedure will be applied for the comparison to highlight above comments: Preparation step. For each pertinent example and sampling function, we simulate 10 replicates of samples of moderate size 1 000, for 30 different values of τ and σ and record the corresponding computation times. This allows us to empirically identify values of τ and σ with similar computing times, that we will call dual values. Under the reasonable assumption that the computation time grows linearly with the number of points to sample, sampling larger sized samples for dual values will take similar times. Selection step. We compute three values σ dual to the following values of τ : the default value [START_REF] Gelman | Bayesian Data Analysis[END_REF] a value large enough to get a high proportion of trajectories whose number of reflections exceeds the upper bound. This value is selected empirically during the preparation step. Simulation step. For each of the couple dual values, we simulate by using each sampling functions 10 samples of size 20 000. For each replicate, the convergence and uniformity diagnostics of Section 4.1 are performed. The results are aggregated first by a mean of the 10 replicates, then either by mean or a maximum of the flows, according to the index. This results into a single score value for each index and each sampling function.

Table 3 and 4 and the four first rows of Table 5 show the aggregated diagnostics results for rlim() and sample points() for Examples 1, 3 and 4. The three first rows correspond to the three choices of dual values. The fourth one corresponds to large values of σ leading to better diagnostic results than the three dual values. The last row of Table 5 shows the comparison with chrrExpSampler() for Example 4, for which the jump lengths τ and σ of rlim() and sample points() are chosen to lead to similar computation times for the three functions.

Comparison of rlim() and sample points() on the trophic model of Example 1

This model is the more complex considered in this paper, with linear inequations involving several flows. The very small default value τ of sample points(), close to 7.8 × 10 -4 , leads to a poor sampling of the polytope, see Table 3. The default value σ of rlim() leads to very satisfactory samples but takes way too long to be computed, about 20 hours for 20 000 solutions.

Above a jump length treshold, precisely τ = 1 here according to the preparation step, the samples returned by sample points() show a drop in quality, a behavior that rlim() does not share.

Table 3 shows the default value τ = 7.8 × 10 -4 , the optimal value according to the statistical indexes τ = 0.7, and a value above the critical threshold, precisely τ = 100 for an easily noticeable effect. The dual jump lengths for rlim() chosen to yield the same computing times are respectively σ = 5 × 10 -5 , σ = 0.03 and σ = 0.2.

The sampling performances of sample points() for τ = 7.8×10 -4 and rlim() for σ = 5 × 10 -5 are not significantly different. Only the burn-in period M and the pilot sample size N RL 0 recommanded by the RL diagnostics really differ, both to the benefit of rlim(). The same is true for τ = 0.7 and σ = 0.03. In all cases, the results of the For τ = 100, the sampling performances of sample points() drop significantly, with much lower quality than rlim() with the dual σ = 0.2. The recommended thinning parameter T returned by the RL diagnostic of 1 for all flows is hard to interpret; it may result from a very large number of duplicates in the sample. Note that rlim() with the dual σ = 0.2 leads to way better sampling performances.

At the cost of a longer computation time, a greater jump length for rlim() always yields better sampling performances according to all indexes. See Table 3, where σ = 1.5 leads to the best results, but at the cost of a computation time of about 2.5 hours.

Comparison of rlim() and sample points() on the urban model of Example 3

Again, we will compare the sampling performances for jump lengths leading to similar computation times on this network of middle complexity. For sample points(), the default τ = 2 600, the optimal value of the preparation step according to all indexes τ = 10 6 , and a larger value leading to duplicates in the returned samples, say τ = 10 7 , lead to choose for rlim() respectively σ = 200, σ = 6 × 10 4 and σ = 2 × 10 5 . The comparison is summarized in Table 4, where the default jump length σ ≈ 3 × 10 6 of rlim() is also presented.

For the default jump length τ = 2 600 and the dual σ = 200, close values are obtained. Both samplings are of very poor quality with a high autocorrelation, I ≈ 80, a number of Geweke rejections of 7, low ESS and RC.

For τ = 10 6 , the autocorrelation drops to I = 16, and both ESS and RC are higher, with especially RC≈ 0.8. The results of the Geweke test are no better than with the default τ . On the other hand, for the dual σ = 6 × 10 4 , rlim() leads to a better mean ESS of 7 035 compared to 2 825 for sample points(), and a better number of Geweke rejections of 4 instead of 7.12, the value of the other indexes being similar.

For τ = 10 7 and σ = 2 × 10 5 , the value of all indexes are clearly better for rlim(). Again, note the hard to interpret low value of the thinning parameter T returned by the RL diagnostic on samples obtained with sample points(). According to all indexes, the sampling performances of sample points() are really bad: absurdly high I = 132, number of Geweke rejections close to 8, and very low ESS and RC. In contrast, rlim() ouputs samples of way better quality even if still not what is to be required.

Finally, for rlim() with the default jump length σ , the sampling performances are really good, and the computation time not unreasonable, about one hour for a sample No parameter alike a jump length is to be fitted for the CHRR random walk of chrrExpSampler(). Thus, in order to compare the sampling performance of the three implementations, we have chosen the jump lengths τ = 1 000 for sample points() and σ = 100 for rlim() leading to computation times similar to chrrExpSampler(). Both rlim() and sample points() clearly outperform chrrExpSampler(). The samples returned by the latter are much more auto-correlated with I = 80.8 versus I ≈ 1, show the highest number of Geweke rejections and very low ESS and RC. This may be due to the CHR exploration method where two following points in the chain only differ by one coordinate. Aside comparison with chrrExpSampler(), we have also compared rlim() and sample points(). Fitted for the same computation time, the two functions return samples of the same order of quality according to all indexes, as shown on rows 1 and 2 of Table 5, until the jump length of sample points() is higher than a threshold and leads to duplicates in the returned sample. Flow values out of their theoretical ranges even happen, and hence the RC for sample points() with τ = 50 000 does not appear on row 3 of Table 5; we fail to figure out a reason for this phenomena. Finally, as shown on row 4, rlim() with a large σ = 1 000 leads to a computation time close to three times longer than σ = 400, with only a negligible improvement in the quality of the returned samples : the jump length σ = 400 seems to be large enough to get the best possible results on this model.

Conclusion

In this paper, we designed a new R package called {samplelim}, and compared it with other existing computing packages for sampling polytopes in linear inverse modeling.

Its sampling function rlim() was compared with the functions xsample() of the R package {limsolve}, sample points() of the R package {volesti}, and chrrExpSampler() of the Matlab toolbox {COBRA}. The comparison was made both in terms of computation time and sample quality. To sum up the comparison results, rlim(): leads to the same quality of samples as xsample() but much faster; is more reliable than sample points() which has a threshold on the jump length leading to a drop in the quality of samples; is more efficient than chrrExpSampler() in terms of computation time, quality of the returned sample and, last but not least, of complexity of the concerned models. Besides the high qualities of its sampling function, the new package {samplelim} keeps the easy to use features of {limsolve}, as well as its annex functions that are important to practitioners.

Most of the time, when information is gathered on the simple elements of a trophic, biochemical or urban metabolic networks, the elements to be considered are very numerous. Then, when building the model network, this number has to be reduced, by aggregating some elements within the same compartment or edge of the network. This is mainly due to the very long computation time associated to the needed number of iterations for stabilizing the results. With the new R package {samplelim}, the very high gain of time and quality will allow the development of these models at a level of aggregation that will be closer to the study scale of the biological processes. Furthermore, coupling networks together may become possible for studying them in a spatial framework.

Fig. 2 :

 2 Fig. 2: Computation time of xsample() (dotted line) and rlim() (continuous line) in Example 2 (time in seconds).

Table 2 :

 2 Aggregated results of the RL diagnostics and Geweke number of rejections of xsample() and rlim() for 10 replicates of 20 000 points in the polytope of Example 2.

  Mean RC of the 28 flows

Fig. 3 :

 3 Fig. 3: Means of the 10 replicates of rlim() in blue and xsample() in green.The categories of compartments (on the left) for Example 2 are described in[START_REF] Caputo | Analysis of trophic networks: an optimisation approach[END_REF]

Table 5 :

 5 Comparison of sample points(), rlim() and chrrExpSampler() on the biochemical network of Example 4 of size 20 000. The RL index I is under the threshold of 5, the number of Geweke rejections close to the expected 2, and both ESS and RC values are satisfying.4.3.3 Comparison of rlim(), sample points(), chrrExpSampler() on the biochemical model of Example 4

Table 1 :

 1 The four considered implementations of MCMC algorithms storing computations from the previous iteration, and including a preprocessing step involving the normal vectors of the facets. A multiphase part is also added, that is a sequence of sampling phases, each leading to a rounding of the original polytope. The efficiency of BiW is thus improved from one phase to the next, and the convergence to the uniform distribution accelerates. Unfortunately, this new algorithm is not yet a part of the CRAN version of {volesti} and therefore we have not been able to consider it for comparison.

	Package	Function	Walk	Language	Community	Input	Ref.
	{limsolve}	xsample()	MiW	R	trophic systems	polytope	[35]
	{COBRA}	chrrExpSampler() CHRR	MatLab	biochemical systems	polytope	[16]
	{volesti}	sample points()	BiW	C++ with R interface	generic	reduced polytope	[3]
	{samplelim}	rlim()	MiW BiW	C++ with R interface	generic+t	polytope

  for the reduced polytope; the value with the highest global sampling performances in the preparation step;

	Algorithm	Jump	Comp.	M	N RL 0	T	I	# Geweke	ESS Range Cov.
	functions	lengths	times (s) (max)	(max) (mean) (mean) rej. (mean) (mean)	(mean)
	sample points() τ = 7.8 × 10 -4 (default)	1.46 13 803 1 115 920	14	78	7.98	139	0.06
	rlim()	σ = 5 × 10 -5	1.74	1 045	939 525	15	80	8.06	194	0.06
	sample points()	τ = 0.7	190	929 1 060 313	15	45	6.41	1 140	0.41
	rlim()	σ = 0.03	180	866	641 974	17	47	5.90	1 588	0.45
	sample points()	τ = 100	1 120 15 178 3 396 166	1	433	7.10	4	0.08
	rlim()	σ = 0.2	1 190	1 329	567 356	14	27	5.64	4 151	0.58
	rlim()	σ = 1.5	9 166	811	609 165	12	8	3.79	9 961	0.72

Table 3 :

 3 Comparison of sample points() and rlim() on the trophic network of Example 1

Table 4 :

 4 Comparison of sample points() and rlim() on the urban network of Example 3 tests are way under what is required for good samples, with in particular, a number of Geweke rejections of 6.

	Algorithm	Jump	Comp.	M	N RL 0	T	I	# Geweke	ESS Range Cov.
	functions	lengths	times (s) (max)	(max) (mean) (mean) rej. (mean) (mean)	(mean)
	sample points()	τ = 2 600 (default)	0.37	558	751 242	12	75	6.96	509	0.30
	rlim()	σ = 200	0.42	755 1 014 509	12	83	6.68	697	0.30
	sample points()	τ = 10 6	78	549	479 050	15	16	7.12	2 825	0.81
	rlim()	σ = 6 × 10 4	70	495	579 383	12	16	4.06	7 035	0.80
	sample points()	τ = 10 7	220	777 1 228 831	2	132	7.87	4	0.10
	rlim()	σ = 2 × 10 5	228	472	433 045	15	14	3.54	10 590	0.85
	rlim()	σ ≈ 3 × 10 6 (default)	3 703	25	31 601	6	2	2.20	16 182	0.98

Version: 1.1.2-7, Date:

2023-09-18

By computation time, we mean CPU-time, i.e., the amount of time the CPU cores of the computer specifically spend on the computations.

https://github.com/pregnault/samplelim Manuscript code: https://github.com/TheoGrente/limcomp


Clément Raguénès
email: clement.raguenes@insa-rennes.fr

Erwan Fourn
email: erwan.fourn@insa-rennes.fr

Cédric Quendo
email: cedric.quendo@univ-brest.fr

Rozenn Allanic
email: rozenn.allanic@univ-brest.fr

Denis Le Berre
email: denis.le-berre@univ-brest.fr

Tunable feeding network for beam steering applications based on chalcogenide material

Keywords: feeding network, tunable devices, beamforming, beam steering, chalcogenide glass, PCM, phase shifter

This abstract presents a new concept of tunable feeding network for antenna arrays. The objective is to illustrate the potential application of a Phase-Change Material (PCM) to antenna array. We will use chalcogenide thin layer to introduce phase shifters in a one-to-four feeding network. Hence, a wide range of states will be achievable, spanning from normal 0°r adiation to a 30°tilted radiation.

I. INTRODUCTION

With the ongoing increase in high data transfer requirements, operating frequencies tend to rise. As the frequency increases, path losses become more significant. To address this issue, radiation efficiency needs to be as optimal as possible. One solution, employed by the recent 5G specification, and futures ones [START_REF] Union | Future technology trends of terrestrial International Mobile Telecommunications systems towards 2030 and beyond[END_REF], is the use of multi-antenna arrays, unlike previous approaches. The main challenge of working with thin beams is to accurately target the reception point. One ability of antennas arrays is to change the radiation direction by applying a reversible phase law to the antennas. The reversibility can be done by making the feeding network reconfigurable. Also, with all the embedded systems constraints, it is interesting to study planar technologies. In this context, our study will work toward a solution to apply a phase law to an array, in a microstrip topology.

Usually, reconfigurable arrays are achieved via electronically tunable components, such as Micro-Electro-Mechanical Systems (MEMS), RF diodes or varactors [START_REF] Yu | Reconfigurable Bidirectional Beam-Steering Aperture With Transmitarray, Reflectarray, and Transmit-Reflect-Array Modes Switching[END_REF], [START_REF] Sundaram | MEMS-Based Electronically Steerable Antenna Array Fabricated Using PCB Technology[END_REF]. These components have already shown their advantages, and their limits, among those there are the electrical consumption, the commutation time, the high frequencies compatibility and the integration. Furthermore, managing each tunable cell's control has become a challenging task due to the growing size of This work has received a French government support granted to the CominLabs excellence laboratory and managed by the National Research Agency in the "Investing for the Future" program under reference ANR-10-LABX-07-01.

the arrays. This challenge is particularly pronounced because of the coupling between the electronic control circuit and the RF circuit [START_REF] Sundaram | MEMS-Based Electronically Steerable Antenna Array Fabricated Using PCB Technology[END_REF], [START_REF] Li | Broadband Electronically Scanned Reflectarray Antenna With Liquid Crystals[END_REF]. As a substitute for these electronic components, certain specialized materials have been under investigation for a number of years. These materials possess the capability to alter certain properties when exposed to a stimulus.

When the property is the electrical conductivity, we talk about Phase-Change Material (PCM). The two most studied ones are Vanadium Dioxide (VO 2 ) [START_REF] Huitema | Highly integrated VO 2 -based tunable antenna for millimeter-wave applications[END_REF] [6] and the chalcogenide family (mostly GeTe and Ge 2 Sb 2 Te 5 ) [START_REF] Mennai | Bistable RF switches using Ge2Sb2Te5 phase change material[END_REF] [START_REF] Singh | Characterization, Optimization, and Fabrication of Phase Change Material Germanium Telluride Based Miniaturized DC-67 GHz RF Switches[END_REF]. Used in planar topologies, these materials can entirely change the circuit's geometry by transitioning the conductivity value from that of a pseudo-insulator to that of a pseudo-conductor with a conductivity ratio of 10 4 to 10 6 .

The main difference between these competing materials is the transition temperatures and states stability. The VO 2 transition between the insulator and the conductive states occurs around 68°C, whereas for chalcogenide materials, temperatures of several hundred degrees are expected. However, to keep VO 2 in the conductive state, a heat source must be maintained, while for chalcogenide glasses, both states are stable. This implies that the source can be switched off once the transition is complete.

Additionally, the selection of the heat source type is crucial. There are two primary types of sources. First, it is possible to use an electrode beneath a PCM layer, heating it through the Joule effect [START_REF] Singh | Characterization, Optimization, and Fabrication of Phase Change Material Germanium Telluride Based Miniaturized DC-67 GHz RF Switches[END_REF]. However, this requires electronic control, which faces the same coupling challenges as electronic components. This is why directly heating the PCM layer with a laser diode can also be interesting. Two approaches can be envisioned: heating all layers with a single laser by diffracting its beam to radiate only specific areas, or utilizing multiple lasers via an optical fiber circuit. Both alternatives offer the advantage of minimizing coupling concerns as there would be no longer any electronic control circuit.

Given these advantages, we present here a proof of concept of a tunable feeding network based on chalcogenide thin films. First, the concept of the tunable phase shifters will be introduced. Then, the application of these phase shifters to the feeding network.

II. DEVELOPMENT

In this paper, our objective is to design a multi-state feeding network with four output ports. First, in the initial scenario, an equal phase is applied to each port. In this case, the radiation direction will be the same as it would be with a single antenna. Then, the targeted phase law will result in a -90°phase difference between each output. Such a phase law results in a 30°beam steering compared to the first scenario. However, other states could be achieved by applying different phase laws, with phase differences less than -90°. For instance, using the notation shown in Fig. 1, applying equal phase shifts to port 1/2 (-90°) and 3/4 (0°) would represent an intermediate state. To model and design our phase shifters, we use filtering theory. The model is based on typical known phase response of first order filters. These filters are known to imply a phase difference up to -90°between its accesses. The 0°to -90°v ariation is led by a special point which is corresponding to the cut-off frequency of the filter. At these specific frequency, the phase should be -45°.

Chalcogenide layers Aluminium lines

More specifically, we chose to use an RC filter for our model. The advantage of this topology lies in its simplicity and the relatively straightforward implementation of the capacitive component in a microstrip line. Indeed, the resistance can be readily associated with the losses of the transmission line, which are related to the characteristic impedance, while the capacitive component can be achieved through local widening of the line. When widening a transmission line, the capacitive behavior is favored relative to its inductive behavior.

To estimate the equivalent value C of the low impedance Z c , which is the impedance of the wide line, it is possible to use the following relation :

C = l Z c v ϕ
where l is the length of the capacitive line, Z c the characteristic low impedance of the line and v ϕ the phase velocity, which depends on the substrate characteristics (basically dielectric permittivity and thickness).

From the values of these two component R and C, the cut-off frequency of the model filter is given by : f 0 = 1/(2πRC). The phase variation ϕ is then described by the following relation :

ϕ = -arctan (ω/ω 0 )
Finally, the relation between the impedance of the capacitive line and the phase difference of the filter is :

Z c = - ωlR tan (ϕ) v ϕ
where ϕ is the phase shift from the RC filter at the pulsation ω = 2πf . For the approximation to be accurate, the length of the line l has to be small compared to the wavelength.

This relation makes us able to compute the impedance Z c to introduce in the transmission line to have the desired phase variation ϕ from the initial line. Then, the phase shift is realized through the introduction of an impedance step in an initial transmission line. The global phase response of this line is composed by the classic phase variation of the line plus the phase response of the impedance step, modeled by the RC filter.

However, this model filter is a low-pass filter. This means the transmitted wave magnitude over the cut-off frequency will decrease quickly. So, from the 0°to -90°phase variation interval, only the 0°to -45°part can be usable. Also, the magnitude level is -3 dB at f 0 , the real usable interval should even be smaller than [-45°, 0°].

The application of these phase shifters to antenna array feeding network is illustrated in the next point.

B. Realization

The widening is realized by deposition of thin chalcogenide layers, that have the ability to change its state, from a pseudoinsulator one to a pseudo-conductor. For every circuits that will be discussed thereafter, we consider a 350 µm thick silicon substrate with high resistivity (> 10 000 Ω.cm) and 1 µm aluminium thick lines. The chosen chalcogenide composition is the nominal composition Ge 2 Sb 2 Te 5 , deposited by RF magnetron sputtering with a thickness of 600 nm. Furthermore, to avoid any chalcogenide film breaks because of the level difference between silicon and aluminium lines, the whole surface is fully planarized by an insulating layer of SiO 2 in combination with dry etch back technics, before the deposition of the chalcogenide layer. The final structure is shown in Fig. 2. ig. 3: Phase response of the feeding network Our chalcogenide layers have a thickness of 600 nm and exhibit a crystalline electrical conductivity of approximately 10 5 S/m. However, this non-ideal conductivity introduces some transmission losses. These characterizations were performed using four-point probe measurements and retro-simulation techniques.

All the chalcogenide parts in the feeding network are shown in yellow in Fig. 1.

C. Results

To handle the initial goal of a 90°phase difference between each outputs, more than one filter have been necessary to implement by output as it can be seen in Fig. 1. By placing multiple impedance steps in series, the phase differences accumulate, allowing us to reach, at 30 GHz, respectively -86°, -84°and -80°between ports 1 and 2, 2 and 3, 3 and 4. The simulated response is shown in Fig. 3. The difference of a few degrees between the ideal -90°and the simulated phase differences is acceptable, considering that the -90°uniform phase differences result in a 30°steering angle, whereas it supposed to be 28°with the simulated values.

However, once the phase differences between the outputs are established, a discrepancy in magnitude becomes apparent. These discrepancies are evident in Fig. 4. The greater the number of steps, the less power is transmitted. It is also worth noting that each step is designed to possess a model cutoff frequency significantly higher than the 30 GHz operational frequency. This implies that the model filter is not intended to attenuate the signal at the operational frequency. Concurrently, the steps are likely to induce some parasitic reflections, contributing to the losses observed in Fig. 4.

Fig. 4: Magnitude response of the feeding network

To address or mitigate this concern, two approaches are currently being investigated. Firstly, introducing gradual transitions between the steps could potentially reduce the parasitic reflections and subsequently lower the losses. These types of gradual variations have been explored in greater depth within this paper [START_REF] Roy | The continuously varying transmission-line technique -application to filter design[END_REF].

Furthermore, some issues happened during the deposition of our circuits, disabling the chalcogenide layers to fully crystallize and reach a sufficient conductivity at the crystalline state. The main reason we study at this point is an unwanted reaction between the SiO 2 and the chalcogenide layers. Indeed, some kind of "bubbles" are appearing inside the chalcogenide layers only above SiO 2 , as Fig. 5 illustrates.

III. CONCLUSION

This paper presents a new type of phase shifter achieved through the use of phase-change material and a modeling method based on filter theory. This method uses the capability to modify the geometry of a transmission line to introduce a first-order RC filter. The phase-change material is employed to create a highly low-impedance region within the original transmission line. Subsequently, filter theory predicts the phase response of such a filter. By operating below its cutoff frequency, the filter exclusively adjusts its phase.

An application of these phase shifters is suggested. It aims to achieve a 30°beam steering through the application of a phase law in the feeding network. Given the objective to reach the -90°phase difference between the four outputs of the array, multiple phase shifter are needed.

Finally, the achievable phase differences between each successive couple of outputs are respectively -86°, -84°and -80°. From this perspective, the initial objective is achieved with an error of only 2°in the beam steering goal, but other issues have been highlighted. Indeed, the power distribution becomes highly uneven and prevents us from expecting a good radiation response. More developments are required to reduce this effect. Among the potential solution, using non uniform transmission line [START_REF] Roy | The continuously varying transmission-line technique -application to filter design[END_REF] to transform our brutal impedance steps to progressive impedance steps could help us to reduce the losses we observed at this point.

From the realization point of view, some studies have to be done before a real measurement, due to some issues between the materials used in this proof of concept.

Fig. 1 :

 1 Fig. 1: Top view of the simulated feeding network

Fig. 2 :

 2 Fig. 2: Illustration of the reason to use SiO 2 to avoid chalcogenide film breaks

Fig. 5 :

 5 Fig. 5: Bubbles apparition illustration.