Introduction
Skin is the largest organ in the body and understanding its mechanical properties is a crucial step in many biomedical applications, from the design of prostheses to surgical intervention Tepole et al. (2015); Neumann (1957). The tissue microstructure is characterized by the presence of semiflexible biopolymer fiber networks such as collagen and elastin, which endow skin with nonlinear and anisotropic mechanical behavior Sherman et al. (2017). The mechanical properties of skin are actually common across many collagenous soft connective tissue Kakaletsis et al. (2021); Meador et al. (2020). Traditionally, the mechanical behavior of skin and other soft tissues has been modelled using expertconstructed constitutive equations Holzapfel (2000); Holzapfel and Ogden (2009); Gasser et al. (2005); Humphrey et al. (1990). In this approach, a functional form describing the main features of the mechanical behavior of a family of materials is constructed first. Then, the free parameters in the equations are fitted to a specific material in the family to obtain a calibrated model. Thus, this approach necessitates choosing a specific functional form for the constitutive model. In turn, inherent restrictions of the functional form can result in poor fitting and poor predictions. Unfortunately, even restricting our attention to skin out of all connective soft tissue, there is currently no consensus on the choice of model that is most suitable in a particular application Shergold et al. (2006); Fung (1993); Delalleau et al. (2007); Lapeer et al. (2010); Joodaki and Panzer (2018); Meador et al. (2020); Lee et al. (2018); Ambrosi et al. (2011); Lee et al. (2018).
A new, emergent approach to material modeling is the use of datadriven methods Ambrosi et al. (2011); Rahmani and Sun (2021); Zhang and Garikipati (2020); Peng et al. (2020). Among them, neural network have been used successfully to describe the mechanical behavior of several materials Zhang and Garikipati (2020); Vlassis and Sun (2021); Lu et al. (2020); Lejeune and Zhao (2021). In this approach, the burden of choosing a suitable analytical representation is no longer present, and since the training of the neural network is not restricted to a parametric family of functions, it often results in more accurate predictions than traditional models Leng et al. (2021); Liu et al. (2020); Reimann et al. (2019). In particular, a common application of machine learning has been to create macroscale constitutive models out of micromechanical simulations Reimann et al. (2019); Lejeune and Zhao (2021).
One potential drawback of this approach is scarcity of high fidelity test data to train the datadriven models. It is often difficult, expensive, and sometimes impossible to gather large amounts of data to train a neural network over the entire input space of deformations. One approach to overcome this problem is to use neural networks in a similar vein as the traditional model development, i.e., training a neural network using data from multiple specimens and multiple loading conditions in order to learn an average response. Then, data of a single material can be fitted by adjusting only a few parameters of the neural network Liu et al. (2020).
The route we follow here lies between the purely datadriven approach and the expert modeling approach. Expert models already include knowledge of physics relevant to soft tissue, observations of both the mechanical behavior and underlying microstructure properties, and intuition from the modeller regarding the main features of the material response. For example, to model skin, we have restricted our attention to the framework of hyperelasticity and used expertdesigned strain energy functions to fit murine and porcine skin data Meador et al. (2020); Han et al. (2021). However, the error in the fits can be undesirable, the parameters nonunique, and the predictions can be highly sensitive to the material parameters, especially outside of the loading regime used during calibration. Here we propose to use the analytical strain energy functions as low fidelity approximations, and complement the data with high fidelity experimental measurements. This approach is based on the recent literature that shows the advantage of multifidelity schemes over single fidelity approaches Lejeune and Zhao (2021); Lu et al. (2020); Bonfiglio et al. (2020); Lee et al. (2018).
Specifically, we train fully connected neural networks that output the strain energy and its derivatives with respect to the isochoric strain invariants, including the effect of anisotropy. The training of the neural network is done with a multifidelity dataset that includes low fidelity data from conventional material models, and high fidelity data from biaxial experiments. Furthermore, the loss function is designed to impose convexity constraints on the strain energy, which ensures physically meaningful elasticity tensors. Finally, we implemented a user material (UMAT) subroutine for the widely used finite element package Abaqus Smith (2009). The UMAT reads all the parameters of a fully connected neural network from the input file and evaluates the neural network to construct the stress and elasticity tensors needed in the finite element simulation. The work shown here will extend the reach of machine learning tools to improve the modeling of soft tissue mechanics, in particular through improved constitutive models ready to be used in commercial finite element codes.
Materials and Methods
Constitutive equations for a hyperelastic material with two families of fibers
In this study we use a form of Helmholtz free energy, , that is a function of the right CauchyGreen deformation tensor,
, and two material direction vectors in the reference configuration,
and . This form of the Helmholtz free energy function allows for greater flexibility in recreating the mechanical behavior of materials where more than one family of fibers is present or even when the orientation of fibers is random, which is usually the case in biological tissues. For soft tissues, which we assume to be nearly incompressible, the additive split into isochoric and volumetric parts is used Holzapfel (2000),(1) 
where is the volume change, and the isochoric strain invariants, , , and are defined as
(2)  
The isochoric right Cauchy Green deformation, , can be defined in terms of the isochoric part of the deformation gradient,
(3)  
(4) 
The second PiolaKirchhoff stress tensor, , follows from the DoyleErickson formula by differentiating the strain energy with respect to and following the procedure outlined by Coleman and Noll Doyle and Ericksen (1956); Coleman and Walter (1974), arriving at
(5) 
(6) 
The following definition of the pressure has been introduced . Additionally, the fictitious second PiolaKirchhoff stress tensor, , is the result from differentiating the isochoric part of the strain energy with respect to the isochoric invariants, i.e. , , and . The full expansion of the fictitious stress tensor is
(7) 
The term fictitious originates from the fact that derivatives with respect to the full right Cauchy Green deformation tensor needs the projection with the fourth order tensor
which relates derivatives with respect to the isochoric part of the deformation to derivatives with respect to the total deformation. If the material under consideration is exactly incompressible, i.e. , then , (5) reduces to
(8) 
and the pressure becomes an unknown Lagrange multiplier field.
The UMAT subroutine requires the computation of the Cauchy stress tensor, . Thus, for completeness, we state the standard push forward operation for the stress
(9) 
The finite element subroutine also requires the computation of the elasticity tensor, nel (2020). For ease of derivation, the material version of the elasticity tensor, , is introduced first,
(10) 
(11) 
The expressions for the volumetric and isochoric parts of the elasticity tensor, and , can be further expanded,
(12) 
(13) 
where the modified pressure term, , has been introduced, as well as the special product noted by () defined as , and an additional fourth order projection tensor ,
Finally, the fictitious elasticity tensor, is obtained from differentiating the fictitious stress tensor with respect to the isochoric part of the deformation tensor,
The full expansion of is available in the Appendix. The only remark needed in the main text is that the tensor requires the second derivatives of the strain energy function with respect to the isochoric invariants: , , , etc. This point will become important in the design in the neural network later on.
As stated above, the elasticity tensor needed in the UMAT subroutine is associated with the deformed configuration. The pushforward operation for the elasticity tensor yields
(14) 
where we have introduced the modified dyadic product defined as . The tensor is related to the Truesdell stress rate; however, Abaqus increments employ the Jaumann stress rate. Therefore, the consistent tangent for Abaqus is not Eq. (14) but rather
(15) 
with the modified dyadic product .
Lastly, we remark that during the training of the neural network, incompressibility can be satisfied exactly. In this case can be determined from boundary conditions. However, in the UMAT we use the following function for ,
(16) 
with the bulk modulus. In this study we set kPa.
Neural network structure and training
We use a fully connected neural network to learn the mechanical behavior of skin. The neural network takes four inputs, the isochoric strain invariants in (2), and produces five outputs, the strain energy, , and its derivatives with respect to the invariants, , , and . Note that the notation is used to denote the predicted values of the neural network. The network architecture is summarized in Table 1.
Training data for the neural network is in the form of stretch, stress and strain energy values. Therefore, the first component of the loss function is simply the comparison of the predicted strain energy, , against the observed . Since the derivatives are direct outputs, a regularization term is added to the loss function to enforce that the predicted derivatives coincide with the derivatives of . Using backpropagation, the derivatives of with respect to the inputs are denoted as , with . Thus, the first component of the loss function is
(17) 
where denotes the training point, out of a total of training points.
The second component of the loss results from comparing the stress computed with the neural network outputs against the observed stress . The stress, defined in Eqs. (8) and (9), is computed based on the neural network outputs , to produce . The loss can then be simply stated as
(18) 
where denotes the Frobenius norm.
To ensure optimal convergence of the finite element framework, the strain energy must be a (poly)convex function of the deformation Ehret and Itskov (2007). In the case of a strain energy function defined in terms of the deformation invariants, the strain energy needs to be a convex function of its inputs. For a function to be convex, it’s Hessian matrix, , must be positive semidefinite. The Hessian matrix is
(19) 
The notation indicates the second derivative of the strain energy computed with the neural network by differentiating the outputs with respect to the input using backpropagation. We impose positivedefiniteness of the Hessian matrix using the principal minor test Prussing (1985). For a matrix to be positive definite, it has to be symmetric and all its leading principal minors, , must be positive. This condition is imposed in terms of two additional loss terms,
(20)  
(21) 
The total loss is a weighted sum of the terms discussed so far,
(22) 
If training data from sources with different fidelities are used, the total loss of the multifidelity (mf) dataset is given as a weighted sum of the losses of the high fidelity (hf) and low fidelity (lf) datasets,
(23) 
The training of the neural network was performed using the Adam optimization algorithm Kingma and Ba (2014)
. The initial learning rate was set to 0.0001. The exponential decay rates for first and second moment estimates,
and, were set to 0.9 and 0.99 respectively. The neural network was trained in 100000 epochs without the use of batching. The training was implemented using Keras
Chollet et al. (2018)with a Tensorflow
Abadi et al. (2016) backend on a workstation with the following specifications: Intel Xeon E51630 3.70 GHz CPU, 16 GB DDR4/2400 MHz random access memory, and Nvidia GeForce GTX 1080 GPU. The values of the weights used in this paper are and .Layer  Number of nodes  Activation function 

Input  4  None 
Hidden layer 1  4  Sigmoid 
Hidden layer 2  8  Sigmoid 
Output  5  Linear 
Synthetic data generation
In the majority of biomedical applications it is difficult to obtain sufficient high fidelity data to train a neural network. The number of measurements might be limited, or the data points may be constrained to a narrow region of the input space. It is then beneficial to make use of low fidelity data if available.
In this study, high fidelity data is in the form of biaxial stressstretch measurements. However, only two or three curves within the fourdimensional input space defined by the invariants is explored. Therefore, we generate synthetic data using the GasserOgdenHolzapfel (GOH) Gasser et al. (2005) material model. The GOH model proposes an isochoric strain energy of the form
(24) 
where is vector denoting the mean fiber direction and parameterized by the angle . The functional forms for the GOH strain energy are
(25)  
(26) 
with the generalized fiber strain
(27) 
The volumetric term that is the same as the one we used to penalize volume changes in our formulation Holzapfel (2000); Gasser et al. (2005); Lee et al. (2018),
(28) 
The derivation of the stress tensor for the GOH strain energy is not repeated here, the interested reader is referred to Lee et al. (2018); Gasser et al. (2005).
Synthetic data with the GOH model is generated by fitting the free parameters to the experimental data using the BFGS optimization algorithm in SciPy Virtanen et al. (2020).
Finite element method implementation
We implemented a general neural network material model in a user material subroutine (UMAT) in the nonlinear finite element package Abaqus. The subroutine was written with minimal assumptions to allow for maximal flexibility. The neural network structure, weights and biases, activation functions, etc. are all imported into the subroutine through the input file.
The subroutine performs the following tasks:

Read in the architecture, weights and biases, activation function types, etc., as a set of material properties.

Perform the forward propagation of the neural network to obtain the predicted strain energy and its first derivatives .

Calculate the second derivatives with backpropagation.

Compute the consistent tangent using Eq. (15).
For the forward propagation, let be the output of layer of the neural network with nodes. Then the output of layer is given as
(29) 
where is the elementwise activation function, is the weights matrix, and
is the biases vector of the
layer of the network. We use the sigmoid activation function in the hidden layers,(30) 
For the derivatives, let be the matrix containing derivatives of the nodes of layer with respect to the inputs of the neural network. Then
(31) 
where is the number of inputs to the neural network and diag denotes a diagonal matrix.
Biaxial stressstretch experiments on porcine and murine skin
We use experimental data from biaxial stressstretch experiments performed on murine Meador et al. (2020) and porcine skin for the training and validation of neural networks. The data is collected in up to 5 different experimental protocols which are defined in Table 2.
Loading  

Offx  0  
Offy  0  
Equibiaxial  0  
Stripx  1  0  
Stripy  1  0 
Results
Performance of the neural network against synthetic data
To test the neural network material model, we first train the network using synthetic data only. We generate eleven curves in the space by first holding while is increased gradually to , with . The values for the stretch are . After reaching the corresponding value, is held constant while is gradually increased (Figure 2A). These loading curves are representative of the the type of test that can be performed experimentally. On the other hand, the neural network takes as inputs the isochoric strain invariants. The invariant space is 4dimensional, but we plot a 3dimensional projection in Figure 2B. We use the GOH material model to generate synthetic stress data points and train the neural network. Various components of the loss are plotted in Figure 2C. The predictions of the trained network are plotted against the training data in Figure 2DF. These results indicate that the neural network is able to recreate expert constitutive models within the training region.
We also test if the neural network performs well outside the training region. We generate three validation datasets. The first validation dataset is built by randomly sampling , to construct a diagonal deformation gradient of biaxial deformations not seen during training. Then, to test predictions under shear, which are not directly part of the training data, we construct a data set of deformation gradients of the form
(32) 
The validation dataset is generated from randomly sampling , , . Lastly, we are interested in the potential of the neural network to extrapolate. An additional validation set is constructed by sampling outside the training region: but ; but ; and and . The errors for the validation datasets are shown in Figure 3. It can be seen that the neural network performs well within the training region but worse toward the boundary of the training region. Note that the stress values are very small in magnitude, for small deformations (see Figure 2DF); thus, even though the absolute error is small for , the relative error can achieve close to 100%. This is not a concern since clearly the stresses at these small deformations are two orders of magnitude smaller compared to stresses at larger deformation. For the third validation set, for example, it can be seen that sampling from and leads to errors on the order of 300%. In this region of the input space the stresses are significant. Thus, the neural network has difficulty extrapolating outside of the training region.
Performance against experimental data: multifidelity data and convexity constraints
Next, we start training the neural network using experimental data. We want to test the effect of using the experimental data alone (sparse high fidelity data), or combining these data with the low fidelity approximation of the GOH model fit (multifidelity data). Concurrently, we want to test if the convexity constraint is required to regularize the fits of the neural network. In Figs. 4 and 5, we show the results for murine skin data and porcine skin data respectively together with error in the predictions which is defined as the average Frobenius norm of the error in stress, .
The first row of Figure 4 corresponds to a neural network that is trained using sparse highfidelity data where the convexity constraints are not imposed. Figure 4A shows the neural network ability to fit the offx and offy data, achieving average errors of 2.415 kPa and 2.246 kPa, respectively. A biaxial test, not used in training, is used to test the predictive capability of the network (Figure 4B). The average error in the validation is 6.137 kPa. Because no convexity constraint is used we can see that it is not satisfied (Figure 4C). Keeping only the highfidelity data but imposing convexity drastically changes the performance. The training loss is poorer (Figure 4E and F), but the validation error improves (Figure 4G). Moreover, the function is convex over the input space (Figure 4H).
The third and fourth rows of Figure 4 show the results of neural networks trained with multifidelity data. It is notable that even though the network of the third row is trained without any convexity constraints, the fact that it is trained on the GOH synthetic data (which is an inherently convex model) helps it achieve better convexity (Figure 4H). Nevertheless, it can be seen that for larger deformations, convexity is lost, as it is evident in the biaxial validation test (Figure 4L). The average error in the validation set for the multifidelity case without convexity constraint is 5.415 kPa (Figure 4K), which is comparable to the sparse high fidelity case without convexity requirements (Figure 4C). The results for the neural network trained with multifidelity data and with the convexity constraint are shown in last row of Figure 4. The training errors are 1.896 kPa and 4.328 kPa (Figure 4M and N), and the validation error is 5.678 kPa (Figure 4O). Since the convexity is imposed, the loss in the convexity is close to zero over the entire input space (Figure 4P). In summary, the convexity constraint is needed, but it can slightly increase the error of the neural network over the training set. This can be reflective of the uncertainties in the experimental data collection, or limitations of the hyperelastic framework. Training with multifidelity (given the convexity constraint) yields the best performance: the error in the training data is comparable to the high fidelity case, but the validation error improves.
In Figure 5 we further study the effects of augmenting the training data and how the neural network differs from relying solely on the expert model. For this we focus on porcine skin. We train two neural networks, one of them is trained with the experimental data only (Figure 5AF), whereas the other is trained on the augmented data (Figure 5GL). In Figure 5
C and D it can be seen that trained only on experimental data, the neural network achieves a low average error of 5.410 kPa, compared to the GOH fit which is 53.164 kPa. Thus, the neural network outperforms the GOH model. This is not surprising since the only task of the neural network is to interpolate the experimental data and satisfy convexity. The contour plots in Figure
5B,D and F shows the difference between the neural network as compared with GOH material model throughout the input space. It can be seen that the two models are very different. Without any data on the rest of the input space, it is unclear if the predictions from the neural network are at all accurate. On the other hand, the GOH model has been developed and trained against thousands of tissue biaxial data. It is reasonable to expect that the GOH model, even though it cannot fit any particular dataset as well as the neural network, it can be trusted to guide the neural network away from the training region. We show that training the neural network on the augmented data, the loss is on average 11.705 kPa against the experimental data (Figure 5C and E), which is still lower with respect to using the GOH model alone. However, as looking at the contours in Figure 5H, J and L we see that the neural network now follows the GOH model closely on the rest of the input space.Therefore, the neural network trained with augmented data is at the very least the best version of the GOH model. It performs better than the GOH material model around the high fidelity data points while approximating the GOH model elsewhere.
The last test of the neural network model is also done with porcine experimental data. In this case we have five different biaxial experiments (see Table 2). We are interested in determining which biaxial tests are the most informative for the neural network material model. Thus, we train the neural network with different combination of experimental data and validate against the rest of the data (Figure 6). We do the same training and testing with the GOH model. In Figure 6A and B we train the material models with only two datasets, and test against the other three. In the training set, as expected from our previous result, the neural network outperforms the GOH model. Here we see that the neural network continues to outperform the GOH model even in the validation dataset. Figure 6C, D show the result of training the models with three of the five biaxial curves, and validating against the remaining two. Again, the neural network obviously outperforms the GOH model in the training set. However, surprisingly, when the neural network is trained with offbiaxial data and tested agains stripbiaxial data, it is now outperformed by the GOH model during the validation step (Figure 6C). In Figure 6E we train against four tests, and validate against the equibiaxial test. The validation error is comparable between the neural network and the GOH model. Yet, the superior performance of the neural network on the training data confirms that datadriven models are a preferable alternative to expert constructed constitutive models.
Finite element method implementation
We show a number of basic finite element simulations to test the capabilities of the neural network material model as a UMAT subroutine for Abaqus. The neural network trained on murine skin data (augmented data with convexity constraint) is defined in the input file. We first consider a rectangular block cm^{3}. Boundary conditions, mesh and results for a uniaxial extension simulation with are depicted in Figure 7A. The result is a homogeneous stress distribution of kPa, , consistent with the results in Fig. 4, confirming that the UMAT subroutine is working as intended.
A shearing simulation is shown in Figure 7B. In this analysis the x surface of the prism is clamped and a displacement boundary condition of is applied on the right surface. The contours of the resulting stress components, and are shown in Figure 7. The Appendix shows a simulation with the GOH fit. As discussed in the previous section, the neural network model with the augmented data is, in a way, the best extension of the GOH model: it retains the expert model features but does not suffer from the constraints of an explicit functional form.
The last simulation in Figure 7 is a torsional loading scenario. In this simulation the x surface of the rectangular prism is clamped and a rotation boundary condition of is imposed on the +x surface. The resulting stresses are presented in Figure 7C. This loading scenario is different from the previous two because it involves significant deformations in the outofplane direction. The UMAT subroutine executes without any problems. The three simulations in Figure 7 showcase the robustness and versatility of our neural network UMAT.
Next, we perform a simulation that is much more closely related to skin biomechanics. Tissue expansion is a widely used technique in reconstructive surgery in which a balloonlike device is inserted and inflated subcutaneously to stretch and grow skin Lee et al. (2018). The domain is a cm^{3} patch of skin modeled with 3200 brick elements. A rectangular expander of dimensions cm^{2} underneath the skin mesh is modeled with the fluid cavity feature in Abaqus. The expander is inflated to 20, 40 and 60 cm^{3} resulting in the principal strain distributions shown in Figure 8. Once again, the simulation converged without issues and the results align with our previous experimental observations of higher deformation at the apex and less toward the periphery of the expander Tepole et al. (2014). The simulation in Figure 8 showcases the ability of our neural network model to be used in realistic finite element simulations through our UMAT.
Note that the simulation in Figure 8 evidences the anisotropy of the model. The fiber directions and are aligned with the Cartesian basis and . The tissue is stiffer in the direction, which is why there is a band of higher stress along that direction in Figure 8. To further showcase the anisotropy in the deformation we also plot the corresponding strains (Figure 9).
Discussion
In this study we propose a neural network material model to replace conventional constitutive equations for nonlinear materials, in particular soft collagenous tissues. The neural network takes isochoric strain invariants as inputs and produces the isochoric strain energy and its derivatives with respect to the invariants as outputs. with this design, objectivity is satisfied a priori. Other efforts in datadriven modeling of materials and structures rely on training directly on the stress data, which requires additional steps to ensure objectivity Ghaboussi et al. (1998); Ghaboussi and Sidarta (1998); Kirchdoerfer and Ortiz (2016). Efforts using invariants or principal stretches have also been shown by others Liu et al. (2020); Mihai et al. (2018); Dabiri et al. (2019), and by us as well but for isotropic materials Leng et al. (2021). The key ideas introduced in this paper are the consideration of anisotropy, training with multifidelity data, convexity constraints, and design of the neural network architecture to compute not only the stress but the consistent tangent needed in finite element simulations.
Training data for the neural network can consist of both expensive (or hard to get) high fidelity data such as results of laboratory experiments, or a combination of high fidelity data supplemented by synthetic data from expert models. A control case shown here is to train the neural network on synthetic data alone. The performance of the neural network with the synthetic data shows that the network can interpolate the expert model. This is not surprising since neural networks are universal approximators Leshno et al. (1993). When trained on sparse high fidelity data alone, imposing convexity requirements does not only ensure a physically admissible model, but also improves the performance over the validation set (see Figure 4). This result reflects the smoothness and convexity of the underlying material behavior. However, even with the convexity regularization, there is no guarantee that the neural network predictions will be accurate outside of the training region. In fact, for the control case in which the neural network is trained on synthetic data alone, predictions outside of the training region are not accurate (see Figure 3).
Unfortunately, high fidelity data of soft tissue mechanics is sparse in most applications. In our previous work, typically only three protocols have been performed: offbiaxial x, offbiaxial y and equibiaxial (see Table 2) Meador et al. (2020). Two other tests, stripbiaxial tests in the x and ydirection are explored here as well. Liu et al. Liu et al. (2020) generated datasets by subjecting tissues to seven biaxial tests. Clearly, more coverage of the input space is always better for datadriven approaches. However, this is a challenge, it requires establishment of multiple repeatable protocols and extensive testing of a individual specimens which can introduce unforeseen uncertainties. Expert models of soft tissue mechanics have evolved over the past few decades and reflect our growing understanding of soft tissues. For example, expert models are often based on microstructure observations Limbert (2019); Kassab and Sacks (2016), are physically meaningful Ehret and Itskov (2007); Holzapfel et al. (2000), and have been carefully designed based on observations of many data Humphrey et al. (1990). On the other hand, expert models can also present limitations, such as nonuniqueness of fit, high sensitivity to parameters, and inability to fit the data due to the inherent constraints of the functional form Tonge et al. (2013). By combining the high fidelity data with an expert model as a low fidelity approximation we aim at getting the best of both, keep datacentric models that can capture the experimental data really well, while at the same time maintaining relatively good performance in regions with scarce high fidelity data.
A strong motivation behind the development of constitutive models of soft tissues is to be able to build predictive finite element simulations to guide device design or treatment planning Baillargeon et al. (2014). Previous work on data driven modeling has fallen short in this regard Ghaboussi et al. (1998); Liu et al. (2020); Cilla et al. (2018). We implemented the neural network model in a UMAT subroutine for Abaqus, a popular finite element package in both academia and industry. The UMAT subroutine code was implemented with maximum flexibility in mind. The definition and all parameters of the neural network are fed into the UMAT through the input file. We showcased finite element simulations with the neural network trained on the murine data. The simulations converged stably without any issues in all loading scenarios considered, from simple deformations to realistic applications such as tissue expansion.
Of course, this work is not without limitations. While it is common to model soft tissues within the hyperelastic framework, other physical phenomena will be included in future work, namely viscoelasticity, interstitial flow, and damage. Additionally, a Bayesian framework is needed to account for the inherent uncertainty in material behavior of biological materials. Nevertheless, we anticipate that the general framework introduced here will open up new avenues in datadriven finite element models that balance highfidelity experimental data with expert knowledge of soft tissue mechanics.
Conclusions
The work presented in this study shows that neural network material models can reliably replace or augment conventional constitutive material models in tissue mechanics analyses. If enough high fidelity data is available, datadriven models can eliminate the burden of choosing a specific functional form and the inherent limitations that come with this choice. However, in most applications, high fidelity data is scarce. Our work demonstrates that a multifidelity approach can leverage expert knowledge in the form of synthetic data, while achieving a better fit to the experimental observations. A strong motivation to develop accurate material models of soft tissue is to build predictive finite element models. We designed the neural network with this application in mind, and implemented a neural network UMAT subroutine for Abaqus, a widely used finite element package.
Data Availability
Murine biaxial stressstrain data are available through Manuel K. Rausch’s dataverse (together with other mechanical raw data): https://dataverse.tdl.org/dataverse/STBML Code associated with this manuscript is available at: https://github.com/abuganza/NN_aniso_UMAT
Acknowledgements
This work was supported by the National Institute of Arthritis and Musculoskeletal and Skin Diseases, National Institute of Health, United States under award R01AR074525 to Adrian Buganza Tepole and the National Science Foundation through awards 1916663 and 1916665 to Manuel K. Rausch and Adrian Buganza Tepole, respectively.
References
 Tepole et al. (2015) A. B. Tepole, M. Gart, C. A. Purnell, A. K. Gosain, E. Kuhl, Multiview stereo analysis reveals anisotropy of prestrain, deformation, and growth in living skin, Biomechanics and modeling in mechanobiology 14 (2015) 1007–1019.
 Neumann (1957) C. G. Neumann, The expansion of an area of skin by progressive distention of a subcutaneous balloon; use of the method for securing skin for subtotal reconstruction of the ear, Plastic Reconstructive Surgery 19 (1957) 124–130.
 Sherman et al. (2017) V. Sherman, Y. Tang, S. Zhao, W. Yang, M. Meyers, Structural characterization and viscoelastic constitutive modeling of skin, Acta Biomaterialia 53 (2017) 460–469.
 Kakaletsis et al. (2021) S. Kakaletsis, W. D. Meador, M. Mathur, G. P. Sugerman, T. Jazwiec, M. Malinowski, E. Lejeune, T. A. Timek, M. K. Rausch, Right ventricular myocardial mechanics: Multimodal deformation, microstructure, modeling, and comparison to the left ventricle, Acta Biomaterialia 123 (2021) 154–166.
 Meador et al. (2020) W. D. Meador, M. Mathur, G. P. Sugerman, T. Jazwiec, M. Malinowski, M. R. Bersi, T. A. Timek, M. K. Rausch, A detailed mechanical and microstructural analysis of ovine tricuspid valve leaflets, Acta biomaterialia 102 (2020) 100–113.
 Holzapfel (2000) G. A. Holzapfel, Nonlinear Solid Mechanics; A Continuum Approach for Engineering, John Wiley & Sons, LTD, 2000.
 Holzapfel and Ogden (2009) G. A. Holzapfel, R. W. Ogden, On planar biaxial tests for anisotropic nonlinearly elastic solids. a continuum mechanical framework, Mathematics and Mechanics of Solids 14 (2009) 474–489.
 Gasser et al. (2005) T. C. Gasser, R. W. Ogden, G. A. Holzapfel, Hyperelastic modelling of arterial layers with distributed collagen fibre orientations, Journal of the royal society interface 3 (2005) 15–35.
 Humphrey et al. (1990) J. D. Humphrey, R. K. Strumpf, F. C. P. Yin, Determination of a constitutive model for passive myocardium: I. a new functional form, Journal of Biomechanical Engineering 112 (1990) 333–339.
 Shergold et al. (2006) O. A. Shergold, N. A. Fleck, D. Radford, The uniaxial stress versus strain response of pig skin and silicone rubber at low and high strain rates, International Journal of Impact Engineering 32 (2006) 1384–1402.
 Fung (1993) Y. C. Fung, Biomechanics; Mechanical Properties of Living Tissues, Springer Science + Business Media, Inc, 1993.
 Delalleau et al. (2007) A. Delalleau, G. Josse, J. M. Lagarde, H. Zahouani, J. M. Bergheau, A nonlinear elastic behavior to identify the mechanical parameters of human skin in vivo, Skin Research and Technology 14 (2007) 152–164.
 Lapeer et al. (2010) R. J. Lapeer, P. D. Gasson, V. Karri, Simulating plastic surgery: From human skin tensile tests, through hyperelastic finite element models to realtime haptics, Progress in Biophysics and Molecular Biology 103 (2010) 208–216.
 Joodaki and Panzer (2018) H. Joodaki, M. B. Panzer, Skin mechanical properties and modeling: A review, Journal of Engineering in Medicine 232 (2018).
 Meador et al. (2020) W. D. Meador, G. P. Sugerman, H. M. Story, A. W. Steifert, M. R. Bersi, A. B. Tepole, M. K. Rausch, The regionaldependent biaxial behavior of young and aged mouse skin: A detailed histomechanical characterization, residual strain analysis, and constitutive model, Acta Biomaterialia 101 (2020) 403–413.
 Lee et al. (2018) T. Lee, S. Y. Turin, A. K. Gosain, I. Bilionis, A. B. Tepole, Propagation of material behavior uncertainty in a nonlinear finite element model of reconstructive surgery, Biomechanics and Modeling in Mechanobiology 17 (2018) 1857–1873.
 Ambrosi et al. (2011) D. Ambrosi, G. Ateshian, E. Arruda, S. Cowin, J. Dumais, A. Goriely, G. Holzapfel, J. Humphrey, R. Kemkemer, E. Kuhl, J. Olberding, L. Taber, K. Garikipati, Perspectives on biological growth and remodeling, Journal of the Mechanics and Physics of Solids 59 (2011) 863–883.
 Lee et al. (2018) T. Lee, S. Y. Turin, A. K. Gosain, I. Bilionis, A. B. Tepole, Propagation of material behavior uncertainty in a nonlinear finite element model of reconstructive surgery, Biomechanics and Modeling in Mechanobiology 17 (2018) 1857–1873.
 Rahmani and Sun (2021) B. Rahmani, W. Sun, A kdtreeaccelerated hybrid datadriven/modelbased approach for poroelasticity problems with multifidelity multiphysics data, Computer Methods in Applied Mechanics and Engineering 382 (2021).
 Zhang and Garikipati (2020) X. Zhang, K. Garikipati, Machine learning materials physics: Multiresolution neural networks learn the free energy and nonlinear elastic response of evolving microstructures, Computer Methods in Applied Mechanics and Engineering 372 (2020) 113362.
 Peng et al. (2020) G. C. Peng, M. Alber, A. B. Tepole, W. R. Cannon, S. De, S. DuraBernal, K. Garikipati, G. Karniadakis, W. W. Lytton, P. Perdikaris, et al., Multiscale modeling meets machine learning: What can we learn?, Archives of Computational Methods in Engineering (2020) 1–21.
 Vlassis and Sun (2021) N. N. Vlassis, W. Sun, Sobolev training of thermodynamicinformed neural networks for interpretable elastoplasticity models with level set hardening, Computer Methods in Applied Mechanics and Engineering 377 (2021).

Lu et al. (2020)
L. Lu, M. Dao, P. Kumar,
U. Ramamurty, G. E. Karniadakis,
S. Suresh,
Extraction of mechanical properties of materials through deep learning from instrumented indentation,
Proceedings of the National Academy of Sciences of the United States of America 117 (2020) 7052–7062. 
Lejeune and Zhao (2021)
E. Lejeune, B. Zhao,
Exploring the potential of transfer learning for metamodels of heterogeneous material deformation,
Journal of the Mechanical Behavior of Biomedical Materials 117 (2021) 104276.  Leng et al. (2021) Y. Leng, S. Calve, A. B. Tepole, Predicting the mechanical properties of fibrin using neural networks trained on discrete fiber network data, arXiv preprint arXiv:2101.11712 (2021).
 Liu et al. (2020) M. Liu, L. Liang, W. Sun, A generic physicsinformed neural networkbased constitutive model for soft biological tissues, Computer Methods in Applied Mechanics and Engineering 372 (2020).
 Reimann et al. (2019) D. Reimann, K. Chandra, N. Vajragupta, T. Glasmachers, P. Junker, A. Hartmaier, et al., Modeling macroscopic material behavior with machine learning algorithms trained by micromechanical simulations, Frontiers in Materials 6 (2019) 181.
 Han et al. (2021) T. Han, T. Lee, J. Ledwon, E. Vaca, S. Turin, A. Kearney, A. Gosain, A. B. Tepole, Bayesian calibration of a computational model of tissue expansion based on a porcine animal model, Manuscript submitted for review. (2021).
 Bonfiglio et al. (2020) L. Bonfiglio, P. Perdikaris, S. Brizzolara, Multifidelity bayesian optimization of swath hull forms, Journal of Ship Research 64 (2020) 154–170.
 Smith (2009) M. Smith, ABAQUS/Standard User’s Manual, Version 6.9, Dassault Systèmes Simulia Corp, United States, 2009.
 Doyle and Ericksen (1956) T. Doyle, J. Ericksen, Nonlinear elasticity, volume 4 of Advances in Applied Mechanics, Elsevier, 1956, pp. 53–115.
 Coleman and Walter (1974) B. D. Coleman, N. Walter, The Foundations of Mechanics and Thermodynamics, Springer Link, 1974.
 nel (2020) How to implement userdefined fiberreinforced hyperelastic materials in finite element software, Journal of the Mechanical Behavior of Biomedical Materials 110 (2020) 103737.
 Ehret and Itskov (2007) A. E. Ehret, M. Itskov, A polyconvex hyperelastic model for fiberreinforced materials in application to soft tissues, Journal of Materials Science 42 (2007) 8853–8863.
 Prussing (1985) J. E. Prussing, The principal minor test for semidefiniteness matrices, Journal of Guidance, Control and Dynamics 9 (1985) 121–122.
 Kingma and Ba (2014) D. P. Kingma, J. Ba, Adam: A method for stochastic optimization, arXiv preprint arXiv:1412.6980 (2014).
 Chollet et al. (2018) F. Chollet, et al., Keras: The python deep learning library, ascl (2018) ascl–1806.
 Abadi et al. (2016) M. Abadi, A. Agarwal, P. Barham, E. Brevdo, Z. Chen, C. Citro, G. S. Corrado, A. Davis, J. Dean, M. Devin, et al., Tensorflow: Largescale machine learning on heterogeneous distributed systems, arXiv preprint arXiv:1603.04467 (2016).
 Virtanen et al. (2020) P. Virtanen, R. Gommers, T. E. Oliphant, M. Haberland, T. Reddy, D. Cournapeau, E. Burovski, P. Peterson, W. Weckesser, J. Bright, S. J. van der Walt, M. Brett, J. Wilson, K. J. Millman, N. Mayorov, A. R. J. Nelson, E. Jones, R. Kern, E. Larson, C. J. Carey, İ. Polat, Y. Feng, E. W. Moore, J. VanderPlas, D. Laxalde, J. Perktold, R. Cimrman, I. Henriksen, E. A. Quintero, C. R. Harris, A. M. Archibald, A. H. Ribeiro, F. Pedregosa, P. van Mulbregt, SciPy 1.0 Contributors, SciPy 1.0: Fundamental Algorithms for Scientific Computing in Python, Nature Methods 17 (2020) 261–272.
 Tepole et al. (2014) A. B. Tepole, M. Gart, A. K. Gosain, E. Kuhl, Characterization of living skin using multiview stereo and isogeometric analysis, Acta biomaterialia 10 (2014) 4822–4831.
 Ghaboussi et al. (1998) J. Ghaboussi, D. A. Pecknold, M. Zhang, R. M. HajAli, Autoprogressive training of neural network constitutive models, International Journal for Numerical Methods in Engineering 42 (1998) 105–126.
 Ghaboussi and Sidarta (1998) J. Ghaboussi, D. Sidarta, New nested adaptive neural networks (nann) for constitutive modeling, Computers and Geotechnics 22 (1998) 29–52.
 Kirchdoerfer and Ortiz (2016) T. Kirchdoerfer, M. Ortiz, Datadriven computational mechanics, Computer Methods in Applied Mechanics and Engineering 304 (2016) 81–101.
 Mihai et al. (2018) L. A. Mihai, T. E. Woolley, A. Gorieli, Stochastic isotropic hyperelastic materials: constitutive calibration and model selection, Proceedings of the Royal Society A 474 (2018).
 Dabiri et al. (2019) Y. Dabiri, A. Van der Velden, K. L. Sack, J. S. Choy, G. S. Kassab, J. M. Guccione, Prediction of left ventricular mechanics using machine learning, Frontiers in Physics 7 (2019) 117.
 Leshno et al. (1993) M. Leshno, V. Y. Lin, A. Pinkus, S. Schocken, Multilayer feedforward networks with a nonpolynomial activation function can approximate any function, Neural networks 6 (1993) 861–867.
 Limbert (2019) G. Limbert, Skin Biophysics: From Experimental Characterisation to Advanced Modelling, volume 22, Springer, 2019.
 Kassab and Sacks (2016) G. S. Kassab, M. S. Sacks, Structurebased mechanics of tissues and organs, Springer, 2016.
 Holzapfel et al. (2000) G. A. Holzapfel, T. C. Gasser, R. W. Ogden, A new constitutive framework for arterial wall mechanics and a comparative study of material models, Journal of elasticity and the physical science of solids 61 (2000) 1–48.
 Tonge et al. (2013) T. K. Tonge, L. M. Voo, T. D. Nguyen, Fullfield bulge test for planar anisotropic tissues: Part ii – a thin shell method for determining material parameters and comparison of two distributed fiber modeling approaches, Acta Biomaterialia 9 (2013) 5926–5942.
 Baillargeon et al. (2014) B. Baillargeon, N. Rebelo, D. D. Fox, R. L. Taylor, E. Kuhl, The living heart project: a robust and integrative simulator for human heart function, European Journal of MechanicsA/Solids 48 (2014) 38–47.
 Cilla et al. (2018) M. Cilla, I. PérezRey, M. A. Martínez, E. Peña, J. Martínez, On the use of machine learning techniques for the mechanical characterization of soft biological tissues, International Journal for Numerical Methods in Biomedical Engineering 34 (2018) e3121. E3121 cnm.3121.
Comments
There are no comments yet.