René Schott
email: schott@loria.fr

G Stacey Staples
email: sstaple@siue.edu

Operator calculus in generalized zeon algebras: theory and application to multi-constrained path problems

Keywords: shortest paths, message routing, operator calculus, semigroup algebras MSC: 68W30, 81R05

Classical approaches to routing problems invariably require construction of trees and the use of heuristics to prevent combinatorial explosion. The operator calculus approach presented herein, however, allows such explicit tree constructions to be avoided. Introduced here is the notion of generalized zeon algebras and their associated operator calculus. The inherent combinatorial properties of generalized zeons make them useful for routing problems by implicitly pruning the underlying tree structures. As an application, an operator calculus approach to multi-constrained path problems is described.

Introduction

Operator calculus (OC) methods on graphs have been developed in a number of earlier works by Schott and Staples [START_REF] Schott | Nilpotent adjacency matrices and random graphs[END_REF][START_REF] Schott | Nilpotent adjacency matrices, random graphs, and quantum random variables[END_REF][START_REF] Staples | Clifford-algebraic random walks on the hypercube[END_REF][START_REF] Staples | A new adjacency matrix for finite graphs[END_REF]. The principal idea underlying the approach is the association of graphs with algebraic structures whose properties reveal information about the associated graphs. By constructing the "nilpotent adjacency matrix" associated with a finite graph, information about self-avoiding structures (paths, cycles, trails, etc.) in the graph are revealed by computing powers of the matrix.

In the operator calculus approach, graded semigroup algebras are generated by "null-square" elements such that properties of the algebra "sieve out" paths. In other words, cycles are removed from consideration automatically.

Consider a directed graph G = (V, E) on n vertices such that associated with each edge (v i , v j ) ∈ E is a vector weight w ij = (w ij1 , . . . , w ijm ) ∈ R m .

The point w * = (w * 1 , . . . , w * m ) ∈ X ⊂ R m is referred to as a Pareto minimum of X if there does not exist w = (w 1 , . . . , w m ) ∈ X such that (∀i) [w i ≤ w * i ], and (1.1) (∃j) [w j < w * j ].

(1.2)

Equivalently, one says that w * is nondominated from below.

Defining the weight of a path in an edge-weighted graph as the sum of vector weights of arcs contained in the path, a Pareto path is then a path whose weight is a Pareto minimum.

For the case m = 1, Dijkstra's algorithm finds all single source minimum paths in a directed graph on n vertices with nonnegative edge weights in O(n 2 ) time [START_REF] Dijkstra | A note on two problems in connexion with graphs[END_REF]. The Bellman-Ford algorithm finds single source minimal paths in digraphs with arbitrary edge weights and runs in O(n |E|) time [START_REF] Bellman | On a routing problem[END_REF][START_REF] Ford | Flows in Networks[END_REF].

In the more general case m ≥ 1, Corley and Moon [START_REF] Corley | Shortest paths in networks with vector weights[END_REF] presented an algorithm for finding all Pareto paths with computational complexity O(mn 2n-3 + mn n ).

The aim of the current work is to find Pareto paths satisfying multiple constraints. Given a vector c = (c 1 , . . . , c m ) ∈ R m , a path is deemed feasible if its vector weight w = (w 1 , . . . , w m ) satisfies

(∀i) [w i ≤ c i ].
Letting P f denote the collection of feasible paths having fixed source v 0 and fixed target v ∞ , the goal is to find a path in P f whose weight is a Pareto minimum. The operator calculus approach described herein can be applied to sieve out the collection of feasible paths and recover all single-source Pareto paths remaining.

Theory: Generalized zeon algebras

Zeon algebras are commutative algebras generated by collections of null-squares, {ζ i : 1 ≤ i ≤ n} with ζ i 2 = 0 for each i. Their combinatorial properties make them useful for a variety of counting properties, as seen in a number of previous works by the current authors. By choosing sufficiently large sets of generators, they can be generalized to algebras whose generators are nilpotent of arbitrary index. The resulting generalized zeon algebras are suitable for a number of combinatorial applications, including multi-constrained routing problems.

Definition 2.1. The n-particle zeon algebra, denoted by Cℓ n nil , is defined as the real abelian algebra generated by the collection {ζ i } (1 ≤ i ≤ n) along with the scalar 1 = ζ 0 subject to the following multiplication rules:

ζ i ζ j = ζ j ζ i for i = j, and
(2.1)

ζ i 2 = 0 for 1 ≤ i ≤ n. (2.2)
A general element u ∈ Cℓ n nil can be expanded as

u = I∈2 [n] u I ζ I , (2.3)
where I ∈ 2 [n] is a subset of [n] = {1, 2, . . . , n} used as a multi-index, u I ∈ R, and

ζ I = ι∈I ζ ι .
Remark 2.2. The zeon algebra Cℓ n nil can be realized as a commutative subalgebra of the Grassmann algebra V over a 2n-dimensional vector space V with orthonormal basis {γ i } by defining

ζ i = γ i γ n+i for each 1 ≤ i ≤ n.
A canonical basis element ζ I is referred to as a blade. The number of elements in the multi-index I is referred to as the grade of the blade ζ I .

The next lemma shows that it is possible to construct elements with arbitrary index of nilpotency within a zeon algebra of sufficiently high dimension.

Lemma 2.3. Let {ζ i : 1 ≤ i ≤ n} be the null-square generators of Cℓ n nil . Then, for any permutation σ ∈ S n and positive integers ℓ

≤ k ≤ n,   k j=1 ζ σ(j)   ℓ = ℓ! I⊆{σ(1),...,σ(k)} |I|=ℓ ζ I . (2.4) Moreover, if ℓ > k, then   k j=1 ζ σ(j)   ℓ = 0. (2.5)
Proof. Since the generators commute, the multinomial theorem applies, with only square-free terms surviving.

Definition 2.4. For positive integer n, let s = (s 1 , . . . , s n ) ∈ N n be an n-tuple of positive integers. Then, the zeon algebra of signature s (or s-zeon algebra), denoted Cℓ s nil , is the real abelian algebra generated by the collection {ν i } (1 ≤ i ≤ n) along with the scalar 1 = ν 0 subject to the following multiplication rules:

ν i ν j = ν j ν i for i = j, and (2.6)
ν i si = 0 for 1 ≤ i ≤ n. (2.7)
For convenience, the following multi-exponent notation is adopted:

ν x := n i=1 ν xi i := ν x1 1 • • • ν xn n .
(2.8)

Letting S = {(x 1 , . . . , x n ) : 0 ≤ x i ≤ s i } ⊂ N 0 n , a general element u ∈ Cℓ s nil can be expanded as u = x∈S u x ν x , (2.9)
where u x ∈ R for each multi-exponent x.

Since the components of signature vectors and arbitrary multi-exponents are nonnegative integers, the 1-norm of such a vector x is simply the sum of the components; that is,

x 1 = n i=1 |x i | = n i=1 x i .
(2.10)

Lemma 2.5. The algebra Cℓ s nil is isomorphic to a subalgebra of the zeon algebra Cℓ s 1 nil .

Proof. For k ∈ {1, . . . , n}, let p(k) denote the k th partial sum

p(k) = k j=1 s j , (2.11)
and define p(0) = 0. In light of Lemma 2.3, the desired isomorphism Cℓ s nil → Cℓ s 1 nil is obtained from the mapping

ν i → si j=1
ζ j+p(i-1) .

(2.12)

Finally, note that the s-zeon algebra is naturally graded according to

Cℓ s nil = s 1 k=0 Cℓ s nil k , (2.13)
where the grade-k part of the algebra is defined by

Cℓ s nil k = span ({ν x : x 1 = k}) . (2.14)
The notation • k extends naturally to elements of Cℓ s nil .

Operator calculus in generalized zeon algebras

The motivation for development of s-zeon operator calculus is based on polynomial operator calculus. To begin, raising and lowering operators are defined naturally in terms of polynomial differentiation and integration operators on basis zeons regarded as polynomials in commuting variables. In this formulation, the generators {ν i } of Cℓ s nil are regarded as variables in the polynomial sense. For any generalized zeon algebra with n generators, let {e i : 1 ≤ i ≤ n} denote standard unit vectors of the form e i := (0, . . . , 1 i th pos. , . . . , 0). Arbitrary multi-exponents are then expressed in the form x = (x 1 , . . . , x n ) = n i=1

x i e i . Definition 3.1. Let s ∈ N 0 n be an arbitrary zeon signature. For 1 ≤ j ≤ n, define the j th s-zeon differentiation operator ∂/∂ν j on Cℓ s nil by linear extension of

∂ ∂ν j ν x = ν x-ej if x j ≥ 1, 0 otherwise. (3.1)
Definition 3.2. The s-zeon integrals are defined by

ν x dν j = ν x+ej if x j < s j -1, 0 otherwise. (3.2)
These polynomial operators induce combinatorial raising and lowering operators by which s-zeon monomials (blades) are raised from grade k to grade k + 1 or lowered from grade k to grade k -1. These raising and lowering operators can also be regarded as creation and annihilation operators in the sense of quantum mechanics.

Definition 3.3. For each 1 ≤ j ≤ n, define the j th raising operator R j by linear extension of

R j ν x = ν x dν j = ν x ν j . (3.3)
Define the j th lowering operator D j by linear extension of

D j ν x = ∂ ∂ν j ν x . (3.4)
Definition 3.4. The j th zeon number operator Λ j is defined on the generalized zeon algebra Cℓ s nil by linear extension of

Λ j (ν x ) := x j ν x . (3.5)
In particular, for arbitrary multi-exponents x, y and scalars α, β,

Λ j (αν x + βν y ) = αx j ν x + β y 1 ν y . (3.6)
Definition 3.5. The dual of the j th zeon number operator, denoted Λ j ⋆ , is defined on the generalized zeon algebra Cℓ s nil by linear extension of

Λ j ⋆ (ν x ) := (1/x j )ν x if x j > 0, 0 otherwise. (3.7)
For arbitrary multi-exponents x, y and scalars α, β,

Λ j ⋆ (αν x + βν y ) = α x j ν x + β y j ν y . (3.8)
An element u ∈ Cℓ s nil is said to be scalar-free if its canonical expansion is of the form

x =0 u x ν x .
(3.9)

Let Cℓ s nil⋆ denote the scalar-free subalgebra of Cℓ s nil⋆ ; that is,

Cℓ s nil⋆ := {u ∈ Cℓ s nil : u is scalar free}. (3.10)
The zeon occupancy operator Λ and its dual

Λ ⋆ are Λ = n j=1 Λ j and Λ ⋆ = n j=1 Λ j ⋆ , respectively.
Lemma 3.6. Letting Cℓ s nil⋆ denote the scalar-free subspace of Cℓ s nil⋆ ,

ΛΛ ⋆ Cℓs nil⋆ = Λ ⋆ Λ Cℓs nil⋆ = I. (3.11)
More specifically, ΛΛ ⋆ : Cℓ s nil → Cℓ s nil⋆ is an orthogonal projection.

Proof. Since the components of multi-exponents are nonnegative integers, writing x = (x 1 , . . . x n ) leads to the component sum as the 1-norm; i.e.,

x 1 = n j=1 x j . Let u = x =0 u x ν x ∈ Cℓ s nil and consider Λ (Λ ⋆ u) = Λ   x =0 u x x 1 ν x   = x =0 u x x 1 x 1 ν x = u.
(3.12)

A nearly identical argument shows Λ ⋆ (Λ u) = u. From the definitions of Λ and Λ ⋆ , it is apparent that for nonzero scalar α, Λ(αν x ) = 0 if and only if x = 0; i.e., αν x = α. The same can be said of Λ ⋆ ; i.e., ker Λ = ker Λ ⋆ = R. As a result, ker ΛΛ ⋆ = ker Λ ⋆ Λ = R.

Given an m-dimensional constraint vector s, a total ordering is induced on the set of multi-exponents x ∈ N 0 m by defining x y if and only if ∃k ≥ 1 such that x i ≤ y i ∀i ≤ k. When any such total ordering is assigned to the multi-exponents, one is able to define minimal elements of Cℓ s nil . This will be useful in subsequent applications in which minimal elements will be associated with optimal solutions. Definition 3.7. Fixing a total ordering of the multi-exponents, define a minimal term of u ∈ Cℓ s nil by

℧u := u x ′ ν x ′ , (3.13)
where x ′ x for all nonzero multi-exponents in the canonical expansion of u.

The s-zeon algebra will be applied in later sections to sieve out paths satisfying multiple constraints. In order to retain identifying information about the paths themselves, another generalization of zeon algebras is considered.

The path algebra RΩ n

For fixed positive integer n, consider the alphabet Σ n := {ω i : 1 ≤ i ≤ n}. For convenience, we adopt the following ordered multi-index notation. In particular, letting u = (u 1 , . . . , u k ) for some k, the notation ω u will be used to denote a sequence (or word) of distinct symbols of the form

ω u := ω u1 ω u2 • • • ω u k . (3.14)
Appending 0 to the set Σ n , multiplication is defined on the words constructed from elements of Σ n by

ω u ω v = ω u .v if u ∩ v = ∅, 0 otherwise, (3.15)
where u .v denotes sequence concatenation. One thereby obtains the noncommutative semigroup Ω n , whose elements are the symbol 0 along with all finite words on distinct generators (i.e., finite sequences of distinct symbols from the alphabet Σ n ). Since there are only n generators, it is clear that the maximum multi-index size of semigroup elements is n. Moreover, these symbols can appear in any order so that the order of the semigroup is

n k=0 n k k! = n k=0 (n) k .
Defining (vector) addition and real scalar multiplication on the semigroup yields the semigroup algebra RΩ n of dimension |Ω n |. This semigroup algebra will be referred to as a path algebra.

Consider the collection of ordered pairs P = {(ω i , ω j ) : i = j ∈ Ω n × Ω n , and note that |P | = n 2n. Imposing an ordering on P , a bijection f :

P → [n 2 -n] is obtained. Any k-subset of [n 2 -n] thereby determines a unique finite word of Ω n : ω u = ω u1 • • • ω u k+1 ↔ {f ((u 1 , u 2 )), . . . , f ((u k , u k+1 ))}. (3.16)
In this way, one obtains a one-to-one semigroup homomorphism φ : Cℓ |P | nil → Ω n onto the words of length two or more in Ω n . Hence, the semigroup algebra RΩ n can be regarded as an extension of a zeon algebra.

Remark 3.8. When the pairs of P are unordered, each k-subset of [n(n -1)/2] determines two finite words of Ω n : ω u , and its reversion ω u ; i.e.,

ω u = ω u1 • • • ω u k+1 ↔ {f ((u 1 , u 2 )), . . . , f ((u k , u k+1 ))} ↔ ω u k+1 • • • ω u1 = ω u .
(3.17)

Operator calculus on graphs

As discussed in previous work, the nilpotent adjacency matrix of a finite graph can be used to sieve out the graph's paths and cycles. The entries of this matrix are the null-square generators of a zeon algebra of appropriate dimension for the graph. The null-square properties of the algebra naturally remove entries corresponding to self-intersecting walks from powers of the nilpotent adjacency matrix.

Of interest in the current work is a method of sieving out paths with multidimensional weights (or costs) simultaneously satisfying a number of constraints. Generalized s-zeon generators will be associated with the graph's edges in such a way that paths whose weights exceed the constraints are zeroed out by the algebra's nilpotent properties.

In particular, extending the nilpotent adjacency matrix construction to Cℓ s nil ⊗ RΩ n allows one to enumerate (list) all paths and cycles satisfying multiple constraints in a finite graph by considering powers of the matrix. The associated tree structure underlying the cycle/path enumeration problem is automatically "pruned" by the inherent properties of the algebra.

The first step is defining a nilpotent adjacency matrix that preserves pathidentifying information. Definition 3.9. Let G = (V, E) be a graph on n vertices, either simple or directed with no multiple edges. Let {ω i }, 1 ≤ i ≤ n denote the null-square, noncommutative generators of RΩ n . Define the path-identifying nilpotent adjacency matrix A associated with G as the n × n matrix

A ij = ω j if (v i , v j ) ∈ E, 0 otherwise. (3.18)
Recalling Dirac notation, the i th row of A is conveniently denoted by v i | A while the j th column is denoted by A |v j . In this way, A is completely determined by

v i |A|v j = ω j if there exists a directed edge v i → v j in G, 0 otherwise, (3.19)
for all vertex pairs (v i , v j ) ∈ E.

Theorem 3.10. Let A be the path-identifying nilpotent adjacency matrix of an n-vertex graph G. For any k > 1 and

1 ≤ i = j ≤ n, ω i v i |A k |v j = k-paths w:vi→vj ω w . (3.20) Moreover, v i |A k |v i = k-cycles w based at vi ω w . (3.21)
More specifically, when i = j, the product of ω i with the entry in row i, column j of A k is a sum of basis blades indexed by k-step paths v i → v j in G. Moreover, entries along the main diagonal of A k are sums of basis blades indexed by the graph's k-cycles.

Proof. The result follows from straightforward mathematical induction on k using properties of the multiplication in RΩ n with the observation that the initial vertex of the walk, v i , is unaccounted for in v i |A k |v j , as seen in (3.18) of the matrix definition. Hence, each term of v i |A k |v j is indexed by the vertex sequence of a k-walk from v i to v j with no repeated vertices, except possibly v i at some intermediate step. Left multiplication by ω i thus sieves out the k-paths.

Considering entries along the main diagonal of A k , note that the final step of a k-cycle based at v i returns to v i so that left multiplication by ω i is not required for cycle enumeration. Note that ω i v i |A k is a row vector whose nonzero entries represent all k-paths with initial vertex v i . Similarly, A|v j is a column vector whose nonzero entries represent 1-paths with terminal vertex v j . Computing the (k + 1)-paths from v i to v j then requires computing ω i v i |A k A|v j . Define the indicator function It follows immediately that the number of multiplications performed in determining ω i v i |A k |v j is bounded above by the number of paths of length at most k -1 having initial vertex v i . Hence, the next corollary is obtained. Corollary 3.11. Given a fixed pair of vertices v 0 and v ∞ , the complexity of enumerating all k-paths from v 0 to v ∞ with the path-identifying nilpotent adjacency matrix is O(n |{paths of length k -1 or less having initial vertex v 0 }|).

χ (v ℓ ,vj )∈E by χ (v ℓ ,vj )∈E = 1 if (v ℓ , v j ) ∈ E, 0 otherwise. ( 3
Remark 3.12. The computational complexity stated above is in terms of basis blade multiplications performed within the algebra. Recall that for disjoint ordered multi-indices p, q, the product ω p ω q = ω p.q is given by sequence concatenation. Hence, some additional polynomial cost is associated with the implementation of the algebra multiplication.

Operator calculus approach to multi-constrained paths

The goal now is to extend the path-identifying nilpotent adjacency matrix approach to include weighted edges. In particular, each edge of the graph will be weighted by an m-tuple of nonnegative integers. In this manner, paths in the graph will have associated m-dimensional additive weights.

Given vectors x = (x 1 , . . . , x m ) and y = (y 1 , . . . , y m ), the notation x ≤ y is taken to mean the following:

x ≤ y ⇔ x i ≤ y i ∀i ∈ {1, . . . , m}.
(3.24)

The strict inequality x < y is analogously defined. Given a finite graph G in which each edge is weighted with an m-tuple of nonnegative integers and a constraint vector c = (c 1 , . . . , c m ) ∈ N 0 m , the multi-constrained path problem is defined as follows.

Definition 3.13. The MCP (Multi-Constrained Path) problem is to find paths

p from v 0 to v ∞ in the graph G such that wt(p) = (vi,vj )∈p wt((v i , v j )) ≤ (c 1 , . . . , c m ) = c. (3.25)
The goal is to find the set of feasible paths P f = {p = (v 0 , . . . , v ∞ ) : wt(p) < c}, i.e., all paths from v 0 to v ∞ that satisfy multiple constraints simultaneously.

An important variant of the MCP problem that is of particular interest is the associated optimization problem. Definition 3.14. Letting P denote the set of feasible paths from v 0 to v ∞ in a weighted graph G, the Multi-Constrained Optimal Path problem (MCOP) is to find a path p = (v 0 , . . . , v ∞ ) ∈ P from v 0 to v ∞ such that wt(p) ≤ wt(q) ∀q ∈ P.

(3.26)

Note that for fixed s ∈ N m , multi-exponents appearing among basis elements ν x ∈ Cℓ s nil must satisfy x < s. Given an m-vector s ∈ N m , multiplication of arbitrary s-zeon blades consequently satisfies

ν x ν y = ν x+y if x + y < s, 0 otherwise. (3.27)
In order to apply the OC approach to problems of identifying paths satisfying multiple constraints (represented by s), the path-identifying nilpotent adjacency matrix will be extended by allowing entries from the algebra Cℓ s nil ⊗ RΩ n . In this approach, a path u = (u 1 , . . . , u m ) of weight x ∈ N 0 m will be represented in Cℓ s nil ⊗Ω n by an element of the form ν x ω u . The concatenation of this path with another path v = (v 1 , . . . , v ℓ ) of weight y is then represented by the product

(ν x ω u )(ν y ω v ) = ν x+y ω u .v if u .
v is a path of multi-weight less than s, 0 otherwise.

(3.28)

Feasible & optimal paths in m-weighted graphs

Often in routing problems, there are costs associated with each edge of a graph. Connections between nodes may require time, energy, money, etc. Given an initial vertex v 0 and terminal vertex v t in a weighted graph, the collection of feasible paths from v 0 to v t refers to all paths whose associated total costs satisfy some predefined constraints. Among these feasible paths, an optimal path can then be chosen. Given a vector of constraints c = (c 1 , . . . , c m ) ∈ N m , properties of the c-zeon algebra Cℓ c nil can be used to sieve out the feasible paths from the collection of all paths. The feasible paths can then be ranked and an optimal path chosen. Definition 3.15. Let c ∈ N m , and let G = (V, E) be a graph on n vertices whose edges (v i , v j ) are multi-weighted by vectors w ij ∈ N 0 m . The c-constrained path-identifying nilpotent adjacency matrix associated with G is the n×n matrix with entries in (Cℓ c nil ) m ⊗ Ω n determined by

Ψ ij = ν wij ω j if (v i , v j ) ∈ E, 0 otherwise. (3.29)
Theorem 3.16. Given a multi-weighted graph G on n vertices with nilpotent multi-weighted adjacency matrix Ψ, a vector of constraints c = (c 1 , . . . , c k ), and a pair of distinct vertices v 0 and v ∞ , the collection of feasible paths v 0 → v ∞ in G is given by

ν 0 ω 0 n ℓ=1 v 0 |Ψ ℓ |v ∞ = paths p:v 0 →v∞ wt(p)<c ν wt(p) ω p . (3.30)
More specifically, feasible paths exist if and only if ν 0 ω 0

n ℓ=1 v 0 |Ψ ℓ |v ∞ is nonzero. For the case v 0 = v ∞ , one has v 0 |Ψ ℓ |v 0 = cycles p:v 0 →v 0 wt(p)<c ν wt(p) ω p . (3.31)
Proof. The result follows from Theorem 3.10 in consideration of combinatorial properties of Cℓ s nil .

Corollary 3.17.

If ν 0 ω 0 n ℓ=1 v 0 |Ψ ℓ |v ∞ = 0, then the optimal path p = (v 0 , . . . , v ∞ )
exists and is given by

℧ ν 0 ω 0 n ℓ=1 v 0 |Ψ ℓ |v ∞ = ν wt(p) ω p . (3.32)
Proof. By Theorem 3.16, the collection of all feasible paths v 0 → v ∞ is given by

ν 0 ω 0 n ℓ=1 v 0 |Ψ ℓ |v ∞ .
By the chosen ordering of paths and definition of ℧, the optimal path is as stated.

Figure 1: Weighted graph on six vertices.

Example 3.18 (Static case). Figure 1 depicts a weighted graph on six vertices. Each edge is assigned a vector of random positive weights in N 4 . The graph is not symmetrically weighted; i.e., weights of vertex pairs vary depending on direction. The constraint vector is c = (40, 27, 30, 40). The weighted pathidentifying adjacency matrix associated with the graph is given by the following: Let P f 0 denote the collection of feasible paths with source v 0 . That is,

0 Ν 1,4, 2,9 Ω 2 0 Ν 6,9,8,1
Ω 4 0 0 Ν 7,8,10,2 Ω 1 0 0 0 0 Ν 1,8,1,10 Ω 6 0 0 0 Ν 5,1 ,8,10 Ω 4 Ν 2,5,1,9 Ω 5 Ν 6,7,10,4 Ω 6 Ν 4
P f 0 = {p = (v 0 , . . . , v * ) : wt(p) < c, v * ∈ V }. (3.33)
Following the approach of Corollary 3.11, the number of multiplications required in computing ν 0 ω 0 v 0 |Ψ ℓ |v ∞ is seen to depend on the number of paths of length ℓ -1 or less having initial vertex v 0 and simultaneously satisfying the constraints represented by c. The following corollary is obtained as an immediate consequence.

Corollary 3.19. Given a fixed pair of vertices v 0 and v ∞ , the complexity of computing the optimal feasible path from v 0 to v ∞ via the operator calculus method is O(n |P f 0 |).

The dynamic multi-constrained path problem

Given a finite set V and positive integer r, an r-dimensional weighting function ϕ is defined on ordered pairs of vertices; i.e., ϕ : V × V → N 0 r . A weighted graph is then defined by the ordered triple G = (V, E, ϕ), where the weight of an arbitrary edge (v i , v j ) ∈ E is the r-vector ϕ((v i , v j )).

Given a fixed collection of vertices V and fixed weighting function ϕ : V × V → N 0 r , a weighted graph process on V is defined as a sequence of weighted finite graphs (G t ) := ((V, E t , ϕ)). Given a pair of vertices v 0 , v ∞ ∈ V and a vector c ∈ N r of constraints, the problem being considered is to recover the set of feasible paths v 0 → v ∞ simultaneously satisfying the constraints c.

Note that any graph sequence (G t ) naturally induces an associated sequence of weighted path-identifying nilpotent adjacency operators (Ψ t ), where each Ψ t has entries in Cℓ c nil ⊗ Ω |V | . It is important to note that a feasible path p from v 0 to v ∞ is not necessarily unique in the dynamic case because the sequence p can be partitioned into steps occurring in different frames of the process. As a consequence, scalar coefficients must be introduced to represent the path multiplicity within the collection P.

Theorem 3.20. The collection of feasible paths of length k ≥ 1 from v 0 → v ∞ requiring m or fewer frames is given by

ν 0 ω v0 0≤ℓ 1 ,...,ℓm ℓ 1 +•••+ℓm=k v 0 |Ψ 1 ℓ1 • • • Ψ m ℓm |v ∞ =
k-paths p=(v0,...,v∞)∈P α p ν wt(p) ω p (3.34) where α p is a scalar coefficient representing the multiplicity of the path p in the collection P.

Proof. Note that k ≥ 1 ensures that at least one of the integers ℓ i is nonzero. For fixed nonnegative integers ℓ 1 , ℓ 2 with ℓ

1 + ℓ 2 = k, ν 0 ω v0 v 0 |Ψ 1 ℓ1 Ψ 2 ℓ2 |v ∞ = ν 0 ω v0 vj =v0 v 0 |Ψ 1 ℓ1 |v j v j |Ψ 2 ℓ2 |v ∞ = vj =v0    ℓ 1 -paths p=(v 0 ,...,v j )∈P In frame 1 α p ν wt(p) ω p       ℓ 2 -paths q=(v j ,...,v∞ )∈P In frame 2 α q ν wt(q) ω q    = k-paths p=(v0,...,v∞)∈P ℓ 1 ,ℓ 2 α p ν wt(p) ω p , (3.35)
where P ℓ1,ℓ2 denotes the collection of feasible paths from v 0 to v ∞ in which ℓ 1 steps occur in frame 1 and ℓ 2 steps occur in frame 2. Proceeding by induction, the result is established for fixed m-tuple of nonnegative integers (ℓ 1 , . . . , ℓ m ) with

ℓ 1 + • • • + ℓ m = k.
Assuming that for positive integer m 0 and fixed non-negative integers ℓ 1 , . . . , ℓ m0 with ℓ

1 + • • • ℓ m0 = k ′ ≤ k, ν 0 ω v0 v 0 |Ψ 1 ℓ1 • • • Ψ m0 ℓm 0 |v ∞ =
k ′ -paths p=(v0,...,v∞)∈P ℓ 1 ,...,ℓm 0 α p ν wt(p) ω p , (3.36) it follows that for ℓ m0+1 = kk ′ , one has where α p denotes multiplicity of path p ∈ P.

ν 0 ω v0 v 0 |Ψ 1 ℓ1 • • • Ψ m0 ℓm 0 Ψ m0+1 ℓm 0 +1 |v ∞ = ν 0 ω v0 vj =v0 v 0 |Ψ 1 ℓ1 • • • Ψ m0 ℓm 0 |v j v j |Ψ m0+1 ℓm 0 +1 |v ∞ = vj =v0     k ′ -
Proof. First, consider the degenerate case ℓ i = 0 for 1 ≤ i ≤ m. In this case, the product of Ψ i 's is the identity operator, and v 0 = v ∞ gives 0 on the left-hand side of the equation. Observing that the maximum path length is n = |V |, the rest follows from Theorem 3.20.

Once the collection of feasible paths is obtained, the optimal path can be selected based on an ordering of the multi-exponents.

Conclusion

The operator calculus approach provides convenient symbolic computational tools for a broad range of combinatorial problems and practical applications. The methods developed here can be applied directly to multi-constrained quality of service (QoS) problems as well as problems related to precomputed routing in store-and-forward satellite constellations [START_REF] Cruz-Sánchez | Operator calculus approach to minimal paths: Precomputed routing in a store-and-forward satellite constellation[END_REF].

  .22)Letting m denote the number of multiplications involved in this computation, v i → v ℓ }χ (v ℓ ,vj )∈E ≤ |{k-paths with source v i }|.(3.23)

  paths p=(v 0 ,...,v j ) in P ℓ 1 ,...,ℓm 0 α p ν wt(p) ω p k ′ )-paths q=(v j ,...,v∞ ) in P 0,...,ℓ m 0 +1 α q ν wt(q) ω q p=(v0,...,v∞)∈Pℓ 1 ,...,ℓ m 0 +1 α p ν wt(p) ω p , (3.37)Hence, the result is established for positive integer m:ν 0 ω v0 v 0 |Ψ 1 ℓ1 • • • Ψ m ℓm |v ∞ = k-paths p=(v0,...,v∞)∈P ℓ 1 ,...,ℓm α p ν wt(p) ω p . (3.38)The proof is thus completed by summing over all such m-tuples.Corollary 3.21. The collection of feasible paths of all lengths from initial vertex v 0 to terminal vertex v ∞ = v 0 requiring m or fewer frames is recovered from the canonical expansion ofν 0 ω v0 0≤ℓ 1 ,...,ℓm ℓ 1 +•••+ℓm ≤n v 0 |Ψ 1 ℓ1 • • • Ψ m ℓm |v ∞ =paths p:v0→v∞ α p ν wt(p) ω p , (3.39)

Corollary 3 .0≤ℓ 1

 31 22. Given a preferential ordering of multi-exponents, the optimal path p = (v 0 , . . . , v ∞ ) from v 0 to v ∞ in the first m frames of the graph sequence (G t ) is given by℧ ,...,ℓm 1≤ℓ 1 +•••+ℓm≤n v 0 |Ψ 1 ℓ1 • • • Ψ m ℓm |v ∞    = α p ν wt(p) ω p , (3.40) provided ν 0 ω v0 v 0 |Ψ 1 ℓ1 • • • Ψ m ℓm |v ∞ = 0for some m-tuple of nonnegative integers ℓ 1 , . . . , ℓ m , not all of which are zero. Here, α p denotes the multiplicity of path p ∈ P.


Pierre Riedinger
email: pierre.riedinger@ensem.inpl-nancy.fr

Irinel-Constantin Morȃrescu
email: constantin.morarescu@ensem.inpl-nancy.fr

A numerical framework for optimal control of switched affine systems with state constraint

Keywords: Singular optimal control, switched systems, complementarity system, numerical method I

In this paper, we address the problem of numerical implementation of optimal control for switched affine systems with state constraint. In order to solve the problem we introduce a relaxed system and establish the connection between the solution of this system and the solution of the initial one. One of the main difficulties is then related to the fact that the optimal solution of the relaxed system is generally singular. We show that, using slack variables, a set of complementarity constraints can be used to take into account the singular nature of the solution. This formulation does not requires a priori knowledge on the structure (regular/singular) of the solution. The optimal control problem is then reformulated as a constraints optimization problem over the Hamiltonian systems and solved via a direct method. In addition, state path constraints are included. Numerical simulations for boost, buck-boost and Flying Capacitor converters, both in continuous and discontinuous conduction mode, illustrate the effectiveness of the proposed methodology.

for which generally no common equilibrium can be defined. In this context, the referred targets, named operating points, are defined as the equilibria of a relaxed system -obtained by expanding the control domain to its convex hull. As a result of the applied control strategy, the average value of the state variables computed on a sliding window coincides with the target [START_REF] Hauroigné | Switched affine systems using sampled-data controllers: robust and guaranteed stabilization[END_REF].

It is worth noting that for switched systems only suboptimal solutions can be defined from the optimal solution of the relaxed systems [START_REF] Patino | Practical optimal state feedback control law for continuous-time switched affine systems with cyclic steady state[END_REF]. Moreover, the optimal solution of the relaxed system is generally singular as deduced from the Pontryagin Maximum Principle (PMP). Precisely, the necessary conditions are inconclusive on the control value when H u = 0, ∀u. This is a well known [START_REF] Robbins | A generalized Legendre-Clebsch condition for the singular case of optimal control[END_REF], [START_REF] Krener | The Higher Order Maximal Principle and its Application to Singular Extremals[END_REF], [START_REF] Michel | Singular Optimal Control ? The State of the Art[END_REF], [START_REF] Bonnard | Singular trajectories and their role in control theory[END_REF] difficulty encountered in optimal control field. The search of analytic solutions is pervasive even for low dimensional systems [START_REF] Patino | Practical optimal state feedback control law for continuous-time switched affine systems with cyclic steady state[END_REF]. In this context only few methods have been proposed to numerically solve such control problems. Indirect shooting method requires a priori knowledge on the structure (regular/singular) of the trajectory [START_REF] Bonnans | Singular Arcs in the Generalized Goddard's Problem[END_REF]. This information can be achieved using direct methods [START_REF] Betts | Practical methods for optimal control using nonlinear programming[END_REF], [START_REF] Riedinger | About solving hybrid control problems, 17th IMACS World Congress[END_REF] but the determination of the control lacks of accuracy. Once the information on the structure of the optimal solution is available, an accurate solution can be obtained if singular arc conditions are imposed for the singular part of the trajectory; see for example the Goddard Rocket Problem in [START_REF] Rao | Algorithm 902: GPOPS, A MATLAB Software for Solving Multiple-Phase Optimal Control Problems Using The Gauss Pseudospectral Method[END_REF].

In order to overcome this inconvenient, we propose to reformulate the control set inclusion deduced from the necessary conditions of PMP as a set of complementarity constraints and to add them to the Hamiltonian system to form a complementarity system [START_REF] Facchinei | Finite-dimensional Variational Inequalities and Complementarity Problems[END_REF]. The multivalued part of the control (corresponding precisely to the singular cases) is then concisely expressed in these complementarity constraints. At this stage, one may view the problem as a Boundary Value Problem (BVP) and try to use a shooting method. Unfortunately, such shooting methods are not appropriate to solve this BVP because of the non uniqueness of the solution.

In this paper we propose to use a direct method [START_REF] Betts | Practical methods for optimal control using nonlinear programming[END_REF] to minimize the cost function subject to the complementarity system. This approach can be classified as mixed direct-indirect method since it uses both the necessary condition of PMP and a large-scale nonlinear programming problem (NLP) arising from the discretization of the cost and of the resulting complementarity systems.

The main advantage of this formulation is that no information on the structure (regular/singular) of the trajectories is a priori required.

To improve the accuracy of numerical results, additional constraints concerning the order of the singularity are taken into account. Finally, the method is extended to path constraints case.

The paper is organized as follows. In Section II one gives the system description and defines the operating points. Section III formulates the optimal control problem derived from PMP. The singular nature of the solution is also detailed and discussed. The complementarity formalism is used to reformulate the optimal control problem in Section IV. To increase the accuracy of the control determination we add constrains related to the order of the singularities. Section (V) contains numerical simulation illustrating the results of the proposed methodology on a variety of power converters.

II. SYSTEM DESCRIPTION

The class of constraint affine systems under consideration can be described by:

ẋ(t) = f (x(t)) + m i=1 u i (t)g i (x(t)) (1)
x(0) = x 0 (2) c(x(t)) ≤ 0 (3)
where

u(t) = (u 1 (t), . . . , u m (t)) ∈ R m is the control law and x(t) ∈ R n represents the state value at time t. The maps f (•) : R n → R n , c(•) : R n → R r and g i (•) : R n → R n , i = 1, . . . , m
are supposed sufficiently smooth. A vector is considered positive/negative if all its components are positive/negative. We note that (3) defines the unilateral state-path constraints.

When u(t) ∈ U = {0, 1} m equation (1) describes a class of Switched System (SS) largely used to express the dynamics of nonsmooth mechanical systems (see [START_REF] Moreau | Unilateral contact and dry friction in finite freedom dynamics[END_REF], [START_REF] Morȃrescu | Passivity-based tracking control of multiconstraint complementarity Lagrangian systems[END_REF] and the reference there in) or electric circuits ( [START_REF] Çamlibel | Stability and controllability of planar bimodal linear complementarity systems, 42nd[END_REF]). On the other hand when u(t) ∈ co(U ) = [0, 1] m where co(U ) stands for the convex hull, equation (1) describes a Relaxed System (RS). The (SS) belongs to the class of nonsmooth systems for which the notion of solution can be properly defined and generalized in the sense of Fillipov [START_REF] Filippov | Differential equations with discontinuous righthand sides[END_REF]. The link between the solutions of (SS) and (RS) is established by a density theorem in infinite time [START_REF] Ingalls | An infinite-time relaxation theorem for differential inclusions[END_REF]. Let us denote L ∞ ([0, +∞), U ) the set of all essentially bounded measurable functions mapping [0, ∞) to U . Then L ∞ ([0, +∞), U ) equipped with the essential supremum norm (i.e. f ∞ = ess sup|f |) becomes a Banach space.

In [START_REF] Ingalls | An infinite-time relaxation theorem for differential inclusions[END_REF] is proven that there exists a switching law u ∈ L ∞ ([0, +∞), U ) such that the trajectory of (RS) is approached as close as desired by the one of (SS). For this reason, the operating points set of (SS) denoted by X ref , is defined as the set of equilibrium points of (RS):

X ref = x ref ∈ R n : f (x ref ) + (4) m i=1 u i,ref g i (x ref ) = 0, u i,ref ∈ [0, 1] . (5)
This set defines the control targets for the state of (SS).

It

III. OPTIMAL CONTROL PROBLEM FORMULATION

In the sequel, we propose a methodology to numerically solve the following optimal control problem for the affine system defined by equations ( 1)-(3):

min u(.) T 0 L(x(t), x ref )dt ( 6
)
where L is the continuously differentiable performance function, x ref stands for a chosen equilibrium point in X ref and T refers to the strictly positive final time and is fixed (T ∈ R + ).

When the (RS) has a bang-bang solution (i.e. ∀t, u(t) ∈ {0, 1} m ), (SS) is also solved.

Otherwise, it exists a non zero measure set (in the Lebesgue sense) T such that the optimal solution takes values u(t) ∈ co(U )\U for all t ∈ T . Although these solutions are not admissible for (SS), Pulse Wise Modulation (PWM), for example, yields an average approximation by high switching frequency [START_REF] Ingalls | An infinite-time relaxation theorem for differential inclusions[END_REF]. Consequently we relax the optimal control problem for (SS) by searching the optimal solutions of (RS) from which suboptimal solutions for (SS) can be derived.

Suppose in a first time that no constraint are taken into account and let inspect the necessary conditions defined by the Pontryagin Maximum Principle (PMP). The Hamiltonian function H is given by (abnormal case is not taken into account):

H(t) = λ ⊤ (t)(f (x) + g(x)u(t)) + L(x(t), x ref ) (7)
where λ(t) ∈ R n depicts the co-state whose dynamic is given by the Hamiltonian system λ

(t) = -H x (t) ( ẋ(t) = H λ (t)).
Optimality condition of PMP that u minimizes the Hamiltonian function leads to the following inclusion,

u i (t) ∈                0 if φ i (t) > 0 [0 1] if φ i (t) = 0 1 if φ i (t) < 0 ∀ i = 1, ..., m (8)
where

φ(t) = H u (t) = λ T (t)g(x(t))
defines the switching function. When u i (t) = 0 or 1, the control is called regular. A singularity arises when a component of φ i (t) vanishes identically on

a time interval [a, b], b > a.
In this case, the PMP is inconclusive concerning the control value

u i (•) on [a, b].
This situation is referred to as singular control [START_REF] Robbins | A generalized Legendre-Clebsch condition for the singular case of optimal control[END_REF], [START_REF] Bryson | Applied Optimal Control[END_REF], [START_REF] Bonnard | Singular trajectories and their role in control theory[END_REF] and it corresponds here to the case where u(t) takes values in co(U )\U, t ∈ [a, b]. We can conclude that solution segments for (RS), which are not admissible for (SS), involve a singular control.

In the sequel, in order to simplify the notation we shall often avoid to emphasize the time dependency of the variables and controller.

H u (x, λ, u) = 0, ∀u ∈ co(U ), ∀t ∈ (a, b). (9)
To determine these arcs, successive time differentiations:

d q dt q H u (x, λ, u) = 0, ∀u ∈ co(U ), ∀t ∈ (a, b), (10)
are required. The differentiation has to be done until at least a component of the control appears.

The smallest positive integer q * such that

∂ ∂u ( d q * dt q * H u (x, λ, u)) = 0, (11)
if it exists, is called the order of singularity [START_REF] Michel | Singular Optimal Control ? The State of the Art[END_REF] (or index of the DAE). As we have shown in [START_REF] Patino | Practical optimal state feedback control law for continuous-time switched affine systems with cyclic steady state[END_REF], for low dimensional systems (i.e. n ≤ 3), the analytic expression of u can be explicitly determined by using the set of algebraic equations [START_REF] Facchinei | Finite-dimensional Variational Inequalities and Complementarity Problems[END_REF] and a second order necessary condition given by the Generalized Legendre-Clebsh Condition [START_REF] Robbins | A generalized Legendre-Clebsch condition for the singular case of optimal control[END_REF], [START_REF] Michel | Singular Optimal Control ? The State of the Art[END_REF], [START_REF] Krener | The Higher Order Maximal Principle and its Application to Singular Extremals[END_REF]:

(-1) ∂ ∂u ( d q * dt q * H u (x, λ, u)) ≥ 0. (12)
On one hand, to determine analytically the optimal solutions including singular arcs is a particularly hard problem and only low dimensional systems (i.e. n=2,3) have been treated [START_REF] Patino | Practical optimal state feedback control law for continuous-time switched affine systems with cyclic steady state[END_REF].

On the other hand, standard indirect numerical methods such as multiple shooting methods are not appropriate to deal with singular arcs without a priori information on the structure of the trajectories. Knowledge on the singular arcs are of few help in determination of the junction time between singular and regular arcs. Furthermore their number may be infinite. In this context, we propose to use direct method and an ad hoc formulation that enable a numerical treatment in presence of singular arcs that do not require a priori information.

IV. OPTIMAL CONTROL USING COMPLEMENTARITY SYSTEM FORMALISM

Complementarity is a useful concept that can be used to describe properly the behavior of nonsmooth systems. Complementarity systems have been largely employed in the analysis and control of mechanical systems [START_REF] Brogliato | Nonsmooth Mechanics[END_REF], [START_REF] Brogliato | Dissipative Systems Analysis and Control. Theory and Applications[END_REF], [START_REF] Morȃrescu | Passivity-based tracking control of multiconstraint complementarity Lagrangian systems[END_REF] and recently they are used to describe the dynamics of electric circuits [START_REF] Çamlibel | Stability and controllability of planar bimodal linear complementarity systems, 42nd[END_REF], [START_REF] Vasca | A new perspective for modeling power electronics converters: complementarity framework[END_REF], [START_REF] Acary | Time-stepping numerical simulation of switched circuits with the nonsmooth dynamical systems approach[END_REF].

Definition 2. A Linear Complementarity Problem (LCP) is a system given by:

           y ≥ 0 Ay + b ≥ 0 y T (Ay + b) = 0 (13)
which is compactly rewritten as

0 ≤ y ⊥ Ay + b ≥ 0 (14)
Such an LCP has a unique solution for all b if and only if A is a P-matrix, i.e. all its principal minors are positive [START_REF] Facchinei | Finite-dimensional Variational Inequalities and Complementarity Problems[END_REF].

Combining complementarity conditions with differential equations we get what usually in the literature is referred to as complementarity systems (see [START_REF] Brogliato | Some perspectives on the analysis and control of complementarity systems[END_REF] and the references therein):

               ẋ(t) = f (x(t), u(t)) y(t) = h(x(t), u(t)) 0 ≤ y(t) ⊥ u(t) ≥ 0 (15)
We also note that in some applications, the complementarity system may contain an external input (forcing term) or the adjoint variable [START_REF] Schumacher | Complementarity systems in optimisation[END_REF]. In the sequel, we propose to reformulate the Hamiltonian system and the necessary condition (8) as a complementarity system.

A. Relay inclusion as a complementarity constraint

Let us recall (see [START_REF] Schumacher | Complementarity systems in optimisation[END_REF]) that the relay-type function ( 8) can be rewritten as a complementarity constraint and an algebraic equation as follows:

0 ≤ s 1 ⊥ u ≥ 0 (16) 0 ≤ s 2 ⊥ (1 -u) ≥ 0 (17) 0 = s 1 -s 2 -φ (18)
where s 1 and s 2 define the positive and the negative part of φ(x, λ). In other words, the components of the two vectors s 1 and s 2 are defined by:

s 1i (t) =        0 if φ i (t) < 0 φ i (t) if φ i (t) ≥ 0 ∀ i = 1, . . . , m s 2i (t) =        -φ i (t) if φ i (t) < 0 0 if φ i (t) ≥ 0
It is important to note that as long as s 1i > 0 or s 2i > 0 the component u i of the control law u is uniquely determined. When s 1i = 0 & s 2i = 0 the uniqueness of u i is lost and the choice of its value is a matter of optimization. We also observe that the situation s 1i • s 2i > 0 cannot occur since this will contradict either [START_REF] Ingalls | An infinite-time relaxation theorem for differential inclusions[END_REF] or [START_REF] Lee | Receding horizon control of switching systems[END_REF].

B. Path constraints

Until now, no consideration has been made on the path constraints (3). Like for singular controls, the index q * of the state constraint is the first time where u appears explicitly in the successive time differentiation of c along the trajectories.This index allows to define the mixed control-state path constraints:

p(x, u) = d q * dt q * c(x) ≤ 0. ( 19
)
and the Lagrangian function:

L(x, u, λ, µ) = H(x, u, λ) + µ T p(x, u) (20)
where µ ≥ 0 represents the vector of Lagrange multipliers.

Theorem 1. If u * is the optimal control for the state constraint affine system (1-3) when the cost is defined by ( 6) then u * satisfies the following constraint optimization problem:

COP : min u(•) T 0 L(x(t), x ref )dt (21)
s.t. ẋ = L λ λ = -L x (22)
x(0) = x 0 λ(T ) = 0 (23)

0 ≤ s 1 ⊥ u ≥ 0 (24) 0 ≤ s 2 ⊥ (1 -u) ≥ 0 (25) 0 = s 1 -s 2 -L u (x, λ, µ) (26)
0 ≤ µ ⊥ (s 1 + s 2 ) ≥ 0 (27) 0 ≤ µ ⊥ -c(x) ≥ 0 (28)
Comments:

• Generally when path constraints are considered, jumps occur on the co-state at junction time τ . They are characterized by λ + (τ ) = λ -(τ ) + α∇c(x(τ )). In fact, α can be taken equals to zeros when the motion follows the border c(x) = 0 for t ≥ τ as proven in the pioneer work of [START_REF] Pontryagin | The Mathematical Theory of Optimal Processes[END_REF].

• When µ equals zeros [START_REF] Patino | Practical optimal state feedback control law for continuous-time switched affine systems with cyclic steady state[END_REF] allows us to retrieve the unconstraint version of PMP. The complementarity relation [START_REF] Robbins | A generalized Legendre-Clebsch condition for the singular case of optimal control[END_REF] assures that µ = 0 when c(x) is not active (i.e. c(x) > 0).

• The constraint c(x) is activated (i.e. c(x) = 0) when µ becomes positive and modifies the co-state dynamics.

• The complementarity condition 0 ≤ µ ⊥ (s 1 + s 2 ) ≥ 0 is necessary to ensure L u = 0 on the constraint c(x) = 0.

It has to be mentioned here two important features of the above formalism. Firstly, this formulation clearly takes into account both singular arcs and path constraints. Secondly, no assumption is made on the structure (singular/regular) of the trajectories. Consequently, we do not requires any initial information on the existence of singular arcs entering in the optimal trajectory.

The optimization problem (COP) can be solve using direct method i.e. by minimizing the cost function [START_REF] Perko | Differential equations and Dynamical systems[END_REF] subject to the complementarity system [START_REF] Raimondo | Min-max model predictive control of nonlinear systems: A unifying overview on stability[END_REF][START_REF] Sira-Ramirez | Sliding mode control of DC-to-DC power converters using reconstructors[END_REF][START_REF] Schumacher | Complementarity systems in optimisation[END_REF][START_REF] Sontag | New characterizations of input-to-state stability[END_REF][START_REF] Vasca | A new perspective for modeling power electronics converters: complementarity framework[END_REF][START_REF] Bonnard | Singular trajectories and their role in control theory[END_REF][START_REF] Robbins | A generalized Legendre-Clebsch condition for the singular case of optimal control[END_REF]. Of course, in order to determine accurate solution, recall that it is important to scaling the problem and to use code that use multiple phase formulation to reduce the sensitivity and that enable grid refinement procedure.

For example, the promising beta version GPOPS Toolbox [START_REF] Garg | A Unified Framework for the Numerical Solution of Optimal Control Problems Using Pseudospectral Methods[END_REF] for Matlab has been used for all the examples given in the next section. From a practical point of view, complementary constraints is achieved by penalization term in the cost.

A practical and significant improvement in minimizing the cost function ( 21) subject to the complementarity system (22-28) can be obtained by including the algebraic constraints [START_REF] Facchinei | Finite-dimensional Variational Inequalities and Complementarity Problems[END_REF] for q = 1, • • • , q * . The last constraint when q = q * is particularly beneficial since the control appears in it explicitly. Of course, as these constraints correspond to the case where the control is singular, these additions are only possible when path constraints are not active (that is when µ = 0 and λ = -H x = -L x ). If the path constraints are active, they are unnecessary since u appears this time in the path constraints.

Theorem 2. If u * is the optimal control of problem defined by cost [START_REF] Çamlibel | Stability and controllability of planar bimodal linear complementarity systems, 42nd[END_REF] and state constraint affine system (1-3) then u * satisfies COP of Theorem 1 and the following additional constraints related to the order of the singularity (10):

0 = s 2q+1 -s 2q+2 - d q dt q φ, q = 2, • • • , q * (29) 0 ≤ p i=1 c i (x)u(1 -u) ⊥ (s 2q+1 + s 2q+2 ) ≥ 0 (30) 0 ≤ s 2q+1 ⊥ s 2q+2 ≥ 0 ( 31
)
where s 2q+1 , s 2q+2 q = 2, • • • , q * are additional slack variables.

It is easy to show from equation ( 29) and ( 31) that s 2q+1 and s 2q+2 are respectively the positive and the negative part of d q dt q φ. Then, [START_REF] Riedinger | About solving hybrid control problems, 17th IMACS World Congress[END_REF] imposes that the i th component of | d q dt q φ| equals zero when u i is singular and the state x does not satisfy one of the constraints given by (3). Following singular arc definition, these relations are fulfilled and their inclusion in (OCP) improves the V. APPLICATION TO POWER CONVERTERS DC-DC converters such as buck, boost, buck-boost, flying capacitor converters enter clearly in the class of constraint affine switched systems. We have applied our numerical approach for all these converters and we present in the sequel two selected cases to illustrate the effectiveness of the proposed procedure.

A. FCC converter

In the case of a few megawatt industrial power applications, the classical power converters have a very high voltage in the switching components (several kilovolts). To compensate this, a new class of power converters called Flying Capacitor Converters (FCC) appeared (or multicellular converter). The structure of the FCC converter reduces the voltage throughout the switches. It is composed of serial connections between semiconductor switching devices and passive storage elements to achieve the target operating voltage as illustrated in Figure 1.

Three switching cells can be isolated, each of them being able to operate in two complementary mode. The behavior of each cell can be described using only one boolean control variable u i ∈ {0, 1} with i = 1, 2, 3. u i = 1 means that the upper switch is closed and the lower switch is open whereas u i = 0 means that the upper switch is open and the lower switch is closed (Figure 1). The state equations of the converter have an affine form given by

      ẋ1 ẋ2 ẋ3       =       0 0 -R L x 3       +       -x 3 C 1 x 3 C 1 0 0 -x 3 C 2 x 3 C 2 x 1 L x 2 -x 1 L E-x 2 L             u 1 u 2 u 3       (32)
where x 1 , x 2 are the voltage on each capacitor and x 3 the load current. The chosen parameters are E = 30 V , R = 10 Ω, C 1 = C 2 = 40µF and L = 10 mF .

The capacitor voltages must be balanced to 2E/3 and E/3, while the demanding load current is fixed to 0.6 A. Then,

x ref = 2E/3 E/3 0.6 T .
The optimization problem consists in minimizing the quadratic criteria:

min u(•) T 0 (x -x ref ) ′ Q(x -x ref )dt.
The weight matrix Q is taken to:

Q =       2 -1 0 -1 1 0 0 0 200      
Figure 2 and 3 show the control and state trajectories starting from the origin. As expected, the control are totally singular and are properly determined using the proposed formulation. In 6) and (1-3) using the same GPOPS toolbox, yields the results provided in Figure 4 and Figure 5. The solver returns a flag indicating that the optimization has finished successfully and that the optimality conditions are satisfied. In fact, we can observe that the trajectories are similar but the control is far from being accurate due to Hamiltonian singularities. All these simulations have been done with the same parameters: 3 phases formulations, 3 mesh iterations, analytic derivative provided and the same mesh tolerance (1e-3).

B. Boost converter

The topology of a boost converter is given in Figure [START_REF] Çamlibel | Stability and controllability of planar bimodal linear complementarity systems, 42nd[END_REF]. In continuous conduction mode, following the position of the switch u, two subsystems ẋ = A i x + B i , i = 1, 2 can be considered

with matrices A 1 =    -R L 0 0 -1 R 0 C   , A 2 =    -R L -1 L 1 C -1 R 0 C   , B 1 =    E L 0    and B 2 = B 1 .
Generally the control strategies proposed in the literature do not take into account the discontinuous mode that may appears when the demanding current in the load is low. In this case, the current in the diode vanishes and produces an autonomous switch. Few results concerning optimal control of hybrid systems are related to the case of autonomous switches. The third subsystem is then described by

A 3 =    0 0 0 -1 R 0 C    , and B 3 =    0 0    .
In fact, it is not necessary to add this last subsystem since it is not difficult to show that this subsystem matches the dynamics of ẋ = u(A

1 x + B 1 ) + (1 -u)(A 2 x + B 2 ) for a control u eq (x) ∈ [0 1].
So, the optimal control problem is solved with path constraint corresponding to Figures 8 and7 show the state trajectories and the optimal control as well as the Lagrange multiplier obtained for the weight matrix Q = diag [START_REF] Acary | Time-stepping numerical simulation of switched circuits with the nonsmooth dynamical systems approach[END_REF][START_REF] Facchinei | Finite-dimensional Variational Inequalities and Complementarity Problems[END_REF]. The time scale is expressed in ms.

The numerical results show that just at the beginning of the simulation the control switches from 1 to 0 before the optimal trajectory reaches the singular arc which happens after approximatively 1ms. When a new reference is applied at time t = 5ms, control switches to zero and the system enters in discontinuous conduction mode when the current vanishes. The value of the Lagrangien multiplier clearly indicates this while the control values met the constraint ẋ1 = 0. The system returns in continuous conduction mode with a short activation of mode 2 at time t = 7.5ms.

The control becomes singular once the reference is reached. Now, Figure 9 compares the precision with which the optimal control is determined according to the formulation adopted to solve the optimal control problem. Starting from the origin, the cost function consists to minimize the quadratic criteria min u(.)

C. buckboost converter

The last example concerns a buck-boost converter and it is chosen to show that this numerical approach can be useful not only to compute optimal open-loop control (in case of predictive control for example) but also to design robust close-loop controller.

In Figure 10 the red circle is the reference and the optimal trajectories starting from different initial conditions are plotted in blue. We can observe that discontinuous conduction mode appears along the negative part of x 2 -axis and that all the trajectories reach the target x ref following a singular arcs. As Figure 10 shown, the state space is partitioning following the control to be used. It can be observed that the singular trajectories, at the junction between the different parts, defined optimal sliding surfaces.

Clearly, there is a partition that defines u as a state feedback u(x). We think that this partition can be achieved using learning machines such as neural networks.

The parameters used in the simulations above are:

E = R = L = C = 1, A 1 =    0 0 0 -1 RC    A 2 =    0 1 L -1 C -1 RC    B 1 =    E L 0    B 2 =    0 0  
 weight matrice Q = diag(q 1 , q 2 ) with q 1 = 1 and q 2 = 10.

VI. CONCLUSION

In this paper we have proposed a numerical framework to solve optimal control problem for affine switched systems. One specific feature of this class is that singular controls have to be considered in order to obtain the optimal solution. As singular controls are determined by an inclusion depending on the state and the co-state, we have reformulated the optimization problem using the augmented (in case of path constraints) Hamiltonian system and a set of slack variables defining a set of complementarity constraints. This constraint augmented optimization problem is then solved with a direct method dedicated to optimal control. The contribution of the proposed methodology is twofold: an accurate optimal control is determined without any additional information concerning the existence and the number of the singular arcs belonging to the optimal trajectory and, the optimal control is determined even in the presence of autonomous switches. The first contribution is an important relaxation since generally in the provided singular control examples in the literature, an a priori knowledge of the number of singular arcs is required (see for example the Godart Rocket example). The second contribution allows to complete the simulation of electric circuits containing switching elements, including the discontinuous conduction mode.

It is noteworthy that a direct transcription of the optimal control problem as a non linear program leads to poor results as the FCC converter example has shown. The effectiveness of our method has been illustrated in several examples.

Definition 1 .

 1 If a component of φ vanishes identically on the time interval [a, b], b > a we say that x belongs to a singular arc for all t ∈ [a, b]. These arcs are characterized by:

Fig. 1 .

 1 Fig. 1. Fying capacitor converter

3 Fig. 2 .

 32 Fig. 2. Optimal singular controls using Hamiltonian and complementarity constraints

3 Fig. 3 .

 33 Fig. 3. Optimal state trajectories using Hamiltonian and complementarity constraints

3 Fig. 4 . 3 Fig. 5 .

 3435 Fig. 4. Numerical results in a naive transcription : optimal controls

Fig. 7 .

 7 Fig. 7. Boost converter: Optimal controls

x 1 ≥ 0 .Fig. 8 .

 08 Fig. 8. Boost converter : State trajectories

T0Fig. 9 .

 9 Fig. 9. Boost converter: Improvement in the numerical control determination; c: direct conscription, b: Hamiltonian like a Complementarity Problem; a: Full Constraint formulation using order of the singularity

Fig. 10 .

 10 Fig. 10. BuckBoost converter: Optimal trajectories in state space. x ref : red circle

June 24, 2011 DRAFT

This work was partially supported by ANR project ArHyCo, Programme "Systèmes Embarqués et Grandes Infrastructures"-ARPEGE, contract number ANR-2008 SEGI 004 01-30011459. The research leading to these results has received funding from the European Community's Seventh Framework Programme (FP7/2007-2013) under grant agreement n • 257462: HYCON2 Network of Excellence "Highly-Complex and Networked Control Systems".


Miriam Jane Smith

Jenny E Higgs

Naomi L Bowers

Dorothy Halliday

Joan Patterson

James Gillespie

Susan M Huson

Simon R Freeman

Simon Lloyd

Scott A Rutherford

A J Wallace

R T Ramsden

Dgr Evans

Professor D Gareth R Evans
email: gareth.evans@cmft.nhs.uk

Cranial meningiomas in 411 NF2 patients with proven gene mutations: Clear positional effect of mutations, but absence of female severity effect on age at onset

Keywords: Meningioma, NF2, Merlin, mutation position effect

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

Introduction

Neurofibromatosis type 2 (NF2) is an autosomal dominant disease characterised by the development of schwannomas, meningiomas and ependymomas, with the great majority of patients having bilateral schwannomas of the eighth cranial nerve [START_REF] Evans | A clinical study of type 2 neurofibromatosis[END_REF][START_REF] Evans | A genetic study of type 2 neurofibromatosis in the United Kingdom. I. Prevalence, mutation rate, fitness, and confirmation of maternal transmission effect on severity[END_REF]. Intracranial meningiomas are the second most common type of tumour found in NF2 patients, although the factors determining which individuals will develop meningiomas are unclear.

Spinal meningiomas may occur in up to one third of NF2 patients [START_REF] Mautner | The neuroimaging and clinical spectrum of neurofibromatosis 2[END_REF][START_REF] Goutagny | Meningiomas and neurofibromatosis[END_REF]. However, they are frequently asymptomatic and many are only detected incidentally when imaging is carried out for an unrelated purpose. Spinal meningiomas are not screened for routinely and they are difficult to distinguish from schwannomas by imaging alone. Therefore, the exact prevalence is unclear.

The NF2 gene is biallelically inactivated in most schwannomas and meningiomas. The spectrum of NF2 mutations consists mainly of truncating (nonsense and frameshift) mutations [START_REF] Ahronowitz | Mutational spectrum of the NF2 gene: a meta-analysis of 12 years of research and diagnostic laboratory findings[END_REF][START_REF] Baser | The distribution of constitutional and somatic mutations in the neurofibromatosis 2 gene[END_REF]. A smaller number of patients have missense mutations, which are associated with a milder disease phenotype [START_REF] Ruttledge | Type of mutation in the neurofibromatosis type 2 gene (NF2) frequently determines severity of disease[END_REF][START_REF] Evans | Genotype/phenotype correlations in type 2 neurofibromatosis (NF2): evidence for more severe disease associated with truncating mutations[END_REF]. Although the NF2 tumour spectrum and molecular findings are reasonably well described, description of disease severity has largely been related to mutation type rather than position within the NF2 gene [START_REF] Ruttledge | Type of mutation in the neurofibromatosis type 2 gene (NF2) frequently determines severity of disease[END_REF][START_REF] Evans | Genotype/phenotype correlations in type 2 neurofibromatosis (NF2): evidence for more severe disease associated with truncating mutations[END_REF][START_REF] Parry | Germ-line mutations in the neurofibromatosis 2 gene: correlations with disease severity and retinal abnormalities[END_REF][START_REF] Kluwe | A missense mutation in the NF2 gene results in moderate and mild clinical phenotypes of neurofibromatosis type 2[END_REF][START_REF] Baser | Genotype-phenotype correlations for nervous system tumors in neurofibromatosis 2: a population-based study[END_REF][START_REF] Selvanathan | Further genotype-phenotype correlations in neurofibromatosis 2[END_REF]. A positional effect has only been described for splice site mutations, with mutations in the 3' end of the gene leading to milder disease [START_REF] Baser | The location of constitutional neurofibromatosis 2 (NF2) splice site mutations is associated with the severity of NF2[END_REF][START_REF] Kluwe | Phenotypic variability associated with 14 splice-site mutations in the NF2 gene[END_REF].

The NF2 gene product, Merlin, is a cytoskeletal protein that provides a functional link between the cell membrane and the actin cytoskeleton [START_REF] Xu | Merlin differentially associates with the microtubule and actin cytoskeleton[END_REF]. The N-terminus contains a FERM (4.1, ezrin, radixin, moesin) domain, which is responsible for binding various interacting proteins [START_REF] Turunen | Structure-function relationships in the ezrin family and the effect of tumor-associated point mutations in neurofibromatosis 2 protein[END_REF]. This can be further divided into three sub-domains. The first sub-domain (amino acids 20-100) has a ubiquitin-like structure, the second sub-domain (amino acids 101-215) has an acyl CoA binding protein-like structure and the third sub-domain (amino acids 216-313) has a structure similar to phosphotyrosine binding protein (PTB), pleckstrin homology domain (PH), or Enabled/VASP Homology 1 (EVH1). These three domains fold into a "cloverleaf" conformation [START_REF] Kang | The structure of the FERM domain of merlin, the neurofibromatosis type 2 gene product[END_REF]. The FERM domain is followed by an α-helical domain (amino acids 314-478) which can fold to bring the C-terminal domain (amino acids 479-595) close to the N-terminus. In its folded conformation Merlin is inhibited from binding actin.

The absence of somatic mutations in exons 14 and 15 of the NF2 gene in sporadic meningiomas [START_REF] Ahronowitz | Mutational spectrum of the NF2 gene: a meta-analysis of 12 years of research and diagnostic laboratory findings[END_REF][START_REF] Baser | The distribution of constitutional and somatic mutations in the neurofibromatosis 2 gene[END_REF] prompted us to investigate whether a positional effect of germline mutations existed for meningiomas. We also felt the evidence base for counselling regarding gender effect on severity of meningioma disease was also relatively under investigated for NF2. Here describe the rates of cranial meningioma onset in a cohort of 411 patients with NF2 in whom a germline or mosaic mutation has been identified and compare severity of disease with relation to the position of the NF2 mutation. We also investigate the effect of gender on meningioma risk.

Materials and Methods

We undertook a retrospective review of 411 patients from our database with a diagnosis of NF2. Patients were referred to the Regional Molecular Genetics Service at St Mary's Hospital in Manchester for mutational analysis. Only individuals with a proven germline or mosaic mutation were included. Mutation analysis of the patients' peripheral blood was undertaken as previously described [START_REF] Moyhuddin | Somatic mosaicism in neurofibromatosis 2: prevalence and risk of disease transmission to offspring[END_REF][START_REF] Wallace | Mutation scanning of the NF2 gene: an improved service based on meta-PCR/sequencing, dosage analysis, and loss of heterozygosity analysis[END_REF]. Briefly, genomic DNA samples were amplified with primers for all 17 exons of the NF2 gene and direct sequencing was used in addition to Multiple Ligation Dependant Probe amplification (MLPA) to identify mutations. Where no mutation was identified, further analysis using DNA of the patient's tumour tissues, where available, was undertaken using a previously described technique [START_REF] Moyhuddin | Somatic mosaicism in neurofibromatosis 2: prevalence and risk of disease transmission to offspring[END_REF].

All patients have had cranial MRI scans to identify meningioma disease presence and number of tumours. Previously operated meningiomas were included in numbers.

Age at initial diagnosis of meningioma was recorded.

Kaplan-Meier analysis was used to assess age at onset of cranial meningioma.

Mutation groups were defined as 1) exons 1-3, which encode the first FERM subdomain of the Merlin protein product; 2) exons 4-6, which encode the second FERM sub-domain; 3) exons 7-9 encoding the third FERM sub-domain; 4) exons 10-13, coding for the α-helical domain and 5) exons 14-15, which encode the C-terminal domain.

Large deletions removing the promoter and exon 1 and whole gene deletions were excluded from the analysis. Multi-exon deletions within the gene were grouped according to the first exon affected by the deletion. Numbers of patients with and without cranial meningiomas were compared using the chi-squared test.

Results

In our cohort of 411 NF2 patients, 215 were male and 196 were female. Sixty six individuals were mosaic. In total 199/411 (48%) NF2 patients had at least one cranial meningioma that was identified on MRI scan or at surgery.

Effect of mutation position

The position and prevalence of each mutation type identified in the NF2 gene and the proportion of individuals in each category with at least one cranial meningioma are presented in Table 1. As expected, the prevalence of cranial meningiomas was lower in individuals with a missense mutation (26%) or splicing mutation (35%), compared with truncating mutations (nonsense or frameshift, 58%). No significant difference was seen between patients with nonsense or frame shift mutations and those with a large deletion.

Mutations within the first two FERM sub-domains had the greatest impact in risk of cranial meningioma (61% and 58% with meningioma, respectively). Mutations occurring in the third FERM sub-domain or within the α-helical domain conferred a lesser risk of cranial meningioma (39% and 45%, respectively). The α-helical domain appears to confer a higher risk than the third FERM domain, although the difference is not significant and is due to a few high risk mutations in the second half of exon 10, at the junction between these two domains. Mutations occurring at junctions between other domains were also seen to have an increased risk (Figure 1b). Patients with mutations in the C-terminal domain were the least likely to develop a cranial meningioma (24%).

The Kaplan-Meier plot in Figure 1 shows a significant difference between the 5' and 3' mutations: exons 1-3, exons 4-6, exons 7-9 and exons 10-13 compared to exons 14-15 (p<0.0001). Cumulative risk of cranial meningioma to age 50 years was 70% for exons 1-3, 81% for exons 4-6, 49% for exons 7-9, 56% for exons 10-13 and 28% for exons 14-15. The hazard ratio (HR) for mutations in exons 7-15 was 0.499 (95% CI 0.376 to 0.663) compared to mutations in exons 1-6 (p<0.0001). A similar hazard ratio of 0.48 (0.27-0.81) was obtained by comparing truncating mutations in exons 14 and 15 with truncating mutations in exons 1-6 (Figure 2). The risk of cranial meningioma up to age 50 years was only 44% for truncating mutations in exons 14 and15 (similar to the risk for all mutation types in these exons) compared to 83% for truncating mutations in exons 1-6 (HR 0.48; p=0.0065). The average age of diagnosis for a cranial meningioma in patients with a mutation in exons 1-13 was 24.8 years, compared to exons 14 and 15 of 34.2 years.

Gender effect

Overall, no gender bias was observed between the number of males (105/215, 49%) and females (94/196, 48%) who developed at least one cranial meningioma (Figure 3a). By 50 years of age approximately 60% of both men and women had developed a cranial meningioma and by 60 years this increased to over 70%. There was no overall significant difference in age of cranial meningioma onset between males and females with the curves closely aligned after 30 years of age. However, these totals disguise an age related effect. There was a significantly greater risk of males developing a cranial meningioma below 20 years of age compared to females (p=0.023). The gender differences for age at onset were sustained across all mutation position groups with cranial meningiomas being detected earlier in males in the first 20 years of life.

Overall, a higher number of cranial meningiomas was found in females than males (p<0.0001). The 105 males with cranial meningiomas had a total of 223 cranial meningiomas (median 2; range 1-10; mean per patient 1.05), whereas the 94 females with cranial meningiomas had a total of 282 cranial meningiomas (median 2; range 1-10; mean per patient 1.42). However, this was mainly due to a small number of females with ≥ 4 cranial meningiomas, four of whom had 10 cranial meningiomas each, compared to only one male with 10 cranial meningiomas.

Effect of somatic mosaicism

We also analysed the effects of somatic mosaic mutations on risk of cranial meningioma formation. The rates of somatic mosaicism have been reported to be approximately 30% of de novo NF2 mutations in patients presenting with bilateral vestibular schwannomas and up to 60% in patients presenting with a de novo mutation and a unilateral vestibular schwannoma [START_REF] Evans | Mosaicism in neurofibromatosis type 2: an update of risk based on uni/bilaterality of vestibular schwannoma at presentation and sensitive mutation analysis including multiple ligation-dependent probe amplification[END_REF]. In our current cohort, 27% of patients with a new mutation were mosaic, comprising 23% of males with a new mutation and 32% of females with a new mutation.

In total, 36/64 (56%) of individuals with mosaic NF2 had at least one cranial meningioma and the mutation position effect can still be seen in somatic mosaic individuals. Of six individuals, mosaic for a mutation in exons 14-15, only one (aged 34) had a cranial meningioma. The ages of the six individuals without meningiomas ranged from 35-63 years.

The risk of cranial meningioma in mosaic patients was found to be greater in women 24/36 (67%) than men 12/28 (43%). Women with somatic mosaic disease were also at a greater risk of developing cranial meningiomas compared to women with nonmosaic NF2 (p=0.0475). This effect was not seen in men, but in women the effect correlated with a significantly higher proportion of mosaic females with truncating mutations and a lower number of mutations in exons 14 and 15 than non-mosaic women.

Discussion

Within the general population, meningiomas account for over 30% of all primary brain and central nervous system tumours (CBTRUS, 2010). NF2 disease is known to confer an increased risk of meningioma and intracranial meningiomas have been reported to occur in approximately 50% of NF2 patients. The molecular genetic factors that determine which patients will develop a meningioma are unclear.

In this study we analysed the germline mutations in a series of 411 patients with NF2, to determine whether the type and position of these mutations within the NF2 gene influence the likelihood of developing cranial meningiomas. As expected, we found that patients with a truncating mutation were more likely to have a cranial meningioma than those with a non-truncating mutation.

In addition to the effect of mutation type, we have demonstrated a clear positional effect of NF2 gene mutations on the risk of developing a cranial meningioma in NF2 disease, with mutations towards the 5' end of the gene being associated with a higher prevalence of cranial meningioma. The highest risk regions for each domain occur at the junctions between domains. We did find that NF2 patients with a constitutional mutation in exon 14 or 15 may develop a cranial meningioma, although the tumours occur much less frequently in these patients and develop at a later age (p<0.001). A similar positional effect on disease severity in NF2 has previously been described for splice site mutations [START_REF] Baser | The location of constitutional neurofibromatosis 2 (NF2) splice site mutations is associated with the severity of NF2[END_REF][START_REF] Kluwe | Phenotypic variability associated with 14 splice-site mutations in the NF2 gene[END_REF]. Although no formal statistical analysis was carried out, interestingly, no meningioma was found in any of the 13 exon 14/15 mutation-bearing individuals [START_REF] Kluwe | Phenotypic variability associated with 14 splice-site mutations in the NF2 gene[END_REF].

The positional effect was seen across all mutation types, with truncating mutations in exons 1-13 conferring a higher risk than those in exons 14 and 15. If all truncating mutations resulted in complete loss of Merlin protein as predicted, we would expect to see no positional effect for this type of mutation. Therefore, our findings provide evidence that at least some nonsense or frameshift mutations must generate mutant proteins, rather than abolishing the protein product through nonsense mediated decay.

In sporadically occurring meningiomas, women are approximately three times as likely to develop a meningioma as men (CBTRUS, 2010). A female preponderance has also been reported previously within NF2 patients [START_REF] Evans | Variation of expression of the gene for type 2 neurofibromatosis: absence of a gender effect on vestibular schwannomas, but confirmation of a preponderance of meningiomas in females[END_REF].

In contrast, several reports have shown that males are more likely to develop a meningioma within the first two decades of life [START_REF] Rushing | Central nervous system meningiomas in the first two decades of life: a clinicopathological analysis of 87 patients[END_REF]. Meningiomas are much less common in childhood than in adults and around 40% of childhood meningiomas occur within the context of NF2.

In our current cohort of 411 NF2 patients, we found no evidence for an overall increased risk of cranial meningioma in females. Overall, there was no difference in risk of cranial meningioma between males and females, although our data support the finding that males are more likely than females to develop cranial meningiomas before 20 years of age. It is possible that a female preponderance after 20 years of age is masked by the increased risk to males from the predominant tumour causing diagnosis of NF2.

Diagnosis with a symptomatic vestibular schwannoma in many NF2 patients gives rise to the simultaneous identification of often asymptomatic meningioma disease.

We did find that women were more likely to develop more than one cranial meningioma than men. In addition, there was an overall increased risk to women with somatic mosaic disease of developing a cranial meningioma than women with nonmosaic NF2.

In conclusion, mutations in exons 14 and 15 confer a reduced risk of cranial meningioma than mutations in exons 1-13. There is a similar risk of cranial meningioma for males and females with NF2 over their lifetime; however men under 20 years of age have a higher incidence of cranial meningioma. This information is of value when counselling patients about their risk of cranial meningioma. The bar graph (b) shows the overall risk of meningioma for mutations in individual exons. 5' regions confer a higher risk than 3' regions and exons encoding the junctions between functional domains confer a higher relative risk than the rest of the domain.

Figure 1 .

 1 Figure 1. The position of mutations in the NF2 gene affects the likelihood of

Figure 2 .

 2 Figure 2. Truncating mutations (nonsense and frameshift) alone also show a

Figure 3 .

 3 Figure 3. There is no overall difference in the likelihood of developing a meningioma

Table 1 .

 1 Effect of mutation type and position on risk of meningioma

	Exon	Proportion with	Missense	Nonsense or	Splice site	In frame deletion
	group	meningioma		frame shift		or large deletion
		Number of patients with	Number of patients with	Number of patients with	Number of patients with	Number of patients with
		meningioma/ total NF2	meningioma/ total NF2	meningioma/ total NF2	meningioma/ total NF2	meningioma/ total NF2
		patients with mutation in	patients with missense	patients with nonsense	patients with splice site	patients with in frame
		this exon group	mutation	or frameshift mutation	mutation	del or large deletion
						mutation

Acknowledgments

We would like to acknowledge John Dawson for his help with statistical analysis. The department of Genetic Medicine is supported by the NIHR Manchester Biomedical Research Centre. We would also like to thank all the patients who provided samples for the study.

Competing interests

The authors have no competing interests to declare.

Licence Statement

"The Corresponding Author has the right to grant on behalf of all authors and does grant on behalf of all authors, an exclusive licence (or non exclusive for government employees) on a worldwide basis to the BMJ Publishing Group Ltd to permit this article (if accepted) to be published in Journal of Medical Genetics and any other BMJPGL products and sublicences such use and exploit all subsidiary rights, as set out in our licence (http://group.bmj.com/products/journals/instructions-for-authors/licenceforms)."


Karen Grønskov

Rebecca L Poole

Johanne M D Hahnemann

Jennifer Thomson

Zeynep Tumer

Karen Brondum-Nielsen

Rinki Murphy

Kirstine Ravn

Linea Melchior

Alma Dedic

Zeynep Tümer

Karen Brøndum-Nielsen

Birgitte Dolmer

Karen Temple

Susanne E Boonen

Deborah Mackay

DELETIONS AND REARRANGEMENTS OF THE H19/IGF2 ENHANCER REGION IN PATIENTS WITH SILVER-RUSSELL SYNDROME AND GROWTH RETARDATION

Keywords: Silver-Russell syndrome, H19 enhancer, deletion, methylation

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

INTRODUCTION

Silver-Russell syndrome (SRS, MIM 180860) is a genetically and clinically heterogeneous condition characterized by intrauterine and postnatal growth retardation, dysmorphic facial features and body asymmetry [START_REF] Abu-Amero | The genetic aetiology of Silver-Russell syndrome[END_REF]. 35-60% of SRS cases show hypomethylation of the H19 imprinting control region (H19-ICR; Fig. 1a) at11p15, which is believed to result in reduced expression of IGF2 (MIM 147470) [START_REF] Eggermann | Epigenetic mutations in 11p15 in Silver-Russell syndrome are restricted to the telomeric imprinting domain[END_REF][START_REF] Gicquel | Epimutation of the telomeric imprinting center region on chromosome 11p15 in Silver-Russell syndrome[END_REF][START_REF] Netchine | 11p15 imprinting center region 1 loss of methylation is a common and specific cause of typical Russell-Silver syndrome: clinical scoring system and epigenetic-phenotypic correlations[END_REF]. The underlying causes of hypomethylation include mosaic maternal uniparental disomy (UPD) of chromosome 11 [START_REF] Bullman | Mosaic maternal uniparental disomy of chromosome 11 in a patient with Silver-Russell syndrome[END_REF], maternal duplication of 11p15 [START_REF] Eggermann | Is maternal duplication of 11p15 associated with Silver-Russell syndrome?[END_REF][START_REF] Eggermann | Chromosome 11p15 duplication in Silver-Russell syndrome due to a maternally inherited translocation t(11;15)[END_REF], and possibly gene mutations of H19-ICR [START_REF] Demars | Analysis of the IGF2/H19 imprinting control region uncovers new genetic defects, including mutations of OCT-binding sequences, in patients with 11p15 fetal growth disorders[END_REF], but is unknown in most cases. In most SRS patients hypomethylation occurs at both H19-ICR and IGF2P0, but is occasionally restricted to H19-ICR or IGF2P0 [START_REF] Bartholdi | Epigenetic mutations of the imprinted IGF2-H19 domain in Silver-Russell syndrome (SRS): results from a large cohort of patients with SRS and SRS-like phenotypes[END_REF].

The H19-ICR is located within the H19 promoter and is normally methylated on the paternally-derived allele. H19 (MIM 103280) encodes a maternally-expressed non-coding RNA while IGF2 encodes the paternally-expressed growth factor IGF-II. In mice the endodermal enhancers of IGF2 are located +10 kb from H19 transcription start site; skeletal muscle and mesodermal enhancers are located at +35 kb and other mesodermal enhancers at +120 kb [10]. The interaction of these enhancer sequences with the promoters of H19 or IGF2 is modulated by CTCF for which the H19-ICR contains seven potential binding sites [11]. CTCF binding to unmethylated H19-ICR prevents interaction between the IGF2 promoter and enhancers, while methylation of the H19-ICR prevents CTCF binding, enabling the enhancers to interact with the IGF2 promoter. Naturally occurring structural anomalies in the enhancer region of H19/IGF2 in patients are very informative in testing the hypothesis of complex imprinting and competing enhancer regulation of H19 and IGF2 and ultimately help to more clearly define the roles of these genes in humans. We report here for the first time SRS patients with chromosomal structural mutations in the H19/IGF2 enhancer region.

METHODS

Patients

Patient 1 was referred for SRS testing at the age of 2 years due to pre-and postnatal growth retardation, characteristic dysmorphic features, congenital anomalies and failure to thrive.

Patient 2 was referred at age 14 weeks with pre-and postnatal growth retardation, typical facial features and severe feeding difficulties. Subsequently his sister (patient 3) presented neonatally with the same clinical features. Patient 4 was identified after karyotyping due to recurrent miscarriages. She had pre-and postnatal growth retardation and atypical diabetes, but there was no clinical suspicion of SRS. Informed consent for the investigations performed and publication was obtained for all patients. Patients 1 and 4 were recruited via the initial diagnostic testing, while patients 2 and 3 were recruited into the research study "Imprinting Disorders Finding Out Why?" (IDFOW), approved by Southampton and South West Hampshire Research Ethics committee 07/H0502/85.

Molecular genetic investigations

Karyotyping was performed on cultured lymphocytes using standard techniques.

Flourescence-in-situ-hybridization (FISH) analysis was performed on metaphases by standard techniques using BAC and fosmids clones. Genomic DNA was extracted from blood leukocytes by standard procedures.

MLPA

MLPA analysis was performed by methylation-specific MLPA kit ME030 version B1 from MRC-Holland (Amsterdam, The Netherlands) following the manufacturer's instructions. 250 ng of DNA was used for analysis.

Pyrosequencing

Genomic DNA (2μg) was bisulphite treated using the EZ-DNA Methylation™ kit, according to manufacturer's instructions (Zymo Research, Orange, CA USA), except that DNA was eluted in 50μl of elution buffer and 50μl of TE pH 8.0.

Pyrosequencing assays of bisulphite-induced C/T polymorphisms [12] for H19 and IGF2P0 were designed using Biotage Assay Design Software™ v1.0.6 (Biotage, Uppsala, Sweden) and performed essentially as previously described [13]. Briefly, primary amplification reactions (50μl, performed in duplicate) contained forward and reverse primers, one of which was biotinylated, Titanium Taq and buffer (Clontech, Mountain View, Ca. USA), and ~20ng bisulphite-treated DNA. After 45 amplification cycles, 20μl single-stranded biotinylated PCR products were pyrosequenced using Pyrosequencing Gold Reagents (Biotage) according to the manufacturer's instructions. The percentage methylation at each CpG site was determined by AQ software. AQ levels >5% were indistinguishable from background fluorescence. In each pyrosequencing assay, a bisulphite treatment control was included in the form of a cytosine not in a CpG dinucleotide; this would be expected to be converted entirely to uracil.

Array CGH analysis

Array CGH analysis (aCGH) used a custom 8 X 60K Agilent array (Agilent Technologies, South Queensferry, UK) designed to densely cover known human imprinted loci. Probes were designed using Agilent Technologies' e-Array design tool with the following filters applied: Tm filter, similarity filter and catalog probe preference. Across the region chr11:1,674,000-2,979,000 (hg18) probes were tightly tiled to provide maximum resolution with an average spacing of 50nt. The remainder of chr11 was represented by Agilent aCGH control probes, which provide a genome-wide backbone to enable efficient normalization.

aCGH was performed according to the manufacturer's instructions. Briefly, 500ng of blood genomic DNA was labeled using the genomic DNA enzymatic labeling kit, hybridized to the array for 24 hours and subsequently washed. Arrays were scanned using an Agilent Technologies' DNA microarray scanner G2539A and data extracted using Agilent

RESULTS

The clinical features of the patients are summarized in Table 1. Cytogenetic analysis of patient 1 showed a normal karyotype, 46,XY while MLPA analysis showed deletion of H19 and part of H19-ICR (Fig. 1a); the peak height of the deleted probes suggested mosaicism, while MLPA analysis of parental DNA was normal, suggesting that the deletion had arisen de novo. Microsatellite analysis showed reduced dosage ratios of D11S4046 (paternal:maternal ratio 1:2.3) consistent with mosaicism (Fig. 1c). Cytogenetic analysis of patients 2 and 3 revealed a pericentric inversion of chromosome 11 (karyotypes 46,XY,inv(11)(p15.5q21)pat and 46,XX,inv(11)(p15.5q21)pat, respectively), inherited in both children from their phenotypically normal father. FISH analysis showed a split signal using BAC clone RP11-534I22, while fosmids W12-1875E5 and W12-1660M14 hybridized to 11p and 11q, respectively (Fig. 1a and Fig. 1d) mapping the 11p15.5 breakpoint to chr11:1,862,216-1,921,555 (hg18).

Patient 4 had a balanced translocation (karyotype 46,XX,t(1;11)(p36.22;p15.5)). FISH analysis mapped the 11p15.5 breakpoint between W12-1875E5 and W12-1660M14 [14].

Paternal origin of the translocation was inferred through karyotyping of the patient's mother, who had a normal karyotype, and the patient's daughter, who had the same translocation as the patient but normal growth.

DNA methylation analysis showed partial hypomethylation of both H19-ICR and IGF2P0 in patient 1. In patients 2, 3 and 4 IGF2P0 was partially hypomethylated while the methylation of H19-ICR was in the normal range (Fig. 1b).

All four patients were analyzed by array CGH (aCGH). In patient 1 aCGH showed a 58.8 kb deletion (chr11:1,918,222-1,977,026 (hg18)) encompassing the endodermal, skeletal muscle specific and part of the mesodermal enhancers, as well as H19 and part of H19-ICR including the seventh CTCF binding site, but not the two repeat blocks (Fig. 1a), though the log intensity difference suggested a mosaic deletion. In patients 2 and 3 aCGH showed an 8.8 kb deletion (chr11:1,918,312-1,927,132 (hg18)) overlapping with the mesodermal enhancer.

Deletions partially overlapping with the above deletions have been registered as rare variations in apparently normal individuals (http://projects.tcag.ca/variation/?source=hg18) but the parental origin of these deletions are unknown. aCGH did not reveal any deletion in patient 4.

DISCUSSION

We present three pedigrees with disruption of the shared H19/IGF2 enhancer region and separation of the mesodermal enhancers from IGF2, one by deletion, one by an inversion breakpoint accompanied by deletion, and one by breakpoint of an apparently balanced translocation. IGF2 expression was reduced in mesodermal tissue in patient 4 [14].

Previously, the small stature of patient 4 was unaccounted for, but we now hypothesize that hypomethylation of IGF2-P0 due to the separation of IGF2 from its enhancers causes growth retardation. The deletion in patient 1 may truncate these enhancers and moreover reduces their distance from IGF2 by 50%, thus potentially affecting the local three-dimensional structure of the DNA and hence its function. We therefore suppose that the mesodermal enhancers are non-functional in all four patients. It is interesting that all four patients are growth-retarded, but only patients 1, 2 and 3, who harbor deletions including the sequence chr11:1,918,312-1,927,132 (hg18), presented with the clinical diagnosis of SRS. One interpretation of this observation is that reduced IGF2 expression from mesodermal enhancers is responsible for only the growth retardation seen in SRS, while other genetic or epigenetic defects cause the other phenotypic traits of SRS. Alternatively, given that SRS is a disorder with marked clinical heterogeneity it may be that patient 4 showed an atypical presentation.

It was striking that all patients showed hypomethylation of IGF2P0, while hypomethylation of the H19-ICR was observed only in patient. This observation suggests that IGF2P0 methylation may not directly depend on H19-ICR methylation in cis, but may instead reflect the promoter activity of IGF2P0 per se, and be established as a consequence of its interaction with distant IGF2 enhancers. Disruption of the enhancer region represents a possible explanation of rare SRS cases with hypomethylation restricted to IGF2P0 [START_REF] Bartholdi | Epigenetic mutations of the imprinted IGF2-H19 domain in Silver-Russell syndrome (SRS): results from a large cohort of patients with SRS and SRS-like phenotypes[END_REF].

The effects of deletions and mutations in the H19 region depend not only on their parent of origin, but also on their localization. Deletions or mutations of H19-ICR described to date are apparently without phenotypic effect when on the paternal allele, whereas in maternal inheritance they are associated with the overgrowth disorder Beckwith-Wiedemann syndrome (MIM 130650) [15,16]. Thus, we hypothesize that a paternal deletion of H19, such as seen in patient 1, may be of less clinical significance than the separation of IGF2 from its enhancers.

For the enhancer region, no data have hitherto been presented in humans; however, mice with paternally inherited deletions were 80% of normal size whereas mice with maternal deletions were of normal size [17]. We predict that detailed characterization of the IGF2 promoter and enhancer elements in humans will uncover further genetic and epigenetic changes that predispose to growth retardation and/or SRS. Moreover, we confirm that IGF2P0 hypomethylation is found in a proportion of SRS cases with normal H19-ICR methylation.

Therefore, in agreement with Bartholdi et al. [START_REF] Bartholdi | Epigenetic mutations of the imprinted IGF2-H19 domain in Silver-Russell syndrome (SRS): results from a large cohort of patients with SRS and SRS-like phenotypes[END_REF] we suggest that IGF2P0 methylation analysis should be included in standard molecular testing for SRS.

(10) Arney KL. H19 and Igf2--enhancing the confusion?

  Primer sequences were: H19 pyrosequencing: forward GTATAGTATATGGGTATTTTTGGAGG, reverse CCATAAATATCCTATTCCCAAATAACC (bio); sequencing GTTTYGGGTTATTTAAGTT. This amplifies the region chr11: 1977625-1977877. IGF2(P0) pyrosequencing: forward TGAGGATGGGTTTTTGTTTGGTAT; reverse TCCTCAATCCACCCAAAATAATAT (bio); sequencing AAAAGTTATTGGATATATAGT.

  Technologies' Feature Extraction software version 10.5.1. Duplications and deletions were mapped and end points determined by Agilent Technologies' DNA Analytics software version 4.0 with the following settings: data was analyzed using the ADM-2 Algorithm with a threshold of 6.0 using the DLR error model, with fuzzy zero turned on and the continuous moving average applied. Web resources Online Mendelian Inheritance in Man (OMIM), http://www.ncbi.nlm.nih.gov/Omim/ UCSC Human Genome Database Build hg18, March 2006, http://www.genome.ucsc.edu Database of Genomic Variants (DGV): http://projects.tcag.ca/variation/?source=hg18

Figure 1 :

 1 Figure 1: genetic organization of the human H19-IGF2 region, and cytogenetic, genetic and

  indicate the proportions of methylated and unmethylated product. c. Electropherogram of microsatellite D11S4046. Middle panel: patient 1; upper panel: father; lower panel: mother. d. FISH analysis of patient 3. Upper panel: W12-1875E5 (red) retained on 11p; lower panel: W12-1660M14 (red) which is transferred to 11q. The centromere probe D11Z2 (green) hybridizes to both normal and derivative chromosomes 11.

ACKNOWLEDGMENTS

MLPA was performed by Judy Boe Rasmussen, The Kennedy Center; FISH was performed by Viv Maloney, Wessex Genetics Service. Jette Bune Rasmussen performed the graphic work.

Café au lait patches Yes No No No Yes No No No Yes No No No No No No No

COMPETING INTERESTS

The authors declare no competing interests