Azary Abboud

Student Member, IEEE, Franck Iutzeler, Member, IEEE Romain Couillet

Member, IEEE Mérouane Debbah

Fellow Member, IEEE Houria Siguerdidjane

Distributed Production-Sharing Optimization and Application to Power Grid Networks

Keywords: Convex optimization of large scale problems, ADMM, randomized methods, smart grids, flow calculations, distributed control

Based on recent works on asynchronous versions of the distributed Alternating Direction Method of Multipliers (ADMM) algorithm, we develop and prove the convergence of a distributed asynchronous method for Production-Sharing Problems over networks. The asynchronous nature of the algorithm allows both for the relaxation of the synchronization constraint often inherent to distributed ADMM-based methods and distributed optimization methods at large, but also allows for random local failures to occur in fully centralized methods. These two considerations motivate the application of the method to the Direct-Current Optimal Power Flow (DC-OPF) problem in power transmission networks. Applied to the DC-OPF, this method leads to an overall network optimal production obtained through a sequence of local computations in subareas of the network (each area waking up randomly while the rest of the network is non-operational) and neighboring data exchanges. In another scenario, the DC-OPF is performed via iterations of a centralized network-wide ADMM method which may contain disconnected nodes (in general with low probability and for a short duration). In both cases, this method still converges and thus provides additional flexibility to classical DC-OPF algorithms. The proposed algorithm, inherently designed for networks of overlapping subareas, is then extended to networks of nonoverlapping areas. Simulations are carried out on the IEEE-30 and IEEE-118 bus test systems which illustrate the convergence, scalability and effectiveness of the proposed algorithms.

I. INTRODUCTION

One of the salient features of power transmission systems is their ability to ensure an optimal transmission throughout the power grid of generated powers at a low economical cost. This naturally leads to the so-called Optimal Power Flow [START_REF] Glover | Power system analysis and design[END_REF] problem, widely studied up until very recently with the emergence of the smart-grid paradigm, see e.g., [START_REF] Giannakis | Monitoring and optimization for power grids: A signal processing perspective[END_REF]- [START_REF] Gan | Optimal power flow in direct current networks[END_REF] for recent works. This problem can be generally cast as that of finding an operational point leading to the least global generation cost subject to power flow equations (such as Kirchhoff and Ohm's laws) and other operational constraints [START_REF] Cain | History of optimal power flow and formulations[END_REF]. The physical constraints make the problem nonlinear and non-convex [START_REF] Wood | Power generation, operation, and control[END_REF], but as far as power transmission is concerned Azary Abboud and Mérouane Debbah are with Laneas Group, Supélec, Gifsur-Yvette, France. Franck Iutzeler is with the ICTEAM Institute, Université Catholique de Louvain, Louvain-la-Neuve, Belgium. Romain Couillet is with the Telecommunication Department, Supélec, Gif-sur-Yvette, France. Houria Siguerdidjane is with the Automatic and Control Department, Supélec, Gifsur-Yvette, France. Email: firstname.lastname@supelec.fr. Part of the material of this paper has been presented at the 39th IEEE International Conference on Acoustics, Speech, and Signal Processing (ICASSP), Florence, May 9, 2014 [START_REF] Abboud | Asynchronous alternating direction method of multipliers applied to the direct-current optimal power flow problem[END_REF]. linear approximations are considered accurate and significantly simplify the problem [START_REF] Gan | Optimal power flow in direct current networks[END_REF].

Our focus is here on the Direct-Current Optimal Power Flow (DC-OPF) problem, which is a linear DC approximation of the original OPF problem [START_REF] Glover | Power system analysis and design[END_REF]. In the DC-OPF formulation, voltage angles differences are considered close to zero and voltage magnitudes are fixed. 1 Embracing the smart-grid paradigm, we seek here to minimize the power generation cost of the grid while taking advantage of its advanced communication characteristics [START_REF] Glover | Power system analysis and design[END_REF], [START_REF] Stott | Review of load-flow calculation methods[END_REF]. A classical method to solve problems of the OPF class relies on a central fusion center which collects data from all the network agents and centrally evaluates the optimal system parameters before communicating these values to the corresponding agents [START_REF] Guerrero | Advanced control architectures for intelligent microgrids? Part II: Power quality, energy storage, and AC/DC microgrids[END_REF]. However, on the one hand, this centralized optimization method is rather prohibitive in large-scale and sparse networks due to the communication burden and the delay induced by the required computations and communication [START_REF] Guerrero | Advanced control architectures for intelligent microgrids? Part II: Power quality, energy storage, and AC/DC microgrids[END_REF], [START_REF] Amin | Toward a smart grid: power delivery for the 21st century[END_REF]. On the other hand, the privacy concerns of the end users or the grid owners in case of network with multiple authorities, and the recent tendency toward multinational electricity markets, make it impossible to efficiently implement the centralized optimization scheme. Besides, several security issues such as single point of failure, bottleneck and agreement on decisions, may arise in case of centralized optimization [START_REF] Guerrero | Advanced control architectures for intelligent microgrids? Part II: Power quality, energy storage, and AC/DC microgrids[END_REF], [START_REF] Kraning | Dynamic network energy management via proximal message passing[END_REF].

To reduce the computational burden and the inherent weakness of a single centralized computation unit, one may leverage on the presence of independent processors in each component of the power grid [START_REF] Amin | Toward a smart grid: power delivery for the 21st century[END_REF] to increase the available computational capacity. This implies having a distributed computing platform, which naturally calls for the intervention of distributed optimization approaches in order to solve the management problems in the power grid. This distributed scheme is specifically preferred when power sources are located at different parts of a large network while being under the control of multiple companies or authorities. We thus aim here to solve the DC-OPF problem by using several processors distributed all over the network. Distributed DC-OPF approaches where multi-utility transmission systems are divided into autonomously managed area were proposed in [START_REF] Kim | Coarse-grained distributed optimal power flow[END_REF], [START_REF] Conejo | A decomposition procedure based on approximate newton directions[END_REF] and in references therein. The distributed computation is a preferable approach that helps in reducing the data exchange overhead and the processing time, leading to a smart optimization of the power system. The pioneering work of Tsitsiklis and Bertsekas [START_REF] Tsitsiklis | Problems in decentralized decision making and computation[END_REF], [START_REF] Tsitsiklis | Distributed asynchronous deterministic and stochastic gradient optimization algorithms[END_REF] in this field has led to extensive studies on distributed optimization techniques and several methods were developed to solve various optimization problems in a distributed manner [START_REF] Duchi | Dual averaging for distributed optimization: convergence analysis and network scaling[END_REF], [START_REF] Nedic | Distributed subgradient methods for multiagent optimization[END_REF].

Several methods can be used to perform distributed optimization, either in a sequential or in a parallel manner, such as Lagrangian relaxation [START_REF] Zhuang | Towards a more rigorous and practical unit commitment by lagrangian relaxation[END_REF], augmented Lagrangian [START_REF] Kim | Coarse-grained distributed optimal power flow[END_REF], approximate Newton directions in conjunction with standard Lagrangian [START_REF] Conejo | A decomposition procedure based on approximate newton directions[END_REF], [START_REF] Luenberger | Linear and nonlinear programming[END_REF] as well as the auxiliary problem principle [START_REF] Losi | On the application of the auxiliary problem principle[END_REF]. These methods are however computationally cumbersome [START_REF] Kalyanmoy | Optimization for engineering design: Algorithms and examples[END_REF], [START_REF] Paudyal | Three-phase distribution opf in smart grids: Optimality versus computational burden[END_REF] as they share the same synchronization that ties the different areas and demand strong coordination (each area optimizes its variables in every iteration and communicates with its adjacent areas before moving to another iteration). The synchronous burden tends to cause degradation in performance mainly because of the inherent communication overhead and the idle processing time [START_REF] Kung | Synchronized and asynchronous parallel algorithms for multiprocessors[END_REF]. Among those methods, the technique known as the Alternating Direction Method of Multipliers (ADMM), while sharing the same performance limitations, has the advantage of converging faster than the aforementioned methods [START_REF] Eckstein | On the Douglas-Rachford splitting method and the proximal point algorithm for maximal monotone operators[END_REF], [START_REF] Boyd | Distributed Optimization and Statistical Learning via the Alternating Direction Method of Multipliers[END_REF], of being easier to implement, and of requiring less time and computational capacity at each processor. The structure of the iterations in this method makes their distribution among the existing processors easier than the previously mentioned methods [START_REF] Boyd | Distributed Optimization and Statistical Learning via the Alternating Direction Method of Multipliers[END_REF].

Recently [START_REF] Iutzeler | Asynchronous Distributed Optimization using a Randomized Alternating Direction Method of Multipliers[END_REF] proposed a fully asynchronous distributed version of the ADMM that is applied sequentially by randomly selecting overlapping areas of the network, the union of which covers the whole network. Despite the asynchrony, by relying on a careful analysis of the proximal splitting method [START_REF] Eckstein | On the Douglas-Rachford splitting method and the proximal point algorithm for maximal monotone operators[END_REF] of which the ADMM is a special instance, it is proved that such an algorithm converges. As opposed to [START_REF] Iutzeler | Asynchronous Distributed Optimization using a Randomized Alternating Direction Method of Multipliers[END_REF], which was not designed for DC-OPF problems specifically, in the present article, we additionally consider that the problem under study (being DC-OPF compliant) is such that the optimization variables are spread across the network (instead of being shared among the nodes) and that these variables have to satisfy some additional box constraints. In power grid terms, the power generated in the network are local quantities only bounded to an indirect user satisfaction constraint, and must be such that maximal capacities of power lines are not exceeded. In [START_REF] Iutzeler | Asynchronous Distributed Optimization using a Randomized Alternating Direction Method of Multipliers[END_REF] the distributed agents seek to find consensus on an (unconstrained) unique parameter minimizing a node-wise aggregate cost. As such, although closely related, our problem formulation differs from [START_REF] Iutzeler | Asynchronous Distributed Optimization using a Randomized Alternating Direction Method of Multipliers[END_REF]. We formulate the problem in a more generic form than the DC-OPF, the Production-Sharing Problem. Each agent has a production ability at a certain cost and can exchange its resources with its neighbors through the links connecting them. This generic formulation makes it applicable to problems involving distributed computations other than power systems (please refer to [START_REF] Abboud | Distributed caching in 5g networks: An alternating direction method of multipliers approach[END_REF] for an example). Because imposing overlapping areas in a power network structure may lead to practical management problems, inspired by the method of passing adjacent variables [START_REF] Conejo | Multi-area coordinated decentralized DC optimal power flow[END_REF] that allows to convert networks of non-overlapping and independent areas into a network of overlapping areas, we extend the approach from [START_REF] Iutzeler | Asynchronous Distributed Optimization using a Randomized Alternating Direction Method of Multipliers[END_REF] to cover this context.

The rest of the article is organized as follows. Section II defines the convex optimization problem with linear constraints for which we provide in Section III the distributed formulation and solution using ADMM. Then, in Section IV, relying on the distributed Production-Sharing optimization w/ ADMM algorithm, we develop the asynchronous distributed Production-Sharing w/ ADMM algorithm. In Section V, we extend our algorithm to the non-overlapping area scenario. Section VI is devoted to the application of the algorithms on the linear DC-OPF problem. The results are presented on the power grid network given by the IEEE-30 bus test system [32], then on the IEEE-118 bus test system [32] in order to verify the scalability of the proposed algorithms. Lastly, Section VII contains our concluding remarks.

Notations

Lowercase (resp. uppercase) boldface characters will denote vectors (resp. matrices). Furthermore, x i denotes the i th entry of vector x and A ij stands for the i th row j th column entry of matrix A.

[A] i corresponds to the i th row of the matrix A. x T and A T denote respectively the transpose of vector x and matrix A. In addition, we will note x = [x 1 ; . . . ; x N ] the vector formed by the vertical concatenation of vectors (x i ) i=1,...,N .

II. PROBLEM FORMULATION

In this section, we will formulate a generic sharing optimization problem which encompasses as a particular case the DC-OPF problem.

A. A Production-Sharing Problem on a Graph

Consider a group of N agents connected by an undirected connected graph of N nodes, G = (V, E), where V = {1, ..., N } is the set of nodes2 and E is the set of bidirectional edges. We will say that i and j are neighbors, or i ∼ j, if {i, j} ∈ E, and we will note N i = {j ∈ V : i ∼ j} the set of neighbors of i. We introduce symmetric edge weights w ij = w ji such that w ij > 0 if and only if {i, j} ∈ E and 0 elsewhere, the matrix W defined by the entries {w ij } is called the weighted adjacency matrix, the diagonal matrix D such that D ii = j∈V w ij is called the degree matrix, and finally L = D-W is called the Laplacian matrix of the graph.

We look at the case where each agent j ∈ V produces and/or consumes some resource (typically electric power) represented by a real value, and the agents can exchange this resource using the links of the aforementioned graph. More precisely, each agent has a production ability in the form of an interval (reduced to {0} for non-producers) and a cost function depending on the quantity produced (which will be assumed convex in the following). The natural problem ensuing is thus to minimize the total production cost under the constraints that i) the demands are satisfied and ii) the network is able to dispatch correctly the resource. This last point can take multiple forms depending on the resource type and will be made explicit when considering the DC-OPF problem next; however, we will write it generically as "the vector of the differences between production and consumption at each agent belongs to the span of the Laplacian matrix".

We write this problem formally as a convex optimization problem 3 with linear constraints:

Production-Sharing Problem min x,y,z∈R N j∈V f j (x j ) subject to ∀j ∈ V, r j ≤ x j ≤ r j y = x -d y = Lz (1)
where x is the resource production vector, d is the demand vector, f j is the production cost function at agent j and [r j , r j ] is its production range, L is the Laplacian of graph G, and finally y and z are intermediate vectors.

In order to be able to perform meaningful derivations, we will make the following standard assumption.

Assumption 1. The functions {f j } j∈V are convex, proper, lower semi-continuous and the graph G is undirected and connected. Furthermore, the set of minimizers of Problem (1) is non-empty.

A simple example: Consider a complete graph with unit edge weights and say that each of the N agents has a unit demand (d = [1, ..., 1] T ). Because the graph is complete L = N I -J where I is the identity matrix and J is filled with ones. This means that the second condition of the above problem simply translates to "the mean of the resource production vector x is one". Now, assume that half of the agents cannot produce any resource (for them, r j = r j = 0 thus x j = 0) while the other half can produce any quantity x (r j = 0, r j = +∞) with cost ax (for every producer j, f j (x j ) = ax j and a > 0); it is immediate to see that each producer will produce the same quantity (as they are indistinguishable in the problem with a complete unweighted graph) and this quantity must be 2 (as only half of the agents produce and the mean of the production vector has to be one).

Now, let us see how this problem can encompass problems encountered in power networks such as the DC-OPF.

B. DC-OPF as a Production-Sharing Problem on a Graph

A power network consists of a set of agents/buses V . Each agent j ∈ V can generate a power p g j with a cost f j (p g j ) in its generation ability range [r j , r j ] (possibly reduced to {0} if the agent is not a generator) as well as a power demand p d j . The agents are linked by transmission lines E, G = (V, E) forming a connected undirected graph. Let θ j be the phase of the current produced by node j ∈ V . In the DC model, the power flow on the transmission line from agent i to j is proportional to their phase difference θ i -θ j through a constant B ij representing the imaginary part of the element in the i th row and j th column of the bus admittance matrix Y [2], [START_REF] Gómez-Expósito | Electric energy systems: analysis and operation[END_REF], [START_REF] Powell | Power System Load Flow Analysis (Professional Engineering)[END_REF]. We also introduce the excess power of node j, p e j = p g j -p d j . Kirchoff's Voltage law and the DC model conditions impose the excess power (which is possibly negative) to be transmitted through the adjacent transmission lines. Formally, for every agent j, this sums up to p e j = i∼j B ji (θ j -θ i ) = ( i∼j B ji )θ j -i∼j B ji θ i . This condition can be rewritten in a vectorial manner as p e = L B θ where L B is the Laplacian matrix defined from the weighted adjacency matrix B as defined previously and θ = (θ 1 , . . . , θ N ) is the phase vector. Let p g = (p g 1 , . . . , p g N ) be the power generation vector and p d = (p d 1 , . . . , p d N ) be the power demand vector. Then, the DC-OPF Production-Sharing Problem can be formulated as follows,

DC-OPF Production-Sharing Problem min p g ,p e ,θ∈R N j∈V f j (p g j ) subject to ∀j ∈ V, r j ≤ p g j ≤ r j p e = p g -p d p e = L B θ (2)
It is straightforward to see that the above problem has the same form as Problem [START_REF] Abboud | Asynchronous alternating direction method of multipliers applied to the direct-current optimal power flow problem[END_REF]. It is important to notice that the matrix B (through its Laplacian) suffices to structure the generated power vector to comply with electrical equalities across the power network.

III. DISTRIBUTED PRODUCTION-SHARING

The Production-Sharing Problem introduced in the previous section is essentially centralized. Indeed, while the underlying graph is present through the Laplacian matrix, the condition y = Lz on the transmission of the excess production implies a coordinated action of the whole network. In this section, we first design an equivalent problem where the network condition is split into overlapping subgraphs and add an additional indicator function in order to ensure the equivalence with the former problem. Then, after noticing that this new problem is well suited for the Alternating Direction Method of Multipliers, we derive a distributed Production-Sharing algorithm solving our problem.

A. Distributed Formulation of the Problem

Starting from Problem (1), as every agent can control its resource generation (and thus its excess) as well as its network intermediate vector (its phase in the DC-OPF case), it is convenient to stack these variables into a vector of size 3N , u = [x; y; z].

With this new variable, it is useful to introduce a function F corresponding to the function to minimize in (1) applied to u and encompassing the range condition:

F : R 3N -→ (-∞, +∞] u =   x y z   -→    j∈V f j (x j ) if ∀j, x j ∈ [r j , r j ] +∞ elsewhere (3)
One can remark that, provided that the functions {f j } j∈V are convex, F is also convex. Now, the rest of the conditions can be written linearly as:

  I -I 0 0 -I L   A   x y z   u =   d 0   d (4)
Hence, Problem (1) rewrites as

min u∈R 3N F (u) subject to Au = d (5)
Let us recall that in our scheme, each agent j has the knowledge of its cost function f j , its production-related variables x j , y j , z j , and the j th row/column of the Laplacian matrix L. It is then straightforward to see that, whereas the upper part of the equality in ( 4) is local to each agent, the rows of the lower part corresponds to the network connections.

In order to derive a distributed algorithm, a common solution is to introduce one variable per condition line per agent and an indicator function ensuring that for each condition line, the sum of the agents related variables is equal to the corresponding value in d. Formally, define the 6N 2 × 3N matrix

à =    diag([A] 1 ) . . . diag([A] 2N )    (6)
and the indicator function4

G : R 6N 2 -→ (-∞, +∞] v -→      0 if ∀i = 1 . . . 2N, 3iN j=3(i-1)N +1 v j = di
+∞ elsewhere, (7) where, diag([A] i ) is the diagonal matrix whose diagonal entries are the elements of the row [A] i .

One can note that vector v has a lot more components than actually needed. Due to the sparsity of the matrix A, several components in v are equal to zero and can be omitted. But, for notation simplicity, the full vector will be kept until we provide the actual algorithm derivation, where we will eliminate unnecessary components and precise the actual number of variables to be updated per agent.

Finally, we obtain the following distributed problem.

Distributed Production-Sharing Problem

min u∈R 3N F (u) + G( Ãu) (8)
By construction, the solutions of this problem allow to derive the solutions of the original problem as mentioned in the following lemma.

B. Application of the ADMM to the Distributed Production-Sharing Problem

Writing an iteration of the ADMM algorithm on Problem (8) leads to the following set of equations

u k+1 = argmin u F (u) + ρ 2 Ãu -v k + λ k ρ 2 (9a) v k+1 = argmin v G(v) + ρ 2 Ãu k+1 -v + λ k ρ 2 (9b)
λ k+1 = λ k + ρ Ãu k+1 -v k+1 (9c)
where ρ > 0 is a free hyper-parameter and λ is the vector of Lagrangian multipliers.

Using the definitions of function F and G, and matrix Ã, we obtain the following algorithm after some derivations reported in Appendix A.

Distributed Production-Sharing optimization w/ ADMM 1) Initialize u and π to the initial values u 0 and π 0 .

2) At iteration k every agent j = 1, . . . , N performs the following steps : a) Update x j , y j and z j

x k+1 j =argmin xj f j (x j ) + ρ 2 x j + π k 1,j ρ -x k j + 1 d 1 (j) x k j -y k j -d j 2 (10a) y k+1 j =y k j + π k 1,j + π k 2,j 2ρ + r k 1,j 2d 1 (j) + r k 2,j 2d 2 (j) (10b) z k+1 j =z k j - 1 N i=1 L 2 ij N i=1 L ij π k 2,i ρ + r k 2,i d 2 (i) (10c) with ∀i = 1, ..., N , d 1 (i) = 2 and d 2 (i) = 2 + |N i | b) Communicate z k+1 j and y k+1 j to neighboring nodes. c) For each constraint i ∈ I j = {i; ∃A ip = 0, p = {j, j + N, j + 2N }}, compute r k+1 1,i = x k+1 i -y k+1 i -d i (11a) r k+1 2,i = -y k+1 i + N j=1 L ij z k+1 j (11b) π k+1 1,i = ρ d 1 (i) r k+1 1,i + π k 1,i (12a)
π k+1 2,i = ρ d 2 (i) r k+1 2,i + π k 2,i (12b)
3) If the stopping criterion is not satisfied, increase k and go to 2). Otherwise, retain x k+1 j , y k+1 j and z k+1 j .

This algorithm is effectively distributed; agent j computes and maintains its personal variables using local information (f j , L j,• ) and exchanges with its neighbors. We can also remark that, contrary to the consensus-based distributed ADMM [START_REF] Schizas | Consensus in ad hoc wsns with noisy linkspart i: Distributed estimation of deterministic signals[END_REF], the relations between the agents are directly linked to the original constraints and not due to a reformulation of an originally centralized problem. The proof of Theorem 1 is straightforward as Distributed Production-Sharing optimization w/ ADMM relies on ADMM applied to a sum of two functions that are proper lower semi-continuous and convex by construction (F is constructed from the (f i ) with Assumption 1; G is a projection onto a linear space). These properties are the sole conditions for the convergence of ADMM (see e.g. [START_REF] Bauschke | Convex analysis and monotone operator theory in Hilbert spaces[END_REF]) and thus for our algorithm.

The convergence speed is expected to be at least O(1/k) as the reasoning from [START_REF] Deng | On the o (1/k) convergence and parallelization of the alternating direction method of multipliers[END_REF] can be adapted; nevertheless, linear convergence speed have also been proved for other ADMMbased algorithms with additional assumptions on F [START_REF] Iutzeler | Explicit convergence rate of a distributed alternating direction method of multipliers[END_REF], [START_REF] Ghadimi | Optimal parameter selection for the alternating direction method of multipliers (admm): quadratic problems[END_REF].

IV. RANDOMIZED DISTRIBUTED PRODUCTION-SHARING

A. Theoretical foundations

In the previous section, we presented an algorithm for solving Problem (1) in a distributed manner over the graph associated with L. However, it is not always possible nor appropriate to compute a full iteration of this algorithm. For example, some agents of the network can randomly fail to exchange or compute their update. Also, it may be faster to solve the problem by taking into account only a random subset of the agents/links at each iteration in the spirit of mini-batch algorithms.

For these reasons, it is interesting to consider a randomized version of the previous algorithm where only some parts of the network are active at a given iteration.

This can be achieved by following the ADMM randomization scheme proposed in [START_REF] Iutzeler | Asynchronous Distributed Optimization using a Randomized Alternating Direction Method of Multipliers[END_REF]. At each iteration k, this method consists of picking a random set of coordinates ξ k (or equivalently a random subset of agents and links), and performing the updates of Eqs. (9b) and (9c) only for the coordinates of ξ k , the other ones being kept at their former value. As for the first update Eq. (9a), only the coordinates needed for the partial update of Eqs. (9b) and (9c) are to be computed.

In this way, provided that the random coordinate selection sequence (ξ k ) k>0 is independent and identically distributed (i.i.d.) and such that the selection probability is positive for every coordinate, the randomized algorithm converges almost surely to a sought solution [START_REF] Iutzeler | Asynchronous Distributed Optimization using a Randomized Alternating Direction Method of Multipliers[END_REF] (see also [START_REF] Bianchi | A stochastic coordinate descent primal-dual algorithm and applications to large-scale composite optimization[END_REF] for refinements).

B. Randomized algorithm

In the case of a network of agents with multiple authorities or when the network is divided into multiple micro-networks, it is advised to use area-based distributed mechanisms to manage efficiently the network. Interestingly, the components of λ and v are linked either to a node or to an edge of the network; thus, updating only a part of these components sums up to performing computations and exchanges only on a subgraph.

Let us decompose the graph G into L connex areas A , l ∈ {1, . . . , L} (see Fig. 1); area A is a graph with vertices V ⊂ V and edges E = {{i, j}; (i, j) ∈ V 2 } ∩ E. Each area A will act as a local processor exchanging data with its closest neighbors when selected by the random sequence (ξ k ) k>0 ; i.e. at times k such that ∈ ξ k+1 . In order to ensure convergence to the sought solution [START_REF] Iutzeler | Asynchronous Distributed Optimization using a Randomized Alternating Direction Method of Multipliers[END_REF], we will assume that i) the areas are connected and encompass the whole graph; and ii) (ξ k ) k>0 is i.i.d. and for all l ∈ {1, . . . , L}, the selection probability for area A is positive:

P[ ∈ ξ 1 ] > 0.
Assumption 2. For any l ∈ {1, ..., L}, let G(V ) (V , E ) be the sub-graph of area A , the following properties are assumed:

1) L l=1 V = V , 2) L l=1 G(V ) is connected.
Furthermore, the area selection sequence (ξ k ) k>0 , valued in the set of the subsets of {1, ..., L}, is independent and identically distributed and such that ∀ ,

P[ ∈ ξ 1 ] > 0.
This means that in order to ensure convergence, the areas must overlap (i.e., have at least an agent in common, see Fig. 1), which can be restrictive. Let us first state the asynchronous algorithm with overlapping areas and extend it to the non-overlapping case in next section.

At each iteration k, applying Eqs. (9a)-(9c) but keeping only the components related to the areas 5 of ξ k+1 in the two latter equations leads to the following algorithm (the calculations are quite similar as the ones in Appendix A).

Asynchronous Distributed Production-Sharing optimization w/ ADMM 1) Initialize u and π to the initial values u 0 and π 0 .

2) At iteration k, a random area A ξ k+1 becomes operational, and each agent j ∈ V ξ k+1 performs the following: a) Update x j , y j and z j using equations (10a), (10b) and (10c) b) Communicate z k+1 j and y k+1 j to the neighboring nodes. c) For each constraint i ∈ I j , compute r k+1 1,i , r k+1 2,i π k+1

1,i and π k+1 2,i using equations (11a) to (12b) respectively.

3) Every non-selected agent j / ∈ V ξ k+1 keeps its former values x k+1 j = x k j , z k+1 j = z k j , ... 4) If the stopping criterion is not satisfied, increase k and go to 2). Otherwise, retain x k+1 j , y k+1 j and z k+1 j .

We remark that this algorithm is effectively randomized in the sense that only the agents of the chosen areas perform computations and exchanges, the others simply keeping their former variables. The performance of this algorithm as well as details about the choice of the area selection sequence will be provided in Section VI. The proof of Theorem 2 is derived from [29, Th. 2] as Asynchronous Distributed Production-Sharing optimization w/ ADMM relies on i) ADMM applied to a sum of two functions that are proper lower semi-continuous and convex by construction (see Theorem 1); and ii) an i.i.d randomization veryfying Assumption 2 in order to match the assumptions of [29, Th. 2]. Again, the convergence speed in the mean square sense is expected to be O(1/k) from [START_REF] Wei | On the o (1= k) convergence of asynchronous distributed alternating direction method of multipliers[END_REF].

As mentioned before, the convergence result is limited to the case where these areas overlap. This is equivalent to having some agents falling under the authority of multiple areas. In power grid networks, this may result in conflicts in decision making as a coordination is thus required between these authorities. In the following section, we extend our algorithms to the case of non-overlapping areas by introducing dummy nodes between the areas. Subsequently, the areas become independent as every agent refers to only one area. 5 It may seem counter-intuitive that the active set of nodes at times k is denoted by ξ k+1 , this index mismatch is justified so that every variable at time k + 1 (denoted • k+1 ) is F k+1 -measurable where F k+1 is the sigmafield induced by (u 0 , π 0 , ξ 1 , ..., ξ k+1 ). Indeed, the only randomness in the algorithm comes from the area selection process and variables at time k + 1 are derived from the last selection variable (ξ k+1 ) in a deterministic manner.

V. EXTENSION TO THE NON-OVERLAPPING AREAS CASE

In this section, we extend our previous result to the case where non-overlapping areas are activated at each iteration. We focus here on the aforementioned Distributed Production -Sharing algorithm but the reasoning below can be easily adapted to a large class of distributed ADMM-based algorithms.

As opposed to Assumption 2, we divide here the graph G (linked to Laplacian L as described in the previous section) into L strictly separated areas A , ∈ {1, . . . , L} where every node belongs to exactly one area. A node will be called a border node if and only if it is connected to a node in a different area, and an inner node otherwise.

Starting from the initial Production-Sharing Problem (1), we aim at decomposing it into L sub-problems, each associated to an area A . Obviously, some constraints couple multiple areas together and thus cannot be assigned to one area. These coupling constraints are related to the connections between border nodes. Inspired by the method of passing adjacent variables [START_REF] Conejo | Multi-area coordinated decentralized DC optimal power flow[END_REF], we propose an approach that transforms our problem so that only the agents of a given area are active at each iteration. For this, we add a dummy node between each pair of connected border nodes of different areas as shown in Fig. 2. These dummy nodes are not associated with any cost function but only serve to rewrite the constraints coupling adjacent areas using a shared variable, that will have to be exchanged as we will see later.

Thus, the changes to be applied to the initial optimization problem are:

• Add chosen dummy nodes;

• Rewrite the coupling constraints accordingly.

Practically, the problem becomes

Non-overlapping Distributed Production-Sharing Problem

min • x, • y, • z∈R • N j∈ • V • f j ( • x j ) subject to ∀j ∈ • V , • r j ≤ • x j ≤ • r j    • y = • x - • d • y = • L • z (13)
where -for all j ∈ V ,

- • V = V • V is the new set of
• L jj = L jj + 1 2 d∼j • L jd , • L jd = -2L ij , • L dd = 4L ij , • L ij = L ij
, if i and j belong to the same area.

With this formulation we can see that Problems (1) and ( 13) are equivalent. Following the reasoning of the previous sections while considering the modified problem, one can derive a new distributed algorithm using ADMM that accounts for the modified functions and matrices. A solution of the original problem of non-overlapping areas can then extracted from a solution of the modified overlapping problem by omitting the entries related to the dummy nodes. One can remark that in the new algorithm, when an area is updated, it only needs local information plus the value of the variables corresponding to the dummy nodes. An iteration is thus now composed of two parts: i) local computations and exchanges in the selected area; ii) communication of border-related values to the adjacent areas.

VI. IMPLEMENTATIONS AND SIMULATIONS

In this section, we implement our algorithms on the DC-OPF problem [START_REF] Glover | Power system analysis and design[END_REF]. Simulations are first carried out on the IEEE-30 bus test system [32] then, on the IEEE-118 bus test system.

In the DC-OPF context defined in II-B, Problem (5) writes

min u∈R 3N F (u) subject to Au = d ( 14
)
where,

A = I -I 0 0 -I L B and d = (p D 1 , . . . , p D N , 0, . . . , 0) T ∈ R 2N .
Note that the cost only depends on the generated power. We consider here quadratic cost functions, more precisely for each generator agent j, we take f j (x j ) = c j x 2 j + c j x j for p g j,min ≤ x j ≤ p g j,max and +∞ otherwise, with the coefficients given in Table I.

In the following, we apply our asynchronous distributed algorithm to the DC-OPF problem using the network of the IEEE-30 bus test system.

A. Overlapping areas

The cost functions (f j ) are defined as above and ρ is set to 1. Eq. (10a) then simplifies to,

x k+1 j = Π [p g j,min ,p g j,max ] 1 2 x k j -c j -π k 1,j + 1 2 (y k j + d j ) 2c j + 1
where Π [p g j,min ,p g j,max ] is the projection onto the interval [p g j,min , p g j,max ]. Firstly, we divide the network representing the IEEE-30 bus test system into L = 3 overlapping areas as per Table II. We compare the synchronous Distributed Production-Sharing algorithm with the asynchronous version in two different settings: i) one randomly chosen area is active at each iteration; and ii) areas A 1 and A 3 are activated together while area A 2 is activated randomly 50% of the time. In the first scenario only one area is activated per iteration. On total, each area is only activated for 33% of the total time. In the second scenario, A 1 and A 3 are activated together, while A 2 is activated 50% of the time.

Secondly, we test the application of the algorithm in the case of nodes failure. That is, the whole network except for random nodes is being updated every iteration. In this scenario, nodes 1 -4, 27 are randomly switched off for 45% of the time and nodes 10 -14, 19, 22, 28 -30 are randomly switched off for 15% of the time. As an example, in a certain iteration when nodes 1 -4 are deactivated, the rest of the network is considered as one area that is currently activated. In another iteration, when nodes 28 -30 are deactivated, the rest of the network, represented by another area, is then activated.

In Fig. 3a, we plot the normalized mean squared deviation (NMSD) to the optimal solution versus the number of iterations. Remark here that while every agent is active in the centralized scheme, only a subset of the agents may be active in the asynchronous version leading to a lower computations per iteration ratio. In Fig. 3b, we check the sufficiency of the power demand in this network by plotting the residual of the power flow constraints.

B. Non-overlapping areas

We take the same IEEE-30 bus system and we divide it into L = 3 non-overlapping areas as given by Table III. We modify this network by introducing 9 dummy nodes on the tie-lines linking two different areas. Using the two scenarios explained in Section VI-A, we plot in Fig. 4a and Fig. 4b the results of applying the Distributed Production-Sharing areas do not overlap. This is mainly because the variables of the shared nodes are being updated more frequently in the overlapping case. Thus, even in the case of a network with multiple authorities, an agreement on the minimum of the global cost can be met using the distributed Production-Sharing w/ ADMM algorithm.

C. IEEE-118 bus test system

We test our algorithms with larger areas using the IEEE-118 bus test system [32] as our network. This network consists of N = 118 buses, 53 generators and 186 branches. We divide it into 3 overlapping areas as in Table IV and we implement the two scenarios described in Section VI-A.

We obtain the plots in Fig. 5a and 5b representing the evolution of the NMSD and the global power demand deficiency (i.e., the residual of the power flow constraints).

These plots illustrate the capability of our algorithms to solve the DC-OPF problem even when the network becomes larger. This feature makes them adequate for large power grid problems.

VII. CONCLUSION

In this paper, we provided distributed algorithms for solving a Production-Sharing Problem. First, a synchronous version based on the ADMM is developed. Then, an asynchronous version is derived; this scheme reduces the number of active agents to a random subset of the original network, allowing for random failures or for deactivation of parts of the network while preserving convergence.

We applied these algorithms to the DC-OPF problem, which can be cast as a Production-Sharing Problem. Simulations were carried on the IEEE-30 and IEEE-118 bus test systems. Multiple scenarios were provided, illustrating the convergence to the optimal network global state despite the presence of random failures and/or inactivated areas.

An interesting topic for future work is to study the asynchronous behavior of the algorithm when the network under study is dynamic. This is the main case when the nodes are equipped with renewable energy sources and storage devices. As for the asynchronous behavior of the algorithm, it would be interesting to study the case where the random coordinate selection sequence is Markovian and not i.i.d. so to model more realistic long-term failures.

APPENDIX A DERIVATION OF THE SYNCHRONOUS ALGORITHM

For every node j ∈ V , let I j = {i : A ip = 0, p = {j, j + N, j + 2N }} be the set of constraints in which node j is involved.

The minimization step (9a) can be written as

u k+1 = argmin x j ,y j ,z j , j=1,...,N N j=1 f j (x j ) + 2N i=1 ρ 2 A ij x j + λ k α ρ -v k α 2 + ρ 2 A ij +N y j + λ k α+N ρ -v k α+N 2 + ρ 2 A ij +2N z j + λ k α+2N ρ -v k α+2N 2 (15)
where α = α(i, j) = 3N (i -1) + j is used as a compact notation to index the values of v and λ corresponding to the element A ij (similarly α+N and α+2N correspond to

A ij +N and A ij +2N ).
This expression is separable into node basis because the component corresponding to the update of x j (respectively, y j and z j ), depends only on the previous values calculated for node j. Thus the iteration on u can be implemented in a distributed manner where each node j ∈ V update its three components x j , y j and z j by solving the following update steps.

x k+1 j =argmin xj f j (x j ) + ρ 2 2N i=1 A ij x j + λ k α ρ -v k α 2 (16a) y k+1 j =argmin yj ρ 2 2N i=1 A ij +N y j + λ k α +N ρ -v k α +N 2 (16b) z k+1 j =argmin zj ρ 2 2N i=1 A ij +2N z j + λ k α +2N ρ -v k α +2N 2 . (16c)
These update steps can be further simplified as follows using the specific structure of matrix A and using compact notations β = β(i, j) = 3N (j -1)+j, γ = γ(i, j) = 3N (j +N -1)+j and ω = ω(i, j) = 3N (i + N -1) + j.

x k+1 j = argmin xj f j (x j ) + ρ 2 N i=1 A ij x j + λ k α ρ -v k α 2 + ρ 2 2N i=N +1 λ k α ρ -v k α 2 = argmin xj f j (x j ) + ρ 2 x j + λ k β ρ -v k β 2 (17)
y k+1 j =argmin yj ρ 2 N i=1 A ij +N y j + λ k α +N ρ -v k α +N 2 + ρ 2 2N i=N +1 A ij +N y j + λ k α +N ρ -v k α +N 2 =argmin yj ρ 2 -y j + λ k β +N ρ -v k β +N 2 + -y j + λ k γ +N ρ -v k γ +N 2 (18)
z k+1 j =argmin zj ρ 2 2N i=1 A ij +2N z j + λ k α +2N ρ -v k α +2N 2 =argmin zj ρ 2 N i=1 L ij z j + λ k ω +2N ρ -v k ω +2N 2 (19)
To further simplify equations ( 17) to [START_REF] Duchi | Dual averaging for distributed optimization: convergence analysis and network scaling[END_REF], we first explicit the update step (9b).

v k+1 =argmin v G(v) + ρ 2 2N i=1 N j=1 A ij x k+1 j -v α + λ α k ρ 2 + A ij +N y k+1 j -v α +N + λ α +N k ρ 2 + A ij +2N z k+1 j -v α +2N + λ α +2N k ρ 2 . ( 20
)
This last equation can be divided into a set of three equations

v k+1 α =argmin vα ρ 2 A ij x k+1 j -v α + λ α k ρ 2 (21a) v k+1 α +N =argmin v α+N ρ 2 A ij +N y k+1 j -v α +N + λ α +N k ρ 2 (21b) v k+1 α +2N =argmin v α+2N ρ 2 A ij +2N z k+1 j -v α +2N + λ α +2N k ρ 2 (21c) subject to N j=1 v α + v α +N + v α +2N = di , ∀i = 1, . . . , 2N.
In order to solve this set of minimization steps we introduce the vector of Lagrangian multipliers π = (π 1 , . . . , π 2N ) T to the set of constraints on the entries of v (i.e., to the set of constraints Ãu = d). After some algebra we obtain,

v k+1 α = 1 ρ (-π k+1 i + λ k α ) + A ij x k+1 j (22a) v k+1 α +N = 1 ρ (-π k+1 i + λ k α +N ) + A ij +N y k+1 j (22b) v k+1 α +2N = 1 ρ (-π k+1 i + λ k α +2N ) + A ij +2N z k+1 j . (22c)
We substitute these results into their corresponding constraints which gives us

π k+1 i = ρ d(i) r i (u k+1 ) + 1 d(i) N j=1 λ k α + λ k α +N + λ k α +2N ( 23
)
where d(i) = 2 + 1 i>N N j =1 1 j ∈Ωi is the degree of the i th constraint (i.e., the count of its nonzero elements), and

r k i =[A] i u k -di is its residual given by r k i =1 i<N N j=1 1 i=j (x k j -y k j ) -d i -1 i>N N j=1 1 i=j +N y k j + N j=1 L i -N j z k j
where 1 i<N = 1 when i < j and 0 otherwise, vice versa for 1 i>N .

As for the dual variable update step (9c), we can substitute it with a set of separated updates of its components λ α . Then, by making use of (22a) we can reduce the computational complexity of the dual update step (9c) by replacing it with a set of updates on π. This is proved as follows,

λ k+1 = λ k + ρ Ãu k+1 -v k+1 ⇒λ k+1 α = λ k α + ρ(A ij x k+1 j -v k+1 α ) (24a)
λ k+1 α +N = λ k α +N + ρ(A ij +N y k+1 j -v k+1 α +N ) (24b)
λ k+1 α +2N = λ k α +2N + ρ(A ij +2N z k+1 j -v k+1 α +2N ) (24c
In the same manner we can prove that λ k+1 α +N = π k+1 i and λ k+1 α +2N = π k+1 i , ∀j = 1, . . . , N. We plug this last result into (22a) (22b) and (22c),

v k+1 α = A ij x k+1 j - r k+1 i d(i) (26a) v k+1 α+N = A ij +N y k+1 j - r k+1 i d(i) (26b) v k+1 α+2N = A ij +2N z k+1 j - r k+1 i d(i) . ( 26c
)
Which can be reduced to

v k+1 α = 1 i=j x k+1 j - r k+1 i d(i) (27a) v k+1 α+N = -1 i=j||i-N =j y k+1 j - r k+1 i d(i) (27b) v k+1 α+2N = 1 i>N L i -N j z k+1 j - r k+1 i d(i) . ( 27c
)
These last updates correspond to the difference between the updated primal variables and the corresponding constraints' residuals. Thus, they can be directly taken into account in the other steps of the distributed optimization algorithm.

As for [START_REF] Losi | On the application of the auxiliary problem principle[END_REF], after some algebra we obtain

π k+1 i = ρ d(i) r k+1 i + π k i . (28)
We adopt a new formulation for the residual r = (r 1 , r 2 ) where, for i = 1, . . . , N

r k 1,i =r k i = N j=1 1 i=j (x k j -y k j ) -d i (29) r k 2,i =r k i +N = - N j=1 1 i=j y k j + N j=1 L ij z k j . (30)
In the same manner, we divide π into π 1 ∈ R N and π 2 ∈ R N , after some algebra we obtain

π k+1 1,i =π k+1 i = ρ d 1 (i) r k+1 1,i + π k 1,i (31)
π k+1 2,i =π k+1 i +N = ρ d 2 (i) r k+1 2,i + π k 2,i (32)
Finally, we can rewrite equations ( 17), ( 18) and ( 19) as Solving the update steps on y j and z j leads to the Distributed Production-Sharing optimization with ADMM listed in Section III.

x k+1 j =argmin x j f j (x j ) + ρ 2 x j + π k 1,j ρ -x k j + 1 d 1 (j) (x k j -y k j -d j ) 2

APPENDIX B NON-OVERLAPPING TO OVERLAPPING ARCHITECTURE

In the original problem, the constraint related to L writes for all j ∈ V , y j = i∈V L ji z i . Furthermore, if j belongs to area , this can be rewritten y j = i∈A L ji z i + k∈V \A L jk z k where the second term represents the coupling relations we want to eliminate thanks to the addition of dummy nodes.

In the new problem, the constraint related to

Lemma 1 .

 1 Let Assumption 1 hold. Then, the solutions of Problem (8) are of the form u = [x ; y ; z ] where (x , y , z ) is a solution of Problem (1).

Theorem 1 .

 1 Let Assumption 1 hold. Then, the sequence (u k ) k>0 = ([x k , y k , z k ]) k>0 generated by Distributed Production-Sharing optimization w/ ADMM converges to u = [x ; y ; z ] where (x , y , z ) is a solution of Problem (1).

Fig. 1 .

 1 Fig. 1. Power grid and its graph presentation divided into 2 overlapping areas.

Theorem 2 .

 2 Let Assumptions 1 and 2 hold. Then, the sequence (u k ) k>0 = ([x k , y k , z k ]) k>0 generated by Asynchronous Distributed Production-Sharing optimization w/ ADMM converges almost surely to u = [x ; y ; z ] where (x , y , z ) is a solution of Problem (1).

Fig. 2 .

 2 Fig. 2. Network of 2 non-overlapping areas converted to overlapping network.

•N

  agents composed of the set of the N original agents V plus the set of the

L

  respectively equal to their original values f j , r j , r j , and d j ;-for all d ∈ from L is a bit more tedious and explained in Appendix B:

Fig. 3 .

 3 Fig.3. Performance of the asynchronous distributed Production-Sharing w/ ADMM algorithm when the IEEE-30 bus system is divided to 3 overlapping areas.

Fig. 4 .

 4 Fig.[START_REF] Bolognani | A distributed control strategy for reactive power compensation in smart microgrids[END_REF]. Performance of the asynchronous distributed Production-Sharing w/ algorithm when the IEEE-30 bus system is divided to 3 overlapping and non-overlapping areas.

Fig. 5 .

 5 Fig. 5. Performance of the asynchronous distributed ADMM algorithm when the network of 118 nodes is divided to 3 overlapping areas.

,

  ∀j = 1, . . . , N.

L

  ij z j -L ij z k j +

L

  For all i, j ∈ V belonging to the same area, we have• L ji = L ji .Noticing that• y d = 0 for all d ∈ • V , we get that 0 = z j , where i and j are the two original nodes between which the dummy node was inserted (and the second equality is due to the fact that L needs to remain symmetric). Node d serves as a dummy node which does not generate or consume, thus • = 0. We assume the symmetry on the edges linking d to i and j, thus• jd = -2L ij , • L dd = 4L ij .

TABLE III IEEE

 III -30 BUS TEST SYSTEM DIVISION INTO 3 NON-OVERLAPPING AREAS

		Area		Nodes	Number of nodes
		A 1		1,2,5-11,17,20,28	12
		A 2		3,4,12-19		9
		A 3		21-27,29,30		9
		0.014		Centralized w/ ADMM
				Overlapping, Asynchronous Distributed ADMM, one area off
		0.012		Non-overlapping, Asynchronous Distributed ADMM, one area off Overlapping, Asynchronous Distributed ADMM, two areas off
				Non-overlapping, Asynchronous Distributed ADMM, two areas off
		0.01
	NMSD	0.008
		0.006
		0.004
		0.002
		0 0	1000	2000	3000	4000	5000	6000	7000	8000
					Iteration k

Although the DC-OPF constitutes a simplified version of the OPF, it is used by several North American Independent System Operators (ISO) for day-to-day operation[START_REF] Giannakis | Monitoring and optimization for power grids: A signal processing perspective[END_REF],[START_REF] Sun | DC optimal power flow formulation and solution using QuadProgJ[END_REF],[START_REF] Stott | DC power flow revisited[END_REF].

We will use the terms nodes and agents interchangeably.

the introduction of the excess variable y is not mandatory but it is kept for clarity purposes.

throughout this paper, we call indicator function, a function which returns 0 when the argument is in some set C and +∞ elsewhere. It is immediate to check that is C is closed and convex, then the indicator of set C is convex and lower semi-continuous.

ACKNOWLEDGMENT This work has been supported by the ERC Starting Grant 305123 MORE (Advanced Mathematical Tools for Complex Network Engineering). The work of F. Iutzeler is also supported by the IAP DYSCO Network (Dynamical Systems, Control, and Optimization).


Fida Azar

Kevin Courtet

Bassil Dekky

Dominique Bonnier

Olivier Dameron

Alain Colige

Vincent Legagneux

Nathalie Theret

Nathalie Théret

Integration of miRNA-regulatory networks in hepatic stellate cells identifies TIMP3 as a key

Keywords: hepatic stellate cells (HSC), extracellular matrix (ECM), The Cancer Genome Atlas hepatic stellate cells, liver fibrosis, miRNA networks, TIMP3 Lay summary

Introduction

Hepatic fibrosis is a dynamic and reversible process of scarring, which leads to excessive deposit of extracellular matrix (ECM) during most chronic liver diseases (CLDs) whatever the etiology including viral hepatitis, alcoholic liver cirrhosis and non-alcoholic fatty liver disease (NAFLD).

Despite recent advances in antifibrotic therapies, [START_REF] Schuppan | Liver fibrosis: Direct antifibrotic agents and targeted therapies[END_REF] CLDs account for 2 million deaths per year. [START_REF] Byass | The global burden of liver disease: a challenge for methods and for public health[END_REF] Understanding the complexity of the underlying molecular mechanisms remains an open challenge for developing new efficient therapies.

Hepatic stellate cells (HSCs) store vitamin A in normal liver but they are also critical players in liver fibrosis. Following injury, they undergo an activation process leading to myofibroblast-like phenotype characterized by proliferation, contractility, chemotaxis and fibrogenic properties. [START_REF] Tsuchida | Mechanisms of hepatic stellate cell activation[END_REF] The molecular mechanisms that regulate HSC activation are driven by changes in microenvironment involving hepatocytes, inflammatory cells and the extracellular matrix (ECM). Such environment generates extracellular stimuli including growth factors and cytokines, cell-cell and cell-matrix interactions that modulate transcriptional gene expression. Epigenetic and post-transcriptional events are also critically involved and the role of microRNAs is now considered central in the regulation of liver fibrosis and HSC activation. [START_REF] Kitano | Hepatic Stellate Cells and microRNAs in Pathogenesis of Accepted Article Liver Fibrosis[END_REF] MicroRNA (miRNAs) are small non coding RNAs interacting with the 3'-untranslated region of mRNAs to promote their degradation and/or repress their translation. 40 miRNAs have been associated with fibrosis in several tissues, including liver, kidney, heart and lung, but more specific tissue signatures were also identified. [START_REF] Vettori | Role of MicroRNAs in Fibrosis[END_REF] In the liver, miRNAs and their potential use as biomarkers or therapeutic tools or targets have been widely documented. [START_REF] Szabo | MicroRNAs in liver disease[END_REF] Regulatory roles of miRNAs have been described in NAFLD, [START_REF] Su | MicroRNAs in the pathogenesis and treatment of progressive liver injury in NAFLD and liver fibrosis[END_REF] viral hepatitis, [START_REF] Sarkar | Hepatitis B Virus Infection, MicroRNAs and Liver Disease[END_REF][START_REF] Shrivastava | MicroRNAs: Role in Hepatitis C Virus pathogenesis[END_REF] alcoholic liver diseases, [START_REF] Xu | MicroRNAs in alcoholic liver disease: Recent advances and future applications[END_REF] cholestatic liver diseases [START_REF] Castro | Cell Death and microRNAs in Cholestatic Liver Diseases: Update on Potential Therapeutic Applications[END_REF] and hepatocellular carcinomas. [START_REF] Chen | Small but Heavy Role: MicroRNAs in Hepatocellular Carcinoma Progression[END_REF] Pathways leading to fibrosis in CLDs are tightly regulated by miRNA-dependent mechanisms. [START_REF] Jiang | The roles of microRNA families in hepatic fibrosis[END_REF] However only 49 genes were reported as targets of down and up-regulated miRNAs in activated HSCs when compared to quiescent cells. [START_REF] Kitano | Hepatic Stellate Cells and microRNAs in Pathogenesis of Accepted Article Liver Fibrosis[END_REF] The present study aimed to develop a whole-genome integrative approach based on miRNAs previously associated with HSC activation in order to identify new targets associated with liver fibrosis.

Methods

Accepted Article

Human Samples

Matching non-tumor liver samples were obtained from patients undergoing surgical hepatectomy for hepatocellular carcinoma (HCC). The patients were 7 women (mean age 61.8±8.6 years) and 57 men (mean age 61.6±8.7 years). Etiology included chronic alcoholism, hepatitis B and C (supplementary table S1). Histological stages of fibrosis were graded according to the METAVIR score: F1, portal fibrosis without septa; F2, portal fibrosis with rare septa; F3, numerous septa without cirrhosis; and F4, cirrhosis. Necroinflammation activity (A) was graded as A0 (absent), A1 (mild), A2 (moderate) or A3 (severe). Liver samples were from the Biological Resources Center (BRC) at Rennes University Hospital. Access to this material was in agreement with French regulations and satisfied the requirements of the local Ethics Committee.

Cell culture

Human hepatocytes and HSCs were isolated from histologically normal specimens of partial hepatectomies from patients undergoing hepatic resection for liver metastases as previously described. [START_REF] Pabic | ADAM12 in human liver cancers: TGF-beta-regulated expression in stellate cells is associated with matrix remodeling[END_REF] HSCs, LX-2 cell line and human Huh7D12 (Huh7) and HepG2 hepatoblastoma cell lines were maintained in Dulbecco's minimal essential medium supplemented with 25 mM D-Glucose, 4 mM L-Glutamine and 10% fetal bovine serum (FBS). LX-2 cells were serum starved for 24H before treatment with IL1- (3ng/ml) or TGF- (5 ng/ml) for 24H in serum free medium.

Mouse model of liver fibrosis

Seven-week-old female C57Bl/6 mice were treated by intraperitoneal injection of CCl4 (Sigma-Aldrich, St. Louis, MO, USA) diluted at 3% v/v in olive oil. Control mice were treated with the vehicle only (olive oil). For acute treatment, a single dose of 0.3 ml/kg of mouse weight was administered to mice. After 4h, 12h, 24h or 7 days, mice were sacrificed. Livers were collected, weighed and treated as previously described. [START_REF] Kesteloot | ADAM metallopeptidase with thrombospondin type 1 motif 2 inactivation reduces the extent and stability of carbon tetrachloride-induced hepatic fibrosis in mice[END_REF] Experiments performed on animals were approved by the local ethics committee (file number 17-1997) in accredited facilities (LA1610002) from the University of Liege.

In situ hybridization

Accepted Article

Formalin fixed, paraffin embedded (FFPE) tissue sections were deparaffinized with xylene (2x5min), rehydrated with decreasing concentrations of ethanol (100-60-30% for 2 minutes each), washed with distilled water and then dried for 10min at room temperature. Localization of mRNA expression in tissue samples was performed using the RNAscope® in situ hybridization (ISH) assay according to the manufacturer's instructions (Advanced Cell Diagnostics Inc., Hayward, Ca, USA). We used RNAscope® mouse probes predesigned by the manufacturer for detection of Mm-Timp3 in the second channel (probe mm-Timp3-C2, Cat No. 471311-C2), and of Mm-Acta2 in the first channel (probe Mm-Acta2-C1, Cat No. 319531). C2 and C1 signals were detected in red (FastRed) and green-blue respectively, using the RNAscope® 2.5 HD Duplex detection kit. Slides were counterstained with 25% hematoxylin and mounted with VectaMount™ Mounting Medium (Vector Laboratories, H-5000).

RNA isolation and gene expression analysis by quantitative real-time PCR

Total RNA and miRNA were extracted either from cells or tissues using QIAzol Lysis Reagent (Cat No. 56302331) and miRNeasy Mini kit (Cat No. 217004) from Qiagen according to the manufacturer's instructions. Reverse transcription was performed from 0.5 µg total RNA using Superscript II (Invitrogen). Real-time quantitative PCR was performed using the qPCR™ Core Kit for SYBR Green I® (EUROGENTEC™) and the ABI Prism 7700 real-time thermo-cycler (Perkin-Elmer, Foster city, CA, USA). Primer pairs for target genes are described in supplementary table S2. For miRNA 200 ng RNAs were reversed-transcribed using the miRCURY LNA miRNA PCR System (Qiagen). The primers were hsa-miR-21-5p (YP0020430), hsa-miR-222-3p (YP0020455), hsa-miR-221-3p (YP00204532), hsa-miR-181b-5p (YP00204530), hsa-miR-17-5p (YP02119304). and hsa-miR-29b-3p (YP00204679) (Qiagen).

Proteomic analyses of ECM from human fibrotic tissue samples

Tissue samples from 18 adjacent non-tumor liver tissues were used for specific extracellular matrix (ECM) extraction. Tissue preparation and ECM protein enrichment were performed according to Dr Naba's protocol. [START_REF] Naba | Enrichment of Extracellular Matrix Proteins from Tissues and Digestion into Peptides for Mass Spectrometry Analysis[END_REF] Briefly, tissue samples were mechanically homogenized and extracted sequentially using the CNMCS compartmental protein extraction kit (Millipore,Billerica, MA). Cytosolic, nuclear, membrane and cytoskeletal proteins were discarded and the final insoluble fraction enriched for ECM proteins was treated for spectrometry analyses. The ECM-

Accepted Article

enriched samples were solubilized in urea, digested with PNGaseF, Lys-C, and trypsin. The resulting peptides were analyzed by LC-MS/MS on Orbitrap Fusion (Thermo Fisher Scientific, San Jose, CA). Mass spectra were interpreted with MaxQuant/Andromeda and MS/MS spectra were searched against the UniprotKB Homo sapiens (CP_Homo-Sapiens) database with fixed modification: carbamidomethyl(C), variable modifications: oxidation (M), oxidation (K), oxidation (P). Peptide identifications were filtered at 1% false discovery rate. Protein digestion and mass spectrometry analyses were performed by the Proteomics Platform of the CHU de Québec Research Center (Quebec, Qc, Canada).

Bioinformatic Tools

CyTargetLinker 17 is a cytoscape plug-in for integrating regulatory interactions in network analysis. We used the Regulatory Interaction Networks (RegINs) provided on the CyTargetLinker website.

RegINs describes networks containing miRNA regulatory interactions that are derived from online interaction databases. In the present study we selected miRTarBase (7.0) that contains MiRNA-Target Interactions (MTIs) experimentally validated. MTIs are viewed as having strong support when they are validated by western blot, qPCR, or reporter assays. Gene annotations were made using Gene Ontology (GO) or DisGeNet, a discovery platform that contains a comprehensive catalogue of genes and variants associated to human diseases.

Statistical analyses

Experimental data are presented as box-and-whisker plots showing distribution of values from at least 3 independent experiments. Significant differences between groups were compared by Mann-Whitney U (two groups) test or Kruskal-Wallis (more than two groups) test or Wilcoxon rank-sum test, when samples are apparied. Statistical calculations were performed with the R software, version 3.1.1. *P < 0.05, **P < 0.01 and ***P < 0.001 were considered significant.

Results

Building miRNA regulatory networks in activated hepatic stellate cells

Accepted Article

Based on the list of miRNA differentially expressed in activated HSCs (supplementary table S3), we searched for all target genes of up-and down-regulated miRNA. Using CyTargetLinker tool, 8243 target genes were identified including 1215 targeted by only up-regulated miRNA, 4862 targeted by only down-regulated miRNAs and 2166 targeted by both up-regulated and downregulated miRNA. In order to reduce the complexity of the graph, we selected miRNA-target interactions (MTI) that were previously experimentally validated according to the criteria from miRTarBase. A miRNA-target interaction is qualified as Strong if experimental evidences including qPCR, western-blot and/or Luciferase reporter assay support MTI. The reduced graph contained 825 target genes including 429 targets of only down-regulated miRNAs and 281 targets of only up-regulated miRNAs (Figure 1). The lists of genes are provided in supplementary table S4 and comparisons between lists are shown as Venn diagrams in supplementary figure S1.

Among the resulting 825 genes, only 81 genes were annotated as liver fibrosis in the DisGeNET database (http://www.disgenet.org/). Among the 49 target genes experimentally validated in activated HSCs 4 (supplementary table S3), 37 were also found by CyTargetLinker using strong criteria, suggesting a good accuracy of the integrative approach (Figure 2). Among these genes, 24 and 11 are targets of down-or upregulated miRNAs, respectively, while 2 are targets of both up-and down-regulated miRNAs. In accordance with the dogma that activated HSCs are responsible for matrix remodeling, downregulated miRNAs led to up-regulation of matrix components (FN1, COL1A1, COL4A1, TNC) and inhibitor of matrix metalloprotease (TIMP2) while up-regulated miRNAs led to down regulation of matrix metalloproteinase (MMP2, MMP9). Similarly, decrease in SMAD7 and increase in both SMAD4 and TGFBR1 expression supported evidences for induced TGF-β signaling and activity in activated HSCs.

In order to validate new target genes predicted by CyTargetLinker, we used transcriptomic data from genes differentially expressed in quiescent and culture-activated HSCs. [START_REF] De Minicis | Gene expression profiles during hepatic stellate cell activation in culture and in vivo[END_REF] For that purpose, we first selected genes showing two-fold changes (>2 and <0.5, for up-and down-regulated genes, respectively) and extracted miRNA-target interactions that matched between prediction and experiment. As shown in Figure 3A and 3B, 34 up-regulated genes and 17 down-regulated genes matched with genes predicted as targets of, respectively, down-regulated and up-regulated miRNAs in activated HSCs. We next focused on these reduced graphs for more in-depth analyzes.

Accepted Article

Down-regulated miRNA networks correlate with extracellular matrix gene overexpression and cell proliferation 22 of the 34 identified up-regulated genes were associated with two major biological processes: extracellular structure organization (GO:0043062) and mitotic cell cycle process (GO:1903047) (supplementary table S5). In accordance with the role of miR29 family members as major regulators of fibrosis, 13 mir29 down-regulation was associated with up-regulation of extracellular matrix genes including SPARC, COL5A1, COL3A1, COL1A1, LOX, LOXL2, TGFB2, ITGB1, DAG1 (Figure 3A). Similarly, IGF1 and LPL which are upregulated during HSC activation were recently identified as targets of mir29. [START_REF] Bouvy-Liivrand | Combinatorial regulation of lipoprotein lipase by microRNAs Accepted Article during mouse adipogenesis[END_REF][START_REF] Kwiecinski | Expression of platelet-derived growth factor-C and insulin-like growth factor I in hepatic stellate cells is inhibited by miR-29[END_REF] Two other mir29 target genes, RNASEL and GPR85 that have not been documented in liver fibrosis to date were also identified. RNASEL is an endoribonuclease implicated in the antiviral response and GPR85 is a member of the G proteincoupled receptor (GPCR) family mainly associated with central nervous system disease.

Beside the miR29-dependent regulatory networks, 8 miRNAs were associated with cell cycle genes (CCNE1, CDK6, CDKN3, CCNA2, CCNE1) and microtubule-regulatory genes (STMN1, KIF22 and AURKB) illustrating regulation of cell proliferation. In agreement with this, miR146a-5p that targets KIF22, CDKN3 and CCNA2 inhibits HSC proliferation. [START_REF] Du | MiR-146a-5p suppresses activation and proliferation of hepatic stellate cells in nonalcoholic fibrosing steatohepatitis through directly targeting Wnt1 and Wnt5a[END_REF] Similarly miR155-5p targeting RAD51 attenuates HSC activation. [START_REF] Dai | MicroRNA-155 attenuates activation of hepatic stellate cell by simultaneously preventing EMT process and ERK1 signalling pathway[END_REF] Interestingly, the reduced graph enlightened the critical role of BIRC5 (alias survivin) which is controlled by four miRNAs (miR165-5p, miR-195-5p, miR-150-5p, miR-355-5p). While the increased expression of BIRC5 was widely documented in chronic liver diseases and HCC, [START_REF] Montorsi | Survivin gene expression in chronic liver disease and hepatocellular carcinoma[END_REF] its regulation in HSCs remained poorly understood. Our study indicates that its increased expression is partly secondary to modifications affecting the abundance of a small number of miRNA.

Identification of a miRNA subnetwork targeting TIMP3

The up-regulated miRNA network was reduced to 9 miRNAs targeting 17 downregulated genes in activated HSCs (Figure 3B, supplementary table S5). Because TGFB1 signal is the driving force of liver fibrosis, it was not surprising to identify down regulation of antagonists of TGFB1 signal including the E3 ubiquitin ligase SMURF1 (responsible for TβRII ubiquitination and degradation in HSC [START_REF] Liu | IQGAP1 suppresses TβRII-mediated myofibroblastic activation and metastatic growth in liver[END_REF] ) and BMP-2 (antagonizing TGF-β1/ROCK-dependent fibrotic signaling [START_REF] Wang | Up-regulation of BMP-2 antagonizes TGF-β1/ROCK-enhanced cardiac fibrotic signalling through activation of Smurf1/Smad6 complex[END_REF] and osteoblast differentiation [START_REF] Spinella-Jaegle | Opposite effects of bone morphogenetic protein-2 and transforming growth factor-beta1 on osteoblast differentiation[END_REF] ). Similarly, TGFBR3 (alias betaglycan) a non-signaling receptor for TGF-β was downregulated during activation of HSCs. [START_REF] Wickert | Corticosteroids stimulate selectively transforming growth factor (TGF)-beta receptor type III expression in transdifferentiating hepatic stellate cells[END_REF]

Accepted Article

The main intriguing observation was the targeting of the Tissue Inhibitor of Metalloproteinases (TIMP3) by five up-regulated miRNAs including miR-21-5p, miR-222-3p, miR-221-3p, miR-181b-5p and miR-17-5p, thereby suggesting a strong negative regulation of TIMP3 in activated HSCs. To our knowledge, the regulation of TIMP3 by miRNAs was not documented in liver fibrosis and HSC to date. Based on literature meta-analysis we observed that miR-17-5p, miR-21-5p, miR-181b-5p and the homologous miR-221/222-3p have already been implicated in the regulation of TIMP3 in other pathologies including cancer (Table 1).

TIMP3 expression is down-regulated in liver fibrosis

In order to explore the regulation of TIMP3 in liver fibrosis, we first measured mRNA levels in adjacent non tumor liver tissues from the whole cohort of 64 patients with HCC. These samples encompass different stages of fibrosis (5 F1, 4 F2, 7 F3, 48 F4) and inflammatory activities (28 A0, 17 A1, 16 A2, 3 A3). We observed a bell-shaped distribution of TIMP3 mRNA levels among fibrosis stages (Figure 4A). Remarkably, TIMP3 mRNA levels were significantly lower in liver fibrosis with highest inflammatory activities (Figure 4B). Because TIMP3 is a secreted protein that binds to extracellular matrix, [START_REF] Brew | The tissue inhibitors of metalloproteinases (TIMPs): an ancient family with structural and functional diversity[END_REF] we used the original approach developed by Naba et al [START_REF] Naba | Enrichment of Extracellular Matrix Proteins from Tissues and Digestion into Peptides for Mass Spectrometry Analysis[END_REF] to analyze its accumulation in ECM-enriched fraction of liver tissues from selected patients with fibrosis at different stages (5 F1, 3 F2, 4 F3 and 6 F4). As shown in Figure 4C, we observed a decrease (although not reaching statistical signification) in TIMP3 protein amounts associated with the highest fibrosis stages (p= 0.08225 between F1 and F4).

To refine our understanding of the regulation of TIMP3 in liver fibrosis, we have investigated its expression in human hepatic cells in culture. TIMP3 was expressed in both in vitro activated human HSCs and the LX-2 cell line stimulated by TGF-β while cultured human hepatocytes and hepatocellular carcinoma cell lines (HepG2 and HuH7) expressed lower levels of TIMP3 (Figure 5A). These observations could appear as counterintuitive since we showed that miRNA network up-regulated in activated HSC is associated with TIMP3 down-regulation (Figure 3B). However, HSCs activation is a complex phenomenon encompassing several phases and depending on the microenvironment. [START_REF] Tsuchida | Mechanisms of hepatic stellate cell activation[END_REF] In agreement with the phenotypic heterogeneity of activated HSCs, cultured-HSCs transdifferentiate into fully mature myofibroblasts while the expression of profibrotic genes is still responsive to TGF-β in LX-2 cells. More recently, using single cell RNA sequencing approach, Krenkel et al [START_REF] Krenkel | Single Cell RNA Sequencing Identifies Subsets of Hepatic Stellate Cells and Myofibroblasts in Liver Fibrosis[END_REF]

described four sub-populations of activated myofibroblasts

Accepted Article

(MFB) in CCl4-fibrotic livers (chronic treatment for three weeks) and further demonstrated that in vitro activation of HSC induced heterogeneity over time.

Since TIMP3 expression was found to depend on the inflammatory context (Figure 4B), we investigated its expression in LX-2 cells stimulated by IL1- a cytokine produced by activated monocytes during inflammatory response. As shown in Figure 5B, IL1- diminished both TIMP3 and TGF-β mRNA levels, while MMP9 expression, used as positive control, was induced as previously reported [START_REF] Robert | Characterization of the MMP/TIMP Imbalance and Collagen Production Induced by IL-1β or TNF-α Release from Human Hepatic Stellate Cells[END_REF] . These observations are therefore consistent with LX-2 cells being in a transient state where gene expression can be modulated in response to both inflammatory cytokines and TGF-β. These responses are obviously disconnected from the complex response to liver injury that occurs in the whole tissue.

We therefore used a CCl4-induced acute liver failure model to investigate Timp3 mRNA levels during the first steps of liver response. Animals were sacrificed at 4-, 12-, 24-and 168-hours after CCl4 injection. At 24h, a decrease in Timp3 expression and a strong increase in IL1- were observed, which perfectly matches with our in vitro studies (Figure 6A). To further explore the regulation of TIMP3 expression in liver fibrosis, we analyzed the expression of miR targeting Timp3 in the same samples. We observed that the expression of miR-17 (p=0.008), miR-181b (P=0.003), miR-222 (p=0.016) and miR-221 (P=0.055) was increased at 24h hour post-injection (Figure 6B). These observations suggest that inflammatory response induced by acute liver injury is associated with downregulation of Timp3 gene expression and upregulation of miRNAs targeting Timp3.

Together our results have shown that five miRNAs upregulated in activated HSC are predicted to target TIMP3, suggesting the existence of a robust mechanism that tightly controls TIMP3 expression. Following liver injury, inflammatory response is associated with HSCs activation and reduced TIMP3 expression as demonstrated by the inverse association of TIMP3 expression with inflammatory activity and IL1-ß expression.

miR-21, miR-221 and miR-222 are markers of poor prognosis in hepatocellular carcinoma

Because TIMP3 expression was recently associated with a favorable survival in HCC, [START_REF] Gu | TIMP-3 expression associates with malignant behaviors and predicts favorable survival in HCC[END_REF] we hypothesized that miRNAs targeting TIMP3 might be inversely associated with survival. Taking advantages of recent studies that integrate miRNA or mRNA expression data from independent sets of HCC [START_REF] Menyhárt | Determining consistent prognostic biomarkers of overall survival and vascular invasion in hepatocellular carcinoma[END_REF][START_REF] Nagy | Validation of miRNA prognostic power in hepatocellular carcinoma using expression data of independent datasets[END_REF] we examined the prognostic value of miR-21, miR-222, mi-R221, miR-181b and Accepted Article miR-17 expression. Using Kaplan-Meier Plotter tool (http://kmplot.com), we selected RNAseq data corresponding to the 379 HCC samples referred as Liver Hepatocellular Carcinoma (LIHC) in The Cancer Genome Atlas (TCGA) database (https://cancergenome.nih.gov/). miR-221, miR-222 and miR-21 expressions are inversely associated with overall survival (Figure 7A), while miR-17 and miR-181b did not show significant association (data not shown). At the opposite, high TIMP3 expression in HCC samples from TCGA database is positively associated with survival. In accordance with the correlation between TIMP3 expression and favorable outcomes, we further showed that downregulation of TIMP3 in HCC is associated with tumor grade and gene signatures from HCC subtypes that include NCIP (The National Cancer Institute Proliferation subtype), SNUR (Seoul National University Recurrence subtype) and RS65, (65-gene risk score subtype) (Figure 7B). To further investigate the correlation between TIMP3 expression and the miRNAs targeting TIMP3 in HCC samples, we selected 290 hepatocellular carcinoma (TCGA samples) corresponding to primary HCC tumors with both miRNA expression data [START_REF] Nagy | Validation of miRNA prognostic power in hepatocellular carcinoma using expression data of independent datasets[END_REF] and gene expression data (RNA-seq normalized read counts) in the "international cancer genome consortium" database (https://icgc.org). Next, we extracted miRNA expression data for these samples from Nagy's study [START_REF] Nagy | Validation of miRNA prognostic power in hepatocellular carcinoma using expression data of independent datasets[END_REF] that integrated several HCC miRNa data sets including TCGA samples and we elaborated a new expression data set to study correlation between TIMP3 and miRs (supplementary table S6).

Spearman's correlation analysis showed that the expression of miR-17 (ρ= -0.29, p=3.0 10 -7 ), miR-21(ρ= -0.25, p=1. 6 10 -5 ) and miR-222 (ρ= -0.20, p=5.0 10 -4 ) was negatively correlated with that of TIMP3. The correlation of TIMP3 expression with miR-181b and miR-221 was negative but not significant (ρ= -0.1, p=0.09). Together these observation support evidence for the implication of miR-17, miR21, miR-181b, miR-221 and miR-222 in regulation of TIMP3 expression in chronic liver diseases.

Discussion

On the basis of miRNA differentially expressed during activation of HSCs, we developed an integrative approach to identify new target genes and regulatory hubs. A major finding from this study was the miRNA subnetwork that targets TIMP3. Unlike other tissue inhibitors of metalloproteases, TIMP3 has been described as a tumor suppressor since decrease in TIMP3 expression is associated with tumor aggressiveness and bad prognosis. [START_REF] Jackson | TIMPs: versatile extracellular regulators in cancer[END_REF] Silencing of TIMP3 has been associated with promoter hypermethylation, loss of heterozygosity or miRNA-dependent regulation. In chronic liver disease, methylation status of TIMP3 is controverted since TIMP3 has Accepted Article been found scarcely or not methylated in tumor, cirrhosis and chronic hepatitis, [START_REF] Lee | Aberrant CpG island hypermethylation along multistep hepatocarcinogenesis[END_REF] while other studies reported association of TIMP3 gene methylation with hepatocellular carcinoma. [START_REF] Lü | Relationship between TIMP-3 expression and promoter methylation of TIMP-3 gene in hepatocellular carcinoma[END_REF] The regulation of TIMP3 expression by miRNAs has been reported in several pathologies (Table 1). However, there are no studies in liver fibrosis yet. Here, we showed that TIMP3 mRNA can be targeted by 5 miRNAs upregulated in in vitro activated HSCs, as compared to quiescent HSCs. In accordance with our study, miR-221-3p, miR-181b-5p and miR-222-3p are up-regulated in cirrhosis, [START_REF] Wojcicka | Next generation sequencing reveals microRNA isoforms in liver cirrhosis and hepatocellular carcinoma[END_REF] miR-21 shows a 2-fold increase in HCV-dependent human fibrosis [START_REF] Marquez | Correlation between microRNA expression levels and clinical parameters associated with chronic hepatitis C viral infection in humans[END_REF] and miR-17-5P was identified as a serum biomarker for HCV-positive cirrhosis. [START_REF] Oksuz | Serum microRNAs; miR-30c-5p, miR-223-3p, miR-302c-3p and miR-17-5p could be used as novel non-invasive biomarkers for HCV-positive cirrhosis and hepatocellular carcinoma[END_REF] Similarly TIMP3 is downregulated in livers from patients with nonalcoholic steatohepatitis and in HSCs treated with advanced glycation end-products which accumulate in patients with diabetes. [START_REF] Jiang | Advanced glycation endproducts induce fibrogenic activity in nonalcoholic steatohepatitis by modulating TNF-α-converting enzyme activity in mice[END_REF] Together these observations suggested that TIMP3 expression was down-regulated in liver fibrosis.

Using HCC RNAseq data from TCGA database, we further showed that TIMP3 expression was associated with overall survival while expression of miR-222, miR-221 and miR-21 predicted bad prognosis. The protective role of TIMP3 expression was also supported by its association with HCC subtypes characterized by low grade and low recurrence (Figure 7B).

When searching in MiRTarBase that contains more than 38x10 4 interactions between 15064 human genes and 2599 miRNAs, only TIMP3 and the tumor suppressor gene PTEN were commonly targeted by the five miRNAs miR-21-5p, miR-81-5p, miR222-3p, miR221-3p and miR17-5p (supplementary figure S2, supplementary table S7). In line with this common regulation, TIMP3 and PTEN were both implicated in miR-221/miR-222-mediated TRAIL resistance [START_REF] Garofalo | miR-221&222 regulate TRAIL resistance and enhance tumorigenicity through PTEN and TIMP3 downregulation[END_REF] and miR-21-5p-mediated doxorubicin resistance. [START_REF] Chen | miR-21-5p confers doxorubicin resistance in gastric cancer cells by targeting PTEN and TIMP3[END_REF] The critical role of PTEN in chronic liver diseases was recently demonstrated in vivo with Pten-null mice being more prone to develop progressive liver fibrosis and HSC activation than their wild type littermates. [START_REF] He | Activation of hepatic stellate cell in Pten null liver injury model[END_REF] In addition, PTEN negatively correlates with liver fibrosis and activated HSCs in rat liver fibrosis. [START_REF] Zheng | Differential expression of PTEN in hepatic tissue and hepatic stellate cells during rat liver fibrosis and its reversal[END_REF] Conversely, over-expression of PTEN inhibited HSCs activation and induced apoptosis. [START_REF] Takashima | The tumor suppressor protein PTEN inhibits rat hepatic stellate cell Accepted Article activation[END_REF] TIMP3 deficiency was previously associated with inflammatory cells accumulation [START_REF] Murthy | Stromal TIMP3 regulates liver lymphocyte populations and provides protection against Th1 T cell-driven autoimmune hepatitis[END_REF] and exacerbated liver injury through enhanced inflammation and IL1-ß expression. [START_REF] Mavilio | A Role for Timp3 in Microbiota-Driven Hepatic Steatosis and Metabolic Dysfunction[END_REF][START_REF] Zhabyeyev | TIMP3 deficiency exacerbates iron overload-mediated cardiomyopathy and liver disease[END_REF] Our observations have consistently shown that IL1-ß and TIMP3 expressions are inversely correlated in the mouse model of CCl4-induced acute liver injury. From a mechanistic point of view, TIMP3 is the major inhibitor of ADAMs enzymes regulating the biodisponibility of cytokines and growth factors activity implicated in inflammatory response. Its decreased expression upon acute response

Accepted Article

to liver injury might therefore contribute to the implementation of the inflammatory response required for activated HSCs-dependent tissue repair. As shown by in situ hybridization patterns (supplementary figure S3), diverse cell types expressed Timp3 with complex dynamics in both healthy and injured regions. The expression dynamics of Timp3 and miRNAs in whole liver are therefore not contributed by the sole hepatic stellate cells. Nevertheless, our observations indicate that IL1-ß inhibits TIMP3 expression in LX-2 cells raising the possibility that repression of this gene is part of an early response of HSCs to inflammation. We have also shown that TGF- induces TIMP3 expression and that IL1- downregulates TGF-ß expression, explaining how it affects also TIMP3 synthesis. A control of TIMP3 production by TGF- might therefore contribute to the maintenance of immune homeostasis, in line with the anti-inflammatory role of TGF- in liver injuries. [START_REF] Schon | Immunomodulatory effects of transforming growth factorβ in the liver[END_REF] Together our results demonstrate that TIMP3 expression can be regulated by miR-21-5p, miR-181b-5p, miR-222-3p, miR-221-3p and miR-17-5p in activated HSC and that TIMP3 expression is inversely correlated to the expression of these miRs, and to the grade and inflammation activity of liver fibrosis. This regulation is associated with IL1--dependent immune response and future challenges will be to integrate the dynamic changes in miRNA networks to propose a model for TIMP3 gene regulation during the development of fibrosis. S3). Menicis's et al [START_REF] De Minicis | Gene expression profiles during hepatic stellate cell activation in culture and in vivo[END_REF] were used to filter miRNA-genes regulatory interaction networks (Figure 1). Statistical significance was calculated using the Wilcoxon rank-sum test (two groups) or the Kruskal-Wallis non-parametric test (more than two groups). HCC subtypes were extracted from [START_REF]Electronic address: wheeler@bcm.edu & Cancer Genome Atlas Research Network. Comprehensive and Integrative Genomic Characterization of Hepatocellular Carcinoma[END_REF] .

Accepted Article

Figure 1 .

 1 Figure 1. miRNA-genes regulatory interaction networks. The list of down-regulated and upregulated miRNA in activated hepatic stellate cells (supplementary table S3) was used as input for CyTargetLinker application. The resulting graphs illustrate the interactions between miRNA and genes. (A) Full graph including all target genes of down-regulated miRNA (blue) and up-regulated miRNA (red). (B) and (C) Graphs were reduced to interactions only supported by strong evidences as defined by CyTargetLinker application either for down-regulated miRNA network (B) or upregulated miRNA network (C). Target genes annotated by the term liver fibrosis in DisGeNet database (Disease id: C0239946) are shown in bold characters

Figure 2 .

 2 Figure 2. Comparison of genes targeted by up-and down-regulated miRNA in activated stellate cells. Genes identified by CyTargetLinker application with strong evidences criteria (green circle, n=825) were compared with target genes of up-(UP-HSC, n=17) and down-regulated miRNA (DOWN-HSC, n=35) previously reported in activated HSCs (supplementary tableS3).

Figure 3 .

 3 Figure 3. miRNA-genes interaction networks of differentially expressed genes between quiescent and cultured HSCs. Genes differentially expressed during activation of HSCs from De

  Resulting reduced networks are shown: (A) Up-regulated miRNA networks. (B) Down-regulated miRNA networks.

Figure 4 .

 4 Figure 4. TIMP3 expression in human liver fibrosis. (A) and (B) Boxplots showing mRNA levels of TIMP3 in liver biopsy samples with various (A) fibrosis stages (5 F1, 4 F2, 7F7 and 48 F4) and (B) inflammatory activities (28 A0, 17 A1, 16 A2 and 3 A3). Statistical significance was calculated using the Wilcoxon rank-sum test. (C) TIMP3 protein levels in liver fibrosis. Specific extraction of extracellular matrix was performed according to Naba's protocol 16 and protein extracts were analyzed by mass spectrometry from human tissue samples at various fibrosis stages (5F1, 3F2, 4F3, 6F4).

Figure 5 .

 5 Figure 5. TIMP3 expression is down-regulated by IL1- in LX-2 cells. (A) Comparative analysis of TIMP3 expression in cultured hepatic cells. TIMP3 and ALB mRNA levels were Accepted Article

Figure 6 .

 6 Figure 6. Timp3 and miR-17-5p miR-181b-5p miR-21-5p miR-221-3p miR-222-3p expression in CCl4-treated mice. Mice were treated by a single dose of CC4 (0.3 ml/kg of mouse weight) and were sacrificed 4-, 12-, 24-and 168-hours post-injection (A) Box plots show mRNA levels of Timp3 and IL1- expression. (B) Box plots show miR-17-5p, miR-21-5p, miR-181b-5p, miR-221-3p and miR-222-3p expression.

Figure 7 .

 7 Figure 7. TIMP3 expression is a marker for good prognosis of hepatocellular carcinoma. (A) Kaplan-Meier survival curves for TIMP3, miR-222, miR-221 and miR-21 expression. Clinical and expression data were extracted from TCGA_LIHC samples (n=364). (B) Box plots show TIMP3 mRNA levels in different HCC grades and subtype from TCGA_LIHC primary HCC samples (n=294): histologic grades (Grades G1, G2 and G3), the Seoul National University recurrence subtype (SNUR, low and high recurrence), the 65-gene risk score subtype (RS65, low and high risk) and The National Cancer Institute proliferation subtype (NCIP, low and high proliferation).

  Figure 1

  Figure 2

  Figure 7

Financial support statement: This work was supported by the Institut National de la Santé et de la Recherche Médicale (INSERM), Université de Rennes 1, the Ligue Nationale Contre le Cancer and the Région Bretagne. BD is recipient of PhD fellowships from the Ligue Nationale Contre le Cancer and Région Bretagne. Acknowledgments We thank Dr. B Turlin and the Biological Resource Centre from Pontchaillou's hospital for providing human tissue samples. We acknowledge C Monseur and the excellent support of the GIGA center at University of Liege, Belgium in which the animal experimentation was carried out.

Accepted Article