1

Molecular archaeology

The genomic history of Ice Age Europeans

Ludovic Orlando

An extensive genomic time-series has been retrieved for 116 humans across Ice Age Europe. The data reveal how climate change drove range contraction and expansion of hunter-gatherer populations developing diverse cultures.

See p. xxx.

For most of the past hundred thousand years, the glacial conditions of the last Ice Age prevailed across Europe. About halfway through this period, the first anatomically modern humans arrived on a continent that had previously been inhabited only by Neanderthals 1,2 . Further groups of hunter-gatherers developed breathtaking art and diverse lithic industries and funerary practices across Europe, and then, around nine thousand years ago (kya), farmers from Anatolia started to settle 3 . Whether pre-farming cultures reflect isolated biological entities or regional differences amongst highly interconnected groups is highly debated. Also unknown is what happened as the climate changed, including during the Last Glacial Maximum (LGM), when large parts of the continent became uninhabitable 4 . On page xxx, Posth et al. 5 provide unprecedented resolution into the genomic makeup of Ice Age European hunter-gatherers, clarifying the biological contours of several cultures, and their relationships. They portray a highly dynamic history, punctuated by range contractions during the LGM, but also expansions that led to local admixture or population replacements as the climate warmed.

By reporting new genomic data for 116 ancient humans and reanalysing data from a further 240 known individuals across 34 countries -including the first western European Gravettian and Solutrean individuals -the authors add substantially to the known genetic atlas of ancient Europe between 35 and 5.2 kya (Fig. 1). To achieve this, they relied on established, cost-effective methodologies that can deliver human sequence data across 1.2 million genomic locations at which single DNA letters are known to show variability between individuals (single nucleotide polymorphisms, or SNPs) 6 -even for specimens that show minimal DNA preservation or are dominated by soil microbes. The authors also restricted their analyses to sequences with evidence of post-mortem DNA damage, to avoid contamination by modern individuals, which was especially pronounced in some specimens. After the LGM, and no later than 17 kya, non-local human groups entered from the Balkans into present-day Italy. They showed Near Eastern genetic affinities that define the 'Villabruna' cluster, and spread further south in small populations that reached Sicily at approximately 14 kya. The implication is that the post-LGM Epigravettian industry did not, as previously thought, develop in situ from earlier local Gravettian groups in Italy -suggesting that the peninsula may not have offered a glacial refugium to humans. The 'Villabruna' influence, however, extended outside Italy, as it was found at El Mirón in Spain around 19 kya. This reveals a southern fringe of Europe acting then as a corridor from the Balkans to Iberia.

Posth et al. also show that, after the LGM, populations expanded further outside their glacial refugia, forming a new genetic blend on the back of the 'Villabruna' and 'Fournol' ancestries. In present-day France, Belgium, Germany and Poland, they formed the 'Goyet-Q2' cluster, which is associated with the Magdalenian lithic industry. The new data also show that it was as late as 14 kya, or thereabouts, that one of the three main genetic ancestries found today in Europe (the WHG, for Western Hunter-Gatherers) was fully formed, as the 'Villabruna' and 'Goyet-Q2' mixed. This group -which Posth et al. now relabel 'Oberkassel' -appears to be genetically homogenous from present-day UK to Poland, indicating an expansion that was possibly related to the abrupt warming of the Bølling-Allerød period (14.8 to 12.9 kya).

  Posth and colleagues' work shows that the Gravettian lithic industry was developed by two main groups that populated Europe during the ten thousand years or so before the LGM. The first group -nicknamed 'Vestonice' in reference to the earliest associated site -stretched across present-day Italy, Czechia and Austria. It descended from earlier groups found further east, including western Russia. The second group, 'Fournol', inhabited western Europe (present-day France and Spain), where it descended from earlier local hunter-gatherers, such as those from the Goyet cave in Belgium (around 35 kya). By 27 kya, mixed ancestries were found at Goyet, indicating an east-to-west range expansion of the 'Vestonice' group right before the LGM. During the LGM (26 to 19 kya), human populations contracted within glacial refugia located at the southern fringe of Europe. In southern France and Iberia, this phenomenon was accompanied by the development of Solutrean flint tools, which Posth et al. reveal to have been produced by people deriving from local 'Fournol' groups. This model of population continuity between western Gravettian and Solutrean tool makers is further supported by the genetic profile of a 23-ky old individual from Andalusia, southern Spain, that another study now reveals 8 .

  Further east, another population (in which individuals of lighter skins and darker eye colours were more frequent at the time) eventually developed the Sidelkino culture of western Russia. Posth et al. detected 'Sidelkino'-related genetic ancestry in Spain, and it increased in the Baltics around the time of the thermal maximum, approximately 8 kya -a period when the 'Oberkassel' genetic influence also extended as far east as the Don-Volga region in Russia. Whether such population movements were the sole product of the then-warmer climate, or a response to farmer groups expanding outwards from Anatolia, remains unknown. Although the spread of Anatolian farmers into Europe kept changing the genetic makeup of local populations, Posth et al. find that individuals derived mainly from 'Oberkassel' ancestry persisted in Germany until the end of the Neolithic at roughly 5.2 kya, in line with a previous study of northern France 9 . This may reflect local communities limiting admixture with incoming farmers but adopting their lifestyle.

Fig. 1 |

 1 Fig. 1 | Map showing the locations of the ancient human individuals analysed. Posth et al. 5 have analysed the genomes of 356 ancient human individuals, including 102 newly reported,

Posth and colleagues have unfolded important population turnovers that have shaped the European genetic landscape in the face of climate change. Their findings caution against simple narratives equating lithic industries and ethnicity. More genetic analyses are, however, needed to gain finer-grained resolution into the formation and expansion timing of several groups (such as the 'Villabruna'), and their possible links to lithic industries presently overlooked (such as the 'Badegoulian'). In the future, sequencing of DNA molecules preserved in cave sediments may help to overcome the limitations of the skeletal genome assemblages available 10 , and even reveal the social structure of underlying human communities 11 , including those first entering Europe.

The European continent has a long tradition of archaeology, where excavation methods were developed from the early nineteenth century. It is therefore not surprising that this continent is also the best characterized archaeogenetically, with extensive genomic time-series now extending from the present day back to the Upper Palaeolithic. The resulting resources and increased understanding of the human past should be commended, but are also an invitation to redouble efforts outside Europe, to avoid developing a Eurocentric vision of human prehistory. For now, the fluidity of ancestries composing the deep genetic history of Europe provides an important lesson: no modern population can claim single origins from the human groups first established on the continent.

Ludovic Orlando is the director of the Centre for Anthropobiology and Genomics of Toulouse (CNRS UMR 5288), Université Paul Sabatier, 37 Allées Jules Guesde, 31000 Toulouse, France. e-mail: ludovic.orlando@univ-tlse3.fr


Student Member, IEEE Ahmed Amghayrir

Noël Tanguy
email: noel.tanguy@univ-brest.fr

Pascale Bréhonnet

Pierre Vilbé

Léon-Claude Calvez

Laguerre-Gram reduced-order modeling

Keywords: Model-order reduction, Laguerre functions, infinite dimensional systems

We present an efficient model reduction procedure based on the Laguerre description of the system to be approximated. Using a one-order operator defined in the Laplace domain we construct a pencil of functions and formulate the problem as the minimization of the L 2 ̟ R + criterion. The use of a weight function in the inner product definition allows a control of the time-error spreading in model reduction procedure. We show how the required Gram matrix can be computed efficiently and prove that the impulse response of the reduced model is also in L 2 ̟ R + . The transfer function approach allows an immediate and promising application in model reduction of infinite dimensional systems. An extension to MIMO systems is also given.

I. INTRODUCTION

Laguerre functions have shown their large potential in numerous applications as in signal analysis and parameter identification [START_REF] Clement | Laguerre functions in signal analysis and parameter identification[END_REF], system identification [START_REF] Wahlberg | System identification using Laguerre models[END_REF], [START_REF] Chou | Continuous-time identification of SISO systems using Laguerre functions[END_REF], approximation of finite or infinitedimensional system [START_REF] Mäkilä | Approximation of stable systems by Laguerre filters[END_REF], [START_REF]Laguerre series approximation of infinite dimensional systems[END_REF], industrial control [START_REF] Dumont | Laguerre-based adaptive control of pH in an industrial bleach plant extraction stage[END_REF],... More recently, starting with a state-space representation of a system, a Laguerre description has been used to compute a reduced-order model of the system. The method is based on Padé approximation of the Laguerre spectra, associated to a Krylov subspace decomposition technique [START_REF] Knockaert | Laguerre-svd reduced-order modeling[END_REF]- [START_REF] Chen | Model reduction in the time-domain using Laguerre polynomials and Krylov methods[END_REF]. In this paper we present a model reduction technique using a transfer-function formalism that allows reduction of irrational transfer functions. The derivative and integral operators usually employed to construct a pencil of functions for the reduction procedure are here replaced by a one-order operator more suited to construct a set of basis functions from the Laguerre description of the system. The model reduction problem is then formulated as a L 2 ̟ R + criterion minimization in which the required Gram matrix can be efficiently computed. The use of a weight function in the inner product definition permits a control of the time-error spreading of the reduced model. Assuming that the original function belongs to L 2 ̟ R + the presented procedure yields a reduced model whose impulse response provably belongs to L 2 ̟ R + . An extension to MIMO systems is also derived.

The paper is organized as follows. From the Laguerre representation of the system we construct in section II a set of linearly independent functions. On the base of this set of functions we formulate in section III the model reduction problem as the minimization of a L 2 ̟ R + criterion and solve it. The belonging to L 2 ̟ R + of the resultant funtion is proved and the extension of the reduction procedure to MIMO systems is given. A numerical example of the model reduction method is then presented in section IV.

II. PENCIL OF FUNCTIONS CONSTRUCTION

The real causal two-parameters Laguerre functions φn (t) related to the Laguerre polynomials Ln (x) by φn (t)

√ γe -αt Ln (γt) , with α > 0, γ > 0 and

Ln (x) e x n! d n dx n
x n e -x , have the following Laplace transforms

φn (s) = √ γ s + α s -γ + α s + α n , n = 0, 1, 2, ... (1)
Defining the weighted inner product of two real causal functions by

f, g ∞ 0 ̟ (t) f (t) g (t) dt, (2)
where the weight function is given by

̟ (t) e -(γ-2α)t , let L 2 ̟ R + denote the Hilbert space of measurable R-valued functions f (t) for which f, f f 2 ̟ < ∞.
It will be noted that a necessary (but not sufficient) condition for f (t) to belong to

L 2 ̟ R + is γ 2 -α > ℜ (-si)
, where (-si) stands for the poles of f (s) [START_REF] Malti | Some results on the convergence of transfer function expansion on Laguerre series[END_REF]. The Laguerre functions constitute an orthonormal basis

in L 2 ̟ R + satisfying φn, φm = δn,m,
with δn,m the Krönecker symbol. In mostly engineering literature the particular choice γ -2α = 0 is made and corresponds to an unit weight function in the inner product definition [START_REF] Wahlberg | System identification using Laguerre models[END_REF]. A different choice of the (γ -2α) value modifies the temporal part of the functions in the inner product and will yield to modify the error-spreading of the impulse response of the reduced model along the temporal axis. Since the Laguerre functions constitute an orthonormal basis in L 2 ̟ R + , the Laplace transform f (s) of any function f ∈ L 2 ̟ R + can be represented by a Laguerre series

f (s) = ∞ n=0 fn φn (s) , (3)
where the Laguerre spectrum {fn} n≥0 is given by fn = f, φn . Using the orthogonality property, the weighted inner product defined in (2) can be replaced by a discrete sum usually more convenient for numerical computation

f, g = ∞ n=0 fngn. (4)
Denoting F (z) the z-transform of the Laguerre spectrum {fn} n≥0 , F (z) can be derived using the bilinear transformation z = s+α s+α-γ which maps the open left plane in s-domain delimited by ℜ (s) < γ 2α into the open unit disk in z-domain

f (s) → F (z) ∞ n=0 fnz -n = √ γz z -1 f γz z -1 -α . ( 5
)
It is worth noting that the inner product (4) may be written in the z-domain as

f, g = 1 2πj |z|=1 F (z) G 1 z dz z . (6)
An important property of the bilinear transformation is the preservation of the order of complexity of transfer functions. As for example if f (s), an irreducible analytic transfer function in some s-domain, is given by f

(s) = p m-1 s m-1 +...+p 1 s+p 0 qms m +...+q 1 s+q 0 therefore F (z) = z(u m-1 z m-1 +...+u 1 z+u 0) vmz m +...+v 1 z+v 0
will be an irreducible, analytic and rational transfer function in corresponding z-domain. The reciprocal transformation giving f (s) from F (z) is

F (z) → f (s) = √ γ s + α F s + α s -γ + α . (7)
Now let us define

f (s) → Λ f (s) s + α s -γ + α f (s) - γ s + α f (γ -α) , (8)
and construct recursively Ωr = f0, f1, ..., fr using

f0 (s) f (s) , fi+1 (s) Λ fi (s)
for i = 0, 1, ..., r -1. Denoting {fi,n} n≥0 the Laguerre spectra of the functions fi (s) and Ξr = {F0, F1, ..., Fr} the set of functions Fi (z) derived from fi (s) via transformation [START_REF]Laguerre series approximation of infinite dimensional systems[END_REF], one obtains the relation linking the successive functions Fi (z) as

Fi+1 (z) = z [Fi (z) -fi,0] (9)
with F0 (z) = F (z). Therefore it follows that the Laguerre spectrum of fi+1 (s) is readily obtained from the Laguerre spectrum of fi (s), and so on

fi+1,n = fi,n+1 = ... = f0,n+i+1 = fn+i+1 (10)
for all n ≥ 0, i ≥ 0. Hence no computation is required and so no error is made when computing the Laguerre spectrum of any member of the set Ωr. Moreover using ( 8) or ( 9) and fi,0 = √ γ fi (γα) = lim z→∞ Fi (z) deduced from (5) one should readily verify that Λ operator preserves the natural frequencies or poles in the Laplace domain. Therefore Ωr constitutes an efficient set of approximating functions for determining a r-order reduced denominator.

III. MODEL REDUCTION

To construct a r-order reduced model of f (s), consider now the approximation of fr (s) using the set Ωr-1. Denoting ê (s) the local error and ai the negatives of the coefficients of the approximation, we can write r i=0 ai fi (s) + ê (s) = 0, where ar = 1.

(11)

Applying transformation [START_REF]Laguerre series approximation of infinite dimensional systems[END_REF] to equation ( 11) yields

r i=0 aiFi (z) + E (z) = 0, (12)
where

E (z) = √ γz z -1 ê γz z -1 -α .
Therefore, in view of (6), we infer that minimizing e, e in Laplace domain to construct an r-order reduced model of f (s) is equivalent to minimizing e, e in z-domain to construct an r-order reduced model of F (z).

Denoting Ψ the Gram matrix constituted of the inner products ψi,j fi, fj for i, j = 0, 1, ..., r -1 and b [ψ0,r ψ1,r ... ψr-1,r] T where T denotes the transpose, then the best r-vector a [a0 a1 ... ar-1] T in the sense of minimizing the error energy e, e in [START_REF] Calvez | Contribution à l'étude des propriétés de la transformation Z et de la transformation de Laguerre[END_REF] can be obtained by solving the linear system

Ψ a = -b. ( 13
)
From ( 4) and [START_REF] Malti | Some results on the convergence of transfer function expansion on Laguerre series[END_REF] the inner products can be expressed exclusively from the Laguerre spectrum of f (s) as

ψi,j = ∞ n=0 fn+ifn+j, (14)
which constitute absolutely convergent series assuming that f (t) belongs to L 2 ̟ R + . Inner products ψi,j are thus replaced by sums, in practice finite, that are simple to calculate. Moreover, cost computation of the Gram matrix will be drastically reduced using the symmetry property ψi,j = ψj,i and the recurrence relation deduced from [START_REF] Nagaoka | Mullis-Roberts-type approximation for continuous-time linear systems[END_REF] ψi,j = ψi-1,j-1fi-1fj-1, i, j = 1, 2, ...

Repeating the use of ( 9) and using [START_REF] Malti | Some results on the convergence of transfer function expansion on Laguerre series[END_REF] yields

Fi (z) = z i F0 (z) - i-1 j=0 fjz i-j . ( 16
)
Substituting ( 16) in ( 12) and solving for F0 (z) yields

F0 (z) = zBr-1 (z) -E (z) Ar (z) (17)
where Br-1 (z) Because E (z) is an error term, assuming that e, e is sufficiently small, equation [START_REF] Morvan | Simplified algorithm for Laguerre approximation of URC networks[END_REF] suggests that

F (z) = zBr-1 (z) Ar (z) (18)
is a r-order candidate model for F (z) = F0 (z). Therefore applying reciprocal transformation [START_REF] Knockaert | Laguerre-svd reduced-order modeling[END_REF] to (18) yields

f (s) = p (s) q (s) = √ γ r i=1 ai i-1 j=0 fj (s + α) i-j-1 (s -γ + α) r-i+j r i=0 ai (s + α) i (s -γ + α) r-i
(19) with ar = 1, a r-model candidate for approximating f (s).

Starting with the Laguerre spectrum {fn} n≥0 of f (s), the model reduction procedure to construct a r-order model of f (s) can be summarized as follows:

1. Form the Gram matrix Ψ and b [ψ0,r ψ1,r ... ψr-1,r] T , by computing the required inner products in the following way:

-compute ψ0,j = ∞ n=0 fnfn+j for j = 0, 1, ...r, -recursively compute ψi,j = ψi-1,j-1fi-1fj-1, for i = 1, 2, ...r -1 and j = i, i + 1, ..., r, -complete the lower triangular part of Ψ using the symmetry property ψi,j = ψj,i.

2. Solve Ψ a =b for a = [a0 a1 ... ar-1] T .

Form the reduced model using (19).

Remarks:

-The Laguerre expansion for some simple transfer functions can be found in [START_REF] Calvez | Contribution à l'étude des propriétés de la transformation Z et de la transformation de Laguerre[END_REF], but usually deducing the infinite Laguerre expansion (3) and computing the Gram matrix Ψ of an irrational transfer function is not a straightforward problem. In practice a N -order truncated Laguerre series (N ≫ r) should be used in the place of (3); assuming that the truncated error is sufficiently small the model reduction will give a good approximation of the original transfer function. It will be noted that the Laguerre spectrum {fn} 0≤n<N of f (s), can always be evaluated by the Weeks method [START_REF] Tricomi | Transformazione di Laplace e polinami di Laguerre[END_REF], [START_REF] Weeks | Numerical inversion of Laplace transforms using Laguerre functions[END_REF].

-It will be noted that the reduced model given by ( 19) preserves the r-first Laguerre coefficients fj of f (s). As the Laguerre spectrum of f (s) is related to the derivatives of f (s) at s = γα [START_REF] Calvez | Contribution à l'étude des propriétés de la transformation Z et de la transformation de Laguerre[END_REF] by

d n f (s) ds n s=γ-α = 1 √ γ n! (-γ) n n j=0 n j (-1) j fj
where n j stands for the binomial coefficients, this implies that the reduced model given by ( 19) preserves the first derivatives of f (s) at s = γα, i.e.

d n f (s) ds n s=γ-α = d n f (s) ds n s=γ-α
, n = 0, 1, ..., r -1.

-Moreover, taking q (s) = r i=0 ai (s + α) i (sγ + α) r-i as the denominator of the reduced model, the choice of the numerator can be optimized in the sense of minimizing the quadratic error f -f , f -f . The problem is then a well-known linear problem which leads to the resolution of f (ν * i + γ -2α) = f (ν * i + γ -2α) for i = 1, 2, ..., r where ν * i stands for the negatives of the complex conjugate of the zeros of q (s).

A. Belonging to L 2

̟ R + consideration Assuming that f (s) is a transfer function of degree r0 ≥ r, Ωr-1 is a set of linearly independent functions and then the Gram matrix Ψ constituted of the inner products ψi,j = fi, fj (i, j = 0, 1, ..., r-1) is a positive definite matrix. Consider now the Lyapunov equation Ψ-

A T ΨA = C where A 0 T I -a represents the companion
form of the reduced system in the state-space representation. Taking into account that a satisfies (13) and using [START_REF] Johnson | A high-q distributed-lumped-active network configuration with zero real-part pole sensitivity[END_REF], the solution of the Lyapunov equation is given by C = c c T + M where c [f0 f1 ... fr-1] T and M 0 0 T 0 ε with ε = e, e . As the quadratic error ε is always non-negative, C matrix is semi-positive definite, therefore Ar (z) cannot have any zero outside the unit circle and thus q (s) cannot have any zero in the right half-plane delimited by ℜ (s) > γ 2α. Furthermore, provided that {A, C} is observable, the reduced model belongs to L 2 ̟ R + [START_REF] Nagaoka | Mullis-Roberts-type approximation for continuous-time linear systems[END_REF]. In the particular case γ = 2α, which corresponds to a unit weighting function in (2), the model-reduction procedure then yields a model proved to be asymptotically stable.

B. Extension to MIMO systems

In MIMO case, the system under consideration is described by a matrix of transfer functions F (s) = f λ,µ (s) . The problem is then to derive a matrix F (s) of transfer functions possessing the same r-order denominator q (s) i.e.

F (s) = f λ,µ (s) = [p λ,µ (s)] q (s) , (20)
in the sense of minimizing the quadratic error ε = λ,µ ε λ,µ where ε λ,µ = e λ,µ , e λ,µ . The solution is given by solving

  λ,µ Ψ λ,µ   a = -   λ,µ b λ,µ   ,
where the Gram matrices Ψ λ,µ and the vectors b λ,µ are computed for each transfer function f λ,µ (s) with the same Laguerre parameters α and γ. Fig. 1. Impulse responses of the original system, of its 4th-order model for (α 1 ; γ 1 ) = (18.5; 37.0) (dashed line), for (α 2 ; γ 2 ) = (8.95; 47.9) (dot line) and of its 4th-order Padé model (dashed-dot line)

Starting with the Laguerre spectrum of each function fλ,µ (s) computed for common values α and γ of the Laguerre parameters, the procedure of MIMO system reduction then follows:

1. Compute the Gram matrices Ψ λ,µ and the vectors b λ,µ relative to the transfer functions f λ,µ (s).

2. Form the Gram matrix Ψ = λ,µ Ψ λ,µ and the vector b = λ,µ b λ,µ then solve Ψ a =b for a = [a0 a1 ... ar-1] T . 3. Construct the reduced model F (s) given in (20) by using ( 19) for all the values of λ and µ.

Note that the numerators p λ,µ (s) could be optimized by solving

f λ,µ (ν * i + γ -2α) = fλ,µ (ν * i + γ -2α)
for i = 1, 2, ..., r where ν * i stands for the negatives of the complex conjugate of the zeros of q (s).

IV. NUMERICAL EXAMPLE

To illustrate the utility of the method we present a typical application of model reduction of an infinite-dimensional system. The system considered is a distributed RC circuit whose irrational transfer function is [START_REF] Johnson | A high-q distributed-lumped-active network configuration with zero real-part pole sensitivity[END_REF], [START_REF] Morvan | Modélisation de circuits et systèmes de dimension infinie[END_REF]

ĝ (s) = 1 1 + sinh( √ RCs) √ RCs with RC = 1.
From ĝ (s), the first N = 150 coefficients of a Laguerre model are recursively computed using [START_REF] Morvan | Simplified algorithm for Laguerre approximation of URC networks[END_REF] for two different weight functions ̟1 (t) = 1 and ̟2 (t) = e -30t . The corresponding Laguerre parameters (α1; γ1) = (18.5; 37.0) and (α2; γ2) = (8.95; 47.9) used for the computation of the truncated Laguerre spectrum of ĝ (s) came from the optimization method described in [START_REF] Morvan | Modélisation de circuits et systèmes de dimension infinie[END_REF], [START_REF] Tanguy | Optimum choice of free parameter in orthonormal approximations[END_REF]. Using the technique presented in this paper the following 4th-order reduced models are derived as g 1 (s) = -0.397s The impulse response of the original transfer function and those of the 4th-order models and Padé approximation are presented Fig. 1 and show the great accuracy of the presented model-reduction procedure. The time zoom near t = 0 underlines the role played by the weight function ̟ (t) = e -(γ-2α)t : a choice of γ > 2α, as used to derive g 2 (s), produces an improvement of the reduced model quality around t = 0 (to a comparative detriment of the quality for t greater). If we define the relative weighted quadratic error by Q g -g, g -g / g, g we obtain Q1 = 2.18 • 10 -4 and Q2 = 6.10 • 10 -4 for the Laguerre-Gram reduced models, and Q P adé = 1.32 • 10 -3 (computed for the weighting function ̟ (t) = 1) for the Padé model. These results confirm the great accuracy of the model reduction procedure.

V. CONCLUSION

An efficient procedure for model reduction of finite or infinite dimensional transfer function described by a Laguerre expansion has been presented. Based on a one-order operator in the Laplace domain is constructed a set of basis functions from which is applied a procedure of quadratic error minimization. The use of a weight function in the inner product definition permits a control of the timeerror spreading of the reduced model. Integrals defining the required inner products are replaced by sums (finite in practice) expressed using the Laguerre spectrum of the original transfer function and can be computed recursively, allowing a numerically convenient calculation. We have shown that the reduced model preserves the first derivatives of the transfer function at a chosen point s = γα and has an impulse response that belongs to L 2 ̟ R + . The procedure has also been extended to MIMO systems. Presented work takes part inter alia of model reduction of infinite dimensional systems. An illustrative example has shown the efficiency of the procedure.

  j-1 and Ar (z) r i=0 aiz i , with ar = 1, are polynomials in z respectively of order r -1 and r.

ACKNOWLEDGMENT

The authors wish to acknowledge the Brittany Region for financial support. They also would like to thank the reviewers for their valuable comments and suggestions.


N Tanguy

L C Calvez

Gram matrix of a Laguerre model: application to model reduction of irrational transfer function

An efficient method for evaluating inner products of repeated integrals and derivatives of a function described by a Laguerre model is presented. An immediate and promising application is model reduction of infinite dimensional transfer functions.

Introduction

The Gram matrix Ψ of a system is very useful in system identification, power density spectrum modeling and model reduction. The Gram matrix contains elements that are inner products of repeated integrals and/or derivatives of the impulse response of the system. In model reduction applications, the elements of this matrix Ψ have to be computed from the mathematical model of the original system. Note that only stable systems will be considered.

Let Ω = {f 1 ,f 2 ,...,f q ,...,f r+1 } denote a set of contiguous real functions where

f q (t) f (t) , f i+1 (t)
df i (t) dt (i = q, q +1, ..., r) ,

f i-1 (t) t ∞ f i (τ ) dτ (i = q, q -1, ..., 2) , (1)
f (t) being the impulse response of the system with transfer function F (s). Clearly, the computation of the Gram matrix involves evaluation of inner products of real functions

ψ i,j f i ,f j = ∞ 0 f i (t) f j (t) dt. (2)
Derivative and integral operators defined in (1) preserve natural frequencies or poles in the * Email address: Noel.Tanguy@univ-brest.fr Laplace domain, therefore Ω constitutes an efficient set of approximation functions for determining a r-order reduced denominator. Moreover, during model-reduction procedure, integral operator will slightly emphasize the lower frequencies, when derivative operator will slightly emphasize the higher frequencies.

Usually f (t) will be the impulse response of the system, but assuming that f (t) is a wellbehaved square-integrable function with Laplace transform F (s), other responses of the system could be considered (for example the step response shifted of its final value).

When systems are described by rational transfer functions, integrals defined by (2) can be efficiently computed in the frequency domain [START_REF] Calvez | Evaluation of scalar products of repeated integrals of a function with rational Laplace transform[END_REF][START_REF] Calvez | Efficient evaluation of model-reduction related integrals via polynomial arithmetic[END_REF][START_REF] Lucas | Evaluation of scalar products of repeated integrals by Routh algorithm[END_REF]. Here we consider systems described by Laguerre models. Indeed the Laguerre functions have shown their large potential in numerous applications as in signal analysis and parameter identification [START_REF] Clement | Laguerre functions in signal analysis and parameter identification[END_REF], system identification [START_REF] Chou | Continuous-time identification of SISO systems using Laguerre functions[END_REF][START_REF] Wahlberg | System identification using Laguerre models[END_REF], approximation of finite or infinite-dimensional system [START_REF] Mäkilä | Approximation of stable systems by Laguerre filters[END_REF][START_REF] Mäkilä | Laguerre series approximation of infinite dimensional systems[END_REF], numerical inversion of the Laplace transform [START_REF] Abate | On the Laguerre method for numerically inverting Laplace transforms[END_REF][START_REF] Tricomi | Transformazione di Laplace e p o l i n a m id iL a g u e r r e ,R .C .A c c a d .N a t[END_REF][START_REF] Ward | The calculation of transients in dynamical systems[END_REF][START_REF] Weeks | Numerical inversion of Laplace transforms using Laguerre functions[END_REF], industrial control [START_REF] Dumont | Laguerre-based adaptative control of pH in an industrial bleach plant extraction stage[END_REF]... This is more particularly in the context of the numerical inversion of infinite-dimensional Laplace transfer functions that we have based our works. That restrictive choice does not withdraw anything from the generality of the presented results.

Let F (s) denote the transfer function of the system described by a Laguerre series

F (s)= ∞ n=0 c n Φ n (s)( 3 )
where {c n } n≥0 is the Laguerre spectrum of the system, Φ n (s) are the Laplace transforms of the Laguerre functions defined by

Φ n (s) √ 2α s + α s -α s + α n (4)
for n =0, 1, 2,... and (-α) is a real pole satisfying α>0. In practice the Laguerre series (3) will be truncated and the Laguerre coefficients c n will be evaluated by a well-known technique based on discrete Fourrier transform computation [START_REF] Abate | On the Laguerre method for numerically inverting Laplace transforms[END_REF][START_REF] Tricomi | Transformazione di Laplace e p o l i n a m id iL a g u e r r e ,R .C .A c c a d .N a t[END_REF][START_REF] Ward | The calculation of transients in dynamical systems[END_REF][START_REF] Weeks | Numerical inversion of Laplace transforms using Laguerre functions[END_REF]. Now, let F i (s) denote the Laplace transform of the function f i (t)a n d{c i,n } n≥0 its Laguerre spectrum. As the Laguerre functions form an orthonormal basis, the inner products defined by ( 2) can be written in Laguerre spectral domain as

ψ i,j = ∞ n=0 c i,n c j,n . (5)
Integrals defined by ( 2) are thus replaced by sums, in practice finite, that are simpler to calculate and numerically more robust. Before evaluating these sums, Laguerre spectra of the repeated derivatives and integrals of f (t) must be computed.

Derivatives and integrals computation

A well-known relation binding the Laplace transform and the z-transform of the Laguerre spectrum {c n } n≥0 is deduced by substituting (4) in (3), then carrying out the transformation s →

α z+1 z-1 ,o n eo b t a i n s Z {c n } n≥0 ∞ n=0 c n z -n = √ 2α z z -1 F α z +1 z -1 (6)
where Z denotes the z-transform operator. This useful relation will enable us to determine the Laguerre spectra of repeated integrals and derivatives of the function f (t).

Let us first consider the derivative operator defined by [START_REF] Calvez | Evaluation of scalar products of repeated integrals of a function with rational Laplace transform[END_REF], in Laplace domain it yields

F i+1 (s)=sF i (s) -f i (+0) . (7)
Substitution of ( 7) in ( 6) leads to

Z {c i+1,n } n≥0 = α z +1 z -1 Z {c i,n } n≥0 - √ 2α z z -1 f i (+0)
which yields the recurrence relation

c i+1,0 = αc i,0 - √ 2αf i (+0) c i+1,n = α (c i,n + c i,n-1 )+c i+1,n-1 (8)
for n =1 , 2, 3 ... and i = q, q +1,...,r.W h e n the initial value f i (+0) is not known it can be computed as follows

f i (+0) = lim s→∞ sF i (s)= √ 2α ∞ n=0 c i,n . (9)
In practice, the Laguerre series is usually truncated at order N . Therefore, and more particularly when the initial value f i (+0) is not readily available, an approximation of the coefficients of the various Laguerre spectra can be computed as follows

c i+1,N -1 = -α c i,N -1 c i+1,n-1 = c i+1,n -α ( c i,n + c i,n-1 ) (10)
for n = N -1,N -2,...,1, i = q, q +1,...,r and where the starting point of the computation is c q,n = c q,n . The approximated Laguerre series of F i (s)i st h e ng i v e nb y

F i (s)= N -1 n=0 c i,n Φ n (s) .
An approximated value of f i (+0) can then be evaluated using [START_REF] Mäkilä | Laguerre series approximation of infinite dimensional systems[END_REF] or ( 9) i.e.

f i (+0) = √ 2α N -1 n=0 c i,n = 1 √ 2α (α c i,0 -c i+1,0 ) . (11)
Note that the procedure [START_REF] Tricomi | Transformazione di Laplace e p o l i n a m id iL a g u e r r e ,R .C .A c c a d .N a t[END_REF] gives a set of functions verifying the strict equality

f i+1 (t)= d f i (t) dt
Integral and derivative operators defined in (1) are reciprocal, therefore the recurrence relations binding the Laguerre spectra of the repeated integrals of the function f (t) are readily obtained from ( 8)

c i-1,0 = 1 α c i,0 + √ 2αf i-1 (+0) c i-1,n = 1 α (c i,n -c i,n-1 ) -c i-1,n-1 (12)
for n =1 , 2, 3,... and i = q, q -1,...,2. Dealing with truncated Laguerre models, and more particularly when the initial value f i-1 (+0) is not readily available, an approximation of the Laguerre coefficients can be recursively computed in the following way

c i-1,N -1 = -1 α c i,N -1 c i-1,n-1 = 1 α ( c i,n -c i,n-1 ) -c i-1,n (13)
for n = N -1,N -2,...,1, i = q, q -1,...,2 and where the starting point of the computation is c q,n = c q,n . Note that the procedure (13) gives a set of functions verifying the strict equality

f i-1 (t)= t ∞ f i (τ ) dτ.

Gram matrix computation

The Gram matrix of repeated integral and derivative functions possesses interesting properties [START_REF] Sévellec | Réduction de l'ordre des systèmes à functions de transfert rationnel[END_REF][START_REF] Vilbé | Suboptimal model reduction via least-square approximation of time-response by its derivatives and integrals[END_REF], that allow reduction of the computing cost of the inner products [START_REF] Chou | Continuous-time identification of SISO systems using Laguerre functions[END_REF]. Indeed integration by parts of (2) yields

ψ i,j =[f i-1 (t) f j (t)] ∞ 0 - ∞ 0 f i-1 (t) f j+1 (t) dt.
Because f (t) is square integrable by assumption, this relation simplifies to

ψ i,j = -f i-1 (+0) f j (+0) -ψ i-1,j+1 (14)
and, in the particular case j = i -1o n eo b t a i n s

ψ i,i-1 = - 1 2 f 2 i-1 (+0) . ( 15
)
The calculation algorithm of the Gram matrix is then the following:

1. Compute the Laguerre spectra of successive integrals and derivatives or their approximations using ( 12) or ( 13) and ( 8) or [START_REF] Tricomi | Transformazione di Laplace e p o l i n a m id iL a g u e r r e ,R .C .A c c a d .N a t[END_REF].

2. Compute the elements of the main diagonal of the Gram matrix using [START_REF] Chou | Continuous-time identification of SISO systems using Laguerre functions[END_REF] or in practice the following truncated sums

ψ i,i ≈ N -1 n=0 c 2 i,n or ψ i,i ≈ N -1 n=0 c 2 i,n .
3. Compute the elements of the first secondary diagonal using (15).

4. Recursively compute the other elements of the lower triangular part of the Gram matrix using [START_REF] Sévellec | Réduction de l'ordre des systèmes à functions de transfert rationnel[END_REF].

5. The upper triangular part of the matrix is readily obtained using the symmetry property ψ i,j = ψ j,i .

Note that if initial values are not readily available, it is then possible to evaluate them using [START_REF] Ward | The calculation of transients in dynamical systems[END_REF]. Remark: It will be noted that the described procedures ( 8) and ( 12) or [START_REF] Tricomi | Transformazione di Laplace e p o l i n a m id iL a g u e r r e ,R .C .A c c a d .N a t[END_REF] and [START_REF] Dumont | Laguerre-based adaptative control of pH in an industrial bleach plant extraction stage[END_REF] give the Laguerre spectra of repeated strict derivatives and integrals. Therefore, in model reduction application, the reduced model of F (s) is provably stable by a Lyapunov method if the Gram matrix is computed on Ω or on Ω= f 1 , f 2 ,..., f q ,..., f r+1 and if the model reduction procedure is based on the approximation of either first or last function of the set by the others [START_REF] Nagaoka | Mullis-Roberts-type approximation for continuous-time linear systems[END_REF].

Example

To illustrate the utility of the method we present a typical application of model reduction of an infinite-dimensional system. The system considered is the celebrated underwater cable whose irrational transfer function is [START_REF] Hsu | Linear approximation of fractional transfer functions of distributed parameter systems[END_REF][START_REF] Morvan | Modélisation de circuits et systèmes de dimension infinie[END_REF] G (s)=e - √ Ks (with K =1). From G (s), the first N = 100 coefficients of a Laguerre model with α =2 .42 is computed using a well-known technique [START_REF] Abate | On the Laguerre method for numerically inverting Laplace transforms[END_REF][START_REF] Tricomi | Transformazione di Laplace e p o l i n a m id iL a g u e r r e ,R .C .A c c a d .N a t[END_REF][START_REF] Ward | The calculation of transients in dynamical systems[END_REF][START_REF] Weeks | Numerical inversion of Laplace transforms using Laguerre functions[END_REF] based on discrete Fourrier transform computation. For ease G (s)= -0.0261s 5 +7.41s 4 -1048s 3 + 63397s 2 + 92844s + 6326 s 6 +98.3s 5 + 4568s 4 + 67984s 3 + 227350s 2 + 132925s + 6686 of experimentation, however, we have coded our own version of the algorithm using a discrete cosine transform. The (not crucial) choice α =2.42 for the computation of the Laguerre spectrum of G (s) is obtained using the optimization method described in [START_REF] Tanguy | Optimum choice of free parameter in orthonormal approximations[END_REF]. In order to derive reduced-order models taking the form of rational functions of order r = 6, Gram matrices of size (r +1)× (r +1) are calculated. We used the method described in [START_REF] Sévellec | Réduction de l'ordre des systèmes à functions de transfert rationnel[END_REF][START_REF] Vilbé | Suboptimal model reduction via least-square approximation of time-response by its derivatives and integrals[END_REF] to derived some reduced models corresponding to the values q =1 , 2, ..., r +1. W e retained the reduced model that minimizes the relative quadratic error defined by

ε ∞ 0 [g (t) -g (t)] 2 dt ∞ 0 [g (t)] 2 dt
where g (t) is the impulse response of the 6thorder model and g (t) is the exact impulse response of G (s)g i v e nb y

g (t)= 1 2K √ π K t 3 2 e -K 4t .
The selected reduced model G (s), given at the top of the page, is based on the Gram matrix constructed with the scalar products of the Laguerre models of G (s), its 4 first derivatives and its 2 first integrals (i.e. q = 3). The corresponding quadratic error ε =5 .24 × 10 -4 confirms the great accuracy of the reduced model.

Conclusion

An efficient procedure for computing the Gram matrix of repeated derivatives and integrals of a function described by a Laguerre model has been presented. Integrals defining the required inner products are replaced by sums (finite in practice) allowing a calculation numerically more convenient. The deduced recurrence relations yield

Figure 1 .

 1 Figure1. Impulse responses of the original system and of its 6th-order model (dashed line)


Trong Bach Nguyen

Laurine Mule Mueni

Tran Nu Thanh

Viet Bui

Huynh Nguyen

Duy Bao

Nguyen Thi

Kim Cuc

Taco Nicolai

Trong Bach Nguyen

Nguyen Thi Kim

Characterization of tuna dark muscle protein isolate

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

viscosity decreased with decreasing of pH, which correlated with the decreasing size of the aggregates.

Keywords: tuna, dark muscle, fish protein isolate, extraction, aggregation

Practical application

Tuna dark muscle (TDM) represents a significant fraction of industrially processed tuna. It is currently considered a low value by-product, because it has poor sensory qualities. However, TDM has a high protein content that does not suffer from sensory defects if it is extracted using the pHshift method. TDM protein isolate therefore has the potential to be a valuable nutritious food ingredient. In order to develop industrial food products with this fish protein isolate it is necessary to characterize it thoroughly.

INTRODUCTION

Proteins can be effectively isolated from fish by-products by dissolving the proteins at high or low pH and subsequently precipitating them at intermediate pH where the net charge of the proteins is close to zero (Tang et al., 2020). Fish protein isolate (FPI) obtained in this manner can be used as nutritious ingredients in industrial food products (Tang et al., 2020;Desai et al., 2019Desai et al., & 2021;;Asad Nawaz et al., 2020). Naturally, in order to optimize the usage of FPI for different types of applications, it is necessary to thoroughly characterize its properties in aqueous solution at different pH.

In industrial processing of tuna fish 50-70% of its body mass is considered as by-products (head, skin, dark meat, viscera, bone, etc) of which about 12% is dark muscle [START_REF] Guérard | Production of tuna waste hydrolysates by a commercial neutral protease preparation[END_REF][START_REF] Herpandi | The tuna fishing industry: A new outlook on fish protein hydrolysates[END_REF]. It has low commercial value, but is rich in nutrients and its extracted protein displays important technological properties [START_REF] Sánchez-Zapata | Quality Characteristics of Dark Muscle from Yellowfin Tuna Thunnus albacares to Its Potential Application in the Food Industry[END_REF]. Tuna dark muscle (TDM) contains approximately 27% protein with high nutritional value, but it also contains hemoglobin and myoglobin that can oxidize and thereby render the muscle dark [START_REF] Chaijan | Physicochemical and gelling properties of shortbodied mackerel (Rastrelliger brachysoma) protein isolate prepared using alkaline-aided process[END_REF][START_REF] Ochoa | Evaluation of the efficacy, safety and acceptability of a fish protein isolate in the nutrition of children under 36 months of age[END_REF][START_REF] Tahergorabi | Fish protein isolate: Development of functional foods with nutraceutical ingredients[END_REF][START_REF] Tian | Nutritional and Digestive Properties of Protein Isolates Extracted from the Muscle of the Common Carp Using pH-Shift Processing[END_REF]. In recent years, protein recovery from fish dark muscle has attracted the interest of both scientists and manufacturers as it contributes to improved economic efficiency and towards zero-waste processing. In which the recovery method and characterization of protein are the main concerns.

The amount of protein that is dissolved from fish muscle depends strongly on the pH [START_REF] Batista | Recovery of proteins from fish waste products by alkaline extraction[END_REF][START_REF] Geirsdottir | Solubility and viscosity of herring (Clupea harengus) proteins as affected by freezing and frozen storage[END_REF]. It was found that at pH > 10 or pH < 3.5 more than 80% of the proteins is soluble, whereas at pH close to the pI  5.2-5.5 it is less than 10% [START_REF] Kristinsson | Fish protein isolate by pH shift[END_REF]. Therefore an efficient method to extract fish protein is to solubilize fish meat in a strong base or acid and subsequently to precipitate the proteins close to pI [START_REF] Davenport | Channel catfish (Ictalurus punctatus) muscle protein isolate performance processed under different acid and alkali pH values[END_REF][START_REF] Foh | Influence of pH shift on functional properties of protein isolated of tilapia (Oreochromis niloticus) muscles and of soy protein isolate[END_REF][START_REF] Kim | New approaches for the effective recovery of fish proteins and their physicochemical characteristics[END_REF][START_REF] Lee | Preparation and characterization of protein isolate from Yellowfin tuna Thunnus albacares roe by isoelectric solubilization/precipitation process[END_REF][START_REF] Surasani | Acid and alkaline solubilization (pH shift) process: a better approach for the utilization of fish processing waste and by-products[END_REF][START_REF] Surasani | Recovery and characterization of proteins from pangas (Pangasius pangasius) processing waste obtained through pH shift processing[END_REF][START_REF] Surasani | Recovery of Proteins from Rohu Processing Waste Using pH Shift Method: Characterization of Isolates[END_REF][START_REF] Undeland | Recovery of functional proteins from herring (Clupea harengus) light muscle by an acid or alkaline solubilization process[END_REF][START_REF] Zhong | Fish Protein Isolates Recovered from Silver Carp (Hypophthalmichthys molitrix) By-Products Using Alkaline pH Solubilization and Precipitation[END_REF])(Zhen-Xing Tang et al., 2020). This method removes the dark color and the lipids [START_REF] Chaijan | Physicochemical and gelling properties of shortbodied mackerel (Rastrelliger brachysoma) protein isolate prepared using alkaline-aided process[END_REF][START_REF] Sánchez-Zapata | Quality Characteristics of Dark Muscle from Yellowfin Tuna Thunnus albacares to Its Potential Application in the Food Industry[END_REF] resulting in a fish protein isolate (FPI) that is useful for food applications. Fish muscle can be dissolved both in strong acid and base, but alkali treatment has advantages in terms of the yield and the strength of gels formed by aqueous FPI solutions [START_REF] Abdollahi | Structural, functional, and sensorial properties of protein isolate produced from salmon, cod, and herring by-products[END_REF][START_REF] Batista | Recovery of proteins from fish waste products by alkaline extraction[END_REF][START_REF] Chaijan | Physicochemical and gelling properties of shortbodied mackerel (Rastrelliger brachysoma) protein isolate prepared using alkaline-aided process[END_REF][START_REF] Davenport | Channel catfish (Ictalurus punctatus) muscle protein isolate performance processed under different acid and alkali pH values[END_REF][START_REF] Kaewdang | Characteristics of Gelatin Extracted from the Swim Bladder of Yellowfin Tuna (Thunnus albacores) as Affected by Alkaline Pretreatments[END_REF][START_REF] Kim | New approaches for the effective recovery of fish proteins and their physicochemical characteristics[END_REF][START_REF] Oliyaei | Effect of Temperature and Alkaline pH on the Physicochemical Properties of the Protein Isolates Extracted from the Whole Ungutted Lanternfish (Benthosema pterotum)[END_REF][START_REF] Zhong | Fish Protein Isolates Recovered from Silver Carp (Hypophthalmichthys molitrix) By-Products Using Alkaline pH Solubilization and Precipitation[END_REF].

The molecular weight of the protein components of FPI and its distribution have been determined by electrophoresis (Brenner et al., 2009a;[START_REF] Foh | Influence of pH shift on functional properties of protein isolated of tilapia (Oreochromis niloticus) muscles and of soy protein isolate[END_REF][START_REF] Geirsdottir | Solubility and viscosity of herring (Clupea harengus) proteins as affected by freezing and frozen storage[END_REF][START_REF] Kristinsson | Recovery and properties of muscle proteins extracted from tilapia (Oreochromis niloticus) light muscle by pH shift processing[END_REF][START_REF] Sun | Effect of pH shifting on conformation and gelation properties of myosin from skeletal muscle of blue round scads (Decapterus maruadsi)[END_REF][START_REF] Tian | Nutritional and Digestive Properties of Protein Isolates Extracted from the Muscle of the Common Carp Using pH-Shift Processing[END_REF][START_REF] Yongsawatdigul | Effects of alkali and acid solubilization on gelation characteristics of rockfish muscle proteins[END_REF].

The results show prominent bands corresponding to myosin heavy (200 kDa) and light (17-23 kD) chains, actin ( 43 kDa) and tropomyosin ( 38 kDa). The secondary structure of tilapia protein isolate (TPI) was characterized by [START_REF] Kobayashi | FT-IR and Raman spectroscopies determine structural changes of tilapia fish protein isolate and surimi under different comminution conditions[END_REF] using Fourier Transform Infrared spectroscopy (FT-IR). The presence of -helical and -sheet structures was demonstrated by absorption bands near 1648 cm -1 and 1660 cm -1 . [START_REF] Brenner | Characterization of fish myosin aggregates using static and dynamic light scattering[END_REF] characterized FPI extracted from cod in dilute aqueous solutions at pH 11 using light scattering techniques. They found that the protein was present in the form of aggregates with a z-average radius of gyration R g = 150 nm and a hydrodynamic radius R h = 180 nm. The FPI also contained about 15 wt% much larger aggregates. The dependence of the scattering intensity on the scattering wave vector showed that the aggregates had a self-similar structure characterized by a fractal dimension of 1.8.

The viscosity of FPI in aqueous solution was investigated by [START_REF] Foh | Influence of pH shift on functional properties of protein isolated of tilapia (Oreochromis niloticus) muscles and of soy protein isolate[END_REF] and [START_REF] Geirsdottir | Solubility and viscosity of herring (Clupea harengus) proteins as affected by freezing and frozen storage[END_REF]. [START_REF] Geirsdottir | Solubility and viscosity of herring (Clupea harengus) proteins as affected by freezing and frozen storage[END_REF] observed that the viscosity of herring protein isolate suspensions increased at pH values far from the pI. [START_REF] Foh | Influence of pH shift on functional properties of protein isolated of tilapia (Oreochromis niloticus) muscles and of soy protein isolate[END_REF] found that at pH 7 the viscosity of tilapia protein isolate increased with increasing protein concentration. Brenner et al. (2009a) reported that reduction of the pH from 11 to between 9.5 and 8.5 led to formation of homogeneous gels, whereas reduction to lower pH values down to 3.8 led to formation of heterogeneous gels showing syneresis.

It is clear from this brief review that FPI extracted from white flesh and mixed fish byproducts has already been studied quite extensively. However, FPI from dark tuna muscle has been rarely investigated. [START_REF] Sánchez-Zapata | Quality Characteristics of Dark Muscle from Yellowfin Tuna Thunnus albacares to Its Potential Application in the Food Industry[END_REF] reported on the iron and myoglobin content of tuna dark muscle as well as the water and oil holding capacity. Here we present an investigation of the amino acid composition, structure, solubility, and physico-chemical properties of FPI extracted from tuna dark muscle at pH 12 as a function of the pH. The objective is to establish if, and to what extent, FPI from tuna dark muscle is different from other types of FPI so as to facilitate its use in food products.

MATERIALS AND METHODS

Materials

The fish protein isolate used for this study was extracted from yellowfin tuna (Thunnus Albacares) dark muscle provided by Hai Vuong Co., LTD (Khanh Hoa province, Vietnam) using the alkaline extraction method described by [START_REF] Undeland | Recovery of functional proteins from herring (Clupea harengus) light muscle by an acid or alkaline solubilization process[END_REF] with slight modifications. Milled dark muscle was washed twice with deionized water (3:1 ratio) at 4 °C to remove the dark color.

Excess water was removed by wrapping the muscle in a filter bag with a pore size of 50 µm and manually pressed. The cleaned muscle was solubilized by adding water at pH 12 at a ratio 1:5. The mixture was homogenized at 4 °C for 2 min at a speed of 3500 rpm using an Ultra-turrax (model: IKA T18, Germany). 2 N NaOH was used to adjust to pH 12 after homogenization. The mixture was stirred for 3 hours at room temperature (20 °C) using an Orbital shaker (model: NB-101M, Korea) and subsequently centrifuged (MF600, Labkorea) at 5000 rpm for 30 min at 4 °C. The supernatant was decanted through a filter bag with a pore-size of 5 µm and the soluble protein was precipitated by reducing the pH to 5.5 using 2 N HCl. The precipitate was thoroughly washed with distilled water at 4 °C. Excess water was removed by manually wrapping and squeezing the protein isolate in a filter bag with a pore size of 5 µm. Finally, the FPI was freeze dried (Lyobeta 35,Spain) and stored at -20 °C. The chemical composition of FPI and TDM is shown in Table 1.

Sample preparation

Stock solutions of FPI were prepared by dissolving the freeze-dried FPI at a concentration of 30 g/L in deionized water at pH 12 by stirring for a few hours at ambient temperature. Solutions at lower concentrations and lower pH were prepared by dilution with deionized water and addition of small amounts of 0.1 N HCl.

Methods

Proximate composition

The composition of tuna dark muscle (TDM) and protein isolate (FPI) extracted from TDM was determined using the procedures of the AOAC (2012). The crude protein content, moisture content, crude fat and ash were determined using the Kjeldahl method, method 928.08, method 948.04 and method 920.153, respectively. Sodium, NaCl, and histamine contents were determined using methods 969.23, 930.23, 997.13 (LOD = 3 ppm), respectively. Heavy metal (As, Hg, Pb, and Cd) contents were determined using methods 986.15, 971.21, 972.23, and 973.34, respectively.

Fourier Transform Infrared (FT-IR) Spectroscopy

FPI powder was characterized using a FT-IR spectrometer (Model Alpha, Bruker, Ettlingen, Germany) equipped with a horizontal attenuated total reflectance (HATR) accessory. The internal reflection crystal was made of zinc selenide, and had a 45° angle of incidence to the IR beam. The sample was mixed with KBr in the ratio of 98:2 (w/w). A total of 32 scans were done in the wavenumber range from 500 to 4000 cm -1 with a resolution of 4 cm -1 . The spectra were collected at room temperature (20 °C) using approximately 0.5 g of the sample that was placed on the surface of the HATR crystal. Analysis of the spectra was carried out using the OPUS 7.0 data collection software program.

Amino acid composition

The amino acid composition of the samples was determined using the method ISO 13903 (2005). In brief, FPI samples were hydrolyzed for 24 hours in 6 N HCl at 110 °C in a sealed vessel.

The samples were then dissolved in deionized water and filtered. About 20 μL hydrolyzed sample or the standard was injected into the amino acid analyzer (Biochrom 30+, Biochrom Ltd., Cambridge, UK) equipped with a high pressure PEEK column (No: H-0563) packed with Ultropac 8 cation exchange resin (No: 132-56). The automated analyzer was configured with post-column ninhydrin derivatization. The ninhydrin derivative of Pro was detected at 440 nm whereas the remaining amino acid derivatives were detected at 570 nm.

Zeta potential

The zeta potential (ζ) of aqueous FPI solutions at 0.1 g/L was determined as a function of the pH by measuring the electrophoretic mobility (U E ) using a zeta potential analyser (Nano Partica SZ-100Z, Horiba Scientific, Japan). The zeta potential was calculated from the electrophoretic mobility by applying Henry's equation [START_REF] Lunardi | Experimental methods in chemical engineering: Zeta potential[END_REF]:

U E = 2ε. ζ. 𝑓(κα) 3η (1)
where, ε is the permittivity, f(α) is the Henry coefficient and η is the dispersion viscosity.

Following the Smoluchowski approximation f(α) was set at 1.5. All measurements were done in triplicate.

Light scattering

Dynamic and static light scattering measurements were conducted with a commercial apparatus (ALV-CGS3, ALV-Langen) operating with a vertically polarized He-Ne laser with wavelength  = 632 nm. The temperature was controlled at 20 °C with a thermostat bath to within  0.2 °C. The measurements were performed on FPI solutions diluted to below 1 g/L in order to avoid the effect of interactions and at angles of observation between 13° and 140°. The relative scattering intensity (I r ), of the solute, calculated as the intensity minus the solvent scattering divided by the scattering intensity of toluene, is related to the weight average molar mass (M W ) and the scattering wave vector (q) dependent structure factor (S(q)) [START_REF] Brown | Dynamic light scattering: the method and some applications[END_REF]:

I r KC = M W S(q) (2)
where K is an optical constant defined as:

K = 4π 2 n S 2  4 N A ( 𝜕n 𝜕C ) 2 ( n tol n S ) 2 1 R tol (3)
with N A is Avogadro's number, ∂n/∂C = 0.189 mL/g the refractive index increment and R tol = 1.35 × 10 -5 1/cm the Rayleigh constant of toluene at 20 °C.

The normalized electric field autocorrelation functions (g 1 (t)) obtained from dynamic light scattering measurements [START_REF] Berne | Laser Spectroscopy: Dynamic Light Scattering[END_REF] were analyzed in terms of a distribution of relaxation times (τ):

   ) ( ) / exp( ) ( ) ( 1    d t A t g (4)
In dilute solution, the relaxation was caused by translational diffusion of the particles. The diffusion coefficient (D) was calculated from the average relaxation rate:

2 1 / q D     (5)
The z-average hydrodynamic radius (R h ) of the solute was calculated from D extrapolated to q = 0 using the Stokes-Einstein equation:

D kT R h  6  (6)
with η the viscosity of the solvent, k Boltzman's constant, and T the absolute temperature.

Confocal Laser Scanning Microscopy

The morphology of FPI solutions was studied with Confocal Laser Scanning Microscopy (CLSM) in the fluorescence with mode utilizing a Zeiss LSM 800 microscope (Oberkochen,

Germany). A water immersion objective lens was used (HCxPL APO 63x NA=1.2 with a theoretical resolution of 0.1 µm in the x-y plane. The proteins were fluorescently labeled by adding 5 ppm rhodamine B that physically bound to the proteins. No effect of labeling was observed at the concentration used. The fluorescent label was excited by a laser at 543 nm and the fluorescent light was detected between 543 and 900 nm. The solutions were inserted between a concave slide and a coverslip and hermetically sealed.

Viscosity measurement

The viscosity was measured at 20 °C as a function of the shear rate using a rheometer (Kinexus Pro 50N, Malvern, England) with a cone and plate geometry (diameter 60 mm; angle 1°).

After loading, the samples were pre-sheared at 100 s -1 for 1 min. The viscosity was determined as a function of the shear rate between 0.1-1000 s -1 .

Sodium dodecyl sulphate-polyacrylamide gel electrophoresis (SDS-PAGE)

FPI solutions at 5 g/L and different pH were mixed with the same volume of an electrophoresis sample buffer (0.25 M Tris (pH 6.8), 2% SDS, 10% glycerol, 100 mM dithiothreitol, 0.2% bromophenol blue) and heated at 95 o C for 10 min. The running gel contained 15% polyacrylamide in 0.5 M Tris-HCl (pH 8.8) and 0.3% SDS, whereas the stacking gel contained 5% polyacrylamide in 0.5 M Tris HCl (pH 6.8) and 0.2% SDS. The electrode buffer contained 0.025 mol/L Tris-HCl, 0.192 M glycine and 0.1% SDS at pH 8.3. SDS-PAGE electrophoresis was performed at 15 mA using a Mini Protean II unit (Bio-Rad Laboratories, Inc., Richmond, CA, USA). Following electrophoresis, the gels were stained with 0.05% Coomassie Blue G250 in 25% v/v ethanol and 10% v/v acetic acid for 4 h. The gel was finally destained with an aqueous solution containing 5% (v/v) methanol and 10% (v/v) acetic acid. Molecular weights of the FPI subunits were estimated using molecular weight standards (Invitrogen, USA) in the range from 3.5 kDa to 260 kDa.

Solubility

Protein solutions at 10 g/L and different pH were centrifuged at 1.3 x 10 4 g for 30 min at 20 °C by a refrigerated centrifuge (Hettich, Mikro 22R, Germany). The protein concentration of the supernatant was determined by UV absorption at 278 nm with a UV-Visible spectrometer (Libra S50 Bio, Cambridge, England) using an extinction coefficient 1.0 L/g.cm. The soluble protein fraction was determined by dividing the protein concentration in the supernatant by the total concentration. All measurements were done in triplicate.

Statistical analysis

Statistical analysis of the data was performed using SPSS software (Version 16.0, Chicago, IL, USA). One-way analysis of variance (ANOVA) by the Student's T-test was used for statistical comparisons of means. Standard deviations and differences at P < 0.05 are indicated.

RESULTS AND DISCUSSION

FPI powder characterization

Amino acid composition

The amino acid composition of tuna dark muscle and fish protein isolate extracted from TDM is given in Table 2. As expected, the total protein content in FPI was higher (90%) than in TDM (73%), but the amino acid composition of the proteins was similar with glutamic acid, aspartic acid, leucine, and lysine as the most abundant amino acids. The fraction of essential amino acids (EAA) was 44% in TMA and 45% in FPI. The amino acid composition was similar to that reported for muscle tissue of yellowfin tuna by [START_REF] Peng | Amino acid and fatty acid composition of the muscle tissue of yellowfin tuna (Thunnus albacares) and bigeye tuna (Thunnus obesus)[END_REF] and for FPI extracted from pangas (Pangasius pangasius) processing waste [START_REF] Surasani | Recovery and characterization of proteins from pangas (Pangasius pangasius) processing waste obtained through pH shift processing[END_REF].

Fourier-Transform Infrared Spectroscopy

As far as we are aware, the only previous FT-IR study of FPI was reported by [START_REF] Kobayashi | FT-IR and Raman spectroscopies determine structural changes of tilapia fish protein isolate and surimi under different comminution conditions[END_REF] for tilapia protein isolate. These authors focused on the amide I (1600 -1700 cm -1 ) and amide II (1500-1600 cm -1 ) regions. They assigned the amide I band with a maximum at 1639 cm -1 to the -sheet structure and the amide II with a maximum at 1550 cm -1 to the -helix structure.

Figure 1 compares the FT-IR spectra of TDM and the FPI powder. Small, but significant differences were observed indicating that the extraction procedure had some influence on the secondary structure of the proteins. The height of the amide II band relative to the amide I band was lower for FPI and a shoulder appeared around 1740 cm -1 , indicating that FPI contained relatively less -helices. The amide II peak for FPI from tuna was higher than for FPI paste from talapia reported by [START_REF] Kobayashi | FT-IR and Raman spectroscopies determine structural changes of tilapia fish protein isolate and surimi under different comminution conditions[END_REF]. The higher α-helix content was perhaps caused by the high concentrations of amino acids (alanine, glutamic acid, lysine and leucine), which generally prefer to adopt an α-helix structure. However, one needs to be careful when comparing these results directly, because the FPI paste investigated by [START_REF] Kobayashi | FT-IR and Raman spectroscopies determine structural changes of tilapia fish protein isolate and surimi under different comminution conditions[END_REF] contained 75% water, whereas the FPI powder studied here contained very little water.

Characterization of tuna dark muscle protein isolate solutions at different pH

Zeta-potential

The zeta potential () of FPI as a function of the pH is shown in figure 2. It changed from strongly positive at pH < 4 to strongly negative at pH ≥ 6. The isoelectric point (pI) was close to pH 5.

Solubility

The solubility of FPI in water was determined as a function of the pH between 2 and 12.

Figure 3 shows that FPI solutions were transparent between pH 6 and pH 12, but became increasingly turbid at lower pH down to pH 5.0 indicating the formation more and larger aggregates. At pH 5.0 a dense precipitate was formed. Protein aggregates also precipitated at pH 4, but the precipitate was less dense. With further decrease of the pH the turbidity decreased indicating that fewer and smaller aggregates were formed.

The fraction of soluble proteins (F pis ) was determined by measuring the protein concentration of the supernatant after centrifugation as described in the section 2.3.9, see figure 4.

At pH ≥ 6.0 and pH ≤ 3.0, F pis was higher than 0.85, but close to the pI, F pis decreased to values as low as 0.12. These results are similar to those reported by [START_REF] Foh | Influence of pH shift on functional properties of protein isolated of tilapia (Oreochromis niloticus) muscles and of soy protein isolate[END_REF] for tilapia protein isolate and [START_REF] Sun | Effect of pH shifting on conformation and gelation properties of myosin from skeletal muscle of blue round scads (Decapterus maruadsi)[END_REF] for myosin extracted from skeletal muscle of blue round scads. Similar observations were also reported for white muscle from cod, rockfish and Alaska pollock (Brenner et al., 2009a;[START_REF] Kim | Negative roles of salt in gelation properties of fish protein isolate[END_REF][START_REF] Yongsawatdigul | Effects of alkali and acid solubilization on gelation characteristics of rockfish muscle proteins[END_REF]. The low solubility close to pI can be attributed to weak electrostatic repulsions between the proteins.

SDS analysis

The protein composition of soluble FPI was analyzed by SDS-PAGE at different pH above pH 6 and below pH 3. The results show numerous bands spread out over the whole range of molecular weights from 15 kDa to 260 kDa independent of the pH, see figure 5. The distinct band at 15 kDa can be attributed to myosin light chains. Indistinct bands at 37 kDa and 42 kDa correspond to tropomyosin and actin, respectively [START_REF] Abdollahi | Structural, functional, and sensorial properties of protein isolate produced from salmon, cod, and herring by-products[END_REF][START_REF] Taktak | Physicochemical, textural, rheological and microstructural properties of protein isolate gels produced from European eel (Anguilla anguilla) by heatinduced gelation process[END_REF].

However, albumins (from the sarcoplasmic fraction) are also typically found in this region [START_REF] Felix | Physicochemical, microstructure and bioactive characterization of gels made from crayfish protein[END_REF]. Remarkably no distinct band corresponding to myosin heavy chains (MHC) was observed at 250 kDa. This could indicate that polymerization of MHC had occurred resulting in larger proteins that could not enter the gel (Y. S. [START_REF] Kim | Negative roles of salt in gelation properties of fish protein isolate[END_REF][START_REF] Nunez-Flores | Importance of salt and temperature in myosin polymerization during surimi gelation[END_REF].

Alternatively, the absence was caused by a high degree of hydrolysis due to the well-known high proteolytic enzyme activity in the dark muscle [START_REF] Abdollahi | Structural, functional, and sensorial properties of protein isolate produced from salmon, cod, and herring by-products[END_REF][START_REF] Kristinsson | Recovery and properties of muscle proteins extracted from tilapia (Oreochromis niloticus) light muscle by pH shift processing[END_REF][START_REF] Zhong | Fish Protein Isolates Recovered from Silver Carp (Hypophthalmichthys molitrix) By-Products Using Alkaline pH Solubilization and Precipitation[END_REF], which would explain the presence of proteins with an almost continuous range of lower molecular weights.

Size & molar mass

The size and molecular weight of FPI in aqueous solution was determined using static and dynamic light scattering at 20 °C. The z-average hydrodynamic radius (R h ) and weight average molar mass (M w ) were determined as described in section 2.3.5 and are shown as a function of the pH in figure 6. Aggregates with R h  230 nm and M w  10 7 g/mol were formed at the highest pH (11 and 12). The aggregate size and molar mass decreased with decreasing pH down to R h  90 nm and M W  1.4 x 10 6 g/mol at pH 7-9. At pH 3 the aggregates were again slightly larger than at pH 6, but at pH 2 the size was close to that at pH 7-9. Most FPI precipitated between pH 3 and 6 and could not be analyzed by light scattering. It is curious that the aggregate size increased strongly between pH 9 and pH 11 even though the zeta potential remained approximately constant. As was mentioned in the Introduction, [START_REF] Brenner | Characterization of fish myosin aggregates using static and dynamic light scattering[END_REF] observed R h = 180 nm for FPI extracted from cod at pH 11.

Microstructure

The microstructure of FPI at 10 g/L was investigated with CLSM at different pH, see figure 7. Between pH 12 and pH 5.8, where the solutions were transparent, the CLSM images were homogeneous showing only a few microscopic particles. Significantly more microscopic particles were seen at pH 5.7 and at lower pH down to pH 3.5 larger protein flocs were formed (Figure 7). At pH 3 and lower the images were again homogeneous (not shown). The larger flocs sedimented during centrifugation and at some pH values even under gravity as was shown in section 3.2.2.

Viscosity

The viscosity of FPI solutions at 10 g/L at different pH was measured as a function of the shear rate, see figure 8. Significant shear thinning was observed at all pH values. At a given shear rate the viscosity was highest at the highest pH values (11 and 12) and decreased strongly with decreasing pH down to pH 6. At pH 3 the viscosity was again higher, but at pH 2 it was close to that in the pH range 6-9. Comparison with the results on the aggregate size shown in section 3.2.4 shows the pH dependence of the viscosity was correlated with that of the average aggregate size -FPI solutions containing larger aggregates being more viscous. The effect of the concentration of FPI solutions is shown in Figure 9 at pH 2 and pH 12. As expected, the viscosity increased with increasing in protein concentration. The increase was more important for the solution at pH 12 that contained larger aggregates. These results can be explained by the fact that larger aggregates take up a larger volume fraction as has been discussed in more detail for whey protein aggregates by [START_REF] Inthavong | Viscosity of mixtures of protein aggregates with different sizes and morphologies[END_REF].

CONCLUSIONS

FPI can be extracted from tuna dark muscle using the pH-shift method. The FPI was rich in essential amino acids and kept the secondary structure of raw material. More than 85% of the FPI proteins could be solubilized at pH ≥ 6 and pH ≤ 3 in the form of aggregates that were significantly larger at pH 10 and 12 than between pH 6 and 9 and below pH 3. At intermediate pH FPI formed large flocs that sedimented under gravity. Soluble FPI solutions were shear thinning with a viscosity that was increased strongly for pH ≥ 10 where the solutions contained larger aggregates. A significant fraction of the myosin heavy chains was broken down into smaller chains. The FPI extracted in this manner has an excellent potential for application in food products for human consumption.

PROXIMATE

Graphical Abstract

FIGURE 2 FIGURE 3

 23 FIGURE 2 Zeta potential of FPI as a function of pH.

FIGURE 4 FIGURE 5

 45 FIGURE 4Fraction of proteins in the supernatant as a function of pH for FPI solutions at 10 g/L.

FIGURE 6

 6 FIGURE 6 Average hydrodynamic radius (circle) and molar mass (triangle) of aggregates of

FIGURE 7 FIGURE 8 FIGURE 9 FIGURE 1 FIGURE 2 FIGURE 3

 789123 FIGURE 7 CLSM images of FPI solutions at 10 g/L at different pH. The images represent 110

FIGURE 4 FIGURE 5

 45 FIGURE 4Fraction of proteins in the supernatant as a function of pH for FPI solutions at 10 g/L.

FIGURE 7 FIGURE 8 FIGURE 9

 789 FIGURE 7 CLSM images of FPI solutions at 10 g/L at different pH. The images represent 110

TABLE 2

 2 Amino acid composition of tuna dark muscle and protein isolate (dry basis).

	Amino Acid	TDM (g/100g)	FPI (g/100g)
	Alanine (Ala)	5.04 ± 0.30 a	5.11 ± 0.62 a
	Arginine (Arg)	4.21 ± 0.14 a	5.63 ± 0.57 b
	Aspartic acid (Asp)	7.21 ± 0.14 a	9.24 ± 0.3 b
	Glutamic acid (Glu)	11.14 ± 0.68 a	16.21 ± 0.63 b
	Glycine (Gly)	3.70 ± 0.20 a	3.38 ± 0.12 a
	Histidine (His) **	3.57 ± 0.27 a	2.40 ± 0.13 b
	Isoleucine (Ile) *	3.95 ± 0.30 a	4.89 ± 0.23 b
	Leucine (Leu) *	6.49 ± 0.24 a	8.29 ± 0.35 b
	Lysine (Lys) *	7.77 ± 0.33 a	9.42 ± 0.33 b
	Methionine (Met) *	2.07 ± 0.07 a	3.00 ± 0.24 b
	Phenylalanine (Phe) *	3.14 ± 0.16 a	4.16 ± 0.26 b
	Proline (Pro)	2.36 ± 0.26 a	2.38 ± 0.25 a
	Serine (Ser)	2.98 ± 0.21 a	3.93 ± 0.17 b
	Threonine (Thr) *	3.07 ± 0.12 a	4.21 ± 0.15 b
	Tyrosine (Tyr *	2.60 ± 0.36 a	3.26 ± 0.11 b
	Valine (Val) *	3.48 ± 0.35 a	5.04 ± 0.22 a
	EAA	32.56	42.26

ACKNOWLEDGEMENTS

This research is funded by Vietnam National Foundation for Science and Technology Development (NAFOSTED) under grant number 106. 99-2018.42.

DATA AVAILABILITY STATEMENT

The data that support the findings of this study are available from the corresponding author upon request.

Authors also authenticate the co-collaborations of the Center of Analytical Services and Experimentation of Ho Chi Minh City (HCMC, Vietnam) for amino acid analysis. And we are also grateful to Prof. Christophe Chassenieux, Prof. Lazhar Benyahia and Dr. Frédéric Renou at IMMM (Le Mans, France) for the valuable discussions.

Many thanks to Nguyen Thi Anh Ky and Tran Thi Hong Van for their assistance in sample preparations.

CONFLICTS OF INTEREST

The authors declared no conflicts of interest.