G A Maugin

Internal Variables and Dissipative Structures

The formalism of interoal state variables is estabfIJshed when gradients of these variables are involved, thus allowing for a spatial localization of dissipative effects that give rise to dissipative structures. Therefore, evolution-diffusion equ ations, rather than usual evolution equations, are obtained. TfÑs leads to a com parison between dissipative structures and local structures exhibited by conser vative systems. Three illustrative cases are bfŁefĦy sketched out: nematic liquid crystals, localization of damage coupled to elasticity or plasticity, and localiz ation of plastic strains in plasticity with hardening. Thanks to the absorption of some tefés in the extra entropy flux, the approach is valid to any order in the gradients and the accompanying heat equation has always the same form, pro viding in fact a possible basis for the observation of some of the dissipative structures of interest.

Introduction

The use of theféodynamical internal variables of state has been advocated as an efficient means of dealing with relaxation and hysteresis effects in complex ther modynamical systems [START_REF] Bampi | Nonequilibrium Thermodynamics -A Hidden Variable Approach[END_REF][START_REF] Maugin | Thermodynamique a variables intemes[END_REF][START_REF] Maugin | Thermodynamics of Hysteresis[END_REF]. However, several recent works [START_REF] Maugin | Vectorial Internal Variables inMagnetoelasticity[END_REF][START_REF] Grmela | Bracket Formulation of Diffusion-convection Equations[END_REF][START_REF] Aifantis | Remarks on Media with Microstructures[END_REF] emphasize that such internal variables should be determined by "complete balance laws" including both a rate term and a divergence term and not according to usual "evolution laws" where the divergence term is absent and with which plasticity, viscoplasticity, damage and fracture in solid mechanics [START_REF] Maugin | Thermomechanics of Plasticity and Fracture[END_REF] and solutions of macromolecules [8-1 O] have accustomed us. That is, although the origin of theféodynamical interpal variables is to be traced back to the description of reacting chemical systems [START_REF]tive Laws and Microstructure[END_REF] and fĦuids with internal degrees of freedom [START_REF] Muller | Thermodynamics[END_REF], where diffusion of the interqal variable is of no practical importance, it seems that many physical phenomena exhibiting a spatial localization or a marked diff usion of dissipative processes are directly amenable by means of spa tially non-uniform internal variaeles. It is the purrose of the present paper to show that the con temporarx framework of irreversible thermodynamics accommodates directly this type of approach. Then the thefêodynamical force associated with the inter nal variable is none other than the (spatial) Erãer-Lagrange derivative of the free energy răth respect to this variable, and the extra entropy flux is shown to include contributions from the flux of internal variables. Of necessity, thermal equations need to be developed simultaneously.

Short illustrations of the general theory are proposed which pertain to nematic liquid crystals, brittle or ductile damage in defofłable solids, and plastic strain localization in plasticity with hardening. Together with nonlinearities the spatial gradients thus introduced allow on e ' to exhibit a general structure which is none other than the one of so-called dissipative structures (see, e.g., Ref. [START_REF] Vasilev | Autowave Processes in Kinetic Systems[END_REF]). These are markedly difÏerent from sofIJtonic structures which rather follow from a de scription by means of thermodynamical observable state variables relevant to conservative systems. Essential difÏerences between these two types of structures are emphasized. Finally, it is noticed that the heat equation in the presence of dissipative structures offers a direct means to appraise the dissipation carried by these structures through experimental techniques such as infrared thermography.

Contents

The general balance laws and basic thermodynamical statements are recalled in Section 2. State laws and the residual dissipation inequality are given in Section 3 when spatial gradients of the internal variables are taken into account. The corre sponding general form of the heat equation is obtained in Section 4. In Section 5 we show how to exploit the dissipation inequality by using dissipation potentials. The ensuing general fofê of dissipative structures is given in Section 6. These structures and solitonic ones are contrasted in Section 7 by emphasizing that these two difÏerent types of localized structures follow from two difÏerent limit processes.

Three examples pertaining to nematic liquid crystals (the internal variable is a vector field), damage coupled to elasticity or plasticity (the internal variable is a scalar damage parameter) and plasticity with hardening (the internal variable is the cumulated plastic strain) illustrate our formal considerations in Section 8. General concluding comments are given in Section 9.

We shall use the notation, and work in the spirit, of our forthcoming book on the thefêomechanics of plasticity [START_REF] Maugin | Thermomechanics of Plasticity and Fracture[END_REF].

General equations

We consider a continuous medium which is a simple material in appearance, conducts heat and, to fix ideas more precisely, we specify that its thermodynami cally recoverable behaviour will be of the elastic type in small strains (this will be relaxed in one of the examples). Irreversible processes, other than heat conduc tion, taking place in the medium will be taken care of through a thermodynam-ical internal variable ex and its spatial non uniformity (e. g., its first gradient �ex).

Here ex denotes an n vector of f'ltn, i.e., it is the ordered list of then independent components of a certain geometrical object, a scalar (n = 1) yielding the simplest case. As our medium apparently is a simple material, the following Euler Cauchy equations of motion hold true as the notion of simple material is stfŁctly related to that of Cauchy's stress ( [START_REF] Maugin | The Principle of Virtual Power in Continuum. Mechanics: Application to Coupled Fields[END_REF]; here T = transpose):

{}V =diva+/, a=aT a • n = T on o Q, a • n = 1' u =-Ii in 0, on oQT, on auu, (2.1) (2.2) (2.3)
wherein e is the matter density, vis the velocity f®eld, a is Cauchy's stress tensor,/ and 1' represent volume and surface prescribed forces, -Ii is a prescribed displace ment f®eld, and n is the unit outward normal to the regular boundary f}Q = f}QT u f}Qu of the open simply connected set Q of 8 3 occupied by the material at time t. In the small strain approximation v � ri, where a superimposed dot denotes the material time derivative hereby reduced to a partial time derivative.

The thermodynamical aspect of the problem is approached through the two laws of thermodynamics here wfŀtten in global form for Q as:

First law of thermodynamics:

K + E = @>( ex t) + Q;
Second law of thermodynamics:

where 1 K = L 2 e v 2 dv , E = L eedv &'(ext) = J f • v dv + J T • v da , Q = -J q • n da, u an an % = J (}17dv, N = -J S • n da . n an (2.4) (2.5) (2.6)
Here Kand E are the total kineti � energy and internal energy, &> < ext> is the power developed by prescribed forces, Q is.the total flux of heat, % is the total entropy, and N is the total flux of entropy. Because of complex dissipative mechanisms occurring in the material we assume that S and q are related by the somewhat unconventional relation (cp. Ref. [START_REF] Muller | Thermodynamics[END_REF])

(2.7)
where e is the thermodynamical temperature ( e > 0, inf e = 0), and k is an extra entropy flux density to be determined in the subsequent procedure.

Accounting for the scalar product of equations (2.1) and (2.2)1 with v in equation (2.4), we are led to the classical energy equation where

E + &> <i > = Q , &><i> = -J a : s dv, a 1 . 8 = (Vu)s = 2 [Vu + (Vu)T].
On localization this yields at all regular points in Q:

t = P <i > - V .
q' P <i >tr : . 8'

(2.8)

(2.9)

(2. 10)

where E = ee is the internal energy per unit volume (we account for the small strain hypothesis), and P <i > is the elementary power developed by the Cauchy stress. The localization of equation (2.5) directly yields the local entropy inequ ality as (S = erJ)

S + V•S>O (2. 11)
or (as e > 0)

es + v • ( e S) -(S • v) e > o .
(2. 12)

Now we combine (2.10) and (2.12) in the usual manner by introducing the Helm holtz free energy per unit volume by 'l'=E-Se .

(2.13)

The following Clausius Duhem (CD) inequality fo llows at once:

-(P + s B) + P(i) + v. (e k) -(S. V) e > 0.

(2.14)

State laws and dissipation inequality

We need to specify the functional dependence of 'I' to proceed with the exploit ation of the CD inequality (2.14). It is only at this point that internal variables enter the picture. We assume (small-strain hypothesis) that s splits additively in an elastic part s e and .an anelastic part s P :

(3.1)

In the spirit of plasticity and viscoplasticity, we suppose that tp depends on sonly through se, and it further depends on e and ex, and possibly on the first gradient of the latter, i. e., Vex (see Appendix for the second gradient). 'I' is assråed to be concave in 8 but it is not necessarily convex in ae and ex, although it may be so in many special cases. Convexity with respect to J7 ex holds good. Thus, symbolically,

a J7 ex .91 = A -J7 • B = -l>'P/l>rx., l> a a - = --J7 •- {> ex a rx.
a J7 rx. Verw often this is split in two parts (resulting thus in stronger conditions):

cpintr = Cl : BP + ,9/ cX > 0, <P1 h = -(S • V )e > 0, (3.11)
. (3.12) where, in some sense, we recognize the different qualitative nature of the two classes of dissipative processes. <P i n tr and <P 1 h are the intrinsic and thermal dissip ations, respectively. The split (3.11)-(3.12) will be assumed in all applications.

At this point we have at hand the thermodynamical state laws (3.8) and the defĥnition of d and the explicit form of S, and the inequalities governing dissipa tive processes, (3.10), or (3.11) and (3.12). The latter are in the standard bilinear form in terms of fluxes and associated forces:

(3.13)
as used in standard irreversible thermodynamics [START_REF] De Groot | Non-equilibrium Thermodynamics[END_REF]. Before exploiting (3.13) we examine the general expression of the heat equation.

Heat equation

The heat equation is none other than equation (2.10) 1 in a disguise. Indeed, on using (2.13), (3.3) and (3.8) in (2.10)1 or, equivalently, just comparing (2.12) and (3.10), we deduce the fĥrst (general) form of the heat equation:

e S + J7 • ( e S) = <Pi n tr,
where the intrinsic dissipation acts like a body source of heat. where the non-negativity of the specific heat C follows from the concavity of 'P with respect toe, (iJ2 'f' /iJ e 2 s 0). -r, together with I and m, are thermoelasticity coeflicients. Equation (4.2) emphasizes the role of body heat source played by both <Pinar and <P,h once a simple constitutive law (such as Fourier's) is assumed for the flux S (and not for q).

Exploiting the dissipation inequality

Leaving aside the case of S for which a simple law Ziegler [START_REF] Ziegler | An Introduction to Thermomechanics[END_REF].

Thermodynamical equilibrium is defined by the condition of spatially uniform temperature and the two conditions:

i,P = 0, d=O.

(5.2)

Therefore, it is natural to postulate the existence of a dissipation potential a positive function which is homogeneous of a cerr¿in degree in a and ci, and is convex in these two variables so that £P and d be given by the derivatives is homogeneous of degree one in tP, then the dissipative mechanical behaviour will be of the plastic type (no strain-rate dependence; no time scale). The homo geneity property is conserved in the Legendre-Fenchel transform so that t'J)r( ;;u a, ex, • . • = up a . s - ;;u s , ex, •.• sP is also homogeneous of degree one in a.

(5. 5 )

If both� and�* are homogeneous of degree two in a or if P , then the complemen tary constitutive equation ( 5 . 3)1 will exhibit a characteristic (relaxation) time while there may exist a convex set limiting the values of a so that we have vis coplasticity (cp. Ref. [START_REF] Maugin | Thermomechanics of Plasticity and Fracture[END_REF]). W e shall comment more thoroughly on the second of these. For instance, if � is quadratic in a. and lJ' is highly nonquadratic in ex while remaining quadratic and convex in Vex, equation (6.1)1 will result in an evolution diff usion equation of the general type

(6.2)
where Df is an operator of diffusion (e. g., 17 2 ) and f is highly nonlinear in ex depending, in addition, on the state of stresses and the field of temperature, these two fĥelds being themselves governed by the evolution-diffusion equations (2.1) and ( 4.2) in the bulk of the material. The system thus obtained for the fields v, u

(such that Ii = v, a simple evolution equation by itself),sP, ex and e clearly offers a formidable potential fer the erĊstence of so-called dissipative structures and auto wave phenomena (cf. Vasilev et al. [START_REF] Vasilev | Autowave Processes in Kinetic Systems[END_REF]) which are so fashionable in contemporary nonlinear physics. The essence of the phenomena here is contained in equations (6.1)1 or (5.3)1 and (5.6h, the fofłer of which deserves greater. attention.

Dissipative and solitonic structures

The udissipative., nature of equation (6.1)1 directly stems from the m anner in which it was obtained and, obviously, from the fact that a is considered as an internal variable. As a consequence of this very formulation no boundary con• dition is a priori associated with (6.1)1 while its diffusion like nature, when 'P depends on J7 a, is evident. Thus accompanying boundary conditions will have to be deduced from (6.1)i by using a ••1ocal .. method, e.g., the so.called pill•box method. However, whether a indeed is an internal variable or a more classical observable variable of state often is a matter of decision at the outset from the part of the scientist, i.e., it depends on the respective scales (most often time scales) at which the processes related to internal and observable variables pro• ceed. Had we considered the new variable a as an additional observable variable, then internal powers associa ted with the generalized velocities ci and P ci would In spite of their common appearance, the philosophies underlying the rAEo equ ations (5.7) and (7.1) are qráte different. Obviously, if time scales are such that relaxarÄon effects dominate over inersial ones (this is pracrÅcally always the case in nematic liqrâd crystals where ex is the director field dcp. de Gennes [START_REF] De Gennes | The Physics of Liquid Crystals[END_REF], Kats and Lebedev [START_REF] Kats | Dynamics of Liquid Crystals[END_REF]), then we can write (7.1) as (7.4) emphasiręng with a scaling factor e the fact that there are inertial effects wfÑch will be discarded. But these will generally be second-order derivatives in time, so that (7.4) may present a singular perturbation in the time domain. On the con• dition that nothing dramatic occurs as e goes to zero, (7.4) will render our equ ation for dissipative strtctures obtained from intefëal-variable theory, i.e., (7.5) On the contrary, in many cases described by equation (7 .1) the problem is almost conservative and dissipation produces only a regular perturbation, so that (7 .1) may be rewritten as (7.6) where the e has been shifted to the dissipative contribution. As e goes to zero, the resulting conservative (albeit, in general, nonlinear and dispersive) equation

(7.7)
will follow without ambiguity.

The two writings (7.4) and (7.6) of (7.1) highlight the essential differences be tween nonlinear evolution-diffusion systems that may exhibit dissipative struc tures and auto-waves as e � 0, and nonlinear conservative systems slightly per turbed by dissipation wfÑch may exhibit pure solitonic structures as e � 0. Structures in nonlinear dissipative systems of the first kind and those in nonlinear conservative systems may have some features in common. Solitary excitation fronts (E. g., motion of a phase-transition interface) and pulses propagating with a stable shape (so-called travelling pulses) are examples of structures which may propagate in both schemes. However, the physical nature and the b'ehaviour of these structures are quite diferent in these two schemes, as one can gather from the two limits obtained in (7.4) and (7.6) as e goes to zero. An example of this is the well known fact that, in the second scheme, colliding solitons penetrate through each other and then recover their individuality up to a possible alter ation in phase, whereas in the case of travelling dissipative s tru<;:tures, colfIJsion of two travelling pulses or fronts usually results in their m�tual annihilation. Fur thermore, while shape and propagation veiocity of stationary travelling pulses are practically unsensitive to i ni tial conditions for dissipative structures, both shape and propagation .velocity of solitons are completely fixed by the initial excitation. This was noticed by Vasilev et al. [START_REF] Vasilev | Autowave Processes in Kinetic Systems[END_REF], p. 15.

• • ..

The above lack of distinction between obseruable and internal natures of some state variables unfortl:lllately exhibits a prevailing ambiguity. For instance, per manent electric dipoles in elastic f erroelectrics are governed by an equation such as (7 .6) in which the second term involves strong nonlinearities, dispersion, and coupling with elasticity. Electroelastic solitons propagate in such a structure [START_REF] Pouget | Solitons and Electroacoustic Interactions in Ferroelectric Crystals-II-Interactions between Solitons and Radiations[END_REF]. In that case the electric polarization is viewed as an observable variable (additional degree of freedom). But to study electromechanical hysteresis in ceramics, an internal variable akin to an electfèc polarization (residual polariz-ation vector, cumulated (scalar) electric polarization) can be introduced [�3]. However, diffusion of this last variable is not yet envisaged, but it could be. Similar remarks hold for the magnetic spin density in elastic ferromagnets (e. g., [START_REF] Maugin | Solitary Waves in Elastic Ferromagnets[END_REF] and [START_REF] Maugin | Thermodynamics of Hysteresis[END_REF]).

It seems that the internal-variable thermodynamical approach is the best suited general framework for the first scheme (dissipative structures). This should come as no surprise if we recall the origin of the concept of internal variables in the kinetics of reactive systems [START_REF]tive Laws and Microstructure[END_REF]. Altogether, however, specific results can only be obtained once the physical and geometrical (tensor order) nature of the relevant internal variables has been specified (see below). Here N, the director velocity, is defined as the Jaumann co-rotational time de rivative of the director field. This is a so-called objective time rate (cf. Ref. [START_REF] Maugin | Continuum Mechanics of Electromagnetic Solids[END_REF], p. 83). As a matter of fact, tP is entirely expressed in terms of objective quantities. The last of (8.9) shows that the fluid considered does not obey Euler's hydrody namics at rest unless the distribution of d's is spatially uniform. In fact, it is not difficult to recognize in (8.1)-(8.10) the equations govefëing the thefêodynamics of nematic liquid crystals [START_REF] Muller | Thermodynamics[END_REF][START_REF] De Gennes | The Physics of Liquid Crystals[END_REF][START_REF] Kats | Dynamics of Liquid Crystals[END_REF], where IJ' is specified as Frank's energy density, the director inertia is neglected, and N is written as a fijnear function of d, thus providing an evolution-diffusion equation for the internal variable d. The evolution equation for the fluid velocity (equarÃon of motion (2.1)1) will be highly nonlinear ind as iiv depends on d. The extra entropy flux (8.10) agrees with the one obtained by Kats and Lebedev ([21], p. 23, equation ( 16)) once N is taken fijnear ind.

The above-given scheme would also apply to semi-dilute solutions of elongated (uncoiled) macromolecules when some kind of orientational order prevails (cf. Maugin and Drouot [START_REF] Maugin | Internal Variables and The Thermodynamics of Macromolecule Solutions[END_REF]). It easily generalizes to the case where the internal variable is a symmetric second-order tensor, the "conformation" of not neces sarily elongated macromolecules [10, 25].

B . Scalar internal variable (localization of d amage)

. .

We return to a soli d -like behaviour. Damage is the • reduction in elastic pr,operties (decrease in Young's modulus) as a result of the irreversible growth of microcav ities and microcracks . .In a simple modelling, following the original idea of Ka chanov [START_REF] Kachanov | Time of Rupture Process under Creep Conditions[END_REF], this can be described by a scalar variable D, such that 0 < D < 1; which measures the decrease in material surface transmitting internal forces (stresses) independently of the orientation of that surface (isotropy hypothesis). It appears that damage is localized and it is subjected to diffusion under the action of stresses. The idea of considering D an as internal variable is not new (cf.

Lemaitre and Chaboche, [START_REF] Lemaitre | Mechanics of Materials[END_REF]).

(a) We first consider damage coupled to elasticity for which plasticity effects•��n be discarded. Hence se = s, and we write tp = 'l' ( s, e, D, VD) . (8.11) The formulation of Sections 2 through 6 applies directly with ix= D. A sensible expression for lJJ which steers the whole problem, is given by

1 lJJ = (1 -D) � (s, e) + 2 v (VD) 2 , (8. 12
)
where � is the thermoelastic energy of the undamaged (virgin) material. Thus, (b) For ductile damage, plasticity with hardening can interfere with damage. Thus we take (8.16) where f3 is the internal variable that accounts for strain hardening (cf. Ref. [START_REF] Maugin | Thermomechanics of Plasticity and Fracture[END_REF]). Assume that the three dissipative processes involved in (8.18) present no time scale, so that the rates {iP, p , fJ} derive from a pseudo-potential of dissipation which is the indicator function of a convex set C whose boundary in (a , !14, d ) space has equation f (a , f!4, d) = 0 (cf. Ref. [START_REF] Maugin | Thermomechanics of Plasticity and Fracture[END_REF]). Then we have the evolution equations

u = o lJJ = c1 -D) o � 0 8 0 8 ' 1 . k = e (V D)D ,
. o f where 1 > 0 if/= 0 andj = 0, and 1 = 0 if/< 0 or/= 0 andj < 0. Fpr instance, we can have (W0, y and u0 are material coefficients)

f (a, !!4, d) = ( 1-D) -1 ll a ll+P4+ 0 --u 0 , [ w, (d )y + l] y+1 W0
from which there follows through (8.19) that and

/3( t) = J (lfP: £P)lf2d t'' p = 1 o 1-D • p -1 if • _ i (J#)Y 8 -W' D -(1 -D) W0 '
where ii is the efective stress (see above).

Furthermore, we need to specify fû. For instance, The quantity h = d K/d eP = 82 Y'a./8(eP)2 is called the hardening modulus (cf.

Ref.

[7], Chapter 5). Writing the consistency condition J 0 to determine i allows one to obtain the following differential equation for eP:

(8.37)

Equations (8.37) and (8.33) 1 and the balance laws of continuum mechanics re called in (2.1)-( 2.3) are none other than the equations that govern the localization of plastic strains and shear-band evolution in metals according to Zbib and Aif antis [START_REF] Zbib | A Gradient-dependent Flow Theory of Plasticity: Application to Metal and Soil Instabilities[END_REF]. Here these equations are given a fully thermodynamical basis which justifÐes why only even-order gradients occur in (8.37) -this comes from the defÐnition of /J//J� -while the fĦow rule (8.32) keeps a priori a very usual form. That is, this approach does not suppose an a priori gradient-dependent flow rule -it is the free energy which is gradient dependent (like in phase transition theory a la Ginzburg-Devonshire). Furthermore, it shows that contrary to what is as sumed by Zbib and Aifantis [START_REF] Zbib | A Gradient-dependent Flow Theory of Plasticity: Application to Metal and Soil Instabilities[END_REF], it is h which is involved in (8.37) and not a more complex expression herr involving both hand de 1/ deP and dc2/ d eP. Had we consi dered sP-dependent c1 and c2, expression (8.37) would have been much more involved, not only through the denominator.

Concludb�g remarks

In spite of their relarÂve complexity and somewhat awkward formulation, the examples briefly treated in Section 8 highlight some of the essential features of our "thermodynamically admissible" approach to dissipative structures. In par ticular, apart from the unconventional form of the entropy flux, it appears that the remainder of the usual framework of irreversible theféodynamics is left unchanged. That is, the standard bilinear form (3.13) using only classical fluxes (e.g., (j) is still valid. Simultaneously, the accompanying heat equation keeps its standard form (4.2). Not unlike the case of conservative structures, most of the relevant physics still is in the definition of the force d whose vanishing valuecp. (5.2h -defines static structures in a way that recalls elementary equilibrium solutions for the order parameter in Landau's theory of phase transition.

The second important point is that evolution-difusion equations seem, with some reason, to take over simple evolution equations and this goes in the direc tion indicated by Grmela [START_REF] Grmela | Bracket Formulation of Diffusion-convection Equations[END_REF] of a general analytical-mechanics formulation of both reversible and irreversible processes. In the course of this development, what made the essence of the notion of internal variables [START_REF] Maugin | Thermodynamique a variables intemes[END_REF], seems to be alienated. That is, while intefëal variables are generally conceived, for the least, as not easily measurable, they were rñewed in all cases as non-controllable. However, with the introduction of diffusion, the notion of accompanying boundary conditions enters the picture and thus control of the said "internal" variables through boundaries appears, contrary to the very notion of being "internal". Hence it seems prefer able to consider equations such as (7.5) as resulting from a simplifÐcation of

(3. 2 )

 2 Through a straightforvard calculation, this allows us to obtain ..P in the form where we have defined the following free energy derivatives:

  space" Euler-Lagrange derivative with respect to rx.. Substituting from (3.3) into (2.14) and accounting for (3.1), we obtain As Sand a are assumed not to depend on 8 and se while the remaining coeffi cients in (3.7) may in general depend on their respective factors, and (3.7) has to remain in one sign for any 8 and s e (see the argument on p. 112 in Ref.[START_REF] Maugin | Continuum Mechanics of Electromagnetic Solids[END_REF]At this point it is astute to select the extra entropy flux k as k =!_ Bri.

8 '

 8 so that (3. 7) reduces to the following residual dissip ati on inequality

(4. 1 )

 1 In the present state of progress we can also use the fact that we know S through (3.8h. Hence, evaluating S and accounting for (3.8) we readily obtain the second form of the heat equations asC B+ J7 •(BS)= <Pin t r + <Pi e 'in which we have defÐned the "thermoelastic" dissipation <P te by

S

  = -K(s;e, 8; O!) • J7 8 (5.1) with A T . (K . A) > "A 2 > 0 for any vector A, K > O. is sufficient, we examine the exploitation of the inequality (3.11) in the light of dissipation potentials in the manner of Rayleigh [1 7] or

• p ofĆ 8 =

 8 o a ' ...A =().@ .w aa • (5.3) With Euler's identity for homogeneous functions the dissipation inequafijty is automatically satisf®ed and, thus, (5.1) and (5.3) are thermodynamically admi ssible. The degree of homogeneity specifies the type of dissipative process consi dered. For instance, iffĆ is homogeneous of degree one in ii, d will not in fact depend on ii, and the behaviour with respect to the variable rx. will be of the "plastic" type, exhibiting no time scale. Similarly, for the more traditional behaviour described by (5.3)1, if the partial Legendre-Fenchel transform fš * defined by fš*(iP, a; . .. ) == Sup [d: £P -fĆ (a, ci; . . . )](5.4) "

  For the time being, on account of (3.5) we record the following possible set of constitutive. equations:

  6. Dissipative structures Equation (5.7) reads {)tp + a� -0 with {)ex a a.tp = tp (se, e; ex, vex) �=�(a, a.; 8e, e, ex, vex).

have( 7 . 3 )

 73 contributed to p(l) and &'m-cp.[START_REF] Maugin | The Principle of Virtual Power in Continuum. Mechanics: Application to Coupled Fields[END_REF]-and a new balance law would govern a with, in general, a kinetic energy, and thus an inertia, associated with ex. Then, after complementing this by a simple study of allied dissipative processes the a• field equation would have read (7.1) where (7.2) is the Lagrangian per unit volume and a a The latter is nothing but the total Euler•Lagrange defŁvative with respect to a. Equation (7.1) is the celebrated equation of Rayleigh ((17], Vol.1, p. 78), used by Maxwell in his famous treatise ((19], Vol. 2, p. 396).

  Vectorial internal variable (liqui d crystals) Let a. = { d i; i = 1, 2, 3} a vector f®eld d called director. Mechanically, the cont inuous medium considered does not present any plasticity-like behaviour. It rather is a fluid of the viscous type, although incompressible. It is easy to see that 8 is replaced by the exact rate-of-strain tensor D = (V v)8 and p ( i ) in (2.10) reads <Iv= <I+ pl, (8.1) where pis the mechanical pressure. Equation (3.2) is replaced by 'P = 'Pc e; d, v d) .

(8. 2 ) 4 )

 24 This free-energy function satisfies a space-isotropy requirement and thus, classi cally, we must check the invariance requirement (obtained by invariance of (8.2) through an infinitesimal rotation of the spatial frame; square brackets indicate anti-symmetrization): Through an elementary computation this allows one to show that equation (3.3) is replaced by 'ii = -s ea : Dd . N + J7 -(B . N)' (8.5) wherein Ni= (DJd); = d; -Qii d i, ll = (Jl v)Aii;i = B <ilkl d k,i > ' subscript A denotes antisymmetfŀzation and parentheses denote sym metrization. Proceeding then as in Section 3 one is led to the following residual dissipation

  <P= dD-(S• V)e ?0.

( 8 . 13 )

 813 If there was no difusion of damage d would simply be the thermoelastic energy of the virgin sample (D = 0). Separating intrinsic and thermal dissipations and assuming linear relations for the dissipative laws, renormalizing the evolution difusion equation for D, and writing the time increment form of the stress-strain relations, we finally obtain the following system of evolution-difusion which summarizes the problem:u =v, i=(Vu)5, (! v = div u, ti =a fJ + (1 -D )(E:i--r B), -r:fJ = (�JUI;, )+ IE v 2 D' CB+ V • (0S) = <Pintr + <Pte'where we have set a= u/(1-D), 'r =osoe'

  and Wais a reference e nergy, -r is a characteristic time an d In is a characteristic length. The stress a is the effective stress of damage theory (cp.[START_REF] Lemaitre | Mechanics of Materials[END_REF], p. 347).

  Only the gradient of cx 2 = D is accounted for (compare Paragraph C below). The thermodynamic laws of state are ofû The intrinsic dissipation reads <l>intr = lT : gP + f!4 iJ + d fJ•

1 KM 1 fû= -( 1 -

 111 D) W (e e ) -13 < 1 + M)/M + -v( Jl D) 2 Mand v are material coefficients. From the consistency conditionj • . 0 (the only case where 1 =t= 0 possibly) it is found that .(ll a ll) (1-D)A.= Kp 1 1M + cro(d/Wo),,'(8.24) where ( ... ) denotes the positive part of its enclosure. Then the "evolution" equation for damage will read with • (ll a ll) D = Kpt!M(Wo/d),, + Uo (8.25) This is strongly coupled to a and s P through equations (8.21) and (8.22)1• Equ ation (8.25) is not, per se, an evolution equation as it contains no time scale. Furthermore, the dynamical behaviour described by (8.21), (8.22)1 and (8.25) is highly singular unless unloading never occurs (for unloading, 1 = 0). But if we are only in a loading phase so that II a II varies monotonically in time, er = II a II can be used as a plastic time and for y � 1, (8.25) will render the following equation: where (linear elasticity) • p -p• O' should be capable of describing the localization of damage (D < 1) during the monotone loading of an elastic-plasrÀc material with hardening. C. Localization of plastic strains an d shear-band evolutionWe consider the free .. energy density(8.28) in elastoplasticity with hardening. Spatial gradients of the scalar internal vafŀ able, up to the second order, are considered (See Appendix). The intrinsic dissip ation will read v. o P a. + VP : oPPcx . . We envisage a rate-independent behaviour. Therefore (compare Paragraph B or Ref.[START_REF] Maugin | Thermomechanics of Plasticity and Fracture[END_REF], Chapters 2 and 3) we hav � the "plastic" evolution equations . af £P =A-, oa(8.31)where 1 is the plastic multiplier and/ is the loading funcrÁon. For the latter we take a direct generalization of the Huber-Mises yield critefèon by writing (this is an expefèmental daträ)f (a, d) = 11a11 + d -u 0 • (8.32) Wehave fŠ >O if/::;:;: Oandf= O, and fŠ = Oif/< Oor/= O andf< O.From(S.31) there follows then that and thus, by eliminating 1 between these two equations t -• a.( t) = J (sP: £P)112 dt ' = e P ' so-called cumulated plastic strain (in fact already introduced as p in (8.21)). We may say that the choice of /has granted a physical meaning to the scalar internal variable a.. If 'Pis further specified as (8.35) we may, for instance, obtain din a fofł such as . . (8.36) where c1 and c2 are constant coefficients [otherwise we should keep if

balance equations (7.4) , although there remains the troublesome question of singular perturbation in the latter as e goes to zero.

Finally, in all practical cases of interest in Section 8, the heat equation ( 4.2) provides, like in the elastic-plastic case (Ref. [7], Chapter 12) , the theoretical support for the observation and measurement of temperature distribution due to localized dissipative structures by means of infrared thermography as intemal vafèable dissipation and the allied "thermoelastic" process provide the heat source in this equation. This is to be further elaborated upon.

A last remark is in order concefëing nonlocality. The introduction of gradients to describe damage effects in Paragraph 8B corresponds to a so-called weak non locality. Other approaches [START_REF] Pigaudier-Cabot | Nonlocal Damage Theory[END_REF]- [START_REF] Saanouni | On the Creep Crack-growth Prediction by a Non local Damage Formulation[END_REF] envisage simple fu nctionals over space to describe these effects, although not in the proper thermodynamical framework. Such a framework for strong nonlocality can be found in Eringen and Mau gin [32] , following many works by A. C. Eringen, a strong advocate of nonlocal continuum mechanics. Whether strong or weak nonlocality is needed to account for localization still is a debated matter (in this regard, see Ref. [START_REF] Maugin | Nonlocal Theories or Gradient like Theories: A Matter of Convenience?[END_REF]).

Appendix

In most of Sections 3 through 8, 'I' is assumed to depend on the first gradient of ex only. The question may be raised about what happens when fÑgher-order gra dients, e.g., the second ones V Vex, are considered? It is a straightforward matter to show that accounting for the second gradient yields �'I' d=-�ex =A-V•B+(V® V) :C, while (3.9) will be replaced by as

With this it is shown that (4.1) is left unchanged as 4> in tr is formally left unaltered, d being still the spatial functional derivative. Equation ( 4.2) is also formally unchanged but for 4>te that will contain a second-gradient of <i, i.e., an additional term n • (V ® V) <i within parentheses with n = "'"-OC/oe.

All conclusions concerning the use and interpretation of ( 4.2) apply also to this more general case. The same proof will apply to gradients of order higher than two as (A.1) will generalize automatically while other terms appearing in the computation will enter the divergence term in (A.3) , and this may always be accommodated ink.


Nazmus Salehin

Patrick P L Tam

Pierre Osteil

PRENET: PREDICTIVE NETWORK FROM ATAC-SEQ DATA

Keywords: ATAC-seq, Footprints, Gene Regulatory Network

Assays for transposase-accessible chromatin sequencing (ATAC-seq) provides an innovative approach to study chromatin status in multiple cell types. Moreover, it is also possible to efficiently extract differentially accessible chromatin (DACs) regions by using state-of-the-art algorithms (e.g. DESeq2) to predict gene activity in specific samples.

Furthermore, it has recently been shown that small dips in sequencing peaks can be attributed to the binding of transcription factors. These dips, also known as footprints, can be used to identify trans-regulating interactions leading to gene expression. Current protocols used to identify footprints (e.g. pyDNAse and HINT-ATAC) have shown limitations resulting in the discovery of many false positive footprints. We generated a novel approach to identify genuine footprints within any given ATAC-seq dataset.

Herein, we developed a new pipeline embedding DACs together with bona fide footprints resulting in the generation of a Predictive gene regulatory Network (PreNet) simply from ATAC-seq data. We further demonstrated that PreNet can be used to unveil meaningful molecular regulatory pathways in a given cell type.

Introduction

Assays for transposase-accessible chromatin sequencing (ATAC-seq) is an efficient and robust technique requiring very little starting material (as few as 50,000 cells), making it a reliable tool for studying restricted biological systems (e.g. embryos). [START_REF] Buenrostro | Transposition of native chromatin for fast and sensitive epigenomic profiling of open chromatin, DNA-binding proteins and nucleosome position[END_REF] Generally, ATAC-seq is used to determine accessible chromatin regions in cells. Recently, it has been shown to be a promising tool to detect transcription factor binding sites (TFBSs) by screening footprints within peaks called from ATAC-seq data. [START_REF] Sung | Genome-wide footprinting: Ready for prime time?[END_REF] Footprints are small dips (10-20 bp) observed in read pile-ups (or peaks) and are indicative of reduced transposase activity that may be associated with TF occupancy.

Combined with RNA sequencing (RNA-seq), these data allow scientists to build gene regulatory networks (GRN) at an unprecedented level of complexity. Nonetheless, RNA-seq data can be difficult to generate because of sample scarcity. This limitation prompted the development of a pipeline that can predict multiple levels of gene regulation only from ATAC-seq data to obtain a comprehensive predictive GRN without further sample preparation (Figure 1A).

Current ATAC-seq pipelines do not extract significantly differentially accessible chromatin (DACs) regions, which are important to determine remodelling of DNA tertiary structure between cell types. For instance, ATAC2GRN [START_REF] Pranzatelli | ATAC2GRN : optimized ATAC-seq and DNase1-seq pipelines for rapid and accurate genome regulatory network inference[END_REF] only provides open chromatin regions based on MACS2 program (Figure 1B -blue area). DiffBind 4 extracts DACs but not specifically located to a promoter and TSS of a given gene (Figure1Bred area). Finally, DAStk 5 can be used to extract footprints using HINT-ATAC [START_REF] Li | Identification of transcription factor binding sites using ATAC-seq[END_REF] program but not within a DACs (Figure 1B yellow area). As we can see on Figure 1B, altogether, the three pipelines aforementioned do not lead to a comprehensive picture of the chromatin status within cells. Thus, we decided to combine both the extraction of DACs specifically located to a promoter (PROM-DACs), predictive of a potentially active gene, together with their associated footprints, predictive of a TF bound to the promoter, to build a two-level Predictive Network (PreNet): TFs regulating genes expression at the promoter. Moreover, PreNet comes with a solution to calculate Fold Changes (FC) and p.values associated to DACs between cell types. In addition, we describe that current footprinting tools, such as pyDNAse [START_REF] Piper | Wellington-bootstrap : differential DNase-seq footprinting identifies cell-type determining transcription factors[END_REF] or HINT-ATAC, [START_REF] Li | Identification of transcription factor binding sites using ATAC-seq[END_REF] produce an important number of false positive footprints which can lead to misinterpretation, hence the devaluation of the data. To prevent this, we developed a False Discovery Rate (FDR) method to eliminate these artefact footprints. Combining both approaches (DACs FC and p.values together with FPs FDR) into one pipeline allowed for the inferring of a more reliable predictive GRN from ATAC-seq data.

To illustrate the utility of PreNet, we used a well-known differentiation model from mouse Embryonic Stem Cells (mESCs) to the Definitive Endoderm (DE). [START_REF] Simon | Functional characterisation of cis -regulatory elements governing dynamic Eomes expression in the early mouse embryo[END_REF] We further compared footprints analysis to ChIP-seq data in similar cell types and benchmark the new FDR filtering against ATAC2GRN by comparing Positive Predictive Value (PPV).

Runtimes against the only available pipeline achieving a similar outcome (ATAC2GRN) further illustrate the increased power of the PreNet pipeline. PreNet can be used with any given set of ATAC-seq data and its respective reference genome. On top of providing a program requiring only ATAC-seq data, we improved the statistical analysis for genuine DACs and FPs leading to an increase of predictiveness as demonstrated here.

Methods

The PreNet pipeline is summarized in Figure 1B.

Availability of data and materials

PreNet can be downloaded from the project code repository https://github.com/ChildrensMedicalResearchInstitute/PreNet. PreNet is provided as a series of snakemake pipelines.

The initial set up for the pipeline requires the tools found in Table 1. Set up for these tools must be completed in the config.snakemake file. The default configuration assumes the tools are present in the current PATH but can be modified to point directly to the tool location.

For input, the snakemake pipeline requires the following (placed in config.snakemake):

1) The trimmed paired-end sequencing files for an ATAC-seq experiment in gzipped FASTQ format, within a "FASTQ" directory.

2) The bowtie2 [START_REF] Langmead | Fast gapped-read alignment with Bowtie 2[END_REF] index for the genome of interest. These may be downloaded from NCBI Genbank or manually created.

3) A list of chromosomes to keep for final analyses along with the size of each of the chromosomes (xxxx.chrom.sizes). The default list contains the non-mitochondrial chromosomes for the mouse genome (chr 1-20, chr X, chr Y). The chrom.sizes file can be found from ENCODE, UCSC or made custom.

4) A blacklist of regions to exclude in BED format. The mouse blacklist (ENCFF547MET.bed) can be found on ENCODE and a similar blacklist exists for the human genome.

5) The regions of interest in SAF format. In this case, promoter regions based on the mm10 reference genome are provided.

6) Transcription factor position-weighted matrix file in MEME format for FIMO [START_REF] Grant | Scanning for occurrences of a given motif[END_REF] allocation. These can be downloaded from a transcription factor database of choice.

Raw data extraction and alignment to reference genome

The pipeline can be started from either raw reads or aligned and filtered files. Within this paper, data from Simon and colleagues 8 (GSE94249) were downloaded from the Gene Expression Omnibus repository and analysed completely using the snakemake pipeline. Paired-end reads were aligned to the GRCm38/mm10 M. musculus reference genome using bowtie2 9 (-k 4 -X 2000 -local -mm), allowing for a maximum of 4 multimapping locations per read. Multimapping reads were then assigned randomly using assign_multimappers.py from the ENCODE project. [START_REF] Davis | The Encyclopedia of DNA elements (ENCODE): Data portal update[END_REF] From the aligned data, unpaired and low-quality reads were excluded using samtools, 11 deduplicated using sambamba [START_REF] Tarasov | Genome analysis Sambamba : fast processing of NGS alignment formats[END_REF] and filtered to remove mitochondrial reads.

Differential Open Chromatin regions analysis of Promoter region

Aligned pairs were converted to pseudo-single ended reads and only the 5'-end was used for differential analysis. This represents the site of Tn5 transposase activity. Regions of open chromatin were identified using MACS2 [START_REF] Zhang | Model-based analysis of ChIP-Seq (MACS)[END_REF] (-q 0.05 -no-model --shift -100 --extsize 200) on pooled replicates. Tn5 transposase events within the promoter (-5 kb to TSS)

and TSS (-1 kb to +l kb) regions (Figure 2A) were estimated using the number of 5'-ends within the region. This was achieved using featureCounts 14 (--read2pos 5) and raw counts were imported into R for analysis with DESeq2. [START_REF] Love | Moderated estimation of fold change and dispersion for RNA-seq data with DESeq2[END_REF] Differentially accessible regions were visualised using deeptools2 (Figure 2C) . 15

Footprinting analysis

Aligned reads were combined to give one file per biological condition. De novo footprint analysis in promoter regions (-5kb to TSS) was performed using HINT-ATAC. [START_REF] Li | Identification of transcription factor binding sites using ATAC-seq[END_REF] Sequences from identified footprints, extracted using bedtools, were assigned to a transcription factor using FIMO 10 from the MEME Suite, utilising the JASPAR 2018 vertebrate database [START_REF] Khan | Update of the open-access database of transcription factor binding profiles and its web framework[END_REF] as a reference for transcription factor (TF) position weight matrices. Regions of predicted TF activity were assigned based on the initial promoter SAF file.

False discovery rate filtering for footprints

Footprints were called using HINT-ATAC (Figure 3A) on all promoter regions and subsequently segregated into those falling either within or outside the called peaks (Figure 3C). HINT-ATAC footprint scores were placed in descending order based on the footprint score. The footprints were filtered based on the proportion of footprints found in inaccessible regions of the genome, outside of ATAC-seq peaks (Figure 3B). For an FDR value of τ, the subset of top footprints is returned such that the proportion of footprints outside of ATAC-seq peaks are equal to τ; these are subsequently filtered for footprints that lie within peaks only. To achieve this, FPs lying outside peaks ("outP") were allocated as false positive footprints and a footprint score cut-off was generated by ranking all FPs, then selecting the lowest score that would provide a final list containing less than 1% of FPs located outside of a peak (in other words, less than 1% of false positive FPs). The goal here is to eliminate FPs in a peak that have a lower score than the 1% FPs detected from a region outside of a peak as they are not more significant, hence meaningful, than a FP found between two single isolated reads.

Benchmarking

Comparison against ATAC2GRN 3 was performed using version 2 of the optimised ATAC2GRN bash pipeline for pairedend experiments. The shell script was modified to run on the benchmarking system. Subsequent allocation of promoter footprints to a transcription factor was done using the PreNet pipeline as ATAC2GRN does not provide an allocation pipeline (only genomic location). The runtime benchmark was performed on a server running Centos 7, Linux kernel version 4.6 and 24 threads running at 2.90GHz. Where applicable, 64-bit versions of the programs have been used.

Running times are given as user time as measured by the shell command time.

ChIP-seq validation

Where peak files are available for pipelines, no further analysis was carried out before comparison (GSM1288315, converted to GRCm38/mm10 M. musculus genome using liftOver [START_REF] Kent | Enabling browsing of large distributed datasets[END_REF] ). Otherwise (GSM1782914), reads were aligned to GRCm38/mm10 genome using Burrow-Wheeler Aligner [START_REF] Li | Fast and accurate long-read alignment with Burrows-Wheeler transform[END_REF] and filtered to remove duplicates, unaligned and blacklist region reads. Peaks were called using MACS2 (-q 0.05). Peak files were then filtered to only include peaks within promoter regions (-5kb to TSS). Regions where transcription factors were assigned to footprints via PreNet were compared to ChIP-seq peaks to determine the Positive Predictive Value (PPV) of PreNet pre-and post-FDR filtering.

Cytoscape GRN generation

Final output from the snakemake pipeline, interactions between FPs filtered by FDR < 0.01 and target DACs, were imported into Cytoscape 21 compatible file format for gene network visualisation. Each GRN comprises of FPs located within a DAC with an open TSS. Nodes with no interactions were excluded from the GRN. Gene ontology search on developmental pathway was conducted using Panther v14. 22

Scalability of PreNet

The snakemake pipeline allow for the allocation of threads for the respective processes. As defaults, threaded processes are allocated 4 threads; therefore, for the processing of 6 files, as in the paper, the pipeline can utilise up to 24 cores.

This can be scaled up or down as required. The pipeline allows for the comparison of any number of biological replicates for 2 conditions.

Results

DACs at promoters are predictive of active genes.

We analysed bulk ATAC-seq dataset from Simon and colleagues [START_REF] Simon | Functional characterisation of cis -regulatory elements governing dynamic Eomes expression in the early mouse embryo[END_REF] on mESC (n=3) differentiation toward DE (n=3).

Firstly, we evaluated whether the promoter regions (Figure 2A and 2B) were differentially opened or closed between the two cell states for both strands of DNA (annotated "+" or "-") (Figure 2C); these regions are annotated as PROM-DACs. We restricted our analysis to promoter regions as there is still strong debate on how to properly annotate enhancer targets (e.g. HiC, methylation…) and how they regulate gene expression. We uncovered 1041 differential opening and 1489 differential closing of PROM-DACs during the mESC to DE differentiation (Figure 2B). This method leads to the detection of false positive regions, such as DACs between two samples that do not present proper enrichment -peaks as defined by MACS2 [START_REF] Zhang | Model-based analysis of ChIP-Seq (MACS)[END_REF] indicating that a large proportion of the DACs was firstly selected from solely closed chromatin regions in both samples.

Genuine PROM-DACs were plotted into a heatmap (Figure 2C). As we can see, selected regions display differential accessibility between both samples after applying our filtering method.

To enable PreNet to select genuine DACs, we then analysed regions surrounding the TSS (+/-1kb) and selected for

those with an open TSS, hence with the ability of being transcribed, using MACS2 peak calling software. We assumed that an active transcription site requires an accessible TSS (TSS-DACs, Figure 2A). This new parameter allowed us to further eliminate 61 (4%) and 182 targets (17%) in the ESC and DE samples, respectively (Figure 2B light green circles). Finally, we calculated the TSS-DAC log2 fold changes (log2FC) from raw reads between the two cell types and plotted it against the PROM-DAC log2FC (Figure 2D). As expected, the correlation between TSS and PROM region log2FC in the two cell states was higher than random (r = 0.65) and significant (p < 2.2×10 -16 ) validating that accessibility of the promoters is linked to that of TSSs. This result is also visible in Figure 2C After annotating the DACs regions to genes, we discovered that pluripotency genes were more accessible in mESCs than in DE cells (e.g. Nr0b1, Dppa2, Dppa4, Zfp42, Prdm5 and Klf5) whilst DE cells displayed a reverse pattern for mesendodermal genes (e.g. Smad3, Fgf8, Pdx1, and Mesp1) (Figure 2D) which are expected to be expressed in both cell types. Example of DACs selected by PreNet are shown in gene tracks for Nr0b1, Dppa2 and Smad3 (Figure 2E).

This result confirms that PreNet selection criteria retain genuine cell-type specific chromatin accessible regions.

PreNet Pipeline Improve the Quality of Footprints Prediction.

To unveil trans-regulatory elements binding on open chromatin regions within the promoter, we extracted FPs from promoter regions using HINT-ATAC (see Methods). Predicted protein-DNA interactions were extracted (2,941,546 for the ESC and 2,978,511 for the DE) (Figure 3A). After investigation we discovered that the majority of these FPs were false positive for two reasons: 1) HINT-ATAC (used in DAStk and ATAC2GRN) does not make the distinction between a proper footprint found in a DACs and dips observed at the edge of a peak or two single-reads separated by a small region (Figure 3B); 2) HINT-ATAC scores FPs based on the flanking reads. That scoring method comes with flaws as FPs could be assigned to closed chromatin regions. More importantly, there is no method allowing to set up a threshold to select genuine FPs within a peak. Indeed, FP are scored but there is no probability associated to that score to predict the presence of these FPs on a promoter in a given cell type.

To overcome these issues, a more stringent method to select FPs was designed based on two criteria:

1) A true positive FP is, by definition, a feature observable within an accessible region (Figure 3B). To visualise the impact of our selection criteria on the FPs list, incremental top ranking FPs against scores were plotted (2,941,546 for the ESC -ESC-FPs are shown as an example on Figure 3A). Then FPs in peaks ("inP") were represented on a similar plot (Figure 3C). Only 482,979 FPs, corresponding to 16% of all FPs previously found in mESC (in DE, 714,562 FPs corresponding to 24% of total FPs (Figure 3G)) are detected within a peak.

Surprisingly, some of the false positive targets display either a very high score or a score close to 0 (Figure 3D).

Both these results raise some concerns regarding the reliability of using HINT-ATAC scoring as the sole method for FPs analysis (Strategy used in ATAC2GRN).

2) Applying a False Discovery Rate (FDR) calculations of FPs using inP versus outP as true positive versus false positive qualification, in order to select in a robust manner more genuine FPs. To image our method, FDR was plotted against cumulative selection of FP (first step = 100 top score FPs, second step = 200 top score FPs, etc…) (Figure 3E) (see Methods). According to our observation in Figure 3D, we can see that in the 100 top score FPs 80% are false positive, similarly to when all FPs are selected (step containing all FPs, n= 2,941,547 in Figure 3E). We selected a step containing a maximum of FPs while retaining less than 1% of FPs. By doing so, 213,271 FPs were shortlisted corresponding to 7% of total ESC-FPs originally found using HINT-ATAC (and 8% for DE-FPs) (Figure 3F andG).

3) Finally, we filtered out the outP's FPs from the FDR corrected list (Figure 3F).

The PreNet tool is believed to generate a more genuine list of FPs from a particular ATAC-seq dataset that remains to be validated.

Validation of PreNet and Genuine GRN Construction

Targets of footprints allocated to SMAD2/3 in DE from PreNet were compared to those found from chromatin immunoprecipitation assay followed by deep sequencing (ChIP-seq) from two independent groups: Wang 2017 (GSM1782914) [START_REF] Wang | The p53 Family Coordinates Wnt and Nodal Inputs in Mesendodermal Differentiation of Embryonic Stem Cells[END_REF] and Yoon 2015 (GSM1288315) [START_REF] Yoon | HEB associates with PRC2 and SMAD2/3 to regulate developmental fates[END_REF] , also in the DE. We compared the putative FP targets associated to SMAD2/3 before and after applying FDR correction to the list of targets from the ChIP-seq experiments giving a

Positive Predictive Value (PPV) that is the proportion of SMAD2/3 allocated footprints that are also found in at least one of the ChIP-seq datasets. Prior to FDR correction, 747 of the 2,794 SMAD2/3 allocated footprints overlap with the ChIPseq dataset, so a PPV of 26.7%. After applying FDR correction, the PPV climbs to 41.7% (584 out of 1,401 footprints) (Figure 3H). According to our expectation, PreNet leads to a better prediction of promoter targets physically bound by SMAD2/3 than using only HINT-ATAC scoring methods.

To further validate the predictiveness of our tool, a GRN was built linking PROM-DACs to FPs across the two cell states (Figure 2C). Only FPs binding PROM-DACs were plotted. We found that the ESC GRN retained 502 genes, which were enriched for "LIF stimulation response", gene set responsible for mESC pluripotency maintenance (Figure 3E). On the other hand, DE cells showed a more restricted network (n = 112 genes) with an enriched gene set linked to "liver development", an endoderm derivative (Figure 3F).

These validation steps (increased PPV and enrichment for meaningful GOs) confirmed that the PROM-DACs selection tools combined with an FP FDR correction embedded in PreNet led to the generation of a more genuine GRN.

Benchmarking PreNet Against Available Pipelines Shows Increased Efficiency.

To increase the impact of PreNet, we performed a run time and PPV comparison as the main criteria for assessing the increased power of PreNet over existing tools using the same datasets. We decided to compare its efficacy against the only available pipeline, ATAC2GRN. [START_REF] Pranzatelli | ATAC2GRN : optimized ATAC-seq and DNase1-seq pipelines for rapid and accurate genome regulatory network inference[END_REF] The user runtime for PreNet was 80 CPU-hours. This process takes 6 FASTQ files as input and provides both differentially accessible regions as well as footprints within the promoter. In comparison, ATAC2GRN took 38.8 and 47.5

CPU-hours for each experiment. Allocating the footprints from the whole genome to a transcription factor using the single-threaded FIMO was stopped at 48h real time. Subset footprints found within the promoter region by ATAC2GRN were allocated using FIMO; this process took 53.22 min and 53.87 min for each condition in addition to the ATAC2GRN pipeline running time. In summary, ATAC2GRN is significantly slower than PreNet.

The predictive value of PreNet was benchmarked using SMAD2/3 ChIP-seq data. ATAC2GRN combined with FIMO allocation resulted in a PPV of 33.3%, wherein 688 out of a total of 2,063 footprints allocated to SMAD2/3 by FIMO overlapped with at least one of the ChIP-seq datasets (Figure 3K). This result is higher than the unfiltered PPV for PreNet (26.7%) but less efficient than FDR correction (41.7%).

Overall, the PreNet pipeline is more efficient at predicting bona fide FPs from ATAC-seq data.

Discussion

Current ATAC-seq analysis pipelines include strong biases leading to annotations of false positive hits for accessible chromatin and potential transcription factors binding on promoter regions. We attempted to solve both issues by firstly restricting the analysis to the promoter regions and intersecting accessible regions with the accessibility of TSS. This analysis pinpointed a close correlation between the TSS and promoter accessibility. We used this filtering to correct for regions that do not have an accessible starting site for transcription, which corresponded to a large proportion of the dataset.

Secondly, we developed a novel approach to improve the specificity of footprints analyses by removing false positive hits using scoring method combined with False Discovery Rate filtering. Then, both gene sets were combined to generate a GRN for each cell type. Genes that could be linked to form a GRN display enrichment for gene ontologies generally associated with each cell type. Finally, we demonstrated that our approach leads to a more robust positive predictive score when compared to the only available tool that performs a similar analysis: ATAC2GRN.

Although, only predictive GRN can be inferred from our pipeline, further validation through wet lab techniques are required, but we believe this approach will help in extracting important information from ATAC-seq datasets by refining putative targets to further validate.

We herein described an innovative method to analyse ATAC-seq data comprehensively while refining the obtained gene sets. This allowed us to extract trans-regulation of gene expression using only a single sequencing method. PreNet could be applied to any ATAC-seq dataset and potentially to single-cell data (not tested here) providing they are compared between clusters of cells. PreNet is believed to enhance ATAC-seq analysis power.

  -, or, in other words, DACs between two closed chromatin regions. To segregate positive regions, we filtered out PROM-DACs on the criterion that at least one of the cell types should display an open chromatin region (i.e., a "peak") based on the threshold applied by MACS2 algorithm. From this analysis, we excluded 28% (422/1489) and 46% (475/1041) of regions in ESC and DE samples respectively (Figure2Bgrey circle),

  with most PROM-DACs having an open TSS (yellow regions surrounding the TSS). Conversely, very few TSSs were inaccessible when the promoter region was open (5 in the ESC and 4 in the DE cells were further removed from the final list) (Figure 2blue circles). PreNet filtering led to the selection of 1001 DACs in ESCs and 855 in DE (Figure 2B and 2D).

Figure 1 :A)Figure 2 :A)

 12 Figure 1: A) Schematic diagram of the PreNet pipeline. B) Flowchart of the PreNet pipeline. Dashed area represents the part of PreNet pipeline used in other available software (ATAC2GRN (blue), DiffBind (red) and DAStk (yellow)). Only, ATAC2GRN was used to benchmark PreNet. DAC: Differentially Accessible Chromatin, FP: Footprint, TSS: Transcription Start Site

Figure 3 :A)

 3 Figure 3: A) All footprints from HINT-ATAC algorithms ranked based on their score in ESC as an example. B) Diagram illustrating annotation of footprints using in the PreNet pipeline. True positives are the footprints found within a peak. False positives, as described herein, are found within closed regions (non-peak-called region) or at the edge of a peak. C) Footprints falling within a peak (top panel) and D) outside a peak (bottom panel) with a maximum score around 10,000. E) False Discovery Rate (%) in cumulative selection of Footprints (Step = 100 FPs). Red dashed lines indicate the threshold used in our study where we selected the highest cumulative step containing a maximum of footprints with less than 1% FPR. F) Footprints after FDR discovery showing a maximum score just above 4,000. G) Summary table of footprints after applying peak allocation and FDR corrections. H) Venn diagram showing intersection of targets of Smad2/3 compared to two ChIP-seq datasets before and after filtering (arrow-head). Black colour indicates FDR correction in footprints selection. Red colour FDR correction overlap with Smad2/3 ChIP-seq data GSM1288315. Blue colour FDR correction overlap with Smad2/3 ChIP-seq data GSM1782914. I) GRN generated from ESC data. Red outline highlights genes related to "Response to LIF". J) GRN generated from DE data. Red outline highlights genes related to "Liver development". K) Venn diagram of ATAC2GRN footprints for comparison to that of PreNet (H). DE: Definitive Endoderm, ESC: Embryonic Stem Cells, FDR: False Discovery Rate, fp: footprints, GRN: Gene Regulatory Network, inP: in a peak, outP: oustide a peak, pe: peaks

Table 1 .

 1 Software requirements for PreNet

	Tool	Version tested	Website
	Samtools 11	1.8	https://www.htslib.org/
			http://bowtie-
	Bowtie2 9	2.2.6	bio.sourceforge.net/bowtie2/index.shtml
	Sambamba 12	0.6.7	https://lomereiter.github.io/sambamba/
	MACS2 13	2.1.1.20160309	https://taoliu.github.io/MACS/
			https://github.com/kundajelab/atac_dnas
		Accessed April	e_pipelines/commits/master/utils/assign
	assign_multimappers.py	3, 2018	_multimappers.py
	featureCounts 14	1.5.3	http://subread.sourceforge.net/
	deeptools 15	2.5.4	https://github.com/deeptools/deepTools
	UCSC kentUtils 16	4	http://genome.ucsc.edu/
	R Statistical Software	3.5.1	https://www.R-project.org/
			https://bioconductor.org/packages/releas
	DESeq2 17	1.22.1	e/bioc/html/DESeq2.html

Acknowledgments

The authors acknowledge the University of Sydney HPC service at The University of Sydney for providing high performance computing resources that have contributed to the research results reported within this paper. Funding bodies: Our work was supported by the Australian Research Council (DP160103651), NS is supported by the Australian Postgraduate Award from University of Sydney and CMRI Scholarship, PO is funded by the Sir Norman Greg fellowship and PPLT is a NHMRC Senior Principal Research Fellow (Grant ID 1110751).

Declaration of interest

The authors declare no competing interests