A Data-Driven Learning Method for Constitutive Modeling: Application to Vascular Hyperelastic Soft Tissues

We address the problem of machine learning of constitutive laws when large experimental deviations are present. This is particularly important in soft living tissue modeling, for instance, where large patient-dependent data is found. We focus on two aspects that complicate the problem, namely, the presence of an important dispersion in the experimental results and the need for a rigorous compliance to thermodynamic settings. To address these difficulties, we propose to use, respectively, Topological Data Analysis techniques and a regression over the so-called General Equation for the Nonequilibrium Reversible-Irreversible Coupling (GENERIC) formalism (M. Grmela and H. Ch. Oettinger, Dynamics and thermodynamics of complex fluids. I. Development of a general formalism. Phys. Rev. E 56, 6620, 1997). This allows us, on one hand, to unveil the true “shape” of the data and, on the other, to guarantee the fulfillment of basic principles such as the conservation of energy and the production of entropy as a consequence of viscous dissipation. Examples are provided over pseudo-experimental and experimental data that demonstrate the feasibility of the proposed approach.


Introduction
Computational (bio-)mechanics is not absent of the data "fever". Even if the so-called data-driven computational mechanics discipline presents distinctive features over what is commonly known as "big data", the possibility of employing raw experimental data to perform simulations has attracted the attention of many researchers recently [1][2][3][4][5][6][7][8][9]. In these approaches, basic equations-those with a higher epistemic value, such as equilibrium or compatibility-are kept, while those which are frequently phenomenological in nature-constitutive equations-are substituted by experimental data, and hence, the name of this family of methods. In a related approach, the so-called equation-free approach substitutes the constitutive law of a material not by experimental data, but by a pseudo-experimental result coming from microscale, first-principle computations [10].
This approach allows us to avoid complex, costly-and often unsuccessful-parameter fitting to obtain the precise form of a constitutive equation. On the other hand, widely accepted constitutive equations have been built under rigorous standards that give rise to the automatic satisfaction of first principles. This is not evident for data-driven approaches. and experimental tests, described in Sections 2.3 and 2.4, respectively. The numerical fitting results of both experiments are shown in Section 3. The paper ends in Section 4 with a discussion of the obtained results on the use of the proposed GENERIC-TDA methodology on the mechanical modeling of biological tissues.

A GENERIC Approach to the Learning Procedure
Recently, W. E and coworkers established a very useful parallelism between DNNs and dynamical systems [33]. Consider a system governed by some state variables z(z 0 , t) : I → S, z ∈ C 1 (0, T], such that its time evolution is established in the form where f is, in general, a nonlinear function-otherwise, the method is of little interest. S represents the phase space of the system (a set of judiciously chosen variables that both describe the energy of the system and can be measured experimentally). I = (0, T] represents the considered time interval. For the selected time horizon T, the f low map is a nonlinear function of z. The dynamical system approach to supervised learning consists in determining f so that the resulting flow map is able to reproduce the experimental data. This parallelism offers the advantage of the vast knowledge developed so far in the field of dynamical systems, which could help us in developing both theoretical insights about DNNs and also to devise alternative routes for developing efficient learning strategies. DNNs can be thought of as discretized dynamical systems, such that every time step corresponds to a layer of the DNN. In this work, we consider viscous-hyperelastic materials. Therefore, the learned function f must satisfy certain well-known principles dictated by thermodynamics. In the hyperelastic case, the right choice for f could arise from Hamiltonian mechanics, i.e., where E represents the Hamiltonian or, in other words, the energy of the system. L(z) represents the so-called Poisson matrix, a skew-symmetric matrix that depends on z, in general, but that results to be constant in many cases. Should the system of interest not be conservative (or Hamiltonian), a new potential needs to be introduced in the formulation-entropy-giving rise to the GENERIC formalism [30] For Equation (2) to represent valid nonequilibrium thermodynamic processes, it must be supplemented with the so-called degeneracy conditions, i.e., If, as stated before, L is skew-symmetric; and choosing M to be symmetric, positive semidefinite, we obtainĖ (z) = ∇E(z) ·ż = ∇E(z) · L(z)∇E(z) + ∇E(z) · M(z)∇S(z) = 0 i.e., one ensures the conservation of energy in closed systems. In turn,Ṡ (z) = ∇s(z) ·ż = ∇S(z) · L(z)∇E(z) + ∇S(z) · M(z)∇S(z) ≥ 0 which is equivalent to the fulfillment of the second principle of thermodynamics. Thus, we notice how, by leveraging the dynamical systems equivalence, we efficiently enforce the conservation of energy and the production of entropy. For an in-depth discussion of the implications of the choice of phase space variables z, we refer the interested reader to our previous works in the field, [31,32]. In general, it is well-known that neither a particular choice of variables is needed, nor a particular scale for the description of the system at hand. The only need is that the variables in the phase space would be able to account for the energy of the system. That is why, in general, a Langevin equation is not suitable for this purpose, see [35] and references therein.
Following the dynamical systems equivalence, the next step consists in the determination, by regression from data, of the form of the elements of the GENERIC description of the dynamics, i.e., L, M, ∇E(z), and ∇S(z). To do so, assume a standard finite difference discretization of the time derivative, where we denote z n+1 = z t+∆t and where L and M are the discrete version of the Poisson and friction operators, respectively. DE and DS represent the discrete gradients. In general, matrix L is constant over the process, while matrix M frequently varies. While there exist machine learning techniques that are able to provide with the precise expressions of the terms involved in a particular PDE from data [2], we pursue a purely numerical route, in which we assume that the elements of the GENERIC formula have a particular manifold structure which is to be unveiled by our method. This is the concept of constitutive manifold that we first stated in our previous works [4,7,36].
The dynamical systems approach to the problem will thus consist of solving the following (possibly constrained) minimization problem within a time interval J ⊆ I: with z meas ⊆ Z, a subset of the total available experimental results. L, M, DE, DS will in general be dependent on z, and therefore the choice of size of z meas (in other words, the number of regressions performed along the time interval, or the number of layers in the DNNs equivalent) will affect the accuracy of the method. An analysis of the influence of this choice was made in [31]. We finally approximate z by employing piecewise polynomials, so that gradient operators can be cast in matrix form, à la finite elements, as where A and B represent the discrete, matrix form of the gradient operators leading to DE and DS, respectively. The GENERIC description of a hyperelastic material is well known [37]. Indeed, for the Hamiltonian, conservative part of the constitutive equation, we have where x = φ(X) ∈ Ω t ⊂ R 3 represents the deformed configuration of the solid, Ω t represents the deformed configuration of the solid at time t, and p ∈ R 3 represents the material momentum density. In this case, The total energy of a hyperelastic body is known to be the sum E = W + K of elastic and kinetic energies. Here, we assume a strain energy density potential w of the form where Ω 0 represents the undeformed configuration of the solid and C represents the right Cauchy-Green deformation tensor. In a general isotropic case, the strain energy density would take the form w = w(X, C, S). If the material under consideration is isotropic hyperelastic, we simply write w = w(C). In turn, the kinetic energy will be Therefore, by numerically identifying the form of the energy potential, one readily observes that the conservative part of the usual constitutive hyperelastic law is found. If, in addition, the material shows a viscous dissipative behavior, the precise form of the entropy potential is to be found.
Once the learning method has been developed, we describe next how the experimental campaign was accomplished.

Treatment of Dispersion and Noise in Data
One of the most important aspects to consider when dealing with soft living tissues is the importance of the dispersion of experimental results. This will be highlighted in Section 3 below. The main objective of the proposed technique is the determination of a constitutive manifold for the term of the GENERIC description of the physical phenomena, this type of dispersion needs to be treated efficiently.
To this end, we suggest to employ Topological Data Analysis (TDA) [38,39]. Behind the proposed method is the assumption of the existence of a constitutive manifold, a concept that we first introduced in [4]. Our set of experimental measurements, in the most general case, is assumed to be composed by D-tuples (D = 9 for this particular case) of the type where U represents the (right) stretch tensor and P the first Piola-Kirchhoff stress tensor. These experimental values are assumed to form a constitutive manifold in a high-dimensional space, see Figure 1. These values are, however, noisy. TDA is nevertheless able to extract the underlying geometry of this manifold, which will later be embedded onto a low-, d-dimensional manifold for interpolation purposes.

Smooth mapping Z(ξ)
Dimensionality reduction Figure 1. Hypothesis about the existence of a constitutive manifold in which the experimental results live. Despite the noise in the data, Topological Data Analysis (TDA) techniques will help us in unveiling the true geometry of the manifold in the high-dimensional setting, which will later be embedded onto a low dimensional space for the ease of computations.
TDA seeks to find the underlying topological structure of data. To this end, it employs a distance parameter R between experimental points. As we make R grow, points at distances lower than R will be connected by edges, triangles, and tetrahedra-in general, by k-simplexes, respectively, in one, two, three, ..., k dimensions. As simplexes appear, they form simplicial complexes. The study of the formation of this simplicial complex structure in the data is precisely the objective of TDA. The optimal R parameter that best describes the topology of the data is found by resorting to the so-called persistence diagrams.
In essence, simplicial homology reflects the number of holes in a given dimension for the simplicial complex. For instance, a chain of edges may close a hole, while the interior space within a tetrahedron is a void. If we record the precise value of R for which a hole or void (in any dimension) appears, and the R value for which it disappears-by the appearance of a filled triangle or tetrahedron, for instance-the resulting plots will indicate which topological features persist more. Those indicate the true topology of the data. For a graphical interpretation of this explanation, see Figure 2.
Under the TDA framework, noisy data produce topological structures with small persistence. The true topology (manifold structure) is described by these structures that persist the most. Once unveiled, the topological structure of data or, equivalently, the right shape of the data manifold, will allow us to interpolate data in the right topological space-in the tangent plane to the manifold at each datum-and not in Euclidean space.  code, we notice that the data set has the topology of a circle, with one single interior hole. In general, those holes or voids that persist the most reflect the persistent topology of the data.

Pseudo-Experimental Data-Learning a Visco-Hyperelastic Response
In this first example, we consider pseudo-experimental (numerical) data coming from a finite element simulation. We consider the same visco-hyperelastic material previously considered in [32], but this time altered with noise. This noise has a standard deviation of 10% of the mean value.
The viscous part of the behavior is modeled after a Prony series expansion of the type With the material as just described, a plane-stress, biaxial experiment is reproduced. This experiment consists of two different loading steps, each one followed by a relaxation step. A total of 50 different data sets have been obtained for this experiment.
A mean GENERIC model. A regression procedure is then accomplished for each one of the 50 different experiments, so as to determine their precise GENERIC expression. With the obtained values, we first compute the mean GENERIC model by simply taking mean values for each one of the GENERIC model components. This "mean" GENERIC model is compared to the noise-free numerical experiment, taken as ground truth.
Extracting the topology of data: GENERIC-TDA model. Instead of just computing the mean values of each term of the GENERIC model, it seems judicious to employ TDA to unveil the topology of data and to determine the final GENERIC model by interpolating from the right neighboring experimental results. To this end, we considered a data set composed by 20 different loading states (all consisting of a load-relaxation-load-relaxation sequence) and the addition of noise (10% sdv) to each one of these processes, so as to obtain 50 different tests for each one of the 20 loading processes. This makes a total of one thousand different tests.
One of these tests (noise-free) is kept as the reference solution. We employ TDA to find the "neighboring" test to this reference solution and, by employing Kriging, to obtain the final numerical values of the GENERIC model. We would like to emphasize that interpolation is made not in the Euclidean space, but in the right tangent space to the data manifold, as unveiled by TDA techniques. Note that we speak of a tangent plane to a noisy manifold, since TDA is able to unveil the topology that persists the most, thus giving an apparent noise-free topology.
Weights for the interpolation are obtained by different Kriging interpolation techniques: Simple, Ordinary, and Local.

Learning the Constitutive Model of Porcine Carotid Tissue
One of the most intricate experimental procedures in the framework of solid mechanics is perhaps that of constitutive modeling of soft living tissues. Here, we will employ data previously obtained and presented in [14] for the constitutive modeling of porcine carotid tissue.
What is remarkable in soft living tissue modeling is, on one hand, the heterogeneity and anisotropy of the tissue; and on the other, the large differences between experimental values found in different specimens. For instance, the behavior of porcine carotid tissue-whose interest is to serve as a proxy of that of humans-differs strongly if the sample is extracted from proximal (i.e., close to the heart) positions of the vessel or if it is extracted from distal positions. Additionally, there is a strong anisotropy regarding circumferential versus longitudinal behavior.
In [14], a traditional fitting procedure was accomplished so as to determine the best fitting model for these data. Taking the mean of the experimental results arising from [14] as the only plausible reference solution-this is standard experimental procedure in the literature-what we did first was to determine by TDA the set of experimental results neighboring this (mean) reference solution. A GENERIC model was then determined for each one of these neighbor results. Three different approaches were then compared: (i) a GENERIC model whose terms are computed as the mean values of each of the GENERIC models for each neighboring curve (the closest ones in the data manifold); and either (ii) ordinary or (iii) local Kriging interpolation techniques among these neighbors of the terms of a new GENERIC model. The result of our approach is a new GENERIC-TDA model whose integration in (pseudo-)time produces a prediction of the tissue behavior.
Experimental tests. To introduce to the reader the most significant details in the experimental models used for our numerical analysis, here we show a brief description of the sample's harvesting and tensile test protocols performed in [14]. The interested reader is referred to this article for a precise description of the experimental campaign.
We consider nine female pigs of 3.5 ± 0.45 months (mean ± SD). The experiments on these swine were approved by the Ethical Committee for Animal Research of the University of Zaragoza. All procedures were carried out in accordance with the Principles of Laboratory Animal Care (86/609/EEC Norm, incorporated into Spanish legislation through the RD 1021/2005).
For each one of the left and right carotids, proximal and distal regions were considered for mechanical testing, as mentioned before. At each location, circumferential and longitudinal strips, approximately 3-mm-wide and 11-mm-long, and 5-mm wide and 15-mm long, respectively, were cut. A total of 14 carotid specimens with 47 and 49 valid tests were performed for the proximal and distal zones, respectively. A minimum number of two test strips along each direction was accomplished for each specimen.
Simple tension tests of the carotid strips were performed in a high-precision-drive Instron Microtester 5548 system, see Figure 3. The procedure is properly described in [14], and consisted of the Instron Microtester 5548 System with two clamps holding the sample. The samples are subjected to the tensile test under a humidity-controlled environment to prevent sample drying.
The applied force was measured with a 5-N load cell with a minimal resolution of 0.001 N, and the axial strain was measured using a noncontact Instron 2663-281 video-extensometer equipped with a high-performance digital camera with a megapixel sensor (0.5 m ± 0.5% ). Different loading and unloading cycles were applied that correspond to 60, 120, and 240 kPa (50%, 100%, and 200% of the estimated physiological stress state in the artery) at 30%/min of strain rate, which can be considered as quasi-static. Therefore, these experiments serve for the hyperelastic modeling of soft tissues, but not for their viscous characterization. The resulting GENERIC model will therefore be purely Hamiltonian under these conditions.
For each one of the fourteen carotids, proximal and distal measurements are obtained, thus giving a total of 28 datasets. Each dataset includes circumferential and longitudinal stresses and stretches, z = {σ c , λ c , σ , λ }. Due to the quasi-static nature of the experiments, time is actually a pseudo-time. In addition to these 14 stress-stretch curves for each location, a fifteenth curve is obtained by computing the mean value of the first 14. This curve will be taken as a sort of reference for comparison purposes, see

Numerical Fitting of the Pseudo-Experimental Data Set
Recalling Section 2.3, Figure 5 shows the "mean" GENERIC model when compared to the noise-free numerical experiment. As can be noticed from this figure, results show a poor accordance to the noise-free version of the data. Constructing a model by just computing the mean of each GENERIC model for noisy data seems not to be a good idea. If we consider it here, it is just because in the experimental framework, phenomenological models are very often obtained after computing means of the available results [14]. Figure 6 shows a displacement comparison among the noise-free sample and the full GENERIC-TDA model with different Kriging interpolation techniques. Additionally, Table 1 shows the obtained 2-norm errors of the mentioned model results.  It is worth noting the high degree of accuracy obtained by employing local Kriging procedures. In combination with TDA, this procedure is able not just to filter the artificial noise added to the data, but to provide a very accurate GENERIC model able to reproduce the visco-hyperelastic model from which pseudo-experimental data was obtained.

Numerical Fitting of Porcine Carotid Tissue
For the mean experimental curves (circumferential and longitudinal) of the distal samples, results are shown in Figure 7. It is worth noting that, in general, it seems not to be a good idea to just compute the mean values of the GENERIC models for each of the neighboring experimental curves. Particularly in the circumferential direction, the deviation from the reference solution is noteworthy. As in the previous section, results provided by Kriging interpolation outperform this approach. Particularly, local Kriging is found to provide the highest accuracy. The predicted behavior is almost indistinguishable from the mean experimental results. but to provide with a very accurate GENERIC model able to reproduce the visco-hyperelastic model from which pseudo-experimental data was obtained.

Numerical Fitting of Porcine Carotid Tissue
For the mean experimental curves (circumferential and longitudinal) of the distal samples, results are shown in Figure 7. It is worth noting that, in general, it seems not to be a good idea to just compute the mean values of the GENERIC models for each of the neighboring experimental curves. Particularly in the circumferential direction, the deviation from the reference solution is noteworthy. As in the previous section, results provided by Kriging interpolation outperform this approach. Particularly, local Kriging is found to provide the highest accuracy. The predicted behavior is almost indistinguishable from the mean experimental results.  With the weights just computed for the distal samples, we constructed a new GENERIC model for the proximal results. Its predictions are shown in Figure 8. Again, by just computing the mean of the GENERIC terms for the neighboring curves does not seem to produce good results. However, With the weights just computed for the distal samples, we constructed a new GENERIC model for the proximal results. Its predictions are shown in Figure 8. Once again, by just computing the mean of the GENERIC terms for the neighboring curves does not seem to produce good results. However, with the Kriging weights computed for the distal samples, results for the proximal samples are equally good. This demonstrates the robustness of the proposed approach.
with the Kriging weights computed for the distal samples, results for the proximal samples are equally good. This demonstrates the robustness of the proposed approach.  It is worth noting that the computational cost of this procedure is by no means high: each sample took on average 2.52 seconds to obtain the corresponding GENERIC model.

Conclusions
In this paper a new methodology for the data-driven learning of constitutive models is proposed. We made an emphasis on those cases in which large experimental deviations are present in the data. By employing first Topological Data Analysis techniques we unveil the shape of the data manifold, It is worth noting that the computational cost of this procedure is by no means high: each sample took on average 2.52 seconds to obtain the corresponding GENERIC model.

Conclusions
In this paper, a new methodology for the data-driven learning of constitutive models is proposed. We made an emphasis on those cases in which large experimental deviations are present in the data. By first employing Topological Data Analysis techniques, we unveil the shape of the data manifold so as to allow us to perform interpolation on the right tangent plane to the manifold. Once the neighboring data are found, a GENERIC expression is found for the material under consideration. In other words, the precise form of the strain energy density and entropy potentials are found. This allows us to predict new loading states to a high degree of accuracy without the need to perform complex parameter fitting procedures to arrive to phenomenological models.
In addition, and in sharp contrast to other existing alternatives, our method is able to guarantee exact (to numerical precision) satisfaction of thermodynamic principles-conservation of energy and positive production of entropy-thanks to the GENERIC formalism.
To the best of our knowledge, no work has been performed in this line applied to biomedical living tissues. Despite the limitation of needing an admissible database to perform the learning process of our method, we strongly believe that the proposed GENERIC-TDA technique can be applied to the numerical fitting of highly nonlinear materials with sound accuracy, as shown in this manuscript. As a proof of concept, our results (developed in both synthetic and real experiments) show the high benefits of using data-driven models for materials simulation in fields where complex physical responses are present. We believe that machine learning methods combined with numerical modeling for biological systems (at any scale) is a very exciting young field with countless challenges and potential usefulness to both biomedical and numerical communities. Funding: This work has been supported by the Spanish Ministry of Economy and Competitiveness through Grant number DPI2017-85139-C2-1-R and by the Regional Government of Aragon and the European Social Fund, research group T24 20R. The support given by ESI Group to F.C. through the ESI Group Chair at ENSAM Paris and tho D.G. and E.C. through the project 2019-0060 "Simulated Reality" is also gratefully acknowledged.