Marius Bozga
email: bozga@imag.fr

Codruta Girlea
email: girlea@imag.fr

Radu Iosif
email: iosif@imag.fr

Codrut ¸a Gîrlea

Iterating Octagons

published or not. The documents may come   L'archive ouverte pluridisciplinaire

Introduction

Counter automata (register machines) are widely investigated models of computation. Since the result of Minsky [START_REF] Minsky | Computation: Finite and Infinite Machines[END_REF] showing Turing-completeness of 2-counter machines, research on counter automata pursued in two directions. The first one is defining subclasses of counter automata for which various decision problems (e.g. reachability, emptiness, boundedness, disjointness, containment, equivalence) are found to be decidable. Examples include Reversal-bounded Counter Machines [START_REF] Ibarra | Reversal-bounded multicounter machines and their decision problems[END_REF], Petri Nets and Vector Addition Systems [START_REF] Reutenauer | Aspects Mathématiques des Réseaux de Petri[END_REF] or Flat Counter Automata [START_REF] Leroux | Flat counter automata almost everywhere! In ATVA[END_REF]. Often, decidability of various problems is achieved by defining the set of reachable configurations in a decidable logic, such as Presburger arithmetic [START_REF] Presburger | Über die Vollstandigkeit eines gewissen Systems der Arithmetik[END_REF]. Such definitions are precise, i.e. no information is lost by the use of over-approximation.

Another, orthogonal, direction of work is concerned with finding sound (but not necessarily complete) answers to the decision problems mentioned above, in a cost-effective way. Such approaches use abstract domains (Polyhedra [START_REF] Cousot | Automatic discovery of linear restraints among variables of a program[END_REF], Octagons [START_REF] Miné | The octagon abstract domain[END_REF], Difference Constraints [START_REF] Bagnara | Data-Flow Analysis for Constraint Logic-Based Languages[END_REF], etc.) and compute fixed points that are over-approximations of the set of reachable configurations.

Both approaches to the analysis of counter automata benefit, in some sense, from algorithms for computing transitive closures of arithmetic relations. If I is the initial set of configurations, and R is the transition relation of the counter automaton, then R * (I) (the image of I through R * ) is the set of reachable configurations, where R * = n≥0 R n is the transitive closure of R. The problem lies essentially in expressing the infinite disjunction from the definition of R * , in a finite way.

In this paper we consider octagonal transition relations that are conjunctions of atoms of the form ±x ± y ≤ c, where x and y are counter values, either at the current step, or at the next step (in which case they are denoted by primed variables), and c ∈ Z is an integer constant. For this class of relations, we prove that the transitive closure is expressible in Presburger arithmetic [START_REF] Presburger | Über die Vollstandigkeit eines gewissen Systems der Arithmetik[END_REF]. This improves the previous result of Comon and Jurski [START_REF] Comon | Multiple counters automata, safety analysis and presburger arithmetic[END_REF], showing that the transitive closure of a difference bound constraint (a conjunction of atoms of the form xy ≤ c, with x, y possible primed counters, and c ∈ Z) is Presburgerdefinable.

We adopt the classical representation of octagonal constraints (or octagons) ϕ(x 1 , . . . , x n ) as difference bound constraints ϕ(y 1 , . . . , y 2n ), where y 2i-1 stands for +x i and y 2i for -x i , and the implicit condition y 2i-1 + y 2i = 0, 1 ≤ i ≤ n. With this convention, [START_REF] Bagnara | An improved tight closure algorithm for integer octagonal constraints[END_REF] provides an algorithm for computing the canonical form of an octagon, by first computing the strongest closure of the corresponding difference bound constraint (using the classical Floyd-Warshall cubic algorithm), and subsequently tightening the constraints of the form y 2iy 2i-1 ≤ c and y 2i-1y 2i ≤ d by adjusting c to 2⌊c/2⌋ and d to 2⌊d/2⌋, respectively1 . We apply this idea to tightening octagonal relations of the form R k , where R is an octagon and k ≥ 1 is the arbitrary number of iterations, obtaining in this way a Presburger formula ψ(k, x 1 , ...,

x n , x ′ 1 , ..., x ′ n ) equivalent to the k-th iteration of ϕ. The transitive closure of R is thus R * = ∃k.ψ.
The main application of this result is that the following problems are decidable, for the class of counter automata with octagonal transition rules:

-reachability: the automaton has a run leading from a configuration in I to a configuration in S, where I and S are Presburger-definable sets of configurations.

-emptiness: the automaton has at least a run starting with all counters set to zero, and leading to a final control location.

In particular, decidability of the reachability problem is useful for the verification of safety properties (e.g. assertion checking) of integer programs, whereas the emptiness problem is a promising approach for the analysis of programs with integer arrays. Indeed, the works of [START_REF] Habermehl | What else is decidable about integer arrays? In FOSSACS[END_REF][START_REF] Bozga | A logic of singly indexed arrays[END_REF] reduce the satisfiability problem for two logics on integer arrays to the emptiness problem of a counter automaton. By enlarging the class of counter automata for which this problem is decidable, we enhance the expressiveness of (decidable) logics of integer arrays. Finally, we have implemented our method in a tool for the analysis of counter automata, called FLATA [10]. In particular, FLATA computes the transitive closure of elementary octagonal cycles, which is used in computing the set of reachable configurations for flat counter automata. We have experimented with a number of test cases and provided several experimental results.

Related Work

The domain of counter machines has been investigated starting with the seminal work of Minsky [START_REF] Minsky | Computation: Finite and Infinite Machines[END_REF]. His result on Turing-completeness of 2-counter machines motivated research on subclasses of counter automata, for which some problems are decidable, such as: Reversal-bounded Counter Machines [START_REF] Ibarra | Reversal-bounded multicounter machines and their decision problems[END_REF], Petri Nets and Vector Addition Systems [START_REF] Reutenauer | Aspects Mathématiques des Réseaux de Petri[END_REF] or Flat Counter Automata [START_REF] Leroux | Flat counter automata almost everywhere! In ATVA[END_REF].

The class of Flat Counter Automata which is closest to our work is the one studied by Comon and Jurski in [START_REF] Comon | Multiple counters automata, safety analysis and presburger arithmetic[END_REF]. In this work, they prove that the transitive closure of a difference bound constraint is Presburger-definable. In [START_REF] Bozga | Flat parametric counter automata[END_REF] we showing how to effectively compute the transitive closure of a difference bound constraint directly as a finite set of linear inequation systems, opening thus the possibility of using SMT tools for the analysis of such models.

Octagonal constraints are investigated in the comprehensive paper of Miné [START_REF] Miné | The octagon abstract domain[END_REF]. Among other results, his paper presents a cubic time tightening algorithm for an octagonal constraint, which is an improvement of the classical algorithm of Harvey and Stuckey [START_REF] Harvey | A unit two variable per inequality integer constraint solver for constraint logic programming[END_REF]. However, this algorithm is not suitable to tightening octagonal constraints of parametric size, as the ones we obtain by iteration. For this reason, we adapt the (newer) tightening algorithm of Bagnara, Hill and Zaffanella [START_REF] Bagnara | An improved tight closure algorithm for integer octagonal constraints[END_REF] to octagons obtained by arbitrary iteration, as they prove that it is sufficient to adjust the constants after the computation of minimal paths (by the Floyd-Warshall algorithm).

Our paper extends to octagons the class of non-deterministic relations for which the transitive closure can be effectively expressed in Presburger arithmetic. This result relies essentially on our method from [START_REF] Bozga | Flat parametric counter automata[END_REF] to computing minimal paths in constraint graphs obtained by iteration, and the idea of [START_REF] Bagnara | An improved tight closure algorithm for integer octagonal constraints[END_REF] for tightening (finite) octagonal relations.

On a different line of work, Boigelot [START_REF] Boigelot | On iterating linear transformations over recognizable sets of integers[END_REF], and Finkel and Leroux [START_REF] Finkel | How to compose presburger-accelerations: Applications to broadcast protocols[END_REF] have studied the computation of transitive closure for affine relations of the form x ′ = Ax + b. Their class of systems differs in nature from ours, in what their transition relations are functional (deterministic), whereas octagons are not. Because of this, their results seem incomparable to ours.

Roadmap The paper is organized as follows. Section 2 introduces difference bound constraints and difference bound relations, recalling a number of results on DBMs, Section 3 defines octagonal constraints, and gives a necessary quantifier elimination result on octagons, while Section 4 is dedicated to our main result, namely computing the transitive closure of an octagonal relation. Section 5 gives implementation results, and finally, Section 6 concludes. Due to space reasons, all proofs are given in [START_REF] Bozga | Iterating octagons[END_REF].

Difference Bounds

In this section we recall several definitions and results on difference bound constraints, and the transitive closure of difference bound relations. In the rest of the paper, Z denotes the set of integer numbers, n > 0 is the number of variables and x = (x 1 , x 2 , . . . , x n ) is the tuple of variables. If ϕ(x) is a logical formula in which x 1 , x 2 , . . . , x n occur free, and v = (v 1 , v 2 , . . . , v n ) is a tuple of values, ϕ[v/x] denotes the formula in which each occurrence of x i has been replaced by v i . We denote by v |= ϕ the fact that ϕ[v/x] is logically equivalent to true. We say that ϕ is consistent if there exists at least one v ∈ Z n such that v |= ϕ, and inconsistent otherwise. If ϕ is a formula, let AP (ϕ) denote the set of atomic propositions in ϕ. If m is a matrix, then m i,j denote the element of m situated at line i and column j.

Difference Bound Constraints

Definition 1. Let x = (x 1 , x 2 , ..., x n ) be a set of variables ranging over Z. Then a formula φ(x) is a difference bound constraint if it is equivalent to a finite conjunction of atomic propositions of the form

x i -x j ≤ α i,j , i = j, 1 ≤ i, j ≤ n, where α i,j ∈ Z.
For instance, xy = 5 is a difference bound constraint, as it is equivalent to xy ≤ 5 ∧ yx ≤ -5. In practice difference bound constraints are represented either as matrices or as graphs, each of these representations being suitable for particular procedures (e.g. closure, iteration). We define these representations and procedures below. Definition 2. Let x = (x 1 , x 2 , ..., x n ) be a set of variables ranging over Z and φ(x) be a difference bound constraint. Then a difference bound matrix (DBM) representing φ is an n × n matrix m such that:

m i,j = α i,j if (x i -x j ≤ α i,j ) ∈ AP (φ) ∞ otherwise For a n × n DBM m we denote by γ(m) = {v ∈ Z n | v i -v j ≤ m i,j , 1 ≤ i, j ≤ n} the set of concretizations of m.
Notice that this is exactly the set of models of the corresponding difference bound constraint ϕ, i.e. γ(m) = {v | v |= ϕ}. A DBM m is said to be consistent if γ(m) = ∅, and inconsistent otherwise. Definition 3. Let x = (x 1 , x 2 , ..., x n ) be a set of variables ranging over Z and φ(x) be a difference bound constraint. Then φ can be represented as a weighted graph G with vertices x 1 , x 2 , ..., x n in which there is an arc with weight α i,j between x i and x j in G if there is a constraint x ix j ≤ α i,j in φ. This graph is also called a constraint graph.

Whenever the graph G is obvious from the context, we will denote by x α -→ y the fact that there exists an edge with weight α from x to y in G. Notice that the DBM m of a difference bound constraint is the incidence matrix of its corresponding constraint graph. This graph will be denoted as G(m) in the following. The three notions (constraint, matrix, graph) are related by the following property : Property 1. Let φ be a difference bound constraint, m the corresponding DBM and G(m) the corresponding graph. Then the following three statements are equivalent :

1. φ is inconsistent 2. G(m) contains at least one negative weight cycle 3. γ(m) = ∅
On one hand, the DBM representation of a difference bound constraint is suitable for computing its normal form, i.e. the most "explicit" formula that has the same set of models. Definition 4. An n × n consistent DBM m is said to be closed if and only if the following hold:

1. m i,i = 0 ∀1 ≤ i ≤ n 2. m i,j ≤ m i,k + m k,j ∀1 ≤ i, j, k ≤ n The shortest path closure of a consistent DBM m is a closed DBM m * such that γ(m) = γ(m * ). It is well-known that, if m is consistent then m * is
unique, and it can be computed from m in time O(n 3 ), by the classical Floyd-Warshall algorithm. Moreover, the following condition holds, for all 1 ≤ i, j ≤ n :

m * i,j = min k-1 l=0 m i l ,i l+1 |x i = x i0 x i1 ...x i k = x j is a path in G(m)
On the other hand, the graph representation of a difference bound constraint is suitable for existential quantifier elimination2 . Concretely, given a difference bound constraint φ(x), the formula ∃x k .φ is a difference bound constraint as well, and its corresponding graph is effectively computable from the graph of φ, as shown by the following property : Property 2. Let x = (x 1 , x 2 , ..., x n ) be a set of variables ranging over Z, φ(x) be a difference bound constraint and m be the corresponding n × n DBM. Then, for any 1 ≤ k ≤ n the formula ∃x k .φ(x) is a difference bound constraint, and its corresponding constraint graph is obtained by erasing the vertex x k together with the incident arcs from the graph G(m * ). Moreover, the DBM of the resulting graph is also closed.

The importance of this result will be made clear in the next section, because it directly implies that the class of difference bound relations is closed under composition. This is the first ingredient of our transitive closure method.

Difference Bound Relations

Definition 5. Let x = (x 1 , x 2 , ..., x n ), x ′ = (x ′ 1 , x ′ 2 , ..., x ′ n ) be sets of variables ranging over Z. A relation R(x, x ′ ) is a difference bound relation if it is equiv- alent to a finite conjunction of terms of the form x i -x j ≤ a i,j ,x ′ i -x j ≤ b i,j , x i -x ′ j ≤ c i,j , or x ′ i -x ′ j ≤ d i,j
, where a i,j , b i,j , c i,j , d i,j ∈ Z, 1 ≤ i, j ≤ n. According to the previous section, a difference bound relation can be represented as a constraint graph with nodes in the set x ∪ x ′ , such that x α -→ y if and only if xy ≤ α ∈ AP (R), for all x, y ∈ x ∪ x ′ . Just like before, the incidence matrix of this graph is the DBM of the relation.

In the following, we denote x (l) = (x

(l) 1 , x (l) 2 , . . . , x (l)
n ), x (≥l) = s≥l x (s) and x (≤l) = s≤l x (s) , for any l ≥ 0. Given a difference bound relation R(x, x ′ ), we define the k-th iteration of R:

R k (x (0) , x (k) ) = ∃x (1) . . . ∃x (k-1) .R(x (0) , x (1) ) ∧ . . . ∧ R(x (k-1) , x (k) ) for all k ≥ 1. The transitive closure of R is defined as R * = ∃k.R k .
Notice that, since existential quantifiers can be eliminated from difference bound constraints (cf. Theorem 2), difference bound relations are closed under composition, and therefore R k is a difference bound relation, for any k ≥ 0.

We now recall a result from [START_REF] Bozga | Flat parametric counter automata[END_REF], namely that the k-th iteration R k of a difference bound relation R is equivalent to a finite set of linear inequation systems. As a result, R * is directly definable in Presburger arithmetic.

The constraint graph G k , representing the matrix of R k (i.e. R(x (0) , x (1) ) ∧ . . . ∧ R(x (k-1) , x (k) )) is the graph composed of k connected copies of the constraint graph of R. x (l) and, for all 1 ≤ i, j ≤ n, for all 0 ≤ l < k:

x (l) i a -→ x (l) j ⇐⇒ x i -x j ≤ a ∈ AP (R) x (l) i a -→ x (l+1) j ⇐⇒ x i -x ′ j ≤ a ∈ AP (R) x (l+1) i a -→ x (l) j ⇐⇒ x ′ i -x j ≤ a ∈ AP (R) x (l+1) i a -→ x (l+1) j ⇐⇒ x ′ i -x ′ j ≤ a ∈ AP (R)
In order to compute R k , for a given k ≥ 1, we must first compute the closure of the constraint graph from definition 6, in order to remove the intermediary nodes x (l) , 0 < l < k. According to Theorem 2, R k is a difference bound constraint, and the graph obtained after eliminating the intermediary variables is the constraint graph of R k . Therefore, R k can be written as a conjunction of the form:

1≤i,j≤n x (0) i -x (0) j ≤ min k {x (0) i ; x (0) j } ∧ 1≤i,j≤n x (0) i -x (k) j ≤ min k {x (0) i ; x (k) j } ∧ 1≤i,j≤n x (k) i -x (k) j ≤ min k {x (k) i ; x (k) j } ∧ 1≤i,j≤n x (k) i -x (0) j ≤ min k {x (k) i ; x (0) j }
where min k {x ; y} ∈ Z ∪ {±∞} denotes the value of the minimal path between nodes x and y in the constraint graph of R k3 . Notice that R k is satisfiable if and only there are no cycles of negative cost within the constraint graph from Definition 6 (cf. property 1). Since this graph is obtained by connecting k copies of R, if there is a negative cycle in the graph, there will also be a negative cycle that goes through x

(0) i , for some 1 ≤ i ≤ n.
The main result of [START_REF] Comon | Multiple counters automata, safety analysis and presburger arithmetic[END_REF][START_REF] Bozga | Flat parametric counter automata[END_REF] is that min k {x

(0) i ; x (0) j }, min k {x (0) i ; x (k) j }, min k {x (k) i ; x (k) j } and min k {x (k) i ; x (0) j } are Presburger definable functions. Example 1. Let R(x, x ′ ) be the following DBM relation, where x = {x 1 , x 2 , x 3 }: R(x, x ′ ) = x 1 -x ′ 1 ≤ -1 ∧x 1 -x 2 ≤ 0 ∧x 2 -x 3 ≤ 7 ∧x ′ 2 -x 2 ≤ -2 ∧x 3 -x ′ 3 ≤ -3
Figure 1 shows the graph representation of R k , for k = 1, 2, 3 . . .. By computing the minimal weight paths in the constraint graph of R k for k = 1, 2, 3, ... we obtain, e.g. :

R (1) (x, x ′ ) = x1 -x2 ≤ 0 ∧ x1 -x3 ≤ 7 ∧ x2 -x3 ≤ 7 ∧ x1 -x ′ 1 ≤ -1 ∧ x1 -x ′ 3 ≤ 4 ∧ x2 -x ′ 3 ≤ 4 ∧ x3 -x ′ 3 ≤ -3 ∧ x ′ 2 -x2 ≤ -2 ∧ x ′ 2 -x3 ≤ 5 ∧ x ′ 2 -x ′ 3 ≤ 2 R (2) (x, x ′ ) = x1 -x2 ≤ -3 ∧ x1 -x3 ≤ 4 ∧ x2 -x3 ≤ 7 ∧ x1 -x ′ 1 ≤ -2 ∧ x1 -x ′ 3 ≤ -2 ∧ x2 -x ′ 3 ≤ 1 ∧ x3 -x ′ 3 ≤ -6 ∧ x ′ 2 -x2 ≤ -4 ∧ x ′ 2 -x3 ≤ 3 ∧ x ′ 2 -x ′ 3 ≤ -3 R (3) (x, x ′ ) = x1 -x2 ≤ -6 ∧ x1 -x3 ≤ 1 ∧ x2 -x3 ≤ 7 ∧ x1 -x ′ 1 ≤ -3 ∧ x1 -x ′ 3 ≤ -8 ∧ x2 -x ′ 3 ≤ -2 ∧ x3 -x ′ 3 ≤ -9 ∧ x ′ 2 -x2 ≤ -6 ∧ x ′ 2 -x3 ≤ 1 ∧ x ′ 2 -x ′ 3 ≤ -8 R (n≥4) (x, x ′ ) = x1 -x2 ≤ -3 -3(n -2) ∧ x1 -x3 ≤ -2 -3(n -4) ∧ x2 -x3 ≤ 7 ∧ x1 -x ′ 1 ≤ -n ∧ x1 -x ′ 3 ≤ -14 -6(n -4) ∧ x2 -x ′ 3 ≤ 1 -3(n -2) ∧ x3 -x ′ 3 ≤ -3n ∧ x ′ 2 -x2 ≤ -2n ∧ x ′ 2 -x3 ≤ 3 -2(n -2) ∧ x ′ 2 -x ′ 3 ≤ -3 -5(n -2)
⊓ ⊔

Octagonal Constraints

This section is dedicated to octagonal constraints. We provide preliminary results that are needed to define transitive closures of octagonal relations, in the next section.

0 7 -1 -1 -1 -1 -1 0 0 0 0 -2 -2 -2 -2 -2 7 7 7 7 -3 -3 -3 -3 -3 . . . . . . x (1) 1 x (0) 2 x (0) 3 x (0) 1 x (2) 1 x (3) 1 x (k-2) 1 x (k) 1 x (k-1) 1 x (1) 2 x
(2) 2

x (k-1) 2

x

(1) 3

x

(2) 3

x

(3) 3 x (k-2) 3 x (k-1) 3 x (k) 2 x (k) 3 x (3) 2 . . . x (k-2) 2 Fig. 1. Graph representation of R k , k = 1, 2, 3, . . . Definition 7. Let x = (x 1 , x 2 , .
.., x n ) be a set of variables ranging over Z. Then a formula φ(x) is an octagonal constraint if it is equivalent to a finite conjunction of terms of the form ±x i ±x j ≤ α i,j , 2x i ≤ β i , or -2x i ≤ δ i , where α i,j , β i , δ i ∈ Z and i = j, 1 ≤ i, j ≤ n.

The name octagon comes from the fact that, in two dimensions, these constraints can be graphically represented by polyhedra with at most eight edges. We represent octagons using the set of variables y = (y 1 , y 2 , . . . , y 2n ), with the convention that y 2i-1 stands for x i and y 2i for -x i , respectively. For instance, the octagonal constraint

x 1 + x 2 = 3 is represented as y 1 -y 4 ≤ 3 ∧ y 2 -y 3 ≤ -3.
If we denote by φ = φ[y 1 /x 1 , y 2 /x 1 , . . . , y 2n-1 /x n , y 2n /x n ], we obtain the following valid entailment: φ(x) → (∃y 2 , y 4 , . . . , y 2n .φ)[x 1 /y 1 , . . . , x n /y 2n-1 ]. Moreover φ(x) ↔ (∃y 2 , y 4 , . . . , y 2n .φ∧

n i=1 y 2i-1 +y 2i = 0)[x 1 /y 1 , . . . , x n /y 2n-1 ].
To handle the y variables in the following, we define ī = i -1, if i is even, and ī = i + 1 if i is odd. Obviously, we have ī = i, for all i ∈ Z, i ≥ 0. The following definition extends the matrix representation from difference bound constraints to octagons: Definition 8. Let x = (x 1 , x 2 , ..., x n ) be a set of variables ranging over Z and φ(x) be an octagonal constraint. Then an octagonal difference bound matrix or octagonal DBM representing φ is an 2n × 2n matrix m such that:

(x i -x j ≤ α i,j ) ∈ AP (φ) ⇐⇒ m 2i-1,2j-1 = m 2j,2i = α i,j (-x i -x j ≤ α i,j ) ∈ AP (φ) ⇐⇒ m 2i,2j-1 = m 2j,2i-1 = α i,j (-x i + x j ≤ α i,j ) ∈ AP (φ) ⇐⇒ m 2i,2j = m 2j-1,2i-1 = α i,j (x i + x j ≤ α i,j ) ∈ AP (φ) ⇐⇒ m 2i-1,2j = m 2j-1,2i = α i,j (2x i ≤ β i ) ∈ AP (φ) ⇐⇒ m 2i-1,2i = β i (-2x i ≤ δ i ) ∈ AP (φ) ⇐⇒ m 2i,2i-1 = δ i
A 2n × 2n octagonal DBM m is said to be coherent if and only if m i,j = mj , ī, for all 1 ≤ i, j ≤ 2n. This property is needed since any constraint x ix j ≤ α, 1 ≤ i, j ≤ n can be represented as both y 2i-1y 2j-1 ≤ α and y 2jy 2i ≤ α. If m is coherent, we denote by

γ Oct (m) = {(v 1 , v 2 , ..., v n ) ∈ Z n | (v 1 , -v 1 , v 2 , -v 2 , ..., v n , -v n ) ∈ γ(m)}
the set of concretizations, i.e. the set of models of the octagonal constraint represented by m. Also, m is said to be consistent if γ Oct (m) = ∅, and inconsistent otherwise. The octagonal graph of an octagonal constraint φ(x), represented as a difference bound constraint φ(y) with DBM m, has vertices y, and an arc labeled by m i,j between y i and y j if and only if m i,j < ∞, 1 ≤ i, j ≤ 2n. Definition 9. A consistent coherent 2n × 2n DBM m in Z is said to be tightly closed if and only if the following hold:

1. m i,i = 0 ∀1 ≤ i ≤ 2n 2. m i, ī is even ∀1 ≤ i ≤ 2n 3. m i,j ≤ m i,k + m k,j ∀1 ≤ i, j, k ≤ 2n 4. m i,j ≤ (m i, ī + mj ,j )/2 ∀1 ≤ i, j ≤ 2n
The intuition behind the last point is the following: since y iyī ≤ m i, ī and yjy j ≤ mj ,j and y i = -yī, y j = -yj if we interpret the DBM m as an octagonal constraint, we have that

y i ≤ ⌊ m i, ī 2 ⌋ and -y j ≤ ⌊ mj ,j 2 ⌋, hence y i -y j ≤ ⌊ m i, ī 2 ⌋ + ⌊ mj ,j
2 ⌋. Therefore, the tightening of an octagonal DBM has to come up with values for m i,j that are smaller than ⌊

m i, ī 2 ⌋ + ⌊ mj ,j 2 ⌋, for all 1 ≤ i, j ≤ 2n.
The tight closure of a 2n × 2n octagonal DBM m is an octagonal DBM m * t such that γ Oct (m) = γ Oct (m * t ), and m * t is tightly closed. By Theorem 7 in [START_REF] Miné | The octagon abstract domain[END_REF], we know that the tight closure of a consistent octagonal DBM exists and is unique.

We now characterize the consistency of octagonal constraints using their DBM representations. The following result, proved in [START_REF] Bagnara | An improved tight closure algorithm for integer octagonal constraints[END_REF], is crucial for the developments of the next section, therefore we cite (a slightly modified version of) it here: Theorem 1 ( [START_REF] Bagnara | An improved tight closure algorithm for integer octagonal constraints[END_REF]). Let m be a 2n × 2n octagonal DBM consistent and coherent, and let m * be its closure. Suppose that

m * i, ī 2 + m * ī,i 2 ≥ 0, for all 1 ≤ i ≤ 2n.
Then γ Oct (m) = ∅ and the DBM m T defined as:

m T i,j = min m * i,j , m * i, ī 2 + m * j,j 2
for all 1 ≤ i, j ≤ 2n, is the tight closure of m.

It follows immediately that an octagonal constraint is consistent if and only if (1) its constraint graph representation does not contain negative weight cycles and moreover, (2) if the sum between the halved weights of the minimal paths from any node x i to xī, and from xī and x i is positive. Notice that the lack of negative weight cycles alone is not enough, as shown by the following example.

Example 2. Let φ(x, y, z, x ′ , y ′ , z ′ ) = (x -y ′ ≤ 1) ∧ (y ′ + x ≤ -2) ∧ (-x + z ′ ≤ 1) ∧ (-z ′ -x ≤ 0
) be an octagonal constraint on integer variables x, y, z and x ′ , y ′ , z ′ . The constraint is inconsistent since it implies (-1 ≤ 2x ≤ -1) and thus has no integer solution. However, its corresponding octagonal DBM does not contain negative weight cycles. The tightening step, i.e. replacing each element m i ī with 2 m i ī 2 , exhibits two negative weight cycles, between +x and -x, and between +z ′ and -z ′ .

⊓ ⊔

We are now ready to give the result of existential quantification over variables occurring within octagons. Namely, we prove that, if φ(x 1 , . . . , x n ) is an octagonal constraint, the formula ∃x k .φ, for some 1 ≤ k ≤ n is again an octagonal constraint, and its representation is effectively computable from the constraint graph of φ.

In the following, we assume w.l.o.g. that k = n (if this is not the case, we proceed to reindexing variables). From now on assume that φ(x 1 , . . . , x n ) is consistent, and let m be its 2n × 2n tightly closed coherent octagonal DBM. We denote by β(m) the 2n -2 × 2n -2 matrix from which the 2n -1 and 2n lines and columns have been eliminated. Notice that these correspond to the x n and -x n terms, respectively. Lemma 1. Let m ∈ Z 2n×2n be the coherent and tightly closed difference bound matrix for φ(x 1 , . . . , x n ). Then β(m) is also coherent and tightly closed.

The following theorem proves that octagons are closed under existential quantification.

Theorem 2. Let x = (x 1 , x 2 , . . . , x n ) be a set of variables ranging over Z, φ(x) a consistent octagonal constraint and m the corresponding 2n × 2n octagonal DBM. Then formula ∃x k .φ(x) is equivalent to erasing vertices y 2k and y 2k-1 together with the incident arcs from the graph G(m * t ). Moreover, the resulting graph is also coherent and tightly closed.

Octagonal Relations

This section is dedicated to our main result, the definition of the transitive closure of an octagonal relation in Presburger arithmetic. An octagonal relation is defined in a similar way to a difference bound relation.

Definition 10. Let x = (x 1 , x 2 , . . . , x n ), x ′ = (x ′ 1 , x ′ 2 , . . . , x ′ n )
be sets of variables ranging over Z. Then an octagonal relation R(x, x ′ ) is a relation that can be written as a finite conjunction of terms of the form

±x i ± x j ≤ a i,j ,±x ′ i ± x j ≤ b i,j , ±x i ± x ′ j ≤ c i,j , ±x ′ i ± x ′ j ≤ d i,j , ±2x i ≤ e i,j or ±2x ′ i ≤ f i,j , where a i,j , b i,j , c i,j , d i,j , e i,j , f i,j ∈ Z and 1 ≤ i, j ≤ n, i = j.
The DBM and graph representation of an octagonal relation are defined as in the previous case of difference bound relations. For the rest of this section, let R(x, x ′ ) be an octagonal relation and R(y, y ′ ) be its difference bound representation. Let G be the graph of R, R k be the k-th iteration of R, m k be its DBM, and G k be the graph corresponding to R k , obtained by connecting k copies of G. The formal definition of G k is similar to Definition 6.

In order to compute R k , we need two ingredients. First, we need to check for consistency of the unfolded relation, that is we need to check that γ Oct (m k ) = ∅. Second, we need to obtain the strongest octagonal constraints between y (0) and y (k) in G k . By Theorem 2, all intermediate vertices y (l) , 0 < l < k can be eliminated from G k , and the result is a tightly closed octagonal graph, whose interpretation is equivalent to R k . By Theorem 1, both points require the computation of the values (m k * ) i (l) ,i (l)

and (m k * ) i (l) ,i (l) , for all 1 ≤ i ≤ 2n and 0 ≤ l ≤ k, where the index i (l) refers to the variable y

(l)
i in the DBM representation of R(y, y ′ ). But these values can now be defined by Presburger formulae, using the results from [START_REF] Comon | Multiple counters automata, safety analysis and presburger arithmetic[END_REF][START_REF] Bozga | Flat parametric counter automata[END_REF]. To understand this point, consider the situation depicted in Figure 2.

0 1 l k -1 i īi 1 ..... i m-2 i m-1 k G l-1 G k-l-1 i 2 Fig. 2. Computing paths of weight (m k * ) i (l) ,i (l)
Assume in the following that G k has no cycles of negative weight. The absence of negative cycles can be checked a-priori using e.g. the method in [START_REF] Bozga | Flat parametric counter automata[END_REF]. Since we are aiming at computing minimal weight paths, it is sufficient to consider acyclic paths only 4 . For a fixed 0 < l < k, an acyclic path between the nodes y i (l) and yī (l) , for some 1 ≤ i ≤ 2n, can be decomposed in at most 2n-1 segments starting and ending in y (l) , but not intersecting with y (l) , other than in the beginning and in the end. Moreover, if the path considered is of minimal weight, these segments are of minimal weight as well. Let 1 ≤ i = i 0 , i 1 , . . . , i m = ī ≤ 2n, be a set of pairwise distinct indices, such that

(m * k ) i (l) ,i (l) = m-1 j=0 (m * k ) i (l) j ,i (l) j+1
As in Figure 2, several of the paths y (l) i

; y

(l) i1 , y (l)
i1 ; y

(l) i2 , . . . , y (l)
im-1 ; y (l) ī

will contain only nodes from the set y (≤l) , whereas the rest will contain only nodes from the set y (≥l) . Notice that, the paths from y (≤l) connect the terminal nodes of G l-1 , whereas the paths from y (≥l) connect the initial nodes of G k-l-1 . Therefore, we have:

-(m * k ) i (l) j ,i (l) j+1
= min l-1 {y

(l) ij ; y (l) ij+1 }, if y (l) ij ; y (l)
ij+1 belongs to y (≤l) , and

-(m * k ) i (l) j ,i (l) j+1
= min k-l-1 {y

(0) ij ; y (0) ij+1 }, if y (l) ij ; y (l)
ij+1 belongs to y (≥l) .

But since the min j {x ; y} functions are definable in Presburger arithmetic [START_REF] Comon | Multiple counters automata, safety analysis and presburger arithmetic[END_REF][START_REF] Bozga | Flat parametric counter automata[END_REF], we obtain that (m * k ) i, ī are Presburger definable as well. Since moreover, integer division can be defined in Presburger as

u 2 = v ⇐⇒ 2v ≤ u ≤ 2v + 1
it is possible to encode the consistency check of Theorem 1 by a Presburger formula, and effectively perform the check required by Theorem 1:

1 ≤ i ≤ 2n 0 ≤ l ≤ k $ m * i (l) ,i (l) 2 % + $ m * i (l) ,i (l) 2 % ≥ 0
The formula used to check consistency of R(x, x ′ ) k , k ≥ 0 is of size O(2 2n log 2n ), where n is the number of variables in x. For more details, the interested reader is pointed to [START_REF] Bozga | Iterating octagons[END_REF].

After quantifier elimination, the strongest octagonal relations are obtained by taking the restriction of the tight closure of m k to y (0) and y (k) , according to Theorem 1 and Theorem 2:

yi (0) -yj (k) ≤ (m * kt ) i (0) ,j (k) = min ( (m * kt ) i (0) ,j (k) , $ (m * kt ) i (0) ,i (0) 2 % + $ (m * k t ) j (k) ,j (k) 2
%) yi (k) -yj (0) ≤ (m * kt ) i (k) ,j (0) = min ( (m * kt ) i (k) ,j (0) , $ (m * kt ) i (k) ,i (k) 2 % + $ (m * k t ) j (0) ,j (0) 2
%) yi (0) -yj (0) ≤ (m * kt ) i (0) ,j (0) = min ( (m * k t ) i (0) ,j (0) , $ (m * kt ) i (0) ,i (0) 2 % + $ (m * k t ) j (0) ,j (0) 2
%) yi (k) -yj (k) ≤ (m * kt ) i (k) ,j (k) = min ( (m * kt ) i (k) ,j (k) , $ (m * k t ) i (k) ,i (k) 2 % + $ (m * kt ) j (k) ,j (k) 2
%)

Since we showed that all weights of the form m * i (l) ,i (l) are definable in Presburger arithmetic, it follows that R k is also Presburger definable. Consequently, the transitive closure of R is the Presburger formula ∃k.R k , leading to the following theorem: Theorem 3. The transitive closure of an integer octagonal relation is Presburger definable.

Example 3. Let R(x, x ′ ) be the following octagonal relation, where 3 shows the graph representation of R k for k = 1, 2, 3, . . .. The n-step closure is:

x = {x 1 , x 2 } R(x, x ′ ) = x 1 + x 2 ≤ 5 ∧ x ′ 1 -x 1 ≤ -2 ∧ x ′ 2 -x 2 ≤ -3 Figure
R n (x, x ′ ) = x 1 + x 2 ≤ 5 ∧ x ′ 1 -x 1 ≤ -2n ∧ x ′ 2 -x 2 ≤ -3n∧ x 1 + x ′ 2 ≤ 5 -3n ∧ x 2 + x ′ 1 ≤ 5 -2n ∧ x ′ 1 + x ′ 2 ≤ 5 -5n ⊓ ⊔ . . . . . . . . . . . . y1 y2 y3 y4 -2 -2 -2 -2 -2 -2 -2 -2 -2 -2 -3 -3 -3 -3 -3 -3 -3 -3 -3 -3 5

Implementation and Experience

We have implemented the method for computing transitive closures of octagonal relations in a tool for the analysis of counter automata, that we are currently developing. Table 1 shows several octagonal relations (first column), the closed form of their iteration (second column), the execution times (in milliseconds) needed to compute the closed forms without consistency checks (third column), and with consistency checks (fourth column). Notice that the first relation is the difference bound relation example from Section 2, for which no consistency check is needed.

R R n time (ms) -check +check x1 -x2 ≤ 0 x2 -x3 ≤ 7 x1 -x ′ 1 ≤ -1 x3 -x ′ 3 ≤ -3 x ′ 2 -x2 ≤ -2 x1 -x2 ≤ 3 -3n x1 -x3 ≤ 10 -3n x2 -x3 ≤ 7 x1 -x ′ 1 ≤ -n x1 -x ′ 3 ≤ 10 -6n x2 -x ′ 3 ≤ 7 -3n x3 -x ′ 3 ≤ -3n x ′ 2 -x2 ≤ -2n x ′ 2 -x3 ≤ 7 -2n x ′ 2 -x ′ 3 ≤ 7 -5n 383 x1 + x ′ 2 ≤ -1 -x2 -x ′ 1 ≤ -2 n ≡2 0 ⇒ x1 -x ′ 1 ≤ -3n/2 n ≡2 0 ⇒ x ′ 2 -x2 ≤ -3n/2 n ≡2 1 ⇒ x1 + x ′ 2 ≤ -1 -3(n -1)/2 n ≡2 1 ⇒ -x2 -x ′ 1 ≤ -2 -3(n -1)/2 75 116 x1 -x ′ 1 ≤ 0 x1 + x ′ 2 ≤ -1 -x2 -x ′ 1 ≤ -2 n ≡2 0 ⇒ x1 -x ′ 1 ≤ -3n/2 n ≡2 1 ⇒ x1 -x ′ 1 ≤ -3(n -1)/2 n ≡2 0 ⇒ x ′ 2 -x2 ≤ -3n/2 n ≡2 1 ∧ n ≥ 3 ⇒ x ′ 2 -x2 ≤ -3 -3(n -3)/2 n ≡2 0 ∧ n ≥ 2 ⇒ x1 + x ′ 2 ≤ -1 -3(n -2)/2 n ≡2 1 ⇒ x1 + x ′ 2 ≤ -1 -3(n -1)/2 n ≡2 0 ∧ n ≥ 2 ⇒ -x2 -x ′ 1 ≤ -2 -3(n -2)/2 n ≡2 1 ⇒ -x2 -x ′ 1 ≤ -2 -3(n -1)/2 81 124 x1 + x2 ≤ 5 x ′ 1 -x1 ≤ -2 x ′ 2 -x2 ≤ -3 x1 + x2 ≤ 5 x ′ 1 -x1 ≤ -2n x ′ 2 -x2 ≤ -3n x1 + x ′ 2 ≤ 5 -3n x2 + x ′ 1 ≤ 5 -2n x ′ 1 + x ′ 2 ≤ 5 -5n 11686 38729
Table 1. Experimental results

The tool, called FLATA, is implemented in Java (TM) , and it is currently available at [10]. The execution times from Table 1 are relative to a Intel (R) Xeon (TM) CPU 3.00GHz equipped with 1 Gb of RAM. As an approximate upper bound, we can in principle handle difference bound relations on ten counters, and/or octagons on five counters, within reasonable execution times (≈ 10 min).

The fourth example from Table 1 takes significantly more time than the previous ones. This happens because in the pre-processing step we need to replace x 1 + x 2 ≤ 5 with x 1x ′ 3 ≤ 5 ∧ x ′ 3 + x 2 ≤ 0, where x 3 is a fresh variable, thus having to deal with a DBM relation with 6 counters.

Conclusions

We have considered the problem of computing transitive closures of octagonal relations, which are finite conjunctions of atomic formulae of the form ±x±y ≤ c, where x and y are possibly primed integer variables. We show that the k-th iteration of such a relation has a Presburger definable closed form. As a consequence, the transitive closure is also Presburger definable. This result enlarges the class of non-deterministic counter automata for which the reachability problem is decidable, to counter automata with octagonal transition relations. This result is expected to have an impact in the fields of software verification, as abstract models of software systems are described using octagons. We have implemented our method in a tool for the analysis of counter automata, and report on a number of experiments.

Definition 6 .

 6 Let x = (x 1 , x 2 , ..., x n ), x ′ = (x ′ 1 , x ′ 2 , ..., x ′ n ) besets of variables ranging over Z and R(x, x ′ ) a DBM relation. The constraint graph of R k has as set of vertices k l=0

Fig. 3 .

 3 Fig. 3. Graph representation of R k , k = 1, 2, 3, . . .

This is needed because we consider the counters to range over integer numbers. If rationals or real number would have been used, this would not have been needed.

In general, difference bound constraints are not closed under universal quantification, however finite disjunctions of difference bound constraints are, since the negation of a difference bound constraint is always equivalent to a finite disjunction of difference constraints.

If min k {x ; y} = -∞ the constraint is logically equivalent to false, and if min{x ; y} = ∞ it is logically equivalent to true.

If a path has a cycle of a positive weight, it cannot be minimal.


Alexandre Cordier

Marcel Klinksiek

Christoph Held

Julien Legros
email: julien.legros@univ-rouen.fr

Sébastien Leveneur
email: sebastien.leveneur@insa-rouen.fr

Biocatalyst and continuous microfluidic reactor for an intensified production of n-butyl levulinate: kinetic model assessment

Keywords: Kinetic modeling, ePC-SAFT, Enzyme, Ping-Pong mechanism

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

Introduction

Lignocellulosic biomass (LCB) is seen as the best alternative to fossil raw materials to make the chemical and fuel industries sustainable [START_REF] Singhvi | Lignocellulosic biomass: Hurdles and challenges in its valorization[END_REF][START_REF] Guragain | Renewable Biomass Utilization: A Way Forward to Establish Sustainable Chemical and Processing Industries[END_REF][START_REF] Zhu | New synthetic approaches to biofuels from lignocellulosic biomass[END_REF]. Compared to first-generation biomass, LCB is not competing with the food sector, avoiding the fuel versus food dilemma. The chemistry of this LCB valorization focuses on the production of platform molecules [START_REF] Lu | Reductive Catalytic Depolymerization of Semi-industrial Wood-Based Lignin[END_REF][START_REF] Schutyser | Chemicals from lignin: an interplay of lignocellulose fractionation, depolymerisation, and upgrading[END_REF][START_REF] Rinaldi | Paving the Way for Lignin Valorisation: Recent Advances in Bioengineering, Biorefining and Catalysis[END_REF] and lignin valorization [START_REF] Lu | Reductive Catalytic Depolymerization of Semi-industrial Wood-Based Lignin[END_REF][START_REF] Schutyser | Chemicals from lignin: an interplay of lignocellulose fractionation, depolymerisation, and upgrading[END_REF][START_REF] Rinaldi | Paving the Way for Lignin Valorisation: Recent Advances in Bioengineering, Biorefining and Catalysis[END_REF].

The platform molecule levulinic acid production has gained much interest these last years [START_REF] Badgujar | A review on catalytic synthesis of energy rich fuel additive levulinate compounds from biomass derived levulinic acid[END_REF][START_REF] Di Menno Di Bucchianico | Production of levulinic acid and alkyl levulinates: a process insight[END_REF][START_REF] Démolis | Synthesis and applications of alkyl levulinates[END_REF]. Levulinic acid or alkyl levulinates are starting materials for producing another platform molecule γ-valerolactone (GVL) [START_REF] Yan | Conversion of levulinic acid and alkyl levulinates into biofuels and high-value chemicals[END_REF][START_REF] Wang | Application of the concept of Linear Free Energy Relationships to the hydrogenation of levulinic acid and its corresponding esters[END_REF][START_REF] Capecci | Solvent effect on the kinetics of the hydrogenation of n-butyl levulinate to γ-valerolactone[END_REF][START_REF] Alonso | Gamma-valerolactone, a sustainable platform molecule derived from lignocellulosic biomass[END_REF][START_REF] Negahdar | Kinetic analysis of the catalytic hydrogenation of alkyl levulinates to γ-valerolactone[END_REF][START_REF] Wright | Development of heterogeneous catalysts for the conversion of levulinic acid to γ-valerolactone[END_REF][START_REF] Capecci | Bayesian Statistics to Elucidate the Kinetics of γ-Valerolactone from n-Butyl Levulinate Hydrogenation over Ru/C[END_REF][START_REF] Delgado | Kinetic model assessment for the synthesis of γ-valerolactone from n-butyl levulinate and levulinic acid hydrogenation over the synergy effect of dual catalysts Ru/C and Amberlite IR-120[END_REF]. Alkyl levulinates can be used directly as blending components for biodiesel or as a fuel oxygenate additives [START_REF] Peng | Catalytic upgrading of renewable furfuryl alcohol to alkyl levulinates using AlCl3 as a facile, efficient, and reusable catalyst[END_REF][START_REF] Peixoto | Production of ethyl levulinate fuel bioadditive from 5-hydroxymethylfurfural over sulfonic acid functionalized biochar catalysts[END_REF], and can also find applications as additives, solvents, and intermediates in fine chemistry [START_REF] Démolis | Synthesis and applications of alkyl levulinates[END_REF]. In the study of Christensen et al., they showed that butyl levulinate (BL) improved conductivity, cold flow properties, and lubricity of diesel fuel and reduced its vapor pressure [START_REF] Christensen | Properties and Performance of Levulinate Esters as Diesel Blend Components[END_REF]. Moreover, it was found that BL remains in solution with diesel down to the fuel cloud point and has more compatibility with elastomers, compared to ethyl levulinate, which tends to separate from diesel at a temperature below 0 °C and results to be more corrosive [START_REF] Christensen | Properties and Performance of Levulinate Esters as Diesel Blend Components[END_REF]. Frigo et al., demonstrated that diesel fuel blended with a mixture of dibutyl ether and BL could reduce particulate emissions without changing engine power efficiency or increasing the NOx emission [START_REF] Frigo | Utilisation of advanced biofuel in CI internal combustion engine[END_REF].

Alkyl levulinates can be produced via the alcoholysis of sugar monomers or LA esterification [START_REF] Badgujar | A review on catalytic synthesis of energy rich fuel additive levulinate compounds from biomass derived levulinic acid[END_REF][START_REF] Di Menno Di Bucchianico | Production of levulinic acid and alkyl levulinates: a process insight[END_REF][START_REF] Di Menno Di Bucchianico | Role of solvent the production of butyl levulinate from fructose[END_REF]. The latter route can be done via homogenous [START_REF] Harwardt | Conceptual Design of a Butyl-levulinate Reactive Distillation Process by Incremental Refinement[END_REF], heterogeneous [START_REF] Russo | Kinetic study of Amberlite IR120 catalyzed acid esterification of levulinic acid with ethanol: From batch to continuous operation[END_REF][START_REF] Tejero | Esterification of levulinic acid with butanol over ion exchange resins[END_REF][START_REF] Peixoto | Acid functionalized coal fly ashes: New solid catalysts for levulinic acid esterification[END_REF][START_REF] Dharne | Esterification of levulinic acid to n-butyl levulinate over heteropolyacid supported on acid-treated clay[END_REF][START_REF] Nandiwale | Esterification of Renewable Levulinic Acid to n-Butyl Levulinate over Modified H-ZSM-5[END_REF][START_REF] Nandiwale | Environmentally benign catalytic process for esterification of renewable levulinic acid to various alkyl levulinates biodiesel[END_REF][START_REF] Maheria | Esterification of Levulinic Acid to n-Butyl Levulinate Over Various Acidic Zeolites[END_REF] or enzymatic catalysis [START_REF] Song | The effect of enzyme loading, alcohol/acid ratio and temperature on the enzymatic esterification of levulinic acid with methanol for methyl levulinate production: A kinetic study[END_REF][START_REF] Yadav | Kinetic Modeling of Immobilized Lipase Catalysis in Synthesis of n-Butyl Levulinate †[END_REF][START_REF] Emel'yanenko | Renewable platform chemicals: Thermochemical study of levulinic acid esters[END_REF][START_REF] Bhavsar | n-Butyl levulinate synthesis using lipase catalysis: comparison of batch reactor versus continuous flow packed bed tubular microreactor[END_REF][START_REF] Zhou | Conversion of levulinic acid into alkyl levulinates: Using lipase immobilized on meso-molding three-dimensional macroporous organosilica as catalyst[END_REF][START_REF] Zhai | Enzymatic synthesis of a novel solid-liquid phase change energy storage material based on levulinic acid and 1,4-butanediol[END_REF][START_REF] Badgujar | Lipase as a green and sustainable material for production of levulinate compounds: State of the art[END_REF]. The use of heterogeneous catalysts such as resins, zeolite or immobilized enzymes should be favored to avoid additional separation stages [START_REF] Harwardt | Conceptual Design of a Butyl-levulinate Reactive Distillation Process by Incremental Refinement[END_REF]. From a chemical engineering viewpoint, using a continuous reactor for biomass valorization should be favored for large production [START_REF] Gérardy | Continuous Flow Upgrading of Selected C2-C6 Platform Chemicals Derived from Biomass[END_REF].

Microreactor or microfluidic technology has raised interest in the scientific community because concentration and temperature gradients are reduced.

Production of n-butyl levulinate from LA esterification over lipase catalysis has been scarcely studied [START_REF] Yadav | Kinetic Modeling of Immobilized Lipase Catalysis in Synthesis of n-Butyl Levulinate †[END_REF][START_REF] Bhavsar | n-Butyl levulinate synthesis using lipase catalysis: comparison of batch reactor versus continuous flow packed bed tubular microreactor[END_REF]. Yadav and Borkar [START_REF] Bhavsar | n-Butyl levulinate synthesis using lipase catalysis: comparison of batch reactor versus continuous flow packed bed tubular microreactor[END_REF] developed a kinetic model of LA esterification by butanol over Novozym®435, without considering the reversibility of this reaction or the temperature effect. Bhavsar and Yadav [START_REF] Bhavsar | n-Butyl levulinate synthesis using lipase catalysis: comparison of batch reactor versus continuous flow packed bed tubular microreactor[END_REF] showed that a continuous packed bed reactor with immobilized Novozym can be used to produce BL.

There is a need to intensify this reaction from an industrial [START_REF] Zhang | Microfluidic reactor with immobilized enzymefrom construction to applications: A review[END_REF] and fundamental standpoint. The use of a microreactor enables to operate in the absence of gradients allowing to use plug-flow model and thus simplifying the kinetic modeling stage [START_REF] Taylor | Modern advancements in continuous-flow aided kinetic analysis[END_REF].

The developed models proposed in the literature for the esterification of carboxylic acid over enzyme do not consider the temperature effect, and the knowledge of this effect on kinetic constants is mandatory for an industrial scaleup. In this manuscript, using microfluidic technology, kinetic models for the synthesis of BL over commercial Novozym®435 were developed and assessed at different temperatures.

Experimental section

Chemicals

All the chemicals were used as provided, without further purification. Butan-1-ol (wt% ≥ 99.9%), levulinic acid (wt% ≥ 99%), (trimethylsilyl)diazomethane (2.0 M in hexane) and the commercial supported Candida Antarctica lipase B (CAL-B), Novozym 435 (5000 U.mg -1 ), were purchased from Sigma-Aldrich. Dichloromethane stabilized by ethanol was purchased from CARLO ERBA Reagents. Deionized water from the Aquadem™ system (Veolia) was used.

Analytical method

The aliquots taken at each residence time (t R ) for the analysis of butyl levulinate concentration were analyzed on a Thermo Scientific™ TRACE™ 1310 GC-FID equipped with an apolar column (DB-5MS, 30 m × 0.250 mm ID × 0.250 µm film thickness). The initial temperature of the analysis method was set at 50 °C for 2 min to reach 250 °C with a temperature rate of 25 °C/min. The aliquots were diluted in dichloromethane, and an excess of trimethylsilyl diazomethane as carboxylic acid scavenger was added (to protect the GC column from corrosion); 1 µL of the resulting solution was injected into the GC.

High field 1 H NMR analyses were performed on a 300 MHz Bruker Spectrospin spectrometer. Chemical shifts (δ) are given with regard to TMS using residual CHCl3 solvent as an internal reference.

Procedure for the biocatalyzed esterification of levulinic acid with butanol in a flow reactor

The Luer-lock syringe, filled with levulinic acid, butanol and water, was connected to a preheating loop (PTFE tubing, ID = 1.59 mm, L = 126 cm) before the packed-bed reactor with an internal diameter of 6.6 mm. The packed-bed reactor was composed of an Omnifit TM column filled with Novozym®435 immerged in a thermostated bath (Fig. 1). Fig. 1. Process for the butyl levulinate synthesis with Omnifit cartridge.

The kinetic measurements were performed. First, the flow system was filled with the solution at a flow rate of 5 mL/min. Then the flow-rate was increased step by step to reach the desired t R . The collection of aliquots began after running the system for 2 mL of reaction product collected, in order to reach the steady-state. The volume of the microreactor is related to the amount of the Novozym®435: 50 mg, 100 mg and 150 mg of supported CAL-B were packed 195 µL, 390 µL and 520 µL, respectively. The kinetic monitoring was performed as described in Table 1. Repeatability was assessed by collecting three samples after the steady-state, and it was found that the standard deviation was lower than 0.05 mol.L -1 . Furthermore, one run was repeated two times (Run 4) on different days to verify the repeatability of the whole system (Fig. S1.1). Fig. S1.1 shows that the repeatability is good. was added, and a 1 H NMR analysis was performed. The obtained spectrum was compared to those of authentic samples of butyl levulinate, n-butanol and levulinic acid, as well as those from PMMA and methyl methacrylate from the literature [START_REF] Hajifatheali | The synthesis of Nmethylbis[2-(dodecylthio)ethyl]amine (SNS) and investigation of its efficiency as new mononuclear catalyst complex in copper-based ATRP[END_REF].

T°C

The deactivation test was carried out at 65 °C with 100 mg of catalyst and inlet levulinic acid concentration of 2.61 mol.L -1 . The outlet concentration of BL was followed with time-on-stream at 65°C.

Results

External and internal mass transfer evaluation

Besides the effect of mass transfer resistance on the kinetics, the presence or absence of flow maldistribution should be determined [START_REF] Doraiswamy | Laboratory Catalytic Reactors[END_REF][START_REF] Commenge | Methodology for multi-scale design of isothermal laminar flow networks[END_REF][START_REF] Hassankiadeh | Dynamic Assessment and Optimization of Catalytic Hydroprocessing Process: Sensitivity Analysis and Practical Tips[END_REF]. According to Doraiswamy and Tajbl [START_REF] Doraiswamy | Laboratory Catalytic Reactors[END_REF], if the ratio reactor diameter on particle diameter is higher than 4, then they conclude that there is a proper liquid distribution with no channeling. In this system, this ratio is higher than 10. Thus, we concluded the absence of flow maldistribution.

To evaluate the influence of both effect, the same methodology presented by Leveneur et al. [START_REF] Leveneur | Green process technology for peroxycarboxylic acids: Estimation of kinetic and dispersion parameters aided by RTD measurements: Green synthesis of peroxycarboxylic acids[END_REF] was applied.

External mass transfer for each experiment was evaluated throughout the coefficient fe (Equation ( 1)) defined by Villermaux [START_REF] Villermaux | Génie de la réaction chimique[END_REF]. If fe is lower than 5%, then the external mass transfer is negligible.

𝑓 𝑒 = 𝑟 𝑂𝑏𝑠 ̅̅̅̅̅ • 𝐿 𝑘 𝐷 • 𝐶 𝑏 (1)
where, L is the ratio particle volume (VP) on the external particle surface (AP), 𝑟 𝑂𝑏𝑠 ̅̅̅̅̅ is the initial observed rate of esterification, 𝑘 𝐷 is the mass transfer coefficient and 𝐶 𝑏 the concentration in the bulk phase. In Equation ( 1), the concentration of LA in the bulk phase was used because butanol is in excess. The mean particle size of Novozym®435 is equal to 0.65 mm, and thus L is equal to 1.08 . 10 -4 m [START_REF] Ravelo | Synthesis of ibuprofen monoglyceride in solventless medium with novozym®435: Kinetic analysis[END_REF].

The mass transfer coefficient 𝑘 𝐷 can be estimated via the Sherwood number (Sh) expressed as

𝑆ℎ = 𝑘 𝐷 • 𝑑 𝑃 ̅̅̅ 𝐷 = 2 + 1.8 • 𝑅𝑒 𝑃 1 2 • 𝑆𝑐 1 3 (2)
where, 𝑑 𝑃 is the mean diameter of the catalyst particle, and Sc stands for the Schmidt number expressed as

𝑆𝑐 = 𝜇 𝑓 𝜌 𝑓 • 𝐷 (3)
D is the molecular diffusion coefficient of LA in butanol calculated by the Wilke-Chang equation [START_REF] Wilke | Correlation of diffusion coefficients in dilute solutions[END_REF]. For instance, the molecular diffusion of LA in butanol at 50°C was found to be 9.02 . 10 -10 m 2 .s -1 . The terms 𝜇 𝑓 and 𝜌 𝑓 represent the viscosity and density of the fluid, i.e., butanol. These physicochemical properties were calculated from Ariba et al.

work [START_REF] Ariba | Physicochemical Properties for the Reaction Systems: Levulinic Acid, Its Esters, and γ-Valerolactone[END_REF]. The fe values were found to be lower than 5% for each experiment showing the absence of external mass transfer.

The internal mass transfer effect was evaluated via the Thieles modulus number 𝜙 𝑆 defined by Equation ( 4) [START_REF] Villermaux | Génie de la réaction chimique[END_REF].

𝜙 𝑆 2 = 𝑟 𝑂𝑏𝑠 ̅̅̅̅̅ • 𝐿 2 𝐷 𝑒 . 𝐶 𝑆 (4)
If 𝜑 𝑆 ′ is lower than 0.1, hence internal mass transfer can be assumed to be negligible.

CS is the LA concentration at the particle surface and in this study CS=CB because there is no external mass transfer. The term De represents the effective diffusion coefficient defined as

𝐷 𝑒 = 𝜀 𝑃 •𝜎 𝜏
• 𝐷, where 𝜀 𝑃 , 𝜎 and 𝜏 represent the porosity, constriction factor and tortuosity of the particle, respectively. From Ravelo et al. [START_REF] Ravelo | Synthesis of ibuprofen monoglyceride in solventless medium with novozym®435: Kinetic analysis[END_REF],

Novozym ® 435 porosity is equal to 0.5. The tortuosity and constriction factor values were fixed to 6 and 1 [START_REF] Chesterfield | Production of biodiesel via ethanolysis of waste cooking oil using immobilised lipase[END_REF]. Based on the 𝜑 𝑆 ′ values of each experiment, the internal mass transfer can be assumed to be negligible.

Novozym 435 stability study

Novozym 435 degradability test (S2) is an essential study because its acrylic support/matrix tends to dissolve in many organic solvents [START_REF] Ortiz | Novozym 435: The "perfect" lipase immobilized biocatalyst?[END_REF][START_REF] Kowalczykiewicz | Rotating bed reactor packed with heterofunctional structured silica-supported lipase. Developing an effective system for the organic solvent and aqueous phase reactions[END_REF][START_REF] Zhao | Migration of reactive trace compounds from Novozym® 435 into organic solvents and ionic liquids[END_REF]. By spectrum comparison of the mixture with a 1 H NMR analysis of different isolated products (butyl levulinate, n-butanol and levulinic acid), the mixture just contains n-butanol, levulinic acid and butyl levulinate and absolutely no traces of PMMA (or associate compound)

were detected. Thus, our measurements are perfectly reliable with no interference from external chemicals.

Fig. 2 shows the BL concentration at the outlet versus time-on-stream. From Fig. 2, one can notice that enzyme deactivation can be neglected during the experiment.

Temperature effect

Fig. 3 shows the effect of temperature on the ratio obtained from Runs 1-4 were compared, because there were carried out in the same operating conditions, except for the reaction temperature (Table 2). As expected, the kinetics of esterification increases with temperature.

These data can also be used to evaluate the effect of mass transfer. The natural logarithm of the initial rate constants versus 1/T was plotted (Fig. S3.1). From Fig. S3.1, the linearity between the natural logarithm of the initial rate constants and 1/T confirms the absence of mass transfer resistance.

Molar ratio effect

The effect of

[𝐵𝑢𝑂𝐻] 𝑖𝑛𝑙𝑒𝑡 [𝐿𝐴] 𝑖𝑛𝑙𝑒𝑡
can affect the thermodynamics and kinetics of esterification [START_REF] Song | The effect of enzyme loading, alcohol/acid ratio and temperature on the enzymatic esterification of levulinic acid with methanol for methyl levulinate production: A kinetic study[END_REF]. Different experiments carried out in similar operating conditions, except the ratio

[𝐵𝑢𝑂𝐻] 𝑖𝑛𝑙𝑒𝑡 [𝐿𝐴] 𝑖𝑛𝑙𝑒𝑡
, were compared (Table 2). One can notice that when

[𝐵𝑢𝑂𝐻] 𝑖𝑛𝑙𝑒𝑡

[𝐿𝐴] 𝑖𝑛𝑙𝑒𝑡 is equal to 3:1 or 5:1, then the reaction rates and equilibrium values are similar. For a

[𝐵𝑢𝑂𝐻] 𝑖𝑛𝑙𝑒𝑡

[𝐿𝐴] 𝑖𝑛𝑙𝑒𝑡 equal to 1:1, there is a deviation when the reaction reaches the equilibrium. Run16_Ratio BuOH:LA 5:1

Evaluation of equilibrium constant

The PC-SAFT equation of state first published by Gross and Sadowski [START_REF] Gross | Perturbed-Chain SAFT: An Equation of State Based on a Perturbation Theory for Chain Molecules[END_REF] expresses the residual Helmholtz energy 𝑎 𝑟𝑒𝑠 as shown in Equation 5.

𝑎 𝑟𝑒𝑠 = 𝑎 ℎ𝑐 + 𝑎 𝑑𝑖𝑠𝑝 + 𝑎 𝑎𝑠𝑠𝑜𝑐 (5)
Thereby the hard-chain (𝑎 ℎ𝑐 ) reference system represents the repulsive interactions between the molecules. The attractive interactions, such as the dispersion (𝑎 𝑑𝑖𝑠𝑝 ), association (𝑎 𝑎𝑠𝑠𝑜𝑐 ), are described as perturbations of the reference system. More details on the modeling procedure and parameters used can be found in Supplementary Information (S4). Thermodynamic modeling of the reaction equilibrium is based on the temperature and pressure-dependent equilibrium constant Kth. It is calculated by the reacting agent concentrations in the equilibrium and activity coefficients according to Equation 6.

Based on Equation 5, we can calculate the temperature and pressure-dependent equilibrium constant Kth.

K th (T, p) = K eq (T, p, x) ⋅ K γ (T, p, x) = ∏ (x i ⋅ γ i ) ν i i (6)
where, Keq is determined from the experimental equilibrium concentrations and the activity coefficients are obtained by PC-SAFT [START_REF] Gross | Perturbed-Chain SAFT: An Equation of State Based on a Perturbation Theory for Chain Molecules[END_REF][START_REF] Lemberg | Predicting the Solvent Effect on Esterification Kinetics[END_REF][START_REF] Altuntepe | Thermodynamics of enzyme-catalyzed esterifications: II. Levulinic acid esterification with short-chain alcohols[END_REF][START_REF] Held | Thermodynamics of Bioreactions[END_REF][START_REF] Tumakaka | Application of the Perturbed-Chain SAFT equation of state to polar systems[END_REF][START_REF] Fuchs | Solubility of Amino Acids: Influence of the pH value and the Addition of Alcoholic Cosolvents on Aqueous Solubility[END_REF]. The activity coefficient of each reactant i in the mixture is calculated from the ratio of the fugacity coefficients in the mixture and of the fugacity coefficient of the pure component.

𝛾 𝑖 = φ i (T,p,𝑥) φ 0i (T,p,x i =1) (7)
The equilibrium constant Kth for each temperature enables the calculation of equilibrium concentrations at different conditions, i.e., molar ratios. The equilibrium constant Kth was calculated based on the experiments with a molar ratio of BuOH:LA 3:1 (Table 3).

Discussion

Ping-Pong Models

Several authors showed that the Ping-Pong Bi-Bi mechanism can be used for the esterification reactions [START_REF] Yadav | Kinetic Modeling of Immobilized Lipase Catalysis in Synthesis of n-Butyl Levulinate †[END_REF][START_REF] Yadav | Synergism between microwave irradiation and enzyme catalysis in transesterification of ethyl-3-phenylpropanoate with nbutanol[END_REF][START_REF] Bansode | Evaluation of reaction parameters and kinetic modelling for Novozym 435 catalysed synthesis of isoamyl butyrate[END_REF][START_REF] Zulkeflee | Modelling of immobilized Candida rugosa lipase catalysed esterification process in batch reactor equipped with temperature and water activity control system[END_REF][START_REF] Mitchell | Looking through a new lens: Expressing the Ping Pong bi bi equation in terms of specificity constants[END_REF][START_REF] Lopresto | Kinetic study on the enzymatic esterification of octanoic acid and hexanol by immobilized Candida antarctica lipase B[END_REF][START_REF] Varma | Kinetics of synthesis of butyl butyrate by esterification and transesterification in supercritical carbon dioxide[END_REF][START_REF] Rizzi | A kinetic study of immobilized lipase catalysing the synthesis of isoamyl acetate by transesterification in n-hexane[END_REF]. Several of them showed that alcohol and carboxylic acid can inhibit the enzyme, but none developed a kinetic model considering the reversibility of this reaction and the temperature effect on the kinetic constants.

Fig. 6 shows the Ping-Pong Bi-Bi mechanism for the esterification of LA by butanol.

Fig. 6. Ping-Pong Bi-Bi and inhibition mechanism for the esterification of LA.

From Varma and Madras study [START_REF] Varma | Kinetics of synthesis of butyl butyrate by esterification and transesterification in supercritical carbon dioxide[END_REF], the rate equation can be derived as

𝑟 𝐸𝑠𝑡𝑒𝑟𝑖𝑓𝑖𝑐𝑎𝑡𝑖𝑜𝑛 = 𝑘 𝑓 • 𝑘 𝑏 • [𝐸] 0 2 • ([𝐿𝐴] • [𝐵𝑢𝑂𝐻] - [𝐵𝐿] • [𝑊] 𝐾 𝑒𝑞 ) 𝐷 ( 8
)
where, [𝐿𝐴], [𝐵𝑢𝑂𝐻], [𝐵𝐿] and [𝑊] are the outlet concentrations of levulinic acid, butanol, butyl levulinate and water, respectively. The term [𝐸] 0 stands for the initial concentration of enzyme. The denominator D is expressed as

E LA E-LA E*-W W E* BuOH E*-BuOH E-BL E BL E*-LA E*-BL + LA + BL E*-W E*-BuOH + W + BuOH Inhibition 𝐷 = 𝑘 𝑏 • [𝐸] 0 • [𝐿𝐴] • [𝐵𝑢𝑂𝐻] + 𝑘 𝑏 • [𝐸] 0 • 𝐾 𝐵𝑢𝑂𝐻 • [𝐿𝐴] • (1 + [𝐿𝐴] 𝐾 𝐼𝐿𝐴 ) +𝑘 𝑏 • [𝐸] 0 • 𝐾 𝐿𝐴 • [𝐵𝑢𝑂𝐻] • (1 + [𝐵𝑢𝑂𝐻] 𝐾 𝐼𝐵𝑢𝑂𝐻 ) + 𝑘 𝑓 • [𝐸] 0 • 𝐾 𝑊 𝐾 𝑒𝑞 • [𝐵𝐿] • (1 + [𝐵𝐿] 𝐾 𝐼𝐵𝐿 ) + 𝑘 𝑓 • [𝐸] 0 • 𝐾 𝐵𝐿 𝐾 𝑒𝑞 • [𝑊] • (1 + [𝑊] 𝐾 𝐼𝑊 ) + 𝑘 𝑓 • [𝐸] 0 𝐾 𝑒𝑞 • [𝑊] • [𝐵𝐿] + 𝑘 𝑓 • 𝐾 𝐵𝐿 • [𝐸] 0 𝐾 𝐼𝐼𝐿𝐴 • 𝐾 𝑒𝑞 • [𝑊] • [𝐿𝐴] + 𝑘 𝑏 • 𝐾 𝐿𝐴 • [𝐸] 0 𝐾 𝐼𝐼𝐵𝐿 • [𝐿𝐴] • [𝐵𝐿] + 𝑘 𝑏 • 𝐾 𝐿𝐴 • [𝐸] 0 𝐾 𝐵𝑢𝑂𝐻-𝑊 • [𝑊] • [𝐵𝑢𝑂𝐻] + 𝑘 𝑏 • 𝐾 𝐵𝑢𝑂𝐻 • [𝐸] 0 𝐾 𝐿𝐴-𝐵𝐿 • [𝐿𝐴] • [𝐵𝐿] (9)
The terms 𝑘 𝑓

𝑟 𝐸𝑠𝑡𝑒𝑟𝑖𝑓𝑖𝑐𝑎𝑡𝑖𝑜𝑛 = (𝑘 𝐿𝐴 • [𝐿𝐴] • 𝑘 𝐵𝑢𝑂𝐻 • [𝐵𝑢𝑂𝐻] -𝑘 𝐵𝐿 • [𝐵𝐿] • 𝑘 𝑊 • [𝑊]) • [𝐸] 0 𝐷 ′ (10)
where, the denominator D' is

𝐷 ′ = 𝑘 𝐿𝐴 • [𝐿𝐴] • (1 + [𝑊] 𝐾 𝐼𝐼𝑊 + [𝐵𝑢𝑂𝐻] 𝐾 𝐵𝑢𝑂𝐻 ) + 𝑘 𝐵𝑢𝑂𝐻 • [𝐵𝑢𝑂𝐻] + 𝑘 𝑊 • [𝑊] +𝑘 𝐵𝐿 • [𝐵𝐿] • (1 + [𝑊] 𝐾 𝑊 + [𝐵𝑢𝑂𝐻] 𝐾 𝐼𝐼𝐵𝑢𝑂𝐻 ) (11)
The terms 𝑘 𝐿𝐴 , 𝑘 𝐵𝑢𝑂𝐻 , 𝑘 𝐵𝐿 and 𝑘 𝑊 are the specific constants of the enzyme for levulinic acid, butanol, butyl levulinate and water, respectively. 𝐾 𝐵𝑢𝑂𝐻 and 𝐾 𝑊 are Michaelistype constants for butanol and water, respectively. 𝐾 𝐼𝐼𝐵𝑢𝑂𝐻 and 𝐾 𝐼𝐼𝑊 are the inhibition constants for butanol and water, respectively.

By considering the binding of BuOH or water with the free enzyme, D' becomes D''

𝐷 ′′ = 𝑘 𝐿𝐴 • [𝐿𝐴] • (1 + [𝑊] 𝐾 𝐼𝐼𝑊 + [𝐵𝑢𝑂𝐻] 𝐾 𝐵𝑢𝑂𝐻 ) + (𝑘 𝐵𝑢𝑂𝐻 • [𝐵𝑢𝑂𝐻] + 𝑘 𝑊 • [𝑊]) • (1 + [𝐵𝑢𝑂𝐻] 𝐾 𝑠𝑖𝐵𝑢𝑂𝐻 ) • (1 + [𝑊] 𝐾 𝑠𝑖𝑊 ) +𝑘 𝐵𝐿 • [𝐵𝐿] • (1 + [𝑊] 𝐾 𝑊 + [𝐵𝑢𝑂𝐻] 𝐾 𝐼𝐼𝐵𝑢𝑂𝐻 ) ( 12
)
where, 𝐾 𝑠𝑖𝐵𝑢𝑂𝐻 and 𝐾 𝑠𝑖𝑊 stand for the constants between BuOH and water with the free enzyme.

According to Mitchell and Krieger [START_REF] Mitchell | Looking through a new lens: Expressing the Ping Pong bi bi equation in terms of specificity constants[END_REF], these equations are mathematically symmetric and less-lumped parameters. Different models were evaluated based on reaction rates developed by Mitchell and Krieger [START_REF] Mitchell | Looking through a new lens: Expressing the Ping Pong bi bi equation in terms of specificity constants[END_REF] and Varma and Madras [START_REF] Varma | Kinetics of synthesis of butyl butyrate by esterification and transesterification in supercritical carbon dioxide[END_REF].

Kinetic modeling

Experiments were performed in isothermal conditions, and internal and external mass transfers were found to be negligible. Plug-flow model was used; thus, material balances for each species can be written as

𝑑[𝐵𝐿] 𝑑𝜏 = 𝑟 𝐸𝑠𝑡𝑒𝑟𝑖𝑓𝑖𝑐𝑎𝑡𝑖𝑜𝑛 (13) 𝑑[𝐵𝑢𝑂𝐻] 𝑑𝜏 = -𝑟 𝐸𝑠𝑡𝑒𝑟𝑖𝑓𝑖𝑐𝑎𝑡𝑖𝑜𝑛 (14) 𝑑[𝐿𝐴] 𝑑𝜏 = -𝑟 𝐸𝑠𝑡𝑒𝑟𝑖𝑓𝑖𝑐𝑎𝑡𝑖𝑜𝑛 (15) 𝑑[𝑊] 𝑑𝜏 = 𝑟 𝐸𝑠𝑡𝑒𝑟𝑖𝑓𝑖𝑐𝑎𝑡𝑖𝑜𝑛 ( 16
)
where, 𝜏 is the space time defined as 𝑉 𝐿

𝑄

, VL and Q are the volume of the liquid in the reactor and volumetric flow-rate, respectively.

Ordinary differential equations ODEs ( 13)-( 16) were solved by the solver DDALPUS algorithm, via a damped Newton method [START_REF] Caracotsios | Sensitivity analysis of initial value problems with mixed odes and algebraic equations[END_REF].

For the non-linear regression, the concentration of BL was used as an observable. The estimation of the different kinetic constants (Equations ( 8)-( 12)) was done via the minimization of the objective function 𝑆(𝜃) expressed as

𝑆(𝜃) = ∑ 𝑤 𝑢 • ([𝐵𝐿] 𝑒𝑥𝑝,𝑢 -[𝐵𝐿] 𝑠𝑖𝑚,𝑢 ) 2 𝑛 𝑢=1 = 𝑆𝑆𝑅 (17)
where, 𝑤 𝑢 is the weigh factor for the experimental value u.

The objective function is expanded as a quadratic function of the parameters around the initial parameter values of the current iteration. The resulting quadratic minimization problem is solved with a modified Gauss-Jordan algorithm within a user-defined feasible region; then, a weak line search is conducted to establish an improved objective value and initial parameter vector for the next iteration. Interval estimates for the individual estimated parameters are then calculated from the final quadratic expansion of the objective function. This minimization is done by the package GREGPLUS to provide optimal parameter estimates with the 95% confidence intervals, expressed by the highest probability density (HPD). GREGPLUS provides the normalized parameter covariance matrix.

The GREGPLUS package and DDAPLUS solver are implemented in the Athena Visual Studio® 14.2 [START_REF] Stewart | Computer-Aided Modeling of Reactive Systems[END_REF] used in this study.

Different models were evaluated based on the Ping-Pong Bi-Bi mechanism developed by the Varma and Madras study [START_REF] Varma | Kinetics of synthesis of butyl butyrate by esterification and transesterification in supercritical carbon dioxide[END_REF] and the Mitchell and Krieger study [START_REF] Mitchell | Looking through a new lens: Expressing the Ping Pong bi bi equation in terms of specificity constants[END_REF]. The term

[𝐸] 0 was expressed by the catalyst loading ρEnzyme, i.e., the mass of catalyst divided by the volume of liquid in the reactor.

The general equation for esterification can be derived as

𝑟 𝐸𝑠𝑡𝑒𝑟𝑖𝑓𝑖𝑐𝑎𝑡𝑖𝑜𝑛 = 𝑘 𝐸𝑠𝑡𝑒𝑟𝑖𝑓𝑖𝑐𝑎𝑡𝑖𝑜𝑛 • 𝜌 𝐸𝑛𝑧𝑦𝑚𝑒 • 1 𝐷 • ([𝐿𝐴] • [𝐵𝑢𝑂𝐻] - [𝐵𝐿] • [𝑊] 𝐾 𝑒𝑞 ) ( 18
)
Different models were assessed, as summarized in Table 4.

Model 1 is the simplest one, by letting the denominator D equal to 1.

Models 2-4 are derived from Varma and Madras. Equation ( 8) was divided by 𝑘 𝑏 .To ease the parameter estimation and avoid division by very low number (close to zero) or high number, the following modification were included

+𝐾 𝐿𝐴 • [𝐵𝑢𝑂𝐻] +𝐾 𝑊 ′′ • [𝐵𝐿] +𝐾 𝐵𝐿 ′′ • [𝑊] +𝐾 𝐿𝑢𝑚𝑝 • [𝑊] • [𝐵𝐿] Model 3 kf . ρEnzyme [𝐿𝐴] • [𝐵𝑢𝑂𝐻] + 𝐾 𝐵𝑢𝑂𝐻 • [𝐿𝐴] +𝐾 𝐿𝐴 • [𝐵𝑢𝑂𝐻] + 𝐾 𝐼𝐵𝑢𝑂𝐻 ′′ • [𝐵𝑢𝑂𝐻] 2 +𝐾 𝑊 ′′ • [𝐵𝐿] +𝐾 𝐵𝐿 ′′ • [𝑊] +𝐾 𝐿𝑢𝑚𝑝 • [𝑊] • [𝐵𝐿] Model 4 kf . ρEnzyme [𝐿𝐴] • [𝐵𝑢𝑂𝐻] + 𝐾 𝐵𝑢𝑂𝐻 • [𝐿𝐴] + 𝐾 𝐼𝐿𝐴 ′′ • [𝐿𝐴] 2 +𝐾 𝐿𝐴 • [𝐵𝑢𝑂𝐻] + 𝐾 𝐼𝐵𝑢𝑂𝐻 ′′ • [𝐵𝑢𝑂𝐻] 2 +𝐾 𝑊 ′′ • [𝐵𝐿] + 𝐾 𝐼𝐵𝐿 ′ • [𝐵𝐿] 2 +𝐾 𝐵𝐿 ′′ • [𝑊] + 𝐾 𝐼𝑊 ′′ • [𝑊] 2 +𝐾 𝐿𝑢𝑚𝑝 • [𝑊] • [𝐵𝐿] + 𝐾 𝐼𝐼𝐿𝐴 ′′′ • [𝑊] • [𝐿𝐴] +𝐾 𝐼𝐼𝐵𝐿 ′′′ • [𝐿𝐴] • [𝐵𝐿] +𝐾 𝐵𝑢𝑂𝐻-𝑊 ′′ • [𝑊] • [𝐵𝑢𝑂𝐻] +𝐾 𝐿𝐴-𝐵𝐿 ′′ • [𝐿𝐴] • [𝐵𝐿] Model 5 kLA . ρEnzyme 𝐾 1 • [𝐿𝐴] + 𝐾 𝐼𝑊 ′′ • [𝑊] • [𝐿𝐴] + 𝐾 𝐵𝑢𝑂𝐻 ′′ • [𝐵𝑢𝑂𝐻] • [𝐿𝐴] +[𝐵𝑢𝑂𝐻] + 𝐾 2 • [𝑊] +𝐾 3 • [𝐵𝐿] + 𝐾 𝑊 ′ • [𝑊] • [𝐵𝐿] + 𝐾 𝐼𝐵𝑢𝑂𝐻 ′ • [𝐵𝑢𝑂𝐻] • [𝐵𝐿] Model 6 kLA . ρEnzyme 𝐾 1 • [𝐿𝐴] + 𝐾 𝐼𝑊 ′′ • [𝑊] • [𝐿𝐴] + 𝐾 𝐵𝑢𝑂𝐻 ′′ • [𝐵𝑢𝑂𝐻] • [𝐿𝐴] +([𝐵𝑢𝑂𝐻] + 𝐾 2 • [𝑊]) • (1 + 𝐾 𝑠𝑖𝐵𝑢𝑂𝐻 ′′ • [𝐵𝑢𝑂𝐻]) • (1 + 𝐾 𝑠𝑖𝑊 ′′ • [𝑊]) +𝐾 3 • [𝐵𝐿] + 𝐾 𝑊 ′ • [𝑊] • [𝐵𝐿] + 𝐾 𝐼𝐵𝑢𝑂𝐻 ′ • [𝐵𝑢𝑂𝐻] • [𝐵𝐿]
To decrease the correlation between the pre-exponential factor and activation energy and ease the parameter estimation stage, the following modified Arrhenius equation was used [START_REF] Buzzi-Ferraris | Planning of experiments and kinetic analysis[END_REF].

𝑘 𝑐 (𝑇) = 𝑒𝑥𝑝 [𝑙𝑛 (𝑘 𝑐 (𝑇 𝑟𝑒𝑓 )) + 𝐸 𝑎 𝑅 • 𝑇 𝑟𝑒𝑓 • (1 - 𝑇 𝑟𝑒𝑓 𝑇 )] (20)
where, Tref is a reference temperature which is the average temperature of the experimental matrix (Table 1).

The following constants were estimated: 𝑙𝑛 (𝑘 𝑐 (𝑇 𝑟𝑒𝑓 )), 𝐸 𝑎 𝑅•𝑇 𝑟𝑒𝑓 , Michaelis-Menten and inhibition constants. Michaelis-Menten and inhibition constants were assumed to be temperature independent.

The effect of the number of estimated parameters on the models was evaluated via the AIC number standing for Akaike Information Criterion [START_REF] Capecci | Bayesian Statistics to Elucidate the Kinetics of γ-Valerolactone from n-Butyl Levulinate Hydrogenation over Ru/C[END_REF][START_REF] Delgado | Kinetic model assessment for the synthesis of γ-valerolactone from n-butyl levulinate and levulinic acid hydrogenation over the synergy effect of dual catalysts Ru/C and Amberlite IR-120[END_REF][START_REF] Mcdonald | Kinetic model discrimination of penicillin G acylase thermal deactivation by non-isothermal continuous activity assay[END_REF].

𝐴𝐼𝐶 = 𝑁𝑢𝑚𝑏𝑒𝑟 𝑜𝑓 𝑖𝑛𝑑𝑒𝑝𝑒𝑛𝑑𝑎𝑛𝑡 𝑒𝑣𝑒𝑛𝑡 • 𝑙𝑛 ( 𝑆𝑆𝑅 𝑛𝑢𝑚𝑏𝑒𝑟 𝑜𝑓 𝑖𝑛𝑑𝑒𝑝𝑒𝑛𝑑𝑎𝑛𝑡 𝑒𝑣𝑒𝑛𝑡 ) +2 • 𝑁𝑢𝑚𝑏𝑒𝑟 𝑜𝑓 𝑒𝑠𝑡𝑖𝑚𝑎𝑡𝑒𝑑 𝑝𝑎𝑟𝑎𝑚𝑒𝑡𝑒𝑟𝑠 (21)

Modeling results

In the first step, preliminary modeling results showed that some parameters tend to zero. Thus, these parameters were discarded:

-For Model 2, 𝐾 𝑊 ′′ , 𝐾 𝐿𝑢𝑚𝑝 and 𝐾 𝐵𝐿 ′′ were discarded in the modeling.

-For Model 3, 𝐾 𝐿𝐴 , 𝐾 𝑊 ′′ and 𝐾 𝐵𝐿 ′′ were not considered.

- By discarding these parameters, the reduced models are displayed in Table 5. Table 6 is a summary of the modeling output for the different models. SSR is the sum of squared residuals, the difference between the experimental and simulated concentrations. AIC values showed that Model 4 is the most probable one (Table 4).

Due to space limitation of the journal, the modeling results of the other models are displayed in Supporting Information (S5). However, the fit of Model 4 to experimental concentration for experiments carried out with a molar ratio LA/BuOH: 1:1 is lower near to the equilibrium than for the other ratio.

This observation is because the equilibrium constant predicted by ePC-SAFT is less reliable for this ratio.

Conclusions

The synthesis of butyl levulinate from the esterification of levulinic acid over Novozym®435, an immobilized enzyme, was investigated in microfluidic technology in isothermal conditions. The enzyme's catalytic activity was found to be stable for 400 minutes. The internal and external mass transfer resistance was found to be negligible.

Thus, a plug-flow model was used to estimate the kinetic constants.

Several kinetic experiments were performed by varying the reaction temperature from 20 to 80 °C, inlet LA concentration from 1.74 to 5.09 mol.L -1 , inlet butanol concentration from 5.09 to 9.07 mol.L -1 , mass of dried enzyme from 50 to 150 mg and residence time from 1 to 30 minutes.

The equilibrium constants were evaluated via the ePC-SAFT equation of state.

We evaluated 6 kinetic models based on power law, the classical Ping-Pong Bi-Bi

Fig. 2 .

 2 Fig. 2. Evolution of experimental BL concentration (blue circle) with error bars (black)

Fig. 3 .

 3 Fig.3. Effect of temperature on the experimental concentration ratio [BL]/[LA]inlet

Figs 4

 4 Figs 4 and 5 show the effect of this ratio on the kinetics, via the normalized ratio [𝐵𝐿][𝐿𝐴] 𝑖𝑛𝑙𝑒𝑡 .

Fig. 4 . 1 Fig. 5 . 35 °C 1 (

 415351 Fig.4. Effect of [𝐵𝑢𝑂𝐻] 𝑖𝑛𝑙𝑒𝑡 [𝐿𝐴] 𝑖𝑛𝑙𝑒𝑡 ratio on experimental concentration ratio [BL]/[LA]inlet at 50°C

Fig. 7

 7 Fig. 7 shows that the residuals, [𝐵𝐿] 𝑒𝑥𝑝,𝑢 -[𝐵𝐿] 𝑠𝑖𝑚,𝑢 , are randomly distributed versus

Fig. 8

 8 Fig.8shows the parity plot, and one can notice that Model 4 can predict the

Fig. 7 .

 7 Fig. 7. Residual plots for Model 4.

Fig. 8 .

 8 Fig. 8. Parity plot for Model 4.

Fig. 9 .

 9 Fig. 9. Fit of Model 4 to the experimental concentrations with 95% prediction

  mechanism and the modified one developed by Mitchell and Krieger. Based on the Akaike information criterion, we found that the classical Ping-Pong model, including the inhibition mechanism by butanol and levulinic acid, can fit the experimental concentrations properly. This model can reasonably predict the experimental concentration by considering the temperature effect on the rate constant.

Table 1 .

 1 Different flow-rates used for the kinetic monitoring.

	Mass of Novozym®435
	(mg)
		50 mg	100 mg	150 mg
	30	6.5 µL/min	13 µL/min	17.33 µL/min
	20	9.75 µL/min	19.5 µL/min	26 µL/min
	15	13 µL/min	26 µL/min	34.66 µL/min
	10	19.5 µL/min	39 µL/min	52 µL/min
	7	27.85 µL/min	55.71 µL/min	74.28 µL/min
	5	39 µL/min	78 µL/min	104 µL/min
	3	65 µL/min	130 µL/min	173.33 µL/min
	2	97.5 µL/min	195 µL/min	260 µL/min
	1	195 µL/min	390 µL/min	520 µL/min
	0	0 µL/min	0 µL/min	0 µL/min

Table 2

 2 shows the experimental matrix used in this study.

	Residence Time
	(min)

Table 2 .

 2 Experimental matrix for esterification.

	123

Inlet (mol/L) Run Temperature (°C) Mass of catalyst (mg) Void volume (μL) BL LA BuOH W Ratio [BuOH]in/[LA]in

	2.4 Stability of Novozym 435
	It is fundamental to evaluate the Novozym 435 stability to know if there are enzymes
	or another product leaching from support denaturation. Two experiments were
	1	20	100	390	0.00 2.61	7.81	0.34	3.00
	2 performed: a degradability test for Novozym 435 in butanol solvent in a batch reactor 35 100 390 0.00 2.61 7.81 0.34 3.00
	3	50	100	390	0.00 2.61	7.81	0.34	3.00
	4 and deactivation in the microfluidic system. 65 100 390	0.00 2.61	7.82	0.34	3.00
	5	80	100	390	0.00 2.61	7.81	0.34	3.00
	6 The degradability test was performed as follows: 100 mg of Novozym 435 were placed 65 150 520 0.00 2.61 7.81 0.34 2.99
	7	65	50	195	0.00 2.61	7.81	0.34	2.99
	8 into 390 µL of a solution containing a [BuOH]inlet/[LA]inlet ratio = 5. The mixture was 65 100 390 0.00 5.09 5.09 0.31 1.00
	9	50	100	390	0.00 5.16	5.16	0.32	1.00
	10 heated at 65°C. After 1 h, the sample was filtered through a 30 µm PTFE frit of the 35 100 390 0.00 5.09 5.09 0.31 1.00
	11	20	100	390	0.00 5.09	5.09	0.31	1.00
	12 Omnifit column (the exact same one that we used for the kinetic study), then CDCl3 65 150 520 0.00 5.09 5.09 0.31 1.00
	13	65	50	195	0.00 5.09	5.09	0.31	1.00
	14	65	100	390	0.00 1.74	8.71	0.32	4.99
	15	50	100	390	0.00 1.77	8.83	0.33	5.00
	16	35	100	390	0.00 1.79	8.95	0.33	4.99
	17	20	100	390	0.00 1.82	9.07	0.34	4.99
	18	65	50	195	0.00 1.74	8.71	0.32	4.99
	19	65	150	520	0.00 1.74	8.71	0.32	4.99
	20	65	100	390	0.00 2.59	7.76	0.70	3.00
	21	50	100	390	0.00 2.59	7.76	0.70	3.00
	22	35	100	390	0.00 2.59	7.76	0.70	3.00
	23	20	100	390	0.00 2.59	7.76	0.70	3.00
	24	65	100	390	0.00 2.46	7.78	1.20	3.17
	25	50	100	390	0.00 2.57	7.69	1.20	3.00
	26	35	100	390	0.00 2.46	7.78	1.20	3.17
	27	20	100	390	0.00 2.46	7.78	1.20	3.17
	124
	125

Table 3 .

 3 Calculated equilibrium constant 𝐾 𝑡ℎ based on the equilibrium concentrations of Runs 1-5 (Table2).

	T / °C	Keq	Kγ	Kth
	20	0.157 3.578	0.56
	35	0.309 2.793	0.67 *
	50	0.356 2.308	0.82
	65	0.491 1.969	0.97
	80	0.695 1.732	1.20
		*interpolated value

  • [𝐸] 0 and 𝑘 𝑓 • [𝐸] 0 represent kinetic rate constants. The terms 𝐾 𝐿𝐴 , 𝐾 𝐵𝐿 , 𝐾 𝐵𝑢𝑂𝐻 and 𝐾 𝑊 are the Michaelis constants for LA, BL, BuOH, and W, respectively. The inhibition constants by levulinic acid and BL are defined by 𝐾 𝐼𝐼𝐿𝐴 and 𝐾

𝐼𝐼𝐵𝐿 . The dissociation constants, representing the dissociation of the inhibitor from the corresponding enzyme-inhibitor, are 𝐾 𝐼𝐿𝐴 , 𝐾 𝐼𝐵𝐿 , 𝐾 𝐼𝐵𝑢𝑂𝐻 and 𝐾 𝐼𝑊 . The adsorption constants 𝐾 𝐵𝑢𝑂𝐻-𝑊 and 𝐾 𝐿𝐴-𝐵𝐿 are lumped constants. Mitchell and Krieger proposed a new rate expression of this Ping-Pong Bi-Bi mechanism [63],

Table 4 .

 4 Kinetic models tested in this study.

	𝐾 𝐼𝐿𝐴 ′′ =	𝐾 𝐵𝑢𝑂𝐻 𝐾 𝐼𝐿𝐴	, 𝐾 𝐼𝐵𝑢𝑂𝐻 ′′	=	𝐾 𝐿𝐴 𝐾 𝐼𝐵𝑢𝑂𝐻	, 𝐾 𝑊 ′′ =	𝐾 𝑊 𝐾 𝑒𝑞	, 𝐾 𝐵𝐿 ′′ =	𝐾 𝐵𝐿 𝐾 𝑒𝑞	, 𝐾 𝐼𝐵𝐿 ′′ =	𝐾 𝑊 ′′ 𝐾 𝐼𝐵𝐿	, 𝐾 𝑊 ′ = 𝑘 𝑓 • 𝐾 𝑊 •	1 𝑘 𝑏	•	1 𝐾 𝑒𝑞	,
										𝐾 𝐼𝑊 ′′ =	′′ 𝐾 𝐵𝐿 𝐾 𝐼𝑊	,
								𝐾 𝐵𝐿 ′ = 𝑘 𝑓 • 𝐾 𝐵𝐿 •	1 𝑘 𝑏	•	1 𝐾 𝑒𝑞
		𝐾 𝐿𝑢𝑚𝑝 =	𝑘 𝑓 𝑘 𝑏 • 𝐾 𝑒𝑞	, 𝐾 𝐼𝐼𝐵𝐿 ′′	=	𝐾 𝐿𝐴 𝐾 𝐼𝐼𝐵𝐿	, 𝐾 𝐵𝑢𝑂𝐻-𝑊 ′′	=	𝐾 𝐿𝐴 𝐾 𝐵𝑢𝑂𝐻-𝑊	𝑎𝑛𝑑 𝐾 𝐿𝐴-𝐵𝐿 ′′	=	𝐾 𝐵𝑢𝑂𝐻 𝐾 𝐿𝐴-𝐵𝐿
	Model 2 ignores the inhibition mechanism, hence 𝐾 𝐼𝐿𝐴 ′′ , 𝐾 𝐼𝐵𝐿 ′′ , 𝐾 𝐼𝑊 ′′ , 𝐾 𝐼𝐵𝑢𝑂𝐻 ′′	, 𝐾 𝐼𝐼𝐿𝐴 ′′ , 𝐾 𝐼𝐼𝐵𝐿 ′′ ,
	𝐾 𝐿𝐴-𝐵𝐿 ′′	, and 𝐾 𝐵𝑢𝑂𝐻-𝑊 ′′	were fixed to zero.
	Model 3 considers the inhibition by butanol and ignores the other inhibition mechanism,
	hence 𝐾 𝐼𝐿𝐴 ′′ , 𝐾 𝐼𝐵𝐿 ′′ , 𝐾 𝐼𝑊 ′′ , 𝐾 𝐼𝐼𝐿𝐴 ′′ , 𝐾 𝐼𝐼𝐵𝐿 ′′ , 𝐾 𝐿𝐴-𝐵𝐿 ′′	, and 𝐾 𝐵𝑢𝑂𝐻-𝑊 ′′	were fixed to zero.
	Model 4 considers all inhibition mechanisms.
	Models 5-6 are derived from Mitchell and Krieger.

  For Model 4, 𝐾 𝐵𝑢𝑂𝐻 , 𝐾 𝐿𝐴 , 𝐾 𝑊 ′′ , 𝐾 𝐵𝐿 ′′ , 𝐾 𝐿𝑢𝑚𝑝 , 𝐾 𝐼𝐼𝐿𝐴 ′′′ , 𝐾 𝐼𝐵𝐿 ′ , 𝐾 𝐼𝑊 ′′ and 𝐾 𝐼𝐵𝐿 ′′′ and 𝐾 𝐿𝐴-𝐵𝐿

				′′	were
	neglected.
	-For Model 5, 𝐾 𝑊 ′′ , 𝐾 𝐼𝑊 ′′ , 𝐾 1 , 𝐾 2 , 𝐾 𝐼𝐵𝑢𝑂𝐻 ′	and 𝐾 3 were neglected.
	-For Model 6, 𝐾 𝐼𝐵𝑢𝑂𝐻 ′′	, 𝐾 𝐼𝑊 ′′ , 𝐾 𝐵𝑢𝑂𝐻 ′′	, 𝐾 2 , 𝐾 𝑊 ′ and 𝐾 3 were neglected.

Table 5 .

 5 Reduced kinetic models for the esterification of levulinic acid over immobilized enzyme.

	Model	Kinetic term		Denominator
	Model 1	kEsterification . ρEnzyme			1
	Model 2	kf . ρEnzyme	[𝐿𝐴] • [𝐵𝑢𝑂𝐻] + 𝐾 𝐵𝑢𝑂𝐻 • [𝐿𝐴]
				+𝐾 𝐿𝐴 • [𝐵𝑢𝑂𝐻]
	Model 3	kf . ρEnzyme	[𝐿𝐴] • [𝐵𝑢𝑂𝐻] + 𝐾 𝐵𝑢𝑂𝐻 • [𝐿𝐴]
			+𝐾 𝐼𝐵𝑢𝑂𝐻 ′′	• [𝐵𝑢𝑂𝐻] 2 + 𝐾 𝐿𝑢𝑚𝑝 • [𝑊] • [𝐵𝐿]
	Model 4	kf . ρEnzyme	[𝐿𝐴] • [𝐵𝑢𝑂𝐻] + 𝐾 𝐼𝐿𝐴 ′′ • [𝐿𝐴] 2
				+𝐾 𝐼𝐵𝑢𝑂𝐻 ′′	• [𝐵𝑢𝑂𝐻] 2
			+𝐾 𝐵𝑢𝑂𝐻-𝑊 ′′	• [𝑊] • [𝐵𝑢𝑂𝐻]
	Model 5	kLA . ρEnzyme		𝐾 𝐵𝑢𝑂𝐻 ′′	• [𝐵𝑢𝑂𝐻] • [𝐿𝐴]
				+[𝐵𝑢𝑂𝐻]
	Model 6	kLA . ρEnzyme	𝐾 1 • [𝐿𝐴] + ([𝐵𝑢𝑂𝐻]) • (1 + 𝐾 𝑠𝑖𝐵𝑢𝑂𝐻 ′′	• [𝐵𝑢𝑂𝐻])
			• (1 + 𝐾 𝑠𝑖𝑊 ′′ • [𝑊])

Table 6 .

 6 Modeling results for each Model.

		Model 1 Model 2 Model 3 Model 4 Model 5	Model 6
	SSR	15.85	15.16	14.83	13.66	15.48	14.03
	Number of estimated parameters	2	4	5	5	3	5
	AIC	-2311.47 -2335.50 -2347.36 -2398.81 -2324.11 -2381.86

Table 7

 7 presents the Normalized parameter covariance matrix for Model 4. According to Toch et al.[START_REF] Toch | A systematic methodology for kinetic modeling of chemical reactions applied to n-hexane hydroisomerization[END_REF], two parameters are correlated if their binary correlation coefficient is higher than 0.95. From Table8, one can notice that the estimated parameters are not correlated.

Table 7 .

 7 Estimated values at Tref = 51°C and statistical data for Model 4.

		2.00
		1.80
		1.60
	(mol.L -1 )	1.40
	CBL_sim
		0.00	0.50	1.00	1.50	2.00	2.50
				CBL_exp (mol.L -1 )

Table 8 .

 8 Normalized parameter covariance matrix for Model 4.

		ln (𝑘 𝑓 (𝑇 𝑟𝑒𝑓 ))	𝐸𝑎 𝑓 𝑅 • 𝑇 𝑟𝑒𝑓	𝐾 𝐼𝐿𝐴 ′′	𝐾 𝐵𝑢𝑂𝐻-𝑊 ′′
	ln (𝑘 𝑓 (𝑇 𝑟𝑒𝑓 ))	1
	𝐸𝑎 𝑓 𝑅 • 𝑇 𝑟𝑒𝑓	0.1	1
	𝐾 𝐼𝐿𝐴 ′′	0.69	0.09	1
	𝐾 𝐵𝑢𝑂𝐻-𝑊 ′′	0.93	0.12	0.53	1

Acknowledgments

The authors thank INSA Rouen Normandy, University of Rouen Normandy, the Centre National de la Recherche Scientifique (CNRS), European Regional Development Fund (ERDF) N• HN0001343, Labex SynOrg (ANR-11-LABX-0029), Carnot Institute I2C, the graduate school for reasearch XL-Chem (ANR-18-EURE-0020 XL CHEM) and Region Normandie for their support. This research was funded, in whole or in part, by the ANR (French National Research Agency) and the DFG (German Research Foundation) through the project MUST (MicroflUidics for Structure-reactivity relationships aided by Thermodynamics & kinetics) [ANR-20-CE92-0002 & Project number 446436621].

Declaration of Competing Interest

The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.