# Automated discovery of interpretable hyperelastic material models for human
brain tissue with EUCLID

Moritz Flaschel  Huitian Yu  Nina Reiter  Jan Hinrichsen  Silvia Budday  Paul
Steinmann  Siddhant Kumar  Laura De Lorenzis  Department of Mechanical and
Process Engineering, ETH Zürich, 8092 Zürich, Switzerland  Department of
Mechanical Engineering, Friedrich-Alexander-Universität of Erlangen–Nürnberg,
91058 Erlangen, Germany  Department of Materials Science and Engineering,
Delft University of Technology, 2628 CD Delft, The Netherlands

###### Abstract

We propose an automated computational algorithm for simultaneous model
selection and parameter identification for the hyperelastic mechanical
characterization of human brain tissue. Following the motive of the recently
proposed computational framework EUCLID (Efficient Unsupervised Constitutive
Law Identitication and Discovery) and in contrast to conventional parameter
calibration methods, we construct an extensive set of candidate hyperelastic
models, i.e., a model library including popular models known from the
literature, and develop a computational strategy for automatically selecting a
model from the library that conforms to the available experimental data while
being represented as an interpretable symbolic mathematical expression. This
computational strategy comprises sparse regression, i.e., a regression problem
that is regularized by a sparsity promoting penalty term that filters out
irrelevant models from the model library, and a clustering method for grouping
together highly correlated and thus redundant features in the model library.
The model selection procedure is driven by labelled data pairs stemming from
mechanical tests under different deformation modes, i.e., uniaxial
compression/tension and simple torsion, and can thus be interpreted as a
supervised counterpart to the originally proposed EUCLID that is informed by
full-field displacement data and global reaction forces. The proposed method
is verified on synthetical data with artificial noise and validated on
experimental data acquired through mechanical tests of human brain specimens,
proving that the method is capable of discovering hyperelastic models that
exhibit both high fitting accuracy to the data as well as concise and thus
interpretable mathematical representations.

###### keywords:

constitutive models, hyperelasticity, brain tissue, interpretable models,
sparse regression

##  1 Introduction

Despite extensive research over the past decades, the characterization of the
mechanical response of biological tissues like those of the human brain
remains an active field with many open questions, see the review article by
Budday et al. (2020). The amount of data is usually limited by ethical
restrictions and by the availability of only few specimens with short
durability (Faber et al., 2022). The general objective is therefore to
leverage these limited data in an efficient manner to discover an adequate
mathematical description of the tissue response. The conventional strategy is
to a priori assume a constitutive model and to calibrate the tunable model
parameters by means of the experimental measurements. However, the a priori
choice of the constitutive model is driven by human experience and intuition;
thus, this classical strategy is highly susceptible to introducing modeling
errors that lead to poor fitting accuracy and/or the need for tedious and
time-consuming iterative model correction procedures. Data-driven methods and
machine learning approaches promise a powerful remedy. In this work, one of
such methods, which falls under the umbrella of methods denoted as EUCLID
(Efficient Unsupervised Constitutive Law Identification and Discovery)
(Flaschel et al., 2021), is used for automatically discovering symbolic
expressions of hyperelastic material models for human brain tissue based on
experimental measurements.

The general idea behind machine-learning-based material modeling is to choose
a model ansatz that exhibits a high expressiveness owing to a vast amount of
tunable parameters, and to calibrate these parameters by leveraging the
available data. Popular examples of versatile machine-learning-based material
models include neural networks (Ghaboussi et al., 1991), splines (Sussman and
Bathe, 2009), Gaussian processes (Frankel et al., 2020), and neural ordinary
differential equations (Tac et al., 2022). In the context of hyperelasticity,
these methods can be used to learn the characteristic strain energy density
function while special attention needs to be attributed to not violating
physical constraints like objectivity or (poly-)convexity of the strain energy
density function (As’ad et al., 2022; Klein et al., 2022a, b; Kalina et al.,
2022; Thakolkaran et al., 2022; Chen and Guilleminot, 2022; Linka and Kuhl,
2022; As’ad and Farhat, 2022; Tac et al., 2022). Until now, the majority of
the proposed machine learning models for hyperelasticy are informed by
artificially generated data, e.g., data generated through simulations at the
microscopic level. An exception is for example the work by Linka et al.
(2021), who train neural networks based on real experimental data for the
characterization of human brain tissue. In contrast to machine-learning-based
approaches stand model-free data-driven methods (Kirchdoerfer and Ortiz, 2016;
Ibañez et al., 2018), which seek to avoid the formulation of a material model
altogether by solving forward problems that are directly informed by the data.

Both the previously mentioned machine-learning-based approaches as well as
model-free data-driven methods have in common that the constitutive behavior
is encoded in a black box that does not allow for physical interpretation,
e.g., it is difficult to interpret the meaning of the many functions and
parameters that compose models encoded in machine learning tools like neural
networks. Furthermore, these methods typically rely on a vast amount of
labeled data pairs that are not available under experimental conditions. The
recognition of these issues motivated the development of the EUCLID method,
which aims to discover interpretable symbolic expressions of material models
using only experimentally available data like displacement fields and global
reaction forces. Initially applied to hyperelasticity (Flaschel et al., 2021;
Joshi et al., 2022) (see also the related and independent work by Wang et al.
(2021)), the framework was later extended to viscoelasticity (Marino et al.,
2023), elastoplasticity (Flaschel et al., 2022), and generalized standard
materials (Flaschel et al., 2023), see Flaschel (2023) for an overview.
Recently, discovering symbolic expressions for material models is gaining more
and more attention in the field, see, e.g., the works by Bomarito et al.
(2021); Kabliman et al. (2021); Park and Cho (2021); Abdusalamov et al.
(2023); Meyer and Ekre (2023).

In this work, EUCLID is used to automatically discover hyperelastic strain
energy density functions for describing the mechanical behavior of human brain
tissue. The idea is to construct a set of candidate material models called the
model library, which can be assembled for example from the vast available
literature on hyperelastic material models (see the non-comprehensive list of
reviewing articles by Boyce and Arruda (2000); Marckmann and Verron (2006);
Steinmann et al. (2012); Dal et al. (2021); He et al. (2021)), and to use a
specifically designed algorithm, based on sparse regression and feature
clustering, to select from the library a material model which is encoded by an
interpretable, simple mathematical expression and which is able to well
capture the mechanical behavior observed in the experimental data. Due to the
wet and shiny surface of brain tissues, displacement field data are difficult
to acquire; thus, against the original philosophy of EUCLID (Flaschel, 2023),
the method in this paper is not informed by full-field displacement and global
reaction force data. Instead, we rely here on labeled data pairs stemming from
uniaxial compression/tension and simple torsion tests of human brain tissue
specimens.

The paper is organized as follows: In Section 2, we construct the model
library in which EUCLID searches for a suitable material model. In Section 2,
the proposed algorithm for model selection is discussed in detail.
Subsequently, the method is tested on artificial data and actual experimental
data in Section 4 and Section 5, respectively, and conclusions are drawn in
Section 6.

Notation: Tensors and matrices may appear in compact or index notation, e.g.,
$\mathbf{F}$ or $F_{ij}$, respectively. In compact notation, first-order
tensors (vectors) and second-order tensors are denoted by bold letters, e.g.,
$\mathbf{F}$. The transpose of a tensor is denoted by $\square^{T}$. When
appearing in index notation, the order of the tensor equals the number of the
indices. If not stated otherwise, the Einstein convention for summation over
repeated indices is used in equations appearing in index notation, e.g.,
$a_{i}b_{i}=\sum_{i}a_{i}b_{i}$. Inner products are denoted by $\cdot$, e.g.,
$\bm{a}\cdot\bm{b}=a_{i}b_{i}$, and outer products by $\otimes$, e.g.,
$\\{\bm{a}\otimes\bm{b}\\}_{ij}=a_{i}b_{j}$. If no operation is indicated
between two tensors, the juxtaposition implies tensor contraction, e.g.,
$\\{\mathbf{F}^{T}\mathbf{F}\\}_{ij}=F_{ki}F_{kj}$. The trace of a tensor is
denoted by $\text{tr}(\square)$, e.g., $\text{tr}(\mathbf{C})=C_{ii}$, the
volumetric part by $\text{vol}(\square)$, the deviatoric part by
$\text{dev}(\square)$, and the determinant by $\text{det}(\square)$.

##  2 Model library

###  2.1 Strain energy density

The constitutive response of a hyperelastic material is completely
characterized by the material strain energy density function. In the spirit of
EUCLID, we construct a material model library, i.e., a large set of potential
candidate material models, by introducing a general parametric ansatz for the
unknown strain energy density. Under the assumption of incompressibility, the
strain energy density $W$ of a hyperelastic material is postulated as
(Holzapfel, 2000)

$W=\tilde{W}(\mathbf{F})-p\cdot(J-1),$ (1)

where $\mathbf{F}$ is the deformation gradient, $p$ is a scalar Lagrange
multiplier that can be physically interpreted as the hydrostatic (or
volumetric) pressure, and $J=\det\mathbf{F}\overset{!}{=}1$ is the determinant
of the deformation gradient. A sufficient condition to satisfy the objectivity
requirement for the material model is that $W$ depends on $\mathbf{F}$ through
the right Cauchy-Green tensor $\mathbf{C}=\mathbf{F}^{T}\mathbf{F}$, i.e.
$\tilde{W}(\mathbf{F})=\hat{W}(\mathbf{C})$. Assuming furthermore isotropic
material behavior, the contribution $\hat{W}$ can be expressed as a function
of the invariants of $\mathbf{C}$, defined as $I_{1}=\text{tr}(\mathbf{C})$,
$I_{2}=\frac{1}{2}(\text{tr}^{2}(\mathbf{C})-\text{tr}(\mathbf{C}^{2}))$ and
$I_{3}=\det(\mathbf{C})$. An equally objective possibility for isotropic
hyperelasticity is to express the strain energy density as a function of the
principal stretches $\lambda_{1}$, $\lambda_{2}$, $\lambda_{3}$ defined as the
eigenvalues of the right stretch tensor $\mathbf{U}$, which in turn is defined
through the polar decomposition $\mathbf{F}=\mathbf{R}\mathbf{U}$ (Holzapfel,
2000). As we seek to avoid a priori assumptions on the material response, we
do not limit ourselves to a strain energy density that depends solely on the
strain invariants or solely on the principal stretches, but instead consider
the general case in which the strain energy density includes a contribution
$\tilde{W}_{I}$ that depends on the strain invariants and a contribution
$\tilde{W}_{\lambda}$ that depends on the principal stretches

$W=\tilde{W}_{I}(I_{1},I_{2},I_{3})+\tilde{W}_{\lambda}(\lambda_{1},\lambda_{2},\lambda_{3})-p\cdot(J-1).$
(2)

By doing so, we obtain a highly general expression for the strain energy
density that encompasses many well-known hyperelastic constitutive models e.g.
of the Mooney-Rivlin or Ogden types.

To introduce a general parametric ansatz for $\tilde{W}_{I}$, we assume that
it can be expressed as a linear combination of a priori chosen nonlinear
feature functions

$\tilde{W}_{I}(I_{1},I_{2},I_{3})=\bm{\theta}_{I}\cdot\bm{Q}_{I}(I_{1},I_{2},I_{3}),$
(3)

where the nonlinear feature functions have been collected in the vector
$\bm{Q}_{I}$ and $\bm{\theta}_{I}$ is a vector of unknown real-valued material
parameters. We construct the feature vector $\bm{Q}_{I}$ such that it contains
features of the generalized Mooney-Rivlin model (Rivlin, 1947) and a
logarithmic feature as it appears in the Gent-Thomas model (Gent and Thomas,
1958)

$\bm{Q}_{I}(I_{1},I_{2},I_{3})=\underbrace{\left[(I_{1}-3)^{m}(I_{2}-3)^{n-m}:n\in\\{1,\dots,N_{\text{Mooney}}\\},m\in\\{0,\dots,n\\}\right]^{T}}_{\text{generalized
Mooney-Rivlin
features}}\oplus\underbrace{\left[\log\left(I_{2}/3\right)\right]}_{\text{logarithmic
feature}},$ (4)

where $\oplus$ denotes vector concatenation and the choice of
$N_{\text{Mooney}}$ dictates the size of the feature vector. This selection of
feature functions was proven in Flaschel et al. (2021) to have a high
approximation power. We choose $N_{\text{Mooney}}=3$ in this work, such that
$\bm{\theta}_{I}$ comprises 10 unknown parameters.

To obtain a general parametric ansatz for the principal stretch dependent
contribution $\tilde{W}_{\lambda}$, we consider the generalized Ogden model
(Ogden, 1972)

$\tilde{W}_{\lambda}(\lambda_{1},\lambda_{2},\lambda_{3})=\sum_{i=1}^{N_{\text{Ogden}}}\frac{2\mu_{i}}{\alpha_{i}^{2}}\left(\lambda_{1}^{\alpha_{i}}+\lambda_{2}^{\alpha_{i}}+\lambda_{3}^{\alpha_{i}}-3\right),$
(5)

where $\mu_{i}$ and $\alpha_{i}$ are unknown real-valued material parameters
and $N_{\text{Ogden}}$ is the number of considered terms in the generalized
Ogden model. In contrast to the model library chosen in Eq. 3, the strain
energy density of the generalized Ogden model depends nonlinearly on the
unknown material parameters. In general, this complicates the inference of the
material parameters from experimental measurements. Therefore, we assume in
the following an a priori fixed set of $N_{\text{Ogden}}$ distinct values of
$\alpha_{i}$. By choosing $N_{\text{Ogden}}$ sufficiently large, it is
expected that this assumption does not significantly restrict the versatility
of the model library. In this work, we choose $N_{\text{Ogden}}=2\cdot 10^{4}$
values of $\alpha_{i}$ evenly distributed between $-100$ and $100$ excluding
zero, i.e., $\alpha_{i}\in\\{-100,-99.99,\dots,-0.01,0.01,\dots,99.99,100\\}$.
With this assumption, $\mu_{i}$ remain the only unknowns in the model ansatz
for $\tilde{W}_{\lambda}$ which can now be written as a linear combination of
nonlinear feature functions

$\tilde{W}_{\lambda}(\lambda_{1},\lambda_{2},\lambda_{3})=\bm{\theta}_{\lambda}\cdot\bm{Q}_{\lambda}(\lambda_{1},\lambda_{2},\lambda_{3}),$
(6)

where the nonlinear feature functions have been collected in the vector
$\bm{Q}_{\lambda}$ with components

$\left\\{Q_{\lambda}(\lambda_{1},\lambda_{2},\lambda_{3})\right\\}_{i}=\lambda_{1}^{\alpha_{i}}+\lambda_{2}^{\alpha_{i}}+\lambda_{3}^{\alpha_{i}}-3,$
(7)

and $\bm{\theta}_{\lambda}$ is a vector of unknown real-valued material
parameters that are related to the parameters $\mu_{i}$ through
$\left\\{\theta_{\lambda}\right\\}_{i}=2\mu_{i}/\alpha_{i}^{2}$.

By defining

$\bm{\theta}=\begin{bmatrix}\bm{\theta}_{I}\\\ \bm{\theta}_{\lambda}\\\
\end{bmatrix},\quad\bm{Q}=\begin{bmatrix}\bm{Q}_{I}(I_{1},I_{2},I_{3})\\\
\bm{Q}_{\lambda}(\lambda_{1},\lambda_{2},\lambda_{3})\\\ \end{bmatrix},$ (8)

the model library for the strain energy density can be written as

$W=\bm{\theta}\cdot\bm{Q}-p\cdot(J-1).$ (9)

To ensure that the strain energy density describes physically meaningful
material behavior, we assume that $\theta_{i}\geq 0$ for all $i$. The non-
negativity of the material parameters is a sufficient (but not necessary)
requirement for stability, see Hartmann (2001).

###  2.2 Stress-strain relation

After having defined the strain energy density, a bijective relation between
the kinematic state and the stress state of the material is obtained through
its differentiation. Specifically, the Piola stress $\mathbf{P}$ is computed
by differentiating $W$ with respect to the deformation gradient

$\mathbf{P}=\frac{\partial
W}{\partial\mathbf{F}}=\frac{\partial\tilde{W}}{\partial\mathbf{F}}-p\mathbf{F}^{-T},$
(10)

where we used $\frac{\partial J}{\partial\mathbf{F}}=J\mathbf{F}^{-T}$. Noting
that
$\tilde{W}=\bm{\theta}\cdot\bm{Q}=\bm{\theta}_{I}\cdot\bm{Q}_{I}+\bm{\theta}_{\lambda}\cdot\bm{Q}_{\lambda}$,
we apply the chain rule to obtain

$\frac{\partial\tilde{W}}{\partial
F_{ij}}=\bm{\theta}_{I}\cdot\left(\frac{\partial\bm{Q}_{I}}{\partial
I_{a}}\frac{\partial I_{a}}{\partial
F_{ij}}\right)+\bm{\theta}_{\lambda}\cdot\left(\frac{\partial\bm{Q}_{\lambda}}{\partial\lambda_{b}}\frac{\partial\lambda_{b}}{\partial
F_{ij}}\right).$ (11)

The derivatives of the strain invariants and of the principal stretches are
provided in A.1 and A.2, respectively. The derivatives of the feature vectors
follow in A.3. The Cauchy stress $\mathbf{T}$ is calculated from the Piola
stress as $\mathbf{T}=\mathbf{P}\mathbf{F}^{T}$.

####  2.2.1 Uniaxial compression and tension

Experimental measurements of specimens under uniaxial compression and tension
deliver labeled data pairs in the form $(u_{\text{UT}},F)$, where
$u_{\text{UT}}$ is the longitudinal displacement at the displacement
controlled end of the specimen and $F$ is the resulting force. Data pairs of
the form $(\lambda_{\text{UT}},P_{11})$, where $\lambda_{\text{UT}}$ is the
longitudinal stretch applied to the specimen and $P_{11}$ is the longitudinal
normal component of the Piola stress, are obtained through

$\lambda_{\text{UT}}=\frac{u_{\text{UT}}+h}{h},\quad P_{11}=\frac{F}{\pi
r_{\text{out}}^{2}},$ (12)

where $h$ and $r_{\text{out}}$ are the height and the outer radius of the
undeformed specimen, respectively. The previously introduced material model
library can be used to derive a relation
$P_{11}(\lambda_{\text{UT}};\bm{\theta})$ that links the experimental inputs
$\lambda_{\text{UT}}$ and outputs $P_{11}$ depending on the material
parameters $\bm{\theta}$.

Applying a longitudinal stretch $\lambda_{\text{UT}}$ to the specimen implies
that $F_{11}=\lambda_{\text{UT}}$ and, due to symmetry, $F_{22}=F_{33}$. As a
result of the incompressibility assumption $\det\mathbf{F}\overset{!}{=}1$,
the deformation gradient under uniaxial compression/tension must hence be

$\displaystyle\mathbf{F}_{\text{UT}}=\begin{bmatrix}\lambda_{\text{UT}}&0&0\\\
0&\frac{1}{\sqrt{\lambda_{\text{UT}}}}&0\\\
0&0&\frac{1}{\sqrt{\lambda_{\text{UT}}}}\\\ \end{bmatrix}.$ (13)

The kinematic state is thus completely determined by the applied longitudinal
stretch.

To obtain the relationship $P_{11}(\lambda_{\text{UT}};\bm{\theta})$, the
unknown hydrostatic pressure $p$ in Eq. 10 needs to be computed. As the
material is allowed to freely deform in the $x_{2}$\- and $x_{3}$-directions,
it is $P_{22}=P_{33}\overset{!}{=}0$. This condition is used along with Eq. 10

$P_{33}=\frac{\partial\tilde{W}}{\partial
F_{33}}-pF^{-1}_{33}\overset{!}{=}0,$ (14)

to find the hydrostatic pressure

$p=\frac{\partial\tilde{W}}{\partial F_{33}}F_{33}.$ (15)

We hence obtain the desired relationship by substituting the pressure in Eq.
10

$P_{11}(\lambda_{\text{UT}};\bm{\theta})=\frac{\partial\tilde{W}}{\partial
F_{11}}-pF^{-1}_{11}=\frac{\partial\tilde{W}}{\partial
F_{11}}-\frac{F_{33}}{F_{11}}\frac{\partial\tilde{W}}{\partial F_{33}}\\\
=\bm{\theta}\cdot\underbrace{\left(\frac{\partial\bm{Q}}{\partial
F_{11}}-\frac{F_{33}}{F_{11}}\frac{\partial\bm{Q}}{\partial
F_{33}}\right)}_{\bm{Q}^{\prime}_{\text{UT}}\left(\lambda_{\text{UT}}\right)},$
(16)

where we denoted all the terms that depend on the longitudinal stretch as
$\bm{Q}^{\prime}_{\text{UT}}(\lambda_{\text{UT}})$. This can be interpreted as
a feature function derivative vector whose scalar product with the material
parameters computes the longitudinal normal component of the Piola stress.
More information for computing
$\bm{Q}^{\prime}_{\text{UT}}(\lambda_{\text{UT}})$ is provided in A.4.

####  2.2.2 Simple torsion

Now we consider a cylindrical specimen in a cylindrical coordinate system
$(r,\vartheta,z)$, where $r$ is the radial coordinate, $\vartheta$ is the
polar angle, and $z$ is the longitudinal coordinate. Assuming the specimen to
undergo simple torsion around the $z$-axis, experimental measurements on
specimens under simple torsion deliver labeled data pairs in the form
$(\phi,M)$, where $\phi$ is the applied angle of deformation at the sheared
end of the specimen and $M$ is the resulting torque. The twist $\psi$ of the
specimen can be deduced from the applied angle of deformation through
$\psi=\phi/h$. The data is transformed into data pairs of the form
$(\tilde{\psi},\tau)$, where $\tilde{\psi}=r_{\text{out}}\psi$ and
$\tau=M/r_{\text{out}}^{3}$ are respectively a normalized twist and normalized
torque.

After normalizing the radial coordinate through $\rho=r/r_{\text{out}}$, the
deformation gradient in the cylindrical coordinate system then reads
(Hartmann, 2001)

$\displaystyle\mathbf{F}_{\text{ST}}=\begin{bmatrix}F_{rr}&F_{r\vartheta}&F_{rz}\\\
F_{\vartheta r}&F_{\vartheta\vartheta}&F_{\vartheta z}\\\
F_{zr}&F_{z\vartheta}&F_{zz}\\\ \end{bmatrix}=\begin{bmatrix}1&0&0\\\
0&1&\rho\tilde{\psi}\\\ 0&0&1\\\ \end{bmatrix},$ (17)

The material model library can be used to derive a relation
$\tau(\tilde{\psi};\bm{\theta})$ that links the experimental inputs
$\tilde{\psi}$ and outputs $\tau$ depending on the material parameters
$\bm{\theta}$.

The normalized torque is given by

$\tau=\int_{0}^{1}\ 2\pi\rho^{2}T_{\vartheta z}\ \mathrm{d}\rho,$ (18)

where the component of the Cauchy stress $T_{\vartheta z}$ is, due to Eq. 17,

$T_{\vartheta z}=P_{\vartheta i}F_{zi}=P_{\vartheta z}.$ (19)

With reference to Eq. 10, we notice that there is no hydrostatic pressure
contribution to $P_{\vartheta z}$ under simple torsion, such that

$T_{\vartheta z}=\frac{\partial{W}}{\partial F_{\vartheta
z}}=\frac{\partial\tilde{W}}{\partial F_{\vartheta
z}}=\bm{\theta}\cdot\frac{\partial\bm{Q}}{\partial F_{\vartheta z}}.$ (20)

We finally obtain the desired relationship $\tau(\tilde{\psi};\bm{\theta})$ by
substituting Eq. 20 in Eq. 18

$\tau(\tilde{\psi};\bm{\theta})=\bm{\theta}\cdot\underbrace{\int_{0}^{1}\
2\pi\rho^{2}\frac{\partial\bm{Q}}{\partial F_{\vartheta z}}\
\mathrm{d}\rho}_{\bm{Q}^{\prime}_{\text{ST}}\left(\tilde{\psi}\right)},$ (21)

where we denoted all the terms that depend on the twist as
$\bm{Q}^{\prime}_{\text{ST}}(\tilde{\psi})$. This can be interpreted as a
feature function derivative vector whose scalar product with the material
parameters computes the normalized torque. More information for computing
$\bm{Q}^{\prime}_{\text{ST}}(\tilde{\psi})$ is provided in A.5.

##  3 Model discovery

After having defined the material model library, i.e., a versatile ansatz for
the hyperelastic strain energy density, we seek an optimal choice of the
material model parameters $\bm{\theta}$. Specifically, in the spirit of EUCLID
- see Flaschel et al. (2021) \- we seek a parameter vector $\bm{\theta}$ which
leads at the same time to a small mismatch between the material model
predictions and the experimental measurements, and to a simple discovered
material model, i.e. one with a small number of terms. Unlike in the original
work on EUCLID (Flaschel et al., 2021), in which the model selection was
driven by a physics-informed optimization problem that was fed by unlabelled
data, in this paper, labelled data pairs (i.e., stress versus strain data
pairs under uniaxial compression/tension, and twist versus torque data pairs
under simple torsion) are available to drive the model selection. Thus, the
method proposed here can be interpreted as a supervised counterpart to the
originally proposed EUCLID. The proposed algorithm for inferring $\bm{\theta}$
is described in the following and a step-by-step overview is provided in B,
see Fig. 10.

###  3.1 Model-data mismatch

First, a quantitative measure of the mismatch between the material model
predictions and the experimental measurements needs to be introduced. In the
case of uniaxial compression/tension, the measurements are provided in the
form of pairs $(\lambda_{\text{UT}}^{(l)},P_{11}^{(l)})$ for a number of load
steps $l=1,\dots,n_{\text{UT}}$. Eq. 16 provides a mapping between the input
and output data, thus leading to $n_{\text{UT}}$ equations

$P_{11}(\lambda_{\text{UT}}^{(l)};\bm{\theta})=P_{11}^{(l)}.$ (22)

Noting that each of the equations above depends linearly on $\bm{\theta}$, we
write them as a linear system of equations

$\mathbf{A}_{\text{UT}}\bm{\theta}=\bm{b}_{\text{UT}},$ (23)

with

$\mathbf{A}_{\text{UT}}=\begin{bmatrix}\left(\bm{Q}^{\prime}_{\text{UT}}(\lambda_{\text{UT}}^{(1)})\right)^{T}\\\
\vdots\\\
\left(\bm{Q}^{\prime}_{\text{UT}}(\lambda_{\text{UT}}^{(n_{\text{UT}})})\right)^{T}\\\
\end{bmatrix},\quad\bm{b}_{\text{UT}}=\begin{bmatrix}P_{11}^{(1)}\\\ \vdots\\\
P_{11}^{(n_{\text{UT}})}\\\ \end{bmatrix}.$ (24)

Analogously, in the case of simple torsion, the measurements are provided in
the form of pairs $(\tilde{\psi}^{(l)},\tau^{(l)})$ for a number of load steps
$l=1,\dots,n_{\text{ST}}$. Eq. 21 provides a mapping between the input and
output data, thus leading to $n_{\text{ST}}$ equations

$\tau(\tilde{\psi}^{(l)};\bm{\theta})=\tau^{(l)}.$ (25)

which can be expressed once again as a linear system of equations

$\mathbf{A}_{\text{ST}}\bm{\theta}=\bm{b}_{\text{ST}},$ (26)

with

$\mathbf{A}_{\text{ST}}=\begin{bmatrix}\left(\bm{Q}^{\prime}_{\text{ST}}(\tilde{\psi}^{(1)})\right)^{T}\\\
\vdots\\\
\left(\bm{Q}^{\prime}_{\text{ST}}(\tilde{\psi}^{(n_{\text{ST}})})\right)^{T}\\\
\end{bmatrix},\quad\bm{b}_{\text{ST}}=\begin{bmatrix}\tau^{(1)}\\\ \vdots\\\
\tau^{(n_{\text{ST}})}\\\ \end{bmatrix}.$ (27)

The linear systems in Eq. 23 and Eq. 26 are concatenated to the system
$\mathbf{A}\bm{\theta}=\bm{b}$ with

$\mathbf{A}=\begin{bmatrix}r_{\text{UT}}\mathbf{A}_{\text{UT}}\\\
r_{\text{ST}}\mathbf{A}_{\text{ST}}\\\ \end{bmatrix},\quad{\bf\it
b}=\begin{bmatrix}r_{\text{UT}}\bm{b}_{\text{UT}}\\\
r_{\text{ST}}\bm{b}_{\text{ST}}\\\ \end{bmatrix},$ (28)

where $r_{\text{UT}}>0$ and $r_{\text{ST}}>0$ scale the contributions of the
uniaxial compression/tension data and the simple torsion data, respectively.
During the experimental investigations conducted in the context of this work
(see, e.g., Section 5), the chosen range of applied deformation in uniaxial
compression/tension and simple torsion led to measurements
$\bm{b}_{\text{UT}}$ and $\bm{b}_{\text{ST}}$, respectively, that were
different in magnitude. Choosing $r_{\text{UT}}=0.3$ and $r_{\text{ST}}=1$ was
found to result in a proper scaling of the contributions from uniaxial
compression/tension and simple torsion, leading to similar maximum absolute
entries in the scaled measurements $r_{\text{UT}}\bm{b}_{\text{UT}}$ and
$r_{\text{ST}}\bm{b}_{\text{ST}}$. Finally, the total mismatch between the
material model prediction and the experimental measurements is quantified by
computing the mean squared error $\mathrm{MSE}$ of the residuals of the linear
system of equations

$\mathrm{MSE}(\bm{\theta})=\frac{1}{n_{\text{UT}}+n_{\text{ST}}}\|\mathbf{A}\bm{\theta}-\bm{b}\|_{2}^{2}.$
(29)

###  3.2 Feature scaling

To obtain a system of equations with dimensionless parameters and system
coefficients, we standardize the system of equations. Defining the mean and
standard deviation of a vector $\bm{x}\in\mathbb{R}^{n}$ as

$\mathrm{Mean}(\bm{x})=\frac{1}{n}\sum_{i}x_{i},\quad\mathrm{Std}(\bm{x})=\sqrt{\frac{1}{n}\sum_{i}\left(x_{i}-\mathrm{Mean}(\bm{x})\right)^{2}},$
(30)

the standardized system of equations reads

$\bar{\mathbf{A}}\bar{\bm{\theta}}=\bar{\bm{b}}$ (31)

with

$\bar{A}_{ij}=\frac{A_{ij}}{\mathrm{Std}(\bm{A}_{j})},\quad\bar{b}_{i}=\frac{b_{i}}{\mathrm{Std}(\bm{b})},\quad\bar{\theta}_{j}=\frac{\mathrm{Std}(\bm{A}_{j})}{\mathrm{Std}(\bm{b})}{\theta}_{j},$
(32)

where we defined $\bm{A}_{j}$ as the $j^{\text{th}}$ column of $\mathbf{A}$.
The mean squared error of the residuals of the standardized linear system of
equations is defined as

$\overline{\mathrm{MSE}}(\bar{\bm{\theta}})=\frac{1}{n_{\text{UT}}+n_{\text{ST}}}\|\bar{\mathbf{A}}\bar{\bm{\theta}}-\bar{\bm{b}}\|_{2}^{2}.$
(33)

###  3.3 Optimization problem

Our objectives are to calibrate the unknown material parameters such that the
material model is in agreement with the experimental data, and to discard
those terms in the model ansatz which have a minor influence on the fitting
accuracy by setting their corresponding parameters to zero. To this end, we
solve an $L_{1}$-regularized optimization problem (also denoted as sparse
regression problem) of the form

$\bar{\bm{\theta}}^{\text{opt}}=\operatornamewithlimits{arg\
min}_{\bar{\bm{\theta}}\geq\bm{0}}\left(\frac{1}{2}\overline{\mathrm{MSE}}(\bar{\bm{\theta}})+\lambda_{p}\|\bar{\bm{\theta}}\|_{1}\right),$
(34)

where $\overline{\mathrm{MSE}}(\bar{\bm{\theta}})$ is the previously
introduced mean squared error function that quantifies the mismatch between
the material model predictions and the experimental data, and
$\|\bar{\bm{\theta}}\|_{1}=\sum_{i}|\bar{\theta}_{i}|$ is an
$L_{1}$-regularization term, also denoted as Lasso (least absolute shrinkage
and selection operator). This term, first used for model selection by Frank
and Friedman (1993) and Tibshirani (1996) and later applied to problems in
dynamics by Brunton et al. (2016), can be interpreted as a convex
approximation of an operator that counts the number of nonzero entries in the
material parameter vector. Thus, adding the $L_{1}$-regularization term to the
minimization problem promotes sparsity of the material parameter vector. In
previous works (Flaschel et al., 2021, 2022), the more general
$L_{p}$-regularization term
$\|\bar{\bm{\theta}}\|_{p}^{p}=\sum_{i}|\bar{\theta}_{i}|^{p}$ with $0<p\leq
1$ was adopted to promote sparsity in the solution vector. The $L_{p}$ term
has the advantage that it converges to the operator that counts the number of
nonzero entries in the material parameter vector as $p$ approaches zero. Thus,
it can be considered as a better measure of sparsity than the
$L_{1}$-regularization term. However, the $L_{p}$ term is non-convex for
$p<1$; this implies that the optimization problem becomes non-convex and thus
multiple solutions are to be expected, corresponding to local minima of the
objective function. In our previous papers (Flaschel et al., 2021, 2022), this
issue was addressed by solving the optimization problem multiple times and
selecting the minimum of the computed local minima, which obviously increases
the computational complexity and cost of the minimization problem. In this
work, we rely on the $L_{1}$ regularization term and thus solve a convex
minimization problem which thus admits a unique solution. The influence of the
regularization is governed by the choice of the weighting factor
$\lambda_{p}\geq 0$. As proposed in Flaschel et al. (2022, 2023); Marino et
al. (2023), an appropriate choice of $\lambda_{p}$ for striking a balance
between fitting accuracy and mathematical complexity of the material model is
obtained through a Pareto analysis (see Section 3.5).

###  3.4 Solver

To solve the problem in Eq. 34, we leverage a coordinate descent algorithm as
proposed by Friedman et al. (2010), which is implemented in the Lasso
subroutine of the open-source Python library sklearn. The subroutine takes
$\bar{\mathbf{A}}$ and $\bar{\bm{b}}$ as well as the value of $\lambda_{p}$ as
input and returns the sparse solution vector $\bar{\bm{\theta}}$ whose entries
are constrained to be greater than or equal to zero. The intercept in the
Lasso subroutine is turned off and the maximum number of iterations is set to
$10^{4}$.

Figure 1: Qualitative illustration of the Pareto analysis for selecting an
appropriate value of $\lambda_{p}$.

###  3.5 Pareto analysis

The solution of the problem in Eq. 34 is highly dependent on the choice of
$\lambda_{p}$. As qualitatively illustrated in Fig. 1, a small value of
$\lambda_{p}$ yields a dense solution vector and therefore a complex
expression for the material model that is in general associated with a high
fitting accuracy, i.e., a small mean squared error. Increasing the value of
$\lambda_{p}$ results in a sparser solution vector and thus a simpler
expression for the material model at the cost of a reduced fitting accuracy
with a larger mean squared error. In the limit of very large values of
$\lambda_{p}$, the optimization problem returns a solution vector with only
zero entries and the mean squared error saturates.

To strike a balance between the accuracy and the complexity of the material
model, a Pareto analysis is applied to the problem. To this end, the problem
in Eq. 34 is solved multiple times using a set of different $\lambda_{p}$
values, leading to a solution vector $\bar{\bm{\theta}}$ for every choice of
$\lambda_{p}$. Here, we choose $41$ values of $\lambda_{p}$ that are evenly
distributed on a logarithmic scale between $10^{-2}$ and $10^{2}$. Our
objective now is to select from all solutions the one that exhibits an optimal
combination of high fitting accuracy and high sparsity. To this end, the mean
squared error and the $L_{1}$-norm of the solution vector are computed for
each solution. The minimum and maximum mean squared errors among all solutions
are denoted as $\overline{\mathrm{MSE}}_{\mathrm{min}}$ and
$\overline{\mathrm{MSE}}_{\mathrm{max}}$, respectively. Further, we define a
threshold mean squared error (see Fig. 1) as

$\overline{\mathrm{MSE}}_{\mathrm{th}}=\overline{\mathrm{MSE}}_{\mathrm{min}}+r_{\lambda}(\overline{\mathrm{MSE}}_{\mathrm{max}}-\overline{\mathrm{MSE}}_{\mathrm{min}}),$
(35)

where $r_{\lambda}$ is a small positive scalar, $0<r_{\lambda}\ll 1$. All
solutions for which the mean squared errors are greater than this threshold
are expected to exhibit a low fitting accuracy and are thus discarded. From
the remaining solutions, i.e., those for which the mean squared error is below
$\overline{\mathrm{MSE}}_{\mathrm{th}}$, we select the sparsest one, i.e., the
solution with the lowest value of the $L_{1}$-norm term (see Fig. 1). In this
way, a sparse solution with a low mean squared error is obtained, whereby the
amount of sparsity is dictated by the user-defined hyperparameter
$r_{\lambda}$. In this work we set $r_{\lambda}=0.02$ and keep it constant
throughout the analyses.

###  3.6 Thresholding

The solution vector $\bar{\bm{\theta}}$ selected in Section 3.5 may contain
values that are close to zero. As these material parameters and their
corresponding features have a vanishing effect on the model response, they are
set to zero. Specifically, we define a threshold value
$\bar{\theta}_{\text{th}}$ (here chosen as $\bar{\theta}_{\text{th}}=0.01$)
and we set $\bar{\theta}_{i}=0$ if $\bar{\theta}_{i}<\bar{\theta}_{\text{th}}$
for all $i$.

###  3.7 Feature clustering

Due to the fine discretization of the chosen range of $\alpha_{i}$ values (see
Section 2), the feature library contains Ogden features with very similar
values of $\alpha_{i}$, which in turn results in a high correlation between
these features. While the solution vector obtained by solving the optimization
problem with Lasso regularization already exhibits a high degree of sparsity,
the sparsity and thus interpretability of the material model can be further
improved by grouping together highly correlated features in the solution. To
this end, we leverage a clustering algorithm that automatically groups
together similar $\alpha_{i}$ values, as qualitatively illustrated in Fig. 2.
For a similar clustering strategy that was used to group together similar
Maxwell elements in viscoelastic constitutive laws discovered with EUCLID, the
interested reader is referred to Marino et al. (2023).

Figure 2: Qualitative illustration of the feature clustering.

In this paper, we adopt the clustering method DBSCAN (Density-Based Spatial
Clustering of Applications with Noise) (Ester et al., 1996) as implemented in
the Python library sklearn. DBSCAN assigns each $\alpha_{i}$ whose
corresponding $\\{\bar{\theta}_{\lambda}\\}_{i}$ is non-zero to a cluster.
This is done in such a way that for each $\alpha_{i}$ in a cluster there
exists at least one value $\alpha_{j}$ (with $i\neq j$) in the same cluster
such that their distance $|\alpha_{i}-\alpha_{j}|$ is smaller than or equal to
a predefined value, here chosen as 0.01. The advantage over the k-means
clustering algorithm used in Marino et al. (2023) is that the needed number of
clusters is obtained automatically. After grouping the $\alpha_{i}$ values
into a set of clusters, an averaged value $\bar{\alpha}$ is computed for each
cluster as

$\bar{\alpha}=\frac{\sum_{k}\alpha_{k}\\{\bar{\theta}_{\lambda}\\}_{k}}{\sum_{k}\\{\bar{\theta}_{\lambda}\\}_{k}},$
(36)

where $k$ ranges over all elements in the considered cluster. Afterwards,
$\bar{\alpha}$ is rounded to the closest original value of $\alpha_{i}$
contained in the material model library (see Section 2). Thus, we identify for
each cluster one value of $\alpha_{i}$ for which the corresponding coefficient
in $\bar{\bm{\theta}}_{\lambda}$ is non-zero, while all other entries in
$\bar{\bm{\theta}}_{\lambda}$ are set to zero.

###  3.8 Final regression without regularization

By solving the sparse regression problem and applying the clustering method, a
small number of relevant features are automatically selected from the
originally large library of candidate features. In a final step, all features
corresponding to parameters that have been identified as zero are disregarded
from the library by removing the corresponding columns from the matrix
$\bar{\mathbf{A}}$, thus obtaining a reduced system matrix
$\bar{\mathbf{A}}_{\text{red}}$. Afterwards, the reduced regression problem is
solved for the remaining non-zero parameters $\bar{\bm{\theta}}_{\text{red}}$,
but this time without introducing any regularization term (i.e., we solve Eq.
34 with $\lambda_{p}=0$).

The motivation behind this last step is twofold. First, while the
regularization term in Eq. 34 forces unneeded parameters to zero, it
simultaneously shrinks the absolute value of the relevant parameters. Solving
the regression problem without regularization for the small number of selected
features leads to material parameters that are not affected by such shrinkage.
And second, after applying the clustering method, most entries in
$\bar{\bm{\theta}}_{\lambda}$ have been identified as zero, while the
remaining non-zero parameters in $\bar{\bm{\theta}}_{\lambda}$ are still
unknown. The final regression problem without regularization is used to
calibrate those non-zero parameters. Finally, the original parameters
$\bm{\theta}$ can be recovered from the scaled parameters $\bar{\bm{\theta}}$
using Eq. 32.

##  4 Numerical verification

###  4.1 Synthetical data

Before applying the proposed method for material model discovery to
experimental data, we verify the method on synthetically generated data
without and with artificially added noise. Such data are obtained by assuming
different hyperelastic material models and generating labeled data pairs under
the assumption of uniaxial compression/tension and simple torsion.
Specifically, we consider the four material models whose strain energy density
functions are provided in Table 1:

  * 1.

first-order Mooney-Rivlin model (MR1)

  * 2.

one-term Ogden model (O1)

  * 3.

two-term Ogden model (O2)

  * 4.

combination of second-order Mooney-Rivlin model and one-term Ogden model
(MR2O1)

The chosen material models are used to generate $n_{\text{UT}}=60$ pairs of
data for uniaxial compression/tension, i.e.,
$(\lambda_{\text{UT}}^{(l)},P_{11}^{(l)})$, where the values of
$\lambda_{\text{UT}}^{(l)}$ are equally spaced between
$\lambda_{\text{UT}}^{(1)}=0.7$ and
$\lambda_{\text{UT}}^{(n_{\text{UT}})}=1.3$, and $n_{\text{ST}}=60$ pairs of
data for simple torsion, i.e., $(\tilde{\psi}^{(l)},\tau^{(l)})$, where the
values of $\tilde{\psi}^{(l)}$ are equally spaced between
$\tilde{\psi}^{(1)}=-1$ and $\tilde{\psi}^{(n_{\text{ST}})}=1$.

To emulate real experimental data, the generated data are perturbed by
artificial noise by considering noisy data pairs
$(\lambda_{\text{UT}}^{(l)},P_{11}^{(l)}+\varepsilon^{(l)}_{\text{UT}})$ and
$(\tilde{\psi}^{(l)},\tau^{(l)}+\varepsilon^{(l)}_{\text{ST}})$, where
$\varepsilon^{(l)}_{\text{UT}}\sim\mathcal{N}\left(0,\sigma^{2}\right)$ and
$\varepsilon^{(l)}_{\text{ST}}\sim\mathcal{N}\left(0,\sigma^{2}\right)$ is
random noise that is drawn independently from a Gaussian distribution with
zero mean and standard deviation $\sigma$. Here, we consider (beside the
noiseless case) two different levels of noise, i.e. $\sigma\in\\{0\
\text{Pa},5\ \text{Pa},10\ \text{Pa}\\}$.

###  4.2 Results and discussion

The synthetically generated data serve as input for the proposed algorithm for
material model discovery. Testing the method on synthetical data has the
advantage that the discovered material models can be compared to the ”ground
truth” (in machine learning jargon), i.e., to the material models that were
assumed during data generation. Table 1 shows the expressions of the
discovered strain energy density functions in comparison with the true ones.
It is observed that the correct mathematical form of the strain energy density
function is recovered for each test case. Table 1 further provides the mean
squared errors that quantify the mismatch between the data and the predictions
of the discovered models. As expected, the $\mathrm{MSE}$ increases for
increasing level of noise.

Table 1: Strain energy density functions of the true and discovered material
models and MSE of predictions.

Benchmarks$\frac{\int}{\int}$ |  Strain energy density $\tilde{W}$ $[\text{Pa}]$ | MSE $[\text{Pa}^{2}]$ |  $\overline{\mathrm{MSE}}$ $[-]$
---|---|---|---
MR1 | Truth$\frac{\int}{\int}$ | $40.00\left(I_{1}-3\right)$ | $+$ | $20.00\left(I_{2}-3\right)$ | - | -
|  $\sigma=0$$\frac{\int}{\int}$ | $40.00\left(I_{1}-3\right)$ | $+$ | $20.00\left(I_{2}-3\right)$ | 0.0000 | 0.0000
|  $\sigma=5\ \text{Pa}$$\frac{\int}{\int}$ | $41.25\left(I_{1}-3\right)$ | $+$ | $17.71\left(I_{2}-3\right)$ | 13.4120 | 0.0022
|  $\sigma=10\ \text{Pa}$$\frac{\int}{\int}$ | $42.51\left(I_{1}-3\right)$ | $+$ | $15.42\left(I_{2}-3\right)$ | 53.6481 | 0.0090
O1 | Truth$\frac{\int}{\int}$ | $2.00\left(\lambda_{1}^{-10.00}+\lambda_{2}^{-10.00}+\lambda_{3}^{-10.00}-3\right)$ |  |  | - | -
|  $\sigma=0$$\frac{\int}{\int}$ | $1.97\left(\lambda_{1}^{-10.03}+\lambda_{2}^{-10.03}+\lambda_{3}^{-10.03}-3\right)$ |  |  | 0.1815 | 0.0000
|  $\sigma=5\ \text{Pa}$$\frac{\int}{\int}$ | $1.94\left(\lambda_{1}^{-10.05}+\lambda_{2}^{-10.05}+\lambda_{3}^{-10.05}-3\right)$ |  |  | 13.5305 | 0.0001
|  $\sigma=10\ \text{Pa}$$\frac{\int}{\int}$ | $1.91\left(\lambda_{1}^{-10.08}+\lambda_{2}^{-10.08}+\lambda_{3}^{-10.08}-3\right)$ |  |  | 53.8606 | 0.0006
O2 | Truth$\frac{\int}{\int}$ | $16.00\left(\lambda_{1}^{-5.00\phantom{0}{}}+\lambda_{2}^{-5.00\phantom{0}{}}+\lambda_{3}^{-5.00\phantom{0}{}}-3\right)$ | $+$ | $\phantom{0}{}8.00\left(\lambda_{1}^{5.00\phantom{-0}{}}+\lambda_{2}^{5.00\phantom{-0}{}}+\lambda_{3}^{5.00\phantom{-0}{}}-3\right)$ | - | -
|  $\sigma=0$$\frac{\int}{\int}$ | $13.50\left(\lambda_{1}^{-5.25\phantom{0}{}}+\lambda_{2}^{-5.25\phantom{0}{}}+\lambda_{3}^{-5.25\phantom{0}{}}-3\right)$ | $+$ | $11.79\left(\lambda_{1}^{4.43\phantom{-0}{}}+\lambda_{2}^{4.43\phantom{-0}{}}+\lambda_{3}^{4.43\phantom{-0}{}}-3\right)$ | 0.2936 | 0.0000
|  $\sigma=5\ \text{Pa}$$\frac{\int}{\int}$ | $13.29\left(\lambda_{1}^{-5.26\phantom{0}{}}+\lambda_{2}^{-5.26\phantom{0}{}}+\lambda_{3}^{-5.26\phantom{0}{}}-3\right)$ | $+$ | $11.79\left(\lambda_{1}^{4.44\phantom{-0}{}}+\lambda_{2}^{4.44\phantom{-0}{}}+\lambda_{3}^{4.44\phantom{-0}{}}-3\right)$ | 13.7349 | 0.0002
|  $\sigma=10\ \text{Pa}$$\frac{\int}{\int}$ | $13.16\left(\lambda_{1}^{-5.26\phantom{0}{}}+\lambda_{2}^{-5.26\phantom{0}{}}+\lambda_{3}^{-5.26\phantom{0}{}}-3\right)$ | $+$ | $11.75\left(\lambda_{1}^{4.45\phantom{-0}{}}+\lambda_{2}^{4.45\phantom{-0}{}}+\lambda_{3}^{4.45\phantom{-0}{}}-3\right)$ | 54.0784 | 0.0007
MR2O1 | Truth$\frac{\int}{\int}$ | $30.00\left(I_{2}-3\right)^{2}$ | $+$ | $\phantom{0}{}2.00\left(\lambda_{1}^{-10.00}+\lambda_{2}^{-10.00}+\lambda_{3}^{-10.00}-3\right)$ | - | -
|  $\sigma=0$$\frac{\int}{\int}$ | $29.25\left(I_{2}-3\right)^{2}$ | $+$ | $\phantom{0}{}1.97\left(\lambda_{1}^{-10.04}+\lambda_{2}^{-10.04}+\lambda_{3}^{-10.04}-3\right)$ | 0.3179 | 0.0000
|  $\sigma=5\ \text{Pa}$$\frac{\int}{\int}$ | $31.69\left(I_{2}-3\right)^{2}$ | $+$ | $\phantom{0}{}1.91\left(\lambda_{1}^{-10.07}+\lambda_{2}^{-10.07}+\lambda_{3}^{-10.07}-3\right)$ | 13.7564 | 0.0002
|  $\sigma=10\ \text{Pa}$$\frac{\int}{\int}$ | $34.16\left(I_{2}-3\right)^{2}$ | $+$ | $\phantom{0}{}1.86\left(\lambda_{1}^{-10.10}+\lambda_{2}^{-10.10}+\lambda_{3}^{-10.10}-3\right)$ | 53.9747 | 0.0007

For the benchmark case MR1, the material coefficients in the strain energy
density function are exactly recovered for the noiseless case. The deviation
between the true and identified parameters increases for increasing noise, as
expected. For the benchmark cases O1 and MR2O1, the true and discovered
parameters are in excellent agreement. However, the material parameters are
not exactly recovered even in the noiseless case. The reason for this is that
computing the averaged values $\bar{\alpha}$ after the feature clustering
(Section 3.7) cannot exactly recover the ground truth exponents in the Ogden
terms, which leads to small deviations in the identified parameters. For the
benchmark case O2, the true and identified parameters exhibit slightly larger
deviations. One explanation for these deviations could be that the model
library comprises many material models and many combinations of material
parameters that are equally suited to describe the material response. The
proposed algorithm does not guarantee to find exactly the same model and
parameters that were assumed during data generation, but may instead provide
an adequate surrogate model that is equally suited to predict the data. We
emphasise that, although the true and identified material parameters show some
deviations, the $\mathrm{MSE}$ for the benchmark case O2 is similar to the
$\mathrm{MSE}$ for the other benchmarks. Thus, the deviation in the model
parameters does not seem to seriously affect the fitting accuracy of the
model.

(a) MR1 under UT

(b) MR1 under ST

(c) O1 under UT

(d) O1 under ST

(e) O2 under UT

(f) O2 under ST

(g) MR2O1 under UT

(h) MR2O1 under ST

Figure 3: Stretch versus stress response under uniaxial compression/tension
and (normalized) twist versus torque response under simple torsion of the
discovered material models in comparison with the synthetical data (noise
level $\sigma=10\ \text{Pa}$).

Fig. 3 illustrates the material response of the discovered material models in
direct comparison with the synthetically generated data for the case with the
highest level of noise, revealing an excellent agreement.

In the following, we discuss in more detail some of the intermediate steps in
the sparsity promoting algorithm for material model discovery (see Steps 3-8
in Fig. 10) in the context of benchmark case MR2O1 with noise level
$\sigma=10\ \text{Pa}$. Fig. 4 illustrates the Pareto analysis, which is
leveraged for the choice of the hyperparameter $\lambda_{p}$ in the sparse
regression problem. It can be seen that the mean squared error
$\overline{\mathrm{MSE}}$ increases and the $L_{1}$-norm of the material
parameter vector $\|\bar{\bm{\theta}}\|_{1}$ decreases for increasing values
of $\lambda_{p}$. For values of $\lambda_{p}$ exceeding approximately $1$, the
material parameter vector obtained from the solution of the sparse regression
problem contains all zeros, so that its $L_{1}$-norm remains constant at zero
and the corresponding $\overline{\mathrm{MSE}}$ saturates. Using the method
described in Section 3.5, the appropriate value of $\lambda_{p}$ is chosen
such that the selected model exhibits both a high fitting accuracy and a low
model complexity. Note also that for values of $\lambda_{p}$ lower than the
chosen one, the $\overline{\mathrm{MSE}}$ is very weakly sensitive to the
specific choice of $\lambda_{p}$.

Figure 4: Illustration of the Pareto analysis for benchmark case MR2O1 (noise
level $\sigma=10\ \text{Pa}$).

In Fig. 5, the effect of the clustering algorithm (see Section 3.7) is
illustrated. Before applying the clustering method, many values of
$\alpha_{i}$ are zero while some non-zero values agglomerate in a small region
near the ground truth value (equal to 10, see Table 1). The clustering
algorithm identifies these values as one cluster. Afterwards, the values of
$\alpha_{i}$ belonging to the cluster are averaged, such that the discovered
material model exhibits only one Ogden term as the ground truth model, and the
feature coefficients are calibrated, as illustrated in Section 3.8.

(a) Before clustering

(b) After clustering and final regression without regularization

Figure 5: Illustration of the feature clustering for benchmark case MR2O1
(noise level $\sigma=10\ \text{Pa}$).

##  5 Experimental validation

###  5.1 Experimental data

In the context of this work, we use multi-modal large-strain mechanical
testing data from 81 cylindrical specimens extracted from three different
human brains. The data acquisition and preprocessing are briefly described in
the following. For details, the reader is referred to Hinrichsen et al.
(2022).

81 cylindrical specimens from different regions of the brain are prepared.
Each specimen is labeled by a brain number and a specimen number as
$\text{HBE}\\_\langle\text{brain~{}no.}\rangle\\_\langle\text{specimen~{}no.}\rangle$.
The specimens have a radius of approximately $r_{\text{out}}\approx 4\
\text{mm}$. The heights of the specimens and the region that they stem from
are recorded in Table 2.

Fig. 6 shows one of the human brains and an exemplary cylindrical sample
before and after mounting it to the rheometer for the mechanical testing. Each
loading mode (compression/tension and torsional shear) consists of three
cycles. Here, we use the data from the third cycle representing the
preconditioned material response as brain tissue was found to show substantial
pre-conditioning behavior (Budday et al., 2017). We limit ourselves to
hyperelastic material models (neglecting poro- and viscoelastic effects). The
strain rates needed to obtain a purely hyperelastic response are (in theory
infinitely) high or low and therefore not feasible in actual experiments.
Thus, the experimental data obtained from the rheometer shows a considerable
hysteresis (see Fig. 6) and we extract the hyperelastic response through a
preprocessing procedure. First, a moving average as well as a low pass filter
is applied to reduce high frequency noises. Afterwards, the hyperelastic
response is approximated by the averaged loading and unloading curves (see
also Budday et al. (2017)). The raw data and the extracted hyperelastic
response are illustrated in Fig. 6. To reduce the dimensionality of the data
while still preserving the characteristic shape of the curve, we reduce the
points per deformation mode to 60.

This finally leads to a number of $n_{\text{UT}}=60$ data pairs
$(\lambda_{\text{UT}}^{(l)},P_{11}^{(l)})$ and $n_{\text{ST}}=60$ data pairs
$(\tilde{\psi}^{(l)},\tau^{(l)})$ acquired from the mechanical tests under
uniaxial compression/tension and simple torsion, respectively. The values of
$\lambda_{\text{UT}}^{(l)}$ are uniformly distributed between
$\lambda_{\text{UT}}^{(1)}=0.85$ and
$\lambda_{\text{UT}}^{(n_{\text{UT}})}=1.15$ and the values of
$\tilde{\psi}^{(l)}$ are uniformly distributed between
$\tilde{\psi}^{(1)}=-0.3$ and $\tilde{\psi}^{(n_{\text{ST}})}=0.3$.

Figure 6: Experimental setup and data preprocessing. Cylindrical samples (b)
are extracted from multiple locations within the human brains (a) and then
mounted to the rheometer (c). The raw data from the rheometer (blue) is
preprocessed to obtain the hyperelastic response (orange) (see (d) and (e)).
This is exemplarily shown for a specimen from the motor cortex (HBE_1_8).
First, the noise is reduced using a moving average as well as a lowpass
filter. Subsequently, the hyperelastic response is obtained by averaging the
loading and unloading curve.

###  5.2 Results and discussion

The proposed algorithm for material model discovery is applied to all
experimentally acquired data sets, keeping all hyperparameters and algorithmic
settings constant. The resulting material parameters are listed in Table 4 and
Table 5. Many material parameters are identified as zero, meaning that the
discovered expressions of the strain energy density functions are concise and
interpretable. The material model discovered for most of the specimens (i.e.,
about $60\%$ of the specimens) is the one-term Ogden model, which has also
previously been used to fit human brain tissue data (Budday et al., 2017,
2020; Hinrichsen et al., 2022). Further models that were identified are the
two-term Ogden model (about $32\%$ of the specimens), and combinations of the
one-term Ogden model with a Mooney-Rivlin term (about $7\%$ of the specimens).

To measure the ability of the discovered material models to fit the given
experimental data, the mean squared errors $\mathrm{MSE}$ and
$\overline{\mathrm{MSE}}$ are provided in Table 4 and Table 5. The scaled mean
squared error $\overline{\mathrm{MSE}}$ ranges from a minimum value of $0.01$
to a maximum value of $0.71$, and the average $\overline{\mathrm{MSE}}$ over
all experimental data sets is $0.17$.

(a) HBE_01_26 under UT

(b) HBE_01_26 under ST

Figure 7: Stress versus stretch response under uniaxial compression/tension
and torque versus (normalized) twist response under simple torsion of a
discovered material model in comparison with the experimental data. The mean
squared error $\overline{\mathrm{MSE}}=0.0091$ is the smallest in the data
set.

Fig. 7 shows the stress versus stretch response under uniaxial
compression/tension and torque versus normalized twist response under simple
torsion of a discovered material model for which the value of
$\overline{\mathrm{MSE}}$ is small in comparison to the other experiments. It
is evident that the discovered material model is very well suited to describe
the given data. Fig. 8 shows the same plots for discovered material models for
which the values of $\overline{\mathrm{MSE}}$ are close to the average value
of $\overline{\mathrm{MSE}}$ over all experiments. The models are in
qualitative agreement with the data. However, especially for the response
under uniaxial compression/tension, the models are not capable to precisely
match the given data. These deviations could possibly be attributed to
modeling assumptions that do not hold true during the mechanical experiments.
For example, the geometries of some of the specimens were not ideally
cylindrical, because the softer specimens deformed under their own weight.
Further, it was assumed throughout the derivations in this paper that the
strain state in the specimens is homogeneous under uniaxial
compression/tension, and it was assumed that the specimens are fixed in
longitudinal direction at the boundaries, but free to move in the transversal
directions. During the experiments, however, both ends of the specimens were
glued, meaning that the cross-sectional area of the specimens remained
unchanged at the glued boundaries. Therefore, the sides of the cylindrical
specimens bulged into a barrel shape, resulting in non-homogeneous strain
states in the specimens. The constraint that the specimens are fixed in both
longitudinal and transversal directions may add additional stiffness to the
experimentally measured stress versus stretch response, which could explain
the underestimated stress by the discovered material models during uniaxial
compression/tension in Fig. 8. Another source for the disagreement between the
discovered models and the data could be that only isotropic material behavior
has been included in the material model library. Thus, possible effects from
material anisotropy in the experimental data cannot be predicted by the
discovered models.

(a) HBE_01_29 under UT

(b) HBE_01_29 under ST

(c) HBE_02_15 under UT

(d) HBE_02_15 under ST

(e) HBE_03_15 under UT

(f) HBE_03_15 under ST

Figure 8: Stress versus stretch response under uniaxial compression/tension
and torque versus (normalized) twist response under simple torsion of a
discovered material model in comparison with the experimental data. The mean
squared errors ($\overline{\mathrm{MSE}}=0.1657$,
$\overline{\mathrm{MSE}}=0.1671$, $\overline{\mathrm{MSE}}=0.1664$ for
HBE_01_29, HBE_02_15, HBE_03_15, respectively) are close to the average
$\overline{\mathrm{MSE}}$ in the data set.

(a) HBE_01_24 under UT

(b) HBE_01_24 under ST

(c) HBE_01_30 under UT

(d) HBE_01_30 under ST

Figure 9: Stress versus stretch response under uniaxial compression/tension
and torque versus (normalized) twist response under simple torsion of
discovered material models in comparison with the experimental data. The mean
squared errors ($\overline{\mathrm{MSE}}=0.7064$,
$\overline{\mathrm{MSE}}=0.6107$ for HBE_01_24, HBE_01_30, respectively) are
the largest in the data set.

Finally, Fig. 9 shows the stretch versus stress and normalized twist versus
torque plots of discovered material models for which the value of
$\overline{\mathrm{MSE}}$ is large in comparison to the other experiments. As
the large mean squared error indicates, the mismatch between the discovered
model and the data is significant. Such large discrepancies are exceptional
for the considered experimental data sets (see Table 4 and Table 5). For
cyclic shear tests, the reason for the large deviation between the model and
the data could be erroneous experimental measurements: the two shown samples
were not loaded symmetrically during cyclic shear tests but were already
twisted before the shear test started. This can happen when a torque is
induced during axial loading (for example due to inhomogeneous microstructure
of the sample). In this case, the rheometer applies a deformation until the
torque goes back to zero. The non-monotonic stress versus stretch response
under uniaxial compression/tension in Fig. 9(a) could be a consequence of a
low signal-to-noise ratio since this sample is very soft. The non-increasing
constant stress plateau under uniaxial compression/tension in Fig. 9(c) could
indicate a failure of the glue that connects the specimen with the testing
machine.

##  6 Conclusion

EUCLID, a recently proposed computational framework for automatically
discovering constitutive models as interpretable symbolic mathematical
expressions, was formulated in this work in a supervised setting and validated
on the mechanical testing data of human brain specimens. It was shown that the
proposed model selection strategy, which comprises sparse regression and
feature clustering, is able to automatically discover interpretable
constitutive models with satisfactory fitting accuracy to the given uniaxial
compression/tension and simple torsion data. The most commonly discovered
model was the one-term Ogden model, which has also previously been widely used
to fit human brain tissue data. The presented results constitute the first
experimental validation of EUCLID. However, several aspects of the proposed
method may be improved or extended in future studies. For example, the method
proposed in this work is informed by labeled data pairs due to the limited
access to full-field displacement data for human brain tissues. The
experimental validation of EUCLID driven by unlabeled data, i.e., full-field
displacement and net reaction force data (see Flaschel (2023)), is the subject
of ongoing research. Moreover, we restricted the model library in which EUCLID
searches for a model to the class of isotropic hyperelastic models. As the
tissue response is in fact not just hyperelastic but affected by dissipative
effects as well as possibly not perfectly isotropic, an important future study
could be to further generalize the model library, such that EUCLID can not
only select a constitutive model for the human brain tissue within a
predefined class of models, but further distinguish between different modeling
classes (such as hyperelasticity, viscoelasticity, elastoplasticity), as shown
by Flaschel et al. (2023) for synthetical data, as well as automatically
detect possible anisotropies and identify the corresponding parameters.

##  Appendix A Derivatives

###  A.1 Derivatives of the strain invariants

The derivatives of the strain invariants with respect to the deformation
gradient are

$\frac{\partial I_{1}}{\partial\mathbf{F}}=2\mathbf{F},\quad\frac{\partial
I_{2}}{\partial\mathbf{F}}=2I_{1}\mathbf{F}-2\mathbf{F}\mathbf{C}.$ (37)

###  A.2 Derivatives of the principal stretches

The derivatives of the principal directions depend on the eigenvectors of
$\mathbf{C}=\mathbf{F}^{T}\mathbf{F}$, denoted by $\bm{N}_{i}$, and the
eigenvectors of $\bm{b}=\mathbf{F}\mathbf{F}^{T}$, denoted by $\bm{n}_{i}$.
Under the assumption
$\lambda_{1}\neq\lambda_{2}\neq\lambda_{3}\neq\lambda_{1}$, it is (see
Holzapfel (2000))

$\frac{\partial\lambda_{1}}{\partial\mathbf{F}}=\bm{n}_{1}\otimes\bm{N}_{1},\quad\frac{\partial\lambda_{2}}{\partial\mathbf{F}}=\bm{n}_{2}\otimes\bm{N}_{2},\quad\frac{\partial\lambda_{3}}{\partial\mathbf{F}}=\bm{n}_{3}\otimes\bm{N}_{3}.$
(38)

###  A.3 Derivatives of the feature vectors

The derivatives of the feature vector $\bm{Q}_{I}$ with respect to the strain
invariants are

$\displaystyle\frac{\partial\bm{Q}_{I}}{\partial I_{1}}$
$\displaystyle=\left[m(I_{1}-3)^{m-1}(I_{2}-3)^{n-m}:n\in\\{1,\dots,N_{\text{Mooney}}\\},m\in\\{0,\dots,n\\}\right]^{T}\oplus\left[~{}0~{}\right],$
(39) $\displaystyle\frac{\partial\bm{Q}_{I}}{\partial I_{2}}$
$\displaystyle=\left[(n-m)(I_{1}-3)^{m}(I_{2}-3)^{n-m-1}:n\in\\{1,\dots,N_{\text{Mooney}}\\},m\in\\{0,\dots,n\\}\right]^{T}\oplus\left[~{}1/I_{2}~{}\right],$

which can be used to calculate the derivatives of the strain energy density
contribution $\tilde{W}_{I}$ with respect to the strain invariants, i.e.,

$\frac{\partial\tilde{W}_{I}}{\partial
I_{a}}=\bm{\theta}_{I}\cdot\frac{\partial\bm{Q}_{I}}{\partial I_{a}}.$ (40)

The components of the feature vector $\bm{Q}_{\lambda}$ are defined as

$\left\\{Q_{\lambda}(\lambda_{1},\lambda_{2},\lambda_{3})\right\\}_{i}=\lambda_{1}^{\alpha_{i}}+\lambda_{2}^{\alpha_{i}}+\lambda_{3}^{\alpha_{i}}-3.$
(41)

Assuming incompressibility $\lambda_{1}\lambda_{2}\lambda_{3}\overset{!}{=}1$,
the components of the feature vector $\bm{Q}_{\lambda}$ can be written as

$\left\\{Q_{\lambda}(\lambda_{1},\lambda_{2})\right\\}_{i}=\lambda_{1}^{\alpha_{i}}+\lambda_{2}^{\alpha_{i}}+\lambda_{1}^{-\alpha_{i}}\lambda_{2}^{-\alpha_{i}}-3.$
(42)

Thus, the derivatives of the components of the feature vector
$\bm{Q}_{\lambda}$ with respect to the principal stretches are

$\displaystyle\frac{\partial\left\\{Q_{\lambda}(\lambda_{1},\lambda_{2})\right\\}_{i}}{\partial\lambda_{1}}$
$\displaystyle=\alpha_{i}\left(\lambda_{1}^{\alpha_{i}-1}-\lambda_{1}^{-\alpha_{i}-1}\lambda_{2}^{-\alpha_{i}}\right),$
(43)
$\displaystyle\frac{\partial\left\\{Q_{\lambda}(\lambda_{1},\lambda_{2})\right\\}_{i}}{\partial\lambda_{2}}$
$\displaystyle=\alpha_{i}\left(\lambda_{2}^{\alpha_{i}-1}-\lambda_{1}^{-\alpha_{i}}\lambda_{2}^{-\alpha_{i}-1}\right),$

which can be used to calculate the derivatives of the strain energy density
contribution $\tilde{W}_{\lambda}$ with respect to the principal stretches,
i.e.,

$\frac{\partial\tilde{W}_{\lambda}}{\partial\lambda_{i}}=\bm{\theta}_{\lambda}\cdot\frac{\partial\bm{Q}_{\lambda}}{\partial\lambda_{i}}.$
(44)

###  A.4 Derivatives for uniaxial compression and tension

In the following, we seek to simplify the expression for
$P_{11}(\lambda_{\text{UT}};\bm{\theta})$ (see Eq. 16) by leveraging the
characteristics of the kinematic state under uniaxial compression and tension.
Considering such kinematic state, the deformation gradient and the resulting
right Cauchy-Green strain tensor read

$\displaystyle\mathbf{F}=\begin{bmatrix}\lambda_{\text{UT}}&0&0\\\
0&\frac{1}{\sqrt{\lambda_{\text{UT}}}}&0\\\
0&0&\frac{1}{\sqrt{\lambda_{\text{UT}}}}\\\
\end{bmatrix},\qquad\mathbf{C}=\begin{bmatrix}\lambda_{\text{UT}}^{2}&0&0\\\
0&\frac{1}{\lambda_{\text{UT}}}&0\\\ 0&0&\frac{1}{\lambda_{\text{UT}}}\\\
\end{bmatrix}.$ (45)

Thus, the strain invariants simplify to

$I_{1}=\lambda_{\text{UT}}^{2}+\frac{2}{\lambda_{\text{UT}}},\quad
I_{2}=2\lambda_{\text{UT}}+\frac{1}{\lambda_{\text{UT}}^{2}},$ (46)

and the principal stretches simplify to

$\lambda_{1}=\lambda_{\text{UT}},\quad\lambda_{2}=\lambda_{3}=\frac{1}{\sqrt{\lambda_{\text{UT}}}}.$
(47)

The simple structure of the kinematic state is used to simplify the feature
vector $\bm{Q}$. The feature vectors $\bm{Q}_{I}$ and $\bm{Q}_{\lambda}$ are

$\bm{Q}_{I}(\lambda_{\text{UT}})=\left[\left(\lambda_{\text{UT}}-3\right)^{m}\left(2\lambda_{\text{UT}}+\frac{1}{\lambda_{\text{UT}}^{2}}-3\right)^{n-m}:n\in\\{1,\dots,N_{\text{Mooney}}\\},m\in\\{0,\dots,n\\}\right]^{T}\oplus\left[\log\left(\frac{2}{3}\lambda_{\text{UT}}+\frac{1}{3\lambda_{\text{UT}}^{2}}\right)\right],$
(48)

and

$\left\\{Q_{\lambda}(\lambda_{\text{UT}})\right\\}_{i}=\lambda_{\text{UT}}^{\alpha_{i}}+2\left(\frac{1}{\sqrt{\lambda_{\text{UT}}}}\right)^{\alpha_{i}}-3,$
(49)

whose derivatives with respect to $\lambda_{\text{UT}}$ are trivial. We
finally obtain the following expression for the feature vector derivative in
Eq. 16

$\bm{Q}^{\prime}_{\text{UT}}(\lambda_{\text{UT}})=\frac{\partial\bm{Q}(\lambda_{\text{UT}})}{\partial\lambda_{\text{UT}}},$
(50)

which after substitution in Eq. 16 results in a simplified expression for
$P_{11}(\lambda_{\text{UT}};\bm{\theta})$.

###  A.5 Derivatives for simple torsion

In the following, we seek to simplify the expression for
$\tau(\tilde{\psi};\bm{\theta})$ (see Eq. 21) by leveraging the
characteristics of the kinematic state under simple torsion. Considering
simple torsion, the deformation gradient and the resulting right Cauchy-Green
strain tensor read

$\displaystyle\mathbf{F}=\begin{bmatrix}1&0&0\\\ 0&1&F_{\vartheta z}\\\
0&0&1\\\ \end{bmatrix},\qquad\mathbf{C}=\begin{bmatrix}1&0&0\\\
0&1&F_{\vartheta z}\\\ 0&F_{\vartheta z}&1+F_{\vartheta z}^{2}\\\
\end{bmatrix},$ (51)

where $F_{\vartheta z}=\rho\tilde{\psi}$. Thus, the strain invariants simplify
to

$I_{1}=I_{2}=F_{\vartheta z}^{2}+3.$ (52)

Defining the auxiliary variable $\bar{C}=1+\frac{1}{2}F_{\vartheta z}^{2}$,
the eigenvalues of $\mathbf{C}$ (which are the squares of the principal
stretches) read

$\lambda_{1}^{2}=\bar{C}-\sqrt{\bar{C}^{2}-1},\quad\lambda_{2}^{2}=1,\quad\lambda_{3}^{2}=\bar{C}+\sqrt{\bar{C}^{2}-1},$
(53)

which results in the principal stretches

$\lambda_{1}=\sqrt{\bar{C}-\sqrt{\bar{C}^{2}-1}},\quad\lambda_{2}=1,\quad\lambda_{3}=\sqrt{\bar{C}+\sqrt{\bar{C}^{2}-1}}.$
(54)

The simple structure of the kinematic state is used to simplify the feature
vector $\bm{Q}$. The feature vector $\bm{Q}_{I}$ is

$\bm{Q}_{I}(F_{\vartheta z})=\left[F_{\vartheta
z}^{2n}:n\in\\{1,\dots,N_{\text{Mooney}}\\},m\in\\{0,\dots,n\\}\right]^{T}\oplus\left[\log\left(F_{\vartheta
z}^{2}/3+1\right)\right],$ (55)

whose derivative with respect to $F_{\vartheta z}$ is trivial, and the feature
vector $\bm{Q}_{\lambda}$ is

$\left\\{Q_{\lambda}(\bar{C})\right\\}_{i}=\left(\sqrt{\bar{C}-\sqrt{\bar{C}^{2}-1}}\right)^{\alpha_{i}}+\left(\sqrt{\bar{C}+\sqrt{\bar{C}^{2}-1}}\right)^{\alpha_{i}}-2,$
(56)

whose derivative with respect to $F_{\vartheta z}$ can be computed using the
chain rule

$\frac{\partial\bm{Q}_{\lambda}}{\partial F_{\vartheta
z}}=\frac{\partial\bm{Q}_{\lambda}}{\partial\bar{C}}\frac{\partial\bar{C}}{\partial
F_{\vartheta z}}=\frac{\partial\bm{Q}_{\lambda}}{\partial\bar{C}}F_{\vartheta
z}.$ (57)

Thus, the derivative $\frac{\partial\bm{Q}}{\partial F_{\vartheta z}}$ is
known and it depends on $\rho$ and $\tilde{\psi}$ through
$\bar{C}=1+\frac{1}{2}F_{\vartheta z}^{2}$ and $F_{\vartheta
z}=\rho\tilde{\psi}$. With reference to Eq. 21, it is

$\bm{Q}^{\prime}_{\text{ST}}\left(\tilde{\psi}\right)=\int_{0}^{1}\
2\pi\rho^{2}\frac{\partial\bm{Q}}{\partial F_{\vartheta z}}\ \mathrm{d}\rho,$
(58)

where the integral over $\rho$ can be computed through numerical quadrature.

##  Appendix B Algorithm

Fig. 10 provides a step-by-step description of the proposed algorithm for
material model discovery.

Step 1: Compute $\mathbf{A}_{\mathrm{UT}}$, $\bm{b}_{\mathrm{UT}}$ and
$\mathbf{A}_{\mathrm{ST}}$, $\bm{b}_{\mathrm{ST}}$, and assemble $\mathbf{A}$
and $\bm{b}$Step 2: Apply feature scaling to obtain $\bar{\mathbf{A}}$ and
$\bar{\bm{b}}$Step 3: Solve the Lasso problem for a set of distinct values of
$\lambda_{p}$Step 4: Discard all solutions from Step 3 with
$\overline{\mathrm{MSE}}(\bar{\bm{\theta}})>\overline{\mathrm{MSE}}_{\mathrm{th}}$Step
5: Select the solution remaining from Step 4 with the smallest value of
$\|\bar{\bm{\theta}}\|_{1}$Step 6: For all $i$, if
$\bar{\theta}_{i}<\bar{\theta}_{\text{th}}$, set $\bar{\theta}_{i}=0$Step 7:
Use the clustering method to group together similar values of $\alpha_{i}$Step
8: Apply regression without regularization to calibrate the non-zero entries
in $\bar{\bm{\theta}}$Step 9: Compute $\bm{\theta}$ from $\bar{\bm{\theta}}$
Figure 10: Step-by-step description of the sparsity promoting algorithm for
material model discovery.

##  Appendix C Specimen information

The corresponding heights and regions of the different specimens are recorded
in Table 2, while the radius of each specimen is $r_{\text{out}}\approx
4\text{mm}$. Table 3 lists the abbreviations of the brain regions and their
meaning (see Hinrichsen et al. (2022)).

Table 2: Heights and regions of the brain tissue specimens. HBE | Region |  $h$ [mm] | HBE | Region |  $h$ [mm] | HBE | Region |  $h$ [mm] | HBE | Region |  $h$ [mm]
---|---|---|---|---|---|---|---|---|---|---|---
01_03 | C | 4.89159 | 02_01 | CR | 2.9083 | 03_01 | BS | 6.1803 | 03_21 | CC | 4.32991
01_04 | BG | 4.44005 | 02_02 | C | 4.36397 | 03_02 | BS | 5.17463 | 03_22 | CC | 3.91166
01_05 | CR | 4.64004 | 02_03 | C | 3.56882 | 03_03 | BS | 6.52776 | 03_23 | CR | 5.91647
01_06 | CR | 4.19198 | 02_04 | CR | 3.70545 | 03_04 | Am | 4.04604 | 03_25 | C | 4.7642
01_07 | BG | 3.86355 | 02_05 | CC | 3.71248 | 03_05 | M | 5.39193 | 03_26 | Am | 4.79861
01_08 | C | 5.1278 | 02_06 | CC | 4.27716 | 03_06 | M | 5.20739 | 03_27 | BG | 5.46595
01_09 | BG | 3.6299 | 02_07 | CC | 3.8597 | 03_07 | CC | 3.7975 | 03_28 | BG | 5.80823
01_10 | CR | 3.79468 | 02_08 | CC | 4.0629 | 03_08 | CC | 4.29526 | 03_29 | CR | 5.61761
01_12 | CC | 4.29779 | 02_09 | M | 5.16654 | 03_09 | BG | 5.64718 | 03_30 | CR | 4.97806
01_13 | CR | 3.89796 | 02_10 | BG | 4.00143 | 03_10 | BG | 4.96097 | 03_31 | C | 4.68306
01_14 | C | 4.59097 | 02_13 | M | 5.28094 | 03_11 | CC | 5.15737 | 03_32 | CR | 4.01945
01_21 | C | 4.42721 | 02_14 | BG | 4.84645 | 03_12 | BG | 5.96805 | 03_33 | CR | 6.07937
01_22 | CR | 3.7167 | 02_15 | BG | 4.69894 | 03_13 | CC | 3.27287 | 03_34 | CR | 4.66727
01_23 | C | 3.78048 | 02_16 | CR | 4.96181 | 03_14 | CB | 6.77129 | 03_35 | CR | 5.81572
01_24 | CR | 2.79804 | 02_17 | C | 4.63 | 03_15 | CB | 6.74769 | 03_36 | C | 3.8775
01_25 | BS | 3.61422 | 02_18 | BG | 3.91647 | 03_16 | CB | 6.40546 | 03_37 | M | 5.80837
01_26 | BS | 4.11383 | 02_19 | C | 4.79345 | 03_17 | CB | 5.18112 | 03_38 | M | 6.11145
01_27 | CB | 4.9045 | 02_20 | BS | 5.41687 | 03_18 | CR | 6.19653 | 03_39 | C | 4.75102
01_28 | CB | 4.6799 | 02_21 | BS | 4.36067 | 03_19 | C | 4.86192 | 03_40 | CR | 5.26147
01_29 | M | 3.121 | 02_22 | CB | 5.37975 | 03_20 | BG | 6.97399 | 03_41 | C | 4.96854
01_30 | M | 4.74689 |  |  |  |  |  |  |  |  |
Table 3: Abbreviations and full names of the brain regions. Am |  Amygdala |  BG |  Basal ganglia |  BS |  Brain stem |  C |  Cortex
---|---|---|---|---|---|---|---
CB |  Cerebellum |  CC |  Corpus callosum |  CR |  Corona radiata |  M |  Midbrain

##  Appendix D Material parameters

The discovered material models and the corresponding parameters are reported
in Table 4 and Table 5. The tables further record different measures for the
fitting accuracy of the discovered models as well as the computational time
needed for the discovery process.

Table 4: Material parameters, mean squared errors and time to solution for the
material model discovery from experimental data - brain no. 01 and 02. Note
that model parameters that are zero throughout all experiments are not shown
in the table. All material parameters are in Pa, except for the dimensionless
parameters $\alpha_{i}$.

| Material parameters$\frac{\int}{\int}$ | MSE | $\overline{\mathrm{MSE}}$ | Time
---|---|---|---|---
HBE | $\mu_{1}$ | $\alpha_{1}$ | $\mu_{2}$ | $\alpha_{2}$ | $I_{2}-3$ | $\left(I_{2}-3\right)^{2}$ | $\left(I_{2}-3\right)\left(I_{1}-3\right)$ | $\left(I_{1}-3\right)^{2}$ | $[\text{Pa}^{2}]$ | $[-]$ | $[s]$
01_03 | 29.85 | -38.84 | - | - | - | - | - | - | 1364.68 | 0.26 | 688
01_04 | 125.13 | -21.14 | - | - | - | - | - | - | 260.16 | 0.08 | 713
01_05 | 0.13 | -68.28 | 6.64 | -33.86 | - | - | - | - | 373.01 | 0.51 | 718
01_06 | - | - | - | - | 33.50 | - | - | - | 189.60 | 0.50 | 638
01_07 | 7.38 | -49.08 | 69.29 | -25.57 | - | - | - | - | 2861.05 | 0.27 | 719
01_08 | 350.45 | -18.80 | - | - | - | - | - | - | 1327.16 | 0.07 | 741
01_09 | 122.55 | -25.75 | - | - | - | - | - | - | 878.70 | 0.14 | 715
01_10 | 1.80 | -57.89 | - | - | - | - | - | - | 1027.67 | 0.43 | 666
01_12 | 33.75 | -20.06 | - | - | - | - | - | - | 16.97 | 0.09 | 680
01_13 | 0.85 | -54.19 | 49.02 | -23.15 | - | - | - | - | 251.27 | 0.18 | 670
01_14 | 274.80 | -21.10 | - | - | - | - | - | - | 606.65 | 0.04 | 680
01_21 | 2.09 | -54.38 | 10.20 | -24.62 | - | - | - | - | 872.84 | 0.49 | 675
01_22 | 2.46 | -45.70 | 45.62 | -25.99 | - | - | - | - | 291.84 | 0.16 | 674
01_23 | 35.76 | -36.88 | - | - | - | - | - | - | 983.44 | 0.23 | 790
01_24 | 16.35 | -20.05 | - | - | - | - | - | - | 140.42 | 0.71 | 684
01_25 | 0.10 | -70.62 | 62.52 | -25.82 | - | - | - | - | 540.38 | 0.18 | 717
01_26 | 207.64 | -24.58 | - | - | - | 392.93 | - | - | 136.29 | 0.01 | 758
01_27 | 2.37 | -52.66 | 64.85 | -24.86 | - | - | - | - | 970.58 | 0.23 | 708
01_28 | 119.72 | -18.48 | - | - | - | 160.40 | - | - | 46.12 | 0.02 | 727
01_29 | 0.01 | -78.58 | 90.36 | -28.71 | - | - | - | - | 1138.84 | 0.17 | 781
01_30 | 55.97 | -28.04 | - | - | - | - | - | - | 3227.89 | 0.61 | 743
02_01 | 6.03 | -51.39 | - | - | - | - | - | - | 1499.17 | 0.35 | 816
02_02 | 1.37 | -51.25 | 76.88 | -22.85 | - | - | - | - | 275.26 | 0.11 | 816
02_03 | 56.66 | -34.54 | - | - | - | - | - | - | 1369.33 | 0.20 | 750
02_04 | 9.11 | -46.83 | - | - | - | - | - | - | 1014.31 | 0.32 | 767
02_05 | 1.06 | -51.48 | 16.30 | -21.06 | - | - | - | - | 64.85 | 0.20 | 700
02_06 | 0.37 | -46.55 | 16.32 | -24.76 | - | - | - | - | 26.25 | 0.18 | 686
02_07 | 53.93 | -18.72 | - | - | - | - | - | - | 20.69 | 0.05 | 710
02_08 | 0.49 | -50.67 | 42.52 | -10.72 | - | - | - | - | 42.58 | 0.18 | 730
02_09 | 134.96 | -25.83 | - | - | - | 280.13 | - | - | 211.21 | 0.03 | 763
02_10 | 7.25 | -47.26 | - | - | - | - | - | - | 808.74 | 0.34 | 692
02_13 | 0.94 | -51.11 | 72.32 | -25.43 | - | - | - | - | 272.57 | 0.10 | 746
02_14 | 92.21 | -24.97 | - | - | - | - | - | - | 166.90 | 0.06 | 766
02_15 | 2.21 | -52.14 | 58.47 | -26.20 | - | - | - | - | 591.73 | 0.17 | 837
02_16 | 3.88 | -46.16 | 19.90 | -28.55 | - | - | - | - | 334.23 | 0.24 | 723
02_17 | 7.69 | -44.42 | 19.19 | -29.34 | - | - | - | - | 1013.33 | 0.34 | 692
02_18 | 88.59 | -22.14 | - | - | - | - | - | - | 92.16 | 0.06 | 732
02_19 | 28.46 | -37.27 | - | - | - | - | - | - | 746.94 | 0.23 | 754
02_20 | 4.99 | -46.94 | 18.22 | -25.82 | - | - | - | - | 363.63 | 0.22 | 784
02_21 | 73.98 | -24.99 | - | - | - | - | - | - | 117.98 | 0.07 | 751
02_22 | 2.96 | -48.92 | 25.56 | -24.22 | - | - | - | - | 298.52 | 0.24 | 731

Table 5: Material parameters, mean squared errors and time to solution for the
material model discovery from experimental data - brain no. 03. Note that
model parameters that are zero throughout all experiments are not shown in the
table. All material parameters are in Pa, except for the dimensionless
parameters $\alpha_{i}$.

| Material parameters$\frac{\int}{\int}$ | MSE | $\overline{\mathrm{MSE}}$ | Time
---|---|---|---|---
HBE | $\mu_{1}$ | $\alpha_{1}$ | $\mu_{2}$ | $\alpha_{2}$ | $I_{2}-3$ | $\left(I_{2}-3\right)^{2}$ | $\left(I_{2}-3\right)\left(I_{1}-3\right)$ | $\left(I_{1}-3\right)^{2}$ | $[\text{Pa}^{2}]$ | $[-]$ | $[s]$
03_01 | 4.13 | -46.12 | 43.08 | -26.43 | - | - | - | - | 548.09 | 0.20 | 703
03_02 | 38.05 | -19.71 | - | - | - | - | - | - | 9.05 | 0.04 | 694
03_03 | 33.32 | -30.51 | - | - | - | - | - | - | 105.64 | 0.11 | 793
03_04 | 4.24 | -43.24 | 63.20 | -23.46 | - | - | - | - | 241.05 | 0.10 | 753
03_05 | 85.82 | -22.43 | - | - | - | - | - | - | 77.84 | 0.05 | 710
03_06 | 71.81 | -28.86 | - | - | - | - | - | - | 123.95 | 0.04 | 717
03_07 | 0.53 | -50.57 | 12.87 | -22.03 | - | - | - | - | 19.29 | 0.16 | 709
03_08 | 0.97 | -48.88 | 12.33 | -17.28 | - | - | - | - | 30.73 | 0.22 | 711
03_09 | 95.08 | -20.61 | - | - | - | - | - | - | 33.18 | 0.02 | 726
03_10 | 22.46 | -30.21 | - | - | - | - | - | - | 58.02 | 0.13 | 776
03_11 | 52.98 | -16.10 | - | - | - | - | - | - | 4.30 | 0.02 | 708
03_12 | 55.48 | -29.03 | - | - | - | - | - | - | 183.62 | 0.09 | 744
03_13 | 8.98 | -11.23 | - | - | - | - | 41.30 | - | 3.30 | 0.04 | 685
03_14 | 5.19 | -41.87 | - | - | - | - | - | - | 86.43 | 0.30 | 767
03_15 | 52.15 | -29.04 | - | - | - | - | - | - | 344.86 | 0.17 | 702
03_16 | 4.14 | -42.41 | - | - | - | - | - | - | 77.03 | 0.33 | 697
03_17 | 22.37 | -34.83 | - | - | - | - | - | - | 212.28 | 0.18 | 700
03_18 | 52.44 | -18.07 | - | - | - | - | - | - | 7.86 | 0.02 | 725
03_19 | 35.71 | -27.60 | - | - | - | - | - | - | 76.98 | 0.11 | 708
03_20 | 69.96 | -24.83 | - | - | - | - | - | - | 181.74 | 0.11 | 719
03_21 | 6.67 | -33.26 | - | - | - | - | - | - | 8.02 | 0.12 | 748
03_22 | 23.35 | -19.87 | - | - | - | - | - | 102.69 | 7.63 | 0.05 | 781
03_23 | 4.50 | -46.22 | 11.66 | -26.94 | - | - | - | - | 295.75 | 0.26 | 737
03_25 | 63.51 | -29.29 | - | - | - | - | - | - | 366.86 | 0.13 | 795
03_26 | 62.21 | -20.92 | - | - | - | - | - | - | 45.17 | 0.06 | 770
03_27 | 68.54 | -23.96 | - | - | - | - | - | - | 70.19 | 0.05 | 764
03_28 | 191.05 | -18.54 | - | - | - | - | - | - | 50.77 | 0.01 | 740
03_29 | 1.84 | -51.56 | 28.07 | -19.34 | - | - | - | - | 388.02 | 0.35 | 693
03_30 | 58.05 | -21.27 | - | - | - | - | - | - | 28.16 | 0.04 | 688
03_31 | 15.53 | -40.39 | - | - | - | - | - | - | 478.73 | 0.25 | 728
03_32 | 32.07 | -25.14 | - | - | - | - | - | - | 23.06 | 0.06 | 730
03_33 | 54.99 | -23.71 | - | - | - | - | - | - | 45.31 | 0.05 | 745
03_34 | 1.21 | -48.04 | 42.59 | -24.79 | - | - | - | - | 115.99 | 0.12 | 742
03_35 | 29.28 | -25.46 | - | - | - | - | - | - | 28.11 | 0.09 | 775
03_36 | 34.02 | -32.20 | - | - | - | - | - | - | 345.65 | 0.21 | 760
03_37 | 105.48 | -23.43 | - | - | - | - | - | - | 131.43 | 0.05 | 766
03_38 | 81.85 | -22.01 | - | - | - | - | - | - | 49.05 | 0.03 | 725
03_39 | 83.47 | -26.77 | - | - | - | - | - | - | 237.35 | 0.08 | 693
03_40 | 42.47 | -23.87 | - | - | - | - | - | - | 22.89 | 0.05 | 692
03_41 | 2.90 | -41.58 | 46.88 | -29.99 | - | - | - | - | 365.95 | 0.13 | 677
|  |  |  |  |  |  |  | Average | 399.57 | 0.17 | 729

## Acknowledgments

MF and LDL acknowledge funding by the Swiss National Science Foundation (SNF)
through grant N. $200021\\_204316$ “Unsupervised data-driven discovery of
material laws”. SB and PS acknowledge funding by the German Research
Foundation (DFG) through grants N. $460333672/\text{CRC}1540/1-2023$
“Exploring Brain Mechanics (EBM)” and BU 3728/1-1. PS acknowledges support
from the European Research Council (ERC) through grant N. $101052785$,
“Configurational Mechanics of Soft Materials (SoftFrac)”.

## Code and data availability

Codes and data related to the EUCLID project are available at https://euclid-
code.github.io/. The codes corresponding to this work are publically available
in the ETH library (Yu et al., 2023) and the mechanical testing data of the
human brain tissue are available upon request.

## CRediT author statement

MF: Conceptualization, Methodology, Software, Investigation, Writing -
Original Draft, Supervision. HY: Methodology, Software, Formal analysis,
Investigation, Writing - Original Draft. NR: Resources, Data Curation,
Investigation, Writing - Original Draft. JH: Resources, Data Curation,
Investigation, Writing - Original Draft. SB: Conceptualization, Resources,
Writing - Review & Editing, Funding acquisition. PS: Conceptualization,
Resources, Writing - Review & Editing, Funding acquisition. SK:
Conceptualization, Methodology, Writing - Review & Editing. LDL:
Conceptualization, Methodology, Writing - Review & Editing, Supervision,
Funding acquisition.

## References

  * Abdusalamov et al. (2023) Abdusalamov, R., Hillgärtner, M., Itskov, M., 2023.  Automatic generation of interpretable hyperelastic material models by symbolic regression.  International Journal for Numerical Methods in Engineering , nme.7203URL: https://onlinelibrary.wiley.com/doi/10.1002/nme.7203, doi:10.1002/nme.7203.
  * As’ad et al. (2022) As’ad, F., Avery, P., Farhat, C., 2022.  A Mechanics-Informed Artificial Neural Network Approach in Data-Driven Constitutive Modeling, in: AIAA SCITECH 2022 Forum, American Institute of Aeronautics and Astronautics, San Diego, CA & Virtual.  URL: https://arc.aiaa.org/doi/10.2514/6.2022-0100, doi:10.2514/6.2022-0100.
  * As’ad and Farhat (2022) As’ad, F., Farhat, C., 2022\.  A Mechanics-Informed Neural Network Framework for Data-Driven Nonlinear Viscoelasticity URL: https://rgdoi.net/10.13140/RG.2.2.21694.36168, doi:10.13140/RG.2.2.21694.36168. publisher: Unpublished.
  * Bomarito et al. (2021) Bomarito, G., Townsend, T., Stewart, K., Esham, K., Emery, J., Hochhalter, J., 2021\.  Development of interpretable, data-driven plasticity models with symbolic regression.  Computers & Structures 252, 106557\.  URL: https://linkinghub.elsevier.com/retrieve/pii/S0045794921000791, doi:10.1016/j.compstruc.2021.106557.
  * Boyce and Arruda (2000) Boyce, M.C., Arruda, E.M., 2000\.  Constitutive Models of Rubber Elasticity: A Review.  Rubber Chemistry and Technology 73, 504–523.  URL: https://meridian.allenpress.com/rct/article/73/3/504/92749/Constitutive-Models-of-Rubber-Elasticity-A-Review, doi:10.5254/1.3547602.
  * Brunton et al. (2016) Brunton, S.L., Proctor, J.L., Kutz, J.N., 2016.  Discovering governing equations from data by sparse identification of nonlinear dynamical systems.  Proceedings of the National Academy of Sciences 113, 3932–3937.  URL: http://www.pnas.org/lookup/doi/10.1073/pnas.1517384113, doi:10.1073/pnas.1517384113.
  * Budday et al. (2020) Budday, S., Ovaert, T.C., Holzapfel, G.A., Steinmann, P., Kuhl, E., 2020.  Fifty Shades of Brain: A Review on the Mechanical Testing and Modeling of Brain Tissue.  Archives of Computational Methods in Engineering 27, 1187–1230.  URL: http://link.springer.com/10.1007/s11831-019-09352-w, doi:10.1007/s11831-019-09352-w.
  * Budday et al. (2017) Budday, S., Sommer, G., Birkl, C., Langkammer, C., Haybaeck, J., Kohnert, J., Bauer, M., Paulsen, F., Steinmann, P., Kuhl, E., Holzapfel, G., 2017.  Mechanical characterization of human brain tissue.  Acta Biomaterialia 48, 319–340.  URL: https://linkinghub.elsevier.com/retrieve/pii/S1742706116305633, doi:10.1016/j.actbio.2016.10.036.
  * Chen and Guilleminot (2022) Chen, P., Guilleminot, J., 2022\.  Polyconvex neural networks for hyperelastic constitutive models: A rectification approach.  Mechanics Research Communications 125, 103993.  URL: https://linkinghub.elsevier.com/retrieve/pii/S0093641322001252, doi:10.1016/j.mechrescom.2022.103993.
  * Dal et al. (2021) Dal, H., Açıkgöz, K., Badienia, Y., 2021.  On the Performance of Isotropic Hyperelastic Constitutive Models for Rubber-Like Materials: A State of the Art Review.  Applied Mechanics Reviews 73, 020802\.  URL: https://asmedigitalcollection.asme.org/appliedmechanicsreviews/article/73/2/020802/1108153/On-the-Performance-of-Isotropic-Hyperelastic, doi:10.1115/1.4050978.
  * Ester et al. (1996) Ester, M., Kriegel, H.P., Xu, X., 1996.  A Density-Based Algorithm for Discovering Clusters in Large Spatial Databases with Noise.  Proceedings of the 2nd International Conference on Knowledge Discovery and Data Mining, Portland , 226–231.
  * Faber et al. (2022) Faber, J., Hinrichsen, J., Greiner, A., Reiter, N., Budday, S., 2022.  Tissue-Scale Biomechanical Testing of Brain Tissue for the Calibration of Nonlinear Material Models.  Current Protocols 2, e381.  URL: https://doi.org/10.1002/cpz1.381, doi:10.1002/cpz1.381. publisher: John Wiley & Sons, Ltd.
  * Flaschel (2023) Flaschel, M., 2023.  Automated Discovery of Material Models in Continuum Solid Mechanics.  Ph.D. thesis. ETH Zurich.  URL: http://hdl.handle.net/20.500.11850/602750, doi:10.3929/ETHZ-B-000602750.
  * Flaschel et al. (2021) Flaschel, M., Kumar, S., De Lorenzis, L., 2021.  Unsupervised discovery of interpretable hyperelastic constitutive laws.  Computer Methods in Applied Mechanics and Engineering 381, 113852.  doi:10.1016/j.cma.2021.113852.
  * Flaschel et al. (2022) Flaschel, M., Kumar, S., De Lorenzis, L., 2022.  Discovering plasticity models without stress data.  npj Computational Materials 8, 91\.  URL: https://www.nature.com/articles/s41524-022-00752-4, doi:10.1038/s41524-022-00752-4.
  * Flaschel et al. (2023) Flaschel, M., Kumar, S., De Lorenzis, L., 2023.  Automated discovery of generalized standard material models with EUCLID.  Computer Methods in Applied Mechanics and Engineering 405, 115867.  URL: https://linkinghub.elsevier.com/retrieve/pii/S0045782522008234, doi:10.1016/j.cma.2022.115867.
  * Frank and Friedman (1993) Frank, I.E., Friedman, J.H., 1993\.  A Statistical View of Some Chemometrics Regression Tools.  Technometrics 35, 109–135.  URL: http://www.tandfonline.com/doi/abs/10.1080/00401706.1993.10485033, doi:10.1080/00401706.1993.10485033.
  * Frankel et al. (2020) Frankel, A.L., Jones, R.E., Swiler, L.P., 2020.  TENSOR BASIS GAUSSIAN PROCESS MODELS OF HYPERELASTIC MATERIALS.  Journal of Machine Learning for Modeling and Computing 1, 1–17.  URL: http://www.dl.begellhouse.com/journals/558048804a15188a,583c4e56625ba94e,651a2e6b0260f708.html, doi:10.1615/JMachLearnModelComput.2020033325.
  * Friedman et al. (2010) Friedman, J., Hastie, T., Tibshirani, R., 2010.  Regularization Paths for Generalized Linear Models via Coordinate Descent.  Journal of Statistical Software 33\.  URL: http://www.jstatsoft.org/v33/i01/, doi:10.18637/jss.v033.i01.
  * Gent and Thomas (1958) Gent, A.N., Thomas, A.G., 1958\.  Forms for the stored (strain) energy function for vulcanized rubber.  Journal of Polymer Science 28, 625–628.  URL: http://doi.wiley.com/10.1002/pol.1958.1202811814, doi:10.1002/pol.1958.1202811814.
  * Ghaboussi et al. (1991) Ghaboussi, J., Garrett, J.H., Wu, X., 1991.  Knowledge‐Based Modeling of Material Behavior with Neural Networks.  Journal of Engineering Mechanics 117, 132–153.  URL: http://ascelibrary.org/doi/10.1061/%28ASCE%290733-9399%281991%29117%3A1%28132%29, doi:10.1061/(ASCE)0733-9399(1991)117:1(132).
  * Hartmann (2001) Hartmann, S., 2001.  Parameter estimation of hyperelasticity relations of generalized polynomial-type with constraint conditions.  International Journal of Solids and Structures 38, 7999–8018.  URL: https://linkinghub.elsevier.com/retrieve/pii/S002076830100018X, doi:10.1016/S0020-7683(01)00018-X.
  * He et al. (2021) He, H., Zhang, Q., Zhang, Y., Chen, J., Zhang, L., Li, F., 2021.  A comparative study of 85 hyperelastic constitutive models for both unfilled rubber and highly filled rubber nanocomposite material.  Nano Materials Science , S2589965121000490URL: https://linkinghub.elsevier.com/retrieve/pii/S2589965121000490, doi:10.1016/j.nanoms.2021.07.003.
  * Hinrichsen et al. (2022) Hinrichsen, J., Reiter, N., Bräuer, L., Paulsen, F., Kaessmair, S., Budday, S., 2022\.  Inverse identification of region-specific hyperelastic material parameters for human brain tissue.  URL: http://biorxiv.org/lookup/doi/10.1101/2022.12.19.521022, doi:10.1101/2022.12.19.521022.
  * Holzapfel (2000) Holzapfel, G.A., 2000.  Nonlinear solid mechanics: a continuum approach for engineering.  Wiley, Chichester ; New York.
  * Ibañez et al. (2018) Ibañez, R., Abisset-Chavanne, E., Aguado, J.V., Gonzalez, D., Cueto, E., Chinesta, F., 2018\.  A Manifold Learning Approach to Data-Driven Computational Elasticity and Inelasticity.  Archives of Computational Methods in Engineering 25, 47–57.  URL: http://link.springer.com/10.1007/s11831-016-9197-9, doi:10.1007/s11831-016-9197-9.
  * Joshi et al. (2022) Joshi, A., Thakolkaran, P., Zheng, Y., Escande, M., Flaschel, M., De Lorenzis, L., Kumar, S., 2022.  Bayesian-EUCLID: Discovering hyperelastic material laws with uncertainties.  Computer Methods in Applied Mechanics and Engineering 398, 115225.  URL: https://linkinghub.elsevier.com/retrieve/pii/S0045782522003681, doi:10.1016/j.cma.2022.115225.
  * Kabliman et al. (2021) Kabliman, E., Kolody, A.H., Kronsteiner, J., Kommenda, M., Kronberger, G., 2021.  Application of symbolic regression for constitutive modeling of plastic deformation.  Applications in Engineering Science 6, 100052.  URL: https://linkinghub.elsevier.com/retrieve/pii/S2666496821000182, doi:10.1016/j.apples.2021.100052.
  * Kalina et al. (2022) Kalina, K.A., Linden, L., Brummund, J., Metsch, P., Kästner, M., 2022.  Automated constitutive modeling of isotropic hyperelasticity based on artificial neural networks.  Computational Mechanics 69, 213–232.  URL: https://link.springer.com/10.1007/s00466-021-02090-6, doi:10.1007/s00466-021-02090-6.
  * Kirchdoerfer and Ortiz (2016) Kirchdoerfer, T., Ortiz, M., 2016\.  Data-driven computational mechanics.  Computer Methods in Applied Mechanics and Engineering 304, 81–101.  URL: https://linkinghub.elsevier.com/retrieve/pii/S0045782516300238, doi:10.1016/j.cma.2016.02.001.
  * Klein et al. (2022a) Klein, D.K., Fernández, M., Martin, R.J., Neff, P., Weeger, O., 2022a.  Polyconvex anisotropic hyperelasticity with neural networks.  Journal of the Mechanics and Physics of Solids 159, 104703.  URL: https://linkinghub.elsevier.com/retrieve/pii/S0022509621003215, doi:10.1016/j.jmps.2021.104703.
  * Klein et al. (2022b) Klein, D.K., Ortigosa, R., Martínez-Frutos, J., Weeger, O., 2022b.  Finite electro-elasticity with physics-augmented neural networks.  Computer Methods in Applied Mechanics and Engineering 400, 115501.  URL: https://linkinghub.elsevier.com/retrieve/pii/S004578252200514X, doi:10.1016/j.cma.2022.115501.
  * Linka and Kuhl (2022) Linka, K., Kuhl, E., 2022.  A new family of Constitutive Artificial Neural Networks towards automated model discovery.  URL: http://arxiv.org/abs/2210.02202. arXiv:2210.02202 [cond-mat].
  * Linka et al. (2021) Linka, K., Reiter, N., Würges, J., Schicht, M., Bräuer, L., Cyron, C.J., Paulsen, F., Budday, S., 2021\.  Unraveling the Local Relation Between Tissue Composition and Human Brain Mechanics Through Machine Learning.  Frontiers in Bioengineering and Biotechnology 9, 704738.  URL: https://www.frontiersin.org/articles/10.3389/fbioe.2021.704738/full, doi:10.3389/fbioe.2021.704738.
  * Marckmann and Verron (2006) Marckmann, G., Verron, E., 2006\.  Comparison of Hyperelastic Models for Rubber-Like Materials.  Rubber Chemistry and Technology 79, 835–858.  URL: https://meridian.allenpress.com/rct/article/79/5/835/93139/Comparison-of-Hyperelastic-Models-for-RubberLike, doi:10.5254/1.3547969.
  * Marino et al. (2023) Marino, E., Flaschel, M., Kumar, S., De Lorenzis, L., 2023\.  Automated identification of linear viscoelastic constitutive laws with EUCLID.  Mechanics of Materials 181, 104643\.  URL: https://linkinghub.elsevier.com/retrieve/pii/S0167663623000893, doi:10.1016/j.mechmat.2023.104643.
  * Meyer and Ekre (2023) Meyer, K.A., Ekre, F., 2023\.  Thermodynamically consistent neural network plasticity modeling and discovery of evolution laws.  URL: https://engrxiv.org/preprint/view/2961, doi:https://doi.org/10.31224/2961.
  * Ogden (1972) Ogden, R.W., 1972.  Large deformation isotropic elasticity – on the correlation of theory and experiment for incompressible rubberlike solids.  Proceedings of the Royal Society of London. Series A, Mathematical and Physical Sciences, , 565–584.
  * Park and Cho (2021) Park, H., Cho, M., 2021.  Multiscale constitutive model using data–driven yield function.  Composites Part B: Engineering 216, 108831.  URL: https://linkinghub.elsevier.com/retrieve/pii/S1359836821002225, doi:10.1016/j.compositesb.2021.108831.
  * Rivlin (1947) Rivlin, R.S., 1947.  Torsion of a Rubber Cylinder.  Journal of Applied Physics 18, 7\.
  * Steinmann et al. (2012) Steinmann, P., Hossain, M., Possart, G., 2012.  Hyperelastic models for rubber-like materials: consistent tangent operators and suitability for Treloar’s data.  Archive of Applied Mechanics 82, 1183–1217.  URL: http://link.springer.com/10.1007/s00419-012-0610-z, doi:10.1007/s00419-012-0610-z.
  * Sussman and Bathe (2009) Sussman, T., Bathe, K.J., 2009\.  A model of incompressible isotropic hyperelastic material behavior using spline interpolations of tension-compression test data.  Communications in Numerical Methods in Engineering 25, 53–63.  URL: http://doi.wiley.com/10.1002/cnm.1105, doi:10.1002/cnm.1105.
  * Tac et al. (2022) Tac, V., Sahli Costabal, F., Tepole, A.B., 2022.  Data-driven tissue mechanics with polyconvex neural ordinary differential equations.  Computer Methods in Applied Mechanics and Engineering 398, 115248.  URL: https://linkinghub.elsevier.com/retrieve/pii/S0045782522003838, doi:10.1016/j.cma.2022.115248.
  * Thakolkaran et al. (2022) Thakolkaran, P., Joshi, A., Zheng, Y., Flaschel, M., De Lorenzis, L., Kumar, S., 2022\.  NN-EUCLID: Deep-learning hyperelasticity without stress data.  Journal of the Mechanics and Physics of Solids 169, 105076.  URL: https://linkinghub.elsevier.com/retrieve/pii/S0022509622002538, doi:10.1016/j.jmps.2022.105076.
  * Tibshirani (1996) Tibshirani, R., 1996.  Regression Shrinkage and Selection via the Lasso.  Journal of the Royal Statistical Society: Series B (Methodological) 58, 267–288.  URL: http://doi.wiley.com/10.1111/j.2517-6161.1996.tb02080.x, doi:10.1111/j.2517-6161.1996.tb02080.x.
  * Wang et al. (2021) Wang, Z., Estrada, J., Arruda, E., Garikipati, K., 2021\.  Inference of deformation mechanisms and constitutive response of soft material surrogates of biological tissue by full-field characterization and data-driven variational system identification.  Journal of the Mechanics and Physics of Solids 153, 104474.  URL: https://linkinghub.elsevier.com/retrieve/pii/S0022509621001459, doi:10.1016/j.jmps.2021.104474.
  * Yu et al. (2023) Yu, H., Flaschel, M., De Lorenzis, L., 2023.  Supplementary software for ”Automated discovery of interpretable hyperelastic material models for human brain tissue with EUCLID”.  ETH Library doi:https://doi.org/10.5905/ethz-1007-638. doi:https://doi.org/10.5905/ethz-1007-638.