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.
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 of actin monomer is determined by the interaction potential , an overall repulsive potential , and a background stochastic force from solvent molecules , where is the vector of all monomer positions, is the gradient with respect to position vector of monomer , and the stochastic forcing is assumed to be Gaussian with indices , running over spatial components, the drag coefficient, the Boltzmann constant, and the temperature. The interaction potential 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 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 a local linearization around the state can be written as
We assume that describes deviation from the reference state such that the time average over simulation time . 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, , , 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 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 of the simulation domain edge length , i.e., . 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, , , with in order to accurately describe a portion of the cytoskeleton of volume μm3, sufficient to determine protrusion force [15]. For a typical human cell volume of 100 μm3, the number of microscopic state variables would be . We seek a reduced description of the mesh within a 1 μm3 control volume through a linear combination , with an orthonormal set of basis vectors. The reduced linear system obtain by projection of (4) onto is with , , and the reduced drag, stiffness matrices, and reduced force, respectively.
The drag and stiffness matrices (both at full resolution , , and reduced resolution , ) are not known explicitly and cannot be readily evaluated through Taylor series expansions such as due to the complex forms of the interaction potentials and 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 that models the solvent stochastic force and a periodic forcing along a single direction as would be produced by an extensional rheometer, (here, denotes a scalar coordinate). The positional data at time steps , , is given by the matrix
Both the drag and stiffness matrices must be symmetric positive definite due to physical considerations ( gives the elastic deformation energy, 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 in which and are orthogonal eigenvector matrices, and the diagonal eigenvalue matrices and have positive components.
The descriptive capability of the reduced model depends crucially on the choice of a basis set that efficiently captures the possible configurations of the actin mesh in as few modes as possible, . To this end, the basis set is chosen from the dominant modes of the singular value decomposition (SVD), . Note that the correlation matrix can be expressed as , and the dominant singular modes 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 , which are simply the coefficients of the linear combination . In some particular cases, it is possible to use straightforward averaging over geometric positions to attach significance to the reduced parameter vector , 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 as , to relate the displacements at the control nodes within dimensions linearly through the element stiffness matrix to the nodal forces . The corresponding formulation for a viscoelastic element is
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, hence , , and the displacements within the element are given as where are the displacements at node of the element, and are form functions with , , coordinates from the element center, . The interpolation (10) can be expressed as with
We now consider the procedure to link the finite element formulation to the data available from the BD simulation. At time step , construct the data set of displacements of all monomers within the mesh . A least squares fit of data to the finite element approximation (12) at time , , could furnish the nodal values , and a subsequent fit to the observed nodal values would provide the element drag and stiffness matrices we seek. However, such a procedure would include all the observed motion, including the stochastic thermal effects. Rather than fitting to , the nodal values are obtained by fitting onto the data projected onto the dominant modes contained in the data set with the projection of onto , . 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 sets , , the nodal displacements are known, and the nodal velocities can be approximated by, e.g., finite difference approximations
As in the case of the full drag, stiffness matrices, similar physical constraints are imposed onto the finite element matrices with , orthogonal matrices and components of diagonal matrices , constrained to be positive.
The columns of and 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 and . The components of and 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.
Consider instances of actin mesh configurations obtained from distinct BD simulations (1), from which the basis sets are obtained through SVDs. The ensemble average basis set captures both variability from thermal forcing and that due to mesh configuration through different initial states . The column of is denoted by can be interpreted as the coefficient set needed to reconstruct an approximation of the microscopic mesh displacement through the finite element interpolation (12) as
A link can now be made to the available data by choosing and computing the solutions , for of the least squares problems
The vectors are not necessarily orthogonal as desired for the columns of , so a -decomposition is computed to obtain . 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 to obtain the velocity modes . 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 modes and are determined, the remaining task to completely define the finite element model is to compute the diagonal matrices and by a least-squares fit over the time history determined by fitting to data sets , i.e., solving the problem
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 of (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 monomers at evenly spaced time steps during the period . 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.
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 . During the first quarter period 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 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 exhibits a different average slope indicating strain hardening.
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) 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 and .
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 , 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.
Data Availability
Simulation data and codes are available and may be directly requested from the corresponding author.
Conflicts of Interest
The authors declare that they have no conflicts of interest.
Authors’ Contributions
BA carried out numerical simulations of the F-actin network, and SM developed the model reduction procedure to define a finite element model.
Acknowledgments
This work was supported by the NSF (CMMI-1068918, DMS-1361375).