The eukaryotic cytoskeleton is a protein fibre network mainly consisting of the semi-flexible biopolymer F-actin, microtubules and intermediate filaments. It is well known to exhibit a pronounced structural polymorphism, which enables intracellular processes such as cell adhesion, cell motility and cell division. We present a computational study on cross-linked networks of semi-flexible polymers, which offers a detailed analysis of the network structure and phase transitions from one morphology to another. We elaborate the morphological differences, their mechanical implications and the order of the observed phase transitions. Finally, we present a perspective on how the information gained in our simulations can be exploited in order to build both flexible and accurate, microstructurally informed, homogenized constitutive models of the cytoskeleton.

1. Introduction

The mechanical properties and the shape of biological cells are governed by the cytoskeleton [1], which plays a crucial role in a variety of intracellular as well as intercellular processes such as cell division, cell migration or mechanosensing. The cytoskeleton is a biopolymer network constituted by actin filaments (F-actin), intermediate filaments and microtubules. These filaments are interconnected by various cross-linking molecules (henceforth abbreviated as linkers). Depending on the current need, cells are able to restructure their cytoskeleton quickly and extensively by means of different linkers [2], which, e.g. leads to the formation of actin stress fibres [3], filopodia [4] and sheet-like structures in lamellipodia [5]. In vitro, other morphologies such as homogeneous-isotropic network structures [6] or cluster networks [7] have been discovered. So far, four fundamentally different networks morphologies have been observed in experiments with semi-flexible biopolymer networks: homogeneous-isotropic networks, bundle networks, cluster networks and lamellar networks (figure 1).

Figure 1. Different F-actin network structures: isotropic-homogeneous actin/heavy-mero-myosin-network (a), actin/fascin-network consisting purely of bundles (b), actin/filamin cluster network (c), lamellipodium (d); images (a)–(c) provided by K. M. Schmoller and A. R. Bausch, Lehrstuhl für Zellbiophysik E27, Technische Universität München, and (d) by M. Vinzenz and J. V. Small, Institute for Molecular Biotechnology of the Austrian Academy of Sciences, Vienna, Austria.

When the morphology changes, the mechanical properties of the network change as well. For example, the dynamic mechanical behaviour of homogeneous-isotropic actin-HMM networks [8] is vastly different from that of actin–fascin bundle networks [9,10].

On the cellular and tissue scale, the cytoskeleton is often treated as an isotropic homogenized continuum employing a pre-fitted constitutive law [11,12]. More complex constitutive models incorporate stress fibre formation [13] and orientation [14]. Other, very simplistic models rely on representations of the subcellular architecture by means of discrete springs rather than a consistent continuum formulation [15]. A comprehensive overview on common modelling approaches is given in [16]. Despite the crucial importance of subcellular structures to cellular mechanics [17,18], and numerous efforts in modelling them computationally (e.g. [19]), all aforementioned approaches neglect the full complexity of the phase transitions which the cytoskeleton has to undergo in order to produce fundamentally different morphologies like actin stress fibres or sheet-like structures in lamellipodia.

To overcome this deficiency, a new class of homogenized constitutive models will be required that is based on detailed information about the microstructure of and transition dynamics between the different possible network morphologies of the cytoskeleton. Notable theoretical efforts [20] explore how cells control the morphology of their cytoskeleton by tuning its thermodynamic equilibrium state via parameters such as the binding energy and concentration of linkers or the filament length. Also, a number of computational studies have dealt with phase transitions [21–23]. A recent one based on a novel, highly efficient finite-element method for Brownian dynamics simulations [24,25] succeeded in establishing an equilibrium phase diagram incorporating all experimentally observed network morphologies in cross-linked, semi-flexible polymer networks [26].

However, the microstructure of and transition dynamics between these morphologies remain poorly understood. In this paper, we examine both in a comprehensive computational study. We interpret our results relying on well-known theoretical concepts and mathematical tools from condensed matter physics like the structure function and the order parameter. This analysis provides the essential geometric and thermodynamic information required to establish more realistic homogenized constitutive relations for the cytoskeleton. The first part of this paper briefly summarizes its computational basis (§2) and key concepts from condensed matter physics for the structural analysis of polymer networks (§3). Subsequently, we provide a short discussion of the equilibrium thermodynamics of cross-linked, semi-flexible polymer networks (§4). In the second part, we present our main findings. Section 5 conducts a detailed examination of the microstructure of different network morphologies and biologically relevant transitions between them. Finally, in §6, we summarize our findings and discuss them in view of their benefit to the formulation of improved homogenized constitutive models for the cytoskeleton.

2. Brownian dynamics finite-element simulations

Here, we intend to provide a brief overview of the Brownian dynamics finite-element approach used in this paper. The details can be found in [24,25].

The essential components of a biopolymer network are filaments, linkers and a fluid, in which the former two are suspended. Filaments may be modelled explicitly as one-dimensional continua, e.g. as Simo–Reissner beams [27] or Kirchhoff beams [28,29]. Linkers are small molecules that diffuse through the network and can establish typically non-covalent chemical bonds to up to two filaments, thereby creating mechanical inter-filament joints. Biological linker species such as fascin, α-actinin or filamin [2] drive morphological changes in biopolymer networks. We think of these linkers as short and very stiff rods and explicitly model them as beams as well. The fluid, in which filaments and linkers are immersed, is only modelled implicitly by means of effective hydrodynamic drag terms and stochastic excitations acting on both species.

This way, we arrive at the present computational framework to study phase transitions in polymer networks. We consider a cubic representative volume element of edge length H which comprises N_f polymer filaments and N_l linkers. Initially, at time t=0, both filaments and linkers are distributed randomly. The dynamics of filaments and linkers is governed by their respective equation of motion, which we will discuss in more detail below. Over the course of time, depending on the density and the properties of filaments and linkers, aggregates emerge, whose formation is driven solely by the Brownian motion and the chemical interaction between the two species. The result of this process of self-organization is one of four morphological archetypes (cf. figure 3), which represent four distinct thermodynamic phases. Subsequent changes (e.g. adding or removing a certain number of linkers) entail transitions between these phases (more details are given in the electronic supplementary material). Thus, by means of appropriate modifications of simulation parameters, our computational model enables a convenient and detailed microscale examination of thermodynamic phase transitions of cross-linked polymer networks. The present section briefly summarizes how to specify the precise equations of motion of filaments and linkers as well as the interactions between them.

We model a filament of length L as a Simo–Reissner beam with a circular cross section, whose centreline is parametrized by a curve parameter s∈[0;L]. By means of the functions x(s,t),θ(s,t), we describe the position and the cross-section orientation pseudo-vector of each point s at time t. For given initial and boundary conditions, we can compute x and θ from the balance of linear and angular momentum

\underset{f_{int}}{\underset{⏟}{f_{el} (s, t)}} - \underset{f_{ext}}{\underset{⏟}{(f_{stoch} (s, t) + f_{visc} (s, t) + f_{e} (x, s, t))}} = 0

2.1a

and

\underset{m_{int}}{\underset{⏟}{m_{el} (s, t) - x^{'} (s, t) \times q_{s} (s, t)}} - \underset{m_{ext}}{\underset{⏟}{(m_{stoch} (s, t) + m_{visc} (s, t) + m_{e} (x, s, t))}} = 0 .

2.1b

Internal contributions (⋅)_int comprise elastic forces and moments f_el and m_el (per unit length), representing the elasticity of the polymer filaments under deformation. Furthermore, q_s is the elastic section force. These quantities can be computed as pointed out in [27,30]. In the microscopic context, inertia is negligible. Instead, filament dynamics is governed by three types of external loads (⋅)_ext: drag forces f_visc and moments m_visc modelling the viscosity of the fluid in which the filaments are immersed, stochastic Brownian forces f_stoch and moments m_stoch modelling random excitations of the filaments by collisions with microscopic particles in the surrounding fluid, and additional external loads (⋅)_e, representing, e.g. contact interactions between filaments. The prime-operator denotes (⋅)′=∂(⋅)/∂s. From here on, we omit arguments, i.e. f_el(s,t) simply becomes f′_el. For simplicity, we make the reasonable and physically consistent assumption that the fluid is Newtonian and can be represented by the 3×3 translational and rotational friction coefficient tensors D_t and D_r. Thus, we can account for viscous drag by

f_{visc} = D_{t} \dot{x}

2.2a

and

m_{visc} = D_{r} \dot{θ} .

2.2b

Superscripted dots mark time derivatives. Elements of these matrices can be drawn from literature (e.g. [31]). Random collisions with surrounding microscopic particles cause a stochastic, Gaussian excitation with zero mean and its variance is linked to the drag coefficients by the so-called fluctuation–dissipation theorem [32], giving

⟨ f_{stoch} (s, t) ⟩ = 0, ⟨ f_{stoch} (s, t) \otimes f_{stoch} (s^{*}, t^{*}) ⟩ = 2 k_{B} T δ_{s s^{*}} δ_{t t^{*}} D_{t}

2.3a

and

⟨ m_{stoch} (s, t) ⟩ = 0, ⟨ m_{stoch} (s, t) \otimes m_{stoch} (s^{*}, t^{*}) ⟩ = 2 k_{B} T δ_{s s^{*}} δ_{t t^{*}} D_{r} .

2.3b

Here, 〈⋅〉 are mean values and δ_ss*:=δ(s−s*) and δ_tt*:=δ(t−t*) are Dirac functions. Temperature T is assumed constant, k_B is the Boltzmann constant. Additional external forces and moments may model, e.g. contact interactions between filaments. For computer simulations, the filament continuum is discretized by beam finite elements. Discrete internal and external forces vectors are computed on the basis of the weak form of (2.1) and a Galerkin discretization.

As shown in figure 2a, in our simulation framework, linkers can be found in three distinct chemical states: free, singly bound and doubly bound. Free linkers are subject to diffusion until they establish a chemical bond to a filament at so-called binding sites on the filaments in accordance to figure 2b, thus becoming singly bound. When binding to a second filament, a singly bound linker becomes doubly bound, i.e. a cross-link (figure 2c). Free linkers are treated as rigid particles whose centre of mass x_l diffuses according to the equation of motion

f_{visc, l} (x_{l}, {\dot{x}}_{l}, t) = f_{stoch, l} (x_{l}, t)

2.4

with the viscous drag force

f_{visc, l} = ζ_{l} {\dot{x}}_{l} .

2.5

In accordance with Stoke’s law [31], the effective friction coefficient of the linker is ζ_l=6πηR_l. The stochastic force exerted on a linker is given by

⟨ f_{stoch, l} ⟩ = 0, ⟨ f_{stoch, l} \otimes f_{stoch, l} ⟩ = 2 k_{B} T ζ_{l} δ_{t t^{'}} I

2.6

with I being the second-order identity tensor. Free linkers of length R_l can establish chemical bonds to binding sites within their reaction volume V_r, which comprises all points in the radius interval [(R_l−ΔR_l)/2;(R_l+ΔR_l)/2] around the linker’s centre of mass. The tolerance ΔR_l accounts for rapid thermal motions and deformations.

Figure 2. Simulated linkers and filaments (green). (a) Chemical states of a linker: free (purple), singly bound (blue) and doublybound (red). (b) The free linker at position x_l is can establish a connection with the binding site at x_A. (c) The singly bound linker at x_A can bind to a second binding site at x_B, creating a cross-link. (Online version in colour.)

Singly bound linkers can bond to binding sites on another filament within the radius interval [R_l−ΔR_l;R_l+ΔR_l] around the binding site to which they are attached. Typically, the molecular structure of linkers and filaments allows chemical bonds only if additional orientational constraints are satisfied. As in [26], we assume that the two filaments to be cross-linked must enclose an angle ϕ∈[ϕ₀−Δϕ;ϕ₀+Δϕ]. Cross-links, i.e. doubly bound linkers, are modelled as beams and are discretized with beam finite elements like the filaments.

Given that the geometrical constraints are met, the reaction kinetics of the filament–linker bond within small time intervals Δt can be modelled as Poisson processes with binding and unbinding probabilities

p_{on} = 1 - \exp (- k_{on} Δ t) and p_{off} = 1 - \exp (- k_{off} Δ t),

2.7

respectively, with rate constants k_on and k_off. Rate constants can be force-dependent but we found that this dependence only has very minor impact on the present problem, making the above simple model a valid choice.

3. Structural analysis of polymer networks

The characterization of the different network morphologies presented in §5 requires certain mathematical and analytical instruments, which we will briefly introduce in the following.

(a) Morphology-dependent local coordinate systems

All inhomogeneous network morphologies (clusters, bundles and lamellae) correspond to geometrical primitives (point, line and plane). Their characterization is thus facilitated by defining local coordinate frames (x₁,x₂,x₃) aligned with the respective geometric primitive in a suitable way. The origins of these coordinate frames are located at the respective morphology’s centre of mass s_l. Cluster networks are analysed using an arbitrarily orientated coordinate frame. For bundle networks, we define the x₁-direction to be orthogonal to the bundle axis. Finally, for the lamellae, the x₁-direction is chosen orthogonal to the plane in which they lie. The respective coordinate frames are shown as insets of figures 4b, 5b and 7b.

(b) Orientational distribution of filaments

The distribution of filament orientations is a characteristic measure for any kind of filamentous network. Using spherical coordinates (φ,ψ), these orientations may be described by means of the orientation density functions

\begin{aligned} ρ (φ) & = ⟨ \frac{1}{\sum_{n = 1}^{N_{f}} L_{n}} \sum_{n = 1}^{N_{f}} \int_{0}^{L_{n}} δ [φ - φ_{n} (s)] d s ⟩ \\ and ρ (ψ) & = ⟨ \frac{1}{\sum_{n = 1}^{N_{f}} L_{n}} \sum_{n = 1}^{N_{f}} \int_{0}^{L_{n}} δ [ψ - ψ_{n} (s)] d s ⟩, \end{aligned}

3.1

where L_n denotes the length of the nth filament, 〈⋅〉 the time average, δ the Dirac-delta function and φ_n(s) and ψ_n(s) the orientation angles of the tangent of the nth filament at the material point s. The azimuth angle φ∈[0;2π[ is measured against the x₂-direction, the polar angle ψ∈[0;π] against the x₁-direction. Functions ρ(φ) and ρ(ψ) are direct measures of the orientational anisotropy of filaments in the network.

(c) Orientation correlation among filaments

Linkers with a preferred binding angle tend to orientationally correlate filaments close to each other. Let the position of the mth binding site on the ith filament be x^[i][m]_b. Then,

I = {(i, m, j, n) : R_{l} - Δ R_{l} \leq ∥ x_{b}^{[i] [m]} - x_{b}^{[j] [n]} ∥ \leq R_{l} + Δ R_{l}, i \neq j}

3.2

is the set of all index tuples belonging to pairs of binding sites on different filaments, which are in the appropriate distance to be connected by a linker. Some of them may indeed be cross-linked, while others are not, e.g. due to the lack of a linker molecule in their vicinity. Let t^[i][m]_b be the unit filament tangent vector at a binding site. Then, the correlation of orientations of neighbouring filaments can be written as the normalized orientation correlation function (OCF)

O (θ) = \frac{1}{N_{O}} \sum_{(i, m, j, n) \in I} δ (θ - \cos^{- 1} | t_{b}^{[i] [m]} \cdot t_{b}^{[j] [n]} |),

3.3

where

N_{O}

produces the unit integral on the interval [0;π/2]. Given a network with perfectly random filament orientations, the linker-mediated interactions between neighbouring filaments are negligible. The distribution function O(θ) is then expected approach

O (θ) = \sin (θ), θ \in [0; π / 2]

due to the geometric preference for angles close to π/2 between randomly oriented straight lines. By contrast, in the case of strong, linker-mediated filament–filament interactions and linkers with a preferred binding angle ϕ, one expects a pronounced maximum at θ≈ϕ.

(d) Filament order parameter

The order parameter of a nematic liquid crystal is commonly given as a scalar value

S = ⟨ P_{2} (\cos ϑ) ⟩ = \frac{1}{2} ⟨ 3 {(\cos ϑ)}^{2} - 1 ⟩,

3.4

representing the ensemble average of the second Legendre polynomial P₂ at a specific point in time and can measure the order of a structure [33]. Usually, ϑ is the enclosed angle between the direction of rigid rod-like molecules and a global director. Whereas S=0 indicates uncorrelated angles, S>0 and S<0 indicate a predominance of small and large ϑ, respectively.

Since we consider assemblies of long semi-flexible polymers, equation (3.4) is altered in order to capture the varying orientation along the filament centreline. Instead to a global direction, we compare the local filament orientation to the orientations of all other filaments depending on their mutual distance d. To this end, the modified order parameter reads

\begin{aligned} S (d) & = \frac{1}{\sum_{i = 1}^{N_{f} - 1} (L_{i} \sum_{j = i + 1}^{N_{f}} L_{j})} ⟨ \sum_{i = 1}^{N_{f} - 1} \sum_{j = i + 1}^{N_{f}} \int_{0}^{L_{i}} \int_{0}^{L_{j}} δ (d - ∥ x_{i} (s_{i}) - x_{j} (s_{j}) ∥) (3 ϑ_{i j}^{2} - 1) {d s}_{i} {d s}_{j} ⟩, \end{aligned}

3.5

where x_(⋅)(s_(⋅)) is the spatial position of the (⋅)th filament and ϑ_ij=|t_i(s_i)⋅t_j(s_j)| with tangents t_(⋅)(s_(⋅)).

(e) Two-point density–density correlation function of linkers

The spatial distribution of the linkers is an important characteristic of the network and is quantified as follows. Let the position vectors of the N_db doubly bound linkers (i.e. cross-links) be x^k_db, k=1,…,N_db with components $x_{db, i}^{k}, i = 1, 2, 3$ . In [33], the number density operator

n (x_{i}) = \sum_{k = 1}^{N_{d b}} δ (x_{i} - x_{db, i}^{k})

3.6

is a measure for the spatial density of doubly bound linkers. The distribution of these linkers can be characterized by the two-point density–density correlation function (DDCF),

C_{i} (x_{i}, {\hat{x}}_{i}) = ⟨ n (x_{i}) n ({\hat{x}}_{i}) ⟩ = ⟨ \sum_{k = 1}^{N_{db}} \sum_{k^{'} = 1}^{N_{db}} δ (x_{i} - x_{i}^{k}) δ ({\hat{x}}_{i} - x_{i}^{k^{'}}) ⟩,

3.7

with non-identical linker positions x^k and x^k′ and two non-identical variable points in space x and

\hat{x}

. In slight variation of

C_{i} (x_{i}, {\hat{x}}_{i})

, we simply evaluate the mutual distance between two doubly bound linkers by means of the normalized two-point DDCF

C_{i} (d_{i}) = \frac{1}{N_{C}} ⟨ \sum_{k = 1}^{N_{db}} \sum_{k^{'} = 1}^{N_{db}} δ (d_{i} - | x_{i}^{k} - x_{i}^{k^{'}} |) ⟩

3.8

with

N_{C}

yielding the unity integral on [0;H].

(f) Structure function

The Fourier transform of the two-point DDCF

I_{x \hat{x}} (q_{i}) = \int e^{- i q_{i} (x_{i} - {\hat{x}}_{i})} C_{i} (x_{i}, {\hat{x}}_{i}) d x_{i} d {\hat{x}}_{i}

3.9

where imaginary unit

i = \sqrt{- 1}

is the so-called structure function. It can often be measured directly in diffraction experiments. The vector quantity q—with the components q_i—describes the difference between incoming and scattered waves. In the present case, we understand it as the frequency of spatial occurrences of detected particles. The structure function of the modified DDCF (3.8) is given by

I_{i} (q_{i}) = \int e^{- i q_{i} d_{i}} C_{i} (d_{i}) d d_{i} .

3.10

(g) Elastic network strain energy

Depending on the beam formulation, the internal elastic energy or strain energy consists of individual terms stemming from axial, bending, shear and torsional deformation. The sum E_int of these contributions of all filaments can be calculated from the known deformations and strain energy functions as described in [27].

(h) Periodic boundary conditions

We ran simulations with filaments and linkers confined in cubes of edge length H and periodic boundary conditions. Note in this case that for the correct evaluation of some of the introduced functions like the OCF, one has to consider not only the actually simulated volume but also its adjacent 26 periodic continuations.

4. Phase diagram of cross-linked semi-flexible polymer networks

In [26], we conducted a vast in silico study of the temporal evolution of the network structure, which led to the constitution of the phase diagram shown in figure 3. Here, the thermal equilibrium configuration is given depending on the preferred binding angle of the linkers ϕ₀ assuming Δϕ=π/16 (cf. §2 on linkers), and the normalized linker concentration

n_{l} = \frac{N_{l}}{N_{b}}

4.1

with N_l being the total number of linkers in the simulated volume and N_b being the total number of binding sites on filaments. Figure 3 depicts the different equilibrium phases: a homogeneous-isotropic phase for low linker concentrations, a bundle phase at high linker concentrations and small preferred binding angles ϕ₀, a cluster phase at intermediate linker concentrations and large preferred binding angles ϕ₀, and a lamellar phase with orthogonal filament order (→squaratic) for ϕ₀ close to π/2 and hexagonal filament order (→hexatic) for ϕ₀ close to π/3. We explored the influence of the different parameters of linkers and filaments, in particular, of the number of filaments N_f (i.e. filament concentrations c_f), their persistence lengths l_p, their axial stiffnesses E_fA_f, and the stiffness and length of the linkers. We found quantitative differences, e.g. as regards the critical linker concentrations n_l,crit at which phase transitions occur, but no qualitative differences. All data presented below were extracted from simulations with edge length H=5 μm of the simulation volume, filament length L_f=4 μm and persistence length l_p ≈ 9.2 μm, filament discretization length h=0.125 μm, linker size R_l=0.1 μm, and size tolerance ΔR_l=0.02 μm. We set the association rate constant to k_on=90 s⁻¹, and the dissociation rate constant to k_off=3 s⁻¹, which translates to a binding energy of E_b ≈ 3.4k_BT. The temperature was set to T=293 K. All simulations were performed with N_f=208 filaments (c_f ≈ 4 μM) except for those described in §5a (for the assessment of the elastic network strain energy) and c of §5, where we used N_f=104 (c_f ≈ 2 μM). A lower filament concentration facilitates the observation of the effect of bundling on the mean elastic strain energy of the filaments and expands the region in phase space where a direct transition from a homogeneous-isotropic to a lamellar network can occur. Other simulation parameters are provided in the electronic supplementary material. Altogether, network evolution was studied for simulated times t_sim ≥ 1000 s, and morphological analyses covered the interval [800 s, 1000 s] for all networks with N_f=208. For networks with a lower filament concentration, we chose t_sim ≥ 2400 s and [2000 s, 2400 s]. All network structures had by then long reached a steady state, which was evidenced by time-invariant DDCF and structure function.

The variation of essential parameters such as L_f, l_p, R_l, and especially the binding site separation h influence critical linker concentrations, thus shifting phase boundaries. However, tuning these quantities within a physiologically sensible range does not lead to the emergence of any additional morphologies. Even with chiral filaments with a helical binding site architecture as in the case of F-actin, the fundamental morphologies remain the same (see the electronic supplementary material).

A phase transition where the derivatives of the free energy function across the phase boundary up to order (n−1) are continuous while the nth derivative is discontinuous is called an nth-order phase transition. For first-order phase transitions (like melting), one typically observes an abrupt change in the microstructure (and thus, e.g. in the order parameter) while second-order phase transitions (as observed, e.g. around the Curie point, which separates paramagnetic and ferromagnetic materials) go along with smooth changes (e.g. of the order parameter).

5. Microstructure and phase transitions

While the interval of possible binding angles ϕ is determined by the molecular structure of a linker and is therefore constant for networks with a specific linker species, the linker concentration n_l can be modulated easily. For example, in the cytoskeleton, the most relevant phase transitions are therefore those along vertical lines in the phase diagram in figure 3: from homogeneous-isotropic gels to the bundle phase (§5a), to the cluster phase (§5b) and to the lamellar phase (§5c), as well as from the cluster phase to the lamellar phase (§5d). The transition from clusters to lamellae can be subdivided into transitions to squaratic and hexatic lamellae. We studied all these phase transitions by gradually increasing the linker concentration. In each case, we chose some (constant) preferred binding angle ϕ₀ with angular tolerance Δϕ=1/16π (having verified robustness of our results against moderate changes in angular tolerance).

(a) Transition from a homogeneous-isotropic gel to the bundle phase

The first major morphological transition results in a bundle network as shown in figure 3. As illustrated in figure 4, the phase transition from a homogeneous-isotropic gel to a bundle network is accompanied by an abrupt change of the network morphology at a critical linker concentration n_l,crit rather than by a continuous transition. We observe bundling for linkers which allow connections between filaments enclosing small angles up to ϕ=π/4, i.e. for parallel linkers. Specifically, we discuss linkers with a preferred binding angle ϕ₀=1/16π. Within the narrow concentration interval n_l∈[0.058;0.062] (i.e. N_l ∈ [400;425]), the DDCF C₁(d₁) in (normal) x₁-direction in figure 4a suddenly changes from a nearly uniform distribution to a peaked distribution with a sharp maximum around zero (and H because of the periodicity of the system). The width of this peak is comparable to the bundle diameter. At concentrations only slightly above n_l,crit, phase separation is observed: a single bundle is surrounded by free filaments devoid of linkers. A further increase of the linker concentration results in almost all filaments available in the simulated volume being aggregated in the bundle. The temporal evolution of such a bundle involves the formation of several thin bundles, which are merged subsequently (some physiological linkers limit bundle growth [34] but only by means of geometric constraints more sophisticated than considered herein. See the electronic supplementary material for details). An increasing linker concentration entails an increasingly long-lived transient state of several bundles, which are only scarcely cross-linked among each other. The corresponding structure functions I₁(q₁) in x₁-direction in figure 4b can be interpreted as follows. The nearly constant two-point density–density correlation in homogeneous-isotropic networks translates into a Dirac function in the frequency space (→ absence of density fluctuations for finite distances). The characteristic dominance of short linker-to-linker distances in bundles (due to their small radius) translates into a broader spectrum in the frequency space with significant linker density variations.

The orientational distribution for homogeneous-isotropic networks with a normalized linker concentration of n_l ≤ 0.058 (N_l=400) is uniform as can be seen in figure 4c. The perfect uniform distribution of azimuth angles of the filament orientations is represented by a circle with radius $1 / \sqrt{π}$ . For n_l>0.062 (N_l=425), the distribution exhibits pronounced maxima in two opposite directions, which is to be interpreted as a strong uniaxial orientational preference. We do not discuss the axial x₂-direction separately as it does not offer additional insight except for the fact that the linker density throughout the bundle is uniform.

During the phase transition, the mean internal elastic energy 〈E_int〉 of the filaments quickly decreases by approximately 5% (figure 4d). The comparatively smooth decrease around n_l,crit is related to the fact that the bundle itself is subject to thermal fluctuations on the scale of the entire structure [10]. It keeps decreasing with increasing linker concentration since cross-links effectively quench thermal fluctuations of the filaments and straighten them. Linker entropy becomes the dominant factor in the total free energy of the system. Note that we employed a reduced filament concentration of N_f=104 filaments (c_f=2μ), which is why the normalized linker concentration n_l is shifted to larger values. We choose a smaller N_f in order to reduce the shielding effect of phase separation, which occurs for the present transition.

Remarkably, we find that both linkers with little orientational freedom studied above (i.e. ϕ₀=π/16 and Δϕ₀=π/16) and linkers with maximal orientational freedom (i.e. Δϕ₀=π/2) induce bundling. In the first case, bundling occurs due to the linkers’ geometrical constraints, which selectively promote the parallel alignment of filaments. In the second case, bundling occurs because arranging filaments in bundles maximizes the number of potential cross-linking sites, representing the state of maximal linker entropy for this kind of filament–linker system. In the latter case, as found in experiments [7,35], linkers create long-lived non-equilibrium states of several interconnected or intertwined bundles (even asterisk-like structures, see the electronic supplementary material) during the process of thermodynamic equilibration as opposed to the essentially unconnected bundles described above. Yet, the final equilibrated state again is one single bundle in the predominant case where the filament length by far exceeds the linker size.

(b) Transition from a homogeneous-isotropic gel to the cluster phase

Next, we examined the morphological transition from a homogeneous-isotropic network to a cluster network (figure 3). Again, we studied the changes in the DDCF, the structure function and the orientational distribution for increasing linker concentrations and for different parameter settings. Typical results are depicted in figure 5 for preferred binding angles ϕ₀=7π/16 and ϕ₀=5π/16, respectively. These angular preferences represent the cases for which linkers at very high concentrations eventually induce squaratic or hexatic order (cf. §5d). At the intermediate linker concentrations investigated here, however, the effect of the angular preference is only moderate.

Figure 5. Transition from homogeneous-isotropic gel to a cluster network with preferred binding angle ϕ₀=7π/16 (a,b,c,g) and ϕ₀=5π/16 (d,e,g,h): DDCF C₁(d₁) (a,d), structure function I₁(q₁) (b,e), and azimuth angle distribution ρ(φ) of the filament orientations (c,f). Normalized linkerconcentrations for ϕ₀= 7π/16: n_l ≈ 0.044 (dark red squares), 0.102 (pink triangles), 0.106 (purple diamonds), 0.117 (blue circles) (N_l=300,700,725,800 linkers); for ϕ₀=5π/16: n_l ≈ 0.058 (light red squares), 0.113 (red triangles), 0.117 (magenta diamonds), 0.120 (violet circles), 0.124 (inverted purple triangles), 0.131 (blue triangles right) (N_l=400,775,800,825,850,900 linkers). The black circle in (c) with radius $1 / \sqrt{π}$ is the uniform distribution. Inset in (b): local coordinate frame. (*g,h*): mean elastic strain energy of the filaments. (Online version in colour.)

Crossing a critical linker concentration n_l,crit induces an abrupt morphological change of the network. For linkers with ϕ₀=7π/16, the change occurs between n_l=0.102 and n_l=0.106, for linkers with ϕ₀=5π/16, we find n_l=0.113 and n_l=0.117. This change is marked by a peaked density–density correlation instead of a hitherto nearly uniform spatial distribution of the linker molecules (cf. figure 5a,d). The width of this peak is comparable to the cluster diameter suggesting a high linker density in the cluster core and a low one in the corona. The non-Dirac structure functions of cluster networks in figure 5d,e exhibit finite-amplitude frequency contributions up to about the 5th and 10th harmonic suggesting increased spatial fluctuations in linker density. The quick monotonous decay of the structure function suggests a smoothly decreasing linker concentration with increasing distance from the cluster core. While the position of filaments and linkers changes significantly during phase transition, the orientation distribution remains largely unaffected, that is, isotropic (cf. figure 5c,f, similar results found for polar angle).

In general, the phase transition from a homogeneous-isotropic gel into a cluster occurs abruptly and in a very similar way for linkers with different preferred binding angles ϕ₀. For smaller preferred binding angles (e.g. ϕ₀=5π/16), the transition is followed by a rather continuous further evolution, which may be interpreted as the beginning of the subsequent (and for small preferred binding angles quite smooth) phase transition from the cluster phase into the lamella phase (cf. §(d)). Here, we observe a continuous flattening of the cluster.

Figure 6 depicts the OCF and the order parameter of a homogeneous-isotropic gel as well as of two clusters with ϕ₀=7π/16 and ϕ₀=5π/16. As expected for homogeneous-isotropic networks, there is no orientation correlation neither on short distances (as captured by the OCF) nor on long distances (as captured by the order parameter). This can be understood from the low number of cross-links and thus very few linker-mediated interactions between filaments, documenting the importance of such interactions for the establishment of orientational correlation and order. While the cluster with ϕ₀=7π/16 exhibits no perceptible orientation correlation, the cluster with ϕ₀=5π/16 displays the onset of an orientational correlation. As consequence, clusters, or rather networks of clusters at least for larger ϕ₀, can be perceived as an isotropic material. In both clusters, linkers impose a certain orientational order over short distances apparent from the positive S(d)>0 for small distances d between filament-binding sites. The OCF reveals that this orientational order originates from a preferred enclosed angle around ϕ₀ in neighbouring filaments.

Upon transition from a homogeneous-isotropic network to cluster networks, we observed an increase in the mean elastic strain energy of the filaments 〈E_int〉 by approximately 7% for both linker types as shown in figure 5g,h. The effect is caused by linkers, which induce additional filament deformation in order to increase the number of available spots for doubly bound linkers, following a rationale similar to the one explaining the straightening of the filaments in the bundle transition (cf. §5a).

(c) Transition from homogeneous-isotropic gel to the lamella phase

For preferred binding angles slightly larger than π/4, a direct phase transition from a homogeneous-isotropic gel to a lamellar network is possible as illustrated in figure 3. The critical linker concentration n_l,crit is found between n_l=0.291 and n_l=0.299 (for a reduced number of filaments N_f=104).

During the phase transition, the DDCFs in figure 7a,c reveal only a moderate filament compaction in in-plane direction but a significant one orthogonal to it (in x₁-direction). This morphological change is reflected even more clearly by the structure functions in figure 7b,d. Within the geometrical limits of the lamella, the linker distribution remains basically uniform. By contrast, perpendicular to the lamellar phase, the linker concentration rapidly and almost instantaneously decreases to very low values. The lamella exhibits a well-defined sixfold symmetry, which leads to a hexagonal distribution of azimuth angles of the filament orientations in figure 7e. The distribution of polar angles (figure 7f) peaks at ψ=π/2 due to the choice of the local coordinate frame.

Interestingly, our simulations revealed that the number of doubly bound linkers increases only slightly during the phase transition itself. When the total linker concentration is increased even further, however, the number of doubly bound linkers increases dramatically indicating a progressive restructuring from a rather disordered lamella towards one with a strict internal hexagonal symmetry. This restructuring explains the ongoing changes in the structure function in figure 7b even beyond n_l,crit.

(d) Transition from clusters to lamellae

Phase transitions from clusters to lamellae can be observed on a rather broad range of preferred binding angles $π / 4 ≲ ϕ_{0} \leq π / 2$ . We will discuss two exemplary cases, ϕ₀=5π/16 and ϕ₀=7π/16 (angle tolerance Δϕ=π/16), which represent the antipoles of this transition. The former case is summarized in figure 8, the latter in figure 9.

Figure 9. Transition from cluster network to a lamellar network with ϕ₀=7π/16. DDCF C₁(d₁) in x₁-direction (a) and its structure function I₁(q₁) (b). DDCF C₂(d₂) (c) and structurefunction I₂(q₂) in x₂-direction (d). Azimuth angle distribution ρ(φ) (e) and polar angle distribution ρ(ψ) (f) of the filament orientations; mean elastic strain energy of the filaments (g). Normalized linker concentrations: n_l ≈ 0.131 (red squares), 0.157 (pink triangles), 0.160 (violet diamonds) and 0.212 (blue circles) (N_l=900,1075,1100,1500 linkers). (Online version in colour.)

Above a critical linker concentration n_l,crit, lamellae with ϕ₀=5π/16 develop a sixfold, hexagonal symmetry (figure 8e), whereas ϕ₀=7π/16 results in a fourfold, tetragonal symmetry (figure 9e). These symmetries are the free energy minima of the system for these ϕ₀ and are characterized by a minimized filament entropy (i.e. thermal fluctuations) due to a maximized linker entropy (i.e. the maximization of the number of binding site pairs eligible for cross-linking). Although both structures are lamellae, network evolution takes different paths. For ϕ₀=5/16, an increasing linker concentration n_l entails a smooth decay of the non-zero modes of the in-plane structure functions (figure 8d) towards the limit of a Dirac function. The linkers distributed in the cluster core are spread out uniformly in the plane of the lamella. By contrast, choosing ϕ₀=7π/16 induces an abrupt quenching of most of the higher modes of the in-plane structure function above n_l,crit (figure 9d). This becomes even more apparent for the out-of-plane structure function shown in figure 9b. We detect a sudden emergence of higher order modes above n_l,crit, evidence of a sudden flattening of the cluster into a lamellar structure. In the case of linkers with ϕ₀=5π/16, this flattening is much less abrupt. Rather, after an initial abrupt but incomplete flattening, we observe a continuous evolution into a lamellar structure. The thickness of both lamellae approaches values on the scale of the linker length for high n_l. The respective evolution of the distribution of filaments orientations in figures 8 and 9 reinforce the general interpretation of our observations.

In summary, we find strong evidence of different linker-dependent modes of morphological transitions. In case of a preferred binding angle ϕ₀=5π/16, increasing linker concentration drives a smooth network evolution. In §5c, this phenomenon has already been linked to an increasing degree of hexagonal order. For linkers, that preferably connect nearly orthogonal filaments (e.g. ϕ₀=7π/16), increasing linker concentration triggers a discontinuous morphological change. From figure 6a, we learn that the cluster geometry imposes a natural preference for large angles between neighbouring filaments. As a consequence, the transition to lamellae requires much less filament reorientation in case of linkers with ϕ₀ ∼ π/2 than with ϕ₀ ∼ π/3. As illustrated in figure 10a, the preferred linker binding angles induce two peaks in the OCFs O(θ) (cf. equation (3.3)), one close to θ=ϕ₀, the other close to θ=0. The first peak reflects the linker-mediated orientational preference of the filaments, while the second peak simply reflects the presence of a large number of filaments in the same direction due to the four- and sixfold symmetry. The internal symmetry of the lamellae furthermore induces a long-range orientation correlation which is reflected in the large values of the order parameter S(d) even at greater filament-binding site distances d (figure 10b).

Figure 10. (a) OCFs O(θ) and (b) the corresponding order parameters S(d) for a homogeneous-isotropic gel (red squares),a lamellar network with ϕ₀=7π/16 and n_l=0.189 (N_l=1300, pink triangles), and a lamellar network with ϕ₀=5π/16 and n_l=0.189 (blue diamonds). Uniform distribution of filament orientations: $O (θ) = \sin θ$ and S(d)=0 (black lines). (Online version in colour.)

Figures 8g and 9g show the evolution of the mean elastic strain energy of the filaments for the hexatic and the squaratic lamella, respectively. The hexatic lamella’s continuously increasing mean elastic strain energy of the filaments from n_l ≈ 0.13 up to n_l ≈ 0.15 result from the gradual filament reorientation. With increasing linker concentration, the characteristic sixfold symmetry is imposed, leading to the straightening of the filaments and the continuous reduction of the mean elastic strain energy. The squaratic lamella exhibits an abrupt drop of the mean elastic strain energy of the filaments by approximately 7% occurring between n_l=0.0157 and n_l=0.0163. This results from the effortless and immediate imposition of orthogonal filament order beyond n_l,crit.

6. Conclusion

Vastly different morphologies like stress fibres or lamellipodia have been observed in the cytoskeleton of biological cells. Nevertheless, in homogenized constitutive models, these different thermodynamical phases and the transitions between them have so far largely been ignored, also owing to a lack of accurate information.

In this paper, we presented the results of a comprehensive computational study of all relevant thermodynamic equilibrium phases and phase transitions in cross-linked, semi-flexible polymer networks like the cytoskeleton. The phases examined can directly be related to morphologies observed in the cytoskeleton. Bundles correspond to actin stress fibres cross-linked by α-actinin [3]. Lamellae with a preferred binding angle ϕ₀=5π/16 correspond to lamellipodia, in which Arp2/3 complex imposes an inter-filament angle of around 70° [4]. The results of our computational study can therefore directly serve as a basis for improved homogenized constitutive laws for cells. Our quantitative examination of the microstructure of the different phases may motivate appropriate strain energy functions in homogenized models. For example, stress fibres (bundles) may be modelled by transversely isotropic strain energy functions, which, e.g. incorporate the viscoelastic specificities of this kind of structure [10]. By contrast, lamellae (as observed, e.g. in lamellipodia) may rather be modelled as constrained mixtures of two or three fibre families arranged in a planar tetragonal or hexagonal symmetry.

Tuning the linker concentrations in their cytoskeleton, cells may induce thermodynamic phase transitions between phases like homogeneous-isotropic networks and hexagonally symmetric sheets. Based on our computational study, we suggest to capture in homogenized models such phase transitions by switching between different strain energy functions (for the different thermodynamics phases) in an appropriate way. To model the transitions realistically, information about their smoothness is required. Our simulations reveal that all relevant phase transitions happen rather abruptly (suggesting first-order phase transitions), which motivates rather discontinuous switching functions between different strain energy functions on the macro-scale. In case of linkers with a preferred binding angle slightly larger than π/4, a modest discontinuous switch towards a cluster or lamellar phase may be followed by an additional smooth evolution of constitutive parameters. For an increasing linker concentration, this smooth evolution may capture the ongoing restructuring within these phases as detected in our simulations by means of the DDCF and structure function. Noting that linkers do not only drive thermodynamic phase transitions and but also entail dramatic mechanical consequences such as passive [36] and active stiffening [37], we suggest the incorporation of the cytoskeleton’s linker concentrations into homogenized constitutive models as thermodynamically motivated internal variables modulating time-dependent changes in the constitutive properties. Note, however, that one has to carefully consider possible problems related to non-convex energy functions as well known from other materials subject to phase transitions (cf. e.g. [38]). As demonstrated above, our Brownian dynamics approach provides a powerful predictive tool capable of driving the formulation of such complex energy functions.

Several of the most important cellular processes like mechanosensing and cell migration are intimately related to vast morphological changes (i.e. phase transitions) in the cytoskeleton that are mostly neglected by current models. There is a rapidly growing interest in modelling cellular mechanics in order to understand these processes better and to potentially alter them. Therefore, we see a pressing need for more realistic, computationally efficient constitutive models of the cytoskeleton. This paper is intended to provide essential micromechanical and thermodynamical information required for the development of such models.

Data accessibility

Simulations were conducted using the institute’s multiphysics code. All presented data are reproducible employing the results of our cited publications.

Authors' contributions

All authors contributed to the design and the conception of the simulations. K.W.M. and C.J.C. devised the computational method and conducted the simulations. All authors participated in writing the paper.

Competing interests

We have no competing interests.

Funding

The research of the first author was funded by the International Graduate School of Science and Engineering (IGSSE) of Technische Universität München (TUM) and the Emmy Noether programme of the German Research Foundation DFG (CY 75/2-1).

Acknowledgements

We thank Robijn Bruinsma (UCLA), Andreas Bausch and Oliver Lieleg (TUM), as well as Andreas Rauch for fruitful discussions.

Footnotes

VIEW FULL TEXT DOWNLOAD PDF

This Issue

Proceedings of the Royal Society A: Mathematical, Physical and Engineering Sciences cover image

October 2015

Volume 471Issue 2182

Article Information

DOI:https://doi.org/10.1098/rspa.2015.0332
Published by:Royal Society
Online ISSN:1471-2946

History:

Manuscript received21/05/2015
Manuscript accepted17/09/2015
Published online08/10/2015
Published in print08/10/2015

License:

Citations and impact

Keywords

Subjects

Your password must have 8 characters or more and contain 3 of the following:

Computational analysis of morphologies and phase transitions of cross-linked, semi-flexible polymer networks

Abstract