A predictive model for protein materials: from macromolecules to macroscopic fibers

TL;DR

This paper introduces a predictive model for the mechanical behavior of protein materials, specifically keratin fibers and spider silks, using limited experimental parameters to accurately simulate their cyclic and damage behaviors.

Contribution

It presents a novel macroscopic model derived from molecular parameters that captures damage, residual strains, and dissipation in protein-based materials.

Findings

Model accurately reproduces cyclic behavior of keratin fibers

Capable of describing dissipation and permanent strains in spider silks

Robust model with high quantitative accuracy

Abstract

We propose a model for the mechanical behavior of protein materials. Based on a limited number of experimental macromolecular parameters (persistence and contour lengths, rate of unfolding dissipation) we obtain the macroscopic behavior of keratin fibers (human, cow, and rabbit hair), taking into account the damage and residual stretches effects which are fundamental in many functions of life. We support our theoretical results by showing that our model is robust and able to reproduce with high quantitive accuracy the cyclic experimental behavior of different keratinous protein materials we tested. We also show the capability of describing, even if with lower precision, the dissipation and permanent strain effects in spider silks.

Figures

Equations

$${\varphi }_{tot}={\varphi }_{tot}(L,L_c)=\kappa \frac{L^2}{L_c-L}+q{L}_c,$$

$$f=\frac{\partial {\varphi }_{tot}(L,L_c)}{\partial L}=\kappa \frac{2L{L}_c-L^2}{(L_c-L{)}^2}.$$

$$L_{un}\equiv L_{max}=\frac{L_c}{\sqrt{\kappa /q}+1},f_{un}=\left(2\sqrt{\kappa /q}+1\right)q.$$

$$L_n=\sqrt{\bar{n}}b=\sqrt{\frac{L_c}{b}}b=\sqrt{L_{c}b},$$

$$\left\{\begin{array}[]{ll}\lambda=\frac{L}{L_{o}},&\mbox{\text@underline{total stretch}},\vspace{0.05 cm}\\ \lambda_{e}=\frac{L}{L_{n}},&\mbox{\text@underline{elastic stretch}},\vspace{0.05 cm}\\ \lambda_{n}=\frac{L_{n}}{L_{o}}&\mbox{\text@underline{permanent stretch}},\vspace{0.05 cm}\\ \lambda_{c}=\frac{L_{c}}{L_{o}}&\mbox{\text@underline{limit extensibility stretch}}.\end{array}\right.$$

$$\begin{array}[]{l}\sigma\!=\!N^{fib}\kappa\frac{\frac{2\lambda}{\bar{\lambda}_{c}(\lambda_{max})}-\left(\frac{\lambda}{\bar{\lambda}_{c}(\lambda_{max})}\right)^{2}}{\left(1-\frac{\lambda}{\bar{\lambda}_{c}(\lambda_{max})}\right)^{2}}\!+\!\mu^{net}\!\left(\lambda\!-\!\frac{1}{\lambda^{2}}\right)+\frac{p}{\lambda\sqrt{\lambda}},\vspace{0.2 cm}\end{array}$$

$${\lambda }_c={\bar{\lambda }}_c({\lambda }_{max})=\left(\sqrt{\frac{\kappa }{\bar{q}}}+1\right){\lambda }_{max}$$

$$\begin{array}[]{ccc}\includegraphics[scale={0.56}]{HUMANHAIR.pdf}&\includegraphics[scale={0.50}]{RABBIT.pdf}&\includegraphics[scale={0.50}]{COW.pdf}\vspace{0 cm}\\ \mbox{a)}&\mbox{b)}&\mbox{c)}\\ \end{array}$$

$$\begin{array}[]{lll}\phi_{tr}&=&\displaystyle-\sum_{i=1}^{n}(Q-J)(1-\chi_{i})-J\sum_{i=1}^{n-1}(1-\chi_{i})(1-\chi_{i+1})\\ &=&\displaystyle-\sum_{i=1}^{n}Q(1-\chi_{i})+Jn_{b_{f}},\end{array}$$

$$\begin{array}[]{l}\phi_{tot}=-k_{B}T\ln[p(L,n_{u},n_{b_{f}})]=k_{B}T\Omega(n_{u},n_{b_{f}})\exp(-\frac{\phi_{e}(L,n_{u})}{k_{B}T})\exp(-\frac{\phi_{tr}(n_{u})}{k_{B}T})\vspace{0.2 cm}\\ =\phi_{e}(L,n_{u})+\phi_{tr}(n_{u})-TS(n_{u},n_{b_{f}}),\end{array}$$

$${\varphi }_{tot}={\varphi }_e(L,n_u)+n_{u}Q+\text{mbox}const.$$

$${\phi }_e={\phi }_e(l,l_c)=\kappa \frac{l^2}{l_c-l},$$

$${\varphi }_{\mathrm{tot}}=\kappa \frac{L^2}{L_c-L}+n_{u}Q.$$

$${\varphi }_{tot}={\varphi }_{tot}(L,L_c)=\kappa \frac{L^2}{L_c-L}+q{L}_c.$$

$$f=\frac{\partial {\varphi }_{tot}(L,L_c)}{\partial L}=\kappa \frac{2L{L}_c-L^2}{(L_c-L{)}^2},$$

$$g=-\frac{\partial {\varphi }_{tot}(L,L_c)}{\partial L_c}=\kappa \frac{L^2}{(L_c-L{)}^2}-q.$$

$$L_{un}\equiv L_{max}=\frac{L_c}{\sqrt{k/q}+1},$$

$$f_{un}=\left(2\sqrt{k/q}+1\right)q.$$

$$L_n=\sqrt{\bar{n}}b=\sqrt{\frac{L_c}{b}}b=\sqrt{L_{c}b},$$

$$L_o=\sqrt{L_c^{o}b}.$$

$$\left\{\begin{array}[]{ll}\lambda=\frac{L}{L_{o}},&\mbox{\text@underline{total stretch}},\vspace{0.05 cm}\\ \lambda_{e}=\frac{L}{L_{n}},&\mbox{\text@underline{elastic stretch}},\vspace{0.05 cm}\\ \lambda_{n}=\frac{L_{n}}{L_{o}}&\mbox{\text@underline{permanent stretch}},\vspace{0.05 cm}\\ \lambda_{c}=\frac{L_{c}}{L_{o}}&\mbox{\text@underline{limit extensibility stretch}}.\end{array}\right.$$

$$f=\kappa \frac{\frac{2\lambda }{{\bar{\lambda }}_c({\lambda }_{max})}-{\left(\frac{\lambda }{{\bar{\lambda }}_c({\lambda }_{max})}\right)}^2}{{\left(1-\frac{\lambda }{{\bar{\lambda }}_c({\lambda }_{max})}\right)}^2}\text{vspace}0.2cm$$

$${\lambda }_c={\bar{\lambda }}_c({\lambda }_{max})=\left(\sqrt{\frac{\kappa }{\bar{q}}}+1\right){\lambda }_{max},$$

Full text

G. Puglisi1, D. De Tommasi1, M.F. Pantano2, N. Pugno3, G. Saccomandi4.

1 Dipartimento di Scienze dell’Ingegneria Civile e dell’Architettura, Politecnico di Bari, Italy

2 Laboratory of Bio-Inspired and Graphene Nanomechanics, Department of Civil, Environmental and Mechanical Engineering, University of Trento, Via Mesiano 77, 38123 Trento, Italy

3 Laboratory of Bio-Inspired and Graphene Nanomechanics, Department of Civil, Environmental and Mechanical Engineering, University of Trento, Via Mesiano 77, 38123 Trento, Italy;

Agenzia Spaziale Italiana, Via del Politecnico, 00133 Roma (Italy);

School of Eng. and Mat. Science, Queen Mary University of London, Mile End Road, London E1 4NS, U.K.

4 Dipartimento di Ingegneria, Università degli Studi di Perugia.

School of Mathematics, Statistics and Applied Mathematics, NUI Galway, University Road, Galway, Ireland

1 Introduction

Experimental analyses (AFM, optical and magnetic tweezers, nanoindentation) clearly show that the outstanding elasticity, toughness, strength, and self-healing properties [15, 2] of protein materials originate from their secondary structure, characterized by the presence of folded (crystal-like) domains, typically in the form of $\alpha$-helices or $\beta$-sheets, which can undergo unfolding as a consequence of an applied displacement. The efficacy of the unfolding mechanism is based on the presence of a large number of non-covalent forces [3], typically hydrogen-bonds, which act as sacrificial joints. Being much weaker than the covalent peptide bonds, these interactions can easily break, causing unfolding phenomena ($\alpha\rightarrow\beta$ transition or $\beta$ domains unfolding). Such unfolding prevents the backbone fracture conciliating the typically conflicting requests for high stiffness and toughness [4]. Indeed, the cooperative strength of non-covalent bonds provides proteins with stability and stiffness before unfolding begins, whereas the increase of the end-to-end length during unfolding allows to accommodate large stretches (20$\%$-40$\%$), resulting in high energy dissipation and toughness.

Here we focus on two different classes of structural protein materials [5] known for their outstanding performance: keratinous materials, which can be found in wool, hair, cells intermediate filaments, epithelial cells and hooves, and silk, a valuable material produced by spiders and silkworms. In the first case the secondary structure is in the form of a $\alpha$-helices undergoing a transition to an almost straight polypeptide chain in the form of $\beta$-sheets [6, 7]; whereas in the second case the protein secondary structure is in the form of (alanine rich) $\beta$-sheet folded domains which undergo unfolding under stretching [8]. Both silk and keratinous materials have attracted significant attention for the design of new bioinspired materials [2, 7, 4, 13, 14, 9, 10, 11], such as nano-polymers [9] and block copolymers [10] in the case of spider silk and hydrogels [11] in the case of keratinous materials. Protein unfolding occurring in the above listed proteins, but also in others, like cytoplasmic IF proteins [15], needs a deep comprehension as it is linked to many biological phenomena (i.e. mechano-transduction and cell motility [16]) and human deceases, such as Alzheimer, type II Diabetes, prion diseases, Parkinson [17].

In this context, as in the case of multiscale dislocation and defect theories in metal plasticity [18], a significant advance of knowledge is strictly related to the availability of models relating the material properties at the macromolecular scale to the macroscopic response of wires and tissues [8, 19, 20, 21, 7, 2]. However, modeling complex phenomena, such as hard domains stability under mechanical loading, different behavior under thermal and mechanical loading [22], asymmetric unfolding/refolding behavior under cyclic loading, still represents a fundamental open task. Molecular Dynamics approaches are limited by the involved time scales [24, 23], whereas the possibility of formally connecting the observed macroscopic rate-independent dissipation with the micro rate-dependent dissipation [25] can be provided by considering the evolution in the multi-valley energy landscapes with three different time scales of loading, thermal vibration, and in-well relaxation [28, 16, 27].

Here, we follow a phenomenological approach and, as in [21] and [12], we consider for the single molecule under assigned elongation a phenomenogical Griffith type approach assuming that the macromolecule unfolds when the energy gain (evaluated as difference of the two force-elongation curves between which the macromolecule ‘jumps’) equals a material parameter, representing the energy ‘dissipated’ in the mechanical unfolding of the crystal. This approach is supported by the experimental observation that the macroscopic toughness ($\sim 10^{2}$MJ/m3) of the considered materials corresponds to the energy of the H-bonds of the unfolded fraction during elongation [30].

The contour length variations associated to the unfolding effect can be similarly evaluated by the macromolecular stretching experiment [11]. Then, based on a (slightly amended) affinity hypothesis [3] and the James and Guth three-chains model [2], following [34] we deduce the macroscopic behavior with damage and residual stretch of protein materials and its analytical relation with few key macromolecular material parameters that can be deduced by AFM experiments.

To show the ability of the model of reproducing the experimental behavior we cyclically tested three different keratin materials and spider dragline silk.

2 Micro-macro model

Macromolecule energy. Following a Griffith type total energy minimization approach [12], as in [6] and [35], we model a protein macromolecule as a lattice of $n$ two-states (folded/unfolded) domains (Fig.1a). We define the state of the $i$-th domain through the variable $\chi_{i}$ with $\chi_{i}=0$ $(=1)$ in the folded (unfolded) state. For simplicity, since we are interested in the unfolding regime, we neglect both the elastic energy and the length of the folded domains. Moreover, we neglect non-local interactions (weak interaction hypothesis). The total internal energy can then be expressed as $\phi_{e}=\sum_{i=1}^{n}\chi_{i}\varphi_{e}(\lambda_{rel}^{i})$, where $\varphi_{e}$ depends on the relative deformation $\lambda^{i}_{rel}=l_{i}/l_{c}$ of the $i$-th unfolded domain, with end-to-end length $l_{i}$, contour length (assumed constant) $l_{c}$ and persistent length $L_{p}$. We also assume the limit extensibility condition: $\lim_{\lambda_{rel}\rightarrow 1}\varphi_{e}(\lambda_{rel})=+\infty$.

Equilibrium requires a constant force $l_{c}\frac{d\varphi_{e}(l_{i}/l_{c})}{dl_{i}}=F$ and, under the hypothesis of convex $\varphi_{e}$, a constant $\lambda_{rel}^{i}=\lambda_{rel}=l_{i}/l_{c}=L/L_{c}$. Here, $L=\sum_{i}\chi_{i}l_{i}$ is the total end-to-end length and $L_{c}=\sum_{i}\chi_{i}l_{c}$ is the total contour length of the unfolded fraction. As a result we have $\phi_{e}=L_{c}\varphi_{e}(\lambda_{\mbox{\tiny{\it rel}}})=L_{c}\varphi_{e}(L/L_{c}).$ Finally, by adopting the elastic energy density proposed in [12] we obtain $\phi_{e}=\phi_{e}(L,L_{c})=\kappa\frac{L^{2}}{L_{c}-L}$, where $\kappa=k_{B}T/(4L_{p})$.

Let us now consider the configurational energy of the different folded/unfolded states. Following [12] consider an Ising type transition energy and a di-block approximation (as supported by MD simulations [13], see Supplementary Material, SM, for details), we end up with a total energy $\phi_{tot}=\kappa\frac{L^{2}}{L_{c}-L}+\frac{L_{c}}{l_{c}}Q$ depending on the single (history dependent) internal variable $L_{c}$, where $Q$ is the unfolding energy of the single hard domain.

In Fig.1b we show the scheme of our minimization procedure, requiring minimization of $\phi_{tot}$ in order to get the transition unfolding length $L_{u}$ (see Fig. 1c and SM for details) corresponding to an increasing number $n_{u}$ of unfolded domains. The system stays in a given configuration (fixed folded fraction) until the elastic energy gain equals the unfolding energy $Q$. This energy is then assumed as ‘dissipated’ by heat. The resulting force-end to end length behavior, shown in Fig.1c, reproduces the experimental behavior in [11] of the single macromolecule.

Continuum macromolecule approximation. Based on the observation of a large number of folded domains (see SM), after introducing the dissipation density $q=Q/l_{c}$, we may consider a continuum limit approximation which allows the chain total energy to be rewritten as

[TABLE]

where now $L_{c}$ can be considered a continuum variable. The equilibrium force in the continuum case is

[TABLE]

Under our Griffith fracture hypothesis of global energy minimization ($\partial\phi_{tot}/\partial L_{c}=0$) we obtain

[TABLE]

As a result the protein unfolding behavior is characterized by a constant stress unfolding plateau, depending only on the persistence length $L_{p}$ and rate of dissipation $q$ as represented in Fig.1d (path OABCD).

Remark. Dashed lines in Fig. 1d) are unloading and reloading paths. The introduction of partial refolding upon unloading would let us to describe also internal hysteric effects. Similarly, it would be possible to introduce rate of loading dependent effects [37] upon considering a different macromolecular behavior [38].

Affinity hypothesis. As classical in polymer elasticity [3], the main hypothesis for the deduction of the macroscopic behavior is the affinity hypothesis which assumes that the macroscopic stretches coincide with the macromolecular ones. We here give an interpretation of the permanent deformations observed in protein materials as a macroscopic counterpart of the variation of the macromolecule natural configuration associated with the unfolding of crystal domains and the availability of new monomers (hidden length). We then extend the affinity hypothesis by identifying both permanent and elastic components of the macromolecular stretch with the corresponding macroscopic counterparts. In such a way we are also able to describe the important permanent stretch effect.

To analytically describe this phenomenon we begin by observing that for a given contour length $L_{c}$ the natural length (zero force) $L_{n}$ of the entropic chain can be expressed, according to a known result of Statistical Mechanics (see *e.g.*​​ [3]), as

[TABLE]

where $b$ is the length of the Kuhn segments and $\bar{n}$ is the number of Kuhn segments of the unfolded chain fraction. Thus, if we denote by $L_{c}^{o}$ the initial value of the contour length, the initial natural end-to-end length is $L_{o}=\sqrt{L_{c}^{o}\,b}.$ Then we may define the following stretch measures (see SM):

[TABLE]

Macrosopic material behavior. Structural proteins are composite materials of hard segments immersed in an amorphous soft fraction of unfolded macromolecules. To describe the complex inter-chains interactions and self avoiding effects, following the additivity hypothesis in polymer elasticity [1], we determine the energy of the real network as the sum of the energy of ideally isolated fibrils (folded molecules in Fig.2b)) whose behavior is schematized in Fig.2a)) plus a network energy term measuring the chains interactions: red springs in Fig.2b). Moreover, according with the classical James and Guth three-chains model, we assume that all chains are aligned along the three principal macroscopic stretch directions. In particular, we suppose that the unfolding protein macromolecules are oriented along the fiber direction, whereas the amorphous soft fraction is supposed to be equally distributed along the three principal directions. For the unfolding chains we assume the continuous limit energy (10), force (11), unfolding forces (13), and stretches (14) deduced above. The network effect is instead modeled by simple Gaussian chains.

Thus, let $(\lambda_{1},\lambda_{2},\lambda_{3})$ be the macroscopic principal stretches, where $\lambda_{1}$ is the stretch in the fiber direction. According to previous hypotheses, the macroscopic energy density can be calculated as $\Phi(\lambda_{i})=N_{fib}\phi_{tot}(\lambda_{1})+\sum_{i=1}^{3}\frac{N_{net}}{3}\phi_{net}(\lambda_{i}).$ Here $N_{fib}$ is the number of macromolecules with unfolding domain, per unit area of fiber section in the reference configuration, and $N_{net}$ is the number of chains, per unit volume, reproducing the real chains network effect. The network energy is modeled as Gaussian $\phi_{net}(\lambda_{i})=\frac{k_{B}T}{2}(\lambda_{i}^{2}-1)$ leading ([3]) to a neo-Hookean network energy $\Phi_{net}=\frac{\mu}{2}(I-3),$ where $I=\sum_{i=1}^{3}\lambda_{i}^{2}$ is the first invariant of the left Cauchy-Green deformation tensor and $\mu=k_{B}TN_{net}/3\,$ represents the shear elastic modulus. On the other hand the energy of the unfolding chains gives a history dependent term depending only on $\lambda_{1}$ stretch. In such a way we describe a fundamental effect with the damage localized along the (maximum elongation) fibre direction. It is worth noticing that this damage anisotropy is crucial for an effective derivation of the macroscopic permanent stretch.

Suppose now that the fiber (assumed incompressible) undergoes a simple extension $\lambda_{1}=\lambda$ and $\lambda_{2}=\lambda_{3}=1/\sqrt{\lambda}$. To take into account the confinement effect we introduce a constant pressure $p$ perpendicular to the fiber skin. Thus, if $\sigma=\sigma_{1}$ is the (Piola/engineering) principal stress in the fiber direction, by using (11) and (14) and by imposing the boundary conditions $\sigma_{2}=\sigma_{3}=p$, we get the stress-stretch relation

[TABLE]

where $\lambda_{max}=\lambda\mbox{ if }\dot{\lambda}_{max}>0,\mbox{ \it{primary loading} }$ and $\lambda_{max}=\mbox{const}\mbox{ if }\dot{\lambda}_{max}=0\mbox{ \it{unloading/reloading}. }$ Here $\lambda_{max}$ is the previously maximum attained strain and, by using (13),

[TABLE]

is the (limit) value of the stretch corresponding to the present value of the contour length.

3 Experimental validation

To prove the efficacy of our model we performed cyclic tensile tests on two different types of material: keratin hair materials and dragline silk.

We first consider keratinous materials, such as human, cow, and rabbit hairs. The tests and the possibility of reproducing the macroscopic behavior of these keratin fibers are shown in Fig. 3.

The macromolecular parameters are coherent with the values reported in the literature. The values of the persistent length match those reported in [7] where $l_{p}=0.4$ nm. The values of the rate of dissipation $q=45\,k_{B}T$/nm agree with the value in [15] from which one deduces that the energy for the mechanical unfolding of a single amino acid is about $100$ kcal/mol, which corresponds to $q=59.6\,k_{B}T$/nm, by assuming an amino acid length of about 1 nm. A higher value $q=180\,k_{B}T$/nm is adopted for the rabbit hair and this is due to the much higher experimental unfolding stress (about 600 MPa) that we here interpret as a particularly high rate of dissipation $q$ corresponding to a higher energy barrier regulating the transition between the folded and unfolded configuration. The number $N=10^{18}$ of macromolecules per unit area corresponds to a coiled-coil cross section of 1 nm as suggested in [7].

It is important to point out that the only ‘macroscopic’ parameters of the model are the pressure $p$ accounting for the cuticle action [40] and the shear modulus $\mu$ of the amorphous fraction. To the knowledge of the authors no value for the confinement pressure $p$ is available, whereas different methods like torsional and traverse extension and torsion tests can be adopted to evaluate $\mu$ (the values considered here are compatible with the results reported in [41]). Despite the simplification of the model it quantitatively well reproduces the stiffnesses, the residual stretches, the damage softening, and the sudden softening corresponding to the reconnection to the primary loading curve (Return Point Memory) of the fiber whose damage and residual stretches are regulated by the previous maximum value of the stretch.

In Fig. 4 we reported the cyclic behavior of a spider silk produced by a cellar spider (Pholcus phalangioides in the inset of the figure). For sample preparation, the spider was left to fall down toward the ground and during its falling it produced its dragline silk, which we collected and glued on a paper frame, resulting in a fiber with 1 cm length and a diameter of about 1 $\mu$m.

The outstanding mechanical properties observed in the macroscopic material behavior is the result of a complex multiscale organization and a five layers skin-core structure [42]. While a complete description of such behavior would ask to take care of this complex structure that extends to the microscale in a hierarchical architecture [43], here we focus, based on the introduced model, on the main source of dissipation, induced by the hard-soft transition [35] of the crystalline fraction, in form of alanine rich $\beta$-sheets, undergoing unfolding under increasing strain [8, 30]. The skin has a confinement effect inducing a prestretch in the aligned macromolecules rich of $\beta$-sheets belonging to the core and counterbalanced by the skin pressure [30]. Comparison of supercontracted vs native silks in [30] determined this pressure to be about $500$ GPa. Specifically in this case the skin undergoes plastic deformation [42] that, following again [30] we model as a constant pressure $p=400$ MPa in the deformed configuration (this corresponds to substituting the term $\frac{p}{\lambda\sqrt{\lambda}}$ in (5) with $\frac{p}{\lambda}$). The matrix contribution is again modeled through a Neo-Hookean material.

Also in this case the values of the macromolecules characteristic parameters are taken from the literature. In fact we considered the values $l_{p}=0.8$ nm and $l_{c}=25$ nm estimated in [30, 44] and [8], the unfolding energy $Q=625\,\,k_{B}T$ agrees with the values estimated in [12]; the number of fibers per unit area $N^{fib}=4.75\times 10^{18}$ (per square meters) reflects the results in [45, 46].

While the model neglects important effects such as the observed presence of crystal domains in the matrix [30, 47], the observation that the behavior at small stretches is typically regulated by both the amorphous fraction and the tertiary structure [11, 43], the importance of hierarchical structure and energy exchange among different scales, Fig.4 shows that it can capture at least qualitatively, but with less quantitative precision, the fundamental phenomena of damage and residual stretches. In particular we may observe the ability of estimating the dissipated energy, that for the considered spider is about 75 MJ/m3 and it is in agreement with the value reported e.g. in [30] ($130-190$ MJ/m${}^{3})$.

Our results show that it is possible to model with high precision the complex mechanical behavior of keratineous protein materials starting from few known molecular properties and, even if with less quantitative precision, also the experimental behavior of spider silk. This represents in our opinion an important step forward in the comprehension of the link between the properties of protein secondary structure and the macroscopic material properties (toughness, stiffness and large deformations prior to fracture, residual stretches), which can be of great importance for better understanding protein materials and for designing new smart bioinspired materials.

Acknowledgments. NMP is supported by the European Research Council (ERC PoC 2015 SILKENE n. 693670, by the European Commission under the Graphene Flagship (WP14 ‘Polymer Composites’, no. 696656) and by the EU’s Horizon 2020 research and innovation programme (NEUROFIBRES, no.732344). DT and GP are supported by funding from Italian Ministry MIUR-PRIN voce COAN 5.50.16.01 code 2015JW9NJT.