Thi Thuy

Van Nguyen
email: ttv.nguyen@uliege.be

Cédric Heuchenne
email: c.heuchenne@ulg.ac.be

Kim Phuc Tran

Thi Thuy Van Nguyen

Kim Phuc Tran
email: kim-phuc.tran@ensait.fr

Anomaly Detection for Compositional Data using VSI MEWMA control chart

Keywords: Compositional data, Markov chain, VSI-MEWMA, control chart, Data Science

come   L'archive ouverte pluridisciplinaire

INTRODUCTION

In the global competitive economy nowadays, an extremely important task for manufacturing companies is to not only offer high-quality products but also reduce waste and increase efficiency in the production processes. The development of advanced technologies in Artificial Intelligence and Data Science fields makes this task more possible, but also more challenging when competing with other companies. Therefore, making a smart decision in manufacturing becomes a crucial task in any production company. In this context, early detection of abnormal products as well as assignable causes to fix the production system as soon as possible is an indispensable part, and Statistical Control Process (SPC) is one of the most effective methods to accomplish this task. Through control charts, SPC helps manufacturing companies monitor product qualities and discover the defects in the production lines. In SPC literature, many studies have been done to design a variety of control charts for monitoring different types of process data, see Montgomery 2013. Among these different data, compositional data (CoDa) are vectors whose components are strictly positive and they often present the proportions, percentages, or frequencies of some whole. Their applications can be found in many domains such as chemical research, econometrics, and the food industry, see [START_REF] Aitchison | The Statistical Analysis of Compositional Data[END_REF]. Due to the constraint on the sum of components of the CoDa vector, it can not be treated as normal data.

In SPC literature, the studies in the control charts for monitoring CoDa data are still limited. In Boyles 1997, a chi-squared C control chart in case p > 3 was proposed in Vives- [START_REF] Vives-Mestres | Signal Interpretation in Hotelling's T 2 Control Chart for Compositional Data[END_REF][START_REF] Tran | Monitoring compositional data using multivariate exponentially weighted moving average scheme[END_REF], the authors proposed a MEWMA-CoDa chart for monitoring CoDa with arbitrary components. This type of control chart was shown to be effective in detecting small to moderate process shift sizes and outperforming its competitor (T 2 -CoDa chart). The influence of measurement errors on the performance of T 2 , MEWMA chart for monitoring CoDa were investigated in [START_REF] Zaidi | Performance of the hotelling t2 control chart for compositional data in the presence of measurement errors[END_REF]Zaidi et al. 2020, respectively. In these control charts, the authors suggested using an isometric log-ration (ilr) transformation to transform CoDa to vector in R p-1 space to handle the constraint of CoDa and the average run length (ARL) to evaluate the performance of proposed control charts.

In the control charts mentioned above, the fixed sampling interval (FSI) was supposed to use. Recently, the design of control charts tends to use variable sampling intervals (VSI). In these charts, the sampling interval between two consecutive samples is allowed to vary due to the value of the current control statistic. Many studies on the VSI control chart have been published so far, see, for example, [START_REF] Castagliola | Monitoring the coefficient of variation using a variable sampling interval control chart[END_REF][START_REF] Nguyen | Monitoring the ratio of two normal variables using variable sampling interval ewma control charts[END_REF], among many others. As we know, the VSI MEWMA control chart for monitoring CoDa has not been used. Consequently, in this study, we propose a VSI MEWMA type control chart for monitoring CoDa, namely VSI MEWMA-CoDa, with arbitrary components based on ilr transformation. The modification of the Markov chain approach proposed by Lee 2009 will be used to compute average time to signal (ATS), criteria to access the performance of VSI control charts.

The rest of this paper is organized as follows: In Section 2, the modeling of CoDa and the suggested isometric logratio transformation are introduced; the VSI MEWMA-CoDa control chart together with the Markov chain approach and optimization procedure to find the optimal parameters are given in Section 3; in Section 4, the performance of the VSI MEWMA-CoDa chart with different scenarios are provided; conclusions and some recommendations for further researches are given in Section 5.

MODELING OF COMPOSITIONAL DATA

By definition, a row vector, x = (x 1 , x 2 , . . . , x p ), is a p-part composition when its components are strictly positive and they carry only relative information, see [START_REF] Aitchison | The Statistical Analysis of Compositional Data[END_REF][START_REF] Pawlowsky-Glahn | Modeling and Analysis of Compositional Data[END_REF]. The relative information here refers only to the proportions between components of the composition, regardless of their numerical values. The sum of the components of x, p i=1

x i , is a constant κ. For instance, κ = 100 refers to measurements in percentage while κ = 1 means that the measurements are proportions. Each composition can be considered as an equivalent class made of proportional factors since the ratios between its components do not change when multiplying it by a positive constant. In this case, if x, y are compositions and x = λy for some constants λ, we say that x, y are compositionally equivalent. To check the equivalency of the two compositions, we can use the closure function C(x), defined as

C(x) = κ • x 1 p i=1 x i , κ • x 2 p i=1 x i , . . . , κ • x p p i=1
x i where κ > 0 is a fixed constant; in this definition, two p-part compositions x, y are compositionally equivalent if C(x) = C(y). The sample space of CoDa is the simplex,

S p = x = (x 1 , x 2 , . . . , x p ) | x i > 0, i = 1, . . . , D; p i=1 x i = κ
In R p space, we can use Euclidean geometry to add vectors or multiply vectors by scalar to obtain their properties or compute their distance. But, due to special structure of CoDa vectors in S p , this geometry can not be applied directly. [START_REF] Aitchison | The Statistical Analysis of Compositional Data[END_REF] introduced the Aitchison geometry, with two operations required for a vector space structure on S p : Perturbation and powering operators. The perturbation ⊕ of x ∈ S p by y ∈ S p (equivalent to the addition in R p ) is defined by

x ⊕ y = C(x 1 y 1 , . . . , x p y p ) ∈ S p and the powering ⊙ of x ∈ S p by a constant α ∈ R (equivalent to the multiplication by a scalar operation in the R p ) is defined by α ⊙ x = C(x α 1 , . . . , x α p ) ∈ S p In practice, CoDa are often transformed to vectors in the Euclidean space to remove its constraints. The center logratio (clr) transformation of vector x ∈ S p , clr(x), is an isometry from S p to a subspace U ⊂ R p , defined by clr(x) = ln

x 1 g m (x)

, ln x 2 g m (x) , . . . , ln x p g m (x) = (ξ 1 , ξ 2 , . . . , ξ p ) where

g m (x) = p i=1 x i 1 p = exp 1 p p i=1 x i
is the geometric mean of the composition and

p i=1 ξ i = 0.
The inverse center log-ratio clr -1 (ξ) recovering [START_REF] Egozcue | Isometric logratio transformations for compositional data analysis[END_REF] showed that the constraint in the component of clr(x) makes singular the clr(x) variancecovariance matrix for random composition. To overcome this drawback, [START_REF] Egozcue | Isometric logratio transformations for compositional data analysis[END_REF] proposed a new transformation which is associated with an orthogonal basis in S p , named isometric log-ratio (irl) transformation. Let e 1 , e 2 , . . . , e p-1 be an orthonormal basis of S p . Any composition x ∈ S p can be expressed as

x from ξ = (ξ 1 , . . . , ξ p ) is clr -1 (ξ) = C(exp(ξ)) = C(exp(ξ 1 ), exp(ξ 2 ), . . . , exp(ξ p )).
x = p-1 i=1 x * i ⊙ e i , x * i = ⟨x, e i ⟩ a = ⟨clr(x), clr(e i )⟩
where ⟨, ⟩ a denotes the Aitchison inner product. Thus, the ilr transformation of x ∈ S p is ilr(x) = x * = (x * 1 , x * 2 , . . . , x * p-1 ). Let B be a (p -1, p) matrix whose i th row is clr(e i ), i = 1, . . . , p -1. This matrix is known as a contrast matrix associated with the orthonormal basis e 1 , e 2 , . . . , e p-1 . The ilr transformation x * of composition x can be computed by

x * = ilr(x) = (x * 1 , . . . , x * p-1 ) = clr(x) • B ⊺ There are many candidates for an orthonormal basis in S p . Egozcue and Pawlowsky-Glahn 2005 proposed a sequential binary partition to define an orthonormal basis. In this basis, e i is defined to be C(e i,1 , . . . , e i,j , . . . , e i,p ) where

e i,j =                exp 1 i(i + 1) if j ≤ i exp - i i + 1 if j = i + 1 1 otherwise
The elements e i,j are called the balancing elements of this basis. Thus, if this type of orthonormal basis is chosen to transform x, i.e., ilr(x) = clr(x) • B ⊺ , the coordinates x * i are called balances and can be obtained by From its ilr coordinate x * , x can be recovered by using the inverse of ilr transformation:

x * i = i i + 1 ln    i j=1 x j 1 i x i+1    .
ilr -1 (x * ) = clr -1 (x * B) = C(exp(x * B))
. Table 1 illustrates the application of ilr transformation in practice for the case p = 4. The first 4 columns present the components of 6 compositions in S 4 and the remaining 3 columns present their corresponding ilr coordinates in R 3 . As can be seen, these ilr coordinates x * i are not constrained any longer. For more detail on CoDa and its properties, see Pawlowsky-Glahn et al. 2015.

VSI MULTIVARIATE EWMA CONTROL

CHART FOR COMPOSITIONAL DATA

VSI MEWMA-CoDa control chart

Let us suppose that, at each sampling period i = 1, 2, . . ., a sample of size n independent p-part composition observations {X i,1 , . . . , X i,n }, X i,j ∈ S p , j = 1, . . . , n is collected, and suppose also that each X i,j , j = 1, . . . , n, follows a multivariate normal distribution N S p (µ, Σ) on the simplex S p , where µ ∈ S p is the center of compositions and Σ is their variance-covariance matrix. Assume that, when the process is in-control, the composition center is µ 0 and when the process is out-of-control, the composition center is µ 1 . The aim of this paper is to design a variable sampling interval MEWMA control chart (denoted by VSI-MEWMA-CoDa) to monitor the center µ of a p-part compositional process. Since CoDa data has a constant constraint on its components, the traditional VSI-MEWMA control chart may not perform well on monitoring this type of data. In [START_REF] Tran | Monitoring compositional data using multivariate exponentially weighted moving average scheme[END_REF], instead of directly monitoring the composition center µ, the authors proposed to monitor the mean vector µ * = ilr(µ) using a FSI MEWMA control chart for the sample mean coordinates vector X * i . In this study, we will apply the idea of [START_REF] Tran | Monitoring compositional data using multivariate exponentially weighted moving average scheme[END_REF] to investigate a VSI-MEWMA control chart for monitoring a compositional process.

Let {X * i,1 , . . . , X * i,n } be the corresponding ilr coordinates of {X i,1 , . . . , X i,n }, i.e. X * i,j = ilr(X i,j ) ∈ R p-1
. Since X i,j follows a multivariate normal distribution N S p (µ, Σ) on S p , its corresponding ilr coordinate X * i,j follows a multivariate normal distribution N R p-1 (µ * , Σ * ) on R p-1 , where µ * = ilr(µ) ∈ R p-1 is the mean vector, Σ * is the (p -1, p -1) variance-covariance matrix of the ilr transformed data. The values of parameters µ * and Σ * depend on the particular choice of matrix B chosen in ilr transformation (see Pawlowsky-Glahn et al. 2015 and section 2). Denote the ilr coordinates of in-control composition center µ 0 and out-of-control composition center µ 1 are µ * 0 and µ * 1 , respectively. The average of n independent p-part compositional observations is defined by

Xi = 1 n ⊙ (X i,1 ⊕ • • • ⊕ X i,n ) then its ilr coordinate X * i is X * i = ilr( Xi ) = 1 n (ilr(X i,1 ) + • • • + ilr(X i,n )) = 1 n (X * i,1 + • • • + X * i,n ) ∈ R p-1 .
We first recall the FSI MEWMA-CoDa control chart proposed by [START_REF] Tran | Monitoring compositional data using multivariate exponentially weighted moving average scheme[END_REF] as follows. Let the MEWMA vector W i be W

i = r( X * i -µ * 0 ) + (1 -r)W i-1 , i = 1, 2, .
. . where Y 0 = 0, r ∈ (0, 1] is a fixed smoothing parameter. In FSI MEWMA-CoDa control chart, [START_REF] Tran | Monitoring compositional data using multivariate exponentially weighted moving average scheme[END_REF] suggested to monitor the statistic

Q i = W ⊺ i Σ -1 Wi W i , i = 1, 2, . . . (1
) where Σ Wi is the variance-covariance matrix of W i . In this work, the asymptotic form of the variance-covariance matrix

Σ Wi Σ Wi = r n(2 -r) Σ *
was used to compute the plotted statistic (and it is also used in our work). An out-of-control signal is issued when

Q i > U CL = H
, where H > 0 is chosen to achieve a specific value of in-control AT S.

In the FSI MEWMA-CoDa control chart, the sampling interval is a fixed constant h F . As for the VSI MEWMA-CoDa control chart, based on the current value of Q i , the time between two successive samples Xi , Xi+1 is allowed to varied. In this chart, the control limit UCL is held the same as in the FSI chart, and an additional warning limit w = U W L (0 < U W L < U CL) is introduced to determine the switch between the long and short sampling intervals: The long sampling intervals h L is used when the control statistic Q 2 i ≤ U W L 2 (safe region) and the short sampling intervals h S is used when U W L < Q 2 i ≤ U CL 2 (warning region). An out-of-control signal is issued when

Q 2 i > U CL 2 .

Markov chain model

Suppose that the occurrence of an assignable cause makes the in-control composition center µ 0 is shifted to µ 1 , or equivalently µ * 0 is shifted to µ * 1 . In this subsection, we will discuss a method based on the Markov chain model to compute the average of the zero-state time to signal (ATS) for the VSI MEWMA-CoDa control chart. Let AT S 0 , AT S 1 denote the ATS when the process runs in-control, and out-of-control, respectively. In comparison with other control charts, it is desirable to design a chart with smaller AT S 1 while their AT S 0 are the same. In the FSI chart, since the sampling interval h F is fixed, we have AT S FSI = h F × ARL FSI . In the VSI chart, since the sampling interval is allowed to vary, the relation between AT S and ARL would be:

AT S VSI = E(h) × ARL VSI . where E(h) denote the average sampling interval. [START_REF] Lowry | A Multivariate Exponentially Weighted Moving Average control chart[END_REF] showed that the performance of a MEWMA-X chart is a function of the n, µ * 0 , µ * 1 and Σ * only through the non-centrality parameter δ where

δ = n(µ * 1 -µ * 0 ) ⊺ (Σ * ) -1 (µ * 1 -µ * 0 ).
Without loss of generality, we can assume n = 1, µ * 0 = 0 (i.e. the in-control composition center is µ

0 = ( 1 p , 1 p , . . . , 1 p
)) and Σ * = I p-1 (the identity matrix in R p-1 ). In this case, the statistic Q i in (1) is modified to

Q i = b ∥W i ∥ 2 2 with b = 2 -r r
. Consequently, the control limits UCL and UWL of VSI MEWMA-CoDa are modified to be U CL = H/b, U W L = w/b To calculate the in-and out-of-control ATS of the VSI MEWMA-X chart, Lee 2009 modified the Markov chain approach proposed by Runger and Prabhu 1996 to approximate its calculation based on the statistic

q i = ∥W i ∥ 2 .
Concerning the in-control case, the one dimensional Markov chain can be used to approximate ATS. In this case, the interval [0, U CL ′ ], where U CL ′ = H/b, is divided into m + 1 sub-intervals (states): the first subinterval has length g 2 and the others have length g, where

g = 2U CL ′ 2m + 1
. The probability of transition from state i to state j, denoted by p(i, j), is given by • for j = 0,

• for i = 0, 1, . . . , m and j = 1, 2, . . . , m, p(i, j) = P (j -0.5)g r 2 < χ 2 (p -1, c) < (j + 0.5)g r
p(i, 0) = P χ 2 (p -1, c) < g 2r 2 .
Let P 1 denote the (m + 1, m + 1) transition probability matrix corresponding to the transient states with the elements p(i, j) then the zero-state in-control AT S of the VSI MEWMA-CoDa control chart is obtained by

AT S = s ⊺ (I m+1 -P 1 ) -1 h,
where s is the (m + 1)-starting probability vector, i.e. s = (1, 0, 0, . . . , 0) ⊺ , h is the (m + 1)-vector of sampling interval with the i th component h i is defined by

h i = h L if ig ≤ U W L h S if ig > U W L .
The expected sampling interval E(h) is calculated by

E(h) = s ⊺ (I m+1 -P 1 ) -1 h s ⊺ (I m+1 -P 1 ) -1 1 m+1
,

where 1 m+1 = (1, 1, . . . , 1) ⊺ is the m + 1 column vector of 1's.

To calculate the zero-state ATS of VSI MEWMA-CoDa chart in the out-of-control case, Lee 2009 modified the two dimensional Markov chain approach which is originally proposed by Runger and Prabhu 1996. In this approach,

W i ∈ R p-1 is partitioned into W i1 ∈ R with mean δ ̸ = 0 and W i2 ∈ R p-2 with zero mean. Then, q i = ∥W i ∥ 2 = W 2 i1 + W ⊺ i2 W i2 .
The transition probability h(i, j) of W i1 from state i to state j is used to analyze the out-of-control component.

Applying the Markov chain-based approach with the number of states of the Markov chain is 2m 1 + 1, for i, j = 1, 2, . . . , 2m 1 + 1, we have

h(i, j) =Φ -U CL ′ + jg 1 -(1 -r)c i r -δ -Φ -U CL ′ + (j -1)g 1 -(1 -r)c i r -δ
where Φ denotes the cumulative standard normal distribution function, c i = -U CL ′ + (i -0.5)g 1 is the center point of state i with the width of each state

g 1 = 2U CL ′ 2m 1 + 1 .
Concerning W i2 component, the transition probability v(i, j) from state i to state j is used to analyze the in-control component. In this case, the Markov chain approach as in in-control case will be applied with p -2 replacing p -1. The control region is partitioned into m 2 + 1 sub-intervals (states) with the width of each states

is g 2 = 2U CL ′ 2m 2 + 1
. The transition probability v(i, j) is given as follows

• for i = 0, 1, 2, . . . , m 2 and j = 1, 2, . . . , m 2 v(i, j) = P (j -0.5)g 2 r 2 < χ 2 (p -2, c) < (j + 0.5)g 2 r 2
,

where c = (1 -r)ig 2 r 2 , • for j = 0, v(i, 0) = P χ 2 (p -2, c) < g 2 2r 2
Let H denote the (2m 1 + 1, 2m 1 + 1) transition probability matrix of W i1 with elements h(i, j), V denote the (m 2 + 1, m 2 + 1) transition probability matrix of ∥Y i2 ∥ 2 with elements v(i, j), and P 2 denote the transition probability matrix of two dimensional Markov chain. Since W i1 and Y i2 are independent, we have P 2 = H ⊗ V, where ⊗ is the Kronecker's matrices product. Matrix P 2 will consist of the transition probabilities of all transient and some absorbing states of the Markov chain.

Let T be the (2m 1 +1, m 2 +1)-matrix with element T (α, β) given by

T(α, β) = 1 if state (α, β) is transient 0 otherwise
and P be the transition probability matrix containing only transient states of the Markov chain. Then, we have P = T(α, β) ⊛ P 2 where symbol ⊛ indicates the elementwise multiplication of matrices.

Let h be the (2m 1 +1)•(m 2 +1) vector of sampling intervals for the bivariate chain. Lee 2009 defined h to be h ⊺ = (1, 0), . . . , (1, m 2 ), (2, 0), . . . , (2, m 2 ), . . . , . . . , (2m 1 + 1, 0), . . . , (2m 1 + 1, m 2 ) with the element h(i, j) defined by

h(i, j) =    h L if a i,j ≤ U W L 2 h S if U W L 2 < a i,j ≤ U CL 2 0 otherwise
where a i,j = (i -(m 1 + 1)) 2 g 2 1 + j 2 g 2 2 . Thus, the zero-sate out-of-control AT S of VSI MEWMA-CoDa control chart is defined by AT S = s ⊺ (I -P) -1 h where s is the initial probability vector with the component corresponding to state (α, β) = (m 1 + 1, 0) is equal to one and all other components are equal to zero. In case [START_REF] Lee | Optimal statistical design of a multivariate EWMA chart based on ARL and MRL[END_REF] showed that the entry corresponding to the component with value equal to 1 of s is the (m(m + 1) + 1)th entry. Concerning the performance of the program used for the computation of the AT S, we follow the recommendation in [START_REF] Tran | Monitoring compositional data using multivariate exponentially weighted moving average scheme[END_REF] and decide to use m 1 = m 2 = 30.

m 1 = m 2 = m,

Optimization procedure

Assume that the fixed sampling interval in FSI control charts is to be a time unit, i.e. h F = 1. Hence, AT S F SI 0 = ARL 0 . In order to evaluate the performances of VSI MEWMA-CoDa with its FSI version, we can compare their out-of-control AT S 1 while constraining the same incontrol values of both AT S 0 and E 0 (h) (average sampling interval). Thus, the VSI MEWMA-CoDa control chart can be designed by finding the optimal combination of parameters that minimize the out-of-control AT S 1 subject to the predefined constraint of AT S 0 and E 0 (h).

In general, a fixed couple (h S , h L ) is typically used, which can be chosen from the suggested list as in the work of [START_REF] Castagliola | Monitoring the coefficient of variation using a variable sampling interval control chart[END_REF]. However, as discussed in the study of [START_REF] Nguyen | Monitoring the ratio of two normal variables using variable sampling interval ewma control charts[END_REF], while h S is quite reasonable to fix, it seems not practical to fix h L due to the fact that when the control statistic falls into the central region, the process is still in safe and the next sampling interval can be flexible to choose if it does not influence the performance of the chart. Based on this reason, we follow the suggestion in [START_REF] Nguyen | Monitoring the ratio of two normal variables using variable sampling interval ewma control charts[END_REF] to fix the proportion between the UCL and UWL values. Let R be the number such that U W L = R • U CL. When the control limit UCL is determined, the warning limit UWL can be computed based on the value of R.

Thus, the optimal design of the VSI MEWMA-CoDa control chart will consist of searching the optimal parameters (r, H, h L ) which minimize the out-of-control AT S 1 for given shift δ subject to constraints in the in-control AT S 0 and E 0

(h) = 1, i.e, (r * , H * , h * L ) = argmin (r,H,h L ) ATS(n, r, H, R, p -1, δ, h L , h S ) subject to the constraint ATS(n, r * , H * , R, p -1, δ = 0, h * L , h S ) = ATS 0 E 0 (h) = 1
By fixing the in-control predefined AT S 0 value, these optimal parameters can be obtained by using the two-steps optimization procedure as follows (1) Find the set of triples (r, H, h L ) such that the incontrol AT S = AT S 0 and E 0 (h) = 1. (2) Among these feasible triples (r, H, h L ), choose (r * , H * ) which provides the smallest out-of-control AT S value for a particular shift δ in vector µ * 0 . As noted in [START_REF] Tran | Monitoring compositional data using multivariate exponentially weighted moving average scheme[END_REF], the value of r must not be too small to avoid unreliable results and the diverging ability in the Markov Chain approach. In this paper, we fix the minimal bound to search for the smoothing parameter r to be 0.05, as recommended in many studies, including [START_REF] Tran | Monitoring compositional data using multivariate exponentially weighted moving average scheme[END_REF].

PERFORMANCE OF THE VSI MEWMA-CODA CONTROL CHART

In this section, we will compare the performance of the VSI MEWMA-CoDa chart with the FSI MEWMA-CoDa chart proposed by [START_REF] Tran | Monitoring compositional data using multivariate exponentially weighted moving average scheme[END_REF]. The comparison will be based on the values of out-of-control ATS 1 while constraining on the same in-control values of both ATS 0 and E 0 (h). To take advantage of the results from the study of [START_REF] Tran | Monitoring compositional data using multivariate exponentially weighted moving average scheme[END_REF], save the calculation costs, and simplify the application in practice, we propose to find the near-optimal values to the VSI MEWMA-CoDa control chart as follows:

• For each optimal couple (r * , H * ) in Table 2 in study of [START_REF] Tran | Monitoring compositional data using multivariate exponentially weighted moving average scheme[END_REF], the value of UWL and h L are chosen to achieve predefined ATS 0 and E 0 (h), • After obtaining UWL and h L , together with the corresponding (r * , H * ), we compute the ATS 1 of VSI MEWMA-CoDa for specific shift sizes δ and compare them with ARL 1 of FSI MEWMA-CoDa chart (Table 3 in [START_REF] Tran | Monitoring compositional data using multivariate exponentially weighted moving average scheme[END_REF].

The procedure to find the near-optimal values is implemented based on following scenarios:

• n = 1, p = 3, ATS 0 = 200, and E 0 (h) = 1;

• δ ∈ {0.25, 0.50, 0.75, 1.00, 1.25, 1.50, 1.75, 2.00};

• h S ∈ {0.1, 0.5}.

The values ARL 1 (FSI column) of MEWMA-CoDa and ATS 1 (VSI columns) for some different scenarios are shown in Table 2. The values of w such that U W L = w/b and h L to obtain the near-optimal value are also provided for each scenario. Some remarks can be drawn from this results as follows

• The VSI MEWMA-CoDa control chart always outperforms the FSI MEWMA-CoDa control chart in detecting the process shifts. For example, when δ = 0.25, h S = 0.1, we have ARL 1 = 64.6 for FSI MEWMA-CoDa chart and ATS 1 = 56.8 for VSI MEWMA-CoDa chart, • The VSI MEWMA-CoDa charts with smaller h S (h S = 0.1) perform better than the ones with larger h S (h S = 0.5). For example, when δ = 0.5, we have ATS 1 = 19.9 in case h S = 0.1 and ATS 1 = 23.5 in case h S = 0.5, • When the shift sizes δ are large (δ ≥ 1.75), the performance of VSI MEWMA-CoDa chart are still better than FSI MEWMA-CoDa chart, but not much. (1.0, 1.8) 12.9 1.00 9.9 (2.9, 1.3) 6.9 (0.9, 1.8) 8.4 1.25 7.1 (1.6, 1.8) 4.9 (0.9, 2.0) 6.3 1.50 5.4 (3.5, 1.2) 3.7 (0.9, 1.9)

CONCLUSION

In this paper, we proposed a VSI MEWMA-CoDa control chart to monitor a normal multivariate random vector defined as the inverse isometric log-ratio of a p-part composition. The optimal procedure to compute the optimal triple (r * , H * , h * L ) and the ATS values of the proposed chart for different shift sizes were presented. We also proposed a method to find the near-optimal values for the VSI MEWMA-CoDa chart to utilize the results in the study of [START_REF] Tran | Monitoring compositional data using multivariate exponentially weighted moving average scheme[END_REF] and reduce the computation costs. The numerical performance comparison between the VSI MEWMA-CoDa chart and standard (FSI) MEWMA-CoDa control chart in terms of ATS 1 (based on the nearoptimal values method) showed that the VSI MEWMA-CoDa chart always outperforms the standard chart. Future research on monitoring CoDa could be concentrated on the extension of the VSI MEWMA-CoDa chart to the VSI MCUSUM-CoDa chart, or investigating the effect of measurement error on these charts. The methods to transform CoDa into normal data before designing these controls charts are also worthy to focus. Due to the wide applications of CoDa in the real-life, the online monitoring of CoDa should be worthy of consideration by researchers in the SPC field.

  type control chart for monitoring CoDa data was proposed. Recently, Vives-Mestres et al. 2014b investigated a T 2 control chart for monitoring CoDa with p = 3 and then Vives-Mestres et al. 2014a extended the work in Vives-Mestres et al. 2014b for individual observations case. Two methods for interpretations of out-of-control signal of individual T 2

2

  where χ 2 (p -1, c) denotes a non central chi-square random variable with p -1 degrees of freedom and non-centrality parameter c = (1 -r)ig r 2 ,

Table 1 .

 1 An example of ilr transformation in S 4

	x 1	x 2	x 3	x 4	x * 1	x * 2	x * 3
	0.10 0.30 0.50 0.10 -0.78 -0.87	0.78
	0.20 0.25 0.20 0.35 -0.16	0.09	-0.42
	0.50 0.10 0.20 0.20	1.14	0.09	0.06
	0.60 0.05 0.05 0.30	1.76	1.01	-0.83
	0.35 0.15 0.10 0.40	0.60	0.68	-0.72
	0.20 0.45 0.05 0.30 -0.57	1.46	-0.52

Table 2 .

 2 Comparison between VSI MEWMA-CoDa and FSI MEWMA-CoDa charts

	δ	FSI	h S = 0.1 (w, h L ) VSI	h S = 0.5 (w, h L ) VSI
	0.25	64.6	(1.7, 1.6)	56.8	(0.7, 2.1)	63.5
	0.50	26.4	(1.7, 1.6)	19.9	(0.9, 1.8)	23.5
	0.75	15.1	(1.6, 1.7)	10.4


Santosh K Misra

Paturu Kondaiah

Santanu Bhattacharya

Didier Boturyn

Pascal Dumy

Co-liposomes comprising a Lipidated Multivalent RGD-Peptide and a Cationic Gemini Cholesterol Induce Selective Gene Transfection in αvβ3 and αvβ5 Integrin Receptor-rich Cancer Cells

HAL is

Introduction

Selective cell targeting has been considered as a major objective in gene therapy to combat different diseases. [START_REF] Barry | [END_REF] The gene carrier species should have specific ligands which can help it to recognize the victim cell, so that it can deliver the cargo to the desired cells only. The concept of cell targeting introduces an effort to recognize natural targets on tumor cell surfaces. Integrin receptors have been recognized as natural targets for such selective transfection. 2 Integrins belong to a major class of αβ heterodimeric transmembrane glycoprotein receptors. [3][4][5] It is known that αvβ3 and αvβ5 integrin are expressed on various cell types such as endothelial cells, 6 osteoclasts, 7 macrophages, 8 platelets, 9 and melanomas. 10 They play a significant role in angiogenesis, vascular intimal thickening and proliferation of malignant tumors. 11 Some important cancer cell lines expressing high  v  3 integrin are A375 (melanoma), M21 (melanoma), SK-MEL-28 (melanoma), M19-MEL (melanoma) while Hela cells (carcinoma) over-expresses v3 as well as v5. 12 However, some cell lines do not express v3 and v5 integrins. These include HEK293 13 293T 14 and NIH3T3 15 cells, which express another integrin receptor as 51 in abundance.

Cells expressing αvβ3 or αvβ5 integrins bind with a broad range of ligands such as fibronectin and vitronectin mainly through the recognition of the ubiquitous triad sequence RGD (Arg-Gly-Asp). 16 The latter has served as a basis for the development and discovery of high and selective integrin-peptide ligands such as c[-RGDfV-], 17 that was exploited for a wide range of applications. [18][19][20][21][22][23] Additionally, the principle of multivalency has then been recognized as an important approach for the design of synthetic ligands. 24 In this context, 50 enhancements of biological activity are obtained from multivalent RGD peptide ligands used to target cell surface receptors such as αvβ3 integrin. [25][26][27] For instance, we have shown that tetrameric RGD-containing scaffolds 28 exhibit advantageous biological properties for imaging [29][30][31] and for targeted drug delivery. 32,33 55 Two strategies may be conceived for achieving cell selective targeting using carrier liposomes. Firstly the covalent conjugation of RGD based molecules to the lipid carrier 34 and secondly the non-covalent loading of RGD based molecules to lipid carriers by physical mixing in optimized molar ratios. 35 Very recently, we 60 prepared RGD-containing lipopeptide to construct fluid supported lipid bilayers (SLB). 36 We then studied RGD-mediated cell adhesion to an SLB surface providing the average ligand spacing to trigger cell adhesion and spreading on a fluid substrate. Herein we present the results of non-covalent RAFT-RGD inclusion in

Results

The concept of cell targeting introduces an effort to recognize natural targets viz. integrin receptors for selective transfection. The multivalency has been established as an important approach for the design of ligands to increase the targeting efficiency. We 20 applied the non-covalent palmitoylated RAFT-RGD inclusion strategy to develop a co-liposomal formulation of GL, DOPE and a palmitoylated RAFT-RGD macromolecule 1 to transfect pEGFP-C3 and PGL-3 in αvβ3 and αvβ5 positive HeLa and negative HEK293 and NIH3T3 cells. This methodology of using 25 lipopeptide 1 was found highly beneficial for introducing high selectivity in gene therapy.

Synthesis of Palmitoyl-RAFT-RGD 4

The synthesis of lipopeptide 1 was produced via solid-phase peptide syntheses (SPPS) according to the methods already 30 developed by our group (Scheme 1). Stable oxime bonds were chosen to append aldehyde-bearing RGD motifs to a cyclodecapeptide scaffold. Building blocks 2 and 3 including respectively the protected amino-oxy function and the palmitoyl group were directly introduced during the SPPS. This strategy 35 decreases the number of steps and the combination of protecting groups required so far for the production of lipopeptide 4. The following head-to-tail cyclization on resin and the subsequent deprotection under acidic conditions furnished the lipopeptide 6. The grafting of glyoxylyl aldehyde containing-RGD targeting 40 elements 7 provided the expected palmitoyl-RAFT-RGD 4 1 in satisfying 15% overall yield.

Optimization of GL: DOPE ratio

Gemini lipid GL5 was screened for the optimization of the lipid: DOPE ratio from 1:0 to 1:5. Here the number of transfected cells 50 (% GFP cells) was higher for the GL: DOPE at ratios 1:3 and 1:5 compared to that in 1:4. However, the amount of DNA transfected to each cell (MFI) was very high when a molar ratio of 1:4 was used. A gemini lipid to DNA ratio of 1:4 was used as the optimal ratio for further studies (Figure S1).

55

Optimization of Palmitoyl-RAFT-RGD 4 inclusion in liposomes

Incorporation of the Palmitoyl-RAFT-RGD 4 1 into cationic liposomes was achieved either by the preparation of co-liposomes or via post-injection method. In the first method, an optimized 60 GL5D was mixed with the lipidated RGD in molar ratios of 2000:1 and 1000:1, respectively. Liposomes were prepared from these mixtures as mentioned earlier and were termed as coliposomal formulations, i.e., Co2000 and Co1000, respectively. A second method was used by injecting an ethanolic solution of the 65 lipidated RGD in preformed GLD suspensions. These formulations have been termed as post-injection formulations, i.e., Po2000 and Po1000, respectively. FACS analysis of the pEGFP-C3 reporter gene transfection showed a distinct advantage of the co-liposomal preparation over the one prepared 70 by the post-injection protocol. Liposome GL5D-RGD was prepared by co-liposomal method showed at least 50% increase in the transfection efficiency compared to the corresponding gemini lipid-DOPE-RGD formulations prepared by the post injection method (Figure S2A).

75

Optimization of GLD: RAFT-RGD molar ratio

Molar ratios of GL5D: RGD were varied from 1:0 to 50:1 for the pEGFP-C3 transfection in HeLa cells. A molar ratio of 100:1 was found to be the best ratio for achieving high and cell-specific transfection. Beyond this ratio, the transfection efficiency the alteration in the liposomal organization and/or aggregation leading to the formation of significantly larger aggregates of the GL5D-RGD formulations at GL5: RGD molar ratio of 50:1 (Figure S2C)

Palmitoyl-RAFT-RGD 4 mediated cell targeting

Optimized GLD and GLD-RGD formulations were used for transfecting three different cell lines, HeLa, HEK293 and NIH3T3 using 0.8 μg pEGFP-C3 per well using the GL/DNA charge ratio (N/P) 0.5 (Figure 2). In HeLa cells, GL5D induced ~50% GFP cell transfection with a MFI of ~25 only while GL5D-10 RGD furnished ~75% GFP cells with a MFI of ~75 at the N/P charge ratio of 0.5. In HEK293 cells, GL5D-RGD gave ~40% GFP cells with a MFI of ~20 at the N/P charge ratio of 0.5. Lipofectamine2000® in presence of 10% FBS (Lipo2000*) was used as positive control in all the cell lines. GL5D-RGD was found to be better than Lipo2000* at least in HeLa and NIH3T3 cells. We further measured the luciferase activity with respect to N/P charge ratios and the relevant data are shown as RLU/mg of protein. In HeLa cells, GL5D-RGD gave around 5-fold higher 65 transfection efficiency compared to that of GL5D at the N/P of 0.5. On the other hand in NIH3T3 cells, the overall transfection was found to be strikingly low and GL5D-RGD showed no improvements in poorer transfection efficiency compared to GL5D. Thus GL5D-RGD was found significantly better 70 transfecting agent in HeLa cells compare to GL5D with high biostatistical significance with p < 0.001 (Figure 3).

Gel electrophoresis

DNA complexation efficiency of the GLD and GLD-RGD formulations were examined at different N/P charge ratios 0.125, 0.25, 0.75 and 1. The GL5D and GL5D-RGD gave around 100% complexation at N/P charge ratio of 0.5. Inclusion of Palmitoyl-110 RAFT-RGD 4 (1) did not however affect the DNA binding efficiency of the GLDs (Figure 4). Miceller SDS mediated decomplexation of DNA from the lipoplexes was analyzed using lipoplexes of N/P ratio 3 followed by incubation with SDS at SDS/GL ratio varying from 1 to 4 (Figure 5). At SDS/GL ratio of GL5D/DNA could release only 40% even at SDS/GL of 5. Stability of the lipoplexes in presence of FBS was measured both using a fixed N/P ratio while varying the % of FBS (v/v) (Figure S3A) and by keeping the % of FBS fixed while altering 30 the N/P charge ratio (Figure S3B). It was found that even 50% of FBS could not dissociate the lipoplexes at N/P charge ratios for the liposomes with or without RGD whereas in presence of 70% FBS, an N/P charge ratio of 2 was enough to protect the DNA from being released.

35

Effect of RGD on the lipoplex delivery efficiency

Overall the transfection efficiency of liposomes was found to be related with the DNA delivery efficiency. The UV absorbance of the undelivered DNA was found to be higher in case of HeLa cells when they were transfected with GL5D compared to GL5D-40 RGD. In NIH3T3 and HEK293, however, it was found to be the opposite (Figure S4).

Cell surface binding

A flow cytometry experiment was performed on both the HeLa and NIH3T3 cells to find out the selective cell surface binding 45 efficiency of the lipoplexes with the Palmitoyl-RAFT-RGD 4 (Figure 7). HeLa cells treated with GL5D-RGD: pyrene, selectively expressed ~95% pyrene positive cells compared to ~70% in case of GL5D: pyrene. In NIH3T3, GL5D-RGD: pyrene as well as GL5D: pyrene both expressed ~65% pyrene positive cells without any selectivity. Presumably the presence of the RAFT-RGD in GL5D-RGD: pyrene assisted selective and highly efficient surface binding to cells containing the integrins receptors, αvβ3 and αvβ5 (HeLa cells). In contrast, in NIH3T3 cells (which are devoid of αvβ3 and αvβ5) the RAFT-RGD could not influence the surface binding percentage to any significant extent.

Confocal Fluorescence Microscopy

Cells were observed under a confocal microscope to visualize the expression of green fluorescent protein (GFP) in the cell cytoplasm. Green fluorescence was clearly visible in the cell cytoplasm, which was obtained via cell transfection using the cationic gemini cholesterol formulations with or without the RAFT-RGD in both HeLa and NIH3T3 cell lines (Figure 8) similar to the cells transfected with the commercially available agents Lipofectamine2000® and Effectene.

Cytotoxicity

The MTT assays of optimized GLDs and GLD-RGDs, with and without DNA were obtained in HeLa, NIH3T3 and HEK293 cell lines (Figure 9). All formulations without DNA were screened for concentrations 0.75 μM to 6 μM while formulations with DNA were screened for N/P charge ratios varying from 0.125 to 1. In the case of the final optimized formulations, the % cell viability was found to be ~80% which dropped to ~60% at higher concentrations and at higher N/P charge ratios, specifically in the case of other non-optimized formulations across each cell line. Importantly an inclusion of addition of RAFT-RGD in the liposomal formulation did not reduce it at all, in any cells, e.g., HeLa, NIH3T3 and HEK293 examined. There are many barriers that impede effective transfection using non-viral gene vectors. Stability of vector/DNA complexes in the bloodstream and specific cellular uptake by the target tissues/organs are some of the main obstacles for in vivo gene delivery. Cationic lipids are well known candidates for gene 50 delivery. 38,39 Liposomal formulations of the cationic cholesterol derivatives are known to achieve DNA transfection in cells with greater efficiency than a number of commercial transfection reagents. These cationic cholesterol based formulations are more efficient than the cationic lipid suspensions based on 55 pseudoglyceryl skeleton in term of gene transfection. 40,41 Moreover gemini lipids have been shown to be more efficient than their monomeric counterparts in terms of their ability to transfect eukaryotic cells. 37 Earlier we have reported serum compatible cholesterol based GLs as efficient transfecting 60 agents. 42,43 But to establish a precisely targeted gene delivery, there was a need to design and develop newer formulations with efficient cell targeting moiety also capable of efficient gene transfection. Accordingly we have utilized the cationic cholesterol derivatives as cytofectins for our studies. The spacer

6 |J. Mater. Chem. B, [2014], [vol], 00-00
This journal is © The Royal Society of Chemistry [year] lengths between the two monomeric units in such gemini lipids play important role in determining the transfection efficacies of the cationic lipids. 38,41 Recognition of specific integrin receptors on cell surfaces by the particular tripeptidic sequence (RGD), allows the binding of 5 RGD loaded liposomal formulations with cell surfaces having abundance of compatible integrin receptors. One can employ different ways of using RGD templates in gene delivery formulations for achieving integrin receptor specific targeting; firstly via a covalent modification of the cationic lipid molecules 10 with RGD based tripeptides 34 or secondly via a non-covalent incorporation of an optimized amount of RGD based lipid soluble molecules. To avoid the complexity of molecular geometry and in order to achieve convenient variation of ratios of RGD to GLs, we decided to proceed with the non-covalent strategy. 15 For the present work, cationic gemini cholesterol (GL5) was used with DOPE at defined proportions to form liposome GL5D, which also contained specified amount of palmitoylated multivalent RGD to prepare GL5D-RGD via co-liposomal and post-injection methods. During the co-liposome formation, 20 presumably the long fatty acid chain in the palmitoylated multivalent RGD ligand has full freedom to be anchored into the hydrophobic region of the vesicular aggregates. Also in case of the palmitoyl-RAFT-RGD 4 , the presence of the sterically demanding cyclic peptide residues might ensure preferential 25 partitioning of these molecules on the outer periphery of the coliposomes. It is certainly the reason for the increase of the transfection efficiency by using the co-liposomal methodology.

Co-liposomal formulation was used for the optimization of the GLD: RGD molar ratio for highest transfection with best 30 selectivity. It was found that 100:1 molar ratio was the best transfecting formulation. Incorporation of RGD templates in molar ratios more or less than 100:1 led to decreases in the transfection efficiency (Figure 2B). Probably, at lower ratios (2000:1 or 1000:1), the number of RGD tags onto the co-35 liposomes was insufficient to achieve maximum while at a higher ratio (50:1) too many RGD tags might have altered the surface properties of the co-liposomes (GLD) in such way that it reduced the transfection efficiency by probably due to the instability and/or aggregation of liposomes with significantly larger 40 hydrodynamic diameter for the GL5D-RGD formulation at GL5: RGD molar ratio 50:1 (Figure S2C)

The overall optimized GL5D-RGD liposomes were then used for transfection of HeLa, HEK293 and NIH3T3 cells using plasmid pEGFP-C3 and PGL-3 vectors. In all the experiments, GL5D, devoid of RGD, were used as negative controls. Transfection of pEGFP-C3 revealed that in HeLa cells GL5D-RGD formulations were significantly better transfectants (Figure 2A) compared to the GL5D (P < 0.001) whereas in NIH3T3 and HEK293 cells, inclusion of RGD did not lead to an 50 enhancement in the transfection efficiency (Figure 2B-C). Similar results were obtained in transfection of PGL-3 plasmids in HeLa and NIH3T3 cells as well. In this instance RGD inclusion in GLDs increased the transfection efficiency by ~5 fold in HeLa cells with P < 0.001 in bio-statistical analysis (Figure 3A). In 55 contrast in NIH3T3 it did not change to any considerable extent (Figure 3B).

An introduction of RAFT-RGD in the GL5D formulation at N/P charge ratio of 0.5 causes significant improvement in HeLa cell selectivity as revealed from the FACS analysis. Similarly the 60 luciferase assay shown in Figure 3 indicates that an introduction of RAFT-RGD in GL5D improves the transfection efficiency at N/P charge ratios of 0.5 in HeLa cells only. It may be further noted here that luciferase assay was performed using PGL-3 plasmid while plasmid pEGFP-C3 was used for the FACS 65 analysis. It is possible that either an improved DNA binding, a better DNA release, an improved serum tolerance or better DNase I sensitivity of the formulations in presence of RGD might have resulted in a high transfection efficiency with specific targeting. 70 But gel electrophoresis experiments performed on formulations with and without RAFT-RGD indicated that the presence of RGD did not enhance the binding of DNA (Figure 4). This also did not enhance the SDS mediated release of DNA (Figure 5) or even afford better DNA protection against DNase I (Figure 6). Even 75 the stability of lipoplexes in presence of serum was not affected by the presence of RGD (Figure S3). These results clearly suggest that there might be another factor, which governs the enhancement of transfection efficiency by the co-liposomes loaded with the RGD. The HeLa cells possess high abundance of 80 αvβ3 and αvβ5 transmembrane integrin receptors while the HEK293 and NIH3T3 cells are devoid of them. As expected, Palmitoyl-RAFT-RGD 4 present in the co-liposomal formulations improves the transfection efficiency in HeLa cells.

In support of the RGD-mediated targeting of the αvβ3 and αvβ5 85 integrins, undelivered lipoplex measurement assay was performed. This indicates that in HeLa cells, GLD-RGD could deliver more lipoplexes in the cytoplasm compared to that of the GLD alone (Figure S4A). In contrast, in NIH3T3 and HEK293 cell lines, the observed pattern was just the opposite. Thus with 90 GLD-RGD, higher amount of lipoplexes were left undelivered. Similarly, GLD-RGD: pyrene was able to tag higher % of cells compared to the GLD: pyrene while NIH3T3 showed no enhancement in surface binding (Figure 7), indicating a selective targeting induced by the RGD peptide (Figure S5).

95

Confocal fluorescence microscopy was performed to find out the GFP expression of the pEGFP-C3 plasmid in HeLa and NIH3T3 cells (Figure 8). It was found that each formulation, with or without RGD along with commercially available Effectene and Lipofectamine2000® (Lipo2000*) (two positive controls), were 100 able to transfect the cells. Confocal fluorescence microscopy data shown here indicate a qualitative expression of GFP. This cannot be used here for any quantitative comparison of different formulations including positive controls, Effectene and Lipofectamine2000®. The results obtained from FACS and 105 luciferase assays are however, quantitative for this purpose. Finally, the cytotoxicity of each formulation was measured using the MTT assay protocol. It confirmed the considerably low toxicity of the formulations with the palmitoylated-RAFT-RGD 4 (Figure 9). Thus, the co-liposomal methodology of RGD tagging 110 on cationic liposomes has been demonstrated to be an efficient and convenient way of cell targeting in the gene delivery experiments.

Conclusions

A new gene delivery formulation has been developed which is 115 composed of a cationic gemini cholesterol (GL5), helper lipid (DOPE) and a palmitoyl-RAFT-RGD 4 template to imbibe the selectivity in liposomes towards cells therefore endowed with abundance of αvβ3 and/or αvβ5 integrins. Stable formulations were developed by a non-covalent co-liposomal mixing of GL, DOPE and palmitoylated RGD-RAFT in different optimized molar ratios. Expression of pEGFP-C3 and PGL-3 in αvβ3 and αvβ5 positive cell (HeLa) was found to be considerably high for the GL5D-RGD formulations compared to other αvβ3 and αvβ5 negative cells (HEK293 and NIH3T3). Similar gel electrophoretic 10 patterns of liposomes with (GL5D-RGD) or without (GL5D) RGD confirmed that the RGD played no role in enhancing the DNA binding, anionic micelle induced DNA release, serum stability or even DNA protection against DNase I to improve the transfection. Measurement of the extent of undelivered lipoplexes and % GLD-RGD: pyrene labeled positive cells support the role of RGD in αvβ3 and αvβ5 mediated cell selectivity by the present gemini lipid based transfectants. The formulations were good transfecting agents in comparison to commercially available Effectene or Lipofectamine2000® (Lipo2000*). As inclusion of 20 RGD did not increase the cytotoxicity of the liposomes and corresponding lipoplexes, it may be highly beneficial to use this formulation for in vivo gene delivery processes. We are currently developing this strategy and in vivo evaluation will be reported in due course.

25

Experimental

Materials and Methods

Gemini lipid (GL), GL5 has been synthesized (Supporting Information) as described elsewhere [37]. Palmitoyl-RAFT-RGD 4 was prepared through a combination of solid and solution-30 phase syntheses according to the methods previously described (Supporting Information) [36]. All reagents, solvents, and chemicals used in this study were of the highest purity available. The solvents were dried prior to use. Column chromatography was performed using 60-120 mesh silica gel. NMR spectra were recorded using a Jeol JNM-300 (300 MHz for 1 H and 75 Hz for 13 C) spectrometer. The chemical shifts (δ) are reported in ppm downfield from the internal standard, TMS, for 1 H-NMR and 13 C-NMR. Mass spectra were recorded on a Kratos PCKompact SEQ V1.

Preparation of lipopeptide

Palmitoylated Multivalent RGD-Peptide 1 was prepared as described in the literature (see the supporting information).

Preparation of liposomes

Liposomes were prepared using a literature protocol. 37 In brief, each of the gemini lipid was mixed with DOPE in an appropriate molar ratio as solution in chloroform. The organic solvent was evaporated under rotary evaporator and finally dried under 60 vacuum. The vacuum dried films of lipid mixtures were subjected to hydration with addition of appropriate amount of buffer, repeated freeze-thaw cycles followed by sonication at 70 o C for 15 min. Each liposomal suspension was found to be transparent and very stable while stored at 4 o C in sterile conditions.

65

Formulations were re-sonicated at room temperature for 5 min before use for each experiment.

Plasmid DNA

Plasmids pEGFP-C3 and PGL-3, which encode for an enhanced green fluorescence protein (GFP) under a CMV promoter and 70 firefly luciferase protein under a SV40 promoter respectively, were amplified in Escherichia coli (DH5α) and purified using Qiagen Midi Prep plasmid purification protocol (Qiagen, Germany). The purity of the plasmid was checked by electrophoresis on 1.0% agarose gel. The concentration of DNA 75 was estimated spectroscopically by measuring the absorption at 260 nm and confirmed by gel electrophoresis. Plasmid preparations showing a value of OD 260 /OD 280 >1.8 were used.

Cell culture

Cells (HeLa, HEK293 and NIH3T3) were cultured in Dulbecco's 80 modified Eagle's medium (DMEM, Sigma) supplemented with 10% fetal bovine serum (FBS) in T25 culture flasks (Nunc, Denmark) at 37 °C in a humidified atmosphere containing 5% CO 2 . Cells were regularly passaged by trypsinization with 0.1% trypsin (EDTA 0.02%, dextrose 0.05%, and trypsin 0.1%) in PBS 85 (pH 7.4).

Transfection procedure

Transfection experiments were performed using literature protocol. 36 In brief, 24 h before transfection, 24-well plates were seeded with 60, 000 cells/well in antibiotic-free media.

Transfection experiments were performed in presence of serum (-FBS+FBS) where lipoplex preparation was carried out in absence of serum but the lipoplexes were incubated with the cells in presence of serum. At the time of transfection, cells were at least 70% confluent. Liposomal formulations and plasmid DNA were 95 diluted in DMEM. Lipoplexes were prepared in DMEM by mixing liposome and DNA in required N/P charge ratio and the total volume was made up to 200 μl. The N/P values represent the ratio of charges on the cationic gemini cholesterol (in mol) to nucleotide base molarity and were calculated by considering the 100 average nucleotide mass of 330. Lipid-DNA complexes were incubated for 30 min at room temperature followed by dilution with 200 μl of 20% FBS containing DMEM (final DNA concentration 12.12 μM). Formulations were added to the wells and incubated for 6 h in optimal condition. After 6h of 105 incubation, old medium was replaced with new medium having 10% FBS. Cells were further incubated for a period of 42h before checking the reporter gene expression. The extent of reporter gene expression was examined by fluorescence microscopy and/or luciferase assay technique and expression was quantified either by flow cytometry analysis or by luminometry. All the experiments were performed in duplicate, and the results presented are the average of at least duplicates of two such independent experiments performed on different days. Lipofectamine2000 was used as positive control using 5 manufacturers protocol in 10% serum, abbreviated as Lipo2000* throughout the manuscript.

Flow cytometry

Transfection efficiency of different liposomal formulations for the reporter gene pEGFP-C3 was examined by performing flow 10 cytometry assay (FACS). 44 Transfected cells were regularly observed under fluorescent microscope for the intensity of fluorescence due to green fluorescence protein (GFP) expressed by the cells. Finally, 42h post transfection, cells were trypsinized by adding 100 μl of 0.1% trypsin and cells were collected in 5% FBS containing DPBS. Duplicate cultures were pooled and analyzed by flow cytometry immediately using a Becton and Dickinson flow cytometer equipped with a fixed laser source at 488 nm. FACS data were analyzed using public domain WinMDI2.8 software.

20

Luciferase assay

Transfection efficiency of various formulations for the reporter gene PGL-3 expression was quantified further by luciferase assay. Transfected cells were assayed using single-luciferase assay kit, provided by Promega (U.S.A.). In a typical assay, after 25 48 h of transfection, the old medium was removed from the wells and the cells were washed twice with 200 μl of PBS and 100 μl of 1x cell lysis buffer was added to each well. Cells were lysed for 30 min in a horizontal rocker at RT. The cell lysate was transferred to centrifuge tubes and centrifuged (4000 rpm, RT) 30 for 2 min. Supernatant was stored at -70 o C until used for assay.

For the assay, 10 μl of supernatant was mixed with 10 μl of luciferase assay substrate (Promega). The lysate and the substrate were both thawed to RT before performing the assay. The luciferase activity was measured in a luminometer (Turner 35 designs, 20/20, Promega, U.S.A.) in standard singleluminescence mode. The measurement was performed for 10 s. A delay of 2 s was given before each measurement. The protein concentration was estimated in each case using Bradford's method with bovine serum albumin as a standard. 45 Comparison 40 of the transfection efficiencies of the individual lipids was presented using data for luciferase assay expressed in relative light units (RLU)/mg of protein or luciferase activity/mg of protein.

Statistical analysis 45

Statistical significance of differences between control and samples were evaluated using two-way ANOVA using GraphPad Prizm 5.0 with Bonferroni analysis wherever applicable. Results were considered statistically significant when the p value was less than 0.05 (* represents p < 0.05; ** represent p < 0.005 and *** 50 represent p < 0.001).

Optimization of Lipid: DOPE ratio

DOPE, a natural lipid, is known to enhance gene transfection efficacy of different cationic liposomal formulations, 46 has been used to prepare high output transfecting agents. Gemini lipid:

55 DOPE ratio was varied from 1:0 to 1:5 of molar ratio to optimize the efficiency of transfection. Efficiency of each liposomal formulation was checked by performing pEGFP-C3, reporter gene transfection on HeLa cells at N/P charge ratio of 0.5 in 10% FBS condition.

60

Optimization of Palmitoyl-RAFT-RGD 4 inclusion methodology

Two different methodologies were used for non-covalent incorporation of the palmitoyl-RAFT-RGD 4 in GLD liposomes. In the first method, optimized amount of GLD was mixed with 65 the RGD-containing scaffold in molar ratios of 2000:1 and 1000:1. Liposomes were then prepared from each of these mixtures as mentioned earlier. These formulations were termed as co-liposomal formulations Co-2000 and Co-1000 respectively. A second method was used via injection of an ethanolic solution of 70 RGD-containing scaffold in pre-formed GLD suspensions in buffered media. Formulations were termed as post-injection formulation Po-2000 and Po-1000 respectively. Efficiency of formulations was checked by performing pEGFP-C3, reporter gene transfection on HeLa cells at N/P charge ratio of 0.5 in 10% 75 FBS condition.

Optimization of GLD: Palmitoyl-RAFT-RGD 4 molar ratio

Palmitoyl-RAFT-RGD 4 was mixed with the GL5D formulation in ratios of GLD: RGD-peptide 1:0, 2000:1, 1000:1, 100:1 and 50:1 respectively. GL5D was optimized as 100:1 molar ratio as best 80 with respect to the transfection efficiency. This optimized ratio was used for further optimization of the N/P charge ratio.

Gel electrophoresis

Complexation and de-complexation ability of the GLD and GLD-RGD at different GL/DNA ratios, 47 DNase I stability of 85 lipoplexes derived from GLD/DNA and GLD-RGD/DNA 48 and the stability of lipoplexes in presence of FBS 49 were all assayed by gel electrophoresis. Lipoplexes were prepared at different GL/DNA charge ratios. After 30 min of incubation, lipoplexes were run on 1% agarose gel to perform gel retardation assay. 90 Lipoplexes were incubated for 15 min in presence of micellar SDS at different SDS: GL molar ratio and run on the 1% agarose gel to perform de-complexation assay whereas DNase I (1 unit) was incubated along with lipoplexes for 30 min at 37 o C and run on the 1% agarose gel to find out the DNA protection efficiency 95 of lipoplexes against DNase I. Stability of lipoplexes in high serum percentage was performed by incubating the lipoplexes with high FBS percentage and run on the 1% agarose gel. We also performed gel electrophoresis for the lipoplexes that remained un-internalized to the cell after 6h of incubation. 100 Observed band intensity proportional to the amount of lipoplexes could not enter inside the cells even after 6h. High intensity band in the electrophoretic gel indicates low uptake of lipoplexes by cells.

Effect of Palmitoyl-RAFT-RGD 4 on lipoplex delivery efficiency

To measure the effect of palmitoyl-RAFT-RGD 4 inclusion on the DNA delivery efficiency of liposomes in each of HeLa, HEK293 and NIH3T3 cells, we performed UV absorption spectral measurements for the lipoplexes 50 which could not enter inside 110 the cells even after 6h of incubation. Experiment performed was similar to general transfection experiment using pEGFP-C3 plasmid complexed with GLD and GLD-RGD. After 6h of incubation, culture medium possessing undelivered lipoplexes was removed. Medium was centrifuged at 13000 rpm for 20 min at 4 o C. Centrifuged lipoplexes were further suspended in water and equally divided into two parts. One part was run on 1% agarose gel and the second part was separated using phenol: choloroform. The absorbance due to DNA was measured using a BioRad UV-spectrometer. Absorbance values (A 260 ) were plotted 10 against individual formulations which indicate the amount of lipoplex left outside the cells. Results were confirmed by gel electrophoresis.

Cell surface binding analysis

Experiment was performed to find out the specificity of palmitoyl-RAFT-RGD 4 loaded liposomes toward the HeLa cell line that over-expresses αvβ3 and αvβ5 integrin receptors. A fluorescent probe, pyrene was included in GLD (GLD: Pyrene molar ratio 100:1) with and without the palmitoyl-RAFT-RGD 4 . After 24 h, cells were incubated with formulations possessing 20 either GLD: Pyrene or GLD-RGD: Pyrene for 2h. Cells were then washed with DPBS and collected for FACS analysis. FACS data were analyzed using WinMDI2.8 software.

Confocal fluorescence microscopy

Transfection of pEGFP-C3 reporter plasmid was visualized under a confocal fluorescent microscope and the experiment was performed as described. 51 Cells were cultured in T25 culture flasks, trypsinized and plated in 12-well plates, having autoclaved glass slips in the wells. Then, cells were plated on glass slips in antibiotic free 10% FBS containing DMEM as medium. Cells 30 were grown for 24 h at 99% humidity, 37 o C temperature and 5% CO 2 condition till cell-monolayer gained ~70% confluency. Further, transfection was performed as mentioned above. Plates were then incubated for a period of 42 h before checking under confocal microscopy. Control transfections were performed in each case by using commercially available transfection agent 'Lipofectamine2000' and 'Effectene' using manufacturers protocol in presence (-FBS+FBS) of serum. After 42 h of incubation, all the media were removed from the wells and the cells were washed with DPBS buffer carefully to remove all the 40 cell debris without disturbing the monolayer of cells. Cells were fixed for 10 min with 1 mL of 4% paraformaldehyde in each well. Cells were washed with DPBS for 3x10 min and kept with 1 mL 0.1% Triton-x-100 for 5 min to increase the membrane permeabilization. Cells were washed with DPBS for 3x10 min.

Glass slips having fixed and permealizable monolayer of cells were taken out from the wells and kept on the glass slides and incubated for 5 min with 1 g/ml of PI (propidium iodide) to specifically stain the nucleus of the cells. Finally cells were washed with DPBS for 3x10 min to remove extra PI to reduce over-staining. Vector shield was used to mount the glass slips upside down on the glass slides which were observed under a confocal microscope (Zeiss) at objective of Plan Apochromat 63x/1.4 Oil DIC and analyzed using a Zeiss LSM image browser.

Cytotoxicity

55

Cytotoxicity of GL5D and GL5D-RGD formulations were determined using 3-(4,5-dimethylthiazole-2-yl)-2,5-diphenyl tetrazolium bromide (MTT) reduction assay in HeLa, HEK293 and NIH3T3 cells. 52 In brief, 15,000 cells/well were plated in 96 well plates. After 24 h, GL5D, GL5D-RGD, GL5D/DNA and 60 GL5D-RGD/DNA formulations were added to the cells in absence of serum. Liposomal formulations GL5D, GL5D-RGD were complexed with 0.2 g of the DNA at various N/P ratios for 30 min to give GL5D/DNA and GL5D-RGD/DNA respectively. After 6 h of incubation, lipoplexes were removed and 200 L of 65 media with 10% FBS was added. After 42 h, 20 L of MTT solution was added and the cells were incubated further for 4 h. Blue formazan crystals were seen in the wells under a microscope. The entire media were removed and 200 L of DMSO was added per well. The absorbance was measured using 70 microtiter plate reader. The % viability was then calculated as = [{A 590 (treated cells)-background]/ [A 590 (untreated cells)-background}] x 100.

Figure 1 .

 1 Figure 1. Molecular structures of the gemini lipid (A) GL5 and (B) a Palmitoyl-RAFT-RGD 4 (1) ligand used in the present investigation.

Scheme 1 .

 1 Scheme 1. Synthesis of 1: (i) 1 equiv. Pd(PPh 3 ) 4 , 100 equiv. PhSiH 3 ; (ii) 2 equiv. PyAOP, DIPEA (pH 8.0); TFA, TIS, H 2 O 45

Figure 2 .

 2 Figure 2. Transfection efficiency of GL5D gemini lipid with and without incorporation of Palmitoyl-RAFT-RGD 4 (100:1) and 55

Figure 3 .

 3 Figure 3. Transfection efficiency of GL5D gemini lipid with and without RGD-RAFT (100:1) in presence of serum (-FBS+FBS) at N/P (lipid/DNA) charge ratio of 0.5. Plasmid DNA (PGL-3) [0.8 g/well] was used for the experiment. Data are expressed as relative light unit/mg of protein as obtained from single luciferase 100

4

  |J. Mater. Chem. B, [2014], [vol], 00-00 This journal is © The Royal Society of Chemistry [year]

Figure 4 .

 4 Figure 4. DNA binding efficiency of GLD and GLD-RGD liposomes. Experiment was performed using 0.2μg pEGFP-C3 5

Figure 5 . 10 ratio 3 .

 5103 Figure 5. SDS mediated release of the DNA from GLD and GLD-RGD based lipoplexes. Experiment was performed using 0.2 μg plasmid DNA per well. Lipoplexes were prepared at N/P 10

Figure 6 .

 6 Figure 6. DNA protection ability of various liposomes against 20

Figure 7 .

 7 Figure 7. Surface binding efficiency of GLD and GLD-RGD formulations in HeLa and NIH3T3 cells. Bar diagrams shows high surface binding of formulations with RAFT-RGD in HeLa cells compared to NIH3T3 cells.

Figure 8 .

 8 Figure 8. Confocal images showing GFP expression of pEGFP-C3 gene transfected to HeLa (A-F) and NIH3T3 (G-L) cells. Liposomal formulations were used at N/P charge ratio 0.5 using 0.8 µg DNA per well. (A) HeLa; (B) + Effectene; (C) + GL5D; (D) + GL5D-RGD; (E) + only pEGFP-C3; (F) + Lipo2000*; (G) NIH3T3; (H) + Effectene; (I) + GL5D; (J) + GL5D-RGD; (K) + only pEGFP-C3; (L) + Lipo2000*.

40

 40

Figure 9 .

 9 Figure 9. Cytotoxicity assay (MTT) performed on HeLa (A and B), NIH3T3 (C and D) AND HEK293 (E and F) cells. Experiment was performed by using 0.2 µg DNA per well in 96 well plates. Experiment was performed in triplicate.Discussion45

  spectrometers. Infra-red (IR) spectra were recorded on a Jasco FT-IR 410 spectrometer using KBr pellets or as neat.Trans-sonic T 460H Elma bath sonicator (40W) was used for sonic dispersal of lipid suspension. Blood serum used in cell culture experiments was obtained from GIBCO. Cell culture flasks and multi-well plates were obtained from NUNC. Dulbecco's modified phosphate buffer saline (DPBS) and Dulbecco's modified eagles medium (DMEM) was obtained from 50 SIGMA. Public domain software, WinMDI was used for the FACS data analysis.

110 8

 8 |J. Mater. Chem. B, [2014], [vol], 00-00 This journal is © The Royal Society of Chemistry [year]

Journal of Materials Chemistry BJournal of Materials Chemistry B Accepted Manuscript

Journal of Materials Chemistry B

Acknowledegment. We thank DST grant (J C Fellowship to SB)


Hoang Trieu

Vy Le

Nelly Pustelnik

Marion Foare

The faster proximal algorithm, the better unfolded deep learning architecture ? The study case of image denoising

Keywords: Image denoising, unfolded proximal algorithms, accelerated methods, deep learning

Deep learning has revolutionized many image processing tasks such as classification or segmentation and, more recently, gives very promising results for solving inverse problems. It however remains a gap between the deep learning blackbox approaches and the more recently unrolled deep learning techniques proposed to bring the physics of the model and standard solving techniques into the network design. In order to understand more precisely the mechanisms, we place ourselves in the framework of the simple study of image denoising and we study four networks designed from unrolled forward-backward iterations in the dual, FISTA in the dual, Chambolle-Pock, and Chambolle-Pock exploiting the strong convexity. Performance and stability obtained with each of these networks will be detailed. A comparison within these approaches, standard penalized likelihood approaches, and the state-of-the-art black-box approach DnCNN is also provided.

I. INTRODUCTION

Deep learning has renewed many fields of image processing tasks such as classification, segmentation and, more recently, inverse problems solving. For several years, there was a gap between standard image processing and neural network procedures as the first one was guided by the physics of the data acquisition and prior knowledge about the object to analyze or recover while the second one was considered as a very efficient prior-free black-box procedure.

Recent works make the bridge between standard image processing techniques and deep neural networks architectures helping in the understanding and the analysis of this complex highly nonlinear tool.

On the one hand, several works have established relations between deep strategies and wavelet transform literature, thus inheriting wavelet properties. Wavelet scattering transforms that cascade wavelet transform with nonlinear modulus and averaging operators, benefit from a mathematical framework explaining important properties of deep convolution networks This work is supported by the ANR (Agence Nationale de la Recherche) from France ANR-19-CE48-0009 Multisc'In and Labex MI-LYON / ANR-10-LABX-0070. We are also grateful for the supports by CBP (Centre Blaise Pascal-ENS lyon) who provide necessary equipments for our numerical experiments.

[1], [START_REF] Oyallon | Generic deep networks with wavelet scatteringe oyallon, s mallat, l sifreinternational conference on learning representations[END_REF] while in [START_REF] Bolcskei | Optimal approximation with sparsely connected deep neural networks[END_REF] the authors prove that affine systems, including several wavelets or frame transforms, can be approximated by deep neural networks with minimal connectivity and memory requirements.

We can also refer to the link between neural networks and nonlocal denoising filtering techniques (e.g., nonlocal means [START_REF] Milanfar | A tour of modern image filtering: New insights and methods, both practical and theoretical[END_REF] or BM3D [START_REF] Dabov | Image denoising by sparse 3D transform-domain collaborative filtering[END_REF]) through their neural tangent kernel [START_REF] Tachella | The neural tangent link between CNN denoisers and non-local filters[END_REF], allowing to avoid the costly training step.

On the other hand, the deep unfolded approaches (cf. LISTA for the pioneering work [START_REF] Gregor | Learning fast approximations of sparse coding[END_REF]) have helped connect penalized likelihood approaches to neural networks. The two strategies feed off each other providing then a better understanding of neural network architecture and a formal framework toward the design of robust neural networks based on explicit links between activation functions and proximity operators [START_REF] Combettes | Lipschitz certificates for layered network structures driven by averaged activation operators[END_REF]. Focus of unfolded architectures -The unfolded deep architectures rely both on an objective function to minimize (e.g., F ) and an algorithmic procedure adapted to the properties of the objective function. For almost any given image processing problem, there exists a massive literature dedicated to the choice of the objective function F and to algorithmic schemes to design a sequence of iterates (x k ) k∈N that converges efficiently to a minimizer x of F . A standard criterion to measure the efficiency of an algorithm is the convergence rate. For instance, handling with strong convexity property into Chambolle-Pock iteration leads to

F (x k )-F ( x) ≤ ζ/k 2 while the standard iterations are reduced to F (x k ) -F ( x) ≤ ζ/k
where ζ is a constant [START_REF] Chambolle | A first-order primal-dual algorithm for convex problems with applications to imaging[END_REF].

However, to the best of our knowledge, when considering unfolded architectures, the impact of accelerated schemes on learning performance has not been studied yet and is the core of our contribution. To better understand the mechanisms with or without employing accelerated schemes and also to facilitate the presentation of the proposed unfolded networks, we handle this question in the context of image denoising focusing on two standard algorithms of the image processing literature: Iterative Soft Thresholding Algorithm (ISTA) [START_REF] Figueiredo | An EM algorithm for waveletbased image restoration[END_REF]- [START_REF] Combettes | Signal recovery by proximal forwardbackward splitting[END_REF] and Chambolle-Pock (CP) [START_REF] Chambolle | A first-order primal-dual algorithm for convex problems with applications to imaging[END_REF]. Contributions -The contributions of this work are: (i) the design of end-to-end trainable unfolded deep architectures for image denoising based on two standard iterative algorithms (ISTA and CP) and their associated accelerated version (FISTA [START_REF] Beck | A fast iterative shrinkage-thresholding algorithm for linear inverse problems[END_REF], [START_REF] Chambolle | On the Convergence of the Iterates of the Fast Iterative Shrinkage/Thresholding Algorithm[END_REF] and CP with strong convexity assumption [9, Algorithm 2]), (ii) the formulation of Deep-(F)ISTA-GD and Deep-(Sc)CP-GD as two neural networks expressed as the combination of weight and activation functions with closed-form expressions. Each of them includes an extra parameter allowing to consider the basic algorithm or its accelerated version, (iii) a study of the performance and robustness of the proposed networks and comparison to standard DnCNN regarding the learning cost and the stability.

The remainder of this paper is organized as follows. Section II provides a brief recall of the considered iterative algorithms and introduces the unfolding procedure of these algorithms. The robustness of the proposed networks is also discussed. Section III is dedicated to the numerical experiments and comparisons with state-of-the-art methods.

II. UNFOLDED DEEP LEARNING ARCHITECTURE

A. Image denoising and associated objective function

Our study focusses on the standard image denoising task, which consists to recover the original image x ∈ R N from its degraded version z = x + n ∈ R N , where n ∼ N (0, δ 2 I N ) denotes an additive white Gaussian noise with a standard deviation δ > 0. Thus, considering a prior imposing sparsity after some linear transform, and according to Bayesian interpretation, maximizing the a posteriori distribution reads:

x = argmin x∈R N F (x) := 1 2 x -z 2 2 + g(Dx), (1)
where D ∈ R |F|×N denotes a linear operator that converts the image in R N to a feature space R |F| with |F| > N , and g denotes a proper, lower-semicontinous, convex function from R |F| to (-∞, +∞]. A standard choice for g is a 1 -norm or a hybrid 1,2 -norm to favor the sparsity of the features. In the standard penalized likelihood approaches, a regularization parameter appears in front of g and the variance of the noise can be added in front of the data fidelity-term. In this work, these parameters are merged within D.

B. Proximal algorithms: (F)ISTA and (Sc)CP

Despite the simplicity of the minimization problem (1), its resolution requires to handle with iterative schemes as no closed form expression exists in a general framework (see a contrario when D is orthonormal, e.g., a wavelet transform). In this work, we focus on two algorithmic schemes: ISTA and CP, as both offer accelerated procedures: FISTA or CP involving strongly convexity (ScCP). (F)ISTA in the dual -(F)ISTA requires to handle the dual formulation of problem (1) to have closed form expression steps, leading to x = z -D u where 1 :

u ∈ Argmin u∈R |F| F (u) := 1 2 D u -z 2 2 + g * (u) (2)
1 The constant terms depending solely on z as been removed.

with g * the Fenchel conjugate of g (e.g. g = • 1 then g * = ι • ∞≤1 the indicator function of the ∞ -ball). The resulting FISTA iterations [START_REF] Beck | A fast iterative shrinkage-thresholding algorithm for linear inverse problems[END_REF] [14] read, for every iteration k,

    u k+1 = prox τ k g * (Id -τ DD )y k + τ Dz y k+1 = (1 + α k )u k+1 -α k u k (3)
where u 1 ∈ R |F| , and

y 1 ∈ R |F| . The sequence (u k ) k∈N converges to u when α k = t k -1 t k+1 and t k+1 = k+a-1 a , a > 2, τ < 1 D 2 and F (u k ) -F ( u) ≤ ζ k 2 .
When α k ≡ 0, these iterations reduce to ISTA. The convergence of the iterates is proved when τ < 2 D 2 for this limit case, and

F (u k ) -F ( u) ≤ ζ k .
(Sc)CP -CP iterations can be directly applied to the minimization problem [START_REF] Bruna | Invariant scattering convolution networks[END_REF]. The data-term being strongly convex of parameter γ = 1, the accelerated CP (ScCP: for Strongly convex CP) [9, Algorithm 2] can be employed, leading to, for every iteration k,

     u k+1 = prox τ k g * u k + τ k D (1 + α k ) x k -α k x k-1 x k+1 = σ k 1+σ k z + 1 1+σ k x k -σ k 1+σ k D u k+1 (4)
where

α k = 1 √ 1+2γσ k , σ k+1 = α k σ k , τ k+1 = τ k α k . The update of x k+1 comes from the proximity operator expression of the 2 -norm square: prox σ k 2 •-z 2 2 x k -σ k D u k+1 = σ k z+x k -σ k D u k+1 1+σ k
. This general framework fits both the ScCP iterations [9, Algorithm 2] and the standard CP iterations [9, Algorithm 1] when γ = 0, σ k ≡ σ, τ k ≡ τ and assuming στ D 2 < 1. The sequence (x k ) k∈N converges to the minimizer of (1) and convergence rates have already been provided in the introduction.

C. Unfolded architecture for (F)ISTA and (Sc)CP

From a training set S = {(x s , z s )|s = 1, . . . , I} where x s denotes a clean image, and z s the associated noisy one, the goal of a deep denoiser f Θ is to learn parameters Θ, to minimize the following standard empirical loss:

Θ ∈ Argmin Θ E(Θ) := 1 I I s=1 x s -f Θ (z s ) 2 2 . (5)
Standard networks with K layers can be written as, ∀z ∈ R N ,

f Θ (z) = η [K] • • • η [1] W [1] z + b [1] • • • + b [K] , (6)
where, for every k ∈ {1, . . . , K}, in the generic learning framework (e.g DnCNN [START_REF] Zhang | Beyond a Gaussian denoiser: Residual learning of deep CNN for image denoising[END_REF]), W [k] denotes a linear transform such as convolutions or pooling to reduce the number of parameters, associated with a bias b [k] , and a nonlinear activation function η [k] (e.g ReLu, sigmoid, HardTanh . . . ). We propose two networks taking the form of ( 6): Deep-(F)ISTA-GD (Network 1) and Deep-(Sc)CP-GD (Network 2). Illustrations of the two networks are provided in Figure 1.

                               W [1] = D [1] 1 (Id |F| -D [1] 1 D [1] 2 )D [1] 1 , b [1] = 0 D [1] 1 z l , η [1] = Id |F| HardTanh λ , W [k] = 0 Id |F| -α k-1 (Id |F| -D [k] 1 D [k] 2 ) (1 + α k-1 )(Id |F| -D [k] 1 D [k] 2 ) , b [k] = 0 D [k] 1 z l , η [k] = Id |F| HardTanh λ , W [K] = 0 -D [k] 2 , b [K] = z l , η [K] = IdN .
Network 2. Deep-(Sc)CP-GD for Gaussian Denoising has the following architecture, for every k ∈ {2, . . . , K -1}:

                       W [1] = IdN 2D [1] 1 , b [1] = 0 0 , η [1] = IdN HardTanh λ , W [k] = 1 1+σ k-1 - σ k-1 1+σ k-1 D [k-1] 2 1+α k 1+σ k-1 D [k] 1 -α k D [k] 1 Id |F| - (1+α k )σ k-1 1+σ k-1 D [k] 1 D [k-1] 2 , b [k] = σ k-1 1+σ k-1 z (1+α k )σ k-1 1+σ k-1 D [k] 1 z , η [k] = IdN HardTanh λ , W [K] = IdN 0 , b [K] = 0, η [K] = IdN .
P λ • * 1 (x) = HardTanh λ (x) = (p i ) 1≤i≤N where p i =      -λ if p i < -λ, λ if p i > λ, p i otherwise.
Proposition 2. We set, for every k ∈ {1, . . . , K}, D

[k] 1 ∈ R |F|×N , D k 2 ∈ R N ×|F| , W [k] , b [k] and η [k] provided by Net- work 1. If g = • 1 , D [k] 1 = τ k D and D [k] 2 = D , u 0 = u 1 = D [1] 1 z l , y 1 = (Id -D [1] 1 D [1] 2 )D [1]
1 z l , then Deep-(F)ISTA-GD network fits the generic (F)ISTA scheme (3).

A similar proposition holds for Deep-(Sc)CP-GD. Proposition 3. We set, for every k ∈ {1, . . . , K},

D [k] 1 ∈ R |F|×N , D [k] 2 ∈ R N ×|F| , W [k] , b [k] , and η [k] provided by Network 2. If g = • 1 , D [k] 1 = τ k D and D [k] 2 = D , u 1 = z l , α 1 = 1, y 1 = D [1]
1 z l , then the Deep-(Sc)CP-GD network fits the generic (Sc)CP scheme (4).

D. Robustness of the newtork

Following [START_REF] Combettes | Lipschitz certificates for layered network structures driven by averaged activation operators[END_REF], Network 1 and Network 2 have a Lipschitz behavior with Lipschitz constant χ = K k=1 W [k] . Thus, given an input z and a perturbation , we can majorize the perturbation on the output via the inequality

f Θ (z + ) -f Θ (z) ≤ χ .
χ can be used as a certificate of the robustness of the network provided that it is tightly estimated. Tighter bound exists but at the price of much more complex computations.

III. EXPERIMENTS

This section provides numerical comparisons between Deep-ISTA-GD, Deep-FISTA-GD, Deep-CP-GD and Deep-ScCP-GD, and aims to illustrate the impact of handling with unfolded deep accelerated schemes. Comparisons with state-of-the-art methods are also provided.

A. Experimental setting

Training and testing datasets -We combined RGB-images obtained from BSD500 and BSD300 [START_REF] Martin | A database of human segmented natural images and its application to evaluating segmentation algorithms and measuring ecological statistics[END_REF] to construct our dataset including 300 training images of size 180 × 180 and 300 testing images of size 320 × 320 both being cropped from the database, as suggested in [START_REF] Zhang | Beyond a Gaussian denoiser: Residual learning of deep CNN for image denoising[END_REF]. In our experiments, no data augmentation like flip or rotation has been done. The images are degraded with a white Gaussian noise with standard deviation δ = {50, 100}. Strategy to compare -We compare four configurations of the unfolded strategies described in Section II:

Deep-ISTA-GD with parameter Θ = {D [k] 1 , D [k] 2 , α [k] ≡ 0}, Deep-FISTA-GD with parameter Θ = {D [k] 1 , D [k] 2 , α [k] }, Deep-CP-GD with parameter Θ = {D [k] 1 , D [k]
2 , σ [k] , α [k] ≡ 1} suited to the case without handling with strong convexity (i.e., γ = 0), Deep-ScCP-GD with parameter Θ = {D

[k] 1 , D [k] 2 , σ [k] , α [k] = 1 √ 1+2σ [k]
} associated with the choice γ = 1. Parameter setting -We train all of our networks on Pytorch with ADAM optimizer [START_REF] Kingma | Adam: A method for stochastic optimization[END_REF] with 500 epochs and a learning rate set to 10 -4 . The batch size is set to 4 in Section III B and set to 10 for Sections III C-D. The impact of the network depth K and the feature size |F| will be evaluated.

B. Standard versus accelerated unfolded schemes

In Figure 2 2[top] displays the performance in terms of average PSNR over the testing dataset w.r.t (K, |F|). The map reported in Figure 2[bottom] displays the value of χ measuring the robustness for each network, where for every k, W [k] is computed with the power method.

We observe on the top row of Fig. 2 that Deep-FISTA-GD achieves better PSNR that Deep-ISTA-GD and that Deep-ScCP-GD achieves better PSNR that Deep-CP-GD for every set of parameters (|F|, K), allowing to conclude that accelerated schemes perform better. Additionally, Deep-ScCP-GD appears to be the most efficient method among the four for each range of parameters. In Fig. 2[bottom], we can observe that all networks have similar behaviours. The deeper is the network, the larger is the parameter χ. This tends to conclude that the deeper is the network the less robust is the network, which might be explained by the huge amount of parameters to learn for deeper networks.

C. Learning behaviour and comparison with DnCNN

In this section, the proposed networks are trained with K = 13 and |F| = 21 being a good compromise between performance and robustness. We compare the proposed networks to DnCNN that appears to be one of the most efficient according to [START_REF] Zhang | Beyond a Gaussian denoiser: Residual learning of deep CNN for image denoising[END_REF]. DnCNN has been re-implemented and retrained with K = 9 and |F| = 13. K and |F| has been selected to obtain a similar number of parameters to train for all networks (i.e. ∼ 14800 parameters). The batch size is set to 10, leading to a gain of almost 1 dB for all methods compared to a training with batch size of size 4.

Training loss

PSNR on test dataset In Fig. 3, we observe that the training losses converge faster for accelerated-based unfolded schemes and that Deep-ScCP-GD reaches a higher PSNR faster than DnCNN. In terms of performance, from the PSNR plot, when δ = 50, only Deep-ScCP-GD provides better (even close) results than DnCNN and we can observe that Deep-ScCP-GD is much more stable than DnCNN (cf. oscillations in the black plot). For a larger value of δ, the benefit of Deep-ScCP-GD compared to DnCNN is observed both in terms of PSNR and stability. In the next section, we only run experiments with Deep-ScCP-GD.

D. Visual comparisons

In Figure 4,we provide visual comparisons between our proposed Deep-ScCP-GD and state-of-the-art methods: TV [START_REF] Rudin | Nonlinear total variation based noise removal algorithms[END_REF], NL-TV [START_REF] Chierchia | A nonlocal structure tensor-based approach for multicomponent image recovery problems[END_REF], DnCNN [START_REF] Zhang | Beyond a Gaussian denoiser: Residual learning of deep CNN for image denoising[END_REF]. The hyperparameters in TV and NL-TV are trained on the training dataset with a multigrid search-parameter strategy described in [START_REF] Le | A proximal based strategy for solving Discrete Mumford-Shah and Ambrosio-Tortorelli models[END_REF]. TV and NL-TV fits the minimization of ( 1

Fig. 1 .

 1 Fig. 1. Architecture of the proposed Deep-(F)ISTA-GD (top) and Deep-(Sc)CP-GD (bottom).

Network 1 .

 1 Deep-(F)ISTA-GD for Gaussian Denoising has the following architecture, for every k ∈ {2, . . . , K -1}:

Proposition 2 (

 2 resp. Proposition 3) establishes the relation between Network 1 (resp. Network 2) and (F)ISTA (resp. (Sc)CP algorithm) described in section II-B. These relations are established when g = • 1 noticing that the proximity operator of the conjugate of the 1 -norm fits the HardTanh activation, a standard activation function used in deep denoiser. Proposition 1. The proximity operator of the conjugate of the 1 -norm scaled by parameter λ > 0 fits the HardTanh activation function, i.e., for every x = (x i ) 1≤i≤N :

  , we provide a comparison between Deep-ISTA-GD, Deep-FISTA-GD, Deep-CP-GD, and Deep-ScCP-GD for different choices of layer number K = {5, 13, 21, 29, 37} and feature size |F| = {13, 21, 29, 37, 45}. The maps reported in Figure

Fig. 2 .

 2 Fig. 2. 1st row: PSNR and 2nd row: χ (exponential scale) on a grid comparison between different choices of depth K and number of features |F| for each model.

Fig. 3 .

 3 Fig. 3. Comparisons between Deep-ISTA-GD , Deep-FISTA-GD , Deep-CP-GD , Deep-ScCP-GD, and DnCNN in terms of Training loss and Averaged PSNR on the testing dataset. (top) δ = 50 (bottom) δ = 100.

Fig. 4 .

 4 Fig. 4. Comparison between TV[START_REF] Rudin | Nonlinear total variation based noise removal algorithms[END_REF], NL-TV[START_REF] Milanfar | A tour of modern image filtering: New insights and methods, both practical and theoretical[END_REF], DnCNN[START_REF] Zhang | Beyond a Gaussian denoiser: Residual learning of deep CNN for image denoising[END_REF] and our proposed network Deep-ScCP-GD on different RGB images degraded with a white Gaussian noise of standard deviation δ.

IV. CONCLUSION

In this work, we proposed two deep unfolded networks for gaussian denoising that can be activated from standard or accelerated schemes. All networks are related to standard proximal algorithms: ISTA, FISTA, CP, and ScCP. Our numerical experiments illustrate the benefit to unroll accelerated schemes when possible. Additionally, the proposed scheme Deep-ScCP-GD is the most efficient from the PSNR performance and stability, and obtain better results than state-of-theart methods. Further experiments on the impact of accelerated schemes for solving more complex inverse problems would be investigated in future works.


Olivier Kimmoun

Hubert Branger

SPACE-TIME MEASUREMENTS OF BREAKING WAVE KINEMATICS AND VOID FRACTION IN THE SURF ZONE

Particle Image Velocimetry (PIV) measurements were performed in a wave tank under regular water waves propagating and breaking on a 1/15 sloping beach. To obtain more accurate measurements, the domain from the initiation of the breaking to the swash zone is splited in fourteen overlapping windows. The full field is reconstructed at each time step by gathering the fourteen synchronized PIV fields. We measured then the complete space-time evolution of velocity field over the whole surf zone. From these data the ensemble-period-average and phase-average components of the flow with their associated fluctuating parts were computed. Moreover using the light intensity of the PIV images and velocity measurements, we estimated the void fraction in each point of the surf zone and at each time step.

INTRODUCTION : EXPERIMENTAL SETUP AND METHOD

Experiments were conducted in the Ecole Centrale wave tank in Marseilles. The glass-windowed tank is 17 m long and 0.65 m width. The water depth in front of the wavemaker is about 0.705 m. A polyvinyl 1/15 sloping beach was mounted on the bottom. The beach, which is about 13 m long, is partly made with transparent glass in order to let the laser light sheet cross vertically the beach from the bottom to the free surface. To characterize the wave field, water elevation was measured with a set of six resistive wave gauges. Velocity field in the water was measured with a PIV technique. In the foam generated during the breaking event, seeds were not clearly visible. Then, the PIV intercorrelation algorithm was not really based on particle displacements, but much more on micro-foam structure and bubble displacements, as explained in [START_REF] Govender | Video-imaged surf zone and roller structures and flow fields[END_REF]. In this experiment, the length of the surf zone was about 3 m long. In order to retrieve the full velocity field all along this distance, we made camera PIV measurements in fourteen different locations ranging from the incipient breaking location up to the swash zone. The size of the PIV images was 37x37 cm 2 . We synchronized PIV image acquisition and wave gauge data using the same time trigs and duration acquisitions for all experiments. This synchronization allowed us to reconstruct, at each time step, a full velocity field covering the entire surf zone. Such experimental measurements and analysis are very time consuming and only one paddle condition for regular waves was conducted. The wave period was T = 1.275 s and the surf similarity parameter was ζ = 0.28. The breaking process starts with a plunging phase, followed by successive splash-up, then the wavefront moves into a roller propagating towards the shore. The PIV acquisition rate was 10 Hz, so with the chosen wave period, the duration of 51 successives PIV velocity field corresponds exactly to 4 wave periods. Consequently all the PIV velocity fields separated by 51 PIV time steps were exactly in phase. By averaging the instantaneous fields corresponding to each phase, the phase-average field and its corresponding fluctuating part were obtained. The different terms that appear in the fluctuating kinetic energy propagation equation can be also estimated [START_REF] Kimmoun | A particle image velocimetry investigation of laboratory surf-zone breaking waves over a sloping beach[END_REF]. We estimated at each time step, the void fraction in each point of the domain by means of light intensity image analysis and vertically-integrated velocity considerations in the domain. We used light scattering theory with a Lambert law to relate void fraction with light intensity.

WAVE KINEMATICS

Phase average vorticity

In order to have a statistical evolution of the periodic flow during a wave cycle, we computed a phase-average of velocities and water profiles by averaging all the instantaneous data that were acquired at the same phase. This procedure is commonly used to study periodic events. The vorticity generated by wave breaking, under the trough-envelope, is of the same order of magnitude as phase speed divided by local depth D. Dimensionless phase-averaged vorticity field Ω/( √ gD/D) is displayed on figure 1-Left during the spash-up event. A large clockwise vortex structure appears first on the front side of the wave crest at the beginning of the plunging phase. This vortex is not localized at the top of the crest, but rather at the toe of the front, just above the mean sea level. The vorticity magnitude increases during the plunging. The first splash-up sequence creates a new vortex ahead of the first one (figure 1-left). Then successive vortices are generated during the following splash-up sequences. The location and intensity of the vortices generated during a wave cycle are presented in figure 2. The first vortex moves on shore with increasing magnitude during the plunging phase, then propagates more slowly downward with a decreasing intensity. The second vortex, initiated by the first splash-up, reaches 30s -1 , moves shoreward and then downward with a decreasing magnitude. The third and fourth vortices moves towards the shoreline along the trough-envelope line. Their magnitude is of order 20s -1 when they propagate near the bore-front, then they slow-down and vanish slowly when they go deeper, dissipated by bottom friction and slowed down by mean return flow.

Fluctuating kinetic energy

An example of the Froude scaled phase-averaged turbulent kinetic energy:

q 2 1/2 / √ gD is presented in figure 1-right.
During the breaking event, we fond that the initiation point for fluctuation generation was located at the top of the front side of the breaking crest. Values of q 2 1/2 increase during the plunging event, and the splash-up phase (figure 1-right).

Then, generated fluctuations are spread downstream and towards the bottom. At this phase, we found again three cells of fluctuating intensities. These cells are located on the edge of the vorticity cells. During the rolling phase, after the splash-up, the non-dimensional velocity scale attached to the bore does not seem to decrease until the bore reaches the shoreline. The fluctuating intensity is maintained by the moving breaking process.

ESTIMATION OF THE SPACE-TIME EVOLUTION OF THE VOID FRACTION

We estimated the void fraction in each point of the domain by means of light intensity image analysis and verticallyintegrated velocity considerations in the domain [START_REF] Kimmoun | A particle image velocimetry investigation of laboratory surf-zone breaking waves over a sloping beach[END_REF]. We used light scattering theory with a Lambert law to relate void fraction with light intensity. An example is presented in figure 3-left. We found an exponential decay of void fraction with depth below the mean sea level, with a very good fit to the Hoque and Aoki (2005) semi-theoretical law [START_REF] Hoque | Distributions of void fraction under breaking waves in the surf zone[END_REF](figure 3-right). It is the first time that an space-time evolution void fraction over all the laboratory surf zone is measured. The airentrainment is large during the first splash-up event with values of void fraction near 0.88 above mean sea level inside the splash-up. Then the void fraction decreases slowly when the wave propagates towards the shore, with values near 0.2-0.3 . Close to the shore, aeration is still present in the front face of the wave, but also at the back face of the bore near the crest.

We acknowledge the support of INTAS grant 06-1000013-9236 and French National program LEFE/IDAO.

Figure 1 .

 1 Figure 1. Splash-up phase : (Left) Non-dimensioned phase-averaged vorticity Ω/ g/D, (Right) Phase-locked Froude normalized fluctuating kinetic energy (q 2 /gD) 1/2

Figure 2 .

 2 Figure 2. Intensity (Left) and position (Right) of the four observed vortices generated during a wave cycle (1 st : red circle, 2 sd : blue square, 3 rd : magenta star, 4 th : green triangle).

Figure 3 .

 3 Figure 3. (Left) Phase-averaged void fraction during the splash up, (Right) Measured vertical distribution of mean void fraction α at different locations (solid red line) and comparison with the semi-theoretical law [3] (dashed blue line)