# Data-Driven Finite Element Models of Passive Filamentary Networks

Authors: Brian Adam; Sorin Mitran
Journal: Complexity
 (2018)

---

## Abstract

A data-driven procedure is introduced to construct finite element models of heterogeneous systems for which an accurate microscopic description is available. A filter to define coarsened finite element nodal values is defined from the principal modes obtained by singular value decomposition of the microscopic data. The resulting finite element nodal values are subsequently used to reconstruct local linearization of the system behavior, defining drag and stiffness matrices for an overdamped system. The procedure is exemplified for an actin mesh described by Brownian dynamics and eight-node cuboid finite elements but is generally applicable with respect to both the microscopic model and the type of finite element approximation. In contrast to standard finite element formulations derived from hypotheses on assumed deformation behavior, the data-driven procedure introduced here is completely determined by the observed behavior be it obtained from simulations or experiment.


---

## Body

## 1. Introduction

Materials with heterogeneous microscopic structure arise in a wide variety of applications and are frequently encountered in nature. For example, the mechanical behavior of a biological cell is largely determined by a network of actin filaments that actively polymerize and depolymerize exerting force on the cellular membrane and leading to cellular motility [1]. Many cellular processes depend on remodeling the actin network configuration leading to diverse mechanical responses [2]. The filamentary actin (F-actin) network exhibits viscoelastic behavior that becomes nonlinear at large deformations [3]. Observed elastic moduli are strongly influenced by dynamically varying cross-linking between network filaments ranging from 0.03 Pa to 300 Pa depending on cross-linker concentration [4]. While the qualitative physics of such networks [4, 5] can identify regions within the phase space of concentrations of various components, detailed and quantitatively correct models of cellular mechanics will probably require numerical simulation. Various approaches to modeling of networks arising in biology are reviewed in [6, 7]. The main difficulty encountered in constructing accurate models is that while the basic behavior of individual components within a cell are fairly well understood, the overall complexity renders such almost first-principle approaches exceedingly expensive. For example, Brownian dynamic methods have been attempted for entire cell simulation [8] or extraction of viscoelastic properties [9], at considerable computational effort. Such studies have identified network stiffening mechanisms [10], the crucial role played by cross-linkers [11] in establishing viscoelastic properties of the network, that compare favorably with the mechanics of networks reconstituted from observation [12].Of particular interest here are finite element models for the actin mesh. One line of research [13] introduces various types of behavior for individual components of the network such as filament elements or cross-linkers. At the other end of the modeling resolution scale are finite element gel continuum [14] models based upon large-deformation mechanics and an assumed constitutive relationship. Related to these approaches, the procedure introduced in this work seeks to use detailed microscopic simulation to obtain a finite element model, since this is the resolution level at which the basic biophysical processes are well defined, and obtain a homogenized model. In contrast to aforementioned approaches, neither is a constitutive relationship at the continuum scale assumed to be known, nor is the goal to construct a finite element model of each constituent element of the network. The approach consists in leveraging microscopic simulation (Brownian dynamics in this work) to extract an effective finite element model, defining the constitutive relationships underlying network response to forcing from available microscopic data. As described in more detail below, the approach combines a succession of data reduction procedures based upon projection in Euclidean spaces. Data from the microscopic model are filtered by the principal modes observed during the motion. The filtered data are used to define nodal deformation values. The principal modes are also reduced to smaller dimension to define effective drag and stiffness matrices describing the linearized viscoelastic behavior of the network. Final model coefficients are determined by fitting to observed behavior over an ensemble of network instances to account for variability in network configuration.


## 2. Methods

The Brownian dynamics (BD) model of [15] is adopted to describe the actin network. In this approach, the velocity vit=xit of actin monomer i is determined by the interaction potential Uix, an overall repulsive potential Uimx, and a background stochastic force from solvent molecules ξ,
(1)γvi=−∇iUix−∇iUimx+ξ,where x=x1…xNT∈ℝm is the vector of all N=m/3 monomer positions, ∇i=∂/∂xi is the gradient with respect to position vector of monomer i, and the stochastic forcing is assumed to be Gaussian
(2)ξtξt′=2γkBTδαβδt−t′,ξt=0,with indices α, β running over spatial components, γ the drag coefficient, kB the Boltzmann constant, and T the temperature. The interaction potential Ui contains terms describing interactions of monomers within filaments (length and angle dependent) and a volume exclusion term of Lennard Jones form. The overall repulsive potential Uim is also of Lennard Jones form. No length and angle interaction terms act upon free actin monomers. For full details consult [15], including the model coefficients used in computations here.The BD (1) form a nonlinear system of stochastic differential equations that accurately describes the processes of actin polymerization that leads to generation of mechanical force. For studies over a small part of a cell, such as in the immediate vicinity of a small portion of the cellular membrane in order to characterize the Brownian ratchet mechanism of cellular protrusion [16–20], direct numerical solution of the BD system (1) is feasible. However, for larger-scale studies of the actin mesh network over an entire cell, the computational cost becomes prohibitive, and a reduced model capturing essential features is required. In particular, it is desirable to obtain a model of the constitutive relationship for the homogenized model informed by the accurate detailed BD model. Though the discussion here is presented for the particular case of the actin mesh network forming the cytoskeleton of a cell, the above situation is generic and also appears in consideration of various gels [5, 21] or the collagen extracellular matrix [22–24] among other applications.From the nonlinear system describing time evolution of the mesh network(3)ẋ=ht,x,ξ,a local linearization around the state x0 can be written as
(4)Dẋ+Kx=f.We assume thatxt describes deviation from the reference state x0 such that the time average over simulation time Txt=T−1∫0Txtdt=x0. Given the nonlinearity and possible history dependence of the network configuration, it can be expected that the drag and stiffness matrices depend on the chosen reference state, Dx0, Kx0, and the analysis presented below would need to be repeated for markedly different reference states. In the example presented here, an ensemble of reference states x0 is generated by simulating growth of a mesh from random placement of actin monomers within the computational domain until the average filament length ℓ reaches some fraction f of the simulation domain edge length L, i.e., f=ℓ/L=0.1. The resulting ensemble exhibits randomness of filament orientation, initial stress produced by filament growth, and cross-linkage formation. The ensemble is representative of newly formed actin meshes, with considerable reconfiguration possible by repeated forcing (e.g., as produced by motile behavior).The microscopic mesh state is characterized by positions, velocities of all monomers,xt, ẋt∈ℝm, with m~O105 in order to accurately describe a portion of the cytoskeleton of volume ~1 μm3, sufficient to determine protrusion force [15]. For a typical human cell volume of 100 μm3, the number of microscopic state variables would be O107. We seek a reduced description of the mesh within a 1 μm3 control volume through a linear combination x=Bpy, with Bp∈ℝm×p an orthonormal set of p basis vectors. The reduced linear system obtain by projection of (4) onto spanB is
(5)Dpẏ+Kpy=g,with Dp=BpTDBp, Kp=BpTKBp, and g=BpTf the reduced drag, stiffness matrices, and reduced force, respectively.The drag and stiffness matrices (both at full resolutionD, K, and reduced resolution Dp, Kp) are not known explicitly and cannot be readily evaluated through Taylor series expansions such as K=∂hx0/∂x due to the complex forms of the interaction potentials Uix and Uimx and the action of the stochastic forcing term ξ. Rather than seeking such an analytical derivation, a data-driven approach is adopted based upon results obtained from numerical simulation with a forcing term
(6)f=axsinωt+ξ,that models the solvent stochastic force ξ and a periodic forcing along a single direction as would be produced by an extensional rheometer, axsinωt (here, x denotes a scalar coordinate). The positional data xti at n time steps ti=iδt, δt=T/n, is given by the matrix
(7)X=x1…xn∈ℝm×n.Both the drag and stiffness matrices must be symmetric positive definite due to physical considerations (xTKx≥0 gives the elastic deformation energy, ẋTDẋ≥0 gives the dissipation energy due to drag of solvent on the monomers), and similar properties are imposed for the reduced matrices through the eigen-decompositions
(8)Dp=VpΛpvVpT,Kp=UpΛpxUpT,in which Vp and Up are orthogonal eigenvector matrices, and the diagonal eigenvalue matrices Λpv and Λpx have positive components.The descriptive capability of the reduced model depends crucially on the choice of a basis setBp that efficiently captures the possible configurations of the actin mesh in as few modes as possible, p≪m. To this end, the basis set is chosen from the dominant modes of the singular value decomposition (SVD), X=BΣST. Note that the correlation matrix C=XXT can be expressed as C=BΣΣTBT, and the p dominant singular modes Bp=b1…bp∈ℝm×p correspond to those eigenmodes of the correlation matrix most readily observed during the simulation.In the reduced model (5), no particular significance is attributed to the components of y, which are simply the coefficients of the linear combination x=Bpy. In some particular cases, it is possible to use straightforward averaging over geometric positions to attach significance to the reduced parameter vector y, e.g., for one-dimensional elements as in [25]. Such an approach is difficult to justify in higher dimensions though.Recall that a standard finite element for an elastic body is expressed asfe=Keue, to relate the displacements ue∈ℝdq at the q control nodes within d dimensions linearly through the element stiffness matrix Ke to the nodal forces fe. The corresponding formulation for a viscoelastic element is
(9)Deu̇e+Keue=fe.The positions of the control nodes can be chosen arbitrarily but are typically taken to be on the boundary of the finite element in order to invoke continuity of deformation between adjacent elements. In the following, we consider the control nodes at the corners of a three-dimensional right rectangular prism, to obtain the commonly used brick element, henced=3, q=8, and the displacements uex∈ℝd within the element are given as
(10)uet,x=∑k=1qNkxuket,where uket∈ℝd are the displacements at node k of the element, and Nkx are form functions
(11)Nkx∈181±x1±y1±z,k=1,…,8,with x, y, z coordinates from the element center, x=x,y,z. The interpolation (10) can be expressed as
(12)uet,x=NTxuet,with
(13)NTx=N1x…Nqx,uet=u1et…uqet.We now consider the procedure to link the finite element formulation to the data available from the BD simulation. At time stepti=iδt, construct the data set
(14)Di=xj0,xjti−xj0,j=1,…,Nof displacements of all monomers within the mesh ui,j≡uti,xj0=xjti−xj0. A least squares fit of data Di to the finite element approximation (12) at time ti, ueti,x, could furnish the nodal values ueti, and a subsequent fit to the observed nodal values would provide the element drag and stiffness matrices De,Ke we seek. However, such a procedure would include all the observed motion, including the stochastic thermal effects. Rather than fitting to Di, the nodal values are obtained by fitting onto the data projected onto the dominant p modes contained in the data set
(15)D¯i=x¯j0,x¯jti−x¯j0,j=1,…,N,with x¯jt the projection of xjt onto spanBp, x¯jt=BpBpTxjt. In essence, the dominant modes obtained from the SVD are used as a filter of the observed microscopic motion.After carrying out the least squares fit to data setsD¯i, i=1,…,n, the nodal displacements ueti are known, and the nodal velocities can be approximated by, e.g., finite difference approximations
(16)veti≡u̇eti≅ueti+1−ueti−12δt.As in the case of the full drag, stiffness matrices, similar physical constraints are imposed onto the finite element matrices(17)De=VeΛveVeT,Ke=UeΛxeUeT,with Ve, Ue∈ℝdq×dq orthogonal matrices and components of diagonal matrices Λve, Λxe constrained to be positive.The columns ofVe and Ue describe the velocity and positional modes of the finite element. In standard finite element analysis of an elastic body, such modes would result from hypotheses on the deformation behavior. Typically, stretch, shear, twist, and bending of the finite element would be of interest, assumed to be independent of another (i.e., orthogonal modes), and used to construct Ve and Ue. The components of Λve and Λxe would result from, say, a Lagrangian formulation of the equations of motion. Such an approach is reasonable for a homogeneous material, but the deformation modes of an actin mesh are likely to be strongly dependent on mesh configuration and can be expected to always couple the intuitively derived stretch, shear, twist, and bending modes. Rather than impose assumed deformation modes, we seek to construct these from an ensemble of actin mesh configurations.Considers instances of actin mesh configurations X1,…,Xs∈ℝm×n obtained from s distinct BD simulations (1), from which the basis sets Bp1,…,Bps∈ℝm×p are obtained through SVDs. The ensemble average basis set
(18)B~p=1s∑l=1sBplcaptures both variability from thermal forcing and that due to mesh configuration through different initial states x10,…,xs0. The kth column of Ue∈ℝdq×dq is denoted by uke∈ℝdq can be interpreted as the coefficient set needed to reconstruct an approximation of the microscopic mesh displacement ukx through the finite element interpolation (12) as
(19)ukx=x−x0≅NTx1uke⋮NTxNuke≡wx,uke.A link can now be made to the available data by choosingp=dq and computing the solutions zke, for k=1,…,p of the p least squares problems
(20)minukeb~k−x0−wx,zke.The vectorszke are not necessarily orthogonal as desired for the columns of Ue, so a QR-decomposition UeR=Ze≡z1e…zpe is computed to obtain Ue. Note that the above procedure homogenizes both variability in mesh configuration and that due to stochastic effects to obtain data-driven deformation modes. An analogous procedure is applied to the velocity data Ẋ to obtain the velocity modes Ve. For a purely elastic system, the velocity modes would be identical to the deformation modes, but for viscoelastic systems, the modes can differ due to phase differences and history dependence (note that even though each monomer is affected by the same type of thermal stochastic forcing ξ, the interaction potentials can induce correlation between the stochastic modes of the overall system).Once the reduced modesUe and Ve are determined, the remaining task to completely define the finite element model is to compute the diagonal matrices Λxe and Λve by a least-squares fit over the time history ueti,u̇eti determined by fitting to data sets D¯i, i.e., solving the problem
(21)minΛxe,Λve∑i=1nVeΛveVeTu̇eti+UeΛxeUeTueti−feti.


## 3. Results

An initial test of the above data-driven finite element construction procedure is carried out for data from seven network instances forced at nondimensional periods ofT=1.57,3.14 (see [15] for details on reference units used to obtain nondimensional quantities). The goal here is to carry out basic verification of the model reduction procedure rather than obtain a fully realistic description of the biomechanics of F-actin networks. The networks contain N=1706 monomers at n=2000 evenly spaced time steps during the period T. A typical sequence of monomer displacement vectors is shown in Figure 1 and is included here mainly to highlight the difficulty in ascertaining dominant deformation types from such data. Though the forcing (6) is applied along a single coordinate direction, the random nature of filamentary orientation and cross-linkage leads to a complex response of the structure.Figure 1
Sequence of actin monomer positions at six equidistant times during a period of the forcing. Notice that some free actin monomers exhibit large displacements, and the overall response to forcing along a single coordinate direction is difficult to ascertain from this data.A stress-strain curve along the forcing direction may be extracted from the microscopic data and compared to that obtained by integrating the finite element model (9). The results are presented in Figure 2. The microscopic model shows considerable variation and hysteresis over the period T. During the first quarter period 0,T/4 corresponding to rising extensional stress applied to the system, the filament network on average shows rising strain but interspersed with strain release events that probably correspond to network reconfiguration. During the next half period T/4,3T/4 of decreasing extensional stress, strain release occurs but with clear plastic behavior due to the viscous terms in the system, and the final quarter period 3T/4,T/2 exhibits a different average slope indicating strain hardening.Figure 2
Comparison of strain-stress curve of microscopic model (plot markers, distinct for each quarter period of the forcing periodT) and the strain stress curve resulting from the finite element model (9).No claim is made within these preliminary results of accurate capturing of true biological F-actin meshes, but the reduced finite element model constructed by seeking periodic solutions of (9) uet=a0+∑k=18akcos2πkt/T+bksin2πkt/T captures the main features of the microscopic simulation to within a relative error of 3.3% for the stress-strain curve. Note that the microscopic data exhibit initial transient behavior, not the long-term periodic behavior that results from repeated application of the oscillatory extensional force (6), and this is observed in the behavior at t≅0 and t≅T.


## 4. Discussion

A general procedure has been introduced to define a finite element model from data obtained by Brownian dynamics simulation of an F-actin network. Though the presentation used data from numerical simulation and exemplified the procedure using an eight-node cuboid finite element, the approach can readily be applied to obtain different types of finite element approximation or use data from experiments. Of particular interest is the extraction of essential behavior from the available data without any assumption on deformation modes or constitutive relationships within the material. For the case of the F-actin network considered here, the derived finite element approximation captures the stress-strain relationship from the microscopic model with a relative error ofεr≤0.04, at negligible computational cost by comparison to Brownian dynamics simulations. Multiple application of the procedure can be used to model nonlinear behavior or parametrized to model active systems concurrently with biochemical reaction models.At this point, the reduced finite element model is deterministic, capturing the average behavior of the filamentary network. The procedure however forms the basis of establishing a stochastic finite element approximation through recognition of the close relationship between the singular value decomposition (SVD) and the Karhunen-Loeve expansion. In this work, only the singular modes from the SVD were used in constructing a deterministic finite element approximation. Such modes correspond to the eigenmodes in a Karhunen-Loeve expansion, while the singular values give the mode variances. Use of the variances would allow formulation of a stochastic finite element, the subject of upcoming work.



---

*Source: 6878265-2018-11-01.xml*