Theory of Weakly Polydisperse Cytoskeleton Filaments.

Cytoskeleton filaments have the extraordinary ability to change conformations dynamically in response to alterations of the number density of actins/tubulin, the number density and type of binding agents, and the electrolyte concentration. This property is crucial for eukaryotic cells to achieve specific biological functions in different cellular compartments. Conventional approaches to biopolymers' solution break down for cytoskeleton filaments because they entail several approximations to treat their polyelectrolyte and mechanical properties. In this article, we introduce a novel density functional theory for polydisperse, semiflexible cytoskeleton filaments. The approach accounts for the equilibrium polymerization kinetics, length and orientation filament distributions, as well as the electrostatic interaction between filaments and the electrolyte. This is essential for cytoskeleton polymerization in different cell compartments generating filaments of different lengths, sometimes long enough to become semiflexible. We characterized the thermodynamics properties of actin filaments in electrolyte aqueous solutions. We calculated the free energy, pressure, chemical potential, and second virial coefficient for each filament conformation. We also calculated the phase diagram of actin filaments' solution and compared with the corresponding results in in vitro experiments.

1. Introduction

Eukaryotic cells can dynamically regulate the biological environment and the polyelectrolyte and mechanical properties of cytoskeleton filaments to achieve specific biological functions as diverse as directional growth, shape, division, plasticity, and migration [1]. For instance, an increase in the number density of G-actin/tubulin and electrolyte concentration can lead to conformation transformations from the orientation-disordered (isotropic) to orientation-ordered (nematic) phase, as well as increasing the filaments’ average length. Additionally, a growth in the number density of binding agents, such as divalent ions or linker proteins, can yield bundling or network conformations [2]. These self-organization behaviors, yet poorly understood, have been observed experimentally. Currently, valuable information on the distribution and type of cytoskeleton conformations in cells is obtained from fluorescence and electron microscopy images [3,4,5,6], whereas confocal microscopy captures their dynamic conformation changes [7,8,9]. However, this information usually provides an incomplete molecular understanding of the interplay between the polydispersity, semiflexibility, polyelectrolyte, and mechanical properties of cytoskeleton filaments on their conformational dynamics, self-organization, and stability. This understanding is crucial to elucidate the biophysical principles underlying fundamental biological functions of eukaryotic cells in normal and pathological conditions, which may vary depending on the cell type and location, gender, age, and inheritance conditions.

A substantial amount of theoretical research has been performed in the field to study the isotropic to nematic phase transformation in macromolecules’ solution. The conventional understanding of the properties of these polyelectrolytes is based on monodisperse (e.g., same filament lengths), mean-field theories, and rod-like cylindrical filament models (e.g., with contour lengths shorter than their persistence length). These methods break down for cytoskeleton filaments because they entail several approximations to treat the inter-filament interactions, electrolytes, and filament structures. Onsager [10] showed that the expansion of the free energy functional up to the second virial coefficient provides accurate results for monodisperse long charged rods in the coexisting isotropic and nematic phases. In the nematic phase, the free energy is expressed as a functional of angular distribution functions of rods, which minimizes the free energy functional. Onsager introduced a trial function for the angular distribution functions depending on a single parameter to increase the efficiency of these calculations. Subsequently, a variety of modifications and extensions of Onsager’s theory were proposed [11,12]. For instance, Odijk introduced a Gaussian-type trial function for the angular distribution. While less accurate than Onsager’s approximation, it overcame some limitations on the numerical calculations. It provides explicit, analytic expressions for the thermodynamic properties that indeed facilitated the phase diagram transition and coexistence analysis [13]. A different approach was proposed to estimate the angular distribution functions for monodisperse, long cytoskeleton filament rods [2]. They were not obtained by minimizing the free energy, but postulating a particular form in each phase with a single variational parameter characterizing the width of the angular distributions.

Additionally, corrections were made to the orientational part of Onsager’s free energy to consider monodisperse macromolecules having contour lengths larger than or the same order as their persistence length [14,15]. Furthermore, later approximations were developed to patch the orientational free energies for rigid and semiflexible macromolecules [13,16]. A different approach was introduced by Sluckin [17], who generalized Onsager’s theory to describe polydisperse rigid rods having the same Gaussian form for the size distribution function in both isotropic and nematic phases. Moreover, Odijk’s ansatz was generalized to the polydisperse case, where the nematic phase onset was formed by rods with lengths larger than the average length of the size distribution in the isotropic phase. Different modifications of the method were proposed to address this shortcoming [18,19,20,21]. In particular, the input size asymmetric distribution functions in the form of Schultz’s and log-normal distributions were considered.

On the other hand, a particular approach was introduced for amphiphilic micellar suspensions [22]. The size distribution function of micellar macromolecules was not considered as an input value, but rather the standard chemical potential that governs the size–angular distribution functions [23]. It was found that the micella average lengths in the coexisting isotropic and nematic phases are different, with larger micella sizes in the nematic phase.

The size distribution function of uncharged polydisperse actin filament rods in the bundling phase was calculated from a free energy that accounts for hard-core filament repulsion and short-range attractions [24]. It was shown that short-range attractions, arising either from linker proteins, depletion-mediated attractions, or polyvalent ions, enhance the tendency of filaments to align parallel to each other, yielding an increase in the average filament length and a decrease in the relative width of the distribution of filament lengths.

These approaches for phase diagram studies consider specific macromolecular properties, neglecting others. For instance, some approaches focused on polyelectrolyte and polydispersity properties only, whereas other approximations accounted for semiflexibility and polyelectrolyte properties, and so on. However, for cytoskeleton filaments, it is imperative to consider the balance and competition between contributions coming from their polydispersity, semiflexibility, polyelectrolyte, and mechanical properties to the total free energy. When accounting for all of these features, one can formulate a more accurate and realistic description of the conformational dynamics, self-organization, and stability properties of cytoskeleton filaments in different cell compartments.

As a first step to face this challenge, we introduce in this article a novel density functional theory for polydisperse, semiflexible cytoskeleton filaments. The approach accounts for the equilibrium polymerization kinetics, length and orientation filament distributions, as well as the electrostatic interaction between filaments and the electrolyte. As a unique feature, the formulation is able to determine critical parameter values governing the isotropic–nematic phase diagram behavior. This approach is essential to study the self-organization behavior of actin filaments taking place in different cell compartments, where the G-actin polymerization and electrolyte conditions may generate filaments with different conformations and lengths long enough to become semiflexible. Specifically, we consider experimental conditions on the actin filaments’ solution where the isotropic–nematic phase diagram transition is due to changes in the G-actin concentration. Additionally, the filament average size for different G-actin concentrations was fixed by changing the gelsolin proteins’ concentration [25,26]. From these special equilibrium conditions, we obtained the size distribution function, whereas we used Sluckin’s trial function to calculate the angular distribution function in the nematic phase. Additionally, we introduced an ansatz for the standard chemical potential excess of actin filaments that results in the asymmetric Schulz distribution function for the actin filaments size. This distribution function agrees with those used in light-scattering experiments on actin filaments [25]. In addition, we generalized the formula for the orientational free energy introduced in [13] to the case of a polydisperse system to account for the filament semiflexibility. Finally, we calculated the isotropic–nematic phase diagrams for a variety of Schulz size distributions, persistence lengths, and concentrations of monovalent ions in the electrolyte solution. We also compared these results with available experimental data [26].

The paper is organized as follows. The theory is described in Section 2; the numerical results are given in Section 3; the discussion is provided in Section 4; the details of the calculations are presented in Appendix A, Appendix B, Appendix C, Appendix D and Appendix E.

2. Materials and Methods

We considered a solution of polydisperse actin filaments, each of them having a length , where is the length of an actin monomer unit and the number of monomer units representing the filament size. Each filament has its own direction in 3D space, which is characterized by a body-angle accounting for the chosen nematic director. In spherical coordinates, we have with the z direction taken along the nematic axis.

The size–angular density distribution function of actin filaments is given by the following expression: where is the total density of the filaments, N is the total filament number of any size, and V represents the volume of the system. Basically, represents the number of filaments of size oriented along the direction . Additionally, is the size distribution function averaged over all angles, and represents the angular distribution functions, which also depend on the filament size . The normalization conditions for the size- and angular-distribution functions and are and respectively. The summation in Equation (2) is performed over all filament sizes , whereas the integration in Equation (3) is carried out over the whole body angle. The size distribution function is characterized by the filament average size , namely the average degree of polymerization, as well as the normalized standard deviation , such as and where

2.1. The Free Energy Density Functional

We define the dimensionless Helmholtz free energy of cytoskeleton filaments f as follows: where F is the free energy, the inverse thermal energy, the Boltzmann constant, T the temperature, the dimensionless ideal gas free energy, the dimensionless energy due to the inter-filament interactions, and the standard chemical potential of -sized filaments.

The expressions for and as functionals of the density distributions and are provided below.

2.1.1. Ideal Gas Free Energy

The ideal gas free energy per volume of the system V can be written in the following form: where is the thermal de Broglie wavelength of -sized filaments, h is the Planck constant, and stands for the -sized filament mass. Substitution of Equation (1) into Equation (8) yields where is the ideal gas free energy per filament (), the monomer thermal de Broglie wavelength, and the filament orientational entropy:

2.1.2. Interaction Free Energy

We write the interaction free energy per volume V in a mean-field fashion, namely where represents the cluster integral

In Equation (12), is the interaction potential between two filaments and the closest distance between them. The expression (11) is a good approximation for the interaction free energy of monodisperse rigid charged rods. In fact, it becomes exact in the limit of infinitely long rods, i.e., for [10].

Substitution of Equation (1) into Equation (11) yields the following expression for as a functional of the densities and :

To calculate the cluster integral in obvious form, we model the interaction potential between two charged rod filaments as follows: where D is the bare rod diameter and the angle between them. The function represents the following repulsive electrostatic inter-filament interaction potential: [27] where

In Equation (16), stands for the inverse Debye length, the modified Bessel function of second kind of first order, the filament linear charge density, the water solvent dielectric permittivity, e the electron charge, and and the ion valency and concentration of species i in solution, respectively. We note that Equation (15) comes from the solution of the linearized Poisson–Boltzmann equation, which is accurate for monovalent ions, namely . Substitution of Equations (14) and (15) into Equation (12) results in the following expression for the cluster integral where and are the effective diameter and the dimensionless function given by Equations (A12) and (A14), respectively. The details of these calculations are presented in Appendix A.

Furthermore, substitution of Equation (17) into Equation (13) provides the following form for the interaction free energy: where is the extension of the second virial coefficient for polydisperse filaments, and the function is given by

We note that the parameter accounts for the orientational averaged contributions coming from the excluded volume and charge density interactions between filaments of size and , as well as the influence of the ionic strength on the second virial coefficient.

It is worth mentioning that the solution of the linearized Poisson–Boltzmann equation in Equations (14)–(16) is valid for infinitely long filaments. Indeed, some correction terms should be considered for short filaments to account for the end effects. However, we considered filaments of average size units of monomers , which correspond to filament lengths in the range of μm only. Thus, the linearized Poisson–Boltzmann equation solution is justified in our study.

2.2. The Distribution Functions

We introduce the Lagrange functional: to find the distribution functions, where f represents the free energy functional given by Equation (7), whereas A and are the Lagrange multipliers accounting for the normalization conditions given by Equations (2) and (3), respectively. For a fixed angular distribution function , the equilibrium size distribution function minimizes the Lagrange functional. Thus, we have

Using the chain rule, the l.h.s. of Equation (21) can be rewritten in the following form: where

An additional relationship for the distribution functions comes from Equation (4). Differentiation on both sides of Equation (4) yields

As a result, Equation (22) can be written as follows:

Substitution of Equation (24) into the r.h.s. of Equation (25) leads to

Consequently, Equations (21) and (26) yield

It is worth noting that -sized filaments can be considered as separated “pseudo-phases”. Additionally, the “pseudo-phases” with all possible sizes are in equilibrium with each other if the equilibrium condition given by Equation (27) is executed. In fact, the value for in Equation (27) does not represent the real chemical potential of -sized filaments since it comes from the Lagrange functional in Equation (23), rather than the Helmholtz free energy f. To calculate , we substitute Equation (20) into Equation (23) and use Equations (3), (7), (9), and (18). We obtain where C is the following constant:

We substitute the expression (28) into the condition of equilibrium (27) to obtain the following equation for the length distribution function : where denotes the standard chemical potential difference between actin units aggregated in a -mer filament and the one in the single dispersed phase , i.e.,

Finally, we substitute the parameter into Equation (30) to obtain the following master equation for :

Similarly, for a fixed distribution function the angular distribution function is obtained by using the variational principle:

Substitution of Equation (20) into Equation (33) and the use of Equations (2), (4), (7), (9), (10), (18), and (19) yield where is a constant that satisfies the following relationship:

To obtain the expression for , we integrate both sides of Equation (34) with respect to and use the normalization conditions (2) and (3). Finally, we replace the expression obtained for into Equation (34) to obtain the following master equation for :

Furthermore, the angular distribution function depends only on the polar angle due to the symmetry of the system. Thus, we have .

To find the size and angular distribution functions and , we solve the expressions (32) and (36) using the successive iterations method. To this end, we chose the monodisperse size distribution as the initial guess for the first iteration step, i.e., , where is a Dirac delta function. Substitution of this expression for into Equation (36) generates the following equation to numerically calculate [28]:

More efficient, but less accurate values for the monodisperse angular distribution function were proposed using some functional forms. For instance, Onsager [10] introduced the trial function . Another commonly used, although slightly less-accurate ansatz is the so-called Odijk’s trial function [13]:

We note that these trial functions depend on a single unknown parameter , which is chosen to minimize the free energy functional f.

In the second iterative step, we substitute the angular distribution function obtained in the first step into Equations (19) and (32)to calculate . Since does not depend on , it follows from Equation (19) that . Additionally, we use in Equation (32) the ansatz: to obtain a size distribution function in the asymmetric Schulz–Zimm form [29]: where z represents the conventional polydispersity parameter, whereas is the constant coming from the normalization condition given by Equation (2). In Appendix B, we show that the polydispersity parameter z is related to the normalized standard deviation of the size distribution function , whereas the parameter y depends on the values and z only, i.e., .

Finally, substitution of Equation (40) into Equation (36) gives an equation to calculate the function . Obtaining the numerical solution of this equation indeed requires a high computational cost because it depends on two variables and . Alternatively, we generalize the monodisperse Odijk’s trial function given by Equation (38) to the polydisperse case. Specifically, we introduce the following parametrization for the polydisperse angular distribution function [17]: where the polydisperse Gaussian parameter depends on the filament size . For the weak polydispersity of the system, the ratio can be considered a small parameter. As a result, the polydisperse Gaussian parameter can be calculated using the following linear expansion around the monodisperse solution: where the two unknown parameters and are chosen to minimize the free energy functional f. In the monodisperse limit, we have that and .

2.3. Free Energy in the Nematic Phase

We substitute the normalization conditions (2) and (3) into Equations (7) and (9) to write the Helmholtz free energy per filament as follows: where and the terms , , represent the orientational, interaction, and mixing free energies, respectively. The expressions for these energies are provided below.

2.3.1. Orientational Free Energy

The orientational free energy in Equation (43) is given by the expression:

Substitution of the expressions for (41) and (42) into Equation (44) yields

Equation (45) represents the generalization of the orientational free energy expression obtained by Odijk for monodisperse rods, i.e., for filaments with contour length , where P is the filament persistence length. Indeed, in the monodisperse limit , the correct limit is obtained for any value of the parameter [13]. For flexible filaments, where , the orientational free energy can be written in the following form [13]:

Substitution of the expressions for (41) and (42) into Equation (46) leads to where and . Certainly, in the monodisperse limit, the correct limit is obtained [11].

To calculate the orientational free energy for any persistence length, we combine these two asymptotic cases for using the following interpolating formula:

The details on these calculations are provided in Appendix C.

2.3.2. Interaction Free Energy

We substitute Equations (40)–(42) into Equations (18) and (19) to obtain where represents the second virial coefficient for polydisperse filaments in the nematic phase and and

The details on these calculations are presented in Appendix D.

Equation (50) represents the generalization of the interaction free energy expression for polydisperse actin filaments. In fact, in the monodisperse limit , the expression for goes to which is the monodisperse distribution function given by Equation (38). Similarly, Equations (50) and (52)–(54) recover the correct expression for the filaments’ interaction free energy in the monodisperse case: with

2.3.3. Mixing Free Energy

The mixing free energy in Equation (43) can be written as:

This expression is not well defined in the monodisperse limit () since the density distribution and the mixing free energy diverges, rather than going to zero. To overcome this shortcoming, we extract the density -independent term in Equation (A45) from the r.h.s. of Equation (58). The resulting renormalized mixing free energy reads

This new expression (59) recovers the correct value in the monodisperse limit, i.e., for (). More details on these calculations are provided in Appendix E.

2.4. Free Energy in the Isotropic Phase

In the isotropic phase, there is no preferential direction for the filament orientations; thus, the orientational distribution function for any filament length, and the orientational free energy becomes

The expression for the interaction free energy is obtained from Equations (4), (18), and (19). We have where is the second virial coefficient for polydisperse filaments in the isotropic phase with

Finally, substitution of Equation (A14) into Equation (63) generates the following expression for h in the isotropic phase where is the elliptic integral given by Equation (A9).

2.5. Isotropic–Nematic Phase at Equilibrium

The isotropic–nematic phase equilibrium in a system of polydisperse macromolecules occurs when the chemical potentials of each species and the osmotic pressures are equal to each other. Specifically, we have [17,19,20]

In the case of the self-aggregation of actin filaments, there is an additional condition for the chemical equilibrium in each thermodynamic phase, namely where is the chemical potential of -sized filaments and the chemical potential of actin monomers. The condition (67) is obtained using a method similar to the one leading to Equation (27). Substitution of Equation (67) into Equation (65) provides a single equilibrium condition for the chemical potentials, i.e.,

To calculate the chemical potential of actin monomers , we write the Gibbs free energy of mixture of actin filaments of all sizes in the following form:

Substitution of Equations (4) and (67) and the use of condition in Equation (69) yield

In Equations (66) and (70), we use the formula for practical calculations of the osmotic pressure.

Using Equation (4), we can also write the average filaments’ length in Equation (70) as follows: where is the G-actin number density. As a result, the average filaments length represents the ratio of G-actin and F-actin concentrations. Additionally, G-actin is usually polymerized in in vitro experiments in the presence of gelsolin, an actin-binding protein known to cap and sever actin filaments [25,26,30]. It was noticed in these experiments that the concentration of F-actin in the solution becomes almost equal to the concentration of gelsolin, i.e., . As a result, the average length of filaments (average degree of polymerization) in these experiments was regulated by the gelsolin concentration. In [26], the variation of the concentration of G-actin was accompanied by a change in the concentration of gelsolin to keep the average length of filaments fixed during the isotropic–nematic phase transition. Using this experimental condition in our theory, the mixing free energy and the term with in the r.h.s. of the expression for free energy (43) are also fixed values and, thus, do not contribute to the phase coexistence properties.

Thus, using Equations (43), (50), and (61), the dimensionless free energy f in the nematic and isotropic phases can be written as follows: where c represents the dimensionless filaments’ density:

The expressions to calculate are given by Equations (49) and (60) for the nematic and isotropic phases, respectively, whereas the corresponding expressions for H are given by Equations (52) and (64).

In the nematic phase, the free energy f (72) is a function of the parameters of the angular distribution function , , and the dimensionless filaments’ density c; thus, . By minimizing f, we obtain the equilibrium parameter functions and for each value of c. Thus, the expression for the nematic free energy at equilibrium becomes . The minimization is carried out using the downhill simplex method in multiple dimensions [31]. Finally, the coexisting dimensionless densities and are obtained from the coexistence conditions (66) and (68), whereas the corresponding actin densities and are obtained from Equations (71) and (73).

In principle, the chemical potential of the non-depleted unimers (G-actin) and the depleted actin concentration are related by Equation (70). We calculate the Helmholtz free energy in Equation (72) as a function of in the isotropic or nematic phase. Subsequently, we substitute the result into the r.h.s. of Equation (70) to obtain the relationship or, vice versa, . In the last case, the chemical potential of unimers can be regarded as an input parameter, which is set by the reservoir condition.

3. Results

In this section, we apply the approach to investigate the isotropic–nematic phase diagram for polydisperse actin filaments in monovalent salt solutions. We considered typical experimental values for key parameters, such as the filament persistence length, the average filament size, the polydispersity parameter, and the electrolyte concentration to elucidate their impact on the orientational and size distributions. We also calculated the thermodynamic properties of the system, including the free energy, pressure, chemical potential and, consequently, the range of G-actin concentrations and filaments’ average lengths leading to conformation transformations from the orientation-disordered (isotropic) to orientation-ordered (nematic) phase. In the numerical calculations, we chose the values for the filaments’ diameter , the monomer units’ length , the actin molar weight kDa, and the linear filament charge density .

We calculated the isotropic and nematic coexisting dimensionless densities for rigid monodisperse rods with length μm for the electrolyte’s concentration M to compare our results with those obtained by Borukhov et al. [2]. Our result was and , whereas the corresponding one in [2] was and (see Figure 3 in [2]). The source of such discrepancy is associated with the cone approximation used in [2] for the angular distribution function.

3.1. Distribution Functions

In this section, we present the analysis on the distribution functions and in the I–N phase coexistence for polydisperse, semiflexible filaments. We used a persistence length μm and an electrolyte concentration M. In Figure 1, we plot the Schulz length distribution function for different average sizes in a weakly polydisperse system with normalized standard deviation ().

Figure 1

Schulz length distribution function as a function of the filament size for different average sizes from Equation (40). The normalized standard deviation is ().

Our results revealed polydisperse distribution functions with lower peaks and broadened distributions for larger values of the average size , while the asymmetric behavior characterizes the different increasing rate lengths of barbed and pointed ends. This asymmetry property plays a crucial role in the phase diagram behavior. Since the normalization integral can be split into the sum of two integrals and (), each of them is independent of the value , and the amount of the filaments with the sizes shorter than is larger than the one with the sizes longer than . In contrast, monodisperse, symmetric distribution functions are represented by delta function distributions .

In Figure 2 and Figure 3, we depicted the I–N phase coexistence values of the parameters and appearing in the angular distribution function for different values of the average length .

Figure 2

The parameter as a function of the average size at the I–N phase coexistence. The normalized standard deviations are (dashed line) and 0.5 (solid line).

Figure 3

The parameter as a function of the average size at the I–N phase equilibrium. The normalized standard deviation is .

It is seen that a larger average length leads to lower parameter values and . In fact, larger filament sizes generate narrower angular distributions, increasing the order of the system. For a given value of , Figure 2 also shows that the parameter decreases with increasing the polydispersity parameter . This is due to the typical broadened size distributions that characterize polydisperse systems, which include short filaments, lowering the order of the system.

In Figure 4, we plot the Gaussian parameter given in Equation (42) as a function of the average filament size . For a fixed size , the parameter decreases with increasing the average length . Additionally, the dependence of the orientational distribution function on the angle for different sizes is plotted in Figure 5.

Figure 4

The Gaussian parameter vs. the average size at I–N phase coexistence for different lengths of filaments .

Figure 5

The angular distribution function as a function of the angle at the I–N phase coexistence for different sizes . The average length is .

It is seen that the increase of size flattens and widens the distribution . On the other hand, the dependence of the total size–angular distribution function on the filament size for several angles is displayed in Figure 6.

Figure 6

The size–angular distribution function as a function of the size at the I–N phase coexistence for different angles . The average length is .

Our results showed that the increase of the angle flattens the total distribution . This is because most filaments in the nematic phase are oriented along the nematic director with . By increasing the angle , the amount of filaments decreases and is directed along the direction with angle .

3.2. Isotropic–Nematic Phase Diagram

In Figure 7, we plot the I–N phase diagram of actin filaments in terms of the coexisting dimensional densities of actin for different average lengths and two values for the normalized standard deviation of Schulz distribution function: (monodisperse) and (weakly polydisperse). In these calculations, we used the values μm and M. In the same figure, we plot the experimental data extracted from Figures 4 and 5 in [26].

Figure 7

The isotropic–nematic phase diagram in terms of coexisting densities of actin for different average lengths of filaments and two values of the normalized standard deviations (dashed lines) and (solid lines). The persistence length is μm, and the electrolyte concentration M. The experimental data (triangles) were retrieved from [26]. The coexisting isotropic and nematic densities are plotted in black and red colors, respectively.

Our results showed isotropic and nematic metastable phases represented by the area enclosed between the red and black curves, whereas the area to the left of the black curves and to the right of the red curves defines the exclusive randomly oriented (isotropic) and parallel oriented (nematic) phases, respectively. Additionally, the red and black curves display an increase in the average filament length with decreasing G-actin concentration and, consequently, the order of the system. For a given G-actin concentration, the nematic phase generates larger average filament lengths as compared to the isotropic phase, whereas higher G-actin concentrations are required in the nematic phase to generate the same average filament length obtained in the isotropic phase. Figure 7 also shows a minor dependence of the coexisting densities on the polydispersity parameter. This is in part due to the experimental condition fixing the same degree of polymerization for both the isotropic and nematic phases. In fact, Figure 7 shows that polydispersity with a larger amount of short filaments slightly reduces both isotropic and nematic coexisting densities as compared to the monodisperse case. Thus, additional factors such as electrostatic effects compete with hard-body interactions to reduce the coexisting densities in the polydisperse case.

We also studied the polydispersity effects on the nematic (anisotropic) order parameter S, which is defined as follows: where is the Legendre polynomial of second order. In the isotropic phase, where is constant, the order parameter S vanishes. In contrast, in a system highly aligned, becomes the Dirac delta function leading to . The I–N phase coexistence values of the nematic order parameter S for different average sizes and two values of the polydispersity parameter are displayed in Figure 8. For both values of , our results revealed a decrease of the order parameter S with increasing the average length . Moreover, the order parameter S is even smaller for polydisperse filaments.

Figure 8

The nematic order parameter S at the I–N phase coexistence for different average lengths . The normalized standard deviation are (dashed line) and 0.5 (solid line).

Overall, we realized that the orientational order of filaments at I–N phase equilibrium is due to the competition between two contributions. In the previous section, we showed that the increase of the average length leads to narrower angular distributions and, thus, higher order in the filaments orientation. On the other hand, Figure 7 displays a decrease in the coexisting nematic density, lowering their order. Since the order parameter S accounts for both contributions and the total order becomes smaller for larger average lengths , we conclude that the contribution coming from the decrease in the coexisting nematic density dominates the order behavior of the filaments over those arising from narrowed angular distributions.

We also considered several filament semiflexibility parameter values since they vary depending on the polymerization buffers, experimental protocols, and techniques. We analyze in Figure 9 the isotropic–nematic phase diagram for two typical values of the persistence length, namely μm and μm, in physiological conditions (electrolyte concentration M).

Figure 9

The isotropic–nematic phase diagram in terms of coexisting densities of actin for different average lengths of filaments and two values of persistence lengths μm (dashed lines) and μm (solid lines). The normalized standard deviation is , and the electrolyte concentration M. The experimental data (triangles) were retrieved from [26]. The coexisting isotropic and nematic densities are plotted in black and red colors, respectively.

For a given value , we found that larger values of the persistence length P decrease both the isotropic and nematic coexisting densities. Additionally, the results for μm give the best fit against the experimental data. This result stems from the impact of the semiflexibility on the orientational free energy . It can be seen from Equation (49) that the increase of the persistence length P leads to the decrease of the orientational term , which boosts the formation of the nematic state and generates smaller values for the coexisting densities.

Finally, we studied the effects of the electrolyte’s concentration on the I–N phase diagram of actin filaments. We plot in Figure 10 the dependence of the effective diameter on the concentration of monovalent ions .

Figure 10

The effective diameter as a function of the monovalent ions’ concentration .

Our results show an exponential decay of with increasing . This result agrees with the formation of an electrical double-layer accumulating many counterions around the filament surface. The lower the electrolyte concentration, the larger its thickness and, consequently, the effective filament diameter are. For salt concentrations larger than M, the effective diameter of filaments is smaller than the diameter of filaments , and the calculations may become inaccurate.

Furthermore, we show in Figure 11 the isotropic–nematic phase diagram for two values of electrolyte concentration M and M. It is seen that the increase in electrolyte concentration leads to the increase in both coexisting isotropic and nematic densities.

Figure 11

The isotropic–nematic phase diagram in terms of coexisting densities of actin for different average lengths of filaments and two values of electrolyte concentrations M (dashed lines) and M (solid lines). The persistence length is μm, and the normalized standard deviation . The experimental data (triangles) were retrieved from [26]. The coexisting isotropic and nematic densities are plotted in black and red colors, respectively.

Clearly, the results for M give the best fit against the experimental data. This behavior is due to a substantial decrease of the effective diameter with increasing the electrolyte concentration (see Figure 10), as well as because the coexisting density is inversely proportional to the effective diameter (see Equation (73)).

In Figure 12, we investigate the I–N phase diagram for a fixed average length in terms of the dependence of the coexisting densities of actin on the concentration of monovalent ions .

Figure 12

The isotropic–nematic phase diagram in terms of coexisting densities of actin for different values of the monovalent ions’ concentrations . The average lengths of filaments are = 800 (solid lines) and 3000 (dashed lines); the persistence length μm, and the normalized standard deviation . The coexisting isotropic and nematic densities are plotted in black and red colors, respectively.

It is seen that the increase of the concentration leads to the increase of both isotropic and nematic coexisting densities . On the other hand, for a fixed value of , the increase of the average size decreases both coexisting densities . Certainly, it is easier for long rods or filaments to form a nematic order rather than for shorter ones. Similarly, for longer filaments (), the I–N phase transition occurs at lower densities compared to short ones () with the same concentration .

Another quantity playing a key role in the I–N phase diagram is the twisting parameter. In some studies [12], the terms with exponential integral in the formulas analogous to Equations (A14) and (64) were dropped for simplicity, such as the electrostatic effects on the phase equilibrium can be described exclusively by the effective diameter and the twisting parameter: which is a measure of a twist of equally charged cylinders, which hinders the formation of the nematic order. On another hand, the increase of the effective diameter promotes the nematic order. The competition between these two effects affects the I–N phase equilibrium in the electrolyte solution [12]. We plot the dependence of the twisting parameter t on the concentration of monovalent ions in Figure 13.

Figure 13

The twisting parameter t as a function of the monovalent ions’ concentration .

For actin filaments, the parameter t decreases with the increase of the concentration . We also rescale the I–N phase diagram in terms of coexisting densities of actin for different values of the twisting parameter t in Figure 14.

Figure 14

The isotropic–nematic phase diagram in terms of coexisting densities of actin for different values of the twisting parameter t. The average lengths of filaments are = 800 (solid lines) and 3000 (dashed lines); the persistence length μm, and the normalized standard deviation . The coexisting isotropic and nematic densities are plotted in black and red colors, respectively.

It is seen that both isotropic and nematic coexisting densities decrease with increasing twisting parameter t, whereas, for a fixed value t, the increase of the average size decreases both coexisting densities . Indeed, according to Figure 10 and Figure 13, the increase of the twisting parameter t decreases the electrolyte concentration and, correspondingly, increases the effective diameter . Furthermore, our results plotted in Figure 14 show that stabilizing effects on the nematic state of larger effective diameters dominate destabilizing twisting effects with larger twisting parameters t [11].

4. Discussion

In the present article, we developed a classical density functional theory to calculate the isotropic–nematic phase diagram for actin filaments in physiological electrolyte solutions. The approach is based on an extension of Onsager’s second-order theory for monodisperse, charged rigid rods [10]. The key ingredients in this work are a unique definition of the orientational free energy and the size and angular distribution functions, which properly account for the polydispersity and semiflexibility of the actin filaments. Unlike many studies on polydisperse rods where the size distribution function is an input parameter [17,19], actin filaments self-aggregate to produce a size distribution function that satisfies the equilibrium condition given by Equation (27). Our results reveal that has the form of Schulz distribution function, which is in accordance with experimental work performed on actin filaments [25,26]. To model an angular distribution function in the nematic phase , we generalized the trial function introduced in [17] for polydisperse filaments, which depends on two parameters and . Additionally, we accounted for the filament semiflexibility by extending the formula for orientational free energy introduced in [13] to the case of a polydisperse system. Finally, we calculated the isotropic–nematic phase diagrams for several values of the normalized standard deviation of Schulz distribution , the persistence lengths of actin filaments P, and the concentrations of monovalent ions . We compared the obtained results with the corresponding experimental data presented in [26]. We found that the set of parameters μm and M gives the best match against the experiment.

In principle, our method can be applied with suitable modifications to other associating particles or exchange colloids that form charged rod-shaped aggregates such as DNA and ionic micelles. In the case of microtubules, the theory should be modified to account for the impact of the lumen on the their electrostatic interactions. It can be also extended to study phase diagrams in confined spaces, as well as under other electrolyte and filament conditions.

For instance, the orientational trial function used in this article is valid for weak polydispersity, i.e., when the normalized standard deviation is small. However, the solution can also be obtained for highly polydisperse systems, while the computational burden dramatically increases.

Additionally, in our study, we followed the conditions of the experiment, where the average length of actin filaments was fixed both in the isotropic and nematic phases [26]. As a result, the last two terms in the expression for free energy (43) did not contribute to the phase equilibrium. In principle, we can also apply our theory for polydisperse systems having different average lengths in the coexisting isotropic and nematic phases. In this case, the last two terms in the equation for free energy (43) cannot be disregarded. Additionally, the approximation introduced in this work for the excess standard chemical potential in Equation (39) should be properly changed/modified if the length distribution of the polydisperse system does not follow Schulz distribution function .

On the other hand, we used Brenner–Parsegian’s formula [27] to approximate the electrostatic part of inter-filament potential in Equation (15). This approximation accounts for water solvent as a dielectric medium characterized by a constant dielectric permittivity. In principle, a more accurate, distance-dependent dielectric medium can be used to capture the high water polarization near the filament surface [32].

Furthermore, actin filaments may form bundles and networks due to the attraction forces between filaments produced by linker proteins or divalent ions in eukaryotic cells. In our method, we accounted only for repulsive interaction between filaments, i.e., the hard core plus electrostatic repulsion in Equation (14). The attraction between filaments can be included, for example, using square-well-type potentials [2,33].

Our approach is also able to describe charged filaments in pathological conditions. For instance, pH changes and G-actin mutations were accounted for from molecular structure models for actin filaments [34].

Overall, more realistic models can be developed for filament–filament and filament–binding protein interactions in constrained spaces (encapsidated polyelectrolytes) [35,36]. For instance, a variety of geometries such as spherical and cylindrical capsids can be used to model different cellular compartments such as dendrites (spines and filopodia), soma, axons, and terminals.