Predicting the Mechanical Properties of Fibrin Using Neural Networks Trained on Discrete Fiber Network Data
PPredicting the Mechanical Properties of Fibrin Using Neural Networks Trainedon Discrete Fiber Network Data
Yue Leng , Sarah Calve , and Adrian B Tepole , Weldon School of Biomedical Engineering, Purdue University, West Lafayette, IN, USA Department of Mechanical Engineering, Colorado University, Boulder, CO, USA School of Mechanical Engineering, Purdue University, West Lafayette, IN, USA
Abstract
Fibrin is a structural protein key for processes such as wound healing and thrombus formation. At the macroscale,fibrin forms a gel and the overall mechanical response is dictated by the mechanics of a microscale fiber network.Hence, accurate description of fibrin gels can be achieved using representative volume elements (RVE) that explicitlymodel the discrete fiber networks of the microscale. These RVE models, however, cannot be e ffi ciently used to modelthe macroscale due to the challenges and computational demands of coupling the micro- and macroscale models.Here, we propose the use of an artificial, fully connected neural network (FCNN) to e ffi ciently capture the behaviorof the RVE models. The FCNN was trained on 1100 discrete fiber networks subjected to 121 biaxial deformations.The stress data from the RVE, together with the total energy on the fibers and the condition of incompressibility ofthe surrounding matrix, were used to determine the derivatives of an unknown strain energy function with respect tothe deformation invariants. During training, the loss function was modified to ensure convexity of the strain energyfunction and symmetry of its Hessian. An invariant representation is compatible with standard finite element codes.Thus, a general FCNN model was coded into a user material subroutine (UMAT) in the commercial software Abaqus.The UMAT implementation takes in the structure and parameters of an arbitrary FCNN as material parameters fromthe input file. The inputs to the FCNN include the first two isochoric invariants of the deformation, and the FCNNoutputs are the derivatives of the strain energy with respect to the isochoric invariants. Second derivatives of the strainenergy function are automatically calculated inside the UMAT through back-propagation. In this work, the FCNNtrained on the discrete fiber network data was used in finite element simulations of fibrin gels using our UMAT.We anticipate that this work will enable further integration of machine learning tools with standard computationalmechanics tools such as finite element analysis. It will also improve computational modeling of biological materialscharacterized by a multiscale structure. Keywords:
Machine Learning, Nonlinear finite elements, Constitutive modeling, Abaqus User Subroutine UMAT,Multiscale modeling
Introduction
It is well known that the biomechanical macroscopic behavior of both native and engineered tissues is largelydetermined by the properties of the underlying microstructural components [1]. Therefore, to allow for a better designof engineered tissues and to increase our fundamental understanding of the mechanical behavior of soft tissue, it isnecessary to investigate this multiscale coupling. Fibrin is an important extracellular matrix (ECM) component in thebody, providing structural integrity to various tissues [2, 3]. It plays a critical role in wound healing and thrombus fate[4, 5]. In adddition, it is also a common sca ff old material used in tissue engineering [6]. The fibrin network structurecan be described by variables such as the fiber diameter and length, the number of branch points, and the volumefraction [7].Due to the inherent multiscale nature of fibrin gels, with the macroscopic dimension described at the centime-ter length scale and the underlying fibrin network defined at the micrometer length scale, our understanding of themacroscale mechanics needs to be informed by the microscale mechanical behavior of the fibrin network. To enable Preprint submitted to ArXiv January 29, 2021 a r X i v : . [ q - b i o . Q M ] J a n he multiscale description, models of discrete fiber networks (DFN) in a representative volume element (RVE) haveoften been proposed [8, 9]. These micromechanical models, based on volume average theory, assumed that the stressand strains at the macroscopic scale are volume averages of the corresponding microscale fields in the RVE, and thusprovide the linkage between the fiber-level mechanics and the tissue-level behaviors [7, 10, 11]. DFN models are ableto predict the mechanical behavior of three-dimensional fibrin networks. The volume-averaging homogenization thatsatisfies the Hill-Mandel condition can then be used as a constitutive model at the macroscale, where the standard bal-ance of linear momentum within continuum mechanics can be solved to determine the mechanical response [12, 13].Moreover, at the macroscale, hyperelasticity is a common assumption that is used to describe soft tissue mechanics[14, 15]. Numerically, the macroscale response can be e ffi ciently modeled within a nonlinear finite element frame-work. The standard coupling between the macroscale and microscale descriptions, however, entails the simulationof individual RVEs associated with each of the integration points of the finite element simulation. This strategy ofnesting DFN models within the finite element framework for large-scale heterogeneous structures is computationallyprohibitive. Not surprisingly, their use has been limited [16, 17].A myriad of strategies for model order reduction have been implemented to predict the mechanical propertiesof multiscale materials that balance computational cost and accuracy [18, 19, 20, 21, 22]. In particular, machinelearning (ML) methods have been extended to various problems of mechanics modeling and multiscale simulations,e.g., modeling macroscopic material behavior by stress-strain curves obtained from micromechanical simulations[17]; approximating the surface response and deriving the macroscopic stress and tangent tensor components [23];capturing the multiscale hydro-mechanical coupling e ff ect of porous media with pores of varying sizes [24], andpredicting the stress-strain behavior of binary composites [25].The neural network method is one example of a ML technique that is capable of extracting meaningful rela-tionships from data given a su ffi cient amount of training data. This method is e ffi cient in automatically discoveringand capturing the underlying complex high-dimensional mapping from the feature vector input to the desired outputwithout the need of manually deriving specific functional forms [26, 27].In this paper, we propose a novel multiscale modeling approach that employs fully connected neural networks(FCNN) to describe the homogenized response of simulations of DFNs. The FCNN is able to predict the mechanicalresponse of various microstructures under di ff erent loads. In this way, the trained FCNN represents a new macroscopicconstitutive relation. To develop the FCNN metamodel, we first generated training data by creating a large number ofDFN microstructures and subjecting the RVEs to various biaxial loading conditions. Derivatives of the strain energywith respect to the strain invariants were obtained by creating a Gaussian process (GP) surrogate for the strain energyfurther constrained on the observed stress data. Using the values of the derivatives of the strain energy as the data,we trained a FCNN model to predict the mechanical properties (e.g., stress and elasticity tensor) of a compositebased on its microstructure and the value of the strain invariants. Notably, the convexity of the strain energy, whichis essential to ensure the model is physically meaningful and numerically stable when the FCNN is used in finiteelement simulations, was considered in the training process by specifically defining the loss function to enforce theseconstraints. Following training, we tested the FCNN by comparing the prediction against the ground truth DFN modelon a validation set. Finally, we implemented the FCNN in a general User Material Subroutine (UMAT) in the popularfinite element package Abaqus.Our work demonstrates that neural networks can be trained by micromechanical simulations, which capture ECMnetwork behavior accurately. The main features of our model are the invariant formulation of the FCNN which satisfiesstress objectivity a priori , the imposition of a convexity constraint on the strain energy function for numerical stability,the e ff ective inference of the nonlinear map, and the implementation of arbitrary FCNN as a UMAT subroutine inAbaqus. We anticipate that this work will enable the widespread use of ML-driven multiscale simulations by reducingthe computational cost without compromising accuracy. Materials and Methods
Discrete fiber network models
In the DFN model, multiple fibers were linked to each other through cross-links, and they were allowed to stretch,compress, or rotate at these points. Bending forces were not considered during fiber deformation. The procedurefor the generation of the networks was in a similar manner to previous e ff orts [7, 28]. Namely, we seeded the three-dimensional RVE with random uniformly distributed nucleation points. Each point gave rise to two segments, which2 igure 1: A microscale network model (A), was used to simulate the fibrin networks (B). grew oppositely along a randomly chosen direction sampled from a uniform distribution on the unit sphere. At eachstep of fiber growth, every segment was checked to decide if it had exceeded the boundaries of the cube or hadcollided with another segment. At the point of the identified intersection or collision, a node was introduced, and thesegment stopped growing. When the growth of all the segments stopped, the algorithm of fibrillar network generationterminated. A sample image of a DFN is depicted in Figure 1A, alongside a confocal image of a 2mg / ml fibrin gel inFigure 1B.Applying forces or deformations to the boundaries of the RVE leads to the deformation of the fibers inside of thenetwork. Individual fibers were considered as hyperelastic and the state of mechanical equilibrium was determinedby finding the minimum of the total energy in the fibers. The energy in a fiber 𝑖 in the network is Ψ ( 𝑖 ) 𝑓 = 𝑘 𝑓 (cid:16) ( 𝜆 ( 𝑖 ) 𝑓 ) − (cid:17) , (1)where the stretch of the fiber 𝑖 is given by 𝜆 ( 𝑖 ) 𝑓 = || x ( 𝑖 ) − x ( 𝑖 ) |||| X ( 𝑖 ) − X ( 𝑖 ) || , (2) X ( 𝑖 ) , X ( 𝑖 ) denote the reference coordinates of the nodes making up the fiber, and x ( 𝑖 ) , x ( 𝑖 ) the location of the nodesafter deformation. The material parameter 𝑘 𝑓 is the sti ff ness of the fiber. The mechanical equilibrium is obtained bysearching for the deformed nodal coordinates of all the fibers such that the total energy (cid:205) 𝑖 Ψ ( 𝑖 ) 𝑓 is minimized. Theboundary conditions are imposed displacements on all nodes at the boundary of the RVE. The homogenized stress onthe RVE is calculated based on the averaging [7], 𝜎 𝑀 = (cid:104) 𝜎 𝜇 (cid:105) = 𝑉 ∫ 𝑅𝑉 𝐸 𝜎 𝜇 . (3) 𝜎 𝑀 denotes the macroscale stress, which is the average of the microscale stress field 𝜎 𝜇 over the RVE volume 𝑉 .For the DFN, the stress can be written in indicial notation as 𝜎 𝑀𝑖 𝑗 = 𝑉 ∑︁ 𝑥 𝑏𝑖 𝑅 𝑏𝑗 (4)where the sum is over the deformed coordinates of the boundary nodes, x 𝑏 , and the corresponding reaction forces R 𝑏 .To generate enough training data for the FCNN, a total of 1100 networks were generated. To explore the di ff erentnetwork geometries, two quantities of interest were identified as the best descriptors for microstructure: the volume3 igure 2: Framework for training a ML-based microscale model. Input data was generated by sampling from microstructure 𝜃, 𝜑 , and deformationparameters 𝜆 𝑥 , 𝜆 𝑦 , which were post-processed to obtain the invariants 𝐼 , 𝐼 . Multiple random geometries were sampled for a given microstructureand deformation. The RVE model was evaluated to generate output stress and strain energy, which were post-processed through Gaussian processregression and optimization to obtain energy derivative outputs (A). A fully connected neural network (FCNN) was trained on a training subsetof the data, and tested on a di ff erent subset of data. FCCN training was further constrained to satisfy convexity and symmetry, requirements for astable finite element implementation (B). fraction 𝜃 , and the fiber diameter 𝜑 . The range of the model parameters was obtained to span similar networks to thosereported in the literature [7, 29]. Namely, the diameter 𝜑 was in the range of [ − ] nm, and the volume fractionwas in the range 𝜃 ∈ [ . , . ] %.Following the generation of the 1100 di ff erent RVE microstructures, data of stress and energy for each of theseRVEs under various deformation regimes were generated. Analogous to mechanical tests performed to characterizesoft material mechanics, the RVEs were subjected to biaxial deformation under plane stress. Eleven stretches in the 𝑥 axis were selected 𝜆 𝑥 = [ . , . , . , . , . , . , . , . , . , . , . ] . Similarly for the 𝑦 axis,we used 𝜆 𝑦 = [ . , . , . , . , . , . , . , . , . , . , . ] . The deformations applied to theRVE were all the possible combinations of 𝜆 𝑥 and 𝜆 𝑦 from the values above, i.e. a total of 121 simulations per RVEwere run for the training of the ML metamodel. The transverse stretch was calculated to maintain incompressibility.While not used for training of the ML metamodel, uniaxial simulations were also explored to compare against previousreports in the literature. Hyperelastic material behavior at the macroscale
In the multiscale formulation, the average stress of the homogenized RVE, 𝜎 𝑀 = (cid:104) 𝜎 𝜇 (cid:105) , is a function of the averagedeformation gradient F 𝑀 = (cid:104) F 𝜇 (cid:105) . In practice, the average deformation gradient known from the macroscale, F 𝑀 , isimposed as a boundary condition to the RVE. Mechanical equilibrium of the RVE then yields the microscale field 𝜎 𝜇 ,which is averaged to compute the corresponding stress 𝜎 𝑀 . In this way, the relationship 𝜎 𝑀 ( F 𝑀 ) is implied althoughnot available explicitly, and it needs instead the evaluation of the RVE model. Under the additional assumption thatthe macroscopic material behavior is that of a nearly incompressible hyperelastic material [30, 31], we propose theexistence of a strain energy function Ψ Ψ = Ψ iso ( ¯ 𝐼 , ¯ 𝐼 ) + Ψ vol ( 𝐽 ) (5)with the standard definition of the isochoric strain invariants and the volume change¯ 𝐼 = tr ( ¯ b ) = ¯ b : I , ¯ 𝐼 = (cid:16) ( tr ( ¯ b )) − tr ( ¯ b ) (cid:17) , 𝐽 = det ( F 𝑀 ) . (6)Furthermore, the invariants in eq. 6 require the definition of the isochoric left Cauchy-Green deformation tensor,¯ b = 𝐽 − / b = 𝐽 − / F 𝑀 ( F 𝑀 ) 𝑇 . (7)4ote that the deformation metrics defined in eqs. (6) and (7) are with respect to the macroscale deformation. Thestrain energy function in (5) is not analytically available, but its existence implies that 𝜎 𝑀 = b 𝜕 Ψ 𝜕 b = P : ¯ 𝜎 + 𝑝 I = ¯ 𝜎 −
13 tr ( ¯ 𝜎 ) + 𝑝 I , (8)where we have introduced the fictitious stress and the pressure,¯ 𝜎 = b 𝜕 Ψ iso 𝜕 ¯ b , 𝑝 = 𝜕 Ψ vol 𝜕𝐽 . (9)Finally, to compute the fictitious stress in (9), the derivatives Ψ ( ¯ 𝐼 , ¯ 𝐼 ) = 𝜕 Ψ iso / 𝜕 ¯ 𝐼 and Ψ ( ¯ 𝐼 , ¯ 𝐼 ) = 𝜕 Ψ iso / 𝜕 ¯ 𝐼 are needed, and, like the strain energy, these derivatives are also functions of the isochoric invariants. Thus, knowledgeof the functions Ψ , Ψ , Ψ vol are enough to compute the stress at the macroscopic level 𝜎 𝑀 . It is remarked againthat these functions are not analytically available. The approach presented here relies on computing the stress andenergy of the fibers from solving the RVE model upon application of the macroscale deformation F 𝑀 as a boundarycondition. Then, a FCNN can be trained to represent the functions Ψ , Ψ that best represent the stress data from theRVE simulations. Strain energy and its derivatives using Gaussian processes and optimization
A Gaussian process (GP) is a stochastic process whose value at any collection of data points can be described witha multivariate normal distribution [32]. Consider a data set comprised of input-output pairs D = { x = ( 𝐼 , 𝐼 ) , Ψ 𝑓 } ,where 𝐼 and 𝐼 are the invariants of the corresponding deformation F 𝑀 applied to the RVE, and Ψ 𝑓 = (cid:205) Ψ ( 𝑖 ) 𝑓 is thetotal strain energy accumulated in the fibers of the RVE. The data set D consists of the vector of inputs and outputsfrom all the training simulations. We used a GP to describe these data. The prior of the GP is fully determined by amean function and a covariance function. A zero-mean prior was considered. For the covariance function we used aradial basis function (RBF) kernel [33]. This kernel is also commonly referred to as the squared exponential kernel.The implementation of the GP in the Python package GPy was used [34]. The hyperparameters of this kernel are thelength scale and signal strength, which were optimized to maximize the likelihood of the observed data pairs. Afterfitting, the posterior GP is denoted Ψ ∗ 𝑓 ( 𝐼 , 𝐼 ) , and it can be used to make predictions at points x ∗ = ( 𝐼 ∗ , 𝐼 ∗ ) that werenot part of the training data. The predicted strain energy Ψ ∗ 𝑓 at a new point x ∗ is a GP with mean function 𝜇 𝑓 = k ( x ∗ , x ) K ( x , x ) − Ψ 𝑓 (10)where k ( x ∗ , x ) is the vector of covariances between the new test point x ∗ and the training data points x . The matrix K ( x , x ) is the covariance between all training inputs. The vector Ψ 𝑓 contains all the observations of the strain energyof the fibers for all the inputs x . The posterior covariance is Σ 𝑓 = K ( x ∗ , x ∗ ) − k ( x ∗ , x ) K ( x , x ) − k ( x ∗ , x ) . (11)Given the posterior GP for the strain energy, Ψ ∗ 𝑓 , the derivatives of the strain energy function with respect to theinvariants are also GPs that can be easily computed from eqs. (12) and (11) [35]. The mean of the gradient of the GPis simply the derivative of the mean function (12), 𝜇 Ψ 𝑖 = k 𝑖 ( x ∗ , x ) K ( x , x ) − Ψ 𝑓 , (12)where k 𝑖 ( x ∗ , x ) is the vector with the derivatives of the kernel 𝜕𝑘 (( 𝐼 ∗ , 𝐼 ∗ , x )/ 𝜕𝐼 ∗ 𝑖 . The covariance of the derivativescan be obtained from the derivative of the kernel function, Σ Ψ 𝑖, Ψ 𝑗 = 𝜕 𝑘 ( x , x (cid:48) ) 𝜕𝐼 𝑖 𝜕𝐼 (cid:48) 𝑗 . (13)However, the GPs for the derivative functions defined with eq. (12) ignore the stress information and rely solelyon the energy information. Thus, we furthered constrained the derivatives of the strain energy based on the stressdata. In the previous section, the stress for a nearly incompressible material was derived.This description was the one5hat we eventually implemented in the finite element subroutine. In the RVE simulations incompressibility could beimposed exactly. For incompressible plane stress behavior the stress can be simpliefied to 𝜎 = b 𝜕 Ψ iso 𝜕 b + 𝑝 I (14)with det ( b ) = det ( F 𝑀 ) =
1, and where 𝑝 becomes a Lagrange multiplier that can be solved for to ensure the planestress condition 𝜎 𝑧𝑧 =
0. Expanding the derivative of the strain energy in (14), 𝜎 = Ψ ( 𝐼 , 𝐼 ) b + Ψ ( 𝐼 , 𝐼 )( 𝐼 b − b ) + 𝑝 I . (15)To summarize, from the set of RVE simulations a GP was constructed over the total strain energy of the fibers as afunction of the invariants of the applied deformation. This GP ignored the stress data. However, a relationship betweenthe stress and the derivatives of the strain energy is available as seen in eq. (15). Thus, an additional optimization wasperformed to correct the derivatives of the strain energy predicted by the GP in order to also satisfy (15). The pressureneeded for the plane stress condition in the RVE was also computed during this optimization step. The optimizationwas carried out with the minimize function in the optimization module of the Python package SciPy [36]. At the endof this step we had all the data as illustrated in Figure 2A. Fully connected neural network
Based on the previous steps, schematized in Figure 2A, the overall input data consisted of the variables { 𝜆 𝑥 , 𝜆 𝑦 , 𝐼 ,𝐼 , 𝜃, 𝜑 } , as well as multiple realizations of the networks even for the same { 𝜃, 𝜑 } . However, not all these inputs areindependent. Ultimately, the goal was to obtain a model that was a function of the deformation invariants. Yet, itwas easier to generate the training data by imposing the principal stretches. Thus, deformations on the RVE were im-posed by sampling the in-plane stretches { 𝜆 𝑥 , 𝜆 𝑦 } which were post-processed to get the invariants { 𝐼 , 𝐼 } . Similarly,the output consisted of the variables { 𝜎, 𝑝, Ψ , Ψ , Ψ } . For a finite element implementation, only prediction of thefunctions { Ψ , Ψ } is needed. Nevertheless, two FCNN were trained. A FCNN to predict the stress tensor directly asa function of the stretches 𝜆 𝑥 , 𝜆 𝑦 was trained. However, this approach imposes the choice of a cartesian coordinatesystem aligned with the principal in-plane directions for the strain and stress tensors. Therefore this network was notexplored in much detail. An alternative approach, favored in this work, is to train the FCNN to predict the derivativesof the strain energy, Ψ , Ψ , as function of the strain invaraints. This second approach was the one implementedin the commercial finite element software Abaqus. In addition to the advantage of an invariant-based formulation,training the FCNN on the strain energy and its derivatives enabled the imposition of meaningful constraints, namelythe convexity of the strain energy. The process of FCNN training is illustrated in Figure 2B.For the FCNN trained to predict Ψ and Ψ , a weighted combination of the mean squared error (MSE) loss, meanabsolute error (MAE) loss, and the mean absolute percentage error (MAPE) loss were used L MSE 𝑖 = 𝑁 𝑑𝑎𝑡𝑎 ∑︁ 𝑛 = (cid:16) Ψ ( 𝑛 ) 𝑖 − ˆ Ψ ( 𝑛 ) 𝑖 (cid:17) L MAE 𝑖 = 𝑁 𝑑𝑎𝑡𝑎 ∑︁ 𝑛 = (cid:12)(cid:12)(cid:12) Ψ ( 𝑛 ) 𝑖 − ˆ Ψ ( 𝑛 ) 𝑖 (cid:12)(cid:12)(cid:12) L MAPE 𝑖 = 𝑁 𝑑𝑎𝑡𝑎 ∑︁ 𝑛 = (cid:12)(cid:12)(cid:12)(cid:12)(cid:12) Ψ ( 𝑛 ) 𝑖 − ˆ Ψ ( 𝑛 ) 𝑖 Ψ ( 𝑛 ) 𝑖 (cid:12)(cid:12)(cid:12)(cid:12)(cid:12) (16)where 𝑖 = , 𝑛 denotes each one of the 𝑁 data pairs {( 𝜑 ( 𝑛 ) , 𝜃 ( 𝑛 ) , 𝐼 ( 𝑛 ) , 𝐼 ( 𝑛 ) ) , ( Ψ ( 𝑛 ) , Ψ ( 𝑛 )) )} , and ˆ Ψ ( 𝑛 ) 𝑖 is theprediction of the strain energy derivative by the FCNN for the corresponding input ( 𝜑 ( 𝑛 ) , 𝜃 ( 𝑛 ) , 𝐼 ( 𝑛 ) , 𝐼 ( 𝑛 ) ) . The lossesin (16) only penalize the error in the prediction of Ψ and Ψ . However, as mentioned previously, an advantage ofassuming a hyperelatic framework is that physically realistic constraints can also be considered during training of theFCNN. Indeed, a major driver for the development of analytical expressions for strain energy potentials of soft tissue6s that convexity can be analyzed in detail [37]. For the FCNN, the Hessian matrix of the strain energy function canbe computed numerically at any of the 𝑛 inputs used for training, H ( 𝑛 ) = (cid:18) Ψ Ψ Ψ Ψ (cid:19) = (cid:169)(cid:173)(cid:171) 𝜕 ˆ Ψ ( 𝑛 ) 𝜕𝐼 𝜕 ˆ Ψ ( 𝑛 ) 𝜕𝐼 𝜕 ˆ Ψ ( 𝑛 ) 𝜕𝐼 𝜕 ˆ Ψ ( 𝑛 ) 𝜕𝐼 (cid:170)(cid:174)(cid:172) . (17)The Hessian contains the second derivatives of the strain energy. On the other hand, since the FCNN output arealready the first derivatives Ψ , Ψ , only one derivative of the FCNN output is required . Computing the derivatives ofthe FCNN output with respect to the inputs is straightforward and follows an algorithm analogous to back-propagation[38]. The first obvious constraint on the FCNN output is that the Hessian must be positive. Therefore, one additionalcomponent of the loss is L sym = 𝑑𝑎𝑡𝑎 ∑︁ 𝑛 = (cid:12)(cid:12)(cid:12)(cid:12)(cid:12) 𝜕 ˆ Ψ ( 𝑛 ) 𝜕𝐼 − 𝜕 ˆ Ψ ( 𝑛 ) 𝜕𝐼 (cid:12)(cid:12)(cid:12)(cid:12)(cid:12) (18)The second constraint is the positive semi-definiteness of H which is enforced by the loss [39] L psd = 𝑁 𝑑𝑎𝑡𝑎 ∑︁ 𝑛 = (cid:12)(cid:12)(cid:12)(cid:12)(cid:12) 𝜕 Ψ ( 𝑛 ) 𝜕𝐼 (cid:12)(cid:12)(cid:12)(cid:12)(cid:12) + 𝑁 𝑑𝑎𝑡𝑎 ∑︁ 𝑛 = (cid:12)(cid:12)(cid:12)(cid:12)(cid:12) 𝜕 Ψ ( 𝑛 ) 𝜕𝐼 𝜕 Ψ ( 𝑛 ) 𝜕𝐼 − 𝜕 Ψ ( 𝑛 ) 𝜕𝐼 𝜕 Ψ ( 𝑛 ) 𝜕𝐼 (cid:12)(cid:12)(cid:12)(cid:12)(cid:12) (19)Hence, the final loss function of the convexity constraints can be obtained by adding the di ff erent losses, each witha given weight. The weights were manually adjusted to obtain desired accuracy and satisfaction of the constraints. Inparticular, the e ff ect of the weights of the symmetry constraint 𝛼 sym and of the positive semi-definiteness constraint 𝛼 psd were investigated in more detail in the Results section.An additional FCNN for the prediction of the principal stresses as a function of the principal stretches was alsoconstructed. However, since in the end this FCNN was not implemented into the finite element code, only the standardMAPE loss was used in training, which is a standard in machine learning [40, 41].The Adam optimization algorithm [42] was used to train the networks and update the FCNN weights duringback-propagation [43]. The initial learning rate was set as 0 . .
9. The training andvalidation data splits were 85% and 15% of the total data points. The total number of observations was 132 , ,
100 networks and the di ff erent deformations. A typical FCNN consists of an input layer, one or morehidden layers, and an output layer. Here we trained and evaluated an FCNN with 3 hidden layers of dimensions8, 8, and 16 respectively, and 1 output layer. The activation function used was ReLu [44](Table 1). The trainingwas implemented using Keras [45] with a TensorFlow [46] back-end on a hardware platform with the followingspecifications: Overclocking 5.0 GHz Intel i9 processor, 32 GB DDR4 / Table 1: NN architectures
Layer (type) Output shape inite element implementation
A general FCNN-based user material subroutine (UMAT) was implemented in the nonlinear finite element pack-age Abaqus Standard (Rhode Island, United States). Rather than hard-coding the specific FCNN that was trained onthe DFN data, the UMAT function was programmed to take in the network architecture and the weights and biasesas material parameters. Thus, the FCNN needs to be specified in the input file. As a result, the UMAT can be usedfor di ff erent FCNN architectures and trained on di ff erent data. In the UMAT function itself, the FCNN is constructedand evaluated based on the information of the input file. This approach gives flexibility. For example, in this pa-per two di ff erent regions in a finite element model were assigned di ff erent FCNNs to model their material behaviorwithout the need for multiple user subroutines. The deformation gradient F at an integration point is passed as oneof the arguments of the UMAT. Based on this deformation gradient, the corresponding isochoric left Cauchy Greendeformation, ¯ b , and its invariants ¯ 𝐼 , ¯ 𝐼 , are computed. Then, using the microstrucure parameters 𝜑 and 𝜃 which aredefined as material properties, together with the network architecture and the weights and biases of the network (alsodefined as material properties), the FCNN evaluates Ψ , Ψ . Evaluation of neural networks is a combination of simpleoperations. Consider a dense layer, layer 𝑖 , of the network. Layer 𝑖 takes in 𝑚 inputs from the previous layer stored invector y 𝑖 − and produces 𝑛 outputs y 𝑖 . In other words, 𝑚 is the number of neurons in layer 𝑖 − 𝑛 is the number ofneurons in the layer 𝑖 . The weights for this layer are arranged in the matrix W 𝑖 which is 𝑚 × 𝑛 . The vector of biasesfor layer 𝑖 is denoted B 𝑖 . Propagation of information in this layer corresponds to the linear map y 𝑖 = W 𝑇𝑖 y 𝑖 − + B 𝑖 . (20)The nonlinearity in the network comes from the activation layers. For the activation layer 𝑖 , the output is simply y 𝑖 = 𝑔 ( y 𝑖 − ) (21)where 𝑔 (•) is the activation function applied element-wise to the vector y 𝑖 − . The ReLu activation functionused here is 𝑔 ( 𝑦 ) = 𝑚𝑎𝑥 ( , 𝑦 ) . In addition to computing the functions Ψ , Ψ , their derivatives with respect to theinvariants ¯ 𝐼 , ¯ 𝐼 are needed for the elasticity tensor, which is essential for the convergence of the Newton iterations ofthe finite element solver. For a dense layer 𝑖 , the input gradient is 𝐽 𝑖 − of dimension 𝑚 × 𝑛 input , with 𝑛 input the numberof inputs to the whole FCNN. The gradient output of a dense layer 𝑖 is J 𝑖 = W 𝑇𝑖 J 𝑖 − . (22)For the activation layer, the gradient output is J 𝑖 = diag ( 𝑔 (cid:48) ( y 𝑖 )) J 𝑖 − . (23)where diag (•) denotes a diagonal matrix, and 𝑔 (cid:48) (•) is the derivative of the activation function. The gradient forthe input layer is initialized as the identity matrix of size 𝑛 input × 𝑛 input . After doing the operations defined in (22)and (23), the gradient output is a matrix of size 𝑛 output × 𝑛 input . For the particular FCNN used in this paper, the finalgradient output contains the Hessian (17) as a sub-matrix.After evaluation of the FCNN, the prediction of Ψ and Ψ are used to compute the isochoric stress based onequations (8) and (9). Additionally, a volumetric strain energy is required, which is not part of the FCNN prediction.Here we opt for a classical expression for the volumetric part of the strain energy Ψ vol = 𝐾 ( 𝐽 − ) (24)where 𝐾 is the bulk modulus and 𝐽 = det ( F ) is the volume change.The consistent tangent can be derived independently of the material model, simply by knowing the second deriva-tives of the strain energy function. The FCNN outputs the prediction of the second derivatives based on the operations822) and (23) . The following terms are used for ease of notation of the consistent tangent 𝛿 = ( Ψ + Ψ + 𝐼 Ψ ) 𝛿 = − ( Ψ + ¯ 𝐼 Ψ ) 𝛿 = Ψ 𝛿 = − Ψ (25)where Ψ , Ψ , Ψ , Ψ are the derivatives predicted by the FCNN that correspond to the entries of the Hessian,eq. (17). The isochoric part of the elasticity tensor is 𝐽 c 𝑚𝑎𝑡ℎ𝑟𝑚𝑖𝑠𝑜 = 𝛿 (cid:18) ¯b ⊗ ¯b −
13 ¯ 𝐼 ( ¯b ⊗ I + I ⊗ ¯b ) +
19 ¯ 𝐼 I ⊗ I (cid:19) + 𝛿 (cid:18) ˆb ⊗ ¯b + ¯b ⊗ ¯b −
13 tr ( ¯b )( ¯b ⊗ I + I ⊗ ¯b ) −
13 ¯ 𝐼 ( ¯b ⊗ I + I ⊗ ¯b ) +
29 ¯ 𝐼 tr ( ¯b ) I ⊗ I (cid:19) + 𝛿 (cid:18) ¯b ⊗ ¯b −
13 tr ( ¯b )( ¯b ⊗ I + I ⊗ ¯b ) + (cid:16) tr ( ¯b ) (cid:17) I ⊗ I (cid:19) + 𝛿 (cid:18) ¯b ¯ ⊗ ¯b − ( ¯b ⊗ I + I ⊗ ¯b ) + ( ¯b : ¯b ) I ⊗ I (cid:19) (26)where the following special dyadic product was used (•) 𝑖 𝑗 ¯ ⊗(◦) 𝑘𝑙 = ( ★ ) 𝑖𝑘 𝑗𝑙 . The volumetric contribution of theelasticity tensor is c vol = (cid:18) 𝑝 + 𝐽 𝜕 𝑝𝜕𝐽 (cid:19) I ⊗ I − 𝑝 I ¯ ⊗ I . (27)The isochoric and volumetric components of the spatial elasticity tensor defined in eqs. (26) and (27) are part ofthe consistent tangent needed for Abaqus’ Newton-Raphson solver. Mechanical equilibrium is solved in Abaqus byincrementally integrating the stress using the Jaumann stress rate. As a result, the consistent tangent needs a correctionwith respect to the standard form of the elasticity tensor [48]. The tangent for Abaqus takes the form c abaqus = c iso + c vol + ( 𝜎 ⊗ I + 𝜎 ⊗ I + I ⊗ 𝜎 + I ⊗ 𝜎 ) (28)in which an additional tensor product is defined as (•) 𝑖 𝑗 ¯ ⊗(◦) 𝑘𝑙 = ( ★ ) 𝑖𝑙 𝑗𝑘 . In summary, the architecture, weights,and biases of the FCNN are passed to the UMAT subroutine as material parameters, and the FCNN is evaluated giventhe current deformation gradient at an integration point to output Ψ , Ψ and their derivatives Ψ , Ψ , Ψ , Ψ .These functions are used to compute the stress 𝜎 and the consistent tangent c abaqus . Results
Response under uniaxial and biaxial loading
We first compared that the results from our DFN implementation matched observations from previous computa-tional models of fiber networks as well as experimental data on fibrin gels. We created 110 networks with a volumefraction 𝜃 ≈ . , 𝜃 ∈ ( . , . ) and fiber diameter 𝜑 = . 𝜃 and 𝜑 because these variablesdescribe the homogenized network properties but do not specify a network uniquely.The 110 RVEs were subjected to uniaxial deformation following similar approaches in the literature [49, 29, 7].. The nonlinear relationship of stress as a function of stretch was observed, with increasing slope at higher stretch9 niaxial test Convergence of DFN simulation Stress vs. Microstructure 𝜆 " 𝜎 " ( k P a ) 𝜎 " ( k P a ) 𝑛 %&’ 𝜑 ( µ m ) 𝜃 ( % ) A B C
Figure 3: Verification of DFN response under uniaxial deformation. In uniaxial tests (A), 110 networks with homogenized microstructure properties 𝜃 ≈ . 𝜑 = 𝜃 , 𝜑 , the random natureof the networks leads to uncertainty in the predicted stress 𝜎 𝑥 , shown as a shaded area (A). The variance in the stress 𝜎 𝑥 was dependent on thenumber of degrees of freedom of the network, defined as the percentage of inner nodes with respect to total number of nodes (B). The stress inuniaxial tests was a function of the homogenized microstructure properties 𝜑 and 𝜃 for a fixed value of fiber sti ff ness, with increasing stress forsmaller fiber diameters and relatively little influence of the volume fraction (C). values. At a stretch of 𝜆 𝑥 = .
2, for example, the mean ( + - S.D.) stress was 1.45 kPa + - 0.33 kPa (Figure 3A). Thestress showed variance, suggesting that even networks with the same homogenized microstructure properties 𝜑, 𝜃 , canhave slightly di ff erent mechanical behavior. To further explore this uncertainty, we plotted the variance in the stress asa function of the relative number of degrees of freedom, i.e. the number of inner nodes with respect to the total numberof nodes (Figure 3B). Decreased variance with higher number of degrees of freedom was observed. Variation of themicrostructure variables 𝜃, 𝜑 , led to consistent trends in the stress response under uniaxial deformation: decreasingfiber diameter led to increase in stress, almost independently of volume fraction (Figure 3C). We note that the sti ff ness 𝑘 𝑓 was kept constant (0.02 kPa), independent of fiber diameter.The response of the networks under biaxial deformations was also explored. Figure 4 shows the response undero ff -biaxial testing keeping 𝜆 𝑦 =
1. For the 110 networks with 𝜃 = . , 𝜑 = Ψ , showed a similar trend to the uniaxial tests. The strain energy increased nonlinearly for increasing stretch,with some variance attributed to the randomness of the 110 networks even though they shared the same homogenizedproperties (Figure 4A). We also explored the e ff ect of the number of degrees of freedom on the variance of thestrain energy. A slight decrease in the variance of Ψ could be found as the number of degrees of freedom increased(Figure 4B), although this was less obvious compared to the stress response. Analyzing the mechanical behaviorwith respect to the microstructure variables, it can be seen that low volume fraction and smaller diameters resulted inhigher energy density (Figure 4B). The derivatives, Ψ and Ψ , showed a slightly nonlinear response with a relativelyconstant variance with respect to the deformation. The variance in the derivatives did not show a strong dependenceon the number of degrees of freedom (Figure 4E,H). Analyzing the contours of the functions Ψ and Ψ , it can be seenthat they depend on the microstructure variables 𝜃 and 𝜑 in a similar way to the stress (Figure 4F,I), which is expectedsince these derivatives appear linearly in the definition of the stress tensor, see eq. (15). Note that for a neo-Hookeanmaterial, the vaue of Ψ will be a constant and the value of Ψ will be zero. For a Mooney-Rivlin material, the valueof both Ψ and Ψ will be non-zero constants. Our results show that the DFN networks are sightly more nonlinearthat these simple but commonly used strain energy functions. Instead of trying to find an analytical expression witha higher degree of nonlinearity to match the data, the next section explains how the FCNN was trained to predict themechanical response without the need of classical strain energy functions. FCNN training and validation
A total of 132,000 data points were used to train the FCNN with the architecture specified in Table 1. Traditionally,a single loss function is used for the optimization of the weights and biases of the FCNN. However, as outlined above,to obtain strain energies that satisfied convexity with respect to the deformation invariants, multiple loss functionswere linearly combined. In particular, there was a trade-o ff between the weights assigned to the loss functions thatdrove the accuracy of the predictions and the weights assigned to the loss functions that drove the convexity andsymmetry constraints. A grid search was performed to select the weights 𝛼 psd and 𝛼 sym , which enforced the positive10 K J / m " ) ( K J / m " ) Biaxial test Convergence of DFN simulation Strain energy density vs. Microstructure
A B C 𝜓 ( k J / m " ) 𝜆 & 𝜑 ( µ m ) 𝜓 ( k J / m " ) 𝑛 *+, 𝜃 ( % ) 𝜓 / ( k J / m " ) 𝜆 & 𝜓 / ( k J / m " ) 𝑛 *+, 𝜑 ( µ m ) 𝜃 ( % ) D E F 𝜓 ( k J / m " ) 𝜓 ( k J / m " ) 𝜆 & 𝑛 *+, 𝜃 ( % ) 𝜑 ( µ m ) G H I
Figure 4: Response of DFNs under o ff -biaxial loading. The strain energy Ψ (A), and its derivatives Ψ , Ψ show a nonlinear response with respectto increasing deformation (D, G). The strain energy variance depends on the number of degrees of freedom (B), while the derivatives are relativelyconstant even as the number of degrees of freedom increases (E, H). In response of changes in microstructure, the strain energy in the fibersincreases with both volume fraction 𝜃 and fiber diameter 𝜑 (C), while the derivatives Ψ and Ψ are mostly insensitive to the volume fraction butincrease with fiber diameter (F, I). ff ect of these parameters. The results for grid search on the accuracy with respect to the change in these weightsis shown in Table 2. This table clearly shows the trade o ff between the accuracy and the convexity of the model.Once 𝛼 psd and 𝛼 sym were large enough, further increase only produced minor changes in the convexity loss. The setof weights 𝛼 psd = 𝛼 sym = was selected since it yielded an overall good validation accuracy while satisfying theconvexity condition. Finally, a model with a total loss of 70.26 was achieved. The MAPE loss for Ψ was 26.05, theMSE loss for Ψ was 0.91, and the MAE loss for Ψ was 0.50. For Ψ , due to predictions close to zero, only theweight associated to the MAE loss was considered, and the value of this loss function for Ψ was 0.38 in the finalmodel. The positive definiteness loss was 1.35E-10, and the symmetry loss 2.32E-05. Table 2: FCNN performance for the grid search over the parameters 𝛼 sym , 𝛼 psd . 𝛼 sym 𝛼 psd (cid:205) L L
MAPE1 L MSE1 L MAE1 L MAE2 L psd1 L sym1
0. 1. 49.41 21.76 0.20 0.33 0.24 4.39 3.441. 1. 57.81 24.59 0.64 0.38 0.28 0.30 0.2410
1. 67.97 31.30 2.20 0.67 0.27 8.39 × − × − × − × − × − × − × − × − × − × − A second FCNN was trained, using the same overall data, but in this case the inputs were the principal stretches 𝜆 𝑥 and 𝜆 𝑦 , and the outputs were the components of the stress tensor 𝜎 . This FCNN was not used in the finite elementimplementation, but it was useful for evaluating the ability of the FCNN approach to capture the mechanical responseof the DFNs. For the FCNN that was trained to predict the stress tensor, only the MAPE loss function was used,as described in the Methods section. No constraints were considered in this case. The final loss for this FCNN was12.72%.As seen in the uniaxial and biaxial response tests (Figures 3 and 4), even for the same microstructure variables 𝜃 and 𝜑 , there is variance in the stress, the strain energy, and the derivatives of the strain energy. The FCNN trained onstress data was first compared to strip-biaxial data from 60 DFNs with the equivalent microstructure 𝜃 = . , 𝜑 = 𝜃 = . , 𝜑 = ff erent DFNs. The loss functions were evaluated with respect to theraw DFN data, i.e. without taking into account the variance in the DFN response. Thus, the FCNN converged naturallytoward the mean response with respect to the microstructure variables even though it was not trained directly on themean data from multiple RVEs with the same 𝜃 , 𝜑 . To explore the error as a function of microstucture, additionaltesting data consisting of 1000 simulations with di ff erent combinations of 𝜃 and 𝜑 were evaluated under strip-biaxialloading. The results are shown in Figure 5B. In this case, the prediction of the FCNN was compared against at least20 fiber networks sharing equivalent 𝜃, 𝜑 pairs that were not a part of the training data.The FCNN that was trained directly on the derivative data Ψ and Ψ was first compared to 100 DFNs with theequivalent microstructure 𝜃 = . , 𝜑 = 𝜃 and 𝜑 , the model predicted Ψ , Ψ with small errors for most networks (Figure 5D,F).12 "" ( k P a ) 𝜆 " Accuracy of FCNN in Biaxial test Error as a function of microstructure 𝜑 ( µ m ) 𝜃 ( % ) 𝜓 * ( k J / m + ) 𝜆 " 𝜑 ( µ m ) 𝜃 ( % ) 𝜓 , ( k J / m + ) 𝜆 " 𝜑 ( µ m ) 𝜃 ( % ) A BC DE F 𝜓 , 𝜓 * 𝜀 . ( % )( k J / m + ) Figure 5: Performance of the FCNN against the testing set. Performance of a neural network for the prediction of stress in strip biaxial loading(solid red line) for showed agreement with the mean (dashed blue line) and confidence interval (shaded blue region) of 60 DFNs with the equivalentmicrostructure microstructure 𝜃 = . , 𝜑 = ff ected by the number of fibers in the network. Loss function values for Ψ and Ψ did not show a distinctive trend as a function of microstructure(D,F). rror as a function of 𝜃 and 𝜓 𝜓 $ % & ’ ( k J / m ) ) 𝜃 ( % ) 𝜓 M A P E ( % ) 𝜓 $ % & ’ ( k J / m ) ) 𝜓 Error as a function of φ and 𝜓 𝜑 ( µm ) Error as a function of 𝜃 and 𝜓 .$%&’ 𝜓 . $ % & ’ ( k J / m ) ) 𝜃 ( % ) 𝜓 . $ % & ’ ( k J / m ) ) Error as a function of φ and 𝜓 .$%&’ 𝜑 ( µm ) 𝜓 . 𝜓 . A BC D 𝜓 M A P E ( % ) 𝜓 . M A E ( k J / m ) 𝜓 . M A E ( k J / m ) Figure 6: Performance of the FCNN as a function of microstructure. Predictions of the derivative function Ψ with respect to the volume fraction 𝜃 (A) and the fiber diameter 𝜑 (B). The error of the second derivative prediction Ψ showed similar trends with respect to 𝜃 and 𝜑 (C,D). B Figure 7: Strain energy density function with respect to 𝐸 and 𝐸 with 𝐸 = 𝜃 = . , 𝜑 = 𝐸 and 𝐸 with 𝐸 = Figure 6 further shows the dependence of the error on the microstructure while also showing the dependence on thetrue value of the derivative. Higher errors were seen on the edge of the plot, namely, for those networks with extremevolume fractions and fiber diameters. While most of the values for Ψ , Ψ were small and concentrated within anarrow range, some outlier values were noticed at the ends of the ranges for volume fraction and fiber diameter(Figure 6). As expected, the FCNN metamodel predictions were better on the region of the output space closer to thecenter of the range spanned by the true Ψ and Ψ .To visualize the convexity of strain energy density function from our FCNN, the strain energy density Ψ was com-puted over the strain space 𝐸 , 𝐸 , i.e. biaxial loading, as well as in the space 𝐸 , 𝐸 , which includes shear. TheFCNN outputs the derivatives Ψ , Ψ which were used to determine the second Piola Kirchho ff stress and integratedover strain curves ∫ S : 𝑑 E . The resulting strain energy function contours computed in this manner were plotted andexamined for a represented microstructure 𝜃 = . 𝜑 = Finite element simulation
Figure 8 shows results from three basic tests done to verify the correct implementation of the FCNN in the UMATsubroutine. We first performed a uniaxial test in which symmetry for X, Y, and Z planes was imposed on the corre-sponding boundaries. X-displacement was applied on one end of the domain to impose an overall stretch of 𝜆 𝑥 = . 𝜑 = . 𝜇 m and 𝜃 = . 𝜎 = .
34 kPa. Secondly, a uniaxial deformation for a higher stretch of 𝜆 𝑥 = . ff ects are visible. The center point of the mesh still behaves correspondingly to uniaxial loading, with astress that aligns with FCNN predictions. This simulation further suggests that the UMAT was implemented correctly.One last test was performed to verify that the simulation converged under di ff erent loading, in particular one that in-volves shearing. In this third simulation, the left end was kept clamped, but a shear deformation was added on theright boundary. All other surfaces were still traction free. The simulation converged also without a problem, and theresulting stress distribution showed a band of shear and normal stress along the diagonal, in agreement with the overallpattern of stress that would be expected from popular strain energy functions such as the neo-Hookean model. Theseresults suggest that the FCNN was implemented correctly in the UMAT, that the second derivatives obtained throughback-propagation are also implemented correctly in the UMAT, and that the potential satisfies convexity leading tothe proper convergence of the simulation. 15 igure 8: Finite element results using the FCNN in the UMAT subroutine. Uniaxial simulation with symmetry boundary conditions on the leftand bottom boundaries leads to a uniform stress field (A). Clamping the left boundary leads to a non-uniform stress distribution; the center of thedomain still behaves as in uniaxial loading (B). Shear simulations lead to a band of shear stress (C). igure 9: Finite element simulation of a rectangular domain with an inner region made out of a material with microstructure 𝜑 = . 𝜇 m and 𝜃 = .
05% and a surrounding material defined by 𝜑 = . 𝜇 m and 𝜃 = . One final example was the modeling of a problem more representative of an application of this kind of materialmodel. We modeled a quarter of a domain in which the center was made out of a di ff erent material with respect to thesurrounding domain. This type of domain could represent a wound or a heterogeneous thrombus depending on therelative material properties. To model only a quarter of the domain, symmetry boundary conditions were used at thebottom and left surfaces, while the top surface has only Y fixed. The right surface was subjected to a uniaxial stretch.The UMAT was used for both domains, but with di ff erent inputs for the microstructure variables. For the surroundingmaterial we used 𝜑 = . 𝜇 m and 𝜃 = . 𝜑 = . 𝜇 m and 𝜃 = . ff er material is higher compared to the inner soft inclusion. Discussion
In this study, we proposed to use a FCNN metamodel to replace a microscale fiber network model. The FCNNwas trained on network microstructures subjected to a wide set of deformations. While the microscale model takesapproximately 1 min to solve for equilibrium of a DFN , the prediction from the FCNN is obtained in less than1msec, highlighting the significant speed-up from using machine learning tools to interpolate physics-based models.Significantly, the FCNN model has a simple structure with a small number of parameters and the metamodel canbe constrained to satisfy physically meaningful equations, yielding an e ffi cient and numerically stable finite elementimplementation. Moreover, the FCNN was trained on many di ff erent networks that span di ff erent microstructures interms of volume fraction 𝜃 and fiber diameter 𝜑 . These constitutive parameters are physically relevant and can bedefined by imaging experimentally-generated fibrin gels, when available. One version of the metamodel was trainedto predict the stress as a function of the principal stretches and performed well in terms of relative error computedwith respect to a testing data set. However, in order to achieve an invariant-based formulation, a FCNN to directlypredict the derivatives of the strain energy with respect to the deformation invariants, Ψ and Ψ , was also trained.The error of the FCNN against the validation data was due in part to the inherent uncertainty in the response ofthe DFNs. Even when the homogenized microstructure variables 𝜃, 𝜑 , are fixed, this does not mean that the networksare fully determined. In practice, the networks were generated by seeding a domain with random nucleation pointsfor fiber formation. Hence, this source of randomness leads to variance in the mechanical response for a given pair 𝜃, 𝜑 . The type of FCNN used here is a deterministic model. When trained on the data from the di ff erent DFNs, theFCNN naturally learns the mean response for a given microstructure. Future work will use Bayesian neural networksthat can capture the variance in the response instead of deterministic predictions [50, 51]. For those samples with17elatively larger error, the networks had extreme values of microstructure parameters. Moreover, some values of Ψ , Ψ in the validation set also had outliers at the edges of the parameter range. These unusual and inconsistent valuescompared to the rest of the training and validation data could be caused by the existence of non-unique solutions tothe problem of fitting a GP and optimizing it to produce the derivative values. Other refinements of the DFN modelcan be considered. In this e ff ort, we opted for a mechanical equilibrium problem in which fibers contribute only dueto stretching but not bending. To make the networks more stable, bending resistance could be considered [52, 53].Di ff erent constitutive models for individual fibers can also be incorporated [54, 9, 55, 1]. Another source of error in theFCNN prediction was the imposition of the convexity and symmetry constraints. The FCNN without this constraintmore closely followed the data, but had small fluctuations in the response that led to loss of positive-definitenesslocally. The constraint loss yielded a FCNN that satisfied these requirements, and while the predictions were wellwithin the confidence intervals of the validation set, they appeared to be less nonlinear than the true values. Therefore,in addition to the Bayesian approach to handle the variance in the response function, networks with a more advancedstructure and activation functions can also be considered in the future to improve the training and testing performancewhile still satisfying the symmetry and convexity requirements.Due to the substantial reduction in computational time and good accuracy in the energy prediction, we embeddedour FCNN as a user material soubroutine in the commercial finite element package Abaqus. The finite elementimplementation showed that the FCNN led to stably convergent simulations. An advantage of our implementationis that the entire FCNN can be passed in a vector of material parameters defined in the input file. In this way, theUMAT does not evaluate the specific FCNN trained with the microscale mechanics data and can actually be used toevaluate any FCNN specified in the input file. For example the recent work by [56] shows that a FCNN can be used tointerpolate heart valve biaxial data. The metamodel in [56] was not used in finite element simulations, however, dueto the generality of our UMAT, the FCNN in [56] could be easily evaluated within our UMAT. All the files for thismanuscript are available online. Conclusions
Our work demonstrates that neural networks can be trained by micromechanical simulations, which capture ECMnetwork behaviors and their relation to microstructural variables such as volume fraction and fiber diameter. Thecomputational e ffi ciency of this kind of machine learning metamodel makes it suitable for implementation withinfinite element simulations. More importantly, the FCNN relies on the data of the microscale model and not on anysort of analytical approximation. Future work includes using the FCNN to predict rate-dependent constitutive relationsand for inverse estimation of microstructure parameters from a combination of mechanical tests and imaging data. Acknowledgements
This work was supported in part by the National Science Foundation under grant No. 1911346-CMMI to PIsTepole and Calve and the Bilsland Dissertation Fellowship to Dr. Yue Leng.
References [1] E. A. Sander, T. Stylianopoulos, R. T. Tranquillo, V. H. Barocas, Image-based biomechanics of collagen-based tissue equivalents, IEEEengineering in medicine and biology magazine 28 (2009) 10–18.[2] V. K. Lai, C. R. Frey, A. M. Kerandi, S. P. Lake, R. T. Tranquillo, V. H. Barocas, Microstructural and mechanical di ff erences between digestedcollagen–fibrin co-gels and pure collagen and fibrin gels, Acta biomaterialia 8 (2012) 4031–4042.[3] A. Mol, M. I. van Lieshout, C. G. Dam-de Veen, S. Neuenschwander, S. P. Hoerstrup, F. P. Baaijens, C. V. Bouten, Fibrin as a cell carrier incardiovascular tissue engineering applications, Biomaterials 26 (2005) 3113–3121.[4] G. SChlag, H. Redl, M. Turnher, H. Dinges, The importance of fibrin in wound repair, in: Fibrin sealant in operative medicine, Springer,1986, pp. 3–12.[5] N. Laurens, P. d. Koolwijk, M. De Maat, Fibrin structure and wound healing, Journal of Thrombosis and Haemostasis 4 (2006) 932–939.[6] Y. Li, H. Meng, Y. Liu, B. P. Lee, Fibrin gel as an injectable biodegradable sca ff old and cell carrier for tissue engineering, The ScientificWorld Journal 2015 (2015).[7] T. Stylianopoulos, V. H. Barocas, Volume-averaging theory for the study of the mechanics of collagen networks, Computer methods inapplied mechanics and engineering 196 (2007) 2981–2990.
8] Y. Li, W. Chen, H. Xu, X. Jin, 3D representative volume element reconstruction of fiber composites via orientation tensor and substructurefeatures, Technical Report, Ford Motor Company, 2016.[9] M. Hadi, E. Sander, J. Ruberti, V. H. Barocas, Simulated remodeling of loaded collagen networks via strain-dependent enzymatic degradationand constant-rate fiber growth, Mechanics of materials 44 (2012) 72–82.[10] B. Agoram, V. H. Barocas, Coupled macroscopic and microscopic scale modeling of fibrillar tissues and tissue equivalents, J. Biomech. Eng.123 (2001) 362–369.[11] N. J. Driessen, C. V. Bouten, F. P. Baaijens, A structural constitutive model for collagenous cardiovascular tissues incorporating the angularfiber distribution (2005).[12] M. G. Geers, V. G. Kouznetsova, W. Brekelmans, Multi-scale computational homogenization: Trends and challenges, Journal of computa-tional and applied mathematics 234 (2010) 2175–2182.[13] J. Fish, R. Fan, Mathematical homogenization of nonperiodic heterogeneous media subjected to large deformation transient loading, Inter-national Journal for numerical methods in engineering 76 (2008) 1044–1064.[14] A. N. Natali, E. L. Carniel, P. G. Pavan, P. Dario, I. Izzo, Hyperelastic models for the analysis of soft tissue mechanics: definition ofconstitutive parameters, in: The First IEEE / RAS-EMBS International Conference on Biomedical Robotics and Biomechatronics, 2006.BioRob 2006., IEEE, pp. 188–191.[15] A. D. Freed, D. R. Einstein, M. S. Sacks, Hypoelastic soft tissues, Acta mechanica 213 (2010) 205–222.[16] M. Alber, A. B. Tepole, W. R. Cannon, S. De, S. Dura-Bernal, K. Garikipati, G. Karniadakis, W. W. Lytton, P. Perdikaris, L. Petzold,et al., Integrating machine learning and multiscale modeling—perspectives, challenges, and opportunities in the biological, biomedical, andbehavioral sciences, NPJ digital medicine 2 (2019) 1–11.[17] D. Reimann, K. Chandra, N. Vajragupta, T. Glasmachers, P. Junker, A. Hartmaier, et al., Modeling macroscopic material behavior withmachine learning algorithms trained by micromechanical simulations, Frontiers in Materials 6 (2019) 181.[18] S. Bhattacharjee, K. Matouˇs, A nonlinear manifold-based reduced order model for multiscale analysis of heterogeneous hyperelastic materi-als, Journal of Computational Physics 313 (2016) 635–653.[19] P. Kerfriden, O. Goury, T. Rabczuk, S. P.-A. Bordas, A partitioned model order reduction approach to rationalise computational expenses innonlinear fracture mechanics, Computer methods in applied mechanics and engineering 256 (2013) 169–188.[20] Z. Liu, J. A. Moore, S. M. Aldousari, H. S. Hedia, S. A. Asiri, W. K. Liu, A statistical descriptor based volume-integral micromechanicsmodel of heterogeneous material with arbitrary inclusion shape, Computational Mechanics 55 (2015) 963–981.[21] J.-C. Michel, P. Suquet, A model-reduction approach in micromechanics of materials preserving the variational structure of constitutiverelations, Journal of the Mechanics and Physics of Solids 90 (2016) 254–285.[22] J. Oliver, M. Caicedo, A. E. Huespe, J. Hern´andez, E. Roubin, Reduced order modeling strategies for computational multiscale fracture,Computer Methods in Applied Mechanics and Engineering 313 (2017) 560–595.[23] B. Le, J. Yvonnet, Q.-C. He, Computational homogenization of nonlinear elastic materials using neural networks, International Journal forNumerical Methods in Engineering 104 (2015) 1061–1084.[24] K. Wang, W. Sun, A multiscale multi-permeability poroplasticity model linked by recursive homogenizations and deep learning, ComputerMethods in Applied Mechanics and Engineering 334 (2018) 337–380.[25] C. Yang, Y. Kim, S. Ryu, G. X. Gu, Prediction of composite microstructure stress-strain curves using convolutional neural networks, Materials& Design 189 (2020) 108509.[26] G. C. Peng, M. Alber, A. B. Tepole, W. R. Cannon, S. De, S. Dura-Bernal, 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.[27] A. B. Tepole, D. Nordsletten, K. Garikipati, E. Kuhl, Special issue on uncertainty quantification, machine learning, and data-driven modelingof biological systems, CMAME 362 (2020) 112832.[28] D. S. Clague, R. J. Phillips, A numerical calculation of the hydraulic permeability of three-dimensional disordered fibrous media, Physics ofFluids 9 (1997) 1562–1572.[29] V. K. Lai, S. P. Lake, C. R. Frey, R. T. Tranquillo, V. H. Barocas, Mechanical behavior of collagen-fibrin co-gels reflects transition fromseries to parallel interactions with increasing collagen content, Journal of biomechanical engineering 134 (2012).[30] L. A. Mihai, A. Goriely, Finite deformation e ff ects in cellular structures with hyperelastic cell walls, International Journal of Solids andStructures 53 (2015) 107–128.[31] L. A. Mihai, H. Wyatt, A. Goriely, A microstructure-based hyperelastic model for open-cell solids, SIAM Journal on Applied Mathematics77 (2017) 1397–1416.[32] M. Seeger, Gaussian processes for machine learning, International journal of neural systems 14 (2004) 69–106.[33] B. Sch¨olkopf, K. Tsuda, J.-P. Vert, Kernel methods in computational biology, MIT press, 2004.[34] S. Group, et al., She ffi eldml / gpy: Gaussian processes framework in python, 2012.[35] A. L. Frankel, R. E. Jones, L. P. Swiler, Tensor basis gaussian process models of hyperelastic materials, Journal of Machine Learning forModeling and Computing 1 (2020).[36] P. Virtanen, R. Gommers, T. E. Oliphant, M. Haberland, T. Reddy, D. Cournapeau, E. Burovski, P. Peterson, W. Weckesser, J. Bright, et al.,Scipy 1.0: fundamental algorithms for scientific computing in python, Nature methods 17 (2020) 261–272.[37] G. Chagnon, M. Rebouah, D. Favier, Hyperelastic energy densities for soft biological tissues: a review, Journal of Elasticity 120 (2015)129–160.[38] V. Avrutskiy, Backpropagation generalized for output derivatives, arXiv preprint arXiv:1712.04185 (2017).[39] J. E. Prussing, The principal minor test for semidefinite matrices, Journal of Guidance, Control, and Dynamics 9 (1986) 121–122.[40] A. Botchkarev, Performance metrics (error measures) in machine learning regression, forecasting and prognostics: Properties and typology,arXiv preprint arXiv:1809.03006 (2018).[41] B. L. Bowerman, R. T. O’Connell, A. B. Koehler, Forecasting, time series, and regression: an applied approach [cd] (2005).[42] D. P. Kingma, J. Ba, Adam: A method for stochastic optimization, arXiv preprint arXiv:1412.6980 (2014).[43] D. E. Rumelhart, G. E. Hinton, R. J. Williams, Learning representations by back-propagating errors, nature 323 (1986) 533–536.
44] V. Nair, G. E. Hinton, Rectified linear units improve restricted boltzmann machines, in: ICML.[45] F. Chollet, et al., Keras: The python deep learning library, ascl (2018) ascl–1806.[46] M. Abadi, A. Agarwal, P. Barham, E. Brevdo, Z. Chen, C. Citro, G. S. Corrado, A. Davis, J. Dean, M. Devin, et al., Tensorflow: Large-scalemachine learning on heterogeneous distributed systems, arXiv preprint arXiv:1603.04467 (2016).[47] N. S. Keskar, D. Mudigere, J. Nocedal, M. Smelyanskiy, P. T. P. Tang, On large-batch training for deep learning: Generalization gap andsharp minima, arXiv preprint arXiv:1609.04836 (2016).[48] A. G. Holzapfel, Nonlinear solid mechanics II, John Wiley & Sons, Inc., 2000.[49] S. Lin, L. Gu, Influence of crosslink density and sti ff ness on mechanical properties of type i collagen gel, Materials 8 (2015) 551–560.[50] V. Mullachery, A. Khera, A. Husain, Bayesian neural networks, arXiv preprint arXiv:1801.07710 (2018).[51] R. M. Neal, Bayesian learning for neural networks, volume 118, Springer Science & Business Media, 2012.[52] A. J. Licup, S. M¨unster, A. Sharma, M. Sheinman, L. M. Jawerth, B. Fabry, D. A. Weitz, F. C. MacKintosh, Stress controls the mechanics ofcollagen networks, Proceedings of the National Academy of Sciences 112 (2015) 9573–9578.[53] L. Zhang, Microstructural modeling of cross-linked fiber network embedded in continuous matrix, Rensselaer Polytechnic Institute, 2013.[54] M. Aghvami, K. L. Billiar, E. A. Sander, Fiber network models predict enhanced cell mechanosensing on fibrous gels, Journal of biome-chanical engineering 138 (2016).[55] R. Raghupathy, V. H. Barocas, A closed-form structural model of planar fibrous tissue mechanics, Journal of biomechanics 42 (2009)1424–1428.[56] M. Liu, L. Liang, W. Sun, A generic physics-informed neural network-based constitutive model for soft biological tissues, Computer Methodsin Applied Mechanics and Engineering 372 (2020) 113402.ness on mechanical properties of type i collagen gel, Materials 8 (2015) 551–560.[50] V. Mullachery, A. Khera, A. Husain, Bayesian neural networks, arXiv preprint arXiv:1801.07710 (2018).[51] R. M. Neal, Bayesian learning for neural networks, volume 118, Springer Science & Business Media, 2012.[52] A. J. Licup, S. M¨unster, A. Sharma, M. Sheinman, L. M. Jawerth, B. Fabry, D. A. Weitz, F. C. MacKintosh, Stress controls the mechanics ofcollagen networks, Proceedings of the National Academy of Sciences 112 (2015) 9573–9578.[53] L. Zhang, Microstructural modeling of cross-linked fiber network embedded in continuous matrix, Rensselaer Polytechnic Institute, 2013.[54] M. Aghvami, K. L. Billiar, E. A. Sander, Fiber network models predict enhanced cell mechanosensing on fibrous gels, Journal of biome-chanical engineering 138 (2016).[55] R. Raghupathy, V. H. Barocas, A closed-form structural model of planar fibrous tissue mechanics, Journal of biomechanics 42 (2009)1424–1428.[56] M. Liu, L. Liang, W. Sun, A generic physics-informed neural network-based constitutive model for soft biological tissues, Computer Methodsin Applied Mechanics and Engineering 372 (2020) 113402.