#  Hierarchical physically based machine learning in material science: the
case study of spider silk

Vincenzo Fazio1, Nicola Maria Pugno1,2∗,
Orazio Giustolisi3, Giuseppe Puglisi3∗∗ 1 Laboratory for Bioinspired, Bionic,
Nano, Meta Materials & Mechanics,
University of Trento, Via Mesiano 77, 38123 Trento, Italy;  2 School of
Engineering and Materials Science, Queen Mary University of London,
Mile End Road, London E1 4NS, U.K.;  3 Department of Civil Environmental Land
Building Engineering and Chemistry, Polytechnic University of Bari, via
Orabona 4, 70125 Bari, Italy<EMAIL_ADDRESS><EMAIL_ADDRESS>
###### Abstract

Multiscale phenomena are characterized by complex structure-function
relationships emerging where entities at different scales aggregate into
structures with unexpected final properties. Mathematical modelling for
multiscale phenomena typically requires the deduction of a set of differential
equations at different scales that predict the macroscopic behavior. The
complexity of these equations and the number of essential parameters make
developing effective, predictive models challenging. To address this
challenge, recent scientific literature has explored the possibility of taking
advantage of the availability of sophisticated numerical techniques in the
artificial intelligence and machine learning fields.

Here, we focus on a fundamental aspect in multiscale phenomena, i.e the
recognition of the hierarchical role of variables. In the spirit of the
considered numerical approaches and especially in determining efficient tools
to deduce effective analytical relations for material modelling, we focus, in
a Pareto front interpretation, on the determination of simple accurate
relations, starting from experimental multiscale analyses. From a physical
point of view, the aim is to deduce information at higher scales from lower
scales data, possibly respecting their hierarchical order. A crucial aspect of
the proposed approach is the deduction of causality relations among the
different variables to be compared with the available theoretical notions and
possibly new interpretations resulting by the data modelling. This result in a
stepwise approximation going from data modelling to theoretical equations and
back to data modelling.

To demonstrate the key advantages of our multiscale numerical approach,
compared to classical, non-physically based data modelling techniques, we
consider the explicit example of spider silk –a material with exceptional
properties hugely analyzed also in the spirit of bioinspiration and strictly
based on an evolutionary hierarchical optimization–. The description of the
complex behavior of such material and the importance of the structure
organization at different scales represent an open problem in material science
and within the field of the design of new bioinspired materials. The
remarkable macroscopic spider silk behavior is indeed the result of
interesting mesostructures arising from the aggregation of amino acids at the
molecular scale. The comparison with recent data modelling results, neglecting
causality and the multiscale character of the examined material, demonstrates
the importance of the search for new data modelling approaches aiming at a
growth of a more deep scientific knowledge in the field.

We argue that, due to the generality of our results, our approach may
represent a proof of concept in many fields where multiscale, hierarchical
differential equations regulate the observed phenomenon.

###### keywords:

Multiscale modelling, Data modelling, Materials Science, Spider silk,
Evolutionary Polynomial Regression approaches.

††journal: Journal

##  1 Introduction

Multiscale models play a crucial role in different fields of theoretical and
applied science, especially due to the increasing possibility of experimental
analyses and technologies working down to the nanoscale such as Atomic Force
Microscopy (AFM), optical tweezers, magnetic tweezers [1], etc. The
technological impact has been incredible, with the possibility of designing
systems at very low scale NEMS, MEMS [2, 3] also with the production of new
smart systems based on a hierarchical organization at different scales:
multiscale metamaterials [4], multiscale gecko inspired adhesive systems [5],
self-cleaning lotus inspired surfaces [6, 7], etc. As a matter of fact, in
different fields a huge experimental literature delivering big data libraries
on hierarchical systems, starting from the nano and micro scales, up to the
macro scale, is now available. The scientific impact of such experimental
results in new numerical/theoretical tools delivering correct physical
interpretation in several fields of impact in engineering [4], medicine [8],
physiology [9], biology [10] and physics [11].

Indeed, the analysis of the now available huge mole of experimental data can
lead to a corresponding increase in the theoretical understanding and
modelling of the resulting physical system, only if adequate numerical
instruments of data modelling are available. We live in the digital age and
the possibility of new instruments such as unprecedented power of calculation
and machine learning techniques opened up the possibility of new incredible
tools of analysis of these multiscale data. On the other hand, as in every
transition, the digital transition brings significant risks and drawbacks if
not deeply analyzed in its possible effects. Thus, machine learning can lead
to a scientific knowledge growth or obscuration, rationalization or
unclearness, access to deeper theoretical models or reliance on purely data
mining approaches.

In this work, we trace a rational way in the direction of deducing new tools
for the modelling of multiscale phenomena based on machine learning techniques
that may lead to real advances in scientific knowledge. Among many different
techniques with the potentiality of significant growth in this perspective,
the proposed approach started from the following observations on the state of
the art in the field. First, we point out that among data-driven techniques,
many of which have been developed in recent years, Artificial Neural Networks
(ANN) and Genetic Programming (GP) are the most adopted to model complex, non-
linear processes including multiscale hierarchical phenomena. Loosely
speaking, ANN uses models consisting of multiple processing elements (neurons)
connected by links of variable weights (parameters) to deduce typically
‘black-box’ representations of the analyzed systems. Learning in ANN involves
adjusting the parameters (weights) of interconnections in a highly
parametrized system.

In few words, the main widely recognized disadvantages of ANN model
construction are the curse of dimensionality, overfitting issues and parameter
estimation [12, 13]. The well-known curse of dimensionality refers to the
exponential increase in the need of parameters when the model input space
grows. This means that the number of connections exponentially raises and in a
such widened space the training set of input becomes more sparse or the amount
of data needed to preserve a constant level of accuracy increases
exponentially. On the other hand, in such a way, ANN acquires greater
flexibility in mapping events with complex structure. However, this leads
often to overfitting problems, that is ANN tends to fit training data too
precisely due to the large number of parameters resulting in the propensity to
generate poor predictions for events not close to the training data set.
Another disadvantage of using ANN is the difficulty of incorporating knowledge
derived from known physical laws into the learning process.

Despite these drawbacks, several significant results in this field have been
reported. Here we recall the very interesting recent fundamental result of the
use of machine learning techniques to predict protein topological
conformations from amino acid sequences with high accuracy [14]. In the field
of multiscale materials modelling, we may recall that Gu et al. [15] have
employed finite elements analysis together with convolutional neural network
algorithms to predict and optimize the toughness of hierarchical composite
systems and validated their results through additive manufacturing and
testing. Recently in [16] the authors adopted ANN to chose among a class of
specific constitutive models depending on the right Cauchy Green deformation
tensor invariants, the model that best reproduces stress-strain behaviors
under different classes of deformation. While the approach is interesting, it
is highly oriented by the specific knowledge of the problem and restricted to
the special case when the class of constitutive laws is already known: i.e.
the stress dependence on the deformation invariants.  As further machine
learning application on materials science, a k-means clustering approach was
employed to predict the behavior of heterogeneous materials under irreversible
processes like inelastic deformations obtaining a data-driven, two-scale model
[17]. Also, Neural Network approaches were employed within the field of the
meta-materials to perform tasks like topological optimization [18], while
Bayesian machine learning was employed in a data-driven design of metamaterial
building blocks based on several design variables [19]. Good predictive
performances were obtained also by neural network methods in linking the
elastic properties of composite materials to their mesoscale structure, in
particular, the three-dimensional microstructure to its effective
(homogenized) properties [20]. For steels, machine learning based
microstructural analysis, property prediction, and properties-to-
microstructure inverse analysis were conducted [21].

The implications are thus fundamental and let us obtain relevant information
for problems that have longly been theoretically unresolved, such as the
recalled long-lasting problem of predicting the protein structures from amino
acid sequences [14]. On the other hand, the main drawback in the perspective
of extending the knowledge for the theoretical modelling of such phenomena is
that ANN leads to “black-box” approaches. There is then a strong limitation on
the ‘operational’ advantages due to the lack of interpretability of the
artificial intelligence results. Some very recent works address this issue
[22, 23], but this is still an open problem [24, 25] due to the intrinsic
nature of the approach, summarized above.

On the other hand, GP is an evolutionary computing method that can generate a
more ‘transparent’ representation of the system. In particular, in this field,
the symbolic regression is a technique proposed by Koza [26] that creates
mathematical expressions to fit a set of data points using the evolutionary
process of the GP. In brief, like all evolutionary computing techniques,
symbolic regression manipulates populations of solutions (in this case
mathematical expressions) using operations analogous to the evolutionary
processes observed in nature. The genetic programming procedure mimics natural
selection through successive generations of solutions improving the fitting of
the data points. As a result, GP allows global exploration of expressions
providing insights into the explicit relations between input and output data.

Here we explore the potential synergies between machine learning and
multiscale modelling to produce robust predictive models that take into
account the underlying physics to handle challenging problems in this field.
In particular, the main purpose is to show that the multiscale character
typically corresponds to a hierarchical organization and such hierarchy should
be crucial in the process of data mining and data modelling, based on the
underlying physical phenomena. The present article, therefore, contributes to
the definition of a physically based data modelling where the expression
“physically based” here wants to highlight that although machine learning has
been effectively used to automate data processing and ensure high accuracy and
repeatability of outcomes of many physical phenomena, the most diffused
techniques typically ignore basic physical principles, which may result in non
effective approaches possibly leading to unphysical results.

More in detail, here we consider the application of a well known symbolic
data-modelling method named Evolutionary Polynomial Regression (EPR) [27, 28],
which enables the discovery of explicit and generalizable equations for the
underlying physical model. Indeed, EPR could be defined as a “gray box”
approach that, unlike typical ANN methods, is structured in such a way that it
is possible to integrate information arising from data modelling with
established physical knowledge of the phenomena. Specifically, an explicit
advantage is the possibility of deducing the independent variables of the
model at each scale and substantial properties of the functional dependencies
among the involved variables. Interestingly, as we explicitly show in this
paper, the use of EPR let us obtain a number of different analytical models,
optimized using a multi-objective method. The choice among these models is
then based on their physical interpretation and the effectiveness of the
model. The new obtained relations can then be used as input for the successive
numerical analysis, thus allowing for a continuous interchange between
physical interpretation and data modelling. Indeed, EPR generates symbolic and
explicit nonlinear equations that include a small number of polynomial
structural parameters.

To analyze the efficiency of the proposed approach in treating complex
multiscale hierarchical phenomena, we here consider the field of constitutive
modelling of complex material behaviors. Many biological examples of
evolutionary material optimization represent the possibility of obtaining
unreached material performances at the macroscale, based on a clever,
hierarchical organization of weak composing materials at the lower scales
[29]. A further enrichment in biological structure is to possibly include
different composing materials [30]. The analytical description of how the
macroscopic performances result from these complex low scale material
organizations is far from being reached and represents an meaningful benchmark
not only for their theoretical interest, but also in the crucial field of
bioinspired material design [31, 32].

Specifically, we here consider the paradigmatic example of spider silk, one of
the most studied natural materials due to its extreme mechanical properties,
particularly its strength and toughness. The availability of increasingly
sophisticated experimental techniques allowed for a deeper understanding of
its complex multiscale, hierarchical material structure. Despite this, many
relevant phenomena governing the strong material history dependence, rate,
temperature, and humidity effects remain unknown. To be specific, our data
modelling analysis is based on recent experimental observations on a large
number of silks from different spider species from all over the world, where
several material properties at different involved scales have been cataloged
for the fist time in a comprehensive database[33].

We refer to previous data modelling results that, while allowing qualitative
interpretation of the results, cannot directly afford new analytic models. As
a result, the availability of such new comprehensive experimental data
represents for sure a fundamental possibility of gaining new insight in this
multiscale model. This requires not only to refer to purely statistical
properties of the available data, as in classical data modelling [33], but to
adopt numerical techniques allowing for the deduction of analytical relations
to be theoretically interpreted and implemented.

This work aims to be general within the framework of a multiscale description
of physical phenomena and the deduction of larger scales properties from the
structures at lower scales. Indeed in the formulation of the specific case
study here analyzed, we have considered three scales starting from the micro
(protein) scale, to the macro scale passing through the meso scale. We
explicitly impose in our approach that these three scales interact with each
other in a hierarchical way. In particular, we consider the three
possibilities of deduction of the meso from the micro properties, a successive
macro from meso and eventually an interesting direct micro to macro deduction.

This is then, in our opinion, a first step toward a more effective adoption of
the new availability of data and data modelling techniques that can be of
fundamental help in several fields of multiscale phenomena when dealing with a
high number of experimental data.

##  2 Evolutionary Polynomial Regression

In this section we give, for the help of the reader, a brief introduction to
the mathematical treatment of numerical optimization problems based on EPR
algorithms. We refer to [27, 28] for a detailed description of the method.

EPR method generates explicit mathematical expressions to fit a set of data
points starting from the symbolic equation

${\bf Y}=\sum_{j=1}^{m}f(\mathbf{X},g(\mathbf{X}),a_{j})+a_{0}$ (1)

where ${\bf Y}$ is the vector of output dependent variables, $\mathbf{X}$ is
the vector of input variables, $f$ is a polynomial function composed of $m$
terms generated by the algorithm, linearly depending on the unknown $a_{j}$
parameters, $g$ is a function defined by the user (in our case we will
consider power laws with each input variable ${\bf X}_{i}$ raised to an
exponent varying among an a priory fixed set of numbers), plus the bias term
$a_{o}$.

Thus in the case of power expressions considered in the following, EPR
technique generates formulae of variable number of polynomial terms,
performing a global search of the expression of symbolic expressions for $f$.
Synthetically, EPR is performed in two steps: a) structure identification and
b) parameter estimation. The first stage entails simultaneously determining
the best ‘arrangement’ of the independent variables and the related exponents.
A multi-objective genetic algorithm termed OPTIMOGA, which stands for
Optimized Multi-Objective Genetic Algorithm, is used to finalize this
optimization. This algorithm has three targets: maximization of the so called
fitness function (in other words a measure of how closely the regression
expression fits the data points), minimization of the number of polynomial
coefficients, and reduction of the number of inputs. Observe that, since the
user defines a priori the set of candidate exponents, the possible negligible
input variables are obtained by including zero as a candidate exponent. This
allows for the fundamental aspect, recalled in the introduction, of
determining the effective independent variables. The values of the parameters
$a_{j}$ are determined in a second stage using the linear Least Squares (LS)
approach, which minimizes the Sum of Squared Errors (SSE). In addition to the
usual LS search, the LS is performed by searching for only positive values
(constrains $a_{j}>0$) to avoiding overfitting, by excluding sequences of
terms with negative/positive $a_{j}$ values that may result from the modelling
of the data noise [34]. Moreover, the uncertainty of the coefficients
($a_{j}$) is evaluated during the search and the distribution of estimated
pseudo-polynomial coefficients is used to eliminate those parameters whose
value is not sufficiently larger than zero [35, 27]. Indeed, it may be argued
that a low coefficient value with respect to the variance of estimates relates
to terms that describe noise rather than the underlying function of the
phenomenon being studied.

As a starting point the candidate independent variables, the general
polynomial structure, functions, exponents, and the maximum number of terms
are assigned based on the starting knowledge of the physical phenomenon. The
exponents can reflect the types of relationships between the inputs and
output. For example, if the vector of candidate exponents is chosen to be
${\bf EX}=[-1,-0.5,0,0.5,1]$, the maximum number of terms is $m=4$ and if the
candidate independent input variables are $k=3$, the polynomial regression
problem is to find a matrix of exponent ${\bf ES}_{4\times 3}$. In a first
stage, an initial population of matrix of exponents is generated. An example
of such a matrix is

${\bf ES}_{4\times 3}=\begin{bmatrix}1&0.5&0\\\ 0&0&1\\\ 0&-0.5&1\\\
-1&0&0.5\end{bmatrix}$ (2)

so that the expression (1) is given as:

${\bf Y}=a_{o}+a_{1}\,{\bf X}_{1}\,{\bf X}_{2}^{0.5}+a_{2}\,{\bf
X}_{3}+a_{3}\,{\bf X}_{2}^{-0.5}\,{\bf X}_{3}+a_{4}\,{\bf X}_{1}^{-1}\,{\bf
X}_{3}^{0.5}$ (3)

The adjustable parameters $a_{j}$ are then computed by minimizing the SSE as a
cost function. It follows the evaluation of the fitness function: if the
termination criterion is satisfied, the output results are shown, otherwise a
new matrix of exponents is generated through Genetic Algorithm (GA) including
crossover, mutation, ranking selection [27] and again the adjustable
parameters are calculated and the fitness function evaluated until the
termination criterion is satisfied.

Interestingly, for a given data set of observations, a regression-based
technique needs to search among an infinite number of possible models to
explain those data. Among different models that equally perform for the
description of a given phenomenon the simplest one is chosen (Occam’s razor or
principle of parsimony). The approach then considers both accuracy and
simplicity, that is it searches for an effective model describing the data as
well as a simple expression i. e. with a small number of inputs and polynomial
coefficients. The small number of constants to be estimated helps to avoid
overfitting problems, especially for small datasets. The equally performing
models are those composing the Pareto dominance front [36, 28] and since EPR
returns the whole set of formulae of the Pareto front, the final choice of the
model among different possible relations is then based on physical
considerations [37]. The selection of a symbolic model through the reading of
a Pareto front of similar but increasingly complex models brings into play the
knowledge of the phenomenon and the insight of the user as opposed to a mere
decision on the acceptance of the single model generated by classical
statistical regression frameworks.

In our opinion, this is a fundamental aspect because lets us understand the
effective functional dependence among the different involved variables. This
is why EPR can be categorized as a grey-box approach because it displays the
link between inputs and outputs with an explicit expression, whose consistency
can be easily analyzed and understood based on physical insight. More
explicitly it allows the possibility of analyzing physical causality among the
variables. In this respect, we observe that GP generates formulae/models for
$f$, coded in tree structures of variable size, performing a global search of
the expression for $f$ as symbolic relationships among ${\bf X}$ while the
parameters $a_{j}$ play a role only in the optimization process. On the other
hand, ANN goal is to map $f$, without focusing on the level of knowledge of
the functional relationships among ${\bf X}$. This is why we argue that EPR
represents a better tool for data-driven knowledge discovery.

##  3 Case study: spider silk

Spider silk is one of the most studied natural materials due to its extreme
mechanical properties, particularly its strength and toughness, which
outperform many high-performance man-made materials. Furthermore, spider silks
are regarded as the foundation of a new class of high-performance fibers in
the context of biomimetics [38, 31]. The availability of increasingly
sophisticated experimental techniques has allowed for a deeper understanding –
both chemically and structurally – of the complex multiscale, hierarchical
material structure at the heart of their notable mechanical behavior over the
last few decades. Despite this, many relevant phenomena governing their
loading history dependence, rate, temperature, and humidity effects remain
unknown, particularly when multiscale effects are taken into account [39].

At the molecular level, spider silks are made up of an amorphous matrix of
oligopeptide chains and pseudo-crystalline regions composed primarily of
polyalanine $\beta$-sheets [40, 41] with dimensions ranging from $1$ to $10$
nm [42], mostly oriented in the fiber direction [43]. The radial cross section
of the fiber is highly organized [44, 45, 41]. Furthermore, the chemical and
structural composition varies according to the different silks produced by the
different glands and, of course, the different species. Here, we focus on the
most performing and extensively studied type of silk known as dragline.

###  3.1 Micro scale

Spider dragline silk fibers (also known as Major Ampullate silk) are
constituted by structural proteins called Spidroins, which are divided into
two major subtypes, MaSp1 and MaSp2. The overall sequence architectures of the
two subtypes are similar, with a highly repetitive core region flanked by
small N-terminal and C-terminal domains (NTD and CTD, respectively). The
repetitive regions, which account for 90$\%$ of the primary structure, are
composed of alternating runs of polyalanine and multiple glycine-rich motifs
arrayed in tandem. Moreover, very recent studies, prompted primarily by
advances in proteomics and sequencing technologies, paint a more complex
picture of dragline silk composition than a simple MaSp1/MaSp2 dichotomy [33].
However, despite the complexity of the composition of the spider silks, here
we only consider the protein MaSp1 and MaSp2 which are widely recognized as
the two main composing the spider silk. From the secondary structure point of
view, the MaSp1 is mainly organized into pseudo-crystalline polyalanine
$\beta$-pleated sheets [44, 46]. On the other hand, the MaSp2 is mainly
constituted by proteins with a proline content preventing the formation of
$\beta$-sheet crystals [41] resulting in a structure with significantly lower
crystallinity and macromolecules with weaker crystal domains, typically in the
form of $\alpha$-helix and $\beta$-turns [41, 47].

We remark that as recognized in polymer mechanics [48] and described also for
the spider silk case in [49], the number of monomers of the macromolecule
(i.e. protein for the silk case) is fundamental for the mechanical behavior of
the material. Based on the fact that (i) the mechanical behavior of the spider
silk material is to be ascribed to the repetitive region features more than
the terminal region of the protein [50], and (ii) the pseudo-crystalline
$\beta$-sheets, mainly present in the MaSp1, are recognized to be the most
impactful feature in determining the exceptional strength of the spider silk
[51], here we consider the following three quantities describing the protein
scale of the silk material:

  * •

length of the repetitive region of the protein MaSp1 in terms of number of
amino acids

  * •

length of the repetitive region of the protein MaSp2 in terms of number of
amino acids

  * •

length of the polyalanine $\beta$-sheet in the protein MaSp1 in terms of
number of alanine amino acid

###  3.2 Meso scale

At the meso scale we consider the proteins’ secondary structure, how
macromolecules are arranged in the fiber and properties regarding the chemical
and structural stability of the polymer. In particular, we analyze the
following material properties:

  * •

Birefringence. It reflects the degree of molecular orientation of silk protein
chains. The birefringence of the dragline silk fiber was calculated from the
retardation value and silk fiber diameter [33].

  * •

Degree of crystallinity. It was calculated based on wide-angle x-ray
scattering (WAXS) analysis [33]. In particular, it was obtained as the ratio
of the total area of the separated crystalline scattering components to that
of the crystalline and amorphous scattering components as resulting from the
1D profile obtained by the two-dimensional (2D) diffraction.

  * •

Degradation temperatures. This quantity gives a measure in the chemical and
structural stability of the silk. In [33] the thermal degradation temperature
has been defined as the temperature that yielded 1% weight losses in the silk
samples. Indeed, heating leads to changes of the molecular weight that in turn
decreases the mass due to the production of gaseous byproducts of the chemical
reactions.

###  3.3 Macro scale

Spider silk is a very interesting material from the point of view of its
mechanical performance at the macroscopic scale. In particular, here we focus
on the material stiffness and strength. The stiffness, of the order of tens of
GPa, is above man-made polymers and at the top among other natural materials.
The strength is even more interesting, being comparable with high strength
steels (order of magnitude of 1 GPa) and with the most performing man-made
composites like the carbon and kevlar reinforced composites [52]. The reason
for these so outstanding properties with respect to standard materials, is not
yet clear, with an relevant role also in the extremely small diameter of
dragline spider silk [53]. For this reason, we also consider the diameter in
the properties at the macro scale. Finally, we address the very significant
role of hydration in the material behavior of spider silks. Indeed, a striking
effect observed in spider silks is the so called supercontraction effect,
addressed, to the knowledge of the author, for the first time in 1977 [54],
that occurs when a spider silk thread is exposed to humidity. Depending on the
silk composition, the experiments show the existence of a Relative Humidity
(RH) threshold beyond which the fiber contracts up to half of its initial
(dry) length. This also results in the possibility of exploiting the
supercontraction in the actuation field [55]. The experimentally observed
contraction depends on several factors, including spider species [56], type of
silk (among the up seven different ones that some spiders can produce [57,
58]), environmental conditions [59] and hydration rate [60]. The quantities we
consider at the macro scale are therefore the following:

  * •

Young’s modulus, obtained from the stress-strain curves determined through
tensile tests of single dragline silk fibers conducted at 25∘C and RH$\approx
50\%$ [33].

  * •

Tensile strength, calculated as the breaking force determined by tensile test
divided by the undeformed cross-sectional areas of the fiber samples
determined by SEM observations [33].

  * •

Diameter, determined by SEM observations [33].

  * •

Maximum supercontraction, calculated as $(L_{0}–L_{f})/L_{0}$, where $L_{0}$
is the length in dry condition and $L_{f}$ in fully wet conditions (RH=100%)
[33].

##  4 Modelling strategy

For all the EPR run the maximum number of terms was set to $3$ to avoid
overfitting and allow the aggregation of the inputs. The chosen set of
candidate exponents is $[-1,-0.5,0,0.5,1]$, adopted to keep the expressions
simple. Moreover, the expressions were optimized with a bias term $a_{o}$
since this element may compensate for the possible lack of relevant inputs in
the model.

As training data, we employed the experimental results recently published in
[33]111Each quantity is considered with the unit of measurement reported in
the original database, namely GPa for Young’s modulus and limit stress, $\mu
m$ for the diameter, ∘C for the thermal degradation temperature, and number of
amino acids for all the micro scale properties, whereas the supercontraction
and the crystallinity are two nondimensionional quantities ranging in (0,1).
We remark that the information regarding the dimension of each term of the
expression found by the EPR algorithm are included within the parameters
$a_{j}$ (see Eqn. (1)) estimated by means the minimization of the SSE.. The
above described quantities are summarized in Table 1 with the corresponding
adopted symbol (the type of font distinguishes the scales):

Micro scale | $a$ | Length of the repetitive region of MaSp1
---|---|---
| $b$ | Length of the repetitive region of MaSp2
| $c$ | Length of the polyalanine $\beta$-sheet in the MaSp1
Meso scale | $A$ | Crystallinity
| $B$ | Birefringence
| $C$ | Thermal degradation temperature (1% loss)
Macro scale | $\mathbb{A}$ | Young’s modulus
| $\mathbb{B}$ | Tensile strength
| $\mathbb{C}$ | Diameter
| $\mathbb{D}$ | Supercontraction
Table 1:  Material properties considered for the data modelling case study
divided by scales.

As a main parameter of accuracy, we report for the different numerical results
the Coefficient of Determination222Here we consider the classical definition
$R^{2}=1-\sum_{i=1}^{N}\frac{(x_{i}^{num}-x_{i}^{exp})^{2}}{(x_{i}^{exp}-\bar{x}^{exp})^{2}}$,
where the $x_{i}^{num}$ are the output variables of the numerical test and
$x_{i}^{exp}$ are the corresponding experimental values, with $i=1,...N$,
where $N$ is the number of experimental observations considered as dependent
variables.. On the other hand, EPR also considers other, not explicitly
reported here, indicators of performance, e.g. the sum of squared errors
(SSE). As a result $R^{2}$ does not necessarily increase as the complexity of
the expressions grows. The physical valence of the expressions found is
discussed, by following [33], through the comparison with the correlation
matrix represented in Fig. 1 obtained by calculating the Pearson correlation
coefficient333It is a measure of linear correlation between two sets of data
$\\{x_{i},i=1,...,n\\}$ and $\\{y_{i},i=1,...,n\\}$ with $n$ the number of
data, defined as
$\rho=\frac{\sum_{i=1}^{n}(x_{i}-\bar{x})(y_{i}-\bar{y})}{\sqrt{\sum_{i=1}^{n}(x_{i}-\bar{x})^{2}}\sqrt{\sum_{i=1}^{n}(y_{i}-\bar{y})^{2}}}$,
where $\bar{x}=\frac{1}{n}\sum_{i=1}^{n}x_{i}$ and
$\bar{y}=\frac{1}{n}\sum_{i=1}^{n}y_{i}$ are the mean values. for all spider
silks. Observe that from the database [33] we considered only the data where
the searched output and the considered input are reported simultaneously.
Thus, since there are some experimental properties missing for some silks of
the database, the number of silks composing the training set is different for
each considered output. This approach allows us to consider for each target
output the maximum number of available information.

Figure 1:  Experimental Pearson correlations among the material properties at
the three scales considering all the silk reported in [33].

The EPR technique has returned a series of equations (expressions) for each
variable searched as output, represented in the figures with increasing
numbers on the horizontal axis in Fig. 2(d,e,f). These models represent the
Pareto front consisting of best formulae considering parsimony (simple
expression) and accuracy in a single formulae space exploration. From the
whole expression set reported in the Appendix, the most suitable equations
(reported in Table 2) were selected considering not only parsimony and
performance, but also analyzing the physical interpretation of the
experimental matrix correlations in Fig. 1).

##  5 Results and Discussion

###  5.1 Meso from micro

Firstly, the meso scale properties have been calculated using all the micro
scale quantities as independent variables (see Fig. 2(a)) according to the
recalled equation $Y=\sum_{j=1}^{m}f(x,a_{j})+a_{0}$. The results of the
accuracy are reported in Fig. 2(d) and the resulting functional dependencies
are reported in Table 2. In Figure 2 we also report the mean of the absolute
value of the Pearson correlation $\overline{|\rho|}$ among each dependent
variable and the independent variables as a measure of the effective existence
of experimental correlation among the considered variables. We remark that the
EPR algorithm, avoiding overfitting, correctly found functional dependence
that exhibits low performances in terms of $R^{2}$ in the case where also the
experimental correlations are low.

As a general qualitative description of the considered numerical tests (see in
particular the variables $A$ and $C$ in Fig. 2(d)), we may typically
distinguish two regimes of the performance curves. In the first regime, the
performance increases rapidly with the number of expressions and thus with the
model’s complexity. In the second regime, the performance curves stabilize in
a saturation band. This indicates an easy way of selecting an optimal model
complexity.

Regarding the selected functional dependence, first, we observe that the
crystallinity $A$ decreases with $b$, in accordance with the general
correlation matrix (Fig. 1). The presence of the bias term is coherent with
the value of $R^{2}=11\%$, since, as recalled before, the bias may compensate
for the lack of relevant inputs in the model.

The birefringence $B$ shows a very low accuracy $R^{2}<5\%$. This can be
interpreted by observing that the experimental result shows a very Pearson
correlation $\overline{|\rho|}=0.07$. We remark that, in this case, the EPR
method avoided data overfitting that could have resulted in more performant,
but physically misleading expressions. We therefore conclude, in this case,
that this meso scale quantity cannot be predicted starting from the considered
micro scale properties and we consider instead the $B$ as an independent
variable to compute the macro scale quantities in the following.

On the other hand, in the case of the Thermal Degradation Temperature $C$, the
EPR found expressions with higher $R^{2}$. In this case, the selected
expression provides a quantitative estimate of the target quantity with a
trend increasing with $a$ and $c$, according to the experimental correlation
matrix (Fig. 1).

Figure 2:  Prediction of material properties using two scales at a time: (a,d) meso from micro experimental properties, (b,e) macro from meso experimental properties, (c,f) macro from micro experimental properties. (a,b,c) Scheme of the strategy to obtain each quantity: solid (Dashed) box indicates experimental (obtained from EPR) quantities. (d,e,f) EPR Model performance in terms of $R^{2}$ plotted against the number of the found model. Scale | Expression | $R^{2}(\%)$ | # model
---|---|---|---
Meso from micro | $A=3.5562\frac{1}{{b}}+0.10262$ | 11 | 2
| $C=3186.7046\frac{1}{{a}}+0.86787{a}{c}^{0.5}+45.2672$ | 23.54 | 4
Macro from meso | $\mathbb{A}=0.091301\frac{B^{0.5}}{{A}}+29.2668{A}$ | 9.01 | 4
| $\mathbb{B}=0.013837\frac{B^{0.5}}{{A}}+0.014276{A}{C}$ | 12.81 | 5
| $\mathbb{C}=0.81928\frac{{A}^{0.5}{C}}{B}$ | 22.37 | 4
Macro from Micro | $\mathbb{D}=+0.061926\frac{{b}}{{c}}+0.0047393$ | 43.35 | 3
Table 2:  Prediction across two scales: selected explicit expressions

###  5.2 Macro-meso

As a second test, the macro properties have been calculated using all the meso
scale quantities as independent variables (see Fig. 2(b)) according to the
usual $\mathbb{Y}=\sum_{j=1}^{m}f(X,a_{j})+a_{0}$. The results of the accuracy
are reported in Fig. 2(e) and the resulting functional dependencies are
reported in Table 2.

Regarding the Young’s modulus $\mathbb{A}$, the chosen expression correctly
reports the monotonic growth with crystallinity $A$, as can be immediately
deduced by comparing the derivative of the expression for $A>0$ with the
experimental correlation matrix.

Regarding the limit stress $\mathbb{B}$, the selected expression correctly
reports the highest experimental correlation, namely the positive one with the
Birefringence $B$.

The expression chosen for the diameter $\mathbb{C}$ has the highest accuracy
among the macro-meso case ($R^{2}=22.37\%$) with a very simple expression
composed of only a single term that includes all three variables at the meso
scale. The correlation is positive for $A$ and $C$ and negative for $B$ in
accordance with the experiments. The relatively high performance of the EPR
method of this particular case corresponds to fairly relevant experimental
correlations $\rho=0.23,-0.39,0.22$ between $\mathbb{C}$ and $A,B,C$,
respectively.

Eventually, we consider the selected expression for the supercontraction
$\mathbb{D}$. In this case, we are not able to produce a good estimate of the
target output from the meso variables ($R^{2}<8\%$). This yet is in agreement
with the general matrix of experimental correlations with
$\overline{|\rho|}=0.07$.

###  5.3 Macro-micro

Eventually, we consider the possibility of direct dependence between macro and
micro variables. Thus, the macro properties have been calculated also using
all the micro scale quantities as independent variables (see Fig. 2(c))
according to equation $\mathbb{Y}=\sum_{j=1}^{m}f(x,a_{j})+a_{0}.$ The results
of the accuracy are reported in Fig. 2(f) and the resulting functional
dependencies are reported in Table 2.

In this case, regarding the Young’s modulus ($\mathbb{A}$) and the limit
stress ($\mathbb{B}$) the values of $R^{2}$ are generally very low and this is
in agreement with the lack of relevant experimental correlations among these
quantities. For the diameter ($\mathbb{C}$), the $R^{2}$ is only slightly
higher, again reflecting lightly higher experimental correlations. On the
other hand, the supercontraction $\mathbb{D}$ is predicted with a relatively
high accuracy ($R^{2}>40\%$), and the selected expression ($R^{2}=43.35\%$)
provides a reasonably precise quantitative estimate of the supercontraction.
Moreover, the expression is very simple and includes the two most relevant
experimental correlations between the supercontraction and the micro scale
properties, i. e. the positive one with $b$ and the negative one with $c$.

By employing this direct macro-micro deduction, from one side we obtain a
relatively precise estimation of the supercontraction property that was
missing from the macro-meso analysis, but from a modelling point of view, we
deduce the possibility of modelling the supercontrction as a macro variable
with a direct functional dependence from the micro ones. Moreover, the low
accuracy in predicting the other macro variables
($\mathbb{A}$,$\mathbb{B}$,$\mathbb{C}$) directly from the micro ones,
enlighten the importance of the meso scale structures in generally determining
the macro properties of the material, as expected from the classical
hierarchical dependence. For the particular case of the spider silk, this
reflects established results in literature pointing out the dependence of the
silk thread macroscopic behavior from the secondary structures of the proteins
[61, 50], here described by the mesoscale variables.

##  6 Theoretical vs experimental correlations

While the objective of this paper is general and mainly related to the
exhibited possibilities of obtaining information on the considered physical
properties, in this section, we show operatively the possibility of comparing
the experimental and theoretical results. In general, we anticipate that the
coefficient of determination of expressions found by the EPR method is
generally low if compared with other frameworks where EPR was applied [62,
63], but this was expected for the study case of spider silks, as in
biological materials a high intrinsic variability for experimental
observations is known [64]. Also, for this particular material, a meaningful
variability of the mechanical property of silks taken from the same individual
under similar conditions is well recognized (see e.g. [65]). Further, the
characteristics of the spider silks have high sensitivity to a large number of
parameters, among which starvation, reeling speed [65] other than the more
expected spider species [56], type of silk (among the up seven different ones
that some spiders can produce [57, 58]), environmental conditions [59] and
hydration conditions [60]. In a very recent work [66], the variability of
spider silk properties has been directly compared to that of carbon fibers,
and significantly higher variability in spider silk in all properties
considered has been reported. For these reasons, even if the $R^{2}$ of the
expressions found by means EPR is generally not as high as other frameworks,
the performance of the data modelling strategies are considered satisfactory.
On the other hand, we remark again that the analysis of the effective
statistical properties of the theoretical results of spider silks is out of
the aim of this paper and is the subject of a forthcoming research of the
authors.

Going to the considered case of spider silks, we are now in the position of
deducing the theoretical meso and macro response based on the only micro
properties. According to the previous reasoning, also the meso variable $B$
(birefringence) is considered here as an independent variable. Moreover,
coherently with the hierarchical assumption of our model, we first deduce the
meso scale variables by the Micro ones and then, based on previous analytical
results, we deduce the macro variables. On the other hand, in the special case
of the supercontraction $\mathbb{D}$, we assume that it directly depends on
micro variables. Notice that all these relevant physical information have been
deduced by previous data modelling.

Figure 3:  Pearson correlations among the material properties at the three
scales. (a) Experimental correlations obtained considering a subset of silks
for which all the analyzed properties are reported simultaneously (b)
Correlations among the material properties obtained from the data modelling
EPR approach (macro and meso) starting from the known micro experimental
properties.

Regarding the experimental data, we consider a subset of the silks analyzed by
[33] and in particular only those for which all the 10 considered properties
(see Table 1) are known simultaneously (the so obtained subset consists of 35
silks). As a possible comparison between the theoretical and experimental data
sets, we consider the Pearson correlation coefficients for each pair of
properties. The results reported in Fig. 3 show a significant correspondence
almost extensible to all the data and a satisfying result in terms of the
values of the correlation coefficients. To get a global comparison we also
adopt a positive definite relative error

$e_{r}=\frac{|e_{a}|}{e_{m}}$ (4)

where $e_{a}=\rho_{t}-\rho_{e}$ is the absolute error, $\rho_{t}$ and
$\rho_{e}$ are the theoretical and experimental Pearson coefficients,
respectively. Here $e_{m}$ is the mean error that since $\rho_{m}$ and
$\rho_{e}$ range in the interval $(-1,1)$, we assume as $e_{m}=1$. The average
value of the relative error by considering all the possible pairs of the
correlation matrix $\overline{e_{r}}=0.33$, with $0<\overline{e_{r}}<2$,
indicates that the functional dependence found by the EPR method reproduces in
a reasonably accurate way the experimental correlations among the properties
of spider silks.

Eventually, as evidenced in [16], an important extension of the proposed
approach would be to consider a Bayesian framework for the uncertainty
quantification in order to compute each output in terms of statistical
distribution with a mean and a confidence interval by also taking into account
the input data variability.

##  7 Conclusions

We showed the possibility of deducing, based on a Genetic Programming
approach, data modelling techniques particularly suitable for the deduction of
analytical models for multiscale problems. Our approach is based on the
Evolutionary Polynomial Regression (EPR) method, which as we showed lets us
deduce models that are both accurate and simple, able to describe the
dependence of macro scale variables from the one at lower scales, in their
hierarchical order. The best performing models are those located on the Pareto
dominance front, which takes into account both accuracy and parsimony and are
returned by the EPR algorithm. The final choice of the model can be then based
on physical considerations.

To explicitly show the possibility of acquiring physical insight in a complex
multiscale problem, and to evidence the key advantages of our multiscale
approach compared to classical, non-physically based techniques, we referred
to the materials science field and in particular to the spider silk: a
biological material with exceptional properties hugely analyzed also in the
spirit of bioinspiration. The choice of this specific case is due to the
observation that such remarkable properties are strictly based on an
evolutionary hierarchical optimization and the macroscopic spider silk
behavior is the result of noteworthy mesostructures emerging from the
aggregation of amino acids at the molecular scale. For this intriguing and
very complex material optimization case many phenomena are still unclear. We
then used this paradigmatic example to show how the presented data modelling
approach can be useful in several directions: determine dependent and
independent variables, indicate their hierarchical organization, deduce
explicit relations among different groups of variables. Furthermore, we showed
that the proposed approach let us overcome the overfitting problem typically
observed in the analysis of big data within the ANN framework diffusely
adopted in this field.

Based on this, new physical knowledge is acquired, that can be used as a
starting point for determining new analytical models, suggesting new
experiments, and applying more focused data modelling analysis. In this sense,
we assume that Machine Learning or Artificial Intelligence can be impactful
for scientific knowledge only if the data modelling approaches are in
continuous synergy with the scientific interpretation of the results. We argue
thus that a new mixed GP - theoretical approach can be a new fruitful approach
in material science, but also in fields as diverse as biology and medicine.

## Declaration of Competing Interest

The authors declare that they have no known competing financial interests or
personal relationships that could have appeared to influence the work reported
in this paper.

## Acknowledgments

Funding: VF and GP have been supported by GNFM (INdAM), GP has been supported
by the Italian Ministry MIUR-PRIN project 2017KL4EF3 and by PNRR, National
Centre for HPC, Big Data and Quantum Computing (CN00000013) - Spoke 5
”Environment and Natural Disasters” and NMP by the European Commission under
the FET Open “Boheme” grant no. 863179 and by the Italian Ministry of
Education MIUR under the PRIN-20177TTP3S.

## References

  * [1] Carlos Bustamante, Jed C Macosko, and Gijs JL Wuite.  Grabbing the cat by the tail: manipulating molecules one by one.  Nature Reviews Molecular Cell Biology, 1(2):130–136, 2000.
  * [2] Bharat Bhushan.  Nanotribology and nanomechanics of mems/nems and biomems/bionems materials and devices.  Microelectronic Engineering, 84(3):387–412, 2007.
  * [3] Jianxiong Zhu, Xinmiao Liu, Qiongfeng Shi, Tianyiyi He, Zhongda Sun, Xinge Guo, Weixin Liu, Othman Bin Sulaiman, Bowei Dong, and Chengkuo Lee.  Development trends and perspectives of future sensors and mems/nems.  Micromachines, 11(1):7, 2019.
  * [4] Qiang Chen and Nicola M. Pugno.  Bio-mimetic mechanisms of natural hierarchical materials: A review.  Journal of the Mechanical Behavior of Biomedical Materials, 19:3–33, 2013.
  * [5] Menon Carlo and Sitti Metin.  A biomimetic climbing robot based on the gecko.  Journal of Bionic Engineering, 3(3):115–125, 2006.
  * [6] Andreas Solga, Zdenek Cerman, Boris F Striffler, Manuel Spaeth, and Wilhelm Barthlott.  The dream of staying clean: Lotus and biomimetic surfaces.  Bioinspiration & biomimetics, 2(4):S126, 2007.
  * [7] Won-Gyu Bae, Hong Nam Kim, Doogon Kim, Suk-Hee Park, Hoon Eui Jeong, and Kahp-Yang Suh.  25th anniversary article: scalable multiscale patterned structures inspired by nature: the role of hierarchy.  Advanced Materials, 26(5):675–700, 2014.
  * [8] Euan A Ashley.  Towards precision medicine.  Nature Reviews Genetics, 17(9):507–522, 2016.
  * [9] Peng Zhang, Jialin Su, and Ulrike Mende.  Cross talk between cardiac myocytes and fibroblasts: from multiscale investigative approaches to mechanisms and functional consequences.  American Journal of Physiology-Heart and Circulatory Physiology, 303(12):H1385–H1396, 2012.
  * [10] Rebecca McLennan, Louise Dyson, Katherine W Prather, Jason A Morrison, Ruth E Baker, Philip K Maini, and Paul M Kulesa.  Multiscale mechanisms of cell migration during development: theory and experiment.  Development, 139(16):2935–2944, 2012.
  * [11] Hantao Ji, William Daughton, Jonathan Jara-Almonte, Ari Le, Adam Stanier, and Jongsoo Yoo.  Magnetic reconnection in the era of exascale computing and multiscale experiments.  Nature Reviews Physics, 4(4):263–282, 2022.
  * [12] S.S. Haykin.  Neural Networks: A Comprehensive Foundation.  International edition. Prentice Hall, 1999.
  * [13] O Giustolisi and D Laucelli.  Improving generalization of artificial neural networks in rainfall-runoff modelling.  Hydrological Sciences Journal, 50(3), 2005.
  * [14] Minkyung Baek, Frank DiMaio, Ivan Anishchenko, Justas Dauparas, Sergey Ovchinnikov, Gyu Rie Lee, Jue Wang, Qian Cong, Lisa N Kinch, R Dustin Schaeffer, et al.  Accurate prediction of protein structures and interactions using a three-track neural network.  Science, 373(6557):871–876, 2021.
  * [15] Grace X Gu, Chun-Teh Chen, Deon J Richmond, and Markus J Buehler.  Bioinspired hierarchical composite design using machine learning: simulation, additive manufacturing, and experiment.  Materials Horizons, 5(5):939–945, 2018.
  * [16] Kevin Linka and Ellen Kuhl.  A new family of constitutive artificial neural networks towards automated model discovery.  Computer Methods in Applied Mechanics and Engineering, 403:115731, 2023.
  * [17] Zeliang Liu, MA Bessa, and Wing Kam Liu.  Self-consistent clustering analysis: an efficient multi-scale scheme for inelastic heterogeneous materials.  Computer Methods in Applied Mechanics and Engineering, 306:319–341, 2016.
  * [18] Hunter T Kollmann, Diab W Abueidda, Seid Koric, Erman Guleryuz, and Nahil A Sobh.  Deep learning for topology optimization of 2d metamaterials.  Materials & Design, 196:109098, 2020.
  * [19] Miguel A Bessa, Piotr Glowacki, and Michael Houlder.  Bayesian machine learning in metamaterial design: Fragile becomes supercompressible.  Advanced Materials, 31(48):1904845, 2019.
  * [20] Ahmet Cecen, Hanjun Dai, Yuksel C Yabansu, Surya R Kalidindi, and Le Song.  Material structure-property linkages using three-dimensional convolutional neural networks.  Acta Materialia, 146:76–84, 2018.
  * [21] Zhi-Lei Wang and Yoshitaka Adachi.  Property prediction and properties-to-microstructure inverse analysis of steels by a machine-learning approach.  Materials Science and Engineering: A, 744:661–670, 2019.
  * [22] W James Murdoch, Chandan Singh, Karl Kumbier, Reza Abbasi-Asl, and Bin Yu.  Definitions, methods, and applications in interpretable machine learning.  Proceedings of the National Academy of Sciences, 116(44):22071–22080, 2019.
  * [23] Mengnan Du, Ninghao Liu, and Xia Hu.  Techniques for interpretable machine learning.  Communications of the ACM, 63(1):68–77, 2019.
  * [24] Finale Doshi-Velez and Been Kim.  Towards a rigorous science of interpretable machine learning.  arXiv preprint arXiv:1702.08608, 2017.
  * [25] Christoph Molnar.  Interpretable machine learning.  Lulu. com, 2020.
  * [26] John R Koza.  Genetic programming, on the programming of computers by means of natural selection. a bradford book.  MIT Press, 1992.
  * [27] Orazio Giustolisi and Dragan A Savic.  A symbolic data-driven technique based on evolutionary polynomial regression.  Journal of Hydroinformatics, 8(3):207–222, 2006.
  * [28] Orazio Giustolisi and DA Savic.  Advances in data-driven analyses and modelling using epr-moga.  Journal of Hydroinformatics, 11(3-4):225–236, 2009.
  * [29] Sinan Keten, Zhiping Xu, Britni Ihle, and Markus J Buehler.  Nanoconfinement controls stiffness, strength and mechanical toughness of $\beta$-sheet crystals in silk.  Nature materials, 9(4):359–367, 2010.
  * [30] Federico Bosia, Tamer Abdalrahman, and Nicola M Pugno.  Investigating the role of hierarchy on the strength of composite materials: evidence of a crucial synergy between hierarchy and material mixing.  Nanoscale, 4(4):1200–1207, 2012.
  * [31] Tina Arndt, Gabriele Greco, Benjamin Schmuck, Jessica Bunz, Olga Shilkova, Juanita Francis, Nicola M Pugno, Kristaps Jaudzems, Andreas Barth, Jan Johansson, et al.  Engineered spider silk proteins for biomimetic spinning of fibers with toughness equal to dragline silks.  Advanced Functional Materials, 32(23):2200986, 2022.
  * [32] Yan Liu, Dan Luo, and Tie Wang.  Hierarchical structures of bone and bioinspired bone tissue engineering.  Small, 12(34):4611–4632, 2016.
  * [33] Kazuharu Arakawa, Nobuaki Kono, Ali D Malay, Ayaka Tateishi, Nao Ifuku, Hiroyasu Masunaga, Ryota Sato, Kousuke Tsuchiya, Rintaro Ohtoshi, Daniel Pedrazzoli, et al.  1000 spider silkomes: Linking sequences to silk physical properties.  Science Advances, 8(41):eabo6043, 2022.
  * [34] Orazio Giustolisi, Angelo Doglioni, Dragan A Savic, and BW Webb.  A multi-model approach to analysis of environmental phenomena.  Environmental Modelling & Software, 22(5):674–682, 2007.
  * [35] ORAZIO GIUSTOLISI and DRAGAN SAVIC.  A novel genetic programming strategy: evolutionary polynomial regression.  In Hydroinformatics: (In 2 Volumes, with CD-ROM), pages 787–794. World Scientific, 2004.
  * [36] Vilfredo Pareto.  Cours d’Economie politique, volume 1.  F. Rouge, 1896.
  * [37] Orazio Giustolisi.  Using genetic programming to determine chezy resistance coefficient in corrugated channels.  Journal of Hydroinformatics, 6(3):157–173, 2004.
  * [38] Gabriele Greco, Tina Arndt, Benjamin Schmuck, Juanita Francis, Fredrik G Bäcklund, Olga Shilkova, Andreas Barth, Nathalie Gonska, Gulaim Seisenbaeva, Vadim Kessler, et al.  Tyrosine residues mediate supercontraction in biomimetic spider silk.  Communications Materials, 2(1):43, 2021.
  * [39] José Pérez-Rigueiro, Manuel Elices, Gustavo R Plaza, and Gustavo V Guinea.  Basic principles in the design of spider silk fibers.  Molecules, 26(6):1794, 2021.
  * [40] Manuel Elices, Gustavo R. Plaza, José Pérez-Rigueiro, and Gustavo V. Guinea.  The hidden link between supercontraction and mechanical behavior of spider silks.  Journal of the Mechanical Behavior of Biomedical Materials, 4(5):658–669, 2011.
  * [41] Alexander Sponner, Wolfram Vater, Shamci Monajembashi, Eberhard Unger, Frank Grosse, and Klaus Weisshart.  Composition and Hierarchical Organisation of a Spider Silk.  PLoS ONE, 2(10):e998, October 2007.
  * [42] Sinan Keten and Markus J. Buehler.  Nanostructure and molecular mechanics of spider dragline silk protein assemblies.  Journal of The Royal Society Interface, 7(53):1709–1721, 2010.
  * [43] Janelle E Jenkins, Sujatha Sampath, Emily Butler, Jihyun Kim, Robert W Henning, Gregory P Holland, and Jeffery L Yarger.  Characterizing the secondary protein structure of black widow dragline silk using solid-state NMR and X-ray diffraction.  Biomacromolecules, 14(10):3472–3483, 2013.
  * [44] SF Li, AJ McGhie, and SL Tang.  New internal structure of spider dragline silk revealed by atomic force microscopy.  Biophysical journal, 66(4):1209–1212, 1994.
  * [45] Lukas Eisoldt, Andrew Smith, and Thomas Scheibel.  Decoding the secrets of spider silk.  Materials Today, 14(3):80–86, 2011.
  * [46] Cameron P. Brown, Jennifer MacLeod, Heinz Amenitsch, Fernando Cacho-Nerin, Harinderjit S. Gill, Andrew J. Price, Enrico Traversa, Silvia Licoccia, and Federico Rosei.  The critical role of water in spider silk and its consequence for protein mechanics.  Nanoscale, 3(9):3805, 2011.
  * [47] Andrea Nova, Sinan Keten, Nicola M. Pugno, Alberto Redaelli, and Markus J. Buehler.  Molecular and Nanostructural Mechanisms of Deformation, Strength and Toughness of Spider Silk Fibrils.  Nano Letters, 10(7):2626–2634, July 2010.
  * [48] Paul J. Flory and Burak Erman.  Theory of elasticity of polymer networks. 3.  Macromolecules, 15(3):800–806, 1982.
  * [49] Vincenzo Fazio, Domenico De Tommasi, Nicola Maria Pugno, and Giuseppe Puglisi.  Spider silks mechanics: Predicting humidity and temperature effects.  Journal of the Mechanics and Physics of Solids, 164:104857, 2022\.
  * [50] Cheryl Y Hayashi, Nichola H Shipley, and Randolph V Lewis.  Hypotheses that correlate the sequence, structure, and mechanical properties of spider silk proteins.  International journal of biological macromolecules, 24(2-3):271–275, 1999.
  * [51] Jeffery L Yarger, Brian R Cherry, and Arjan Van Der Vaart.  Uncovering the structure–function relationship in spider silk.  Nature Reviews Materials, 3(3):1–11, 2018.
  * [52] Michael F Ashby and Kara Johnson.  Materials and design: the art and science of material selection in product design.  Butterworth-Heinemann, 2013.
  * [53] D Porter, J Guan, and F Vollrath.  Spider silk: super material or thin fibre?  Advanced Materials, 25(9):1275–1279, 2013.
  * [54] Robert W. Work.  Dimensions, birefringences, and force-elongation behavior of major and minor ampullate silk fibers from orb-web-spinning spiders—the effects of wetting on these properties.  Textile Research Journal, 47(10):650–662, 1977.
  * [55] Vincenzo Fazio, Nicola Maria Pugno, and Giuseppe Puglisi.  “water to the ropes”: A predictive model for the supercontraction stress of spider silks.  Extreme Mechanics Letters, 61:102010, 2023.
  * [56] Cecilia Boutry and Todd Alan Blackledge.  Evolution of supercontraction in spider silk: structure–function relationship from tarantulas to orb-weavers.  Journal of Experimental Biology, 213(20):3505–3514, 2010.
  * [57] Fritz Vollrath.  Spider webs and silks.  Scientific American, 266(3):70–77, 1992.
  * [58] JM Gosline, CC Pollak, PA Guerette, A Cheng, ME DeMont, and MW Denny.  Elastomeric network models for the frame and viscid silks from the orb web of the spider araneus diadematus.  1994\.
  * [59] Gustavo R Plaza, Gustavo V Guinea, José Pérez-Rigueiro, and Manuel Elices.  Thermo-hygro-mechanical behavior of spider dragline silk: Glassy and rubbery states.  Journal of Polymer Science Part B: Polymer Physics, 44(6):994–999, 2006.
  * [60] Ingi Agnarsson, Cecilia Boutry, Shing-Chung Wong, Avinash Baji, Ali Dhinojwala, Andrew T Sensenig, and Todd A Blackledge.  Supercontraction forces in spider dragline silk depend on hydration rate.  Zoology, 112(5):325–331, 2009.
  * [61] Kenjiro Yazawa, Ali D Malay, Hiroyasu Masunaga, Y Norma-Rashid, and Keiji Numata.  Simultaneous effect of strain rate and humidity on the structure and mechanical behavior of spider silk.  Communications Materials, 1(1):10, 2020.
  * [62] Luigi Berardi, Orazio Giustolisi, Zoran Kapelan, and DA Savic.  Development of pipe deterioration models for water distribution systems using epr.  Journal of Hydroinformatics, 10(2):113–126, 2008.
  * [63] Laura Enríquez, Daniele Biagio Laucelli, Juan Saldarriaga, and Orazio Giustolisi.  Encapsulating the chlorine transport in water distribution networks by a symbolic formula using epr.  PREPRINT ?
  * [64] Douglas Cook, Margaret Julias, and Eric Nauman.  Biological variability in biomechanical engineering research: Significance and meta-analysis of current modeling practices.  Journal of biomechanics, 47(6):1241–1250, 2014.
  * [65] Bo Madsen, Zheng Zhong Shao, and Fritz Vollrath.  Variability in the mechanical properties of spider silks on three levels: interspecific, intraspecific and intraindividual.  International journal of biological macromolecules, 24(2-3):301–306, 1999.
  * [66] Gabriele Greco, Hamideh Mirbaha, Benjamin Schmuck, Anna Rising, and Nicola M Pugno.  Artificial and natural silk materials have high mechanical property variability regardless of sample size.  Scientific Reports, 12(1):3507, 2022.

##  Appendix A EPR Expressions

###  A.1 Meso-micro

####  A.1.1 Cristallinity

$\displaystyle A$ $\displaystyle=0.19253$ (.1) $\displaystyle A$
$\displaystyle=3.5562\frac{1}{{b}}+0.10262$ (.2) $\displaystyle A$
$\displaystyle=3.9403\frac{1}{{b}}+0.0097339{c}+0.013178$ (.3) $\displaystyle
A$ $\displaystyle=0.0010868{b}+2.0441\frac{{c}^{0.5}}{{b}}$ (.4)
$\displaystyle A$
$\displaystyle=0.007611\frac{{b}}{{a}^{0.5}}+1.961\frac{{c}^{0.5}}{{b}}$ (.5)
$\displaystyle A$
$\displaystyle=0.0077062\frac{{b}}{{a}^{0.5}}+2.4415\frac{1}{{b}}+0.38037\frac{{c}}{{b}}$
(.6) $\displaystyle A$
$\displaystyle=0.0079125\frac{{b}}{{a}^{0.5}}+10.8129\frac{1}{{c}{b}}+0.50353\frac{{c}}{{b}}$
(.7)

####  A.1.2 Birefringence

$\displaystyle B$ $\displaystyle=45.0698$ (.1) $\displaystyle B$
$\displaystyle=88.7638\frac{1}{{b}^{0.5}}+31.2338$ (.2) $\displaystyle B$
$\displaystyle=188.8929\frac{1}{{a}}+54.9397\frac{1}{{b}^{0.5}}+31.4271$ (.3)
$\displaystyle B$
$\displaystyle=571.6574\frac{1}{{a}}+1.7581\frac{{a}}{{b}^{0.5}}+19.1995$ (.4)
$\displaystyle B$
$\displaystyle=120.3924\frac{{b}^{0.5}}{{a}}+18.2244\frac{{a}}{{b}}+7.3676$
(.5) $\displaystyle B$
$\displaystyle=135.0812\frac{{b}^{0.5}}{{a}}+19.3865\frac{{a}}{{b}}+0.036473{a}+2.329$
(.6) $\displaystyle B$
$\displaystyle=141.865\frac{{b}^{0.5}}{{a}}+17.6973\frac{{a}}{{b}}+0.69564\frac{{a}}{{b}^{0.5}}$
(.7)

####  A.1.3 Thermal degradation temperature

$\displaystyle C$ $\displaystyle=226.8169$ (.1) $\displaystyle C$
$\displaystyle=3.1984{c}+201.1538$ (.2) $\displaystyle C$
$\displaystyle=785.1137\frac{1}{{c}}+13.9926{c}+13.4973$ (.3) $\displaystyle
C$ $\displaystyle=3186.7046\frac{1}{{a}}+0.86787{a}{c}^{0.5}+45.2672$ (.4)
$\displaystyle C$
$\displaystyle=2482.7659\frac{1}{{a}}+383.6809\frac{1}{{c}}+0.27232{a}{c}+25.2964$
(.5)

###  A.2 Macro-meso

####  A.2.1 Young’s Modulus

$\displaystyle\mathbb{A}$ $\displaystyle=9.4674$ (.1)
$\displaystyle\mathbb{A}$ $\displaystyle=0.86198B^{0.5}+3.7949$ (.2)
$\displaystyle\mathbb{A}$ $\displaystyle=0.60871\frac{1}{{A}}+29.105{A}$ (.3)
$\displaystyle\mathbb{A}$
$\displaystyle=0.091301\frac{B^{0.5}}{{A}}+29.2668{A}$ (.4)
$\displaystyle\mathbb{A}$
$\displaystyle=1.3608\frac{B^{0.5}}{{A}{C}^{0.5}}+29.2732{A}$ (.5)
$\displaystyle\mathbb{A}$
$\displaystyle=20.1089\frac{B^{0.5}}{{A}{C}}+1.9687{A}{C}^{0.5}$ (.6)
$\displaystyle\mathbb{A}$
$\displaystyle=2.3655\frac{B}{{A}{C}}+0.098189\frac{1}{{A}}+2.0067{A}{C}^{0.5}$
(.7) $\displaystyle\mathbb{A}$
$\displaystyle=27.51\frac{1}{{A}{C}}+2.2479\frac{B}{{A}{C}}+2.0041{A}{C}^{0.5}$
(.8) $\displaystyle\mathbb{A}$
$\displaystyle=10.3522\frac{B^{0.5}}{{A}{C}}+1.3585\frac{B}{{A}{C}}+1.9953{A}{C}^{0.5}$
(.9)

####  A.2.2 Limit Stress

$\displaystyle\mathbb{B}$ $\displaystyle=1$ (.1) $\displaystyle\mathbb{B}$
$\displaystyle=0.14344B^{0.5}+0.2349$ (.2) $\displaystyle\mathbb{B}$
$\displaystyle=0.09374\frac{1}{{A}}+3.1305{A}$ (.3) $\displaystyle\mathbb{B}$
$\displaystyle=0.013904\frac{B^{0.5}}{{A}}+3.1765{A}$ (.4)
$\displaystyle\mathbb{B}$
$\displaystyle=0.013837\frac{B^{0.5}}{{A}}+0.014276{A}{C}$ (.5)
$\displaystyle\mathbb{B}$
$\displaystyle=0.011517\frac{B^{0.5}}{{A}}+0.0017796{A}{C}B^{0.5}+0.20612$
(.6) $\displaystyle\mathbb{B}$
$\displaystyle=0.021735\frac{1}{{A}}+0.0090613\frac{B^{0.5}}{{A}}+0.001894{A}{C}B^{0.5}+0.14047$
(.7) $\displaystyle\mathbb{B}$
$\displaystyle=0.012666\frac{B^{0.5}}{{A}}+0.0048689{A}{C}+0.001261{A}{C}B^{0.5}+0.09517$
(.8)

####  A.2.3 Diameter

$\displaystyle\mathbb{C}$ $\displaystyle=2$ (.1) $\displaystyle\mathbb{C}$
$\displaystyle=81.9474\frac{1}{B}$ (.2) $\displaystyle\mathbb{C}$
$\displaystyle=177.4203\frac{{A}^{0.5}}{B}+0.051737$ (.3)
$\displaystyle\mathbb{C}$ $\displaystyle=0.81928\frac{{A}^{0.5}{C}}{B}$ (.4)
$\displaystyle\mathbb{C}$
$\displaystyle=0.00021001{C}+0.80165\frac{{A}^{0.5}{C}}{B}$ (.5)
$\displaystyle\mathbb{C}$
$\displaystyle=0.0037544\frac{{C}}{B^{0.5}}+0.76892\frac{{A}^{0.5}{C}}{B}$
(.6) $\displaystyle\mathbb{C}$
$\displaystyle=0.011782\frac{{C}}{B}+0.0025068\frac{{C}}{B^{0.5}}+0.76003\frac{{A}^{0.5}{C}}{B}$
(.7) $\displaystyle\mathbb{C}$
$\displaystyle=0.028813\frac{{C}}{B}+0.69464\frac{{A}^{0.5}{C}}{B}+0.010027\frac{{A}^{0.5}{C}}{B^{0.5}}$
(.8)

####  A.2.4 Supercontraction

$\displaystyle\mathbb{D}$ $\displaystyle=0.31695$ (.1)
$\displaystyle\mathbb{D}$ $\displaystyle=1.4429\frac{1}{B^{0.5}}+0.093513$
(.2) $\displaystyle\mathbb{D}$
$\displaystyle=11.8488\frac{{A}^{0.5}}{B}+0.18064$ (.3)
$\displaystyle\mathbb{D}$
$\displaystyle=0.00085048{C}+10.8458\frac{{A}^{0.5}}{B}$ (.4)
$\displaystyle\mathbb{D}$
$\displaystyle=0.00086571{C}+158.9386\frac{{A}^{0.5}}{{C}^{0.5}B}$ (.5)
$\displaystyle\mathbb{D}$
$\displaystyle=0.00050622{C}+3.8693e-05{C}B^{0.5}+12.6449\frac{{A}^{0.5}}{B}$
(.6) $\displaystyle\mathbb{D}$
$\displaystyle=0.00067818{C}+2.1548e-05{C}B^{0.5}+172.7934\frac{{A}^{0.5}}{{C}^{0.5}B}$
(.7)

###  A.3 Macro-Micro

####  A.3.1 Young’s Modulus

$\displaystyle\mathbb{A}$ $\displaystyle=8.5718$ (.1)
$\displaystyle\mathbb{A}$ $\displaystyle=+150.4503\frac{1}{{b}}+4.8595$ (.2)
$\displaystyle\mathbb{A}$
$\displaystyle=+103.8759\frac{1}{{a}}+93.8966\frac{1}{{b}}+3.4673$ (.3)
$\displaystyle\mathbb{A}$
$\displaystyle=+18.4325\frac{{c}}{{a}}+36.55\frac{1}{{c}}$ (.4)
$\displaystyle\mathbb{A}$
$\displaystyle=+111.3405\frac{{c}}{{a}{b}^{0.5}}+37.8293\frac{1}{{c}}$ (.5)
$\displaystyle\mathbb{A}$
$\displaystyle=+712.237\frac{{c}}{{a}{b}}+5.6688\frac{{b}^{0.5}}{{c}}$ (.6)
$\displaystyle\mathbb{A}$
$\displaystyle=+120.4316\frac{{c}}{{a}{b}^{0.5}}+21.7225\frac{1}{{c}}+2.2107\frac{{a}^{0.5}}{{c}}$
(.7) $\displaystyle\mathbb{A}$
$\displaystyle=+121.6927\frac{{c}}{{a}{b}^{0.5}}+0.40972\frac{{b}}{{c}}+114.29\frac{{a}^{0.5}}{{c}{b}}$
(.8)

####  A.3.2 Limit Stress

$\displaystyle\mathbb{B}$ $\displaystyle=1.2533$ (.1)
$\displaystyle\mathbb{B}$ $\displaystyle=+17.8802\frac{1}{{a}}+0.77347$ (.2)
$\displaystyle\mathbb{B}$
$\displaystyle=+5.1279\frac{1}{{a}^{0.5}}+3.14\frac{1}{{c}}+0.025086$ (.3)
$\displaystyle\mathbb{B}$
$\displaystyle=+146.7895\frac{1}{{a}^{0.5}{b}}+0.014999{b}+0$ (.4)
$\displaystyle\mathbb{B}$
$\displaystyle=+0.26238\frac{{b}}{{a}^{0.5}{c}^{0.5}}+24.7876\frac{1}{{b}}+0$
(.5) $\displaystyle\mathbb{B}$
$\displaystyle=+0.29999\frac{{b}}{{a}}+0.056394\frac{{b}}{{c}}+24.7595\frac{1}{{b}}+0$
(.6) $\displaystyle\mathbb{B}$
$\displaystyle=+132.185\frac{1}{{a}^{0.5}{b}}+0.052774\frac{{b}}{{a}^{0.5}}+0.056459\frac{{b}}{{c}}+0.037522$
(.7) $\displaystyle\mathbb{B}$
$\displaystyle=+0.44217\frac{{b}}{{a}^{0.5}{c}}+49.6034\frac{{c}^{0.5}}{{a}^{0.5}{b}}+0.054095\frac{{b}}{{c}}$
(.8)

####  A.3.3 Diameter

$\displaystyle\mathbb{C}$ $\displaystyle=1.5415$ (.1)
$\displaystyle\mathbb{C}$ $\displaystyle=+44.5973\frac{1}{{b}}+0.43669$ (.2)
$\displaystyle\mathbb{C}$ $\displaystyle=+1030.4476\frac{1}{{a}{b}}+0.82562$
(.3) $\displaystyle\mathbb{C}$
$\displaystyle=+1367.643\frac{1}{{a}{b}}+0.013763{b}$ (.4)
$\displaystyle\mathbb{C}$
$\displaystyle=+8702.7135\frac{1}{{a}{c}{b}}+0.01363{b}+0.17322$ (.5)
$\displaystyle\mathbb{C}$
$\displaystyle=+9507.752\frac{1}{{a}{c}{b}}+0.0019587{c}{b}$ (.6)
$\displaystyle\mathbb{C}$
$\displaystyle=+9491.6515\frac{1}{{a}{c}{b}}+0.001894{c}{b}+7.7482e-05{a}{c}$
(.7) $\displaystyle\mathbb{C}$
$\displaystyle=+9609.6043\frac{1}{{a}{c}{b}}+0.0014633{c}{b}+7.4705e-05{a}^{0.5}{c}{b}$
(.8)

####  A.3.4 Supercontraction

$\displaystyle\mathbb{D}$ $\displaystyle=+0.32479$ (.1)
$\displaystyle\mathbb{D}$ $\displaystyle=+0.0055013{b}+0.097111$ (.2)
$\displaystyle\mathbb{D}$ $\displaystyle=+0.061926\frac{{b}}{{c}}+0.0047393$
(.3) $\displaystyle\mathbb{D}$
$\displaystyle=+0.53648\frac{1}{{c}}+0.050578\frac{{b}}{{c}}$ (.4)
$\displaystyle\mathbb{D}$
$\displaystyle=+0.8009\frac{1}{{c}}+0.0072816\frac{{a}^{0.5}{b}}{{c}}$ (.5)
$\displaystyle\mathbb{D}$
$\displaystyle=+15.9524\frac{1}{{c}{b}}+0.008755\frac{{a}^{0.5}{b}}{{c}}$ (.6)
$\displaystyle\mathbb{D}$
$\displaystyle=+430.4092\frac{1}{{a}{c}{b}}+0.0090635\frac{{a}^{0.5}{b}}{{c}}$
(.7)