Variational approximation methods for long-range force transmission in biopolymer gels

Haiqin Wang(王海欽) and Xinpeng Xu(徐新鵬)

1Physics Program,Guangdong Technion-Israel Institute of Technology,241 Daxue Road,Shantou 515063,China

2Technion-Israel Institute of Technology,Haifa 3200003,Israel

Keywords: biopolymer gels,cell-cell communications,force transmission,variational methods

1. Introduction

Cells in animal tissues are surrounded by extracellular matrix(ECM).[1]Physical,topological,and biochemical compositions of ECMs are not only complex, tissue-specific, but also markedly heterogeneous.[1]However,in mostin vitroexperiments,only one or two types of constituting proteins(such as collagen, fibrin, and elastin) are extracted from ECMs to form cross-linked biopolymer gels.[1–3]A striking feature of these reconstituted biopolymer gels is their asymmetric elastic response to extension (or shear) and compression.[4–6]They stiffen under small shear stresses(or small strains about 5%–10%)with their elastic moduli increasing as the 3/2-power law of the applied stress. In contrast,biopolymer gels soften when compressed upon a small amount where they almost lose resistance to shear stress completely.[7,8]Such nonlinear elasticity of fibrous biopolymer gels can be attributed to the microstructural nonlinearities of the constituent semiflexible filaments which stiffen (due to inextensibility) under extension and soften(due to buckling)under compression.[4,9–12]

When cells are embedded in such fibrous biopolymer gels, they pull on the gel and produce non-equilibrium forces by myosin motors that consume adenosine triphosphate(ATP).Recent experiments have shown that the displacements and structure changes induced by cells adhered to nonlinear biopolymer gels can reach a distance of tens of cell diameters away,[13–17]in contrast to the distance of several cell diameters reached by displacements induced by cells on linear synthetic gels. Such phenomena support long-range cellcell mechanical communication, a process that can mechanically couple distant cells and coordinate processes such as capillary sprouting[17]and synchronous beating.[18]The longrange transmission of cellular forces is usually attributed to the unique nonlinear mechanics and fibrous nature of the biopolymer gels.[2,4,12,14,19–23]Particularly, the decay of cell-induced displacement has been measured quantitatively in anin vitroexperiment for a fibroblast spheroid embedded and contracting in a three-dimensional (3D) fibrin gel[16](see Fig. 1(a)).The cell-induced displacement is found(as shown in Fig.1(b))to decay over distance by a power-law scaling asu～r-0.52,which is much more slowly than that induced by cells in linear elastic medium(scaling asu～r-2).In this work,we show that these experimental observations and measurements can be explained by a direct variational analysis based on the threechain continuum model of biopolymer gels[12,19]and the variational principle of minimum free energy(MFEVP).

Fig. 1. Slow power-law decay of the displacements induced by contracting cells embedded in a 3D biopolymer gel measured in experiments. (a)The displacement vector fields (color quivers) induced by an individual fibroblast spheroid in a 3D fibrin gel. (b) The decay of the displacement magnitude over the distance r away from the cell center is averaged for multiple time points. Each curve is for a different cell. An effective near-field power-law u ～r-0.52 is obtained, which decays much slower than that in linear isotropic medium with u ～r-2. Reproduced from Notbohm[16] with permission from Royal Society.

The rest of this paper is organized as follows: We introduce general variational methods in Section 2 for static problems in elastic materials.In Section 3,we apply the variational methods to study the transmission of internal forces applied by a spherically contracting cell in a 3D biopolymer gel. In Section 4,we summarize our major results and give some remarks.

2. Variational methods for elastic materials

Variational principles have been widely used in the continuum modeling of soft and biological matter such as the variational principle of minimum free energy (MFEVP) for static problems[24]and Onsager’s variational principle(OVP)for dynamic problems.[25–30]In this work, we apply MFEVP to study the transmission of internal cellular forces in elastic biopolymer gels.

2.1. Variational formulation for the continuum modeling of elastic biopolymer gels

We consider a general elastic material,in which the total free energy functional can be written as

The explicit form of the deformation energy density?edepends on the structure, interactions, and properties of the constituents of the elastic material. Particularly for biopolymer gels composed of crosslinked stiff biopolymer filaments,they show highly nonlinear elastic properties as schematically shown in Figs. 2(a1)–2(a2):[4–8]They stiffen in the stretched direction and soften in the compressed direction due to the inextensibility and microbuckling of individual filaments, respectively. Several energy forms have been proposed to to model these unique nonlinear elasticity from some continuum models analogous to rubber elasticity,[4,11,12,19,32,33]e.g.,1-chain sphere models, 3- and 8-chain cubic lattice models,and 4-chain tetrahedra model. In these models, the macroscopic elastic deformation energy of the gel is obtained by adding the free energy of individual blocks (or material elements)that are deformed affinely. Here we employ the simple three-chain model that are easier to handle analytically and are found to best fit the experimental data for the elastic deformation of various biopolymer gels at small deformation.[12,19]In a typical three-chain model for cross-linked biopolymer gels as shown in Figs.2(b1)–2(b2),a primitive cubic block is constructed with lattice points representing the cross-linking sites,and the three-chain segments between sites at the block edges are aligned along the three principal directions of deformation (denoted by ?xiwithi= 1,2,3). The elasticity of the chain segments represents the emergent response of the constituent polymers to applied forces. In the three-chain model for biopolymer gels,it is usually simply assumed that the chain segments have the same nonlinear elasticity as single semiflexible biopolymers. In this case,the deformation energy densityFetakes the following form:[12,19,33]

in terms of the three principal strain componentsεi(i=1,2,3)and for the special case ofε1＞0 andε2,ε3＜0. Here ?K ≡K0-2μ0/3 is the modified bulk modulus[12,19]withμ0andK0being the linear shear and bulk moduli, respectively,which are related to Young’s modulus and Poisson’s ratio byμ0=E0/2(1+ν0),andK0=E0/3(1-2ν0). The dimensionless parameter ?ρ2,3is given by ?ρ2,3=1-(1-ρ(ε2,3))(1+εb/ε2,3)2withρ(ε)=ρ0+(1-ρ0)Θ(ε+εb),Θbeing the Heaviside step function,and 0≤ρ0?1.Note that the stretchstiffening and the compressive-softening nonlinear elasticity of biopolymer gels are captured in the energy densityFeby the terms,1/(1-ε1/εs)(in the stretched principal direction withε1＞0), and ?ρ2,3with a small softening ratioρ0(in the compressed principal directions withε2,ε3＜0),respectively. The nonlinear elasticity becomes significant when the strain magnitude is comparable to critical strainsεs(stretch)andεb(compression),respectively.In experiments,εsandεbare measured(in fibrin gels[4–8]) to be around 10% and 1%, respectively.Furthermore, fromσei=?Fe/?εiwe obtain the three principal stress components:

Fig.2. Schematic illustrations of the continuum elastic model of biopolymer gels in both rectangular geometry(for homogeneous deformation)and spherical geometry (for inhomogeneous deformation) at reference states [(a1), (b1), (c1)] and deformed states [(a2), (b2), (c2)]. [(a1),(a2)]Deformation of individual biopolymer filaments in the stretched direction(stiffening due to inextensibility along the ?x1 direction)and the compressed direction(softening due to buckling along the ?x2 direction),respectively.[(b1),(b2)]The physical picture of the three-chain model of a 3D biopolymer gels along three principal directions. The gel stiffens along the stretched ?x1 direction and softens along the compressed ?x2,3 directions. [(c1),(c2)]A spherically contracting cell of radius Rc is embedded in an infinite 3D biopolymer gel in spherical geometry. The cell generates a displacement-uc at its boundary with uc ＞0 for contractile cells. We compare the theoretical results with the experimental data in the range of Rc ≤r ＜Rm.

2.2. Variational methods of approximation: Ritz method

Now, we have shown that variational principles such as MFEVP provide an equivalent(and more convenient)applications of vectorial governing(force balance or Euler–Lagrange)equations. However, variational principles should not be regarded only as a mathematical substitute or reformulation of force balance equations. They also provide some powerful variational methods of finding approximate solutions to these equations, e.g., the Ritz method and the least-squares method.[34]In these variational methods, some trial solutions to the problem are assumed where the state variable functions are taken as combinations of some simple functions with a small number of adjustable parameters. Then the total free energy that is a functional of state variables can be integrated over space and reduces to a function of these adjustable parameters. Correspondingly,the functional minimization of the free energy with respect to state variables reduces to a function minimization with respect to a small number of parameters, which gives approximate solutions of the static problem and determines the equilibrium properties of the system. Such methods simplify the static boundary value problem significantly. They bypass the derivation of the governing equilibrium equations and go directly from the variational statement to an approximate solution of the static problem.These simplified solution methods are, therefore, called direct variational methods or variational methods of approximation.[34]

Note that the trial solutions can be either completely empirical arising from experiences gained in systematic numerical analysis or experimental measurements,[28,35]or assumed to be linear combinations of a finite set of basis functions such as algebraic and trigonometric polynomials.[34]The latter choice of trial solutions is usually known as the Ritz method,in which the trial solution can be approximated to arbitrary accuracy by a suitable linear combination of a sufficiently large set of basis functions. However, we would like to emphasize that no matter what forms the trial solutions are assumed to be, they must satisfy the specified essential boundary conditions (unnecessary to satisfy the natural boundary conditions explicitly,because they are included intrinsically in the variational statement).

To be specific in elastic materials,the only state variable is the displacement fieldu(r)and the total free energy functional is?t=?t[u(r)].The trial approximate solution ofu(r)can be completely empirical, for example, taking the form ofU(r;c),which satisfies the specified essential boundary conditions and is parameterized byNas yet unknown (independent) parametersc=(c1,c2,...,cN). In addition, in the Ritz method,we seek a more explicit approximationU(r;c)in the form of finite series expansion of a set of continuous,linearly independentNbase functions that satisfy the specified essential boundary conditions of the problem. Note that,in general,the accuracy of either form of trial solutionsU(r;c) can not be knownin priori; it is known only after we make the calculations using the trial solution and compare the results with experiments and/or simulations/computations.This,therefore,requires more efforts made to choose a reasonable trial solution, which is chosen not arbitrarily but should be based on enough knowledge and deep understanding of the system and the problem that are gained from abundant experimental measurements and numerical simulations/computations.

Substituting either form of the approximate trial solutionsU(r;c) into the total free energy functional?t[u(r)],we obtain(after carrying out the integration with respect tor):?t=?t(c). The parameterscare then determined by minimizing?twith respect toc:

which represents a set ofNlinear equations amongc1,c2,...,cN,whose solution together with the above trial solution yields the approximate solutionU(r). This completes the description of the variational methods of approximation.

3. Transmission of forces induced by a spherically contracting cell in biopolymer gels

Now we use the above variational methods to study the decay of displacementsu(r)induced by a spherically contracting cell that is well adhered to a 3D nonlinear elastic biopolymer gel.[12,19,36]In the experiments for a fibroblast cell contracting in a 3D fibrin gel[16](as shown in Fig.1),the decay of the cell-induced displacement over distance is found to follow a power-law scaling asu～r-0.52,which is much more slowly than that induced by cells in linear elastic medium (scaling asu～r-2). In this section, we show that these experimental observations can be well understood and explained by the three-chain model of biopolymers and direct approximate variational analysis.

An adherent cell can apply active contractions to their surrounding matrix. In this sense, the cell can be modeled as a contractile force dipole.[37]In the simplest case, we here consider the decay of displacements that are induced by a spherically contracting cell in a 3D infinite extracellular matrix[19,38,39]as performed in experiments and schematically shown in Fig. 1(a). In this case, the active cell contraction can be characterized by a boundary condition of fixed radial displacement-ucat the cell boundaryr=Rcas shown in Fig.1(a),i.e.,

withuc＞0 for contractile cells. We would like to point out that a fibrous extracellular matrix with strong crosslinkers and pore sizes much smaller than the cell dimensions can be modeled as a continuum elastic material as described by the energetic models presented in Section 2. Then,in any matrix element of size much larger than the average pore size of the matrix gel but much smaller than the cell dimensions,the matrix deformation can be assumed to be homogeneous and affine as schematically shown in Figs. 2(c1) and 2(c2). Therefore, the cell-matrix mechanical interaction yields an elastic boundary value problem with inhomogeneous deformation in the extracellular biopolymer gel that depends on the radial coordinate,r,in a spherical geometry. In this case,the principal directions of the strain and stress tensor of each gel element are along the radial and the two perpendicular angular directions. This spherical symmetry significantly simplifies our analysis without losing any important physical mechanisms. In spherical coordinates(r,θ,φ),the displacement vector is then given byu=u(r)?r, the non-zero components of the strain tensorεijare

3.1. Force transmission in biopolymer gels under linear limits

We first calculate the decay of displacements or force transmission in the following two linear limits.

Force transmission in linear isotropic limit:ε1/εs?1 and|ε2,3|/εb?1. In this limit, the elastic energy density in Eq.(3)reduces to

which takes the same form ofFefor linear isotropic materials[12,13]if substitutingμ0=E0/2(1+ν0) andK0=E0/3(1-2ν0). To use the Ritz-type variational method to calculate the decay of cell-induced displacements in such biopolymer gels under the linear isotropic limits, we first need to find a good trial solution. From the theory of linear elasticity,[31]we know that the Green functions of point forces and force dipoles follow power-law scaling,which is like that of point charges and electric dipoles in dielectric medium. It is, therefore, reasonable to take the power-law form of empirical trial solution,u=-uc(Rc/r)n, which satisfies the essential boundary conditions atr=Rcandr →∞. Note that here in the empirical trial solution, there is only one undetermined parametern; however,if we know very little about the physics of the system and the problem, we need to use the original Ritz method,in which we need to take a trial solution in the form of finite series expansion of a set of continuous,linearly independent base functions(such as polynomial or sinusoidal/cosine functions). In that case, there can be many undetermined parameters and the approximate solutions will depend on the number of base functions used in the trial solution.

3.2. Force transmission in nonlinear biopolymer gels

We now calculate the decay of displacements in nonlinear elastic biopolymer gels.

Fig. 3. An approximating power-law decay of cell-induced displacements in nonlinear gels of compressive-softening but no stretchstiffening. (a) Elastic deformation energy F as a function of the power-law exponent n. Several different nonlinearity parameters Ab ≡uc/Rcεb have been chosen.The black dots in the curves denote the optimized exponent n* that corresponds to minimum energy. Here we take ν0 =0.01, uc/Rc =0.2, ρ0 =0.1. (b) The optimized exponent n* is plotted as a function of Ab for different Poisson’s ratios. Here we also take uc/Rc=0.2 and ρ0=0.1.

Here we do the integral over a finite matrix region,Rc≤r ≤Rm, takeRm=10Rc＞Rb, and we plot the energy?(n) in Fig.3(a)for different nonlinearity parameters,Ab=uc/Rcεb.We then minimize?(n)numerically with respect tonand obtain the optimized power-law exponentn*as a function ofεb,ν0, andρ0. In Fig. 3(b),n*is plotted as a function ofAbfor different Poisson’s ratios. In the limit ofAb→∞, the biopolymer gels behave simply as linear anisotropic material and the correspondingn*are marked by open symbols. Asν0increases to the incompressible limitν0=0.5,n*approaches to its value in linear isotropic materials,n*=2.

Fig. 4. An approximating power-law decay of cell-induced displacements in nonlinear gels with both compressive-softening and stretchstiffening. (a)Elastic deformation energy F as a function of the powerlaw exponent n. Several different stiffening nonlinearity parameters As ≡uc/Rcεs have been chosen. The black dots in the curves denote the optimized exponent n* that corresponds to minimum energy. Here we take ν0=0.4,uc/Rc=0.2,ρ0=0.1,and Ab=2.0(with εb=0.1).(b)The optimized exponent n* is plotted as a function of As for different buckling nonlinearity parameter Ab. Here we also take uc/Rc=0.2 and ρ0=0.1.

Force transmission in biopolymer gels with both stretch-stiffening and compressive-softening. In general,one must consider the full elastic energy densityFein the form of Eq. (3), which takes into account of both stretchstiffening and compressive softening. In this case, we still take the approximate empirical trial solution,u=-uc(Rc/r)n,and numerically integrate the deformation energy functional?[u(r)] over a finite matrix region,Rc≤r ≤Rm, and obtain an energy function ofn:?(n;εs,εb,ν0,ρ0). In Fig.4(a),we plot the energy?(n)for different nonlinearity parameters,As=uc/Rcεs, where we takeRm=10Rc＞Rbas measured in experiments.[16]We then minimize?(n)numerically with respect tonand obtain the optimized power-law exponentn*as a function ofεs,εb,ν0, andρ0. In Fig. 4(b), we plotn*as a function ofAsfor variousAb. In the limit ofAs→0,the biopolymer gels behave in the same way as nonlinear gels with only compressive softening(as shown in Fig.3)and the correspondingn*are marked by open symbols in Fig.4(b).

In addition,we find that the slow power-law decay of displacements induced by contracting cells,u ≈-uc(Rc/r)0.52,measured in experiments by Notbohmet al.[16]can be quantitatively explained by the full three-chain model of biopolymer gels and the above variational approximations. For this purpose,we first estimate the order of magnitude of all the model parameters that have clear physical meaning and have been measured in manyin vitroexperiments of cells and biopolymer gels. For example, in fibrin gels, the critical strains,εs(for the presence of stretch stiffening)andεb(for the presence of compressive softening),are measured to be around 10%[4–6]and 1%,[7,8]respectively. Poisson’s ratio of fibrin gels ranges from 0.1 to 0.4.[12,41]Therefore, for a typical cell of contraction strengthuc/Rc～10%, we obtain the nonlinear parametersAs≡uc/Rcεs～1 andAb≡uc/Rcεb～10. Based on the experimental observation on the power-law decay of cellinduced displacements and the understanding on magnitudes of the parameters,we then take the approximate trial solution to be the power-law form ofu=-uc(Rc/r)n, fixν0=0.1,and choose the small gel-softening ratioρ0to be 0.1. We find that the experimental data of cell-induced displacements and the effective power-law exponentn ≈0.52 from Ref.[16]can be fitted very well(see Fig.5)using the three-chain model by takingAs=1.76 andAb=33.3, which are in the same order of magnitude as that estimated above. Furthermore, we have also numerically solved the equilibrium equation(9)derived from the full three-chain model of biopolymer gels, in which all the parameters take the same values as that in our previous approximate variational analysis. The numerical results in Fig.5 also show very good agreement with the experimental data. We can, therefore, conclude that the decay of cell-induced displacements in biopolymer gels can be well explained by the three-chain model of biopolymers,and the slow power-law decay with an exponent～0.52 measured in experiments is an effective power-law that can be further understood by direct approximate variational analysis.

Fig. 5. The decay of displacement induced by the spherically contracting cell in nonlinear biopolymer gels with both compressive softening and stretch stiffening.The approximating power law u=-uc(Rc/r)-0.52(shown by dashed lines) measured in experiments[16] is fitted here by choosing ν0 =0.1,ρ0 =0.1,Ab =33.3, and As =1.76, which lie in the reasonable range of these parameters measured in experiments. The corresponding numerical solution of the equilibrium equation (9) is shown by solid curves.Open symbols with different colors are taken from experiments shown in Fig.1(b)by Notbohm et al.[16]

4. Conclusions

Variational methods have been widely used in the modeling and analysis of soft and biological matter. We utilize the variational methods that are based on the variational principle of minimum free energy (MFEVP) to study the slow-decay of the displacements induced by a spherically contracting cell in a 3D biopolymer gel. We employ the three-chain model to describe the nonlinear(stretch-stiffening and compressivesoftening)elasticity of biopolymer gels. We firstly show that in the linear limits, the classical scaling laws for the decay of cell-induced displacements in linear isotropic and linear anisotropic elastic medium can both be obtained exactly by the Ritz-type variational methods. We then show that in general nonlinear biopolymer gels,when relevant physical parameters take reasonable values as that measured in other separate experiments,we can use variational methods to reproduce the scaling lawu～r-0.52that is measured experimentally for the decay of displacements induced by fibroblast cell spheroid contracting in 3D fibrin gels. This work evidences the validity of continuum modeling in describing fibrous biopolymer gels and deepens our understanding of long-range force transmission in biopolymer gels and extracellular matrix that is essential for efficient matrix-mediated cell-cell communications.

Acknowledgments

X. X. is supported by the National Science Foundation for Young Scientists of China (Grant No. 12004082),Guangdong Province Universities and Colleges Pearl River Scholar Funded Scheme (2019), 2020 Li Ka Shing Foundation Cross-Disciplinary Research (Grant No. 2020LKSFG08A), Provincial Science Foundation of Guangdong (Grant No. 2019A1515110809), Guangdong Basic and Applied Basic Research Foundation (Grant No.2020B1515310005), Featured Innovative Projects(Grant No.2018KTSCX282),and Youth Talent Innovative Platforms(Grant No. 2018KQNCX318) in Universities in Guangdong Province.