domain: spi phys math info

text:
Pierre Vuillemin
email: pierre.vuillemin@onera.fr

Adrien Maillard
email: adrien.maillard@jpl.nasa.gov

Charles Poussot-Vassal
email: charles.poussot-vassal@onera.fr

Optimal Modal Truncation

Keywords: Model approximation, modal truncation, mixed optimisation

This paper revisits the modal truncation from an optimisation point of view. In particular, the concept of dominant poles is formulated with respect to different systems norms as the solution of the associated optimal modal truncation problem. The latter is reformulated as an equivalent convex integer or mixed-integer program. Numerical examples highlight the concept and optimisation approach.

Introduction

Large-scale dynamical models often arise in the industry due to the inherent complexity of the systems or phenomena to be studied and the complexity induced by the processes and tools used for their modelling (e.g. Finite Element Methods, etc.). The dimension of these dynamical models then translates into a high numerical and computational burden that can prevent from performing simulation, analysis, control or optimisation. Model approximation is meant to alleviate the issue by building a much smaller model catching the main dynamics of the initial one and that could be used instead. In this context, this article is aimed at improving a standard linear model approximation method, the modal truncation, by adding considerations based on usual systems norms minimisation.

Let us consider consider a Linear Time-Invariant (LTI) dynamical model H represented by its state-space realisation, ẋ(t) = Ax(t) + Bu(t) y(t) = Cx(t) + Du(t)

where u(t) ∈ R nu is the control input, x(t) ∈ R n the internal state, y(t) ∈ R ny the output and A, B, C, D are real matrices of adequate dimensions. Generally speaking, the objective of model approximation consists in finding a LTI model Ĥ described by ẋ(t) = Âx(t) + Bu(t)

ŷ(t) = Ĉ x(t) + Du(t) (2)
where the number of input and outputs remains unchanged but the dimension of the state is decreased, i.e. x(t) ∈ R r with r n, and such that the input to output behaviour of Ĥ is close to the one of H in some sense.

Submitted. Preprint version September 2020. Funding. Part of this research was carried out at the Jet Propulsion Laboratory, California Institute of Technology, under a contract with the National Aeronautics and Space Administration (80NM0018D0004).

For sake of simplicity the state-space representations (1) and ( 2) of H and Ĥ are used indistinctly from their transfer functions, defined as follows,

H : C \ ρ(A) → C ny×nu s → H(s) = C(sI -A) -1 B + D (3)
where ρ(A) = {λ ∈ C/det(λI -A) = 0} is the set that contains the eigenvalues of A.

Linear model approximation has been widely studied and several methods are now available. See e.g. [START_REF] Antoulas | Approximation of large-scale dynamical systems[END_REF] for an overview of classical approaches and [START_REF] Benner | Model reduction and approximation: theory and algorithms[END_REF][START_REF] Antoulas | Interpolatory methods for model reduction[END_REF] for a more recent treatment of the topic. Among the classical approaches, the modal truncation consists in projecting the large-scale model H onto its dominant eigenspace. While the method is not the most efficient in general, it remains widely used in practice by engineers due to its conceptual simplicity and the fact that it preserves some poles of the initial model which are quantities of physical interest. Besides, it may be efficient enough to produce faithful reduced-order models, especially with poorly damped systems such as flexible structures.

In this context, the objective of this paper is to revisit the modal truncation from an optimisation point of view. In particular, usual systems norms are used to define specific dominant poles so that the resulting reduced-model is the best among all the possible models obtained by modal truncation with respect to the chosen norm. It is shown how this problem translates naturally into (binary) Integer Programming (IP) or Mixed Integer Programming (MIP). As the resulting problems share some similarities, various norms are considered: the H 2 norm, its time and frequency-limited versions and the H ∞ norm.

The mandatory concepts and tools in systems theory and approximation are recalled in section 2, especially elements concerning the norms of systems that are considered in the article and about the modal approximation algorithm. Then in section 3, the latter is reformulated into different optimisation problems depending on the considered norm. Numerical applications are presented in section 4 to highlight the concepts of the approach. Finally, concluding remarks are exposed in section 5 together with some insights into possible extensions of this work.

Note that throughout this article, the following hypothesis are considered:

• The model H is assumed to be stable, i.e. ρ(A) is contained within the open left complex plane C -. Indeed, should it have unstable poles, then they should be kept in the reduced-order model Ĥ anyway to have a bounded input-output approximation error.

• The standard modal truncation approach requires the full eigenvalue decomposition of the matrix A which involves O(n 3 ) dense linear algebra operations. Therefore the dimension n of the state should remain moderate in practice. For iterative dominant eigenspaces computation, see [START_REF] Rommes | Efficient Computation of Multivariable Transfer Function Dominant Poles Using Subspace Acceleration[END_REF] and references therein.

• Similarly, one assumes that A has only semi-simple eigenvalues. While the approach could theoretically be extended in presence of Jordan blocks, their computation is ill-conditioned and could hardly be achieved in practice for non-trivial cases.

It should also be noted that a matrix E ∈ R n×n multiplying ẋ in (1) could be considered without affecting the reminder of the article. Indeed, should the matrix E be non-singular, then it could be inverted to fall back on a system of the form [START_REF] Antoulas | Approximation of large-scale dynamical systems[END_REF]. Otherwise, the associated transfer function would contain a polynomial part which elements of order larger or equal to 1 should be kept for the truncation error to be finite, similarly to the unstable case.

Notations Let us denote by R the set of real numbers, R + the subset of positive numbers, C the set of complex numbers, S n ++ the set of symmetric and positive definite matrices of size n, L 2 (I) the Lebesgues space of square integrable functions on I (R if not specified). Given a complex valued matrix M , M T denotes its transpose, M * its conjugate and M H its conjugate transpose, M F is its Frobenius norm, [M ] i,j is its i, j-th element, tr(M ) is its trace, vec(M ) is the vector obtained by concatenating the columns of M . Considering z ∈ C, Re(z) and Im(z) are its real and imaginary parts, respectively, ı = √ -1 is the imaginary unit. Floor and ceiling functions are denoted by • and • , respectively.

Preliminaries

Diagonal canonical form

This Diagonal Canonical Form (DCF), given in proposition 1, appears naturally throughout the modal approximation and allows to simplify the expression of usual systems norms as recalled in the next section. See remark 1 for its practical computation.

Proposition 1 (Diagonal Canonical Form). Let us denote by I n = {i = 1, . . . , n} the set of integers ranging from 1 to n. Provided that the dynamic matrix A of H in (1) has semi-simple eigenvalues λ i ∈ C (i ∈ I n ), then the transfer function associated with the system can be decomposed as follows,

H(s) = D + i∈I n Φ i s -λ i , (4)
where Φ i ∈ C ny×nu is the residue associated with the pole λ i .

Remark 1. From a practical point of view, the decomposition (4) is obtained by computing the full eigenvalue decomposition, AX = X∆

where ∆ = diag(λ 1 , . . . , λ n ) and X ∈ C n×n contains the right eigenvectors of A. Then by applying the change of variable x(t) = Xξ(t), the dynamic (1) becomes,

ξ(t) = X -1 AX ∆ ξ(t) + X -1 B B ∆ u(t), y(t) = CX C ∆ ξ(t) + Du(t), (6)
and the residues are thus given, for i ∈ I n , as

Φ i = C ∆ e i e T i B ∆ = c i b T i . (7)
Note that the residues are written as the outer product of two vectors c i ∈ C ny and b i ∈ C nu corresponding to the columns (resp. rows) of C ∆ (resp. B ∆ ) thus ensuring that rank(Φ i ) = 1 and that the right hand side of (4) is indeed of order n.

Systems norms

This section is aimed at recalling some key elements concerning norms of LTI systems that are useful for this article. For a more exhaustive introduction to the topic, interested readers may refer to [15, chap.4] and references therein.

In particular, the H 2 -norm, its frequency and time limited counterparts are given in definition 1, definition 2 and definition 3, respectively. Provided the considered model is in DCF as in equation (4), these norms then have a simplified expressions as detailed in proposition 2, proposition 3 and proposition 4. In addition, the definition of the H ∞ -norm is recalled in definition 4.

Definition 1 (H 2 -norm). Considering a stable and strictly proper LTI model H as in equation (1), its H 2 -norm is defined in the frequency domain as,

H 2 1 2π ∞ -∞ H(ıν) 2 F dν. (8)
Proposition 2. Assume that the stable and strictly proper LTI model H is in DCF as in equation (4), then its H 2 -norm can be computed as follows,

H 2 2 = i∈I n tr(Φ i H(-λ i ) T ) = i,k∈I n tr(Φ i Φ T k ) -λ i -λ k . (9)
The H 2 -norm can be related to the time-domain in various ways. In particular, in the context of model reduction, let us consider the approximation error model E = H -Ĥ, then for any input signal of bounded energy u ∈ L 2 (R), the worst-case output error between the two models y -ŷ ∞ is upper bounded as follows

y -ŷ ∞ ≤ H -Ĥ 2 u 2 . ( 10
)
Definition 2 (Frequency-limited H 2 -norm). The frequency-limited H 2 -norm [START_REF] Gawronski | Advanced structural dynamics and active control of structures[END_REF], denoted H 2,ω -norm, is defined by restricting the frequency interval to [-ω, ω] in (8) so that

H 2,ω 1 2π ω -ω H(ıν) 2 F dν. ( 11
)
Proposition 3. The frequency-limited H 2 -norm of a LTI model H in DCF form can be be computed similarly to (9) as follows,

H 2 2,ω = ω π D 2 F - 2 π i∈I n tr(ΦiH(-λi) T )atan ω λi , (12)
where atan(z) = 1 2j (log(1 + ız) -log(1 -ız)) and log(z) is the principal value of the complex logarithm for z = 0.

Proof. See [START_REF] Vuillemin | Poles residues descent algorithm for optimal frequency-limited H 2 model approximation[END_REF].

Note that H 2,ω is only a semi-norm when considering the whole Lebesgue space L 2 (ıR) but it is a norm for rational functions as considered here.

Let h(t) denotes the impulse response of H, corresponding to the inverse Laplace transform of H(s), i.e.

h(t) L -1 (H)(t) = Ce At B + Dδ(t), (13)
where δ is the Dirac impulse. Due to Parseval's equality, for a stable and strictly proper model H, the H 2 -norm can also be computed in time-domain as follows,

H 2 = h 2 ∞ -∞ h(t) 2 F dt. ( 14
)
For a stable system H, h(t) = 0 for t < 0, therefore the integral in ( 14) can be restricted to R + . In addition, by restricting even further the integration interval to [0, τ ], one can define the time-limited H 2 -norm [START_REF] Gawronski | Advanced structural dynamics and active control of structures[END_REF][START_REF] Goyal | Time-limited H 2 -optimal model order reduction[END_REF] as detailed in definition 3.

Definition 3 (h 2,τ -norm). Considering a stable and strictly proper LTI model H as in (1) with impulse response h(t), its time-limited H 2 -norm, denoted h 2,τ -norm here, is defined for τ > 0 as

h 2 2,τ τ 0 h(t) 2 F dt. ( 15
)
Proposition 4. The time-limited h 2 -norm of a LTI model H in DCF form can be computed similarly to (9) as follows,

h 2 2,τ = i,k∈I n tr(Φ i Φ T k ) λ i + λ k e (λi+λ k )τ -1 . (16)
Proof. The DCF (4) with D = 0 enables to re-write the impulse response as

h(t) = i∈I n Φ i e λit , (17)
which naturally leads to the expression (16) after integration.

Again, h 2,τ is only a semi-norm for the whole space of square integrable functions L 2 (R), but it is a norm for the impulse response functions associated with rational functions.

A time-domain bound of the error similar to the one available with the H 2 -norm ( 10) can be derived [START_REF] Goyal | Time-limited H 2 -optimal model order reduction[END_REF].

Definition 4 (H ∞ -norm). Considering a stable LTI model H as in (1), its H ∞ -norm is defined as follows, H ∞ sup ν∈R σ 1 (H(ıν)), ( 18
)
where σ 1 (H(ıν)) is the largest singular value of the transfer matrix.

The H ∞ -norm represents the worst amplification gain of the system and is a widely used measure of robustness. Similarly to the H 2 -norm, within the context of model reduction, it enables to bound the L 2 gain of the approximation error,

y -ŷ 2 ≤ H -Ĥ ∞ u 2 . ( 19
)
Note that the computation of the H ∞ -norm either requires an iterative bisection procedure or the resolution of a Semi Definite Program (SDP) [START_REF] Scherer | Multiobjective Output-Feedback Control via LMI Optimization[END_REF].

Reminder on the modal truncation

Modal truncation consists in keeping only r elements from the decomposition (4) to form the reducedorder model. In particular, by defining the subset I n r ⊂ I n containing r unique elements from

I n , Ĥ(s) = D + r i=1 Φi s -λi = D + i∈I n r Φ i s -λ i . (20)
Note that for Ĥ to have a real realisation, the retained complex eigenvalues of H must be selected together with their complex-conjugate pair. This implies that some combinations are not allowed to form I n r . In standard modal truncation, D is chosen equal to D. Modal truncation then boils down to select the dominant poles-residues couples that should be included within I n r . Dominant poles may be defined in various ways. Below, the definition based on the bound of the H ∞ -norm of the approximation error that is generally considered in the reduction literature is recalled.

Approximation error and bounds

Let us denote by E = I n \ I n r the set of discarded indexes. The approximation error E between H and Ĥ is then naturally given as,

E(s) = H(s) -Ĥ(s) = D -D + i∈E Φ i s -λ i . ( 21
)
For its H 2 norm to be bounded, its direct feedthrough must be zero, i.e. D -D must be zero. In that case, E 2 is readily obtained by considering the specific formulation of the norm for transfer functions with such structure ( 9),

E 2 2 = i∈E tr(Φ i E(-λ i ) T ). ( 22
)
For the H ∞ -norm, only an upper bound has been derived based on the triangular inequality. For D = D, it states that,

E ∞ ≤ i∈E Φ i 2 |Re(λ i )| . ( 23
)
Usual criterion for determining I n r Dominant poles are generally defined as the poles λ i that have the largest ratio

Φ i 2 |Re(λ i )| . ( 24
)
Such a choice to fill the set I n r enables to minimise the H ∞ bound (23) of the approximation error. Still, it does not make the resulting reduced-order model optimal with respect to the H ∞ -norm and we shall see in section 3 that optimality considerations allow to characterise dominant poles in a more generic way.

Optimal modal truncation

The main idea here consists in formulating the modal truncation method as an optimisation problem. The latter is stated formally in problem 1. Dominant poles-residues are then defined in definition 5 as the elements associated to the corresponding optimal solution. Problem 1 (Optimal modal truncation). Let us consider a stable n-th order LTI dynamical model H in DCF (4), a reduction order 0 < r < n and the set I n r of indexes containing the elements to be kept within reduced-order model Ĥ as in (20). Considering in addition some system norm • (e.g. H 2 , H ∞ , etc.), the optimal modal truncation problem can then be formally stated as min

I n r H -Ĥ s.t.
Ĥ given by (20) Ĥ has real coefficients (25)

Definition 5 (Dominant poles-residues). Suppose that În r solves the optimal modal approximation problem (25), then the set of r dominant poles-residues with respect to the associated system norm is defined as

Λ( În r ) = {λ i , Φ i } i∈ În r (26)
As illustrated in example 1, solving problem 1 boils down to select the r poles-residues amongst n that minimise the error. Note that as complex poles must come by pairs, the exact number of possible unique combinations depends on their number within the initial model H as detailed in proposition 5.

Proposition 5. Considering a n-th order LTI model H with n c pairs of complex poles and n r = n-2n c real poles, then the number of possible combinations for r-th order modal truncation (r < n) with a real state-space realisation is 0

H 1 H 2 H 3 H 2 H 1 H 3 H 3 H 1 H 2
κ(n r , n c , r) =      nr r if n c = 0 nc Ne if n r = 0 Ne i=Ns nc i nr r-2i if n r > 0, n c > 0 ( 27
)
with N s = max(0, r-nr

2

) and N e = min(n c , r 2 ). Proof. Case 1 (n c = 0) : if there is no complex pole in the model, n = n r , the binomial coefficient counts how many sets of size r formed of n r real poles exist. Case 2 (n r = 0) : if there is no real pole in the model, n = n c and we have to count as previously. As the set of poles in an r-modal truncation is of size r and the complex poles come by pairs, the quotient N resulting from the division of r by 2 represents the maximum number of complex poles in an r-modal truncation. Case 3 (n c > 0, n r > 0) : in between N s and N e , the formula uses the same principle and sums over the different configurations that the set of poles can take when there is a combination of real and complex poles. N s and N e are set so binomial coefficients are always defined. When i = N s , there is no or the minimum number (depending the number of real poles) of complex pole in the truncation. When i = N e and r is even, the set of poles is only made of complex poles (r -2i = 0) and the first term is equal to the second case. When i = N e and r is odd, the set contains only N e complex poles and one real pole. The right term equals nr 1 = n r which effectively counts the number of possibilities to fill the single slots in the sets of complex poles.

While the number of combinations κ grows slower than n r as soon as the model H has some complex eigenvalues, it is still large enough to prevent from scanning exhaustively the decision tree.

Example 1. Let us consider the following third order model

H(s) = H1(s) + H2(s) + H3(s) = 1 s + 1 + 1 s + 3 + 2 s + 5 , (28)
that should be reduced to r = 2. In this simple case and considering e.g. the H 2 -norm, one can evaluate all the possible combinations. figure 1 shows a tree where each node represents the subsystem which is added in the reduced-order model. For instance, the leftmost branch leads to

H 1 + H 2 .
The approximation error is also displayed next to each node. Note that this tree contains redundant branches (in grey) due to commutativity of the sum. As only real poles are considered, at the end, there is n r = κ(3, 0, 2) = 3 unique possible combinations and the optimal model Ĥ is clearly given by

Ĥ(s) = H 1 (s) + H 3 (s), (29)
meaning that {-1, 1} and {-5, 2} are the dominant poles-residues for the H 2 -norm.

To reformulate problem 1 in a practical way, let us consider the following parametrization of the reduced-order model,

H α (s) = D α + i∈I n α i Φ i s -λ i , (30)
where α i = e T i α ∈ {0, 1} are binary variables acting as activation variables and D α ∈ R ny×nu . With this parametrization, the order of the reduced model may be enforced by the constraint

1 T α = r, (31)
where 1 ∈ R n is a vector full of ones. In addition, to ensure that H α has a real realisation, the complex conjugate pairs of poles must be kept together. This translates into additional linear constraints between the α i . Indeed, let us consider the set of n c complex conjugate pairs indexes,

I C = {{i, j} ∈ I n /λ i ∈ C, Im(λ i ) > 0, λ i = λj }, (32)
and the matrix M ∈ R nc×n which contains a row for each couple {i, j} ∈ I C such that,

[M ] •,i = -[M ] •,j = 1. ( 33
)
Then the realness constraint is

M α = 0. ( 34
)
By coupling the parametrization (30) with the constraints (31), (34) and the binary constraint, problem 1 may be reformulated as a binary optimisation problem as stated in problem 2.

Problem 2 (Binary formulation of optimal modal approximation). Considering the parametrization (30) for the reduced-order model, problem 1 is equivalent to the following problem,

min α H -H α s.t. α ∈ {0, 1} n 1 T α = r M α = 0 (35)
In the following section, problem 2 is specified for H 2 -norm, its frequency/time limited counterparts and the H ∞ -norm.

In H 2 -norm

Considering the framework introduced in problem 2, the optimal H 2 modal truncation problem can be recasted as a convex binary quadratic problem as stated in theorem 1.

Theorem 1 (Optimal H 2 modal truncation). Considering the notations of problem 2, the r-th order optimal H 2 modal truncation model H α is such that D α = D and α is the solution of the following convex quadratic binary problem,

min α (1 -α) T Q(1 -α) s.t. α ∈ {0, 1} n 1 T α = r M α = 0 (36)
where Q ∈ C n×n is a hermitian matrix which entries are given as,

[Q] i,j = tr(Φ i Φ H j ) -λ i -λ * j , i, j = 1, . . . , n (37)
Proof. For the H 2 -norm of the approximation error E α = H -H α to be finite, D α is constrained to be equal to D and may be discarded in the sequel. In that case, E α is given as

E α (s) = H(s) -H α (s) = i∈I n (1 -α i ) Φ i s -λ i , (38)
and from equation ( 9), the quadratic nature of the error w.r.t. α appears,

E α 2 2 = i,k∈I n (1 -α i ) tr(Φ i Φ T k ) -λ i -λ k (1 -α k ). ( 39
)
As each complex pole-residue pair in H comes with its complex conjugate, the sums in the right-hand side of (39) may be reordered so that {λ k , Φ T k } is replaced by their conjugate. The H 2 error can then be rewritten as

E α 2 2 = (1 -α) T Q(1 -α). (40)
Coupling the objective function (40) with the constraints (31), (34) and the binary constraint leads to the optimisation problem (36). Its objective is a squared-norm and is therefore strictly convex. The equality constraints of problem (36) are linear and thus convex. Therefore, it is a binary convex quadratic problem.

As the relaxation of the binary problem (36) is convex, efficient branch and bounds algorithms (see e.g. [START_REF] Bliek1ú | Solving mixed-integer quadratic programming problems with IBM-CPLEX: a progress report[END_REF]) can be used to solve the overall optimal H 2 truncation problem.

Yet, state of the art solvers may not handle the fact that Q is complex. However, as each complex element comes with its conjugate in the sum (40) (the overall sum is real), Q can be replaced with Q = 1 2 (Q + Q H ) which leads to the same objective function as long as M α = 0.

About the initialisation While existing general purpose solvers are perfectly able to determine a feasible starting point, providing a meaningful feasible initial solution may help to prune rapidly some parts of the tree.

As highlighted in proposition 6, the H 2 -norm of the approximation error is upper bounded by the sum of each subsystems norms. Therefore, this suggests to initially select the pole-residue pairs with largest associated H 2 -norm (42). Again, complex conjugate pairs must be kept together meaning that corner cases have to be dealt with. In particular, if one pole remains to be selected but the next largest couple is complex, then either discard it until a real pole is reached or discard the last real pole selected to get the complex pair. Alternatively in those cases, decrease or increase r by one. The initialisation process is highlighted in example 2. Proposition 6. The H 2 -norm of the approximation error E α is bounded by the sum of the discarded subsystems norms,

E α 2 ≤ i∈E H i 2 , (41)
where the individuals subsystems norms are given as

H i 2 = Φ i F -2Re(λ i ) . ( 42
)
Proof. Applying the triangular inequality to the norm of the approximation error E α directly leads to the result.

Example 2. Considering the trivial model of example 1, the H 2 criterion (42) indicates,

H 1 2 = 1/ √ 2 ≈ 0.71, H 2 2 = 1/ √ 6 ≈ 0.41, H 3 2 = 2/ √ 10 ≈ 0.63. ( 43
)
This suggests H 1 + H 3 as reduced model, which is indeed the optimal one in that case. Note that the H ∞ criterion (24) leads to the same ordering here but this could be otherwise for multiple inputs multiple outputs cases.

In H 2,ω -norm

Similarly to the H 2 -case, the optimal H 2,ω modal truncation problem can be recasted as an optimisation problem as stated in Theorem 2. The main difference lies in the fact that for ω < ∞, the frequency-limited H 2 -norm remains finite even when there is a direct feedthrough. Therefore, unlike in the H 2 case, D α ∈ R ny×nu remains a free design variable in the parametrisation (30) and the resulting optimisation problem is thus a mixed convex quadratic program.

Theorem 2 (Optimal H 2,ω modal truncation). Considering the framework of problem 2, let us define

ᾱ = 1 -α ∈ R n , m = n + n y n u , d = vec(D) ∈ R nynu , d α = vec(D α ) ∈ R nynu , a ω (λ i ) = atan( ω λi ) and x α = [ᾱ T , d T α ] T ∈ R m .
Then, the r-th order optimal H 2,ω modal truncation model H α is obtained through x α , solution of the following strictly convex mixed quadratic problem,

min α,dα x T α Q ω x α + c T ω x α s.t. α ∈ {0, 1} n d α ∈ R nynu 1 T α = r M α = 0 (44)
where Q ω ∈ C m×m is a hermitian matrix defined as

Q ω = Q ω 1 2 U ω 1 2 U H ω ω π I nynu ( 45
)
where the top left block is related to the matrix Q (37), for i, j = 1, . . . , n

[Q ω ] i,j = - 1 π (a ω (λ i ) + a ω (λ * k )) [Q] i,j (46)
and the off-diagonal term is given, for i = 1, . . . , n, as

e T i U ω = - 2 π a ω (λ i )vec(Φ i ) T . ( 47
)
Additionally, the linear term is given by

c ω = - 2 π a ω (λ 1 )vec(Φ 1 ) T • • • a ω (λ n )vec(Φ n ) T ωd T T . ( 48
)
Proof. Based on equation ( 12), the norm of the approximation error can be written as three components,

E α 2 2,ω = E 1,ω + E 2,ω + E 3,ω = ω π D e 2 F . . . - 2 π i∈I n (1 -α i )tr(Φ i D T e )atan( ω λ i ) . . . - 2 π i,k∈I n (1 -α i ) tr(Φ i Φ T k )atan( ω λi ) -λ i -λ k (1 -α k ), (49)
where

D e = D -D α .
Similarly to the H 2 case, the error (49) is also quadratic with respect to the optimisation variables. Indeed, first, note that the terms within the sums of E 3,ω can be rearranged so that tr(Φ

i Φ T k )atan( ω λi ) is replaced by tr(Φ i Φ H k ) 1 2 (atan( ω λi ) + atan( ω λ * k
)). Therefore, E 3,ω is similar to (40) and can be written as

E 3,ω = (1 -α) T Q ω (1 -α), (50)
where Q ω is given by equation ( 46). Then, by using vectorisation, E 1,ω can be transformed as follows,

E 1,ω = ω π tr(D T D) -2tr(D T D α ) + tr(D T α D α ) = ω π d T d -2d T d α + d T α d α .
(51)

Finally, using the same vectorisation process, the last element divides into a quadratic part and a linear part,

E 2,ω = (1 -α) T (U ω d α + f ω ), ( 52
)
where, for i = 1, . . . , n, e

T i f ω = -2 π tr(Φ i D T )atan( ω λi ), e T i U ω = -2 π atan( ω λi )vec(Φ i ) T . ( 53
)
By considering ᾱ and stacking it with d α , the final structure of the approximation error appears,

E 2 2,ω = ᾱT d T α T Q ω 1 2 U ω 1 2 U T ω ω π I nynu ᾱ d α + f T ω -2 ω π d T ᾱ d α + ω π d T d (54)
The constant part can be discarded and the constraints remain the same as in problem 2 with the additional optimisation variables D α ∈ R ny×nu . As • 2,ω is a norm for rational functions, the objective is strictly convex making the overall optimisation problem a (strictly) convex mixed quadratic program.

About the initialisation Similarly to the H 2 case, the frequency-limited norm of the approximation error can be upper bounded as shown in proposition 7. Subsystems with highest individual norm may be selected initially. Note that unlike the H 2 -norm case which is parameter free, here, some subsystems may become more relevant depending on the value of the frequency bound ω as illustrated in example 3.

Proposition 7. The H 2,ω -norm of the truncation error E α is bounded by the sum of the discarded subsystems norms,

E α 2,ω ≤ i∈E H i 2,ω , (55)
where the individuals subsystems norms are given as

H i 2,ω = Φ i F -2Re(λ i ) - 2 π Re(atan( ω λ i )). ( 56
)
Proof. The triangular inequality applied to the approximation error leads to the bound. The H 2,ωnorm of each individual system is obtained by expanding H(-λ i ) in [START_REF] Rommes | Efficient Computation of Multivariable Transfer Function Dominant Poles Using Subspace Acceleration[END_REF] and reordering the terms in the sums to pair complex conjugate components. To reduce the system to an order 2, there is only κ(2, 1, 2) = 2 solutions in the decision tree. Each coincides either with H 1 or with H 2 . The H 2,ω -norm of the approximation errors are computed for varying values of the frequency bound ω ranging from 10 -2 to 10 3 and are reported in figure 2 (top) together with the value of the heuristic criterion (56) associated with each mode (bottom). One can see that H 1 is dominant (the error is lower) for low values of ω while H 2 becomes dominant after 100 rad/s. Besides, as shown by the bottom figure in that simple case, the sorting criterion (56) is coherent with the optimal result. An uncertainty area appears just before 100 rad/s where the criteria for the complex eigenvalues crosses one of the real pole but not the other one.

In h 2,τ -norm

As the H 2 case, the direct feedthrough is here constrained to be equal to D so that the resulting optimal h 2,τ modal truncation problem reduces to a convex binary quadratic problem as stated in theorem 3.

Theorem 3 (Optimal h 2,τ modal truncation). Considering the notations of problem 2, the r-th order optimal h 2,τ modal truncation model H α is such that D α = D and α is the solution of the following convex quadratic binary problem,

min α (1 -α) T Q τ (1 -α) s.t. α ∈ {0, 1} n 1 T α = r M α = 0 (58)
where Q τ ∈ C n×n is a hermitian matrix which entries are given, for i, j = 1, . . . , n, as,

[Q τ ] i,j = (1 -e (λi+λ * j )τ )[Q] i,j . (59)
Proof. As in the H 2 case, the direct feedthrough D α of the reduced-order model H α must be equal to D. The error is then the same as in equation ( 38) which, combined with the poles-residues expression of the h 2,τ -norm (16) leads to

eα 2 2,τ = i,k∈I n (1 -αi) tr(ΦiΦ T k ) λi + λ k (e (λ i +λ k )τ -1)(1 -α k ), (60)
where e α (t) = L -1 (E α )(t). Again, the norm of the approximation error exhibits a quadratic structure and may be reformulated as

e α 2 2,τ = (1 -α) T Q τ (1 -α), (61)
where the entries of the matrix Q τ are given, after reordering of the elements, by (59).

About the initialisation In proposition 8, an upper bound of the approximation error that can motivate the selection of the initial poles-residues is presented.

Proposition 8. The h 2,τ -norm of the approximation error L -1 (E α ) = e α is bounded by the sum of the discarded subsystems norms,

e α 2,τ ≤ i∈E h i 2,τ , (62)
where the individuals subsystems norms are given as

h i 2,τ = Φ i F -2Re(λ i ) 1 -e 2Re(λi)τ . (63)

In H ∞ -norm

The structure of the H ∞ -norm makes the associated optimal modal truncation problem quite different from the previous ones. Indeed, considering problem 2 with the H ∞ -norm leads to a convex, albeit non-smooth, mixed optimisation problem. While a general purpose convex optimizer may be used to solve the convex relaxation in a branch and bound process, the specific of the H ∞ -norm can be exploited to derive an alternative equivalent problem which may be more direct to address. The latter in detailed in theorem 4.

Theorem 4 (Optimal H ∞ modal truncation). Let us consider the complex diagonal form of the approximation error E α ,

A e = ∆, B e = B ∆ , C e = C ∆ ∆ α , D e = D -D α , (64)
with ∆ α = I n -diag(α 1 , . . . , α n ). Then the r-th order optimal H ∞ modal truncation model H α is given by the solution of the following mixed SDP,

min γ,P,α,Dα γ s.t. α ∈ {0, 1} n D α ∈ R ny×nu γ ∈ R + P ∈ S n ++ 1 T α = r M α = 0 M e (α) ≺ 0 (65)
where the last constraint is a Linear Matrix Inequality (LMI) characterised by the matrix

Proof. As its objective is convex, by adding a slack variable γ ≥ 0, the H ∞ modal approximation problem can be re-written equivalently (see [START_REF] Boyd | Convex Optimization[END_REF]) as

min γ,α,Dα γ s.t. α ∈ {0, 1} n D α ∈ R ny×nu γ ∈ R + 1 T α = r M α = 0 E α ∞ ≤ γ (67)
The last constraint can then be transformed even further. Indeed, let (∆, B ∆ , C ∆ , D) denote the complex realisation associated with the DCF of H so that H(s) = C ∆ (sI n -∆) -1 B ∆ + D. Equation (64) then represents a complex realisation associated with the approximation error E α = H -H α .

By using the Bounded Real Lemma (see e.g. [START_REF] Scherer | Multiobjective Output-Feedback Control via LMI Optimization[END_REF]), the constraints E α ∞ ≤ γ can be traded by considering additional slack variables as a symmetric and positive definite matrix P ∈ S n ++ such that the following matrix inequality is satisfied,

The inequality can further be transformed by a Schur complement leading to the LMI M e (α) ≺ 0 where M e (α) given by equation (66).

The convex relaxation of the binary constraints in (65) leads to a SDP which can be easily formulated using a modelling framework such as YALMIP [START_REF] Löfberg | Yalmip : A toolbox for modeling and optimization in matlab[END_REF] and solved by associated solvers such as SeDuMi [START_REF] Polik | SeDuMi: a package for conic optimization[END_REF]. Still, while a SDP is convex, it remains difficult to solve, especially as the dimension of the problem increases. For these reasons, the optimal H ∞ modal truncation is restricted to models of moderate dimension n.

The LMI matrix (66) is complex-valued. Should the SDP solver only handle real matrices, then an equivalent real-valued matrix should be considered as detailed in proposition 9.

Proposition 9. Let us consider the set of indexes associated with real poles I R = {i ∈ I n , λ i ∈ R} and the set I C defined in (32). Let us define the unitary transformation matrix T ∈ C n×n such that

T i,i = 1, i ∈ I R , (69)
and such that for {i, j} ∈ I C ,

T i,i = T j,i = 1/ √ 2, T i,j = -T j,j = ı/ √ 2. ( 70
)
Then, assuming that the constraint M α = 0 is satisfied, the complex realisation (64) can be replaced by the following real realisation,

Ãe = T H A e T, Be = T H B e , Ce = C ∆ T ∆ α . (71)
Proof. The real matrices Ãe , Be and C ∆ T correspond to the standard real-valued realisation associated with the DCF obtained by combining complex conjugate elements. The diagonal matrix ∆ α acts as a filter on the output matrix which may be applied before or after the transformation to real form provided that the complex pairs are kept together. Indeed, as long as this is the case, i.e. that the constraint (34) is satisfied, the matrices ∆ α and T commutes. More specifically, as ∆ α is idempotent, ∆ α T = T ∆ α is equivalent to ∆ α T (I -∆ α ) = 0. Looking at each entry leads to α i T i,j (1 -α j ) = 0 for i, j ∈ I n which is satisfied considering the zero elements of T and the dependencies between the α i for complex poles. Therefore, Ce = C ∆ ∆ α T = C ∆ T ∆ α which leads to the result.

Numerical illustrations

In subsection 4.1, an academic example is used to highlight the different nature of poles dominance depending on the considered norm. Then, a more realistic example is considered in subsection 4.2 for H 2 dominant modes identification.

Academic example

To highlight the differences between the different norms within the context of modal truncation, let us consider the following 8-th order academic model,

H(s) = H 1 (s) + H 2 (s) + H 3 (s) + H 4 (s), (72)
where the subsystems are given as follows,

H 1 (s) = 25 s 2 /10
The impulse responses of H and its components are plotted in figure 3. The gains of their frequency responses are plotted in figure 4. To reduce this model to an order r = 2 by modal truncation, 1 with ω = 0.1 and T = 0.2. Note that no direct feedthrough is considered here.

In this simple case, the results could, in a sense, be inferred from the gain diagram in figure 4:

• H 1 has visually an important mean contribution and should therefore play an important role in the H 2 -norm.

• H 2 has a large contribution only in low frequency below 0.1rad/s and should therefore be important w.r.t. the H 2,ω -norm.

• H 3 contains the fastest modes and should therefore be dominant at the beginning of the impulse response thus dominating the h 2,T -norm when T = 0.2. This is validated by the impulse responses in figure 3.

• H 4 has the highest gain and should therefore be dominant w.r.t the H ∞ -norm.

However, as illustrated in the next example, such a heuristic analysis is no longer tractable on realistic models and the systematic approach detailed in this work then shows its benefits.

Application to dominant modes selection

In this example, let us consider the ISS2 model from [START_REF] Leibfritz | Description of the benchmark examples in COMP l e ib 1.0[END_REF]. It is a 270-th order model with 3 inputs and outputs. It has only complex eigenvalues and the number of combinations for modal approximation grows as κ(0, 135, r) where r should be chosen even. Assuming that we are looking for the r = 10 dominant poles, then there are more than 10 8 possible combinations.

The H 2 optimal modal truncation problem is solved with a simple branch and bound algorithm starting from the initialisation scheme suggested in Section 3.1. The set of initial poles ρ(A) are plotted together with the r H 2 dominant ones in figure 5. The singular values associated with the initial and reduced-order transfer functions are plotted in figure 6.

The initialisation heuristic choice turns out to be almost optimal in that case as only 3 nodes lead to an improvement of the H 2 error. This is not surprising since the model represents a highly flexible structure which dominant modes are, in a sense, easily distinguishable due to the large magnitude of some of the associated residues. Indeed, as shown by figure 6, dominant modes are mainly associated with peaks in the frequency-domain response. However, this may not be as clear in general, especially when the model has some well damped dynamics.

It is clear from figure 5 that dominant poles are not necessarily the ones with lowest natural frequency and that the associated residues play and important role. Indeed, the dominance must be understood in and input-output sense through the lens of the actuators, sensors and eigenvectors. The pole dominance as defined here is therefore associated with and input-output setting. This article revisits the well-known modal truncation technique from an optimisation point of view.

In particular, dominant poles are defined as the solution of the associated optimal modal truncation problem with respect to different systems norms. The latter reduces to a convex integer or mixed integer program depending on the considered norm: H 2 -norm and its frequency/time limited variants, H ∞ -norm.

The same approach may be developed for other norms or quantities of interest. For instance, one may be interested in knowing which poles are dominant in response to some specific input other than an impulse such as a step, a combination of sine, etc. This can give meaningful insights to determine the main dynamics in an automatic way for control design purposes.

In this article, the choice has been made to stick to the original modal truncation framework, i.e. the poles and the associated residues are kept fixed (excepted the term D). Therefore, the resulting model is only optimal among all the models with same residues and poles. To improve its matching with the initial large-scale model, the residues (one side in the MIMO case) should be considered as free variables. To avoid the multiplication of the binary variables with the residues in the objective, the optimisation problem may then be modified with the big-M technique such that each binary variable behaves as an activation variable through the constraint Φ i F ≤ α i M where M is large enough. The problems then remain convex and may be solved similarly.

More generally, the underlying concept of this article consists in determining the dominant components of the additive decomposition of a LTI model. While the diagonal canonical form is considered here, the idea may be generalised to other decompositions. An interesting candidate is the block diagonal form (see e.g. the discussion in [START_REF] Moler | Nineteen dubious ways to compute the exponential of a matrix[END_REF]). Indeed, the latter preserve the Jordan blocks and can thus be safely applied to systems with multiple eigenvalues. However this induces changes for the computation of the norms that must be integrated to the associated optimisation problems.

Figure 2 :

 2 Figure 2: H 2,ω -norm of the approximation error (top) and heuristic sorting criteria associated with each eigenvalue (bottom).

M

  e (α) =   A T e P + P A e P B e C T e

  A T e P + P A e + C T e C e P B e + C T e D e D T e D e -γ 2 I nu ≺ 0.

Figure 3 :

 3 Figure 3: Impulse responses of H and its components with the time axis in log scale.

Figure 5 :Figure 6 :

 56 Figure 5: Initial set of poles ρ(A) and subset Λ containing the 10 H 2 dominant poles.

  Example 3. Let us consider the 4-th order model H = H 1 + H 2 with

	H 1 (s) =	2.2 s + 0.1	+	1.2 s + 0.2
				(57)
	H 2 (s) =	1.2 s 2 /10 4 + 0.02s/100 + 1

Table 1 :

 1 Norms of the approximation errors H -H i , i = 1, . . . , 4 only κ(2, 3, 2) = 4 combinations are possibles. These combinations are actually associated with each subsystem H i .The different norms of the approximation errors are computed and reported in table

	100
	50
	0
	-50
	-100
	-150
	-200
	-250
	10 -2	10 -1	10 0	10 1	10 2	10 3	10 4
		Figure 4: Gain of the frequency responses of H and its components.


Marc Lahaye
email: marc.lahaye@inrae.fr

Tiphaine Gauclain

Jade Bakiepa

Xavier Falourd

Jean-Michel Le-Quere

Rémi Bauduin

O1.6 -Role of apple flesh mechanical properties in juice production

Keywords: Cell wall, cellulose, pectin, DMA, NMR

à la diffusion de documents scientifiques 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.

As turgor pressure is the major contributor to apple flesh firmness, mechanical properties were measured by Dynamic Mechanical Analysis (DMA) on fresh and plasmolyzed tissue. Groups of varieties were identified on the basis of the impact of plasmolysis on the residual viscoelastic properties (storage modulus E' and tan ). Insights in the cell wall structure and organization was obtained by solid state 13 C NMR spectroscopy using a variable contact time CPMAS (VCT-CPMAS) experiment and modeling of signals intensity with regard to magnetization transfer time.

The T1 relaxation time of the proton attached to C4 of cellulose taken as a probe of cellulose organization showed a clear correlation with E' of fresh apple flesh which indicated that tightly packed cellulose fibrils in the cell wall are unfavorable for flesh firmness. The magnetization transfer time from neighboring proton (mainly from water) to pectin methyl ester protons also demonstrated a correlation with E' indicating that pectin hydration plays a major role on fruit viscoelastic properties. Both DMA and NMR results are in support of a poro-viscoelastic type mechanical properties of apple flesh in which water in the porous cell wall network contribute to the fruit firmness [1]. Relation between cider apple juice yield obtained by an experimental press with viscoelastic mechanical properties of fresh and plasmolyzed fruit collected at different levels of maturity showed that relation could be obtained only from fresh fruit E'.

If there is a global positive trend between fruit firmness and juice yield, better relations were obtained according to varieties. Thus, mechanical properties of apple flesh represent a weak indicator of juice yield and other dominant genetically dependent structural factors are likely involved.

This study was supported in part by the interregional Bretagne, Normandie and Pays de la Loire research program Optipress-2. Part of the work was realized on the UR 1268 BIA BIBS instrumental platform (IBiSA, Phenome-Emphasis-FR ANR-11-INBS-0012).


Natural Deep Eutectic solvents pre-treatments ease pectin water extraction from apple pomace

Mingrui Chen, Marc Lahaye

CC:Glycerol (G) were tested alone or in sequence (first. CC:G and second. CC:LA) as a pre-treatment of apple pomace to ease water extraction of pectin. Extraction yields, sugar composition, degree of methyl and acetyl esterification and molecular weight distribution were characterized to evaluate the impact of these pre-treatments. The process was adjusted to limit cell wall polysaccharide degradation and losses allowing for the later recovery of other valuable polymers (hemicellulose, cellulose) from the residues.

Results show that combining CC:LA and CC:G pretreatment improves extraction yield (21.1, pomace dw%) compared to classical methods (the highest yield 17.5%). Moreover, CC:G enhanced the later CC:LA pretreatment ability by 4.5%. Overall, the sequential NADES pretreatment facilitated extraction of a large amount of high methyl esterified pectin (65.7% of pomace uronic acids). This sequential pretreatment lead to the recovery of pectin with higher rhamnose/uronic acid ratio compared to classical extracts which indicated a higher content in rhamnogalacturonan I structural domains.

Thus, NADES pretreatment succeeded in increasing the extraction yield of pectin by hot water extraction and offer a "green" alternative to classical harsh chemical extraction. The impact of the chemical characteristics of pectin enriched in rhamnogalacturonan domains remains to be evaluated with regard to their functional properties.

Keywords: Biorefinery; Cell wall polysaccharides; Valorization of food waste Acknowledgements: MC is beneficiary of a PhD scholarship from the China Scholarship Council

3 -Food systems and sustainability of F&V processing O3.3 -Natural Deep Eutectic solvents pre-treatments ease pectin water extraction from apple pomace

	Mingrui CHEN, Marc LAHAYE*
	INRAE, UR BIA, F-44316, Nantes, France
	*Correspondence: marc.lahaye@inrae.fr
	As a widely produced fruit in the world, apple once processed into a variety of products generates a
	large amount of pomace that is a traditional source of pectin. Industrial pectin extraction processes
	classically use harsh chemicals that pose environmental concerns. Recently, natural deep eutectic
	solvents (NADES) formed by hydrogen bonding acceptor and hydrogen bonding donor metabolites
	emerged as innovative solvents for extracting diverse molecules from biomass and were explored for
	extracting polysaccharides from food related wastes.

In the present study, as an alternative to classical extraction (water extraction, mineral acid extraction, sequential extraction with chelating agent), two NADES: Choline chloride (CC):Lactic acid (LA),


Song Tao

Chen Changyou

Wang Pingping
email: wangpp@mail.iee.ac.cn.

Long-Fei Wu

Biocompatible analysis of marine magnetotactic ovoid strain MO-1 for in vivo application *

Keywords: magnetotactic bacteria, motility, cell interaction, median lethal dose

Magnetotactic bacteria are capable of biosynthesizing magnetic nanoparticles, also called magnetosomes, and swimming along magnetic field lines. The abilities endow the whole cells of magnetotactic bacteria with such applications as targeted therapy and manipulation of microrobots. We have shown that the intact marine magnetotactic bacteria MO-1 kill efficiently antibiotic-resistant pathogen Staphylococcus aureus in vivo, but the biocompatibility of this marine bacteria is unknown. In this study, the strain MO-1 was chosen to analyze its biocompatibility and potential for biomedicine applications. Results showed that MO-1 cells could be guided at 37 °C under an external magnetic field and swim in the blood plasma and urine. They could keep active locomotivity within 40 minutes in the plasma and urine, although their velocity slowed down. When incubated with human cells, magnetotactic bacteria MO-1 had no obvious effects on cellular viability at low dose, while the cell toxicity increased with the augmentation of the quantity of the MO-1 cells added. In the in vivo experiments, the median lethal dose of magnetotactic bacteria MO-1 in rats was determined to be 7.9 × 10 10 cells/kg. These results provided the foundation for the biocompatibility and safety evaluations of magnetotactic bacteria MO-1 and suggested that they could be basically used in clinical targeted therapy.

INTRODUCTION

Magnetotactic bacteria are a group of phylogenetically diverse and Gram-negative microorganisms with the ability of synthesizing magnetite (Fe 3 O 4 ) or greigite (Fe 3 S 4 ) nanoparticles in a genetically controlled manner (R.P. [START_REF] Blakemore | magnetotactic bacteria[END_REF][START_REF] Bazylinski | Magnetosome formation in prokaryotes[END_REF][START_REF] Simmons | South-seeking magnetotactic bacteria in the Northern Hemisphere[END_REF]Faivre and Schüler, 2008). These nanoparticles, also called magnetosomes, are arranged in a chain-like manner, which helps magnetotactic bacteria to swim along magnetic lines [START_REF] Lefèvre | Evolution of the bacterial organelle responsible for magnetotaxis[END_REF]. Magnetotactic bacteria also have the flagella that enable them to swim autonomously and efficiently [START_REF] Zhang | Flagella and Swimming Behavior of Marine Magnetotactic Bacteria[END_REF]. On the basis of different types of swimming behavior, magnetotactic bacteria can be divided into two main categories: one is polar strain that swims persistently in one direction along the magnetic field, and the other one is axial strain that swims in either direction along magnetic field lines [START_REF] Bazylinski | Magnetosome formation in prokaryotes[END_REF]. The features of magnetosomes and autonomous movements along the magnetic field endow them with rapid responsibility to external magnetic field.

As a biosynthetic magnetic nanoparticle, magnetosome is coated by a biomembrane, which could be easily decorated for in vitro uses, such as DNA extraction and immunoassays [START_REF] Matsunaga | Biotechnological application of nano-scale engineered bacterial magnetic particles[END_REF]. Magnetosomes could also be used in vivo as a contrast agent [START_REF] Hu | Comparison of the 1H NMR Relaxation Enhancement Produced by Bacterial Magnetosomes and Synthetic Iron Oxide Nanoparticles for Potential Use as MR Molecular Probes[END_REF][START_REF] Heinke | MPS and MRI efficacy of magnetosomes from wild-type and mutant bacterial strains[END_REF], drug carriers [START_REF] Sun | In vitro and in vivo antitumor effects of doxorubicin loaded with bacterial magnetosomes (DBMs) on H22 cells: the magnetic bio-nanoparticles as drug carriers[END_REF][START_REF] Alphandéry | Chains of magnetosomes with controlled endotoxin release and partial tumor occupation induce full destruction of intracranial U87-Luc glioma in mice under the application of an alternating magnetic field[END_REF], and a tumor photo-/magneto-thermal transmitter [START_REF] Alphandéry | Chains of Magnetosomes Extracted from AMB-1 Magnetotactic Bacteria for Application in Alternative Magnetic Field Cancer Therapy[END_REF][START_REF] Alphandéry | Preparation of chains of magnetosomes, isolated from Magnetospirillum magneticum strain AMB-1 magnetotactic bacteria, yielding efficient treatment of tumors using magnetic hyperthermia[END_REF]Chen et al., 2016). However, they are not capable of autonomous swimming and do not target the deep area in the body even when an external magnetic field is applied, thereby greatly limiting their applications.

Applications of the whole cell of magnetotactic bacteria have recently attracted considerable attention. The whole magnetotactic bacteria could not only be guided to target the deep area but also perform as energy converters, which are conducive to microrobot construction and targeted therapy [START_REF] Fdez-Gubieda | Magnetotactic bacteria for cancer therapy[END_REF]. Magnetotactic bacteria are suitable actuators of bacterial microrobot due to their autonomous swimming ability. Magnetosopirrillum gryphiswalense (MSR-1), an axial strain, was constructed as a guidable delivery vehicle for cargo transport and release [START_REF] Alsaiari | Magnetotactic bacterial cages as safe and smart gene delivery vehicles[END_REF]. It was also integrated with ciprofloxacin-loaded microtubes to build controllable microrobots, and the robots could target an infectious biofilm under the guidance of an external magnetic field and then release the antibiotics [START_REF] Stanton | Magnetotactic Bacteria Powered Biohybrids Target E. coli Biofilms[END_REF]. A polar magnetotactic bacterium may be more suitable as an actuator than the above axial strain because it swims persistently in one direction along the magnetic field. A polar marine magnetotactic coccus, designated strain MC-1 is one of the best actuators for building microrobots. MC-1 bacteria pushed a microsphere and swam along the magnetic field, thus showing a great potential in manipulated movement under large weight loading [START_REF] Martel | Controlled manipulation and actuation of micro-objects with magnetotactic bacteria[END_REF]. Further, MC-1 cells were used to transport drug-loaded nanoliposomes to target the tumor in the blood under the control of magnetic field [START_REF] Taherkhani | Covalent Binding of Nanoliposomes to the Surface of Magnetotactic Bacteria for the Synthesis of Self-Propelled Therapeutic Agents[END_REF][START_REF] Felfoul | Magneto-aerotactic bacteria deliver drug-containing nanoliposomes to tumour hypoxic regions[END_REF]. Another polar magnetotactic bacterium, magnetotactic ovoid strain MO-1, is also a highly suitable candidate for constructing microrobots. Magnetotactic ovoid strain MO-1 was isolated from the Mediterranean Sea, and it is a kind of chemoautotrophic bacterium [START_REF] Lefèvre | Isolation and characterization of a magnetotactic bacterial culture from the Mediterranean Sea[END_REF]. MO-1-derived microrobots were constructed by attaching microbeads to the cell surface, and the microrobots could swim in the microfluidic channel under the control of magnetic field [START_REF] Ma | Construction and operation of a microrobot based on magnetotactic bacteria in a microfluidic chip[END_REF]. MO-1 cells have also been fabricated as a separating microrobot to target separate pathogens and then kill them using different external magnetic fields [START_REF] Chen | Construction of a microrobot system using magnetotactic bacteria for the separation of Staphylococcus aureus[END_REF]Chen et al., 2017). By using magnetosomes in the cell body, intact magnetotactic bacteria could produce positive magnetic resonance imaging contrast for visualizing mouse tumor xenografts, thereby providing a potential tool for improved cancer imaging visualization in preclinical and translational studies [START_REF] Benoit | Visualizing implanted tumors in mice with magnetic resonance imaging using magnetotactic bacteria[END_REF]. The whole cell of magnetotactic bacteria could also transmit the magnetic field energy into the thermal energy through the intracellular nanoparticles. When applied in the alternating magnetic field, they produced heat and then strongly affected cancer cell proliferation on clinical conditions [START_REF] Gandia | Unlocking the Potential of Magnetotactic Bacteria as Magnetic Hyperthermia Agents[END_REF]. Magnetotactic bacterium-based magnetic hyperthermia was also used to target and kill pathogenic cells, thus could potentially treat infectious diseases (Chen et al., 2016). Aside from thermal energy, magnetotactic bacteria could transmit the magnetic field energy to the mechanical force, which could inhibit the viability of pathogen (Chen et al., 2017).

Although the whole magnetotactic bacteria show many potentials in the construction of microrobots and targeted therapy, their biocompatibility and biosecurity should be considered seriously. On the one hand, when magnetotactic bacteria enter into the body fluid, such as blood and urine, for drug delivery and targeted treatment, their status and swimming ability require sufficient evaluation. When incubated in the blood, the terminal velocity of magnetococcus marinus strain MC-1 decayed for a duration of 40 minutes [START_REF] Martel | Flagellated Magnetotactic Bacteria as Controlled MRI-trackable Propulsion and Steering Systems for Medical Nanorobots Operating in the Human Microvasculature[END_REF]. On the other hand, magnetotactic bacteria, as foreign substances, are certainly attacked by the immune system in the human body, and the toxicity should be assessed. As for magneto-ovoid strain MO-1, which has been used for the construction of bacterial microrobots and targeted therapy, the biocompatible studies in the blood and urine are still unknown until now, and they need to be carried out.

In the present study, magneto-ovoid strain MO-1 was used to investigate its biocompatibility for in vivo applications. The swimming ability of MO-1 bacteria in the blood and urine was first assessed. The interaction of bacteria and human cells was examined to show the biocompatibility and the cellular safe dose. Finally, the median lethal dose (LD 50 ) of magnetotactic bacteria in rats was determined for in vivo biomedical utilization.

MATERIAL AND METHOD

Material

Chemical regents containing NaCl, MgCl 2 , Na 2 SO 4 , CaCl 2 , KCl, NH 4 Cl, NaHCO 3 , K 2 HPO 4 , and NaNO 3 were purchased from BBI Co., Ltd. Succinic acid, L-cysteine, tryntone, and yeast extract were bought from AMRESCO Commercial Finance, LLC. MTT was purchased from Sigma-Aldrich. The mineral solution and vitamin solution for magnetotactic ovoid strain MO-1 were prepared in accordance with previous reports [START_REF] Frankel | Magneto-aerotaxis in marine coccoid bacteria[END_REF]. Ferric quinate was also prepared. Artificial urine was purchased from Qingdao Jisskang Biotechnology Co., Ltd. (China), and its pH was adjusted to 7.5 by using NaOH solution.

Animal

Adult male Sprague-Dawley (SD) rats weighing 180-200 g were used in the LD 50 experiments and purchased from Beijing HFK Bio-Technology Co., Ltd (China). Blood was acquired after the rats were anesthetized, and then the plasma was obtained by centrifuging the rat blood.

Bacterial culture and transmission electron microscopy (TEM)

The magnetotactic ovoid strain MO-1 was grown at room temperature in EMS2 medium as previously reported (Lefevre, et al., 2009). The EMS2 medium was first added with (per L of deionized water; pH 6.8) 19.45 g NaCl, 5.9 g MgCl 2 , 3.24 g Na 2 SO 4 , 1.8 g CaCl 2 , 0.55 g KCl, 0.2 g NH 4 Cl, 1.26 g NaHCO 3 , 0.6 g L-cysteine, 3 g HEPES, and 0.35 g agarose. The medium was then autoclaved at 120°C for 20 min. After the medium was cooled to approximately 45 °C, the following solutions were added (per L) to previously sterile-filtered stock solutions in the following order: 0.5 mL vitamin solution, 0.5 mL ferric quinate (0.4 mM), 3 mL sodium thiosulfate (40% w/v), 4 mL dipotassium hydrogen phosphate (1.86% w/v), and 5 mL modified Wolfe's mineral solution. Ultrathin copper grids were immersed in the bacterial suspensions to observe the MO-1 cell morphology, and then the grids were dried in a clean bench after retrieval for observation using a transmission electron microscope (Hitachi H-7650B; Hitachi, Japan).

Speed analysis of magnetotactic bacteria in plasma and urine

Speed analysis of magnetotactic bacteria MO-1 in the plasma and urine was carried out in a microfluidic chip. The microfluidic channel contains four corners, and its plane geometry was shown in the previous report [START_REF] Ma | Construction and operation of a microrobot based on magnetotactic bacteria in a microfluidic chip[END_REF]. The width and depth of the channel were 200 μm and 20 μm, respectively. The rat plasma or artificial urine was pushed into a microfluidic channel.

After one of the micropools in the chip was blocked, the MO-1 solution was injected into the other micropool. A magnetic field was applied to guide the MO-1 cells in the microfluidic channel. The movement of MO-1 swarm in the plasma or artificial urine was observed under 40× magnification in the dark-field mode. The swimming behavior was recorded using a Canon 500D camera that recorded 1920 × 1080 pixel images at 20 frames per second. The experiments were conducted at 37 °C to mimic the in vivo biomedical applications of human physiological temperature and repeated thrice. The temperature controlling instrument of Live Cell Imaging System was used to keep the microfluidic chip with MO-1 cell samples at 37 °C. All the blood and urine solutions were at human physiological pH. The average swimming velocities were calculated using the MTrackJ plugin for ImageJ [START_REF] Murat | Opposite and Coordinated Rotation of Amphitrichous Flagella Governs Oriented Swimming and Reversals in a Magnetotactic Spirillum[END_REF][START_REF] Zhang | Measurement of Free-Swimming Motility and Magnetotactic Behavior of Magnetococcus massalia Strain MO-1[END_REF].

Cell culture and MTT assay

Human immortalized epidermal cell HaCaT (No. 3111C0001CCC000373) and human cervical carcinoma cell Hela (No. 3111C0001CCC000011; China Infrastructure of Cell Line Resources, Beijing, China) were maintained in Dulbecco's modified Eagle's medium (DMEM, Gibco/Life Technologies, Grand Island, NY, USA) supplemented with 10% fetal bovine serum (Australian origin, Hyclone Laboratories, Logan, UT, USA), 100 IU/mL penicillin (Sigma-Aldrich, St. Louis, MO, USA), and 100 μg/mL streptomycin (Sigma-Aldrich) at 37 °C in 5% CO 2 . The HaCaT and Hela cells were seeded at 10,000 cells/well in 100 μL DMEM in a 96-well plate and incubated for 24 h. The medium was then replaced with the magnetotactic bacteria MO-1 working solution for the succeeding experiments. The magnetotactic bacteria MO-1 working solution (per well in 100 μL DMEM) contains 0 (control sample), 2 × 10 5 , 4 × 10 5 , 8 × 10 5 , 1.6 × 10 6 , 3.3 × 10 6 , 6.5 × 10 6 , 1.3 × 10 7 , 2.5 × 10 7 , and 5 × 10 7 bacteria.

After incubation was performed for 24 h, the viability of the HaCaT and Hela cells was analyzed by MTT method. In brief, the cell medium was replaced with 100 μL DMEM without phenol red (Gibco) containing MTT agent and incubated at 37°C in the dark for 4 h. The medium was discarded, and 100 μL of dimethyl sulfoxide (Sigma-Aldrich) was added to the well. After shaking for 15 min, the cell plates were centrifuged at 3000 rpm for 5 min, and the supernatant was transferred to a new 96-well plate. The absorbance was measured at 570 nm with a reference at 630 nm by using a microplate absorbance reader (168-1130 model, Bio-Rad Laboratories, Hercules, CA, USA). The relative cell viability (%) relative to control was calculated using the following equation: [A(test sample)/A(control sample)] × 100% (where A is absorbance). The experiments were repeated thrice.

LD 50 of magnetotactic bacteria in rats

The magnetotactic ovoid strain MO-1 in exponential growth was collected using a magnetic block and suspended in phosphate buffer solution (PBS, pH7.0). Then, 16 SD rats were randomly assigned to four groups (n= 4), and 5 × 10 8 , 5 × 10 9 , 5 × 10 10 and 5 × 10 11 MO-1 cells in 3 mL PBS were injected into four groups of rats intravenously. The animals in each group were observed for 1 week. Afterwards, Behrens-Kärber method was used to calculate the LD 50 and 95% confidence interval (CI) [START_REF] Zhang | Median effective effect-site concentration of sufentanil for wake-up test in adolescents undergoing surgery: a randomized trial[END_REF]:

  1 50 lg 0.5 m LD X i p        , 50 lg LD pq Si n   , 95% CI:     50

  

where X m is the logarithm of the maximum dose, i is the difference in the logarithm of adjacent concentrations, p and q are the lethal success and failure rates in each group, and n is the number of rats in each group.

RESULT AND DISCUSSION

Characterization of magnetotactic bacteria MO-1

Magneto-ovoid strain MO-1 is a marine species affiliated with zeta proteobacterial strains, and it grew at aerobic-anaerobic interfaces in the culture flasks. The TEM results showed that magnetotactic bacteria MO-1 owned two flagellar bundles on the long hemisphere of the cell (Fig. 1a, black arrows). They have an autonomous swimming ability toward north pole in a magnetic field. TEM images showed that they also synthesized bio-membrane-coated magnetosomes with an average number of 11.1 ± 2.94 (Fig. 1b), which assisted them to feel and move along the magnetic fields. The morphologies of most MO-1 magnetosome particles have an elongated cubo-octahedral shape, with the length of 77.06 ± 11.1 nm, width of 71.54 ± 10.59 nm respectively (n=149). The features of fast autonomous movement and magnetotaxis may lead them to reach the deep site inside the organism body for drug delivery and targeted therapy.

Speed of magnetotactic bacteria MO-1 in rat plasma and artificial urine

Movement ability is an important characteristic for MO-1 cells to perform a task. Generally, in vivo application scenarios in biomedical areas are always accompanied with the body liquid environment, such as blood and urine. Here, the magnetotactic velocity of magnetotactic bacteria MO-1 in the blood plasma and artificial urine was assessed at the human physiological temperature of 37 °C to mimic the in vivo biological application scenarios. The rat plasma was acquired, and the MO-1 cells and the plasma were injected into the microfluidic chip. At 37 °C, when a magnetic field with a strength of 2.2 mT was applied, MO-1 cells could swim along the magnetic field in the blood plasma liquid, the fresh MO-1 cells swam at the beginning speed of about 300 μm/s, while the swimming speed decreased as the time pasted (Fig. 2a and2b). In the first 4 minutes, the MO-1 cells swam at the velocity of above 200 μm/s and then travelled at more than 100 μm/s in the following 20 minutes. In the next time, the MO-1 cells continued to swim at a lower speed below 100 μm/s. During the entire process, the MO-1 cells were capable of responding to the external magnetic field guidance at 37 °C (Fig. 2b). The experimental results showed that magnetotactic bacteria MO-1 could survive and keep their swimming motility in the plasma at the human physiological temperature for at least 45 minutes. When artificial urine was used as another model of body fluid, the state of magnetotactic bacteria MO-1 at the human physiological temperature (37 °C) was similar to that in the blood plasma (Fig. 3a). The MO-1 cells could be navigated by an external magnetic field and swim in the urine at a speed of above 200 μm/s during the first 6 minutes. The bacteria over time continued swimming at the speed of more than 130 μm/s for 20 minutes. Finally, their speed further decreased in the next time but still moved in the urine at a lower speed below 100 μm/s (Fig. 3a and3b). The swimming experimental results confirmed that the MO-1 cells could survive and swim for at least 40 minutes in the urine at physiological temperature. Targeted therapy or drug delivery in the human body mediated by magnetotactic bacteria probably need the ability of these bacteria to swim autonomously in the body fluid. The results of this study demonstrated that magnetotactic bacteria MO-1 are capable of surviving and swimming in the blood and urine. According to the average speed and surviving time, the distance that magnetotactic bacteria MO-1 swam in the blood for 45 minutes or in the urine for 40 minutes at physiological temperature could even reach up to approximately 40 cm. Considering the size of human organ and the principle of proximal intervention with external carriers, the magnetotactic bacteria MO-1 swimming distance is enough to satisfy the requirements of drug delivery and targeted therapy. The blood results are also similar to those when another magnetotactic bacteria of MC-1 were used [START_REF] Martel | Flagellated Magnetotactic Bacteria as Controlled MRI-trackable Propulsion and Steering Systems for Medical Nanorobots Operating in the Human Microvasculature[END_REF]. Other factors, such as blood flow velocity and vascular branch, were not considered in the present study. However, it's been reported that magnetotaxis enabled MTB to overcome counter-directional flow and migrate effectively through porous micromodels (Rismani Yazi, et al., 2017;Rismani Yazi, et al., 2018). As to the scene that the velocity of MO-1 cells actually decreased as the incubation time increased, the main reason is possibly the osmotic pressure induced by different salt concentrations. For example, the NaCl concentrations in the blood and urine are 8.5 and 4.2 g/L respectively, much lower than those (19.45 g/L) needed by MO-1 to grow. Consistently, it has also been implied that the velocity of the bacteria could also be controlled by varying the concentration of Na + in the medium [START_REF] Lefèvre | Isolation and characterization of a magnetotactic bacterial culture from the Mediterranean Sea[END_REF]. Besides, a magnetic field with a higher strength could slow down the speed of magnetotactic bacteria [START_REF] Yang | Dynamic Model and Motion Mechanism of Magnetotactic Bacteria with Two Lateral Flagellar Bundles[END_REF]. It's also suggested from the survival and movement abilities that magnetotactic bacteria do not proliferate spontaneously inside an animal nor cause severe infection and pathogenicity. Thus, their safety for the in vivo application is much better than that of other pathogenic bacteria.

Interaction of magnetotactic bacteria MO-1 with human HaCaT and Hela cells

Although the infection and pathogenicity of magnetotactic bacteria as a therapeutic tool may be negligible according to the analysis, the dose effect at the cellular level and the in vivo toxicity of magnetotactic bacteria MO-1 remained to be studied. In the present study, the biocompatibility and biological safety dose of magnetotactic bacteria MO-1 at the cellular level were examined via MTT assay. HaCaT and Hela cells (1 × 10 4 cells) were treated with magnetotactic bacteria MO-1 at different bacterial cell concentrations ranging from 0 to 5 × 10 7 bacteria/well in 100 μL DMEM in a 96-well plate for 24 h. As shown in Fig. 4a, the incubation of MO-1 cells at the concentration of less than 3.3 × 10 6 bacteria/well had no effect on the viability of the HaCaT cells. The cellular viability decreased by approximately 40% at the concentration of 6.5 × 10 6 bacteria/well and continued to decrease as the number of MO-1 cells was increased to 5 × 10 7 bacteria/well (Fig. 4a).

As for Hela cells, co-culture of MO-1 cells did not induce cell toxicity at the bacterial concentration of less than 1.6 × 10 6 bacteria/well (Fig. 4b). However, when the MO-1 concentration was at 3.3 × 10 6 bacteria/well, the viability of Hela cells decreased significantly to 70% and continued to decrease with the increased concentration of magnetotactic bacteria MO-1.

A potential reason for the increased cytotoxicity may be that increased magnetotactic bacteria touched with cells and then were phagocytosed, leading to increases in oxidative stress and reactive oxygen species inside the cell [START_REF] Wang | Magnetic nanoparticles trigger cell proliferation arrest of neuro-2a cells and ROS-mediated endoplasmic reticulum stress response[END_REF][START_REF] Wang | The interaction of bacterial magnetosomes and human liver cancer cells in vitro[END_REF].

Previous reports showed that 6.5 × 10 6 bacteria of magnetococcus marinus strain MC-1 incubated with 1 × 10 5 murine J774A.1 and NIH/3T3 cells in a 96-well plate induced a huge decrease in cellular viability to approximately 30% [START_REF] Taherkhani | Covalent Binding of Nanoliposomes to the Surface of Magnetotactic Bacteria for the Synthesis of Self-Propelled Therapeutic Agents[END_REF]. The decrease in cellular viability caused by magnetotactic bacteria MO-1 in the present study was much smaller than the results of strain MC-1. The difference may be attributed to the different cell lines and added bacterial species. The method of coating nanoliposome on the surface of MC-1 bacteria was proposed to reduce the toxicity of magnetotactic bacteria [START_REF] Taherkhani | Covalent Binding of Nanoliposomes to the Surface of Magnetotactic Bacteria for the Synthesis of Self-Propelled Therapeutic Agents[END_REF].

Determination of LD 50 of magnetotactic bacteria in rats

From an application perspective, the effect of magnetotactic bacteria on the whole human body must be considered. Different doses of magnetotactic bacteria MO-1 were injected into rats in each test group to acquire the LD 50 , and the rats were observed persistently for 1 week. As shown in Fig. 5, the survival rate decreased as the doses of magnetotactic bacteria MO-1 were increased. The experimental animals did not die in the groups of 5 × 10 8 bacteria/3 mL, one died in the 5 × 10 9 bacteria/3 mL group, three died in the 5 × 10 10 bacteria/3 mL group, and all animals died in the 5 × 10 11 bacteria/3 mL group. Behrens-Kärber method was used to calculate the LD 50 of the rats (0.2 kg) at 7.9 × 10 10 bacteria/kg, with a 95% CI of 1.99 × 10 10 -3.155 × 10 11 bacteria/kg. According to a guidance (U.S. Food and Drug Administration, 2005), the dose per unit body weight of an adult (60 kg) is equivalent to 0.16 of that of rats (0.2 kg). Therefore, the LD 50 of an adult human (60 kg) was calculated to be 1.264 × 10 10 MO-1 cells/kg. With the abilities of magnetic nanoparticle synthesis and autonomous swimming, magnetotactic bacteria have been used as contrast agents and in vivo drug delivery. As a contrast agent, magnetospirillum magneticum AMB-1 was intravenously injected into mice with 1 × 10 9 cells [START_REF] Benoit | Visualizing implanted tumors in mice with magnetic resonance imaging using magnetotactic bacteria[END_REF]. The number of magnetotactic bacteria MC-1 intravenously injected into rats for drug delivery was reported to be 1 × 10 8 cells [START_REF] Felfoul | Magneto-aerotactic bacteria deliver drug-containing nanoliposomes to tumour hypoxic regions[END_REF]. Thus, the actual quantity of magnetotactic bacteria for in vivo uses was much less than its LD 50 (7.9 × 10 10 MO-1 cells/kg in rats) obtained in the present study. Of course, different species of magnetotactic bacteria may have different LD 50 values, thereby needing further investigations. Regardless of the species of magnetotactic bacteria, after entering the body from the blood vessels or urethra, these bacteria may interact with the initially encountered organ or tissue and then interact with different kinds of cells, especially immune cells. Whether magnetotactic bacteria induced an intense immune response or not depends mainly on the immunogenicity and quantity of the bacteria. In any case, as a kind of Gram-negative bacteria, the endotoxins of magnetotactic bacteria MO-1 should be paid with sufficient attention when used in biomedical applications [START_REF] Alphandéry | Applications of magnetotactic bacteria and magnetosome for cancer treatment: A review emphasizing on practical and mechanistic aspects[END_REF]. On one hand, agents such as poly-L-lysine could be modified to reduce or eradicate the endotoxins of magnetotactic bacteria (Hamdous, et al., 2017;Mandawala, et al., 2017). On the other hand, rats could actually overcome some endotoxins of MO-1 cells according to the study.

Nevertheless, the endotoxins of magnetotactic bacteria / magnetosomes released from bacterial cells into the animal body should be considered seriously. Taken together, on the basis of the LD 50 obtained in the present study, the dose of MO-1 cells may be safe and it could satisfy the requirement of clinical tumor therapy.

CONCLUSIONS

In summary, the biocompatibility of magnetotactic bacteria was determined in this study. The MO-1 cells were able to swim along the magnetic field lines in the blood and urine and survive at 37 °C for approximately 40 minutes. When the HaCaT and Hela cells were incubated with magnetotactic bacteria MO-1, their cellular effect was dependent on the does and within the range of below 3.3 × 10 7 MO-1 cells/mL, which was well tolerated by the human cells. In the animal experiments, an LD 50 of 7.9 × 10 10 bacterial cells/kg in rats was determined. These data provided the key to understand the biocompatibility and safety of magnetotactic bacteria to a certain extent and suggested that magnetotactic ovoid strain MO-1 could basically be used in targeted therapy.

However, the distribution of whole bacteria in the body should be further investigated.

ACKNOWLEDGEMENT

We would like to thank Dr. Haitao Chen for the assistance in analyzing the swimming speed of magnetotactic bacteria MO-1 in the plasma and urine by using the MTrackJ plugin for ImageJ.

Fig. 1

 1 Fig. 1 Magnetotactic ovoid strain MO-1. a Magnetotactic ovoid strain MO-1. Black arrows indicated the two

Fig. 2

 2 Fig. 2 Speed of magnetotactic bacteria MO-1 in rat plasma under external magnetic field. The swimming

Fig. 3

 3 Fig. 3 Speed of magnetotactic bacteria MO-1 in artificial urine under external magnetic field. The

Fig. 4

 4 Fig. 4 Cell viability of HaCaT (a) and Hela (b) cells after interaction with magnetotactic bacteria MO-1.

Fig. 5

 5 Fig. 5 Determination of LD50 of magnetotactic bacteria in rats.

* Supported by the National Natural Science Foundation of China (No. 51937011, No. 51907192) and the National Key Research and Development Program of China (Grant NO. 2017YFC0108501).