Fabrice Thalmann

A schematic model for molecular affinity and binding with Ising variables

Keywords: 87, 10, Vg Biological information -82, 35, Gh Polymers on surfaces; adhesion -68, 43, De Statistical mechanics of adsorbates

After discussing the relevance of statistical physics in molecular recognition processes, we present a schematic model for ligand-receptor association based on an Ising chain. We discuss the possible behaviors of the affinity when the stiffness of the ligand increases. We also consider the case of flexible receptors. A variety of interesting behaviors is obtained, including some affinity modulation upon bond hardening of softening. The affinity of a ligand for its receptor is shown to depend on the details of their rigidity profile, and we question the possibility of encoding information in the rigidities as well as in the shape. An exhaustive study of the selectivity of patterns with length n < 8 is carried on. Connection with other spin models, in particular spin glasses is mentioned in conclusion.

Introduction

Molecular biology and biochemistry differ from chemistry and physics by the very high specificity of the interactions and the processes that they aim at describing. Words like functions or shapes are used instead of molecules, atoms, forces or fields. These biomolecules are perceived as capable of processing information (mutual recognition) and performing dedicated actions (switches or catalytic reactions). It takes only two orders of magnitude in size, from 0.1 nm to 10 nm, to abandon a world of thermal chaos and vibrations, and to enter a world of dedicated agents and reliable procedures [START_REF] Alberts | Molecular Biology of the Cell[END_REF].

It was suggested by E. Fisher that selective interaction and binding between molecules were primarily a consequence of their complementary shapes, which has since been known as the lock and key paradigm. These specific interactions between highly complementary moieties accounts for the specialized and efficient action of enzymes and often explains how drugs work at the molecular level. The lock and key paradigm was then extended to include final conformational changes that may happen upon binding, a mechanism known as induced fit [START_REF] Koshland | [END_REF]. Since a receptor and its ligand have two complementary shapes, any structural change results in decreasing the binding properties of the pair.

The tremendous increase of known 3d molecular structures (NMR and X-ray scattering) [3,4] and the evergrowing computing capacities has made of the search for complementary ligand-receptor interaction a very intense and competitive research field. Numerical approaches ba-sed on the lock and key principle are commonly known as docking algorithms. The issue of these strategies depend on how successfully are Coulombic forces, hydrogen bonding, solvation properties, electronic densities. . . accounted for [START_REF] Israelachvili | Intermolecular and Surface Forces[END_REF][START_REF] Leckband | [END_REF].

As the size of the interacting bodies increases, it is natural to question whether statistical physics still plays a role in these highly optimized recognition processes. In cells, most molecular interactions are subtly balanced in order to achieve reversibility and to prevent the occurrence of irreversible aggregation. Entropic contributions may help to achieve this balance.

One of the current goal of computational biochemistry, still out of reach, is to predict realistic ligand binding stabilization energies with an accuracy of about 1 kcal.mol -1 (1.6 k B T ), compared with experimental measurements [7]. The quenching of conformational degrees of freedom upon binding, and the subsequent entropic changes must be accounted for when computing thermodynamical association constants with such accuracy.

The immune system provides unrivaled cases of specific mutual recognition. For instance, antibodies are specialized proteins which recognize and bind in an extremely specific and accurate manner to foreign bodies invading a living organism. A recent bioinformatic study of interactions between T-cells and the major histocompatibility complex (MHC) supports the view that selective interactions between peptides may owe more to a delicate balance among many weak additive interactions, rather than to a strong complementary and exclusive mutual binding [8].

Also striking is the phenomenon of allostery. Some proteins have their function subordinated to the presence of an effector. An historically famous example is the transcription regulation of Lac-operon, for which it was demonstrated by Monod and Jacob that the fabrication of the lactose digesting enzymes in e.coli was conditioned to the presence of a significant amount of lactose in the environment of the bacteria [9]. Modern biology teaches us that in the absence of lactose, the protein Lac-repressor binds to a stretch of dna, and prevents the expression of the genes downstream. When lactose molecules bind to Lac-repressor, the protein shape changes, and it looses its ability of binding dna, enabling the expression of the genes under its control. The dna and lactose binding sites are located on distinct regions of the repressor protein, suggesting an action at distance caused by the presence of lactose molecules. Recently, the idea of a simple conformational change of allosteric molecules has been challenged. By studying a schematic mechanical model of Lacrepressor, R. Hawkins and T.C.B MacLeish estimated the contribution of internal, vibrational degrees of freedom, i.e. a change in protein stiffness induced by the ligand [10]. They concluded that positive or negative binding entropy changes ∆∆S were associated to changes (hardening or softening) in the effective spring constants used in their mechanical model of repressor proteins. It is precisely this idea of stiffness and entropic modulation of the binding site efficiency that motivates the present work.

There is, as a matter of fact, a deep and formal connection between statistical physics, recognition, binding, and information theory. When a ligand binds selectively to a patterned substrate, it accomplishes some kind of decoding and reads a piece of information conveyed by its target, the more conspicuous case being the association of complementary base pairs in dna-dna or dna-rna duplexes which is the cornerstone of genetic information processing. In the following Sections, we intend to show that a simple Ising spin chain can be turned into an elementary model for the binding of a ligand molecule onto a patterned receptor, in the limit where both thermal fluctuations and internal entropic degrees of freedom are relevant. Within this framework, one can tune the interactions between ligand and binding site, as well as the internal stiffness of the ligand, and we also consider flexible binding sites. The binary nature of the receptors makes them natural information carriers. After defining the affinity and selectivity of a ligand for its receptor in this situation, we discuss how much dependent are the affinities and selectivities on the stiffness parameters. Then, we proceed by giving examples of non monotonic behaviors of the affinities with increasing stiffnesses. We exhibit a case of decreased affinity upon local stiffening of the ligand, reminiscent from Hawkins and McLeish results. We show that there are affinity biases between receptors with similar shapes but different local rigidities. We finally perform an exhaustive comparison of all pairs of patterns up to a length of 8 monomers, and discuss their intrinsic ability to reliably encode information, which is found to decrease with their length. In the following sections, the word receptor will be used with the same meaning as "binding site", i.e. an object the size of the ligand that binds to it. This linguistic shortcut must not occult the fact that in many realistic cases, the receptor is a much bigger object than the ligand, and the binding site only constitutes a subpart of it. This schematic model does not pretend to accurately describe a realistic experimental situation. However, despite its simplicity, it already displays a fairly rich phenomenology which may find a counterpart in some real cases.

Using Ising-like models for modeling selective binding is not a new proposal. Indeed, Schmid, Behringer and coworkers introduced and performed intensive statistical studies of models with very similar Hamiltonians [11][12][13]. Their model describe the contact between hydrophobic and polar patches belonging to too opposite and complementary binding sites, and is presented as a coarse-grained approach to hydrophilic-hydrophobic interactions that are central to both protein folding and protein-protein interactions. Their spin variables describe the contact between the two rigid moieties and account for short range local rearrangements of the coarse-grained residues. A coupling J between "spin variables" is interpreted as a cooperativity property, while we associate ours to a stiffness parameter. They model 2d small rectangular patches while we discuss more schematic 1d interacting chains; our attempt to encode information in the stiffness pattern in addition to the spatial shape is not found in their work. The model introduced is Section 3 is similar to their approach, but the general approach of Section 4 is original.

Statistical physics of surface-bound receptors binding

Let us consider a number of ligands L in contact with a single binding site M linked to a rigid surface. A balance is reached resulting in an equilibrium association constant K a L,M for the exchange:

L free + M free ↔ L • M bound , (1)
from which the probability p L,M of observing a ligand bound to M obeys:

K a L,M [L free ] C 0 (1 -p L•M ) = p L•M , (2)
with brackets [.] denoting molar concentrations, and C 0 a molar standard reference concentration, making the equilibrium association constant dimensionless, for instance C 0 = 1mol.L -3 . This is a generalization of the law of mass action that one would write for the complexation equilibrium of two species L, M in solution.

K a L,M [L free ].[M free ] = [L • M bound ].C 0 . (3)
The equilibrium constant depends a lot on the solvation (hydration) of L and M, the ionic content of the solution and all the details of close range interactions between species.

As our current goal is to emphasize the role of the internal degrees of freedom of L, we disregard all solvent related interactions by making a kind of ideal solution assumption, and considering that all ligand-receptor interactions are short ranged. This assumption on the solvent is equivalent to saying that all bound and unbound conformations of the ligand receptor pair have the same solvation free-energy. The mutual affinity of such a pair comes entirely from shape and stiffness considerations.

It becomes possible to express the equilibrium constant K a L,M as a partition function ratio :

K a L,M = N A v 8π 2 V 0 Z L•M bound Z L free Z M free . ( 4
) = N A v 8π 2 V 0 C a L,M (5)
Equation ( 4) is a particular instance of the equilibrium association constant derived and presented as eq. ( 13) in ref. [14]. N A is the Avogadro number and V 0 is the volume occupied by one mole in the reference concentration state C 0 . Z L free stands for the Boltzmann-Gibbs sum over all the internal conformations of L, with fixed orientation and center of mass. Z M free is a similar configuration integral over the internal conformations of M when L and M are apart and not interacting. We assume that the configuration integral of the bound complex L • M reduces to a product vZ L•M bound , where Z L•M bound represents the sum over all internal configurations of L in contact with M (with fixed centers of mass and orientations) and the volume v corresponds to all the positions of the center of mass of the ligand L relative to the receptor M that are considered as forming a bound state. In the context of bimolecular associations, v should range up to a few angstrom cube (10 -30 m 3 ). Note that for rigid receptors Z M free equals 1 by construction.

The association constant is obtained by setting p L•M = 1/2 in equation [START_REF] Koshland | [END_REF]. One can then argue that the partition sum per free ligand

8π 2 [L free ]N A Z L free Z M free (6)
equals the partition sum of the bound complex

vZ L•M,bound . (7)
In writing these expressions, we neglected the specific volume change occurring during the binding process. In other words, we assume that the Gibbs ∆G and Helmholtz ∆F thermodynamic quantities coincide. Equation (4) links the thermodynamical affinity constant that can be experimentally determined and the configuration integrals Z L•M,bound , Z Lfree and Z Mfree that are the main focus of this work. Equation (4) can be also written as:

∆G 0 ≃ ∆F 0 = ∆U -T ∆S c , (8)
with

∆F (0) = -k B T ln(K a L,M
) and where ∆U designates the average enthalpic change of moieties L and M upon binding, caused by their change in conformation and mutual interaction, while ∆S c represents the corresponding change in configurational entropy. This expression stands a particular case of

∆G 0 = ∆U + ∆W -T ∆S c , (9)
demonstrated in refs. [7,14,15], where ∆W accounts for the contribution of an explicit solvent to the formation of the ligand receptor pair. In aqueous solutions, both hydrogen bonding and hydration forces originate from specific interaction with the solvent. Strictly speaking, solvent mediated interactions are associated with ∆W in eq. ( 9). However, it is to some extent possible to take them into account by means of an effective hamiltonian and treat these interactions as if they were part of the direct interaction term ∆U . We define the affinity of L for M as C (a)

L,M (eq. 5). For a given number of ligands L and binding sites M, the affinity controls the fraction of bound pairs of molecules (adsorption isotherms) and an increased affinity leads to an increase of bound pairs. However, when many patterns M 1 , M 2 . . . compete for the same ligand, the affinity is not a good assessment of how exclusive is the binding of L for a given M. This is why one needs a relative selectivity parameter for comparing the respective behaviors of a ligand for a target (matching pattern M) and for a decoy (mismatching pattern W):

S r (L, W) = Z L•M bound Z L•W bound • Z W free Z M free , (10)
or equivalently S r = C a L,M /C a L,W . We will be also interested in the absolute selectivity, when comparing the affinity of the matching pattern with the affinity of a complete set of decoys:

S a (L) = Z L•Mbound Z Mfree • 1   W =M Z L•Wbound Z Wfree   , (11)
if the ligand is given a choice between all the possible patterns, or

S ′ a (L) = inf W =M Z L•Mbound Z Mfree • Z Wfree Z L•Wbound , (12)
if only one decoy W is present at a time. The inf operator of eq. ( 12) runs over all the possible patterns W that are different from the matching pattern M, and picks up the best competitor of M, i.e. the one for which the affinity difference is the lowest. Eq. ( 11) is better suited for assessing the affinity of a ligand for its matching pattern, while all the other possible competitors W = M are simultaneously present. This definition of the selectivity depends crucially on the set P of allowed patterns W. Changing this set P results in changing the selectivity parameters S a,r . Using a subset P ′ of P naturally leads to higher selectivities. In practice, a large number of poor affinity decoys can eventually beat a good ligand-receptor pair [16].

A minimal model for flexible ligand and rigid receptor

We now introduce a model containing the basic components of the above discussion: matching, internal degrees of freedom, stiffness, information, and we search for phenomenon such as stiffness dependent affinity and matchingdecoding.

The ligand L is modeled as an articulated chain of n monomers. Each bead i is allowed to occupy only two positions labelled by a binary variable s i = ±1 and the ligand can adopt 2 n distinct internal conformations. Ligand stiffness is enforced by means of next nearest neighbors couplings J i,i+1 = ±J, in the spirit of the original spinglass model by Edwards and Anderson [17,[START_REF] Mézard | Spin glass theory and beyond[END_REF].

The signs of the couplings {J i,j } define the ground state shape (native shape) of L, up to a trivial two-fold degeneracy, while the moduli J = |J i,j | describe the energetic cost associated to bending distortion of the ligand (strain). We assume that ligands have a well defined orientation, say from left to right, and we do not consider the possibility of left-right reversal. This can be justified considering that biopolymers (proteins, nucleic acids) always display such an orientation along the chain. Finally we also disregard the possibility of lateral shift between ligand and receptors, as occurring for instance in the hybridization of dna oligomers. This physically relevant situation increases significantly the combinatorics of the association and efficient algorithms have been designed to tackle these alignment problems [START_REF] Durbin | Biological sequence analysis[END_REF]. In this work we purposely focus more on the thermodynamical stability of two molecules that are prepositioned in the right conformation, or which possess a unique, non-degenerated optimal relative conformation. The outcome of this schematic model is expected to be relevant on length scales of about 1 nm between consecutive "monomers", large enough to invoke some coarsegraining of the underlying molecules, but small enough to preserve the importance of entropic, conformational degrees of freedom.

We now introduce a symbolic notation to describe the ground state of these flexible molecules. For that purpose, one must distinguish between open, free end molecules and cyclic, closed end molecules. Open ligands with n monomers have 2 n internal configurations. There are 2 n-1 possible ground states, each one of them being doubly degenerated due to up-down symmetry. To fully describe the shape and ground state of a ligand, we denote by the letters u and d the position, respectively up or down, of the first monomer on the left. We then associate a p to each antiferromagnetic coupling constant J > 0, and a m to each ferromagnetic one J < 0. A symbol (o) is added at the end to signal an open chain. The ligand L represented on the left of Figure 1 is thus ascribed to the symbol L =u/ppp(o). The corresponding Ising configuration reads s 1 = +1, s 2 = -1, s 3 = +1, s 4 = -1, and the coupling constants are J 12 = J 23 = J 34 = J.

Cyclic ligands with n monomers have spin s 1 and s n coupled with a term J n1 s 1 s n . Cyclization is not really justified by the initial molecular association problem, it is here just a convenient way to simplify calculations and get rid of boundary effects, giving to all monomers the same importance and simplifying the interpretation of results. Cyclic ligands of length n have the same number of configurations as open ligands, but it requires n coupling constants to fully determine their ground state. Cyclic ligands with an odd number of antiferromagnetic p couplings are "frustrated", meaning that no ligand configuration can satisfy simultaneously all the constraints imposed by the J i,j . The ground state conformation is then at least four-fold degenerated. Cyclic ligands with a even number of p have a well defined two-fold degenerated ground state. For instance, the unfrustrated cyclic ligand with the same shape as represented on the left of Figure 1 is coded as L ′ =u/pppp(c), where the suffix (c) reminds of the cyclic character of the molecule. Ligand L ′′ =u/pppm(c) is frustrated with no well defined ground shape. In this schematic approach, the use of open or cyclic ligands is essentially a matter of convenience, as they generate qualitatively similar results.

In the same way, the patterned receptor (binding site) is represented by a sequence of binary values b i , 1 ≤ i ≤ n, each one taking a value ±1. When the ligand comes in contact with the receptor, it gains a negative stabilizing energy -A whenever the monomer position s i and the receptor value b i match. We do not give a penalty to a mismatch situation b i = s i , but this could just be done by shifting negatively of the total configurational energy. The coupling constant A represents a short range interaction, possibly mimicking hydrogen bonding or hydrophobic patches. In both cases, the effective contact parameter A may depend on temperature. The resulting total "Hamiltonian" describing the ligand and the receptor in close contact is:

H r {s i } = n ′ i=1 J i,i+1 s i s i+1 -A n i=1 δ sibi . (13)
The sum runs until n ′ = n-1 for open chains, and n ′ = n, with coupling J n,n+1 = J n,1 for cyclic chains. As δ sb = (1 + sb)/2 for Ising variables, the Hamiltonian:

H r {s i } = n i=1 J i,i+1 s i s i+1 - A 2 n i=1 s i b i - nA 2 . ( 14
)
assumes the form of a random field Ising spin glass with both quench bond disorder J ij and quenched random magnetic field -Ab i /2. However, contrary to usual disordered systems studies, we do not perform here any average over the quenched random bonds, as these couplings contain the relevant information. The related partition function, expressed with the inverse Boltzmann factor β, is

Z L•M bound = {si=±1} exp nβA 2 + n ′ i=1 (-βJ i,i+1 s i s i+1 )+ βA 2 n i=1 (s i b i ) . (15)
Meanwhile, Z M free = 1 and Z L,free is a special instance of ( 15) with A = 0:

Z Lfree = {si=±1} exp n ′ i=1 (-βJ i,i+1 s i s i+1 ) , (16) = 2 n [cosh(βJ) n ± sinh(βJ) n ] ,
result valid for cyclic chains, with a sign + without bond frustration, and a signotherwise. Finally, from the binding free energy ∆F , defined as:

∆F = -k B T ln Z L•M bound Z Lfree Z Mfree , (17)
one deduces the partial enthalpic ∆U and entropic -T ∆S c contributions by numerically differentiating with respect to β.

∆U = ∂(β∆F ) ∂β ; ( 18
)
-T ∆S c = ∆F -∆U , ( 19
)
where it is assumed that A does not depend on temperature (enthalpic contribution).

To describe the shape of rigid receptors, it suffices, in principle, to enumerate the values b i . However, anticipating the case of flexible receptors that will be considered in the coming section, we use for receptors the same convention as for ligands, i.e. a first letter u or d, followed by a list of couplings p and m. Receptors with matching ground state are called M, W being associated with those with mismatching ground states.

Let us illustrate the preceding section with L=u/pppp(c), M=u/pppp and W=u/pmpm. The coupling constants of

L are {J 1,2 = J 2,3 = J 3,4 = J 4,1 = J}, the matching mo- tif corresponds to b 1 = b 3 = 1, b 2 = b 4 = -1 and the mismatching motif to b 1 = b 4 = 1, b 2 = b 3 = -1.
The calculation of the partition functions by enumeration of the 16 configurations of L, or with a transfer matrix method gives:

Z L free = 12 + 2e 4βJ + 2e -4βJ ; Z LM bound = e 4βJ (1 + e 4βA ) + 4(e βA + e 2βA (20)
+e 3βA ) + 2e -4βJ e 2βA ;

Z LW bound = 2e 4βJ e 2βA + 4e βA + 2e 2βA +4e 3βA + 2e 2βA e -4βJ + e 4βA + 1.

From now on, we assume that k B T sets the energy scale, and introduce the dimensionless coupling constants a = exp(βA/2), j = exp(βJ). Affinity and selectivity are rational fractions of a and j.

C a L,M (a, j) = (1 + a 8 )j 4 + 4(a 2 + a 4 + a 6 ) + 2a 4 j -4 12 + 2j 4 + 2j -4 , (21)
and S r (a, j) = (1 + a 8 )j 4 + 4(a 2 + a 4 + a 6 ) + 2a 4 j -4 2a 4 j 4 + 4a 2 + 2a 4 + 4a 6 + 2a 4 j -4 + 1 + a 8 .

(22) We are interested in assessing the role of the stiffness parameter j. In Figure 2, we observe that the affinity of the ligand L for the matching pattern M increases monotonically with j. Figure 3 represents the variation with j of the thermodynamic potentials ∆F , ∆U and -T ∆S c .

At the opposite, the stiffness j reduces the affinity of L for the mismatching pattern W, as represented in Figure 4, with the thermodynamic functions shown in Figure 5. In addition, one notices that the special case j = 1 represents a soft ligand which can adapt to any pattern. Quite naturally, the selectivity between W and M is 1 for j = 1 and tend towards a finite value for j → ∞ (Figure 6). The maximal selectivity depends on the short range contact parameter a and is reached as soon as j ≥ a.

We conclude that stiffness is always favorable when the shape of a ligand and a receptor agree, but becomes unfavorable is a mismatch is present. When ligand and receptor shapes almost agree but not perfectly, there must be an optimal compromise between a very soft ligand j = 1, which precludes any selectivity at all, and a very hard ligand j ≫ 1 which excessively penalizes the mismatches.

Flexible receptors

The next step is to consider flexible receptors M: we expect then soft ligands to beat stiff ligands, as they will better fit the various configurations of M. In our model, M and L play a dual role and it becomes possible to treat ligand and binding site (receptor) on the same footing, by inserting coupling constants K between the "spins" b i . This situation arises when two molecules with similar weight and structure bind together (Figure 7).

H r {s i , b i } = n ′ i=1 J i,i+1 s i s i+1 -A n i=1 δ sibi + n ′ i=1 K i,i+1 b i b i+1 . ( 23
)
The problem can be naturally solved for cyclic ligands and receptors with 4x4 transfer matrices. When all the coupling constants have same absolute magnitude, we define K i,i+1 = η i |K|, k = exp(β|K|), J i,i+1 = ǫ i |J|, j = exp(β|J|) where η i and ǫ i are ±1, to find:

Z L,M bound = a n Tr n i=1 T (ǫi,ηi) (a,j,k) ; ( 24
)
Z L free Z M free = 4 n [cosh(βJ) n ± sinh(βJ) n ] (25) •[cosh(βK) n ± sinh(βK) n ] ;
(the sign ± depending on the bond frustration along the chains) with noncommuting matrices defined as:

T (ǫ,η) (a,j,k) =    aj ǫ k η j -ǫ k η j ǫ k -η aj -ǫ k -η j -ǫ k η a -1 j ǫ k η a -1 j -ǫ k -η j ǫ k -η j ǫ k -η a -1 j -ǫ k -η a -1 j ǫ k η j -ǫ k η aj -ǫ k -η j ǫ k -η j -ǫ k η aj ǫ k η    .
(26) In computing the trace, one actually performs a summation over the 4 n internal configurations. Thermodynamic quantities ∆F , ∆S c and ∆U , given by equations ( 17), ( 18) and [START_REF] Durbin | Biological sequence analysis[END_REF], are then obtained by numerically differentiating the transfer matrix results. One notices that the transfer matrix is entirely built from dimensionless parameters a, j and k. The numerical differentiation of β∆F with respect to β assumes that A, J and K do not depend on temperature, leading to da/dβ = β -1 a ln(a), dj/dβ = β -1 j ln(j) and dk/dβ = β -1 k ln(k). As a result, the thermodynamic quantities derived from equations ( 17), ( 18) and ( 19) are automatically expressed in units k B T . Non purely enthalpic contributions to the contact energy parameter A could also be included in this numerical scheme by using a different prescription for the derivative da/dβ. The external random field b i which was quenched for rigid receptors is now annealed (the random bonds still quenched), and we checked that for large values of K (namely k = exp(βK) ≥ 5) the result for a flexible ligand receptor pair tends to the predictions for the rigid receptor. One can easily convince oneself that there is no difference between a short rigid receptor and a short flexible receptor with doubly degenerated ground state (i.e. no frustration) for which the condition k ≫ j ≫ 1 holds. In this limit, Z M,free ≃ 2, a factor which also appears in Z L.M,bound , leaving C a L,M unchanged. In practice, we regarded k = 10 as sufficient to reach the rigid situation. This corresponds to an energy gap of 3 k B T between the receptor ground state and its first distorted state. Indeed, all the results regarding rigid receptors presented in this study were actually obtained by setting k to large values such as k = 10 and applying eq. (24).

To calculate the partition function of a flexible, open, ligand receptor pair, one replaces the last matrix T (ǫ,η) (a,j,k) by a matrix T (a,1,1) representing the freely oscillating ends. Formula (24) becomes

Z L,M bound = a n Tr T (a,1,1) n-1 i=1 T (ǫi,ηi) (a,j,k) ; ( 27
)
The transfer matrix formalism can be modified if one wishes to pick-up a particular bond coupling J i,i+1 and assign to it a value different from the usual J. In our numerical implementation of the transfer matrix product, we use a 5-letters alphabet {p, m, P, M, .} to describe a ligand pattern L.

-Characters p and m are respectively used for ǫ = 1 and ǫ = -1. -Character . denotes a vanishing coupling constant J = 0 (j = 1).

-Characters P and M represent ǫ = 1 and ǫ = -1, but with a different (usually larger) magnitude of J or K.

Symbols P and M code for some localized hardening of ligand and receptors. Symbol . loosely connects two adjacent and stiffer domains of L. More complex scenarios can be considered, but all are subject to the same limitation, which is that these transfer matrices can deal only with nearest neighbor couplings.

Ligand stiffness and affinity

We now provide a list of examples illustrating various behaviors. All curves represent the affinity variations when the ligand stiffness is increased from j = 1 (soft ligand) till j = 8, i.e. a 2 k B T activation barrier associated to local conformational change (spin reversal).

-Uphill nonmonotonic affinity:

Ligand L =u/ppp.ppp(o) vs M =u/ppmmppm (Fig- ure 8)
The ligand L is made of two loosely connected adjacent domains ppp. The resulting affinity is nonmonotonic, first increasing, then decreasing (Figure 9, a feature previously seen in Figure 4). This behavior emerges from a competition between the matching subdomain, favored by large values of the stiffness parameter j and the mismatching subdomain whose affinity decreases with j. In this particular case the mismatching domain forces the overall affinity to decrease below its starting point, a maximum being reached around j ≃ 1.6. Such a behavior illustrates the concept of optimal stiffness, where the association of a ligand with a (more) rigid receptor requires some tuning of its average rigidity.

-Downhill nonmonotonic affinity: L =u/ppppppp(o) vs W =u/ppmmppm (Figure 10) For a = 2 and k = 10 (rigid pattern), the affinity curve is not monotonic with j, first decreasing, then increasing (Figure 11). The affinity is minimal for a certain value j ≃ 1.6. This behavior in enhanced when a is increased and reduced when the receptor stiffness k is reduced, as seen in Figure 12. The interest of this situation is to be the exact opposite of the preceding situation, with a local affinity minimum for j ≃ 1.6.

-Unbinding upon local hardening: Ligand L =u/ppPpp(o) vs M =u/ppmmp (Figure 13)

The receptor shows a slight mismatch located under a region in L which is locally more rigid (stiffness J > J) than the average coupling. Increasing J while keeping J constant should hamper the capacity of L to accommodate the mismatch, and lead to lower affinity. This is precisely what is seen in Figure 14. The effect is seen for open chains (j = 1.5, j = 2.) and cyclic chains (j = 2). It is reduced if the average chain stiffness is larger (j = 3). If one admits that the stiff bond P is under control of an external agent or effector, the issue is that the affinity of the ligand for its target decreases even though their shape are not altered. Figures 15 and16 show respectively the energetic and entropic contributions to the changes in affinity associated with increasing the stiffness parameter j. One can read from these data the associated relative entropy change ∆∆S = ∆S(j) -∆S(j). Ref [10] suggests that a positive ∆∆S can be associated to the decreased affinity of lac-repressor for dna in the presence of lactose. This change in ∆∆S arises from increasing one spring constant and decreasing another one, in an harmonic model of lac-repressor. In our case, because our model is not harmonic, changes in ∆∆U and ∆∆S cannot be separated and both contribute to changing the affinity. Note that if the hard bond P was located on top of a matching subdomain, the opposite behavior of the affinity with J would occur. For instance the affinity of the pair u/ppPpm(o)-u/ppppm increases with large values of J (not shown).

To sum up, by combining matching and mismatching subdomains and bonds of adjustable stiffness, it is possible to obtain a variety of nonmonotonic behaviors of a ligand affinity for its target. It turns out to be possible to control the mutual affinity by acting on the rigidity of selected bonds, mimicking possibly the influence of a cofactor (or effector) involved in some allosteric mechanism.

Selectivity and information decoding

Let us consider two persons A and B willing to communicate, and in possession of a number N of specific ligandreceptors pairs. A sends a ligand i to B, and B brings this (still unknown to him) ligand in contact with all the receptors available in his library. By using an analytical tool (fluorescence, quartz crystal microbalance, resonant surface plasmon absorption. . . ) B determines the label i of the ligand he has received. By referring to a preestablished codebook, shared with A, B determines the content of the message. Information is conveyed to B by means of selective binding [START_REF] Lehn | Supramolecular Chemistry: concepts and perspectives[END_REF].

Information is coded in the shape of molecules, and especially biomolecules. When the biomolecule is in solution, this information about shape and chemical composition is made available to other molecules and a few of them will be able to selectively bind to it. Any selective adhesion process can be considered as information reading, or to be more precise, information decoding.

There are a number of issues raised by this medium of communication. One may be concerned by the reliability of the message transmission, directly connected to the absolute selectivity properties of the ligand-receptor pairs S a . One may also asks oneself whether there are physical bounds to the minimal binding energy or contact area necessary to ensure a transmission error rate lower than a predefined threshold value [START_REF]Maxwell's Demon 2: Entropy, Classical and Quantum Information[END_REF][START_REF] Jones | Information and Coding Theory[END_REF].

Coding and decoding information is every day's concern for radio engineers, and many coding schemes have been developed in order to safely carry information around.

The simplified model discussed in the previous section establishes a convenient connection with the realm of digital information, as ligand and receptors are already described in terms of binary chunks of information. Results are shown in Figure 18. It seems that the selectivity bias is indeed lower for d H = 1 than it is for d H = 2. One expects two qualitatively different limits for the selectivity dependence in the number of mismatches. If the stiffness J is larger than the contact energy A, the selectivity is likely to correlate better with the Hamming distance d ′ H between stiffness patterns: one can flip an entire interval of spins by changing only two characters in the stiffness pattern. In the opposite limit, the distance between configurations d H should be dominant. Ref. [23] deals in details with related issues.

Selectivity of finite length symbols

We calculated the absolute selectivity (eq. 11) among the subset of all n-long ligand and receptor symbols, to determine the true amount of information that such a family of ligands may carry.

To be more precise, each ground state conformation L among the 2 n-1 possible ones (due to up-down degeneracy) was compared with its matching counterpart M along with all the 2 n-1 -1 other competitors W, to define a shape-dependent gap of selectivity S a (L), eq. ( 11). In this procedure, all bonds have the same stiffness j and the receptors are considered rigid, k = 10. This shapedependent gap of selectivity was subsequently minimized with respect to all different possible L, resulting in a quantity Gap(a, j, n) characteristic of this family of ligands. In particular, Gap(a, j, n) provides a lower bound for the occurrence of false positive if a ligand is put in contact with an equal number of all possible receptors. Using eq. ( 12) instead of eq. ( 11) leads to a different quantity Gap ′ (a, j, n) related to the noise to signal ratio when transmitting a message using a library of compounds as suggested earlier in this section.

Finally, the lower bound Gap(a, j, n) has to be optimized with respect to the internal parameters a and j, in order to assess, in a context independent manner, the intrinsic coding capacity of a family of ligands with given length n, defining in this way GAP(n). The resulting formulas are:

Gap(a, j, n) = inf ligands L of length n S a (L) ; (28)
GAP(n) = sup a,j
Gap(a, j, n)

= sup a,j inf ligands L of length n S a (L) .

The inf operator defines the selectivity bias of a matching receptor competing against a league of all other receptors. The sup operator maximizes Gap with respect to a and j.

In practice, cyclic chains were preferred to open chains, as they confer to all monomers the same importance. This optimization was performed by exhaustively scanning a rectangular grid of values (a, j) ∈ [1, 9]×[1, 9] with step 0.1, k being kept equal to 10. Our findings are rendered as a 2d plot of Gap(a, j, n = 8) vs (a, j) (Figure 19), and shows that the gap increases with both a and j, but grows only marginally beyond a, j ∼ 2. Thus, we took as representative the maximum obtained for (a = 8.8, j = 8.8) for estimating GAP(n). This arguably constitutes only a rough estimate of the sup a,j operator of eq. ( 29). Similar results were observed for all the values of n = 3, . . . 8 considered in this work. Values of ln[GAP(n)] are reported in Table 1. ln[GAP] decreases almost linearly with n (Table 1) as the number of decoys grows exponentially. Extrapolating to large n, GAP(n) vanishes for n ≥ n * ≃ 20, length for which, in the presence of an equal number of all receptors, a ligand has more chance to bind a mismatching receptor than its own complementary one.

If one uses eq. ( 12) as an alternative definition of GAP, one observes a saturation of this quantity to a n independent value, close to e 4.33 ≃ 76 (for a = 8.8, j = 8.8, k = 10, which are quite large values). This means that, irrespective of their length, any matching ligand receptor pair L -M has an affinity larger by a factor 76 than any other mismatching pair L -W. Selectivity S ′ a is not sensitive to the growing number of decoys (or "complexity") with n.

One learns from these results that if one wishes to preserve a minimal selectivity while increasing the ligand length n, one must necessarily select a subset of shapes as valid ligands and disregard the other possibilities. This is tantamount to introducing "redundancy" in the "coding scheme". Note that the numerical value of Table 1 depends on the details of the a,j maximization and on the choice of k.

Sensitivity to stiffness profiles

The following step is to investigate the importance of stiffness patterns compared with shape patterns. Is it possible to encode information in the rigidity profile ? This situation is motivated, for instance, by the bioinformatic study of Sacquin-Mora et al. who investigated correlations between mechanical properties and binding location in proteins [24]. It would be indeed very interesting to determine whether self assembling pairs possess some kind of tactile sense and are sensitive to their respective local surface rigidity.

We consider four identically shaped ligands and receptors with patterns L =u/PPppPP(c) (ligand), u/PPppPP, u/PPPPpp,u/PpPpPp and u/ppPPpp (receptors), where p stands for a soft coupling and P for a stiff coupling (Figure 20).

The corresponding affinities are shown in Figure 21. The graph shows that the best affinity is observed when both patterns coincide (self) and the lowest when the stiffness pattern is opposite to the ligand (opposite). The mismatch u/PPPPpp is almost optimal and the intermediate case u/PpPpPp gives an intermediate value. The thermodynamic study shows that in both situations, the entropic contribution is dominated by the enthalpic contribution (Figures 23, 22 and 24). The stiffness patterns act on both enthalpic and entropic terms and the issue of the competition is not simple to predict. For instance, ligand L =u/PpPpPp(c) has a stronger affinity for W =u/PPppPP than for its self pattern M (not shown).

If these results show that the affinity is sensitive to the stiffnesses, encoding information in stiffness profiles does not seem as straightforward as it is for shapes. Conveying information by means of selective binding, as sketched earlier in this section, would imply to design a quadruplet L 1 , L 2 , M 3 , M 4 of ligands and receptors with identical shapes but distinct rigidity profiles, with affinities obeying:

C a L1,M1 ∼ C a L2,M2 ; C a L1,M2 ∼ C a L2,M1 ; C a L1,M1 , C a L2,M2 ≫ C a L1,M2 , C a L2,M1 . (30)
We could not come up with such a quadruplet. These quadruplets may not exist, or require longer lengths that the ones considered in this work (say n ≥ 8). Stiffness profiles modulate the recognition process, but we were not able to prove that they could be substituted to shape profiles. We believe it is still an open issue to know for sure if stiffness profiles alone can bear some information.

Conclusion and Perspectives

It is clear that a ligand with too few degrees of freedom cannot achieve good selectivity. In the top of Figure 25, a fine receptor is in contact with a coarse ligand, which averages out the details of the receptor. If one changes one bit of the receptor, the binding properties of the ligand are only marginally altered, and the selectivity ratio stays close to 1. At the opposite, a ligand with many monomers in contact with each element of the receptor, will have a greater tolerance to the fluctuations of a particular monomer (Figure 25, bottom). Binding is enhanced by the addition of every monomer contribution. This case is a close analogous to the so-called "repetition code" in information theory. The repetition code consists in repeating many times every single bit of information to ensure the safe transmission of a code word (the message). It is a greedy procedure, as the length of the message is increased by the same factor. Here, a greater selectivity is achieved, at the expense of a greater complexity of the molecule L which is thrice as long. As a rule, repetition is the easiest way one can come up with to amplify the trends observed in Figures 9 and12. By glueing together copies of the same patterns, one automatically enhances the characteristics of the ligand-receptor repeated unit.

To summarize a comparison with the results presented in [11][12][13], we can say that we focussed mostly on interesting characteristic trends exhibited by a few selected pairs of ligand and receptors, while the authors of [11][12][13] favor global and averaged trends running over the whole set of possible patterns. Our definition of the Gap and GAP indicators differ from the free energy difference ∆F which serve as a criterion in their work. Our transfer matrix approach treats exactly the binding statistics of a given pair of ligand and receptor, without need of Monte-Carlo sampling, nor mean-field or large J approximations that are required in their 2d approach. Finally we believe that our model is the first one that makes it possible to investigate the role of local stiffness modulation, and the possibility of stiffness encoding of information.

The connection between error-correcting transmission codes and spin systems was recognized by N. Sourlas [25][26][START_REF] Nishimori | Statistical Physics of Spin Glasses and Information Processing[END_REF]. In particular, it was shown that the usual binary parity checks b 1 ⊕ b 2 ⊕ . . . ⊕ b p which involves the sum of p binary digits modulo 2 was in fact equivalent to coupling p spins s 1 s 2 . . . s p , with s i = ±1, leading to a formal connection between information theory and p-spin glasses. Our model is currently restricted to nearest neighbors coupling, and cannot account for long range couplings. Information redundancy is thus limited to simple repetition codes. It would be interesting to investigate how a second layer of spins and couplings could be added in order to better enforce robustness with respect to single bit mismatches. One also notices that a 2d generalization of the ligand shape would indeed be equivalent to a genuine 2d Edwards-Anderson spin glass. Spin glasses are well known for their long-lived or metastable states [START_REF] Mézard | Spin glass theory and beyond[END_REF]28]. Each one of these states can be put in correspondence with a matching random-field representing a different receptor.

Another interesting connection between spin systems and pattern recognition, is the Superparamagnetic clustering of data [29,30]. Wiseman, Blatt and Domanyi showed that it was possible to train a two dimensional array of Potts spins in order to recognize picture features and patterns (2d inhomogeneous distributions of points). This work could provide hints on how to train a 2d elastic network for shape and stiffness recognition.

The nonmonotonic behaviors of the affinity, or the decrease of affinity upon local stiffening of the chain are still modest, showing only a cut by half in the case illustrated in Figure (14). One may be interesting in finding stronger effects by hardening more than a single bond.

The schematic model introduced and studied here is intended to guide us towards more realistic examples, such as simple molecules that could be designed and investigated with the help of coarse grained or all atom numerical models. This, we believe, should be the next step to endeavor. 10). Sr is directly related to the probability of the ligand to bind the decoy instead of the matching pattern. As expected Sr = 1 when j = 1 because a completely soft ligand binds equally well any receptor (non specific binding).

Here a = 3. Because the stiff bond of the ligand lies on a mismatching region, its hardening leads to a decrease in the affinity. Assuming that an effector molecule is the cause of this bond stiffening, as sketched in Figure 13, one is left with a situation that reminds the proposed mechanism of Hawkins and McLeish. At the opposite, when the stiff bond lies on a matching region, its hardening leads to an increase of the affinity (not shown). 13 and14). The vertical dashed lines indicate the points where all bond stiffnesses are identical (j = j). One can infer from these curves the relative variation of configurational entropies ∆∆S caused by a change in j. For instance ∆∆S is positive (binding is favored) for j ≥ j = 1.5. The trend is reversed for j ≥ j = 3., and unsettled for j ≥ j = 2. Length

6. 1

 1 Dependence in the number of mismatches I now consider the selectivity of L =u/pppp(c), in contact with altered configurations W 1 =u/pmmp and W 2 =u/pmpm. Both altered patterns differ from the original one respectively by one and two values of the spin s i (Hamming distances d H = 1 and d H = 2, cf Figure 17). However the Hamming distance relative to the stiffness pattern is d ′ H = 2 for both W 1 and W 2 and is not directly related to the Hamming distance of configurations.

Fig. 1 .

 1 Fig. 1. Picture of the ground state of a ligand molecule L =u/ppp(o) and a matching receptor M =u/ppp. With 4 monomers, the flexible ligand explores 2 4 = 16 configuration states, while the rigid receptor does not change its conformation. On the right, a mismatching pattern W =u/pmp is shown.

Fig. 2 .

 2 Fig. 2. Affinity C a of open L =u/ppp(o) and cyclic L ′ =u/pppp(c) with their rigid matching receptor M, (Figure 1 on the left) as a function of the ligand stiffness parameter j, with a = 3. Here stiffness is shown to enhance the affinity of the pair. A ligand whose native state matches the binding site of a receptor is said to be pre-organized for matching.

Fig. 3 .Fig. 4 .

 34 Fig. 3. Thermodynamics of u/ppp(o) on its matching receptor, a = 3, expressed in units β -1 = kT

Fig. 5 .

 5 Fig. 5. Thermodynamics of L =u/ppp(o) on receptor W =u/pmp, a = 3, expressed in units β -1 = kBT .

Fig. 6 .

 6 Fig.6. Selectivity Sr between ligand L =u/ppp(o), L ′ =u/pppp(c), their matching rigid receptor M =u/ppp and a mismatching receptor W =u/pmp, as a function of the ligand stiffness j (see eq. 10). Sr is directly related to the probability of the ligand to bind the decoy instead of the matching pattern. As expected Sr = 1 when j = 1 because a completely soft ligand binds equally well any receptor (non specific binding). Here a = 3.

Fig. 7 .

 7 Fig. 7. Ground state shape of flexible cyclic ligand L =u/pppp(c) and flexible receptor M =u/pppp. The 4 monomers ligand and receptor are treated on the same footing, and the partition sum encompasses 4 4 = 256 configuration states, which is done by tracing the fourth power of a transfer matrix.

Fig. 8 .

 8 Fig. 8. Ground state shape of ligand L =u/ppp.ppp(o) and motif M =u/ppmmppm, involved in Figure (9). The dot stands for a very weak coupling constant between monomers 4 and 5.

Fig. 9 .

 9 Fig. 9. Nonmonotic behavior of the affinity C a vs j for L =u/ppp.ppp(o), M =u/ppmmppm, k = 10. The effect is more pronounced for larger values of a.

Fig. 10 .

 10 Fig. 10. Ground state shape of flexible ligand L =u/ppppppp(o), and rigid receptor W =u/pmmpmmp involved in Figure (11).

Fig. 11 .

 11 Fig. 11. Nonmonotic downhill behavior of the affinity C a vs j for L =u/ppppppp(o), M =pmmpmmp, k = 10. The effect is more pronounced when a increases.

Fig. 12 .EFig. 13 .

 1213 Fig.12. Nonmonotic behavior of the affinity C a vs j for L =u/ppppppp(o), M =u/pmmpmmp, a = 2. The effect disappears when the flexibility k of the receptor decreases and the limit case k = 1 corresponds to a very soft, patternless receptor.

Fig. 14 .

 14 Fig.14. Affinity C a as a function of the stiffness exp(βJ ) of the hard link, with j = exp(βJ) = 2 and a = exp(βA/2) = 3. Because the stiff bond of the ligand lies on a mismatching region, its hardening leads to a decrease in the affinity. Assuming that an effector molecule is the cause of this bond stiffening, as sketched in Figure13, one is left with a situation that reminds the proposed mechanism of Hawkins and McLeish. At the opposite, when the stiff bond lies on a matching region, its hardening leads to an increase of the affinity (not shown).

Fig. 15 .

 15 Fig.15. Configurational entropy change -T ∆Sc upon binding, when the stiffness of the bond linking bead 3 and 4 is increased (Figures13 and 14). The vertical dashed lines indicate the points where all bond stiffnesses are identical (j = j). One can infer from these curves the relative variation of configurational entropies ∆∆S caused by a change in j. For instance ∆∆S is positive (binding is favored) for j ≥ j = 1.5. The trend is reversed for j ≥ j = 3., and unsettled for j ≥ j = 2.

Fig. 16 .

 16 Fig.16. Energy (enthalpy) change ∆U upon binding, when the stiffness of the bond linking bead 3 and 4 is increased (Figs13 and 14). The vertical dashed lines indicate the points where all bond stiffnesses are identical (j = j). The relative energy change ∆∆U is positive (binding is unfavored) for j = 1.5 and j = 2, but slightly negative for j = 3.

Fig. 17 .

 17 Fig. 17. Ligand u/pppp(c) (top left) with a matching receptor u/pppp(top right), mismatching receptors u/pmmp (bottom left) and u/pmpm (bottom right). The Hamming distances in terms of configuration are respectively dH = 0, 1 and 2, while in terms of stiffness one has d ′ H = 0 and d ′ H = 2.

Fig. 18 .

 18 Fig.[START_REF] Mézard | Spin glass theory and beyond[END_REF]. Selectivity Sr for Hamming distance dH = 1 (straight lines) and dH = 2 (circles). The selectivity of the pair u/pppp(c)-u/pmpm is higher than the one of u/pppp(c)-u/pmmp, which means that u/pppp(c) binds u/pmmp (dH = 1) better than u/pmpm (dH = 2). Other parameters are a = exp(βA/2) = 3 and k = exp(βK) = 10.

Fig. 19 .

 19 Fig.[START_REF] Durbin | Biological sequence analysis[END_REF]. Three dimensional graph of the relative selectivity ln[Gap(a, j, 8)] as a function of a and j, for patterns with length 8. Selectivity increases slowly with a and j. Receptors are rigid, k = 10.

Fig. 20 .

 20 Fig. 20. A ligand and four receptors with same shape but different rigidity patterns are represented. The ligand is L =u/PPppPP(c), made of two rigid bonds, two soft bonds and then two rigid bonds (top). The receptors are (a) middle left u/PPppPP, (b) middle right u/PPPPpp, (c) bottom left u/PpPpPp, (d) bottom right u/ppPPpp.

Fig. 21 .

 21 Fig. 21. Affinities C a for the four pairs introduced in Figure 20 as a function of j (stiffness of the soft bonds of the ligand, symbol p,m), see text for details. Parameters are a = 2, k = 2, k = 10 (receptor hard bonds with symbol P,M), j = 10 (ligand hard bonds with symbol P,M).

Fig. 22 .

 22 Fig.[START_REF] Jones | Information and Coding Theory[END_REF]. Free energy differences in the case of Figure20and Figure21result from a competition between energy (enthalpy) and entropy, where the energy term is larger by a factor 2 compared with the entropic one.

Fig. 23 .

 23 Fig.23. Energy (enthalpy) differences in the case of Figure20and Figure21dominates the thermodynamics of binding.

Fig. 24 .

 24 Fig.24. Differences in entropy -T ∆Sc. The entropy loss is minimal for identical stiffness profile and maximal for opposite stiffness profile.

Fig. 25 .

 25 Fig. 25. A coarse molecule and a fine receptor (top), a fine molecule and a coarse receptor (bottom).

Table 1 .

 1 Table of selectivities for patterns with increasing sizes. See text for details.

		n ln(GAP)
	2	3.65
	3	3.22
	4	2.91
	5	2.67
	6	2.47
	7	2.30
	8	2.15

Acknowledgements

The author thanks Carlos Marques for discussions on this topic, and is indebted to the referees for many meaningful comments.


Simon R Eugster
email: eugster@inm.uni-stuttgart.de

Francesco Dell'isola

Roberto Fedele

Pierre Seppecher

Piola Transformations in Second-Gradient Continua

Keywords: Continuum Mechanics, Second-Gradient Continua, Piola Transformation, Principle of Virtual Work, Lagrangian Formulation, Eulerian Formulation

Second-gradient continua are defined as those continua whose internal virtual work functionals depend on the first and second-gradient of the virtual displacement. These functionals can be represented either in Lagrangian (referential) or Eulerian (spatial) description thus defining respectively the Piola-Lagrange as well as the Cauchy-Euler stress and double-stress. In this paper, we deduce the Piola transformation formulae, i.e., those relationships between all Lagrangian and Eulerian fields relevant for the formulation of the Principle of Virtual Work. In particular, we derive the Piola transformations of stress and double-stress as well as the Piola transformations for external virtual work functionals compatible with second-gradient internal work functionals. The latter transformations contain in fact the Piola transformations of the contact surface and line forces as well as the contact surface double-forces.

Introduction

The equations of motion for three-dimensional continua are naturally deduced from the Principle of Virtual Work by using the so-called Lagrangian description. The chosen space of configurations consists of the set of admissible placements, which are functions that map the position of a material particle in the reference configuration into its position in the current configuration. The reference configuration is the domain of the considered scalar, vector and tensor fields in Lagrangian description. In some applications of continuum mechanics, as for instance in fluid mechanics, it is more useful to consider fields defined on the current configuration. Also when dealing with deformable solids, externally applied forces are described more naturally as vector fields defined in the current configurations of the bodies to which they are applied. It is thus of the greatest importance to to be able to use the so-called Eulerian description too. Obviously, Lagrangian and Eulerian descriptions must be fully equivalent. Indeed, Gabrio Piola established in his fundamental works [START_REF]The complete works of Gabrio Piola: Volume I: Commented English translation[END_REF][START_REF] Dell'isola | The complete works of Gabrio Piola: Volume II: Commented English translation[END_REF] the formulas determining the relationships between Eulerian and Lagrangian descriptions. The obtained Piola transformations allow for the deduction of the equations the internal, external and inertial virtual work functionals should not change their value whether they are expressed in Lagrangian or in Eulerian form. In other words, the Piola transformation can be obtained by regarding the placement as a change of variables in the integral expression of 35 the involved virtual work functionals.

Piola transformations in classical first gradient theory are well known and widely used. Similar transformations applicable for second-gradient continua are not a trivial generalization. In this paper, we determine the Pi-40 ola transformations of stress and double-stress introduced in the Lagrangian configuration into their Eulerian counterparts. Moreover, we give the transformation formulae of the compatible Eulerian external virtual work functionals, which determine the relationships between the Eulerian and Lagrangian external line forces, surface forces and surface double-forces. We surprisingly find that Eulerian double-forces not only contribute to the expression of Lagrangian double-forces but also to Lagrangian line and surface forces. It has to be explicitly remarked that, to our knowledge, in the literature on second-gradient continua one can find either formulations in Eulerian description [START_REF] Toupin | Elastic materials with couple-stresses[END_REF][START_REF] Germain | La méthodes des puissances virtuelles en mécanique des milieux continus -1ère Partie, Théorie du second gradient[END_REF][START_REF] Germain | The method of virtual power in the mechanics of continuous media, I: Second-gradient theory[END_REF][START_REF] Smelser | La théorie du second gradient et la capillarité[END_REF][START_REF] Casal | Relation entre l'équation de l'énergie et l'équation du mouvement en théorie de Korteweg de la capillarité[END_REF] or in Lagrangian description [START_REF] Auffray | Analytical continuum mechanics à la Hamilton-Piola least action principle for second gradient continua and capillary fluids[END_REF][START_REF] Seppecher | Microscopic interpretation of strain-gradient and generalized continuum models[END_REF][START_REF] Giorgio | A study about the impact of the topological arrangement of fibers on fiber-reinforced composites: Some guidelines aiming at the development of new ultra-stiff and ultra-soft metamaterials[END_REF][START_REF] Giorgio | Lattice shells composed of two families of curved Kirchhoff rods: an archetypal example, topology optimization of a cycloidal metamaterial[END_REF][START_REF] Giorgio | In-depth gaze at the astonishing mechanical behavior of bone: A review for designing bio-inspired hierarchical metamaterials[END_REF][START_REF] Scerrato | Bio-inspired design of a 280 porous resorbable scaffold for bone reconstruction: A preliminary study[END_REF][START_REF] Placidi | A review on 2D models for the description of pantographic fabrics[END_REF][START_REF] Placidi | A strain gradient variational approach to damage: a comparison with damage gradient models and numerical results[END_REF][START_REF] Turco | Equilibrium paths of Hencky pantographic beams in a three-point bending problem[END_REF], but the relation between them seems not to be available. For instance generalization to N -th gradient continua [START_REF] Dell'isola | How contact interactions may depend on the shape of Cauchy cuts in Nth gradient con-295 tinua: approach "à la D'Alembert[END_REF][START_REF] Dell'isola | Least action and virtual work principles for the formulation of generalized continuum models[END_REF] needs only the Lagrangian description. Situations are similar when homogenization methods are used to obtain the deformation energy of materials belonging to the class of higher gradient continua, see among others [START_REF] Barchiesi | Large inplane elastic deformations of bi-pantographic fabrics: asymptotic homogenization and experimental validation[END_REF][START_REF] Barchiesi | Twodimensional continua capable of large elastic extension in two independent directions: Asymptotic homogenization, numerical simulations and experimental evidence[END_REF][START_REF] Dell'isola | Large deformations of planar extensible beams and pantographic lattices: heuristic homogenization, experimental and numerical examples of equilibrium[END_REF][START_REF] Alibert | Truss modular beams with deformation energy depending on higher displacment gradient[END_REF].

Lagrangian Kinematics

The physical space is modeled as a three-dimensional Euclidean vector space E 3 with the inner product denoted by ⟨•, •⟩ E 3 . We assume the reference configuration of the considered body B to be a subset Ω ⊂ E 3 , which is sufficiently regular to perform all the required calculations, see Figure 1 for a particular example of such a configuration. The topological boundary of Ω is denoted by ∂Ω. The boundary ∂Ω is the union of a finite number of twodimensional orientable surfaces with boundary, called faces of Ω. The faces are oriented in accordance with their corresponding outward-pointing unit normal fields N . Each of the faces' boundary curves are assumed to have a piecewise continuous tangent T as well as an outward-pointing unit normal B that is tangent to the face. The union of all boundary curves is denoted by ∂∂Ω. Remark that each curve constituting the boundaries of the faces, which are called edges of Ω, must be regarded as part of the boundary of both intersecting faces. Hence, considering Figure 1, the outward-pointing unit normal B to ∂∂Ω takes two values on the edge curve Γ. These are denoted by B + and B -depending on whether they are tangent to the face Σ + or Σ -, respectively. Edges of Ω are assumed to concur in a finite number of wedges together with a finite number of other edges. See [START_REF] Dell'isola | How contact interactions may depend on the shape of Cauchy cuts in Nth gradient con-295 tinua: approach "à la D'Alembert[END_REF] for more details about the differential-geometric assumptions accepted here.

The motion of the body B is defined as the mapping Π : Ω × R → E 3 , (X, t) → x = Π(X, t), which is a differentiable parametrization of suitably regular placement maps Π t = Π(•, t) : Ω → E 3 , which are assumed to be one to one. At some fixed time instant t, the image ω = Π t (Ω) ⊂ E 3 is called the current configuration of the body B and represents the spatial points occupied by the body B in its deformed state. We use the upper-case X ∈ Ω and the lower-case x ∈ ω to denote points in the reference and current configuration, respectively. Although working in E 3 , clever choices of representation allow for a precise classification of the appearing objects. For that reason, we use an arbitrary skew basis (G 1 , G 2 , G 3 ) to represent the referential points as linear combination X = X A G A . Note that we apply Einstein's summation convention, which implies summation over upper and lower indices that appear twice in a term. Evaluation of the inner product for all combinations of these base vectors yields the components of the associated bilinear form, which are referred to as referential metric components G AB = ⟨G A , G B ⟩ E 3 . The points in the current configuration are represented with respect to an alternative basis (g 1 , g 2 , g 3 ), which induces the current metric components g ij = ⟨g i , g j ⟩ E 3 . Henceforth, we refer to these two bases as referential and current basis and follow the convention to use upper-and lower-case letters also for the appearing indices. Hence, the placement map Π t is represented in such a way that

x = x i g i = Π t (X) = Π i t (X)g i . (1)
The components of the first gradient F = ∇Π t and the second gradient F = ∇F = ∇(∇Π t ) of the placement map are then given as

F i A := (∇Π t ) i A = ∂Π i t ∂X A , F i AB := (∇F ) i AB = ∂ 2 Π i t ∂X A ∂X B , (2)
both of which are functions of (X, t).

In order to check whether the current configuration minimizes a certain functional, the placement map Π t can be embedded in a variational family of mappings Πt : R × Ω → E 3 , (ε, X) → x = Πt (ε, X), such that Π t (X) = Πt (0, X). The variation of the placement map is then defined as

δΠ t (X) := ∂ Πi t ∂ε (0, X)g i =: δΠ i (X, t)g i (3)
inducing a virtual displacement field

δΠ : Ω × R → E 3 , (X, t) → δΠ i (X, t)g i , (4)
which geometrically can be seen as a time-dependent vector field along the placement Π t . Note that we use the notational convention to denote Lagrangian fields, which are defined on Ω, with uppercase letters. Geometrically similar fields are the velocity and acceleration fields defined as

V (X, t) := ∂Π i ∂t (X, t)g i , A(X, t) := ∂ 2 Π i ∂t 2 (X, t)g i . ( 5
)
Since the variation of the placement map induces a change in the first and the second-gradient, by the symmetry of the second derivatives, one readily sees the relations

δF i A = ∂ ∂ε ∂ Πi t ∂X A ε=0 = ∂δΠ i ∂X A , δF i AB = ∂ ∂ε ∂ 2 Πi t ∂X A ∂X B ε=0 = ∂ 2 δΠ i ∂X A ∂X B . (6)

Virtual Work Principle in Lagrangian Form

The basic assumption in a variational formulation of continuum mechanics as founded by d'Alembert and Lagrange consists in postulating that the motion of every continuum can be characterized by suitably choosing the constitutive equations for the three linear continuous virtual work functionals δW int Ω , δW ext Ω , δW dyn Ω and by assuming that, at every instant of time t, the virtual work identity

δW tot Ω (δΠ) := (δW int Ω + δW ext Ω + δW dyn Ω )(δΠ) = 0 (7)
holds for every virtual displacement field δΠ. The functional δW int Ω captures the virtual work expended in the internal interactions among parts of the considered body Ω on every virtual deformation process. The functional δW ext Ω allows for the calculation of the virtual work expended in the interactions of the body with its external world. Finally, δW dyn Ω gives the inertial virtual work expended on virtual displacements.

For second-gradient continua, the internal virtual work of Ω can always be represented as the integral of a volume density in the form1

δW int Ω (δΠ) = - Ω P A i δF i A + P AB i δF i AB , (8)
where P A i and P AB i are the components being work dual to the first and second gradient of the virtual displacement. We call them the Piola-Lagrange stress P and Piola-Lagrange double-stress P.

Let ρ 0 = ρ 0 (X) be the reference density of the continuum with mass m = Ω ρ 0 . We incorporate inertial effects by postulating the following virtual work contribution

δW dyn Ω (δΠ) = - Ω ρ 0 g ij A j δΠ i , (9)
where A j are the components of the acceleration. The external virtual work functional specifies the interactions between the considered continuum and its external world. Postulating the Principle of Virtual Work [START_REF] Casal | Relation entre l'équation de l'énergie et l'équation du mouvement en théorie de Korteweg de la capillarité[END_REF] implies that the class of possible external virtual work functionals for a continuum is limited once its internal and inertial virtual work functionals are specified. This class is called the set of compatible external virtual work functionals. To be more specific, if one assumes an external virtual work contribution that cannot be balanced either 110 by the internal or the inertial virtual work contributions, the Principle of Virtual Work tells us that these contributions must vanish. As external virtual work functionals include boundary interactions concentrated on ∂Ω, the determination of their compatibility requires a process of 115 integration by parts.

The internal virtual work [START_REF] Auffray | Analytical continuum mechanics à la Hamilton-Piola least action principle for second gradient continua and capillary fluids[END_REF] is in fact a representation of a second-order distribution in the form of (A.5). As proven in Appendix A, the internal virtual work functional can also be represented as

δW int Ω (δΠ) = Ω ∂ P A i ∂X A δΠ i - ∂Ω P A i N A δΠ i + ∂Ω M ∥ C L ∂ ∂X C P AB i N B M ∥ L A δΠ i - ∂Ω P AB i N A N B ∂δΠ i ∂X C N C - ∂∂Ω P AB i B A N B δΠ i , (10)
where M ∥ denotes the projector to the tangent spaces of ∂Ω and

P A i = P A i - ∂P AB i ∂X B . (11)
Since the inertial virtual work is provided only by a volume integral, the representation [START_REF] Giorgio | A study about the impact of the topological arrangement of fibers on fiber-reinforced composites: Some guidelines aiming at the development of new ultra-stiff and ultra-soft metamaterials[END_REF] together with the virtual work principle [START_REF] Casal | Relation entre l'équation de l'énergie et l'équation du mouvement en théorie de Korteweg de la capillarité[END_REF] implies that the external virtual work functionals must be of the form

δW ext Ω (δΠ) = Ω F Ω i δΠ i + ∂Ω F ∂Ω i δΠ i + ∂Ω D ∂Ω i ∂δΠ i ∂X C N C + ∂∂Ω F ∂∂Ω i δΠ i . (12)
In this expression, the co-vector fields F Ω , F ∂Ω and F ∂∂Ω are dual to virtual displacements and are, due to their integration domain, forces per unit reference volume, surface and line, respectively. Moreover, an additional surface 120 density field D ∂Ω appears, which, following the nomenclature used by Germain [START_REF] Germain | La méthodes des puissances virtuelles en mécanique des milieux continus -1ère Partie, Théorie du second gradient[END_REF][START_REF] Germain | The method of virtual power in the mechanics of continuous media, I: Second-gradient theory[END_REF], is called surface density of double-forces. These contributions are dual to the normal derivative of the virtual displacement. Using ( 9) -( 12) in [START_REF] Casal | Relation entre l'équation de l'énergie et l'équation du mouvement en théorie de Korteweg de la capillarité[END_REF], for an unconstrained continuum, the Principle of Virtual Work implies the equations of motion given by the partial differential equations

∂ ∂X A P A i - ∂P AB i ∂X B + F Ω i = ρ 0 g ij A j in Ω . (13)
Additionally, also the boundary conditions follow from the virtual work principle. These are on the faces

F ∂Ω i = P A i N A -M ∥ C L ∂ ∂X C P AB i N B M ∥ L A D ∂Ω i = P AB i N A N B on ∂Ω (14)
and on the edges

F ∂∂Ω i = P AB i B A N B + + P AB i B A N B -on ∂∂Ω . (15)
Recalling (A.12) from Appendix A, the symbols () ± de-125 note the limits on the curves constituting ∂∂Ω from the faces ± of the quantities in the brackets.

Virtual Work Principle in Eulerian Form

From the Lagrangian point of view considered so far, all scalar, vector and tensor valued functions depend on referential points X and time t. Since the placement Π t : Ω → ω is invertible, we can use also the Eulerian description, where the spatial points x = Π t (X) are regarded to be the independent variable. The inverse function π t = Π -1 t : ω → Ω is written with lower case letters, as it will be done henceforth for every map with ω as its domain. Moreover, we introduce the inverse map π :

ω × R → Ω × R defined as (x, t) → (X, t) = π(x, t) := (π t (x), t) = (Π -1 t (x), t) . ( 16
)
This change of parameterization induces the spatial virtual displacement δx, the spatial velocity v and the spatial acceleration a, all of which are defined via their Lagrangian counterpart as

δx := δΠ • π , v := V • π , a := A • π . (17)
According to our terminology, [START_REF] Dell'isola | How contact interactions may depend on the shape of Cauchy cuts in Nth gradient con-295 tinua: approach "à la D'Alembert[END_REF] are the Piola transformations of the virtual displacement, velocity and acceleration. The components of the first and second-gradient of the spatial virtual displacement are abbreviated as

δd i j := ∂δx i ∂x j , δd i jk := ∂δx i ∂x j ∂x k . ( 18
)
Lagrangian and Eulerian formulations describe the same physical phenomena. For this reason, the value of the virtual work expressions should not be affected by a change of variables. Therefore

δW int ω (δx) := δW int Ω (δΠ) , δW ext ω (δx) := δW ext Ω (δΠ) , δW dyn ω (δx) := δW dyn Ω (δΠ) . (19)
Consequently, the virtual work equality holds also in the Eulerian description:

δW tot ω (δx) := (δW int ω + δW ext ω + δW dyn ω )(δx) = 0 . ( 20
)
Since after a simple change of variables, the internal virtual work integral expression is still a representation of a second-order distribution, it can be written in the form

δW int ω (δx) = - ω c j i δd i j + c jk i δd i jk , (21)
where c j i and c jk i are the components of the work conjugates to the first and second gradient of the spatial virtual 130 displacement. We call them the Cauchy-Euler stress c and the Cauchy-Euler double-stress c.

The inertial virtual work is expected to be of the form

δW dyn ω (δx) = - ω ρg ij a j δx i , (22)
where ρ = ρ(x, t) denotes the current mass density of the continuum resulting in the current inertial force density

f dyn,ω i = -ρg ij a j .
135

To characterize the compatible external virtual work functional, the same integration by parts procedure can be applied as in the Lagrangian formulation. Defining

cj i = c j i - ∂c jk i ∂x k , (23)
in Appendix A it is proven that the Eulerian internal virtual work functional has the following representation:

δW int ω (δx) = ω ∂c j i ∂x j δx i - ∂ω cj i n j δx i + ∂ω m ∥ c l ∂ ∂x c c jk i n k m ∥ l j δx i - ∂ω c jk i n j n k ∂δx i ∂x c n c - ∂∂ω c jk i b j n k δx i . (24)
Note that n and b denote the outward-pointing unit normals to the faces ∂ω and edges ∂∂ω of the current configuration and m ∥ the projector to the tangent spaces of ∂ω. Consequently, the compatible external virtual work functionals must be of the form

δW ext ω (δx) = ω f ω i δx i + ∂ω f ∂ω i δx i + ∂ω d ∂ω i ∂δx i ∂x c n c + ∂∂ω f ∂∂ω i δx i , (25)
where the co-vector fields f ω , f ∂ω and f ∂∂ω are forces per unit current volume, surface and line, respectively. Also in the Eulerian framework, there appears a surface density of double-forces d ∂ω i , which is a density per unit current surface and which is dual to the normal gradient with respect 140 to the current normal vector.

The Principle of Virtual Work then implies for an unconstrained continuum the governing equations of motion ∂ ∂x j c j i -

∂c jk i ∂x k + f ω i = ρg ij a j in ω , (26)
together with the boundary conditions on the faces

f ∂ω i = cj i n j -m ∥ c l ∂ ∂x c c jk i n k m ∥ l j d ∂ω i = c jk i n j n k on ∂ω (27)
and on the edges

f ∂∂ω i = c jk i b j n k + + c jk i b j n k -on ∂∂ω . (28)

Piola Transformations of Mass Density, Stress and Double-Stress

In the previous section, we have introduced as Eulerian dual quantities the Cauchy-Euler stresses together with the current inertial force, the current external forces and current double-forces. The Piola transformation problem raises the following question: "Which are the relationships between the Lagrangian and Eulerian stresses and double-stresses, inertial forces as well as external forces and double-forces such that the identities ( 19) are satisfied?".

We present in this section the Piola transformation between the referential and current mass density as well as the Piola transformations between the Piola-Lagrange and the Cauchy-Euler stresses and double-stresses.

Let Φ be a time-dependent Lagrangian field with domain Ω that is related to the corresponding Eulerian field ϕ with domain ω by Φ(X, t) = ϕ(Π t (X), t). Recalling (2), the chain rule implies that the gradients of the Lagrangian and Eulerian fields are connected by

∂Φ ∂X A (X, t) = ∂ϕ ∂x j (Π t (X), t)F j A (X, t) . ( 29
)
Since this relation can also be written as

∂Φ ∂X A (π(x, t)) = ∂ϕ ∂x j (x, t)F j A (π(x, t)) , (30)
we will drop the arguments in what follows. The indication of the arguments may be complemented by the reader. Using this convention together with (2), we obtain, by taking once more the gradient of (29), the expression

∂ 2 Φ ∂X A ∂X B = ∂ 2 ϕ ∂x j ∂x k F j A F k B + ∂ϕ ∂x j F j AB . (31)
Consequently, the gradients of the Lagrangian and Eulerian virtual displacement fields are related by

δF i A = δd i j F j A , δF i AB = δd i jk F j A F k B + δd i j F j AB . (32)
In a similar way, the time derivative of the Lagrangian field Φ together with (5) leads to

∂Φ ∂t (X, t) = ∂ϕ ∂x j (Π t (X), t)V j (X, t) + ∂ϕ ∂t (Π t (X), t) (33)
As done in (30), by using V j (π(x, t)) = v j (x, t), this relation can be written in spatial coordinates as

∂Φ ∂t (π(x, t)) = ∂ϕ ∂x j (x, t)v j (x, t) + ∂ϕ ∂t (x, t) , (34)
which then serves as definition of the material time derivative d dt ϕ(x, t) := ∂ϕ ∂x j (x, t)v j (x, t)

+ ∂ϕ ∂t (x, t) . ( 35
)
This is why the spatial acceleration field a can be expressed as the material time derivative of the spatial velocity field, which in formulas reads as

a i (x, t) = A i (π(x, t)) = ∂V i ∂t (π(x, t)) = d dt v i (x, t) . ( 36
)
To obtain the transformation formulae between referential and spatial fields, we will take the Lagrangian virtual work expressions ( 8) and ( 9) and change the integration variable from the referential coordinates X to the spatial 160 coordinates x resulting in an integral expression defined on the domain ω.

Transforming the right hand side of ( 8), by the change of variables given by π t (x), whose Jacobian determinant is equal to J -1 (π(x, t)) = det(F (π(x, t))) -1 , we get with the above introduced abuse of notation

δW int Ω (δΠ) = - ω J -1 P A i δF i A + P AB i δF i AB ( 37
)
The integration domain ω makes it clear that all capitalized Lagrangian fields must be composed with π. Next, we insert the relations (32) and compare the outcome with the Eulerian virtual work expression [START_REF] Dell'isola | Large deformations of planar extensible beams and pantographic lattices: heuristic homogenization, experimental and numerical examples of equilibrium[END_REF] in the following way

δW int Ω (δΠ) = - ω J -1 P A i δd i j F j A - ω J -1 P AB i (δd i jk F j A F k B + δd i j F j AB ) = - ω J -1 P A i F j A + P AB i F j AB δd i j - ω J -1 P AB i F j A F k B δd i jk = - ω c j i δd i j + c jk i δd i jk = δW int ω (δx) . (38)
Comparison of the last and second last line leads to the Piola transformation of stress and double-stress

c j i = J -1 P A i F j A + P AB i F j AB , c jk i = J -1 P AB i F j A F k B , (39)
which can be written in a more synthetic way as

2 c = J -1 P • F t + P : F t , c = J -1 P : (F t ⊗F t ) . (40)
Note that with our notation, either the left hand sides have to be composed with Π or the right hand sides with π. The reader will remark that, by simply assuming P = 0, we get

165 the Piola transformation formula c = (J -1 P • F t ) • π for first gradient continua.
Transforming in the same way as before the Lagrangian inertial virtual work expression, we get

δW dyn Ω (δΠ) = - ω ρ 0 g ij A j δΠ i J -1 . (41)
Using the definition of the spatial virtual displacement and acceleration from [START_REF] Dell'isola | How contact interactions may depend on the shape of Cauchy cuts in Nth gradient con-295 tinua: approach "à la D'Alembert[END_REF], we can further manipulate the expression to

δW dyn Ω (δΠ) = - ω ρ 0 J -1 g ij a j δx i = - ω ρg ij a j δx i = δW dyn ω (δx) , (42)
where we recognize in the last step the relation between referential and current mass density. This Piola transformation of the density expressed in Lagrangian coordinates is generally referred to as referential equation of continuity and reads as

ρ 0 (X) = J(X, t)ρ(Π(X, t), t) . (43)
The spatial equation of continuity

d dt ρ(x, t) + ρ(x, t) ∂v i ∂x i = 0 (44)
follows immediately by carrying out the following computations. First, take the time derivative of (43), then insert the identity ∂J/∂t = J∂v i /∂x i and, last, express the result in spatial coordinates with the spatial velocity field and utilize the material time derivative (35).

Piola Transformations of external forces and double-forces

The final relations left to show are the Piola transformations of the forces and double-forces appearing in the Lagrangian and Eulerian external virtual work functionals [START_REF] Giorgio | In-depth gaze at the astonishing mechanical behavior of bone: A review for designing bio-inspired hierarchical metamaterials[END_REF] and (25). We proceed in a similar way as just done for the stresses and double-stresses. However, some further unexpected difficulties arise. Carrying out the change of variable given by the inverse of the placement, the Eulerian normal derivative to ∂ω of the virtual displacement does not transform into the Lagrangian normal derivative to ∂Ω. As a consequence, the Eulerian virtual work functional of double-forces, once transformed into Lagrangian coordinates followed by a Lagrangian surface integration by parts, is composed of three virtual work functionals: a first one to be identified with the Lagrangian virtual work functional of double-forces, a second one contributing to Lagrangian surface forces and the last one contributing to Lagrangian edge forces.

Using the inverse of the right Cauchy-Green strain C -1 with the components

(C -1 ) AB = (F -1 ) A i g ij (F -1 ) B j , (45)
as well as

K A = M ∥ A C (C -1 ) CB N B = (F -1 ) A i g ij (F -1 ) B j -∥F -T • N ∥ 2 G AB N B , (46)
we get after long but straightforward calculations the following transformation formulae for the external forces:

F Ω i = Jf ω i , F ∂Ω i = ∥JF -T • N ∥f ∂ω i -M ∥ B A ∂ ∂X B (Jd ∂ω i K A ) , F ∂∂Ω i = ∥F • T ∥f ∂∂ω i + J(C -1 ) AB B A N B d ∂ω i + + J(C -1 ) AB B A N B d ∂ω i -. (47)
Using definition (A.3), we can write the Piola transformations of the external forces in direct notation as

F Ω = Jf ω , F ∂Ω = ∥JF -T • N ∥f ∂ω -div ∂Ω ∥ Jd ∂ω ⊗ K , (48)
F ∂∂Ω = ∥F • T ∥f ∂∂ω + J(B • C -1 • N )d ∂ω + + J(B • C -1 • N )d ∂ω -. (49)
Finally, the Piola transformation of the surface doubleforce, once in index notation and once in direct notation, reads as

D ∂Ω i = J∥F -T • N ∥ 2 d ∂ω i , D ∂Ω = J∥F -T • N ∥ 2 d ∂ω . ( 50
)
Appendix A. Equivalent form for second-order distributions

Let us consider a regular manifold V embedded in an ndimensional Euclidean vector space E n and the idempotent projectors m ∥ and m ⊥ on its tangent and normal spaces. When V has co-dimension one, and if n denotes its unit normal vector, we have

m ⊥ γ α = n γ n α , m ∥ γ α = δ γ α -n γ n α . (A.1)
Given a vector field w defined in the neighborhood of V, the divergence theorem for Riemannian submanifolds with boundaries is stated as

V m ∥ γ α ∂ ∂x γ m ∥ α β w β = ∂V m ∥ α β w β m ∥ γ α b γ = ∂V w β b β , (A.2)
where ∂V denotes the boundary of V and where the unit vector b is tangent to V and normal to ∂V. Defining div V ∥ by setting for all smooth fields ϕ

div V ∥ (ϕ) α := m ∥ γ α ∂ ∂x γ (ϕ) , (A.3)
the divergence theorem reads as

V div V ∥ (m ∥ • w) = ∂V w • b . (A.4)
In accordance with the theory of distribution, both the virtual work expression in Lagrangian and Eulerian description can be considered as distributions D represented in the form

D(ϕ) = v s α ∂ϕ ∂y α + s αβ ∂ 2 ϕ ∂y α ∂y β , (A.5)
where the derivatives of the test functions ϕ are taken with respect to the coordinates y α of a 3-dimensional Euclidean space. Note that the index i appearing in both ( 8) and ( 21) does not play any role in the present considerations and is therefore omitted. Since the order of a distribution can be defined as the sum of the highest order derivative of the test function and the co-dimension of the integration domain in the domain of the test functions, D is a second-order distribution. The symbol v denotes the generic integration domain which satisfies the same continuity requirements as discussed in Section 2 for the reference configuration Ω. The faces of the subset v are denoted by ∂v and come along with the outward-pointing unit normal field n. The symbol ∂∂v denotes the edges on which the outward-pointing unit normals b to the boundaries ∂∂v are defined. Moreover, these normals are in the tangent planes to the faces constituting ∂v.

Using the product rule in the second integrand of (A.5), we can write

D(ϕ) = v s α - ∂s αβ ∂y β ∂ϕ ∂y α + v ∂ ∂y β s αβ ∂ϕ ∂y α .
With the abbreviation sα = s α -∂s αβ /∂y β and applying the product rule for the first integrand, we end up with

D(ϕ) = v ∂ ∂y α (s α ϕ) - v ∂s α ∂y α ϕ + v ∂ ∂y β s αβ ∂ϕ ∂y α .
Using the divergence theorem for the first and the third term and introducing the zeroth-and first-order distributions

D 0 v (ϕ) = - v ∂s α ∂y α ϕ , D 0 ∂v (ϕ) = ∂v sα n α ϕ , (A.6)
the distribution (A.5) can be written in the form

D(ϕ) = D 0 v (ϕ) + D 0 ∂v (ϕ) + ∂v s αβ ∂ϕ ∂y α n β . (A.7)
The last term here, is the only expression in which still derivatives of ϕ appear. Therefore, we will manipulate this term further by using the projectors (A.1) for the faces ∂v The distribution D I ∂v is a second-order transverse distribution involving the normal derivative of the test function (∂ϕ/∂y γ )n γ and cannot be reduced any further.

Due to the idempotence of the projector m ∥ and by applying once more the product rule, we can manipulate the first term in the last line of (A.8) in the following way Note, to obtain D 0 ∂∂v , the divergence theorem (A.2) has been applied leading to a line integral along the edges of v. We used here a notational convention in the expression of the integral. As depicted in Figure 1, we observe that an edge γ is the intersection of two subsurfaces σ + and σ -, say. Hence, γ is traversed twice: once with the surface normal n -, edge normal b -and the limit (s -) αβ approached from the surface σ -, as well as once with the corresponding n + , b + and (s + ) αβ . Consequently, if we denote each of the edge curves by γ i for i = 1, . . . , n e , then the integral expression of the first equality in (A.11) reads as 3 This equivalence can be applied to the Lagrangian internal virtual work functional [START_REF] Auffray | Analytical continuum mechanics à la Hamilton-Piola least action principle for second gradient continua and capillary fluids[END_REF], when carrying out the -P A i , s αβ → -P AB i , sα → -P A i , ϕ → δΠ i , y α → X A , v → Ω, b → B and n → N with the corresponding projector m ∥ → M ∥ . Remark the now appearing implicit summation over the index i, which has no influence on the derivations presented above. As a consequence of (A.13), the Lagrangian internal virtual work functional can then also be represented by [START_REF] Giorgio | A study about the impact of the topological arrangement of fibers on fiber-reinforced composites: Some guidelines aiming at the development of new ultra-stiff and ultra-soft metamaterials[END_REF]. The following replacements are required in (A.6), (A.9) and (A.11) to obtain for the Eulerian internal virtual work functional the representation (24): s α → -c j i , s αβ → -c jk i , sα → -c j i , ϕ → δx i , y α → x j as well as v → ω.

Figure 1 :

 1 Figure 1: Reference and current configuration of a body with a piecewise continuous boundary surface and piecewise continuous boundary curves.

∂v s αβ ∂ϕ ∂y α n β = ∂v s αβ ∂ϕ ∂y γ n β δ γ α = ∂v s αβ ∂ϕ

  ∂y γ n β m ∥ γ α + m ⊥ γ α = ∂v s αβ ∂ϕ ∂y γ n β m ∥ γ α + D I ∂v (ϕ) , (A.8)where we have introduced the distributionD I ∂v (ϕ) = ∂v s αβ n α n β ∂ϕ ∂y γ n γ .(A.9)

  γ s αβ n β m ∥ λ α ϕ = D 0 ∂∂v (ϕ) + D0 ∂v (ϕ) .(A.10)In the last step, we have introduced the second-and firstorder distributionsD 0 ∂∂v (ϕ) = ∂∂v s αβ n β b α ϕ ,

s

  αβ n β b α + + s αβ n β b α -ϕ . (A.12)In conclusion, from the point of view of the theory of distributions, the second-order distribution D(ϕ) from (A.5) can equivalently be represented asD = D 0 v + (D 0 ∂v + D0 ∂v ) + D I ∂v + D 0 ∂∂v . (A.13)

Without fixing a particular constitutive law for P and P, the sign of the virtual work expression is irrelevant.

Note the non-standard operators required for the higher order tensors. Using the vectors a, b ∈ E

, the co-vectors α, β ∈ (E * ) 3 and the third order tensor A ∈ ⊗ 1 2 E 3 , these operators can be defined by the equalities (α ⊗ a)⊗(β ⊗ b) = (α ⊗ β) ⊗ (a ⊗ b) and α • A : (a ⊗ b) = (a ⊗ b) : A t • α. The remaining operators are standard and can be recovered by comparing (40) with (39).

following replacements in (A.6), (A.9) and (A.11): s α → 3 It is also common to introduce b± = t×n ± as the vector product between tangent vector to the curve γ and the surface normal n, see[START_REF] Toupin | Elastic materials with couple-stresses[END_REF][START_REF] Germain | The method of virtual power in the mechanics of continuous media, I: Second-gradient theory[END_REF] 24]. In that case b -= -band the sign changes in the second term of (A.12).

Acknowledgment. The paper was published with the financial support of the Ministry of Higher Education and Science of the Russian Federation as part of the program of the Moscow Center for Fundamental and Applied Mathematics under the agreement No. 075-15-2019-1621.