Complexity
Volume 2018 (2018), Article ID 6878265, 7 pages
https://doi.org/10.1155/2018/6878265

Research Article

Data-Driven Finite Element Models of Passive Filamentary Networks

Brian Adam and Sorin Mitran

Department of Mathematics, University of North Carolina at Chapel Hill, Chapel Hill, NC 27599-3250, USA

Correspondence should be addressed to Sorin Mitran; mitran@unc.edu

Received 6 July 2018; Accepted 20 September 2018; Published 1 November 2018

Academic Editor: Francisco Chinesta

Copyright © 2018 Brian Adam and Sorin Mitran. This is an open access article distributed under the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

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  μm³, sufficient to determine protrusion force [15]. For a typical human cell volume of 100 μm³, the number of microscopic state variables would be . We seek a reduced description of the mesh within a 1 μm³ 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.

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 . 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.

Figure 2: Comparison of strain-stress curve of microscopic model (plot markers, distinct for each quarter period of the forcing period ) 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) 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).

References

L. Blanchoin, R. Boujemaa-Paterski, C. Sykes, and J. Plastino, “Actin dynamics, architecture, and mechanics in cell motility,” Physiological Reviews, vol. 94, no. 1, pp. 235–263, 2014.

A. S. Abhilash, B. M. Baker, B. Trappmann, C. S. Chen, and V. B. Shenoy, “Remodeling of fibrous extracellular matrices by contractile cells: predictions from discrete fiber network simulations,” Biophysical Journal, vol. 107, no. 8, pp. 1829–1840, 2014.

J. Stricker, T. Falzone, and M. L. Gardel, “Mechanics of the F-actin cytoskeleton,” Journal of Biomechanics, vol. 43, no. 1, pp. 9–14, 2010.

M. L. Gardel, J. H. Shin, F. MacKintosh, L. Mahadevan, P. Matsudaira, and D. A. Weitz, “Elastic behavior of cross-linked and bundled actin networks,” Science, vol. 304, no. 5675, pp. 1301–1305, 2004.

C. Storm, J. J. Pastore, F. C. MacKintosh, T. C. Lubensky, and P. A. Janmey, “Nonlinear elasticity in biological gels,” Nature, vol. 435, no. 7039, pp. 191–194, 2005.

C. . P. Broedersz and F. . C. MacKintosh, “Modeling semiflexible polymer networks,” Reviews of Modern Physics, vol. 86, no. 3, pp. 995–1036, 2014.

F. Meng and E. M. Terentjev, “Theory of semiflexible filaments and networks,” Polymers, vol. 9, no. 12, 2017.

M. Klann and H. Koeppl, “Spatial simulations in systems biology: from molecules to cells,” International Journal of Molecular Sciences, vol. 13, no. 6, pp. 7798–7827, 2012.

T. Kim, W. Hwang, H. Lee, and R. D. Kamm, “Computational analysis of viscoelastic properties of crosslinked actin networks,” PLoS Computational Biology, vol. 5, no. 7, article e1000439, 2009.

G. Zagar, P. R. Onck, and E. van der Giessen, “Two fundamental mechanisms govern the stiffening of cross-linked networks,” Biophysical Journal, vol. 108, no. 6, pp. 1470–1479, 2015.

O. Lieleg, K. M. Schmoller, M. M. A. E. Claessens, and A. R. Bausch, “Cytoskeletal polymer networks: viscoelastic properties are determined by the microscopic interaction potential of cross-links,” Biophysical Journal, vol. 96, no. 11, pp. 4725–4732, 2009.

M. H. Jensen, E. J. Morris, and D. A. Weitz, “Mechanics and dynamics of reconstituted cytoskeletal systems,” Biochimica et Biophysica Acta (BBA) - Molecular Cell Research, vol. 1853, no. 11, pp. 3038–3042, 2015.

C. J. Cyron, K. W. Müller, A. R. Bausch, and W. A. Wall, “Micromechanical simulations of biopolymer networks with finite elements,” Journal of Computational Physics, vol. 244, pp. 236–251, 2013.

J. Zhang, X. Zhao, Z. Suo, and H. Jiang, “A finite element method for transient analysis of concurrent large deformation and mass transport in gels,” Journal of Applied Physics, vol. 105, no. 9, 2009.

Y. Inoue, T. Deji, and T. Adachi, “Brownian dynamics simulation study on force–velocity relation in actin-based membrane protrusion,” Computational Particle Mechanics, vol. 2, no. 4, pp. 329–337, 2015.

G. G. Borisy and T. M. Svitkina, “Actin machinery: pushing the envelope,” Current Opinion in Cell Biology, vol. 12, no. 1, pp. 104–112, 2000.

K. Keren, Z. Pincus, G. M. Allen et al., “Mechanism of shape determination in motile cells,” Nature, vol. 453, no. 7194, pp. 475–480, 2008.

A. Mogilner and G. Oster, “Cell motility driven by actin polymerization,” Biophysical Journal, vol. 71, no. 6, pp. 3030–3045, 1996.

A. Mogilner and G. Oster, “Force generation by actin polymerization ii: the elastic ratchet and tethered filaments,” Biophysical Journal, vol. 84, no. 3, pp. 1591–1605, 2003.

A. Mogilner and B. Rubinstein, “The physics of filopodial protrusion,” Biophysical Journal, vol. 89, no. 2, pp. 782–795, 2005.

J. Gernandt, G. Frenning, W. Richtering, and P. Hansson, “A model describing the internal structure of core/shell hydrogels,” Soft Matter, vol. 7, no. 21, pp. 10327–10338, 2011.

J. M. Holzwarth and P. X. Ma, “Biomimetic nanofibrous scaffolds for bone tissue engineering,” Biomaterials, vol. 32, no. 36, pp. 9622–9629, 2011.

P. X. Ma and R. Y. Zhang, “Synthetic nano-scale fibrous extracellular matrix,” Journal of Biomedical Materials Research, vol. 46, no. 1, pp. 60–72, 1999.

T. A. Ulrich, A. Jain, K. Tanner, J. L. MacKay, and S. Kumar, “Probing cellular mechanobiology in three-dimensional culture with collagen–agarose matrices,” Biomaterials, vol. 31, no. 7, pp. 1875–1884, 2010.

Y. Feng and S. Mitran, “Data-driven reduced-order model of microtubule mechanics,” Cytoskeleton, vol. 75, no. 2, pp. 45–60, 2018.