Anna Zhukova
email: anna.zhukova@inria.fr

David James Sherman
email: david.sherman@inria.fr

Knowledge-based generalization of metabolic models

Keywords: metabolic modeling, generalization, genome-scale reconstruction

Background: Genome-scale metabolic model reconstruction is a complicated process beginning with (semi-)automatic inference of the reactions participating in the organism's metabolism, followed by many iterations of network analysis and improvement. Despite advances in automatic model inference and analysis tools, reconstruction may still miss some reactions or add erroneous ones. Consequently a human expert's analysis of the model will continue to play an important role in all the iterations of the reconstruction process. This analysis is hampered by the size of the genome-scale models (typically thousands of reactions), which makes it hard for a human to understand them.

Results: To aid human experts in curating and analyzing metabolic models, we have developed a method for knowledge-based generalization that provides a higher-level view of a metabolic model, masking its inessential details while presenting its essential structure. The method groups biochemical species in the model into semantically equivalent classes based on the ChEBI ontology, identifies reactions that become equivalent with respect to the generalized species, and factors those reactions into generalized reactions.

Conclusions: Generalization allows curators to quickly identify divergences from the expected structure of the model, such as alternative paths or missing reactions, that are the priority targets for further curation. We have applied our method to genome-scale yeast metabolic models and shown that it improves understanding by helping to identify both specificities and potential errors.

Introduction

Genome-scale metabolic models are complex networks that describe the thousands of reactions and molecular species that participate in an organism's metabolism. The complexity of these networks makes it difficult for human curators to understand, analyze, and verify them, since the individual reactions that require their attention are hidden among the reactions that are correctly described. The priority targets are gaps in pathways, and reactions that are in some sense unusual, because they represent shortcuts or alternative paths. Efficient curation of genome-scale models requires analysis and exploration tools that synthesize high-level views of the network and focus curator attention on these small subsets of unusual reactions.

Curation is performed after the automatic inference of the draft metabolic model. Inference methods combine databases of reactions and pathways with genomic information and existing models for similar organisms [START_REF] Swainston | The SuBliMinaL Toolbox: automating steps in the reconstruction of metabolic networks[END_REF]. Genomic data for the new organism is compared to the data of the reference organism, to find genomic evidence such as the presence of catalysing enzymes for the reactions conserved in the new organism. Starting from the inference of a draft model, the model refinement process includes several iterations of model analysis, error detection, and improvement [START_REF] Thiele | A protocol for generating a high-quality genome-scale metabolic reconstruction[END_REF].

The models produced at each iteration are intended for computer simulation, and so describe all the reactions thought to participate in the organism's metabolism. Although automatic model inference tools and genomic comparison methods are becoming steadily more sophisticated, they may still leave gaps in the model or add erroneous reactions. Curation by human experts is necessary.

Much of the complexity of the reaction network comes from biochemically similar reactions that operate on slightly different substrates. For example, in the peroxisome compartment of Yarrowia lypolitica model (MODEL1111190000 [START_REF] Loira | A genome-scale metabolic model of the lipid-accumulating yeast Yarrowia lipolytica[END_REF]) six acetyl-CoA oxidase reactions are present, transforming fatty acyl-CoAs differing in their carbon chain length (decanoyl-CoA, lauroyl-CoA, etc.) into the corresponding unsaturated fatty acyl-CoAs. There are also several similar reactions for other steps of the β-oxidation of fatty acids pathway [START_REF] Metzler | Biochemistry: The Chemical Reactions of Living Cells[END_REF]. Although all of these details are needed for accurate computer simulation, and are common to many models, not all of them are interesting for a curator. It is instead the differences from the common pattern that demand attention. They may be errors in the model, such as missing steps or erroneous connections between pathways, or they may be organism-specific differences such as alternative pathways that are biologically interesting.

To aid human understanding of genome-scale models, while keeping the details needed for a computer simulation, we propose a 3-level zoomable approach:

• The most abstract level represents compartmentalization of the model, and focusses on such questions as: Are all the compartments present? Are they well connected by transport reactions?

• The second level shows the modules inside of each of the compartments. The questions to be addressed on this level include: Are all the essential processes present? Is the structure of each process correct?

Is there any organism-specific adaptation of the structure?

• The most detailed level is intended for computer simulation and represents the inner structure of each of the modules with all the species, reactions and their kinetics, stoichiometry and constraints.

The two abstract levels are intended for a human expert, and the last one for the computer.

In this study we focus on the second level of abstraction, that represents the modules inside compartments.

A fair amount of work has been done on identifying reusable modules. These approaches can be divided into two groups: series and parallel.Aseries approach operates on chains of reactions, and generalizes them as a series, consequently hiding the structure of the network. An example of a series approach is representing the network as a set of metabolic pathways (KEGG [START_REF] Kanehisa | KEGG for integration and interpretation of large-scale molecular data sets[END_REF], MetaCyC [START_REF] Caspi | The MetaCyc database of metabolic pathways and enzymes and the BioCyc collection of pathway/genome databases[END_REF]), that can be further divided, for example, into reaction modules (conserved sequences of reactions along the metabolic pathways) [START_REF] Muto | Modular Architecture of Metabolic Pathways Revealed by Conserved Sequences of Reactions[END_REF].

The other type of approach operates on reactions that are parallel, keeping the steps and preserving the general view of the network. An example of this approach is grouping reactions based on EC (Enzyme Commission) numbers [START_REF] Tohsato | A multiple alignment algorithm for metabolic pathway analysis using enzyme hierarchy[END_REF]. The drawback of this approach is that it is not applicable to networks with no EC number assigned or reactions with no catalysing enzymes identified. We have developed another parallel -reaction method for knowledge-based generalization of metabolic models, which does not depend on enzyme information. It provides a higher-level view of a model while keeping its essential structure and omitting the details.

Mathematical basis

Basic definitions

We represent a metabolic model M as a pair of two sets: a set S of biochemical species, and a set R of reactions between them: We represent each reaction r ∈ R as a pair of sets of species: its reactants and products. A chemical reaction may be represented by a balanced chemical equation, showing the formulae of the reactants and products, and the changes that take place [START_REF] Clugston | Advanced Chemistry (Advanced Science)[END_REF]. This definition leads to restriction (1) that all the species participating in the reaction must be different.

M = S,
r = {s (rs) 1 ,...,s (rs)
k }, {s

,...,s

(ps) l } ∈ R ⊂ 2 S × 2 S ,
where s

(rs) 1 = ... = s (rs) k = s (ps) 1 = ... = s (ps) l (1)
To perform the model generalization, we define an equivalence operation ∼ on the species set, and group species into equivalence classes: [s] ∼ = {s ∈ S|s ∼ s}.

Species equivalence imposes reaction equivalence: two reactions are equivalent if their corresponding reactant and product species sets are pairwise equivalent.

∀r, r ∈ R r = {s (rs) 1 ,...,s (rs)
k }, {s

(ps) 1 ,...,s (ps) l } , r = {s (rs) 1 ,...,s (rs)
k }, {s

(ps) 1 ,...,s (ps) l } r ∼ r ⇐⇒ ∧ k = k, l = l, ∀i 0 ≤ i ≤ k ∃! ĩ 0 ≤ ĩ ≤ k : s (rs) i ∼ s(rs) ĩ , ∀j 0 ≤ j ≤ l ∃! j 0 ≤ j ≤ l : s (ps) j ∼ s(ps) j .
Equivalent reactions are factored together into a generalized reaction that operates with generalized species (i.e. species equivalence classes):

[r] ∼ = {[s (rs) 1 ] ∼ ,...,[s (rs) k ] ∼ }, {[s (ps) 1 ] ∼ ,...,[s (ps) l ] ∼ } .
In order to maintain the number of distinct species participating in a reaction, restriction (1'), analogous to restriction (1), must be satisfied:

[s

(rs) 1 ] ∼ = ... =[s (rs) k ] ∼ =[s (ps) 1 ] ∼ = ... =[s (ps) l ] ∼ (1 )
In order to avoid creation of paths in the generalized model that are not based on the evidence from the initial model, we introduce restriction (2): Species that do not participate in any pair of equivalent reactions and do not have any common equivalent species must not be grouped together.

∀s, s ∈ Ss∼ s ⇒∨ ∃r ∼ r ∈ R : s ∈ reactants(r) ∧ s ∈ reactants(r) ∃r ∼ r ∈ R : s ∈ products(r) ∧ s ∈ products(r) ∃ ṡ ∈ S : s ∼ ṡ ∧ ṡ ∼s.
(2)

Note that restriction (2) can be reformulated as maximizing the number of species equivalence classes while keeping the reaction equivalence classes unchanged.

The generalized model M/ ∼ is a pair of generalized species and reaction sets (quotient sets):

M/ ∼ = S/ ∼,R/ ∼ -generalized model, S/ ∼ = {[s 1 ] ∼ ,...,[s ñ] ∼ } -quotient species set, R/ ∼ = {[r 1 ] ∼ ,...,[r m] ∼ } -quotient reaction set.
The generalized model is a zoom out of the initial model. It provides a higher-level view by including less species and reactions, but more generic ones. For example, 3-oxodecanoyl-CoA, 3-oxolauroyl-CoA, and 3-oxohexanoyl-CoA species of the initial model can be generalized into oxo-fatty acyl-CoA.

Every reaction of the generalized model corresponds to at least one reaction of the initial model, having the same topology (number of distinct reactant and product species) and operating on species that can be zoomed out into those participating in the generalized reaction. The appropriate level of abstraction is defined with respect to the initial model as the most general one that satisfies restrictions (1') and ( 2).

The method and restrictions are described in figure 1.

Specific and ubiquitous species

We say that a ubiquitous species is one that participates in many reactions (more than some threshold), such as water, hydrogen, oxygen, etc. Grouping such species would increase the number of reactions each group would participate in, beyond the sharing already common to most of the models, and decrease readability;

in fact, during visualisation these species are often even duplicated to improve readability [START_REF] Rohn | VANTED v2: a framework for systems biology applications[END_REF]. Consequently we do not generalize ubiquitous species. In the generalized model each forms a trivial equivalence class:

S (ub) = {s (ub) 1 ,...,s (ub)
n }⊂S : ∀i [s

(ub) i ] ∼ = {s ub i }
Specific species are the others, which we divide into non-trivial equivalence classes and generalize accordingly.

Model generalization problem

Problem 1 Given a metabolic model M = S, S (ub) ⊂ S, R that describes n species (including n ≤ n ubiquitous ones) and m reactions, find an equivalence operation ∼ that obeys restrictions (1') and ( 2), and minimizes the number of reaction equivalence classes R/ ∼.

We will solve model generalization problem 1 in three steps:

1. Define the most general equivalence operation ∼ (having minimal number of species equivalence classes S/∼), that does not take into account the restrictions;

2. Modify the current equivalence operation to satisfy the restriction (1');

3. Modify the current equivalence operation to satisfy the restriction (2).

Step 1. Equivalence operation ∼.

Definition 2 Given a model M = S, S (ub) ⊂ S, R : S = n, S (ub) =n ≤ n, R = m, we define an equivalence operation ∼ on the species set S as forming n +1 equivalence classes in the quotient set S/∼: one for each of the ubiquitous species, and one for all the other species: ub) .

∀s (ub) ∈ S (ub) [s (ub) ] ∼ = {s (ub) }, ∀s, s ∈ S\S (ub) [s] ∼ =[s] ∼ = S\S (
Lemma 1 For any equivalence operation ∼ on the model M = S, S (ub) ⊂ S, R , the corresponding quotient species set S/ ∼ and quotient reaction set R/ ∼ are partitions of, respectively, the quotient species set S/∼ and the quotient reaction set R/∼ induced by ∼:

∀ equivalence operation ∼ defined on S, S (ub) ,R ∀s ∈ S [s] ∼ ⊂ [s] ∼ ∀r ∈ R [r] ∼ ⊂ [r] ∼ 9
n among them being ubiquitous, and m reactions.

Result:∼ -equivalence operation described in Lemma 1, M/∼ = S/∼,S (ub) (ub) ] ∼ ←{s (ub) } ; // map s (ub) to its equivalence class end S (ub) /∼←{[s (ub) ] ∼|s (ub) ∈ S (ub) };

/* Generalize specific species */ for s ∈ S\S (ub) do

[s] ∼ ← S\S (ub) ; // map s to its equivalence class end S/∼←S (ub) /∼∪{S\S (ub) };

/* Generalize reactions */ // map a reaction to its generalized version that operates with generalized species

gen ← λr. {[s] ∼|s ∈ reactants(r)}, {[s] ∼|s ∈ products(r)} ; for r ∈ R do [r] ∼ ←{r ∈ R|gen(r)=gen(r)}; end R/∼←{[r] ∼|r ∈ R}; return ∼, S/∼,S (ub) /∼,R/∼ Algorithm 1: Compute∼
The Compute∼ algorithm forms the equivalence classes for ubiquitous and then specific species as in Definition 2 and then computes the generalized reactions.

Step 2. Stoichiometry preserving property obedience

Algorithm

We start with the given equivalence operation ∼ 0 =∼, and iteratively improve it, until the stoichiometry preserving property (1') is obeyed. We denote the equivalence operation obtained at the i-th iteration step as ∼ i .

At each iteration, if there exists a species equivalence class that violates the stoichiometry preserving property (1'), i.e.:

∃s =s ∈ S, r ∈ R : s ∈ species(r) ∧ s ∈ species(r) ∧ [s] ∼ i =[s] ∼ i , we partition this species equivalence class [s] ∼ i =[ s] ∼ i into two: [s] ∼ i+1 ∨ [s] ∼ i+1 =[ s] ∼ i =[ s] ∼ i to form
a new approximation ∼ i+1 of the equivalence operation. When no species equivalence class violating the restriction (1') can be found, the current equivalence operation is returned as result.

At each iteration one equivalence species class is partitioned. In the worst case, the equality operation = (each species is equivalent only to itself) will be achieved. As it obeys restriction (1'), the process will terminate.

M/∼ = S/∼,S (ub) /∼⊂S/∼,R/∼ -corresponding generalized model.

Result

Species equivalence class partition

In a species equivalence class that violates the restriction (1') there are usually only a few conflicts present, and multiple solutions of the partition problem exist.

Species ontology

In order to make the choice of the species equivalence classes biologically meaningful, we use an ontology that describes hierarchical is a relationships (more specific to more general) between biochemical species.

Definition 3 A term t is a model term if it corresponds to a specific species in the metabolic model.

We assume that no two model terms are connected by a descendant-ancestor (more specific-more general) relationship in the ontology; otherwise, we mark the ancestor term ubiquitous:

∀t, T ∈ terms (∃ species(t),species(T ) ∈ S ∧ t ∈ descendants(T ) ⇒ t = T ).
We iteratively remove all the leaf terms that are not model terms from the ontology, so that all the model terms become leaves, and all the leaves become model terms.

For each species equivalence class that needs to be partitioned, we first find the least common ancestor T of the ontological terms corresponding to its species. If the ontology allows for multiple inheritance, and there are several such least common ancestors, we pick the first one. Then we look among the T -th descendant terms for those that are compatible (to avoid multiple inheritance). ∀r ∈ R {species(t),species( t)} ⊂ species(r).

(1 )

To do so, we first exclude all the terms that violate the stoichiometry preserving property (1"). We thus obtain an exact set cover problem.

Problem 4 (Set cover) Given a set X and a collection of its finite subsets Ψ, such that S∈Ψ S = X, find a minimum-size subset Π ⊂ Ψ whose members cover all of X: S∈Π S = S∈Ψ S = X.

Remark 1 Set cover is NP-complete [START_REF] Karp | Reducibility Among Combinatorial Problems[END_REF].

Problem 5 (Exact set cover) As in Set cover problem, except that here the sets used in the cover are not allowed to intersect.

Remark 2 Exact cover is NP-complete [START_REF] Goldreich | Computational Complexity: A Conceptual Perspective[END_REF].

Exact set cover applied to ontological terms

Each ontological term t defines a set S(t) of its descendant leaf terms (including t if it is a leaf). The instance consists of a set X of the model terms of interest, and a collection Ψ of all sets defined by their common ancestor T , its descendant terms, and their relative complements with respect to X: ∀S ∈ Ψ X\S ∈ Ψ, excluding all the sets that violate the stoichiometry preserving property (1"). We look for a minimum-size exact cover of X.

Note, that in this case an exact cover always exists, e.g. the one formed by all the leaf terms.

Choice of the ontology

We assume that any term that violates property (1") is removed from the ontology. Note that the term T is also removed.

If the ontology has no multiple inheritance, i.e. ∀S, S ∈ Ψ S ∩ S = ∅⇒S ⊆ S ∨ S ⊆ S, the problem becomes trivial: the set of the root terms forms the solution. The size of the solution, though, depends on the characteristics of the ontology, e.g. for a completely flat ontology (i.e., with no relationships) the solution consists of singleton equivalence classes.

If multiple inheritance is allowed, any Ψ ⊆ 2 X becomes possible, and the problem becomes NP -complete.

We use the ChEBI ontology (de Matos et al., 2010) of chemical compounds, as it is de facto a standard for species annotation in metabolic models. ChEBI consists of three main branches: chemical entity, role, and subatomic particle.T h echemical entity branch describes terms useful for annotation of biochemical species in a metabolic model.

The level of detail in the ChEBI hierarchy is not uniform: some sub-branches are more developed than others, so equally precise terms may be placed unequally deep in the hierarchical tree. For example, both hydrogen peroxide (CHEBI:16240) and decanoyl-CoA (CHEBI:28493) terms describe precise chemical molecules; but hydrogen peroxide is only 5 terms away from the chemical entity in the ChEBI hierarchy, while decanoyl-CoA is 11 terms away.

Besides that, different types of classification are combined together in the hierarchical tree, leading to multiple inheritance. For example, in the fatty-acid (CHEBI:35366) sub-branch, several classification types are present, including:

• classification based on the length of the carbon chain:

short-chain fatty acid (CHEBI:26666): 2-4 carbons;

medium-chain fatty acid (CHEBI:59554): 6-12 carbons;

etc.

• classification based on the presence of double bonds in the carbon chain:

saturated fatty acid (CHEBI:26607): no double bonds;

unsaturated fatty acid (CHEBI:27208): one or more double bonds;

• classification based on substituent groups:

hydroxy fatty acid (CHEBI:24654): one or more hydroxy substituents;

oxo fatty acid (CHEBI:59644): at least one aldehydic or ketonic group;

etc.

Moreover, using only hierarchical relationships in the ChEBI ontology is not always enough. Examples

show, that similar reactions can happen to the acid and the base in a conjugate acid-base pair. A conjugate acid-base pair is two species, one an acid and one a base, that differ from each other through the loss or gain of a proton [START_REF] Stoker | General, Organic, and Biological Chemistry[END_REF]. For instance, in the Rhea database of chemical reactions [START_REF] Alcántara | Rhea-a manually curated resource of biochemical reactions[END_REF] We extend Ψ so that it is closed under the operation of relative complement: ∀S, S ∈ Ψ S\ S ∈ Ψ. This allows for solving the set cover problem instead of the exact cover one: As Ψ is closed under the operation of complement intersection, we can obtain an exact set cover C from any set cover C = {S 1 ,S 2 ,...,S m } by replacing its elements with their relative complements with the previous elements

of C: C = {S 1 ,S 2 \S 1 ,...,S m \ m-1 i=1 S i }.
To approximate the solution of the set cover problem, we use a greedy algorithm.

Greedy Algorithm

Among the available subset candidates S i ∈ Ψ, pick the one of the largest size and add it to the resulting set cover Π. Repeat this operation until all elements of X are covered.

Data: X -set of interest, Ψ ⊆ 2 X -set of subsets of X Result:Π⊆ Ψ -set cover of X Π ←∅; // resultant cover while X = ∅ do // select S ∈ Ψ that covers maximum elements of X S (max) ← max(Ψ,criterion = λS. (S ∩ X)); Ψ ← Ψ\{S (max) }; X ← X\S (max) ; Π ← Π ∪{S (max) };
end return Π Algorithm 3: GreedySetCover Greedy set cover is a polynomial time approximation algorithm that achieves an approximation ratio of H( X), where H(n)i st h en-th harmonic number: [START_REF] Chvatal | A Greedy Heuristic for the Set-Covering Problem[END_REF]. It is the best possible polynomial time approximation algorithm for set cover, under plausible complexity assumptions [START_REF] Feige | A threshold of ln n for approximating set cover[END_REF].

H(n)= n i=1 1 i ≤ ln n +1
Step 3. Species equivalence class number maximization Problem 6 Given an equivalence operation ∼ defined on a metabolic model M = S, S (ub) ⊂ S, R ,fi n d an equivalence operation ∼ that obeys restriction (2) and does not change the reaction equivalence classes: R/ ∼= R/ ∼.

Algorithm

To satisfy restriction (2) we associate each species s in the initial model with a pair of sets of reaction equivalence classes in the quotient reaction set R/∼, induced by reactions where it participates as a reactant or product.

s → R (rs) s = {[r (rs) 1 ] ∼ ,...,[r (rs) o ] ∼ },R (ps) s = {[r (ps) 1 ] ∼ ,...,[r (ps) t ] ∼ } .
Definition 7 Given an equivalence operation ∼ defined on a metabolic model M = S, S (ub) ⊂ S, R ,w e define an equivalence operation ∼ as forming a separate species equivalence class for each of the ubiquitous species, and putting ∼-equivalent specific species that intersect in their product or reactant reaction classes in the same equivalence class:

∀s (ub) ∈ S (ub) ,s ∈ Ss (ub) ∼s ⇐⇒ s (ub) = s, ∀s, s ∈ S\S (ub) s ∼s ⇐⇒ ∧ s ∼ s (R (rs) s ∩ R (rs) s = ∅) ∨ (R (ps) s ∩ R (ps) s = ∅) ∨ (∃ ṡ ∈ S : s ∼ ṡ ∧ ṡ ∼s).
Any further partition of the quotient species set would imply the partition of the quotient reaction set.

Hence the number of species equivalence classes is maximal for the current number of reaction equivalence classes, and restriction (2) is satisfied.

Data: ∼ -equivalence operation defined on a metabolic model M = S, S (ub) ⊂ S, R , M/∼ = S/∼,S (ub) /∼⊂S/∼,R/∼ -corresponding generalized model.

Result

Applications

To illustrate the model generalization method we show its application to the genome-scale metabolic network of the lipid-accumulating yeast Yarrowia lipolytica (MODEL1111190000 [START_REF] Loira | A genome-scale metabolic model of the lipid-accumulating yeast Yarrowia lipolytica[END_REF]). The generalized model is attached as additional files 1, 2. The generalization has created 100 non-trivial quotient species and 217 non-trivial quotient reactions, and reduced the total number of species from 1847 to 1072 and of reactions from 2002 to 893.

The generalization method shows the best performance if the model is well-annotated with ChEBI. The species lacking ChEBI annotations are forced to form trivial quotient species in the generalized model as there is no evidence of their biochemical similarity with any other species in ChEBI. For 430 species in the Y. lipolytica model no appropriate ChEBI annotation was found, thus they could not be grouped with other species.

Peroxisome is an example of a well-annotated compartment in the Y. lipolytica model: only two species have no ChEBI annotations: YLR043C disulphide and TRX1. The generalization process reduced the number of reactions in peroxisome from 65 to 27. Figures 2 and 3 represent the peroxisome before and after the generalization and were produced using Tulip graph visualisation tool [START_REF] Auber | Tulip A Huge Graph Visualization Framework[END_REF].

The model before the grouping of equivalent reactions and species into generalized ones is shown on figure 2: different colors correspond to different equivalence classes. The same color code is used in figure 3 representing the generalized model that operates with quotient species and reactions. For example, the violet unsaturated FA-CoA node is a quotient of 8 species: hexadec-2-enoyl- CoA,and trans-tetradec-2-enoyl-CoA (colored violet in figure 2). In a similar manner, the light-green acCoA oxidase quotient reaction, that converts fatty acyl-CoA (yellow) into unsaturated FA-CoA (violet), generalizes 6 corresponding light-green reactions of the initial model (figure 2).

The generalized model describes the β-oxidation of fatty acids pathway [START_REF] Metzler | Biochemistry: The Chemical Reactions of Living Cells[END_REF] The more precise model is needed for simulation, while the more general one is clearer to a human, and reveals the main properties of the model. For example, the generalized model highlights the fact that there is a particularity concerning C24:0-CoA (tetracosanoyl-CoA) (red, inside the cycle in figure 3): there exists a "short-cut" reaction, producing it directly from another fatty acyl-CoA (yellow), avoiding the usual four-reaction beta-oxidation chain, used for other fatty acyl-CoAs.

Another application of model generalization is metabolic model comparison. The generalization brings the models to the same level of abstraction and highlights the differences such as gaps. Examples can be found in [START_REF] Zhukova | Knowledge-based generalization of metabolic networks: a practical study[END_REF].

Discussion

We have developed a method that provides a semantically zoomed-out view of a metabolic model, that keeps its essential structure but hides the details.

We have implemented our method as a Python program, that is available for download from http://metamogen.gforge.inria.fr. It takes an SBML model as an input, annotates its species with ChEBI terms (if the annotations are not present in the model) and generalizes it. It produces two SBML files, as an output. The first output file contains the generalized model. The second output file uses groups [START_REF] Hucka | Groups Proposal[END_REF] extension of SBML, and contains the initial model plus the groups representing quotient species and reaction sets.

We have applied our method to genome-scale metabolic models of yeasts. We have illustrated it here on the lipid-accumulating yeast Y. lipolytica and have shown that generalization helps finding gaps and peculiarities, as well as compresses the information stored in the model, which can be used for model visualisation and model comparison. In the example, the chain of β-oxidation of fatty acids reactions in the constitutive peroxisome of Y. lipolytica is generalized to a cycle of reactions, highlighting the alternative path for C24:0-CoA (tetracosanoyl-CoA).

Currently the generalization method depends on the ChEBI ontology. It cannot generalize species that lack ChEBI annotations. In future work we will overcome this limitation.

The method zooms out a model to the most general level of abstraction that is consisted with the model structure, i.e. does not violate the restrictions (1') and (2). It remains to be seen whether there are intermediate levels of abstraction that can be useful for model analysis. In particular it may be interesting to define the maximal generalization for a group of organisms, in order to highlight the specific differences of the individual models with respect to a common generalization.

  R -m o d e l , S = {s 1 ,...,s n } -s p e c i e ss e t , R = {r 1 ,...,r m } -reaction set.

Problem 2

 2 Given an equivalence operation ∼ defined on a metabolic model M = S, S (ub) ⊂ S, R find an equivalence operation ∼ that obeys restriction (1') and induces a quotient species set S/ ∼ of minimal size S/ ∼, such that S/ ∼ is a partition of the quotient species set S/∼ induced by ∼, i.e., ∀s ∈ S [s] ∼ ⊂ [s] ∼ .

  :∼ -equivalence operation described in Problem 3, M/∼ = S/ ∼,S (ub) / ∼⊂S/ ∼,R/∼corresponding generalized model. S/ ∼←S/ ∼ ; // resultant quotient species set S/ ∼⊂2 S S (ub) / ∼←S (ub) / ∼ ; // resultant quotient ubiquitous species set S (ub) / ∼⊂2 S (ub) R/ ∼←∅; // resultant quotient reaction set R/ ∼⊂2 R ∼ ←∼ ; // resultant equivalence operation ∼ : S ∪ R → S/ ∼∪R/ ∼ /* Partition quotient species that do not obey restriction (1') */ for S (gen) ∈{ S(gen) ∈ S/ ∼|∃s =s ∈ S(gen) ,r ∈ R : s ∈ species(r) ∧ s ∈ species(r)} do Π=Partition(S (gen) ); S/ ∼←Π ∪ S/ ∼\{S (gen) } ; // Update S/ ∼ for S(gen) ∈ Π do for s ∈ S(gen) do [s] ∼ ← S(gen) ; ← λr. {[s] ∼|s ∈ reactants(r)}, {[s] ∼|s ∈ products(r)} // map a reaction to its generalized version that operates with generalized species for r ∈ R do [r] ∼ ←{r ∈ R|gen(r)=gen(r)}; end R/ ∼←{∼(r)|r ∈ R}; return ∼, S/ ∼,S (ub) / ∼,R/∼ Algorithm 2: PreserveStoichiometry 13

  :∼ -equivalence operation described in Problem 2, M/∼ = S/ ∼,S (ub) / ∼⊂S/ ∼,R/∼corresponding generalized model. S/ ∼←∅; // resultant quotient species set S/ ∼⊂2 S S (ub) / ∼←S (ub) / ∼ ; // resultant quotient ubiquitous species set S (ub) / ∼⊂2 S (ub) R/ ∼←R/ ∼ ; // resultant quotient reaction set R/ ∼⊂2 R ∼ ←∼ ; // resultant equivalence operation ∼ : S ∪ R → S/ ∼∪R/ ∼ /* Update specific species generalization */ // Map a species to a set of its ∼-equivalent species that participate in ∼-equivalent reactions r sim ← λs.{s ∼ s|∃r, r ∈ R : s ∈ reactants(r) ∧ s ∈ reactants(r) ∧ r ∼ r}; p sim ← λs.{s ∼ s|∃r, r ∈ R : s ∈ products(r) ∧ s ∈ products(r) ∧ r ∼ r}; sim ← λs.r sim(s) ∪ p sim(s); S/ ∼←S (ub) / ∼∪{sim(s)|s ∈ S\S (ub) }; // Merge all quotient species sets that intersect while ∃S (gen) = S(gen) ∈ S/ ∼ : S (gen) ∩ S(gen) = ∅ do S/ ∼←(S/ ∼\{S (gen) , S(gen) }) ∪{S (gen) ∪ S(gen) }; end for S (gen) ∈ S/ ∼ do for s ∈ S (gen) do [s] ∼ ← S (gen) ; // map s to its equivalence class end end return ∼, S/ ∼,S (ub) / ∼,R/∼ Algorithm 4: Maximize

  happening inside the Y. lypolitica peroxisome in a generic way: as a transformation of fatty acyl-CoA (yellow) into unsaturated FA -C o A (violet), then into hydroxy FA-CoA (green), 3-oxo FA-CoA (magenta), and back to fatty acyl-CoA (with a shorter carbon chain); while the initial model describes the same process in more details, specifying those reactions for each of the fatty acyl-CoA species present in the organisms' cell (e.g. decanoyl-CoA, dodecanoyl-CoA, etc.). That is why the beta-oxidation chain of the reactions in the initial model, transforming step-by-step the fatty-acyl-CoA with the longest carbon chain into the one with the shortest chain, in the generalized model appears as a cycle (generalizing all the fatty-acyl-CoAs into one species, regardless the chain-length).

Figure 1 :

 1 Figure 1: Model generalization method and restrictions. a) Generalization first groups the species into equivalence classes, and then factors them into generalized species. The reaction equivalence classes and factoring are inferred from the species classes. b) Restriction (1'). The top part shows a correct generalization that obeys the restriction. Two bottom parts show generalizations that would change the reaction stoichiometry, and thus are not allowed. c) Restriction (2). The top part shows a correct generalization that obeys the restriction. The bottom part violates the restriction as there is no evidence in the model (i.e. no equivalent reaction) of the species b 3 belonging to the same equivalence class as b 1 and b 2 .

Figure 2 :

 2 Figure2: Peroxisome of the Y. lypolitica model (MODEL1111190000 ). Species are represented as circular nodes, and the reactions as square ones, connected by edges to their reactants/products, according to SBGN notation[START_REF] Moodie | Systems Biology Graphical Notation: Process Description language Level 1[END_REF]. Ubiquitous species are of smaller size and colored gray. Specific species are divided into six non-trivial equivalence classes, and colored accordingly (violet, light-blue, yellow, green, light-green, magenta). The specific species that form trivial equivalence classes are all colored red. Reactions are divided into fifteen non-trivial equivalence classes, also represented by different colors. Reactions that form trivial equivalence classes are all colored blue. The size of the model does not allow for readable species labels, so they are omitted.

Figure 3 :

 3 Figure 3: The generalization of the peroxisome of the Y. lypolitica model (see figure 2). The generalized model operates on quotient species and reactions. The number given in parentheses and the size of each node indicates how many entities it generalizes.

  Definition 1 The model generalization process groups chemical species present in the model into equivalence classes, and merges each class into a generalized chemical species. Reactions that involve same generalized chemical species are then factored together into a generalized reaction.By applying the model generalization process, we can build a simplified model that focusses on the high level relationships. The simplified model can be further divided into pathways.

  Definition 4 Terms t 1 ,...,t k are compatible if and only if their descendant model terms do not intersect:t 1 ,...,t k are compatible ⇐⇒ ∀ i = j ∈{1,...,k} descendants(t i ) ∩ descendants(t k ) ∩ leaves(T )=∅.Given a term T , find a compatible term set among its descendants, such that it has minimal size, covers all the T -th descendant leaf terms, and satisfies the stoichiometry preserving property (1"):

	Problem 3

? t 1 ,...,t k ∈ descendants(T ): ∧ k = k min , t 1 ,...,t k are compatible, leaves(T ) ⊂ descendants(t 1 ) ∪ ...∪ descendants(t k ), ∀i = j ∈{1,...,k}∀t ∈ leaves(t i ), t ∈ leaves(t j )

Acknowledgements

The authors would like to thank Drs. Romain Bourqui and Antoine Lambert of the LaBRI MABioVis team for advice on graph layout.

AZ was supported by a CORDI-S doctoral fellowship from Inria.

Complete algorithm

Data: M = S, S (ub) ⊂ S, R : S = n, S (ub) =n ≤ n, R = m -metabolic model describing n species, n among them being ubiquitous, and m reactions.

Result: ∼ -approximation of the equivalence operation described in Problem 0, M/ ∼= S/ ∼,S (ub) / ∼⊂ S/ ∼,R/ ∼ -corresponding generalized model.

Algorithm 5: Compute∼

Author Disclosure Statement

No competing financial interests exist.

Additional Files

Additional file 1 -The generalized SBML level 2 version 4 model of Yarrowia lipolytica (derived from MODEL1111190000 )

The model contains generalized reactions and species.

Additional file 2 -The SBML level 3 version 1 model of Yarrowia lipolytica (derived from MODEL1111190000 ) with groups extension representing equivalent species and reactions

The model contains all the elements (reactions, species, etc.) of the initial MODEL1111190000 model, and is enriched with groups representing species and reaction equivalence information.


Xiaolin Zeng

How Vertex reinforced jump process arises naturally

.

Introduction

One of the most remarkable results in probabilistic symmetries is the de Finetti's theorem [START_REF] De | Funzione caratteristica di un fenomeno aleatorio[END_REF], which states that the law of any exchangeable sequence valued in a finite state space is in fact a mixture of i.i.d. sequences. This theorem has a geometrical interpretation via Choquet's theorem. More precisely, the subspace of exchangeable probabilities forms a convex, and those probabilities given by i.i.d. sequences are exactly the extreme points of the convex [START_REF] Aldous | Exchangeability and related topics[END_REF].

In the 1920s, W.E. Johnson [START_REF] Sandy | We johnson's" sufficientness" postulate[END_REF] conjectured that, under some technical conditions, if a process X n is exchangeable and P(X n+1 = i|X 0 , . . . , X n ) depends only on the number of times i occurs and the total steps n, then X n is nothing but the famous Polya urn: drawing balls uniformly from an urn and put back one additional ball with same color as the drawn one. This is a process with linear reinforcement. In term of random walk, the natural counterpart of Polya urn is the edge reinforced random walk (ERRW). Diaconis conjectured that this process have the same characterization as Polya urn. In [START_REF] Silke | How edge-reinforced random walk arises naturally[END_REF] S.W.W.Rolles have shown that both conjectures are true under technical conditions.

The vertex reinforced jump process (VRJP) is a linearly reinforced process in continuous time. In a recent paper, Sabot and Tarres [START_REF] Sabot | Edge-reinforced random walk, vertexreinforced jump process and the supersymmetric hyperbolic sigma model[END_REF] have shown that ERRW is a mixture of VRJP, which indicates that the VRJP are building blocks of ERRW, thus should share a similar characterization. This paper gives this characterization (Theorem 4), as a counterpart of Rolles' result; namely, the only continuous time process which is partially exchangeable and the transition probability depends only on neighbor local times is VRJP, under technical conditions.

Let us first recall the definition of ERRW, let G = (V, E) be a locally finite undirected graph without direct loops (edges with one endpoint). Let Z n denote the location of the random process at time n. Let a e > 0, e ∈ E. For n ∈ N, define w n (e), the weight of edge e at time n, by if {Z n , v} ∈ E 0 otherwise. Now let us introduce some definitions before stating Rolles' result. Again G = (V, E) is a locally finite undirected graph without direct loops, with its vertex set V and edge set E. Denote i ∼ j if {i, j} ∈ E. Following Rolles, we call (Z n ) n≥0 a nearest neighbor random walk on G, if it is a discrete time random process (not necessarily Markov) such that successive positions are neighbors.

An admissible path of the random walk is a sequence of vertices of G, denoted π = (v 0 , v 1 , . . . , v n ) such that consecutive vertices are neighbors. The number of visits to vertex i of path π is denoted

N i (π) := #{k : v k = i, k = 0, . . . , n};
Similarly, the number of transition counts in the path π of an oriented edge e = (i, j) is denoted

N e (π) = N i,j (π) := #{k : v k = i, v k+1 = j, k = 0, . . . , n -1}.
Two paths ξ, η are said to be equivalent and denoted ξ ∼ η, if ξ and η start at the same state and the transition counts from i to j of any pair (i, j) are equal for ξ and η, i.e. N i,j (ξ) = N i,j (η) for all (i, j).

Remarks 1. Two equivalent paths necessarily end at the same vertex.

Definition 1. A nearest neighbor random walk is partially exchangeable if any two equivalent paths have the same probability.

Theorem 1 (Diaconis & Freedman [START_REF] Diaconis | de finetti's theorem for markov chains[END_REF]). Let Z n be a recurrent random walk (i.e. with probability one it returns to Z 0 infinitely often), then Z is a mixture of Markov chains if and only if it is partially exchangeable. Moreover, the mixing measure is uniquely determined.

As it turns out that edge reinforced random walk is a mixture of reversible Markov chains, Rolles introduced the following more restrictive notion of partial exchangeability: for π = (v 0 , . . . , v n ) and e = (i, j) let Ñe (π) := #{k :

v k = i, v k+1 = j or v k = j, v k+1 = i, k = 0, . . . , n -1}.
Definition 2. A nearest neighbor random walk is partially exchangeable in a reversible sense if it satisfies the following: for any two paths ξ, η, if Ñe (ξ) = Ñe (η) for all e ∈ E, then ξ and η have the same probability.

In [START_REF] Silke | How edge-reinforced random walk arises naturally[END_REF] Theorem 1.1, Rolles proved that if a nearest neighbor random walk is recurrent and partially exchangeable in a reversible sense, then it is a mixture of reversible Markov chain.

Rolles' main result in [START_REF] Silke | How edge-reinforced random walk arises naturally[END_REF] states that, if G = (V, E) is a strongly connected graph and Z n is a nearest neighbor random walk on G such that the following assumptions are satisfied:

1. Z is partially exchangeable in a reversible sense (Definition 2).

2. For all v ∈ V and e ∈ E there exists a function f v,e taking values in [0, 1] such that for all n ≥ 0

P(Z n+1 = v|F n ) = f Zn,e (N Zn (Z 0 , . . . , Z n ), ÑZn,v (Z 0 , . . . , Z n )).
Then Z is an edge reinforced random walk or a Markov chain under some technical conditions (c.f. [START_REF] Silke | How edge-reinforced random walk arises naturally[END_REF] for precision).

Next we define the vertex reinforced jump process X t . Assign positive weights (W e ) e∈E to the edges, the process X t starts at time 0 at some vertex i 0 , if X is at vertex i ∈ V at time t, then, conditioned on the past, the process jumps to a neighbor j of i with rate W i,j (1 + l j (t)), where for e = {i, j}, W i,j = W e and l j (t) is the local time of vertex j at time t:

l j (t) := t 0 1 Xs=j ds.
Theorem 2 (Sabot & Tarres [START_REF] Sabot | Edge-reinforced random walk, vertexreinforced jump process and the supersymmetric hyperbolic sigma model[END_REF]). The ERRW Z n with weights (a e ) is equal in law to the discrete time process associated with a VRJP X t in random independent weights W e ∼ Gamma(a e , 1) And finally, the VRJP X t turns out to be partially exchangeable within a time scale (c.f. next section for the definition of partial exchangeability in continuous times). Let

D(s) = i∈V (l i (s) 2 + 2l i (s)),
then the process Y t = X D -1 (t) is a mixture of Markov processes with an explicit mixing measure, in addition, the mixing measure turns out to be related to a σ-model introduced by Zirnbauer, c.f. [START_REF] Sabot | Edge-reinforced random walk, vertexreinforced jump process and the supersymmetric hyperbolic sigma model[END_REF] Theorem 2.

In this paper we give a counterpart of Rolles' result for VRJP, namely we characterize exchangeable jump processes with local rate functions.

Definitions and results

Definition 3. We call (X t ) t≥0 a nearest neighbor jump process on G, if it is a random process which is right continuous without explosion, and each jump is from some vertex i to one of its neighbors j (i.e. i ∼ j). Definition 4. A nearest neighbor jump process X t is a mixture of Markov jump processes if there exists a probability measure µ on Markov jump processes such that L(X t ) = L(Y t )µ(dY ), where L denotes the law of respective processes. If for µ a.s. the Markov processes are reversible, then the process X t is a mixture of reversible Markov processes.

Freedman introduced the notion of partial exchangeability in continuous time in [START_REF] Da Freedman | De finetti's theorem in continuous time[END_REF].

Definition 5 (Freedman). A continuous process X t is partially exchangeable if for each h > 0, the law of {X nh ; n = 1, 2, • • • } satisfies the following property: for any two paths ξ = (ξ 0 , . . . , ξ l ), η = (η 0 , . . . , η l ) such that ξ ∼ η , P(X 0 = ξ 0 , . . . , X lh = ξ l ) = P(X 0 = η 0 , . . . , X lh = η l ).

We recall the de Finetti's theorem in continuous time showed by Freedman [START_REF] Da Freedman | De finetti's theorem in continuous time[END_REF]. Theorem 3. Let X t be a continuous time process starting at i 0 ∈ G, X t is mixture of Markov jump processes if 1. X t has no fixed points of discontinuity, more precisely, for every t, if t n → t, then P(X tn → X t ) = 1;

2. X t is recurrent;

3. X t is partially exchangeable.

Our main theorem is:

Theorem 4. Let X t be a nearest neighbor jump process on G satisfying the following assumptions:

1. For all i ∈ V , there exists C 2 diffeomorphisms h i such that X is partially exchangeable within the time scale D(s) = i∈V h i (l i (s));

2. G is strongly connected (i.e. any two adjacent vertices are in a cycle);

3. The process, at vertex i at time t, jumps to a neighbor j of i with rate f i,j (l j (t)) for some continuous functions f i,j

Then X is a vertex reinforced jump process within time scale, i.e. there exists another time scale D such that X D-1 (t) is a vertex reinforced jump process.

Remarks 2. In fact, the hypothesis of Theorem 4 implies that the functions f i,j (x) are necessarily of the form W i,j x + ϕ j .

Remarks 3. Note that we do not a priori require f i,j = f j,i , i.e. there is no assumption of reversibility for X t ; however the VRJP is a mixture of reversible Markov jump processes within time change.

Remarks 4. Concerning the third assumption, we cannot prove the result with rate f i,j (l i , l j ), but the case where f i,j (l i , l j ) = f i (l i )f j (l j ) can be treated. In fact, by applying a time change, the process with rate function of the form f i (l i )f j (l j ) can be reduced to our theorem.

In section 3, we introduce an equivalent notion of partial exchangeability and, as an example, we give a different proof of partial exchangeability of VRJP within a time scale. Section 4 contains the proof of Theorem 4.

3 The two notions of partial exchangeability

Partial exchangeability, infinitesimal point of view

Consider a nearest neighbor jump process on G satisfying the third assumption of Theorem 4. As we have assumed regularity on the trajectory of the process (c.f. Definition 3), to describe the law of our process, it is enough to describe the probability of the following events:

σ = {X [0,t 1 [ = i 0 , X [t 1 ,t 2 [ = i 1 , X [t 2 ,t 3 [ = i 2 , . . . , X [t n-1 ,tn[ = i n-1 , X [tn,t] = i n },
which will be denoted

σ : i 0 t 1 -→ i 1 t 2 -t 1 ---→ i 2 . . . i n-1 tn-t n-1 -----→ i n t-tn ---→
in the sequel and we call such an event a trajectory.

It turns out that when the jump rate is a continuous function of local times, the law of our process can be characterized by some function, which will be called density in the sequel. In fact, for the study of certain history depending random processes, we have the following lemma: Lemma 1. If (X t ) is a jump process with jump rate depending only on local times and the current position of the random walker, i.e. there exists functions f i,j (l) such that conditioned on the past, X t jumps from i to j at rate f i,j (l(t)), and, moreover, f i,j (l(t)) does not depend on the variable l i (t). Then there exists functions d σ , such that for all bounded measurable functions Φ defined on the trajectories,

E(Φ(X u , u ≤ t)) = n≥1 i 0 ,...,in d σ Φ(σ)dt + d i 0 t - → Φ(i 0 t - →)
where d σ = exp(-

t 0 j∼Xs f Xs,j (l(s))ds) n k=1 f i k-1 ,i k (l(t k )) and d i 0 t - → = P(X s = i 0 , 0 ≤ s ≤ t).
Remarks 5. We believe that Lemma 1 still hold when f i,j (l) depends on l i (t). In fact, if we can find a time changed process such that its jump rates do not depend on l i (t), it is immediate by re-applying the inverse time change that Lemma 1 holds in the general cases.

Proof. As f i,j (l(t)) does not depend on l i (t), the holding time of X t at i is exponentially distributed with rate

j∼i f i,j (l(t))
and the probability of jumping from i to j is

p(i, j) := f i,j (l(t)) k∼i f i,k (l(t))
.

Moreover, the process up to time t is characterized by the events

i 0 s 1 -→ i 1 s 2 -→ • • • sn -→ i n s n+1 ---→, s 1 , . . . , s n+1 > 0, n+1 i=1 s i ≤ t. For 1 ≤ k ≤ n + 1, denote t k = s 1 + • • • + s k , P(X t follows the trajectory i 0 s 1 -→ i 1 s 2 -→ • • • sn -→ i n s n+1 ---→, s k > 0, s k ≤ t) = tn≤t n k=1   p(i k-1 , i k ) exp( j∼i k-1 f i k-1 ,j (l(t k-1 ))s k ) • j∼i k-1 f i k-1 ,j (l(t k-1 ))   P(s n+1 > t -t n )ds = t 1 <t 2 <•••<tn<t exp(- t 0 j∼Xs f Xs,j (l(s))ds) n k=1 f i k-1 ,i k (l(t k-1 ))dt, with ds = ds 1 • • • ds n , dt = dt 1 • • • dt n .
Now the lemma follows by distinguishing different trajectories.

Definition 6. We say that X t admits a density if the assumptions in Lemma 1 are satisfied, and we denote its density as d σ .

Let us now give another definition of partial exchangeability for continuous time processes in terms of density. Define two trajectories σ and τ to be equivalent and denoted σ ∼ τ , if their discrete chain strings are equivalent and the local times are equal at each vertex. Formally,

Definition 7. Let σ = i 0 t 1 -→ i 1 t 2 -t 1 ---→ i 2 • • • i n-1 tn-t n-1 -----→ i n t-tn ---→, τ = j 0 s 1 -→ j 1 s 2 -s 1 ----→ j 2 • • • j n-1 sn-s n-1 -----→ j n t-sn ---→ .
Then σ and τ are equivalent if and only if

∀i ∈ V, l σ i (t) = l τ i (t) ∀i, j N i,j (σ) = N i,j (τ ).
where N i,j (σ) denotes the number of jumps from i to j in σ, i.e. N i,j (σ) = N i,j ((i 0 , . . . , i n )), and l σ i (t) = Proof. Suppose that the process X t is partially exchangeable in density, let h > 0, consider the event I = {X 0 = i 0 , X h = i 1 , . . . , X nh = i n }, let (j 0 = i 0 , j 1 , . . . , j n ) be an equivalent string of (i 0 , . . . , i n ), and J = {X 0 = j 0 , X h = j 1 , . . . , X nh = j n }. We construct a bijection T which maps trajectories of I to those of J. As (i 0 , . . . , i n ), (j 0 , . . . , j n ) are equivalent, for any pair of neighbors (i, j), there are exactly a same number of transition counts from i to j. Let us define T to be the transformation which is a permutation of the time segmentations [lh, (l + 1)h) of size h; which, for any k, moves the kth transition i kth --→ j of I to the kth transition i kth --→ j of J, and leaving the last time segmentation [nh, ∞) invariant. Figure 1 illustrates an example of such application.

0 1 0 2 1 0 2 1 0 1 Figure 1: The transformation T for I = {X 0 = 0, X h = 1, X 2h = 0, X 3h = 2, X 4h = 1} and J = {X 0 = 0, X h = 2, X 2h = 1, X 3h = 0, X 4h = 1}. Let σ = k 0 s 1 -→ k 1 s 2 -→ k 2 • • • k N -1 s N --→ k N s N +1 ---→
be one trajectory of the event I, we check that

T (σ) = k ′ 0 s ′ 1 -→ k ′ 1 s ′ 2 -→ k ′ 2 • • • k ′ N -1 s ′ N --→ k ′ N s ′ N +1 ---→
is a trajectory of the event J, and that T is one-one and on-to (c.f. Figure 2). If we fix the total number of jumps N and the discrete trajectory (k 0 , k 1 , . . . , k N ), then T can be though as a substitution of integration. Thus

0 1 0 2 1 0 2 1 0 1 σ T (σ) u1 u2 u3 u4 u5 u6 u7 u8 u9 u1 u2 u3 u4 u5 u6 u7 u8 u9 s 1 s 2 s 3 s 4 s 5 s 6 s ′ 1 s ′ 2 s ′ 3 s ′ 4 s ′ 5 s ′ 6
P(I) = N k 0 ,k 1 ,...,k N 1 s 1 ,...,s N +1 ∈I(N,k 0 ,...,k N ) d σ ds 1 • • • ds N +1 = N k ′ 0 ,k ′ 1 ,...,k ′ N 1 s ′ 1 ,...,s ′ N +1 ∈I ′ (N,k ′ 0 ,...,k ′ N ) d T (σ) ds ′ 1 • • • ds ′ N +1 = P(J),
where I(N, k 0 , . . . , k N ) is the subset of R N +1 defined as the set of (s 1 , . . . , s N +1 ) such that the event k 0

s 1 -→ k 1 s 2 -→ • • • k N s N +1
---→ is in I; and I ′ (N, k ′ 0 , . . . , k ′ N ) is its image by applying T ; see Figure 2 for a concrete example. As T preserves local times and the numbers of transition counts, these two integrals are whence equal.

Proposition 2. If a jump process is partially exchangeable in the sense of Definition 5, and its jump rate is a continuous function of local times, then it is also partially exchangeable in the sense of Definition 8.

Proof. Let X t denote such a process, for h > 0, consider the σ-algebra

F h = σ(X nh , n ≥ 0), let F 0 = σ(∪ h>0 F h ) and F = σ(X t , t ≥ 0).
As in [START_REF] Da Freedman | De finetti's theorem in continuous time[END_REF], we only consider h running through the binary rationals. Note that F 0 = F thanks to the right continuity of the trajectories.

Let

σ = i 0 t 1 -→ i 1 t 2 -t 1 ---→ i 2 • • • i n t-tn
---→ be a trajectory with n jumps (say n ≥ 1 to avoid triviality). Let {X (h) ∼ σ/h} denotes the event

{X 0 = σ 0 , X h = σ h , . . . , X N h = σ N h , with N = ⌊t/h⌋}.

It turns out that

d σ = lim h→0 P(X (h) ∼ σ/h)h -n .
In fact, let Ψ = 1 X (h) ∼σ/h , by definition of d σ ,

E(Ψ(X u , u ≤ t)) = P(X (h) ∼ σ/h) = k≥1 i 1 ,...,i k d τ Ψ(τ )dt 1 • • • dt k (1)
where

τ = i 0 t 1 -→ i 1 t 2 -t 1 ---→ i 2 • • • i k-1 t k -t k-1 -----→ i k t-t k ---→ .
When h is small enough, the sum in (1) must be over k ≥ n, and we have

P(X (h) ∼ σ/h) = P 1 + P 2 .
where for some p k , k = 1, . . . , n depending on h P 1 = P((X u ) 0≤u≤t makes n jumps at times s 1 , . . . , s n with s k ∈ (p k h, (p k + 1)h] and the trajectory is i 0 , . . . , i n ) P 2 = P((X u ) 0≤u≤t makes more than n + 1 jumps and X (h) ∼ σ/h)

Note that the jump rates are bounded from both below and above, and any holding time in the event of P 2 must be in an interval of length lesser than 2h, whence the probability of making n + l (l ≥ 1) jumps following the trajectory σ/h is smaller than the probability of n + l independent exponential variables (of parameter C) each smaller than 2h, where C is an upper bound of the jump rates. Whence

P 2 ≤ l≥1 (P(cst ≤ Exp(C) < cst + 2h)) n+l ≤ l≥1 (P(Exp(C) < 2h)) n+l = O(h n+1 ).
Thus P 2 can be dropped when taking the limit. In addition,

P 1 = (pn+1)h pnh • • • (p 1 +1)h p 1 h d σ dt 1 • • • dt n ,
note that here d σ depends only on t 1 , . . . , t n and it is an absolutely integrable function, by Lebesgue differentiation theorem (Theorem 1.6.19 [START_REF] Tao | An introduction to measure theory[END_REF]) lim h→0 P 1 /h n = d σ . Now let

d σ d τ P(X ih ∼ σ/h)h -n P(X ih ∼ τ /h)h -n
σ ∼ τ , when h is sufficiently small, proceeding as in the diagram shows that d σ = d τ .

Example: VRJP is partially exchangeable within a time change

Recall that Y s = X D -1 (s) with D(s) = i∈V (l i (s) 2 + 2l i (s)), It turns out that we can write down the density of the trajectory σ of the (time changed) VRJP process Y (For convenience, write s n+1 for s in the sequel). The density of

σ := i 0 s 1 -→ i 1 s 2 -s 1 ----→ i 2 • • • i n-1 sn-s n-1 -----→ i n s-sn ---→ is (c.f. [11]), denoting S i (t) = t 0 1 Yu=i du the local time of Y , d σ =( 1 2 ) n n k=1 W i k-1 ,i k i∈V,i =in 1 1 + S i (s) • exp (- i∼j W i,j ( (S i (s) + 1)(S j (s) + 1) -1)), (2)
which clearly depends only on final local times and transition counts, thus by Proposition 1, Y is partially exchangeable. On finite graph it is rather easy to prove that the VRJP is recurrent (for example, using a representation of VRJP by time changed Poisson point process as in [START_REF] Sabot | Edge-reinforced random walk, vertexreinforced jump process and the supersymmetric hyperbolic sigma model[END_REF], and then use an argument as in [START_REF] Davis | Reinforced random walk[END_REF] or [START_REF] Sellke | Reinforced random walk on the d-dimensional integer lattice[END_REF]). Therefore, Y is a mixture of Markov jump processes.

For convenient, we include a proof of this in the sequel (after the proof of Proposition 3), since the mechanisms of the proof enlightens the proof of the main theorem.

Proof of Theorem 4.1 Computation of densities

Let X be a nearest neighbor jump process on G satisfying the assumptions of Theorem 4, in particular, recall the time scale

D(s) = i∈V h i (l i (s)). ( 3
)
Let l i (t) be the local time of the process X at vertex i at time t. Let us denote the process after time change to be

Y t = X D -1 (t) , (4)
let S i (s) = s 0 1 Yu=i du (5)
denote the local time of Y . Consider the trajectory

σ : i 0 t 1 -→ i 1 t 2 -t 1 ---→ i 2 • • • i n-1 tn-t n-1 -----→ i n t-tn ---→ (6)
where 0 < t 1 < • • • < t n < t, after applying the time change, the corresponding trajec-

tory for Y is σ Y : i 0 s 1 -→ i 1 s 2 -s 1 ----→ i 2 • • • i n-1 sn-s n-1 -----→ i n s-sn ---→
where s k = D(t k ).

Proposition 3. With the same settings as in equations (3) (4) (5) ( 6), the density of the trajectory σ Y for Y is

d Y σ = exp   - s 0 j∼Yv f Yv,j (h -1 j (S j (v))) h ′ Yv (h -1 Yv (S Yv (v))) dv   n k=1 f i k-1 ,i k (h -1 i k (S i k (s k-1 ))) h ′ i k-1 (h -1 i k-1 (S i k-1 (s k )))
.

Proof. Applying Lemma 1 to the process X,

d σ = exp   - t 0 j∼Xu f Xu,j (l j (u))du   n k=1 f i k-1 ,i k (l i k (t k-1 )).
Recall that in (3) we assumed that h i : R + → R + are diffeomorphisms satisfying h i (0) = 0.

Next we compute the density for the same trajectory σ but for the process Y s = X D -1 (s) , as we have S i (D(s)) = h i (l i (s)), derivation leads to

S i (D(s)) ′ = D ′ (s)1 Y D(s) =i = h ′ i (l i (s))1 Xs=i . Hence (D -1 (t)) ′ = 1 D ′ (D -1 (t)) = 1 h ′ Y t -• h -1 Y t -(S Y t -(t)) , l i k (t k-1 ) = h -1 i k (S i k (D(t k-1 ))) = h -1 Ys k (S Ys k (s k-1 )).
Substituting s = D(t), we have

d Y σ = exp (- s 0 j∼Yv f Yv,j (h -1 j (S j (v))) h ′ Yv (h -1 Yv (S Yv (v))) dv) n k=1 f i k-1 ,i k (h -1 i k (S i k (s k-1 ))) h ′ i k-1 (h -1 i k-1 (S i k-1 (s k )))
.

Back to the partial exchangeability of VRJP

Proof. Apply the previous proposition to VRJP, where f i,j (l j ) = W i,j (1+l j ) and

h i (l i ) = l 2 i + 2l i . The density d Y σ is 1 2 n exp   - s 0 j∼Yu W Yu,j S j (u) + 1 2 S Yu (u) + 1 du   n k=1 W i k-1 ,i k S i k (s k-1 ) + 1 S i k-1 (s k ) + 1 .
As our trajectory is left continuous without explosion, starting at i 0 , if we calculate the product through the trajectory, by telescopic simplification, it results that the product reduces to

i∈V i =in 1 S i (s) + 1 n k=1 W i k-1 ,i k .
To compute the integral inside the exponential, it is enough to note that, in the expression:

i∼j W i,j ( (S i (s) + 1)(S j (s) + 1) -1),
the local times S i (s), i ∈ V of the process Y only vary (linearly) with s when the process is at i, i.e., when Y t = i. Therefore, the derivative of the above expression with respect to s equals to

j∼Ys W Ys,j S j (s) + 1 2 S Ys (s) + 1
which is what we integrate inside the exponential. Whence ( 2) is proved, and expression (2) depends only on final local times and transition counts, the result hence follows.

Determination of time change h

In the sequel we work with the time changed process Y , to simplify notations, we will write d σ for d Y σ when it does not lead to any confusion. By Proposition 3, the density of certain trajectory contains an exponential term and a product term, let us denote

d σ = exp(-σ) • σ, with      σ = s 0 j∼Yv f Yv ,j (h -1 j (S j (v))) h ′ Yv (h -1 Yv (S Yv (v))) dv σ = n k=1 f i k-1 ,i k (h -1 i k (S Ys k (s k-1 ))) h ′ i k-1 (h -1 i k-1 (S Ys k-1 (s k )))
where the exponential term stems from those exponential waiting times, and the product term corresponds to the probability of the discrete chain. The heuristics of the proof in this subsection is the following: as we assumed partial exchangeability, if we consider two equivalent trajectories, then their densities share the same expression, by comparing them we can hence deduce certain equalities involving f i,j and h i etc. It turns out that these equalities determine h i s then f i,j s.

The following fact is simple but important, suppose that at time s, the random walker arrives at i 0 , each vertex i has accumulated local time l i := S i (s); then it jumps to i 1 after an amount of time t, by Proposition 3, the density has acquired a multiplicative factor exp

  - s+t s j∼i 0 f i 0 ,j • h -1 j (l j ) h ′ i 0 • h -1 i 0 (l i 0 + v) dv   • f i 0 ,i 1 • h -1 i 1 (l i 1 ) h ′ i 0 • h -1 i 0 (l i 0 + t) . ( 7
)
This fact is in constant use in the sequel, when we explicit the density of certain trajectory.

Lemma 2. Let σ = i 0 s 1 -→ i 1 s 2 -s 1 ----→ i 2 • • • i n-1 sn-s n-1 -----→ i n s-sn
---→ be a trajectory, then

σ = σ + σ where σ = s 0 j∈σ,j∼Yv f Yv,j (h -1 j (S j (v))) h ′ Yv (h -1 Yv (S Yv (v))) dv, σ = s 0 j / ∈σ,j∼Yv f Yv,j (h -1 j (S j (v))) h ′ Yv (h -1 Yv (S Yv (v)))
dv

and if τ is such that τ ∼ σ, then σ = τ .
Proof. Note that for j / ∈ σ, S j (u) = 0 for all u ≤ s. Let Ĥi be the primitive of

1 h ′ i • h -1 i such that Ĥi (0) = 0, σ = j / ∈σ s 0 1 Yv∼j f Yv,j (0) h ′ Yv (h -1 Yv (S Yv (v))) dv = j / ∈σ,i∈σ,j∼i f i,j (0)
s 0 1 Yv=i h ′ i (h -1 i (S i (v))) dv = j / ∈σ,i∈σ,j∼i f i,j (0) Ĥi (S i (s))
which depends only on final local times, thus if τ ∼ σ, then τ = σ.

In the sequel cst denotes some constant, which can vary from line to line.

Lemma 3. If the process X admits such a time change D which makes it partially exchangeable in density, then for any i ∼ j, there exists some constants λ i,j such that

f i,j (x) = λ i,j h ′ j (x), ∀x ≥ 0. (8)
Proof. Let ǫ > 0, consider the following two trajectories for the process Y :

σ = i ǫ - → j ǫ - → i t - → j s - → i • - → τ = i t - → j s - → i ǫ - → j ǫ - → i • - →
Note that σ and τ have the same transition counts and the final local times on vertex i, j are respectively equal. Thus the densities of these trajectories are a.s. equal by partial exchangeability. By Lemma 2,

d σ = σ • exp( σ + σ),
where

     σ = f i,j •h -1 j (0) h ′ i •h -1 i (ǫ) • f j,i •h -1 i (ǫ) h ′ j •h -1 j (ǫ) • f i,j •h -1 j (ǫ) h ′ i •h -1 i (ǫ+t) • f j,i •h -1 i (ǫ+t) h ′ j •h -1 j (ǫ+s) σ = ǫ 0 f i,j •h -1 j (0) h ′ i •h -1 i (v) dv + ǫ 0 f j,i •h -1 i (ǫ) h ′ j •h -1 j (v) dv + t 0 f i,j •h -1 j (ǫ) h ′ i •h -1 i (ǫ+v) dv + s 0 f j,i •h -1 i (ǫ+t) h ′ j •h -1 j (ǫ+v) dv. d τ = τ • exp( τ + τ ),
where

     τ = f i,j •h -1 j (0) h ′ i •h -1 i (t) • f j,i •h -1 i (t) h ′ j •h -1 j (s) • f i,j •h -1 j (s) h ′ i •h -1 i (t+ǫ) • f j,i •h -1 i (ǫ+t) h ′ j •h -1 j (ǫ+s) τ = t 0 f i,j •h -1 j (0) h ′ i •h -1 i (v) dv + s 0 f j,i •h -1 i (t) h ′ j •h -1 j (v) dv + ǫ 0 f i,j •h -1 j (s) h ′ i •h -1 i (t+v) dv + ǫ 0 f j,i •h -1 i (ǫ+t) h ′ j •h -1 j (s+v) dv;
We do not explicit σ and τ as they cancel when we compare these expressions (c.f. Lemma 2). Letting ǫ → 0 yields that exp( σ) = exp( τ ); therefore σ = τ , i.e.

∀s, t, f i,j • h -1 j (s) h ′ j • h -1 j (s) • f j,i • h -1 i (t) h ′ i • h -1 i (t) = cst.
Now fix t, let s vary, whence

∀s, f i,j • h -1 j (s) = cst • h ′ j • h -1 j (s),
and let λ i,j denotes this constant, as h -1 j is a diffeomorphism, its range is R + , which allows us to conclude.

The next lemma states in some sense that the exponential part and the product part appearing in the density of a trajectory can be treated separately. Lemma 4. Let σ, τ be two trajectories, and denote

d σ = exp( σ) • σ, d τ = exp( τ ) • τ, if σ ∼ τ , then σ = τ .
Proof. We have S Ys k (s k ) = S Ys k (s k-1 ), thus Lemma 3 yields that

f i k-1 ,i k •h -1 i k (S Ys k (s k-1 )) = λ i k-1 ,i k h ′ i k • h -1 i k (S Ys k (s k )). Whence the product part is σ = n k=1 f i k-1 ,i k (h -1 i k (S Ys k (s k-1 ))) h ′ i k-1 (h -1 i k-1 (S Ys k-1 (s k ))) = n k=1 λ i k-1 ,i k i =i 0 h ′ i • h -1 i (0) i =in h ′ i • h -1 i (S i (s))
, and the last term depends only on the transition counts and final local times.

Lemma 5. Let H i = h ′ i • h -1 i , then for some constant A i (recall that h i is assumed C 2 diffeomorphism), (H 2 i ) ′ = A i and if i ∼ j, then λ i,j A j = λ j,i A i .
Remarks 6. The latest equality tells that the process is reversible. However, we did not assume the reversibility of the process, but vertex reinforced jump processes are reversible (as a mixture of reversible Markov jump process), so are the edge reinforced random walks. In contrast, directed edge reinforced random walks are mixtures of non reversible Markov chains, with independent Dirichlet environments. We can hence expect that the reversibility is a consequence of a non oriented linear reinforcement (where linearity leads to partial exchangeability). Proof. Recall that we have assumed that the graph is strongly connected, i.e. if i, j are two adjacent vertices, there exists a shortest cycle

i n-1 i n i 1 i 2 i n-2 i n-1 i n i 1 i 2 i n-2 σ τ
i 1 ∼ i 2 ∼ i 3 • • • ∼ i n ∼ i 1 with i 1 = i, i n = j
and the i k s are distinct and n ≥ 2. Let (i 1 = i, i 2 , i 3 , . . . , i n = j) be a cycle as described, consider the trajectories (c.f. Figure 3)

σ = i 1 r 1 -→ i n r 2 -→ i 1 s 1 -→ i 2 s 2 -→ i 3 • • • i n-2 s n-2 ---→ i n-1 s n-1 ---→ i n τ = i 1 r 1 -→ i 2 s 2 -→ i 3 • • • i n-2 s n-2 ---→ i n-1 s n-1 ---→ i n r 2 -→ i 1 s 1 -→ i n .
As σ ∼ τ , by Lemma 4 and Lemma 2, σ = τ . Also let

σ ′ = i 1 r 1 -→ i n r 2 -→ i 1 s 1 -→ i 2 s 2 -→ i 1 τ ′ = i 1 r 1 -→ i 2 s 2 -→ i 1 s 1 -→ i n r 2 -→ i 1 ,
thus σ′ = τ ′ . We are going to compute explicitly σ, τ etc, using [START_REF] Da Freedman | De finetti's theorem in continuous time[END_REF], let s =

r 1 +r 2 +s 1 +• • •+s n-1 and recall that Ĥi is the primitive of 1 h ′ i • h -1 i such that Ĥi (0) = 0. σ = (i,j)∈σ 2 ,i∼j λ i,j s 0 1 Yv=i h ′ j • h -1 j (S j (v)) h ′ i • h -1 i (S i (v)) dv = λ i 1 ,i 2 H i 2 (0)( Ĥi 1 (r 1 + s 1 ) -Ĥi 1 (0)) + λ i 2 ,i 1 H i 1 (r 1 + s 1 )( Ĥi 2 (s 2 ) -Ĥi 2 (0)) + λ i 1 ,in H in (0)( Ĥi 1 (r 1 ) -Ĥi 1 (0)) + H in (r 2 )( Ĥi 1 (r 1 + s 1 ) -Ĥi 1 (r 1 )) + λ in,i 1 H i 1 (r 1 )( Ĥin (r 2 ) -Ĥin (0)) + λ in,i n-1 H i n-1 (0)( Ĥin (r 2 ) -Ĥin (0)) + λ i n-1 ,in H in (r 2 )( Ĥi n-1 (s n-1 ) -Ĥi n-1 (0)) + ∆
where ∆ is defined as follows: let

Q k := H i k (0)( Ĥi k-1 (s i k-1 ) -Ĥi k-1 (0)) and Q ′ k := H i k (s k )( Ĥi k+1 (s i k+1 ) -Ĥi k+1 (0)), ∆ = n-1 k=3 λ i k-1 ,i k Q k + λ i k ,i k-1 Q ′ k-1 .
For τ we have: τ = (i,j)∈τ 2 ,i∼j λ i,j s 0

1 Yv=i h ′ j • h -1 j (S j (v)) h ′ i • h -1 i (S i (v)) dv = λ i 1 ,i 2 H i 2 (
0)( Ĥi 1 (r 1 ) -Ĥi 1 (0)) + H i 2 (s 2 )( Ĥi 1 (r 1 + s 1 ) -Ĥi 1 (r 1 )) + λ i 2 ,i 1 H i 1 (r 1 )( Ĥi 2 (s 2 ) -Ĥi 2 (0)) + λ i 1 ,in H in (0)( Ĥi 1 (r 1 ) -Ĥi 1 (0)) + H in (r 2 )( Ĥi 1 (r 1 + s 1 ) -Ĥi 1 (r 1 ))

+ λ in,i 1 H i 1 (r 1 )( Ĥin (r 2 ) -Ĥin (0)) + λ in,i n-1 H i n-1 (s n-1 )( Ĥin (r 2 ) -Ĥin (0)) + λ i n-1 ,in H in (0)( Ĥi n-1 (s n-1 ) -Ĥi n-1 (0)) + ∆ with the same ∆. Also σ′ = λ i 1 ,i 2 H i 2 (0)( Ĥi 1 (r 1 ) -Ĥi 1 (0)) + H i 2 (0)( Ĥi 1 (r 1 + s 1 ) -Ĥi 1 (r 1 )) + λ i 2 ,i 1 H i 1 (r 1 + s 1 )( Ĥi 2 (s 2 ) -Ĥi 2 (0)) + λ i 1 ,in H in (0)( Ĥi 1 (r 1 ) -Ĥi 1 (0)) + H in (r 2 )( Ĥi 1 (r 1 + s 1 ) -Ĥi 1 (r 1 ))

+ λ in,i 1 H i 1 (r 1 )( Ĥin (r 2 ) -Ĥin (0)) τ ′ = λ i 1 ,i 2 H i 2 (0)( Ĥi 1 (r 1 ) -Ĥi 1 (0)) + H i 2 (s 2 )( Ĥi 1 (r 1 + s 1 ) -Ĥi 1 (r 1 ))

+ λ i 2 ,i 1 H i 1 (r 1 )( Ĥi 2 (s 2 ) -Ĥi 2 (0))
+ λ i 1 ,in H in (0)( Ĥi 1 (r 1 ) -Ĥi 1 (0)) + H in (0)( Ĥi 1 (r 1 + s 1 ) -Ĥi 1 (r 1 ))

+ λ in,i 1 H i 1 (r 1 + s 1 )( Ĥin (r 2 ) -Ĥin (0)).

Recall that σ -σ ′ = τ -τ ′ , which leads to λ in,i n-1 H i n-1 (0)( Hin (r 2 ) -Hin (0)) + λ i n-1 ,in H in (r 2 )( Ĥi n-1 (s n-1 ) -Ĥi n-1 (0)) = λ i 1 ,in (H in (r 2 ) -H in (0))( Ĥi 1 (r 1 + s 1 ) -Ĥi 1 (r 1 ))

+ λ in,i 1 (H i 1 (r 1 ) -H i 1 (r 1 + s 1 ))( Ĥin (r 2 ) -Ĥin (0)) + λ in,i n-1 H i n-1 (s n-1 )( Hin (r 2 ) -Hin (0)) + λ i n-1 ,in H in (0)( Ĥi n-1 (s n-1 ) -Ĥi n-1 (0)) letting s n-1 → 0 leads to λ i 1 ,in (H in (r 2 ) -H in (0))( Ĥi 1 (r 1 + s 1 ) -Ĥi 1 (r 1 )) = λ in,i 1 (H i 1 (r 1 + s 1 ) -H i 1 (r 1 ))( Ĥin (r 2 ) -Ĥin (0)) as i 1 , i n , r 2 , s 1 , r 1 are arbitrary, divide the formula by r 2 s 1 and let r 2 , s 1 go to zero leads to λ i 1 ,in H ′ in (0) Ĥ′ i 1 (r 1 ) = λ in,i 1 H ′ i 1 (r 1 ) Ĥ′ in (0), finally note that Ĥ′ i = 1/H i , thus λ i 1 ,in (H 2 in ) ′ (0) = λ in,i 1 (H 2 i 1 ) ′ (r 1 ). Lemma 6. For all i ∼ j, let W i,j = λ i,j A j /2 = λ j,i A i /2, there exists constant ϕ j depends only on j, such that f i,j (x) = W i,j x + ϕ j .

Proof. As (H 2 j (s)) ′ = A j , there exists B j such that H 2 j (s) = A j s + B j , therefore f i,j • h -1 j (s) = λ i,j H j (s) = λ i,j A j s + B j .

On the other hand, (h

-1 j ) ′ (s) = 1 √ A j s+B j
, thus for some C j ,

h -1 j (s) = 2 A j A j s + B j + C j .
f i,j (h -1 j (s)) = f i,j ( 2 A j A j s + B j + C j ) = λ i,j A j s + B j , which leads to f i,j (x) = W i,j x + ϕ j , where ϕ j is some constant depends only on j. Applying the time change

D(s) = i l i (s) -ϕ i ϕ i ,
the resulting process will be of jump rate W i,j ϕ i ϕ j (1 + T j (t))

where T j (t) is the local time for the time changed process Z t = X D-1 (t) .

w 0

 0 (e) = a e for all e ∈ E , w n+1 (e) = w n (e) + 1 for e = {Z n , Z n+1 } ∈ E, w n (e) for e ∈ E \ {{Z n , Z n+1 }}. Let P (a) v 0 denote the probability of the ERRW on G starting at v 0 with initial weights a = (a e ) e∈E . Then P (a)v 0 is defined byZ 0 = v 0 , P (a)v 0 -a.s.,P (a) v 0 (Z n+1 = v|Z 0 , . . . , Z n ) =wn({Zn,v})e,Zn∈e wn(e)

Figure 2 :

 2 Figure 2: An example of σ and T (σ).

Figure 3 :

 3 Figure 3: the trajectories σ and τ in Lemma 5.

  Definition 8. A continuous time nearest neighbor jump process is said to be partially exchangeable in density if the densities are equal for any two equivalent trajectories. Or equivalently, the density depends only on final local times and the transition counts.

t 0 1 σs=i ds denotes the local time.

3.2 Equivalence of the two notions

It turns out that in the case of nearest neighbor jump process with continuous jump rate functions, the notion of partial exchangeability in Definition 5 and in Definition 8 are equivalent.

Proposition 1. If a continuous time nearest neighbor jump process is partially exchangeable in the sense of Definition 8, then it is partially exchangeable in the sense of Definition 5.

Acknowledgments: I would like to thank Christophe Sabot for his constant support in this project. I would like to thank an anonymous referee for carefully reading the paper and providing corrections.


Michel Paques

Florence Rossant

Nicolas Lermé

Chahira Miloudi

Caroline Kulcsár

José-Alain Sahel

Kevin Loquin

Laurent Mugnier

Isabelle Bloch

Edouard Koch

Caroline Kulcsar

Adaptive optics imaging of retinal microstructures: image processing for medical applications (Orale)

Keywords: adaptive optics retinal imaging, agerelated macular degeneration, photoreceptors, arterial hypertension

teaching and research institutions in France or abroad, or from public or private research centers.

INTRODUCTION

Adaptive optics (AO) is an optoelectronic technique that compensates for optical aberrations by the mean of a deformable mirror. AO imaging was first developed in astronomy and later applied to fundus images. AO improves the lateral resolution of fundus images, which can reach the micrometric scale. AO-based fundus cameras allow in vivo visualization of microstructures of the retina, making the retina the only human tissue in which subcellular structures may be routinely observed in vivo. High resolution imaging of photoreceptors therefore holds promises for cellular-scale diagnosis and guidance of therapy. An increasing number of clinical applications are identified for AO imaging, hence it is progressively translating from the lab to the clinics and will undoubtedly become a routine procedure for many retinal diseases in the next years.

The main applications of AO in ophthalmology span from photoreceptor to retinal pigment epithelial cells and vessels, each of them being affected by specific diseases. The first photoreceptor images using AO were obtained more than 15 years ago [START_REF] Liang | Supernormal vision and high-resolution retinal imaging through adaptive optics[END_REF]. However, there are still some obstacles that prevent the widespread use of AO in clinical routine. One relates to the integration of information from AO into medical decisions. Although AO fundus imaging technology and image processing is reaching technological maturity, clinical interpretation of the complex features of diseased retina indeed remains challenging because of a number of factors. Among these factors are the level of pigmentation of the fundus, the transparency of the retina, the presence of other sources of light reflection in diseased retina, the spatial and temporal variability of photoreceptor reflectance, and the variable orientation of photoreceptor outer segments. Other issues are related to the fundus image technology: scanning laser and flood illumination show some differences in term of resulting imaging which are not yet fully understood. A crucial factor is also that, in order to optimize the management of these complex images by ophthalmologists, carefully designed structure-specific softwares for image processing and analysis are necessary .

We will show here that there are some clinical applications of AO in public health that already offers promising opportunities for clinicians. Most of this paper relates to flood illumination-based AO.

MEDICAL APPLICATIONS OF ADAPTIVE OPTICS IMAGING

Photoreceptors

Historically, photoreceptors were the first microstructures specifically identified by AO imaging. The diameter of the largest cones reaches ~3µm, making them visible under the aspect of a mosaic of hyperreflective dots (Figure 1). Rods and the smaller foveal cones became more recently accessible with the increasing resolution of AO imaging systems. This allowed documenting the loss of photoreceptors in diseases such as retinitis pigmentosa. A precise counting of residual photoreceptors may be an essential biomarker of the effect of biotherapies such as gene or cell therapy. Image processing has brought major improvements to AO imaging. Cone counting algorithms have now reached a Supported by Agence Nationale de la Recherche (ANR-09-TECS-009 and ANR-12-TECS-0015-03), the Foundation Fighting Blindness and the french Ministry of Research (CIFRE131145A10). high degree of reliability. One of the most reliable and accurate approach is based on a recursive construction of thresholded components, when the seeds of the recursions are the regional maxima of the input image [START_REF] Loquin | Photoreceptor detection in in-vivo adpative optics images of the retina: towards a simple interactive tool for the physician[END_REF][START_REF] Loquin | Automatic Photoreceptor Detection in In-Vivo Adaptive Optics Retinal Images: Statistical Validation[END_REF].

These automatic analysis algorithms benefits from advanced preprocessing methods. As flood-illumination systems usually produce noisy images making microscopic structures hardly visible, a commonly used solution is to register these images and average them to increase the signal-to-noise ratio. Myopic deconvolution of AO images has improved the quality of flood-illumination AO images [START_REF] Blanco | Marginal blind deconvolution of adaptive optics retinal images[END_REF] since it enables to extract cone mosaics from raw (i.e. non-averaged) images (Figure 2), hence avoiding artifacts related to photoreceptor scintillation and/or subsequent image processing. Sub-pixellic registration, developed by one of [START_REF] Kulcsár | Robust processing of images sequences produced by an adaptive optics retinal camera[END_REF] also had a significant impact on the resulting image quality, through the suppression of peripheral rotational artifacts (Figure 3) and hence enlarging the field over which photoreceptors can be counted. Cone photoreceptors have a strong directional absorbance and reflectance selectivity, a property referred to as the Stiles-Crawford effect (SCE) [START_REF] Stiles | The luminous efficiency of rays entering the eye pupil at different points[END_REF][START_REF] Gao | Measuring directionality of the retinal reflection with a Shack-Hartmann wavefront sensor[END_REF][START_REF] Westheimer | Directional sensitivity of the retina: 75 years of Stiles-Crawford effect[END_REF]. This may be observed on AO image by the strong variations of the cone mosaic under different illumination incidence. Preliminary work from our group suggests that the analysis of the SCE, which requires the fusion of images taken at various illumination angles, may provide a novel biomarker related to the healing process of photoreceptors.

Arterial hypertension

As an easily accessible part of the microcirculation, the retina provides a convenient mean for in vivo evaluation of the consequences of arterial hypertension. Several epidemiological studies have shown that the severity of hypertensive retinopathy as observed on fundus photographs is indeed correlated with the incidence of arterial hypertension or of end-organ damage such as coronary diseases, age-related cognitive decline and lacunar stroke. By AO imaging, the blood column of vessels appears as dark stripes 50-100 µm wide with an axial light reflection (Figure 3). A linear structure along both sides of the blood column of arteries is visible, which most likely corresponds to parietal structures [START_REF] Chui | The use of forward scatter to improve retinal vascular imaging with an adaptive optics scanning laser ophthalmoscope[END_REF].

Supported by Agence Nationale de la Recherche (ANR-09-TECS-009 and ANR-12-TECS-0015-03), the Foundation Fighting Blindness and the french Ministry of Research (CIFRE131145A10). Advanced image processing methods were proposed by our group, allowing automatical or semi-automatical segmentation of parietal structures of arterioles [START_REF] Lermé | Coupled Parallel Snakes For Segmenting Healthy and Pathological Retinal Arteries in Adaptive Optics Images[END_REF][START_REF] Lermé | Segmentation of Retinal Arteries in Adaptive Optics Images[END_REF]. Our approach relies on three main steps: (i) detection and accurate localization of the vessel axial reflection, (ii) presegmentation through a tracking procedure relying mostly on gradient information, (iii) accurate delineation of the arterial walls, based on a new active contour model integrating parallelism information [START_REF] Lermé | Coupled Parallel Snakes For Segmenting Healthy and Pathological Retinal Arteries in Adaptive Optics Images[END_REF]. The proposed energy functional, that defines the active contour model, makes evolve independently four curves towards the artery wall contours assuming that they are almost parallel to the axial reflection. In [START_REF] Lermé | Segmentation of Retinal Arteries in Adaptive Optics Images[END_REF], other structural features were also modelled and integrated, in order to improve the robustness against low contrasted walls and morphological deformations that occur along vessels in case of pathologies. For this, new coupled energy terms were added in the energy functional, modeling symetry properties on either side of the axial reflection. In this new model, the four curves do no evolve independently any more but jointly. The strength of parallelism and symetry constraints is controlled through weighting parameters that have been optimized experimentally so as to cover the largest variety of images and pathologies. However, in case of failure, the clinician can adjust them manually (semi-automatic mode). From the segmentation result, graphic representations of morphometric parameters of a given vessel segment can be generated, which facilitates medical analysis of vascular morphometry. Figure 4 shows an example of segmentation (a), along with several measurements: diameters and wall thckness (b) and wall-to-lumen ratio (c). The latest morphometric biomarker (wall-to-lumen ratio, WLR) represents the ratio between the wall thickness and the vessel internal diameter (lumen). It is calculated all along the artery segment, from measurements averaged on a sliding window (50μm wide in this case). In a cohort of healthy and hypertensive subjects, using these tools we found that the thickness of the wall of arterioles was positively correlated with systolic and diastolic blood pressure. In multivariate analysis taking into account age, lumen diameter, and systolic and diastolic blood pressure, only diastolic blood pressure and lumen diameter were found independently correlated to the WLR [START_REF] Koch | Supported by Agence Nationale de la Recherche (ANR-09-TECS-009 and ANR-12-TECS-0015-03), the Foundation Fighting Blindness and the french Ministry of Research (CIFRE131145A10)[END_REF]. We could also observe and measure changes of vascular morphometry in focal lesions, which provided insights into their pathophysiology.

AO imaging is currently the only technique allowing a direct view of vascular parietal structures. AO imaging of retinal arterioles therefore offers a unique opportunity to explore microvascular changes related to age, diabetes or arterial hypertension in vivo in humans at a near-histology level. Longitudinal explorations of hypertensive patients are underway, addressing questions such as the predictive value of such measures on the incidence of end-organ damage, and the effect of treatment on microvascular parameters. Phenotyping of retinal vessels by AO imaging may contribute to a better management and understanding of end-organ damage, especially in the brain given the functional and anatomical similarities of the retinal and cerebral circulations. Stratification of end-organ damage risk may be improved by biomarkers issued from AO imaging.

Age-related macular degeneration

Geographic atrophy (GA), the atrophic form of age-related macular degeneration, is one of the most frequent causes of vision loss in industrialized countries. The disease is characterized by the primary alteration of a single layer of cells, the retinal pigment epithelium (RPE). Accurate characterization of the progression of atrophic lesions Supported by Agence Nationale de la Recherche (ANR-09-TECS-009 and ANR-12-TECS-0015-03), the Foundation Fighting Blindness and the french Ministry of Research (CIFRE131145A10).

during GA is important for the estimation of its long-term prognosis and consequently for the evaluation of therapeutic results. It is likely that improved imaging of small GA lesions may allow evaluating the effect of treatment at earlier stages of the disease, where such treatment may be more susceptible to show efficacy.

AO imaging improves the resolution of the changes affecting the RPE [START_REF] Gocho | Adaptive Optics Imaging of Atrophic Geography[END_REF] (figure 5). Time-lapse imaging revealed for the first time the migration of these cells, presumably inflammatory cells. Such dynamic changes can be observed within as well as outside atrophic areas. Time-lapse AO imaging also allowed detection of very small spots of atrophy and follow-up of their progression. AO imaging is therefore of interest to detect the emergence of atrophic spots, to document the progression of GA at a small temporal and spatial scale, and also to investigate the dynamic process underlying the redistribution of (presumed) inflammatory cells. A complex, dynamic process of redistribution of melanin-loaded cells indeed precedes and accompanies the emergence and progression of GA. Further studies correlating AO to conventional en face imaging and OCT would help to better define the place of AO imaging in the management of GA patients.

Vascular inflammation

Recently, we observed that AO may reveal paravascular cellular infiltrates in eyes with vascular inflammation [START_REF] Errera | Retinal vasculitis imaging by adaptive optics[END_REF], which are not otherwise detectable. This offers the opportunity of diagnosing vascular inflammation and also monitoring the therapy of such patients in a noninvasive way. Since vascular inflammation occurs in many general diseases, including multiple sclerosis, detection of paravascular infiltrates may of interest in a number of neurological and/or inflammatory diseases.

CONCLUSION

Up to now, the main application of AO fundus imaging has been photoreceptor detection and counting; this somewhat overshadowed the exploration of other retinal structures.

With a simple procedure applicable in a routine setting, AO imaging can document microscopic features of common ophthalmological and general diseases and will possibly become a reference tool for their monitoring. Such medical applications will rely heavily on image processing and analysis, for which specific applications are needed for each structure of interest, and possibly of each disease.

Figure 1 .

 1 Figure 1. Imaging of the normal cone mosaic by adaptive optics (AO) flood illumination imaging (bar, 300µm).

Figure 2 .

 2 Figure 2. increased contrast of the cone mosaic on raw image resulting from myopic deconvolution, which enables cone counting to be performed.

Figure 3 .

 3 Figure 3. Improvement of AO images by subpixellic registration (bar, 300µm).

Figure 3 .

 3 Figure 3. Imaging of the wall of an arteriole (between arrows in the magnification) by AO flood illumination imaging (bar, 300µm).

Figure 4 .

 4 Figure 4. Segmentation of a retinal artery segment in AO imaging (a) and extracted measurements (µm): internal diameter (lumen, in red), external diameter (magenta) and wall thickness (brown) (b); WLR calculated from the previous measures (c).

Figure 5 .

 5 Figure 5. Illustration of pigment redistribution in a case of geographic atrophy. Note the multiple black dots in the AO image, presumably representing inflammatory cells.

Supported by Agence Nationale de la Recherche (ANR-09-TECS-009 and ANR-12-TECS-0015-03), the Foundation Fighting Blindness and the french Ministry of Research (CIFRE131145A10).