## Abstract

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

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

**and**

*x***from the balance of linear and angular momentum**

*θ**a*and

*b*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

*a*and

*b*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

*a*and

*b*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 2*a*, 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 2*b*, thus becoming singly bound. When binding to a second filament, a singly bound linker becomes doubly bound, i.e. a cross-link (figure 2*c*). Free linkers are treated as rigid particles whose centre of mass *x*_{l} diffuses according to the equation of motion
*ζ*_{l}=6*πηR*_{l}. The stochastic force exerted on a linker is given by
** 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.

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 *ϕ*∈[*ϕ*_{0}−Δ*ϕ*;*ϕ*_{0}+Δ*ϕ*]. 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
*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*_{1},*x*_{2},*x*_{3}) 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*_{1}-direction to be orthogonal to the bundle axis. Finally, for the lamellae, the *x*_{1}-direction is chosen orthogonal to the plane in which they lie. The respective coordinate frames are shown as insets of figures 4*b*, 5*b* and 7*b*.

### (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
*L*_{n} denotes the length of the *n*th filament, 〈⋅〉 the time average, *δ* the Dirac-delta function and *φ*_{n}(*s*) and *ψ*_{n}(*s*) the orientation angles of the tangent of the *n*th filament at the material point *s*. The azimuth angle *φ*∈[0;2*π*[ is measured against the *x*_{2}-direction, the polar angle *ψ*∈[0;*π*] against the *x*_{1}-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 *m*th binding site on the *i*th filament be *x*^{[i][m]}_{b}. Then,
*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)
*π*/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 *π*/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
*P*_{2} 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
*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*^{k} and *x*^{k′} and two non-identical variable points in space ** x** and

*H*].

### (f) Structure function

The Fourier transform of the two-point DDCF
** 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

### (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 *ϕ*_{0} assuming Δ*ϕ*=*π*/16 (cf. §2 on linkers), and the normalized linker concentration

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 *ϕ*_{0}, a *cluster* phase at intermediate linker concentrations and large preferred binding angles *ϕ*_{0}, and a *lamellar* phase with orthogonal filament order (→squaratic) for *ϕ*_{0} close to *π*/2 and hexagonal filament order (→hexatic) for *ϕ*_{0} 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*_{f}*A*_{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^{−1}, and the dissociation rate constant to *k*_{off}=3 s^{−1}, which translates to a binding energy of *E*_{b} ≈ 3.4*k*_{B}*T*. 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 *n*th derivative is discontinuous is called an *n*th-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 *ϕ*_{0} 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 *ϕ*_{0}=1/16*π*. Within the narrow concentration interval *n*_{l}∈[0.058;0.062] (i.e. *N*_{l} ∈ [400;425]), the DDCF *C*_{1}(*d*_{1}) in (normal) *x*_{1}-direction in figure 4*a* 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*_{1}(*q*_{1}) in *x*_{1}-direction in figure 4*b* 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 4*c*. The perfect uniform distribution of azimuth angles of the filament orientations is represented by a circle with radius *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*_{2}-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 4*d*). 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. *ϕ*_{0}=*π*/16 and Δ*ϕ*_{0}=*π*/16) and linkers with maximal orientational freedom (i.e. Δ*ϕ*_{0}=*π*/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 *ϕ*_{0}=7*π*/16 and *ϕ*_{0}=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.

Crossing a critical linker concentration *n*_{l,crit} induces an abrupt morphological change of the network. For linkers with *ϕ*_{0}=7*π*/16, the change occurs between *n*_{l}=0.102 and *n*_{l}=0.106, for linkers with *ϕ*_{0}=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 5*a*,*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 5*d*,*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 5*c*,*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 *ϕ*_{0}. For smaller preferred binding angles (e.g. *ϕ*_{0}=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 *ϕ*_{0}=7*π*/16 and *ϕ*_{0}=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 *ϕ*_{0}=7*π*/16 exhibits no perceptible orientation correlation, the cluster with *ϕ*_{0}=5*π*/16 displays the onset of an orientational correlation. As consequence, clusters, or rather networks of clusters at least for larger *ϕ*_{0}, 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 *ϕ*_{0} 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 5*g*,*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 7*a,c* reveal only a moderate filament compaction in in-plane direction but a significant one orthogonal to it (in *x*_{1}-direction). This morphological change is reflected even more clearly by the structure functions in figure 7*b,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 7*e*. The distribution of polar angles (figure 7*f*) 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 7*b* 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 *ϕ*_{0}=5*π*/16 and *ϕ*_{0}=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.

Above a critical linker concentration *n*_{l,crit}, lamellae with *ϕ*_{0}=5*π*/16 develop a sixfold, hexagonal symmetry (figure 8*e*), whereas *ϕ*_{0}=7*π*/16 results in a fourfold, tetragonal symmetry (figure 9*e*). These symmetries are the free energy minima of the system for these *ϕ*_{0} 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 *ϕ*_{0}=5/16, an increasing linker concentration *n*_{l} entails a smooth decay of the non-zero modes of the in-plane structure functions (figure 8*d*) 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 *ϕ*_{0}=7*π*/16 induces an abrupt quenching of most of the higher modes of the in-plane structure function above *n*_{l,crit} (figure 9*d*). This becomes even more apparent for the out-of-plane structure function shown in figure 9*b*. 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 *ϕ*_{0}=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 *ϕ*_{0}=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. *ϕ*_{0}=7*π*/16), increasing linker concentration triggers a discontinuous morphological change. From figure 6*a*, 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 *ϕ*_{0} ∼ *π*/2 than with *ϕ*_{0} ∼ *π*/3. As illustrated in figure 10*a*, the preferred linker binding angles induce two peaks in the OCFs *O*(*θ*) (cf. equation (3.3)), one close to *θ*=*ϕ*_{0}, 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 10*b*).

Figures 8*g* and 9*g* 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 *ϕ*_{0}=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.

- Received May 21, 2015.
- Accepted September 17, 2015.

- © 2015 The Author(s)