# Numerical Modeling of Force-Stiffness Response of Cross-Linked Actin Networks Using Tensegrity Systems

Authors: Xian Xu; Yafeng Wang; Yaozhi Luo
Journal: Mathematical Problems in Engineering
 (2015)

---

## Abstract

A three-dimensional tensegrity structure is used as a computational model for cross-linked actin networks. The postbuckling behavior of the members under compression is considered and the constitutive relation of the postbuckling members is modeled as a second-order polynomial. A numerical scheme incorporating the equivalent constitution of the postbuckling members is used to predict the structural response of the tensegrity model under compression loads. The numerical simulation shows that the stiffness of the tensegrity structure nonlinearly increases before member buckling and abruptly decreases to a lower level as soon as members buckle. This result qualitatively mimics the experimentally observed stiffness to compression stress response of cross-linked actin networks. In order to take member length variety into account, a large number of simulations with the length of buckling members varying in the given range are also carried out. It is found that the mean response of the simulations using different buckling member length exhibits more resemblance to the experimental observation.


---

## Body

## 1. Introduction

Tensegrity is a special class of pin-jointed assemblies, whose stability is provided by the self-stress state between tensioned elements and compressed elements. Sculptures containing the essential characteristics of tensegrity were first created by Snelson in 1948, while the term “tensegrity” was first used by Fuller in 1961 to name a class of cable-bar structures in a patent [1, 2]. Since it was invented a half century ago, tensegrity has attracted much attention from both researchers and practicers. Architectural application was believed as a promising area for tensegrity applications. But after many endeavors towards this direction, applications of this structural system in real projects are very few. The limited use of the system can be attributed to its inherently complex and flexible nature. However, it inspired some innovative structural systems such as cable domes [3] and suspended-domes [4], which have attained great success in architectural applications. Currently, new applications of tensegrity structures are taking place in the field of “smart” structures [5–10], locomotive robots [11], and cellular mechanics [12–15], due to their adjustability and resemblances to cellular structures. A comprehensive review and discussion on the development and applications of tensegrity can be found in Motro [16], Skelton and de Oliveira [17], and Sultan [18].Tensegrity-based models have been used to successfully predict the mechanical behaviors of whole cells [12, 13, 19–22], the erythrocyte membrane [23], viruses [24], actin stress fibers [25], and actin networks [26]. Postbuckling behavior of compression elements especially has been considered to explain the linear stiffening behavior of living cells [27, 28]. Stretch force [27] and shear force [28] were considered and emphasis was put on the linear stiffening behavior exhibited after the buckling of compression elements. The behavior before and during the buckling of compression elements was short of analysis. The application of tensegrity structures as mechanical models of living cells is based on the assumption that the cytoskeleton enjoys tensegrity structures with the microtubules and microfilaments acting as compression elements and tension elements, respectively. Since living cells are initially tensed, prestress is always included in the tensegrity model of the cytoskeleton. Though tensegrity models with zero prestress were also examined in the previous studies, they were treated as comparison cases to the ones with prestress [27, 28]. Xu and Luo [26] used a planar 3-member tensegrity structure with a buckling member to explain the nonlinear stiffening and reversible softening behavior exhibited by dendritic actin networks. Compression load which has not been considered in the previous studies is applied to the tensegrity model to mimic the compression imposing on the actin networks [26]. Since the average stiffness of dendritic actin networksin vitro was found to be independent of prestressing by myosin II motors [29], no prestress was introduced into the tensegrity model of actin networks [26].Based on the previous study [26], a more general tensegrity model in three dimensions is proposed in this paper. Instead of the analytical method used for the planar model, a numerical approach is used to predict the structural response of the 3D tensegrity model and the constitutive relation of the postbuckling members is modeled as a second-order polynomial. Numerical results show that the stiffness to compression force response of the proposed model qualitative corresponds to the behavior exhibited by the actin networks under compression. The work presented in this paper confirms the finding reported in the previous study [26] and extends the tensegrity model from a conceptual 2D system to a more general 3D system. The analytical method used in the previous study [26] is only fit for simple tensegrity system with a very small number of members. The numerical approach presented in this paper is a general scheme which is able to handle more complex tensegrity system with a large number of members and can be used in the future study of using more sophisticated tensegrity model to obtain quantitative correspondence to the experimental observation.


## 2. Methods

### 2.1. Derivation of Stiffness

Here, the basic stiffness formulation described in Guest [30] is used. The tangent stiffness matrix K relates to the nodal displacement, written as a vector u, and the applied load at each node, written as a vector P, as(1)
K
u
=
Pwhich is well known as the equilibrium equation in displacement form. The equilibrium of the structural system also can be described in force form as(2)
B
T
t
=
P
,where B is the compatibility matrix of the structural system and t is the internal force matrix [31]. The compatibility matrix B can be determined from kinematic relationships e
=
e
[
u
] by using the principle of virtual displacements [28](3)
B
=
∂
e
∂
u
,where e is the elongation matrix of the members.The tangent stiffness of the system can be obtained by(4)
K
=
∂
P
∂
u
.Substituting (2) and (3) into (4) yields(5)
K
=
B
T
C
B
+
D
u
,where(6)
C
=
∂
t
∂
e
,
D
=
∂
B
T
∂
u
t
.The tangent stiffness matrix also can be written as [32, 33](7)
K
=
A
G
^
Α
T
+
S
,where A
=
B
T; A
G
^
Α
T is called the “material stiffness matrix”; and S is called the “geometric stiffness matrix.” G
^ is a diagonal matrix of modified axial stiffness, with an entry for each member i (1
≤
i
≤
m, and m is the number of members of the structural system),(8)
g
^
i
=
g
i
1
-
ε
i
,where g
i
=
A
i
E
i
/
l
i is the axial stiffness of member i and ε
i
=
t
i
/
(
A
i
E
i
) is the nominal strain of member i. A
i, E
i, and l
i are defined as the cross-sectional area, Young’s modulus, and length of member i, respectively.The geometric stiffness matrixS can be written as the Kronecker product of stress matrix Ω and a 3D identity matrix I
; that is,(9)
S
=
Ω
⊗
I
.The components of matrix Ω are given by(10)
Ω
i
j
=
-
t
^
i
,
j
=
-
t
^
j
,
i
if
i
≠
j
,
i
,
j
a member
∑
k
≠
i
t
^
i
k
if
i
=
j
0
if there is no connection between nodes
i
and
j
,where t
^
i
,
j is the tension coefficient (
t
/
l
) in the member connecting nodes i and j.


### 2.2. Constitutive Relation

A constitutive relation as that used by Volokh et al. [28] for the strut of their tensegrity models for the cell cytoskeleton is adopted for all the members of the tensegrity system considered in this paper. We do not intentionally sort the members of a tensegrity system as cables or struts, because actin networks consist of one type of elements: actin filaments. Each member in the tensegrity system is modeled as a pinned elastic column (Figure 1). The rest length of the elastic column is notated as L and the bending stiffness of it is notated as E
I, where E is Young’s modulus of actin filaments and I is the moment of inertia of the sectional area of an actin filament. When it is tensioned, it behaves linearly in accordance with Hooke’s law. When it is compressed, it also behaves linearly in accordance with Hooke’s law as long as the compression in it is smaller than the Euler critical force P
c
r. Once the compression in the member reaches P
c
r, the member begins to buckle. The postbuckling behavior of a member is assumed to be in accordance with the solution of the “elastica” problem [34](11)
P
P
c
r
=
4
π
2
K
sin
⁡
α
2
2
,
(12)
l
α
=
L
4
π
P
c
r
P
E
sin
⁡
α
2
-
1
,where P
c
r
=
π
2
E
I
/
L
2; l is the chord length after the buckling (Figure 1); K
[
·
] and E
[
·
] are complete elliptic integrals of the first and second kind, respectively; and α is the end slope of the deformed bar (Figure 1).Figure 1
Postbuckling deformation of an elastic bar pinned at the endpoints under axial forceT; l is the chord length; α is the end slope of the deformed bar.The nominal strain of the bar during postbuckling is given by(13)
ε
=
l
-
L
L
.Substituting (11) into (12) and then substituting (12) into (13), we have(14)
ε
=
2
E
sin
⁡
α
/
2
K
sin
⁡
α
/
2
-
1
.Equation (11) can be reexpressed as(15)
σ
σ
c
r
=
4
π
2
K
sin
⁡
α
2
2
,where σ is the nominal stresses during the postbuckling and σ
c
r
=
P
c
r
/
A. Equations (14) and (15) can be regarded as parametric equations describing the normalized constitutive relation of a bar during postbuckling. Least square method is employed to find a polynomial expression for the relation of σ
/
σ
c
r
~
ε. A second-order polynomial which is sufficient to give a good accuracy [28] is used as a seed equation for the constitutive relation. It yields that(16)
σ
σ
c
r
=
a
ε
2
+
b
ε
+
c
,where a
=
-
0.8227, b
=
0.3029, and c
=
-
0.9944. As shown by Figure 2, (16) is fit well with the solution of the “elastica” problem. More details on the approximating of the postbuckling behavior by polynomial can be found in Volokh et al. [28]. And the full constitutive relation for a bar before and after buckling is(17)
σ
=
E
ε
,
σ
≤
σ
c
r
σ
c
r
a
ε
2
+
b
ε
+
c
,
σ
>
σ
c
r
.Figure 2
Second-order polynomial fitting of normalized constitutive relation of a bar during postbuckling.


### 2.3. Numerical Simulations

Based on the above analysis, it is obvious that the tangent stiffnessK depends on the nodal displacement u. Since the member elongation e depends on the member constitutive relationships, the tangent stiffness K also depends on the member constitutive relationships in accordance with (6). Hence, (1) is nonlinear. To solve such a nonlinear equation, iterations are introduced. Equation (1) can be rewritten in an increment form as(18)
K
i
δ
u
i
=
p
-
r
i
,
u
i
+
1
=
u
i
+
δ
u
i
,
u
0
=
0
,where K
i can be estimated by substituting u
i into (5) and r
i
=
B
t
[
u
i
] is the residual force. Then it can be iteratively solved by a Newton-Raphson scheme.Here we consider a 30-member tensegrity structure which has been intensively investigated as a model of the cytoskeleton. Three nodes in the plane ofx
y are jointly fixed to a rigid base, and another three nodes at z
=
L
3
/
2 are subjected to forces directed to the negative direction of the z-axis. Under the given boundary condition and loads, 6 members of it are in compression; 21 members of it are in tension; and the remainder are in zero stresses, as shown by Figure 3. The rest lengths of the compression members and the tension members are L and L
6
/
4, respectively. A commercial FEM package ANSYS 9.0 is used to carry out the simulation. Each member of the tensegrity is mimicked by a LINK8 element. The nonlinear elastic constitutive relation described in (17) is input to the FEM model by user defined material. Incremental loading is used and the embedded Newton-Raphson solver is chosen to solve the iteration equations at every load step.Figure 3
30-member tensegrity model (thick lines: members in compression; thin lines: members in tension; and dashed lines: members with zero stress).


### 2.4. Scales and Mechanical Properties of Actin Filaments

Experimentally obtained scales and mechanic properties for actin filaments are used in the simulations. By analyzing the thermally driven fluctuations in their shape, the bending stiffness and Young’s modulus of actin filaments are estimated to beE
I
=
7.29
×
10
-
26 Nm2 and E
=
2.6 Gpa [35]. The cross section area of the actin filament is approximately A
=
18.85 nm2, based on the assumption that the cross section of the actin filament has an ellipse shape with major and minor radii of 3 nm and 2 nm [35]. The length of filaments in dendritic actin networks is reported to be less than 1 μm [36]. Without loss of generality, values in the range of 0.5–1 μm are used for the representative member length L in our calculations.


## 2.1. Derivation of Stiffness

Here, the basic stiffness formulation described in Guest [30] is used. The tangent stiffness matrix K relates to the nodal displacement, written as a vector u, and the applied load at each node, written as a vector P, as(1)
K
u
=
Pwhich is well known as the equilibrium equation in displacement form. The equilibrium of the structural system also can be described in force form as(2)
B
T
t
=
P
,where B is the compatibility matrix of the structural system and t is the internal force matrix [31]. The compatibility matrix B can be determined from kinematic relationships e
=
e
[
u
] by using the principle of virtual displacements [28](3)
B
=
∂
e
∂
u
,where e is the elongation matrix of the members.The tangent stiffness of the system can be obtained by(4)
K
=
∂
P
∂
u
.Substituting (2) and (3) into (4) yields(5)
K
=
B
T
C
B
+
D
u
,where(6)
C
=
∂
t
∂
e
,
D
=
∂
B
T
∂
u
t
.The tangent stiffness matrix also can be written as [32, 33](7)
K
=
A
G
^
Α
T
+
S
,where A
=
B
T; A
G
^
Α
T is called the “material stiffness matrix”; and S is called the “geometric stiffness matrix.” G
^ is a diagonal matrix of modified axial stiffness, with an entry for each member i (1
≤
i
≤
m, and m is the number of members of the structural system),(8)
g
^
i
=
g
i
1
-
ε
i
,where g
i
=
A
i
E
i
/
l
i is the axial stiffness of member i and ε
i
=
t
i
/
(
A
i
E
i
) is the nominal strain of member i. A
i, E
i, and l
i are defined as the cross-sectional area, Young’s modulus, and length of member i, respectively.The geometric stiffness matrixS can be written as the Kronecker product of stress matrix Ω and a 3D identity matrix I
; that is,(9)
S
=
Ω
⊗
I
.The components of matrix Ω are given by(10)
Ω
i
j
=
-
t
^
i
,
j
=
-
t
^
j
,
i
if
i
≠
j
,
i
,
j
a member
∑
k
≠
i
t
^
i
k
if
i
=
j
0
if there is no connection between nodes
i
and
j
,where t
^
i
,
j is the tension coefficient (
t
/
l
) in the member connecting nodes i and j.


## 2.2. Constitutive Relation

A constitutive relation as that used by Volokh et al. [28] for the strut of their tensegrity models for the cell cytoskeleton is adopted for all the members of the tensegrity system considered in this paper. We do not intentionally sort the members of a tensegrity system as cables or struts, because actin networks consist of one type of elements: actin filaments. Each member in the tensegrity system is modeled as a pinned elastic column (Figure 1). The rest length of the elastic column is notated as L and the bending stiffness of it is notated as E
I, where E is Young’s modulus of actin filaments and I is the moment of inertia of the sectional area of an actin filament. When it is tensioned, it behaves linearly in accordance with Hooke’s law. When it is compressed, it also behaves linearly in accordance with Hooke’s law as long as the compression in it is smaller than the Euler critical force P
c
r. Once the compression in the member reaches P
c
r, the member begins to buckle. The postbuckling behavior of a member is assumed to be in accordance with the solution of the “elastica” problem [34](11)
P
P
c
r
=
4
π
2
K
sin
⁡
α
2
2
,
(12)
l
α
=
L
4
π
P
c
r
P
E
sin
⁡
α
2
-
1
,where P
c
r
=
π
2
E
I
/
L
2; l is the chord length after the buckling (Figure 1); K
[
·
] and E
[
·
] are complete elliptic integrals of the first and second kind, respectively; and α is the end slope of the deformed bar (Figure 1).Figure 1
Postbuckling deformation of an elastic bar pinned at the endpoints under axial forceT; l is the chord length; α is the end slope of the deformed bar.The nominal strain of the bar during postbuckling is given by(13)
ε
=
l
-
L
L
.Substituting (11) into (12) and then substituting (12) into (13), we have(14)
ε
=
2
E
sin
⁡
α
/
2
K
sin
⁡
α
/
2
-
1
.Equation (11) can be reexpressed as(15)
σ
σ
c
r
=
4
π
2
K
sin
⁡
α
2
2
,where σ is the nominal stresses during the postbuckling and σ
c
r
=
P
c
r
/
A. Equations (14) and (15) can be regarded as parametric equations describing the normalized constitutive relation of a bar during postbuckling. Least square method is employed to find a polynomial expression for the relation of σ
/
σ
c
r
~
ε. A second-order polynomial which is sufficient to give a good accuracy [28] is used as a seed equation for the constitutive relation. It yields that(16)
σ
σ
c
r
=
a
ε
2
+
b
ε
+
c
,where a
=
-
0.8227, b
=
0.3029, and c
=
-
0.9944. As shown by Figure 2, (16) is fit well with the solution of the “elastica” problem. More details on the approximating of the postbuckling behavior by polynomial can be found in Volokh et al. [28]. And the full constitutive relation for a bar before and after buckling is(17)
σ
=
E
ε
,
σ
≤
σ
c
r
σ
c
r
a
ε
2
+
b
ε
+
c
,
σ
>
σ
c
r
.Figure 2
Second-order polynomial fitting of normalized constitutive relation of a bar during postbuckling.


## 2.3. Numerical Simulations

Based on the above analysis, it is obvious that the tangent stiffnessK depends on the nodal displacement u. Since the member elongation e depends on the member constitutive relationships, the tangent stiffness K also depends on the member constitutive relationships in accordance with (6). Hence, (1) is nonlinear. To solve such a nonlinear equation, iterations are introduced. Equation (1) can be rewritten in an increment form as(18)
K
i
δ
u
i
=
p
-
r
i
,
u
i
+
1
=
u
i
+
δ
u
i
,
u
0
=
0
,where K
i can be estimated by substituting u
i into (5) and r
i
=
B
t
[
u
i
] is the residual force. Then it can be iteratively solved by a Newton-Raphson scheme.Here we consider a 30-member tensegrity structure which has been intensively investigated as a model of the cytoskeleton. Three nodes in the plane ofx
y are jointly fixed to a rigid base, and another three nodes at z
=
L
3
/
2 are subjected to forces directed to the negative direction of the z-axis. Under the given boundary condition and loads, 6 members of it are in compression; 21 members of it are in tension; and the remainder are in zero stresses, as shown by Figure 3. The rest lengths of the compression members and the tension members are L and L
6
/
4, respectively. A commercial FEM package ANSYS 9.0 is used to carry out the simulation. Each member of the tensegrity is mimicked by a LINK8 element. The nonlinear elastic constitutive relation described in (17) is input to the FEM model by user defined material. Incremental loading is used and the embedded Newton-Raphson solver is chosen to solve the iteration equations at every load step.Figure 3
30-member tensegrity model (thick lines: members in compression; thin lines: members in tension; and dashed lines: members with zero stress).


## 2.4. Scales and Mechanical Properties of Actin Filaments

Experimentally obtained scales and mechanic properties for actin filaments are used in the simulations. By analyzing the thermally driven fluctuations in their shape, the bending stiffness and Young’s modulus of actin filaments are estimated to beE
I
=
7.29
×
10
-
26 Nm2 and E
=
2.6 Gpa [35]. The cross section area of the actin filament is approximately A
=
18.85 nm2, based on the assumption that the cross section of the actin filament has an ellipse shape with major and minor radii of 3 nm and 2 nm [35]. The length of filaments in dendritic actin networks is reported to be less than 1 μm [36]. Without loss of generality, values in the range of 0.5–1 μm are used for the representative member length L in our calculations.


## 3. Results

Six typical values in whichL = 1.0, 0.9, 0.8, 0.75, 0.7, 0.6, and 0.5 μm for the member rest length are considered. The relationships of structural stiffness versus the external force for the abovementioned three cases are given in Figure 4. It is clearly shown that all the tensegrity models with different rest length of members present a nonlinear stiffening response at the beginning region and then suddenly enter a process of stress softening. The turnover point between the nonlinear stiffening process and the stress softening process is corresponding to the Euler buckling of the member under compression. Since shorter members have larger Euler critical buckling forces, it is reasonable to expect that models with shorter members have larger stiffness than those with longer members. This has been verified by the simulation results (Figure 4). And it is also shown that the magnitude of the rest member length does not affect the shape of K versus P curve but the rise and span of it.Figure 4
Numerical solution of structural stiffnessK versus compression force P relationship using typical rest length for members in compression.Comparing Figure4 to the reported stiffness versus stress data of actin networksin vitro [29], it is notable that the softening process exhibited by the tensegrity models covers a smaller regime than that observed in actin networks. The possible reason for this is that in the tensegrity model we considered the members in compression are identical to each other and buckle simultaneously under the given load, while in the actin networks there are numerous filaments with various lengths under compression, which buckle sequently as the load increases. It seems that the tensegrity models considered above only mimic a basic piece of structure of the actin network.Therefore, a complex tensegrity system with various member lengths should be a more accurate model for the actin network. But constructing a tensegrity structure with a large number of members of various lengths is an unsettled problem as yet. Considering that the behavior of modular tensegrity systems is similar to the single cells [37], we assume that the behavior of a complex tensegrity system with various member lengths is the superposition of the behavior of single cells with various member lengths. 500 calculations are conducted for the 30-member tensegrity structures with L uniformly varied in the range of 0.5–1.0 μm. Then the means of the structural stiffness obtained from each of the 500 calculations are computed. In the plot of the mean structural stiffness versus compression force relationship the softening process has been extended to a regime much larger than that exhibited by the single tensegrity cells, as shown by Figure 5. It also shows that the compression force at the end of the softening regime is about twice that at the end of the stiffening regime, which is in good agreement with the corresponding observation of the actin network.Figure 5
Mean structural stiffnessK versus compression force P relationship for 500 30-member tensegrity cells with rest member length L uniformly ranging from 0.5 μm to 1.0 μm.


## 4. Conclusions

This paper presents a numerical study on simulating the mechanical response of cross-linked actin networks by a 3D tensegrity model. A general numerical scheme considering the postbuckling behavior of members under compression is used to carry out the simulation. The simulation result shows that the nonlinear stiffening and softening behavior exhibited by cross-linked actin networks has been qualitatively predicted by the proposed tensegrity model. Due to the symmetry of the model, the members under compression are identical to each other and buckle simultaneously and thus the stiffness of the tensegrity model abruptly decreases to a lower level as soon as members buckle. To take member length variety into account, a large number of simulations with the member length varying in the given range are carried out. It is found that the mean response of the simulations exhibits more resemblance to the experimental observation. This finding indicates that cross-linked actin networks may work as a large-scale tensegrity system with a large number of members in various lengths. It is also worth noting that the work presented in this paper is still a qualitative study and no quantitative correspondence has been reached yet. To construct a large-scale and sophisticated tensegrity model that could obtain quantitative mimic to the cross-linked actin network will be the future work of the authors.



---

*Source: 182918-2015-12-22.xml*