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.


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 the isotropic-nematic phase diagram behavior. This approach is essential to study the selforganization 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 Appendices A-E.

Materials and Methods
We considered a solution of polydisperse actin filaments, each of them having a length L = l m ν, where l m 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 dω = sin θdθdϕ with the z direction taken along the nematic axis.
The size-angular density distribution function of actin filaments ρ ν (ω) is given by the following expression: ρ ν (ω) = ρn ν η ν (ω), (1) where ρ = N/V 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, n ν 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 n ν and η ν (ω) are 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 n ν is characterized by the filament average size < ν >, namely the average degree of polymerization, as well as the normalized standard deviation σ, such as and where

The Free Energy Density Functional
We define the dimensionless Helmholtz free energy of cytoskeleton filaments f as follows: where F is the free energy, β = 1/k B T the inverse thermal energy, k B the Boltzmann constant, T the temperature, f ig the dimensionless ideal gas free energy, f int the dimensionless energy due to the inter-filament interactions, and µ (0) ν the standard chemical potential of ν-sized filaments.
The expressions for f ig and f int as functionals of the density distributions n ν and η ν (ω) are provided below.

Ideal Gas Free Energy f ig
The ideal gas free energy per volume of the system V can be written in the following form: where Λ ν = h √ 2πm ν k B T is the thermal de Broglie wavelength of ν-sized filaments, h is the Planck constant, and m ν = νm 1 stands for the ν-sized filament mass. Substitution of Equation (1) into Equation (8) yields where f ig = βF ig /N is the ideal gas free energy per filament (N = ρV), Λ 1 = h √ 2πm 1 k B T the monomer thermal de Broglie wavelength, and s (or) ν the filament orientational entropy:

Interaction Free Energy f int
We write the interaction free energy per volume V in a mean-field fashion, namely where B ν 1 ν 2 (ω 1 , ω 2 ) represents the cluster integral In Equation (12), w 12 is the interaction potential between two filaments and r 12 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 L/D → ∞ [10].
Substitution of Equation (1) into Equation (11) yields the following expression for F int as a functional of the densities n ν and η ν (ω): To calculate the cluster integral B ν 1 ν 2 (ω 1 , ω 2 ) in obvious form, we model the interaction potential between two charged rod filaments w 12 (r 12 , ω 1 , ω 2 ) as follows: where D is the bare rod diameter and γ 12 the angle between them. The function Γ(r 12 , γ 12 ) represents the following repulsive electrostatic inter-filament interaction potential: [27] Γ(r 12 , where In Equation (16), 2 stands for the inverse Debye length, K 1 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 ξ i and ρ i 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 |ξ i | = 1. Substitution of Equations (14) and (15) into Equation (12) results in the following expression for the cluster integral B ν 1 ,ν 2 (ω 1 , ω 2 ) where D e f f and B * (sin γ 12 ) 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 h ν 1 ν 2 is given by We note that the parameter h ν 1 ν 2 accounts for the orientational averaged contributions coming from the excluded volume and charge density interactions between filaments of size ν 1 and ν 2 , 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 450 ÷ 5000, which correspond to filament lengths in the range of 1.2 ÷ 13.5 µm only. Thus, the linearized Poisson-Boltzmann equation solution is justified in our study.

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 B ν 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 n ν 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 {n ν } 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 "pseudophases". 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 f 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 n ν : where ∆µ Finally, we substitute the parameter y ≡ − ln n 1 C into Equation (30) to obtain the following master equation for n ν : Similarly, for a fixed distribution function n ν , 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 E ν is a constant that satisfies the following relationship: To obtain the expression for E ν , we integrate both sides of Equation (34) with respect to ω 1 and use the normalization conditions (2) and (3). Finally, we replace the expression obtained for E ν 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 n ν 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., n ν = δ(ν− < ν >), where δ is a Dirac delta function. Substitution of this expression for n ν 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 η(θ) = α cosh(α cos θ)/(4π sinh α). 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 n ν . Since η(θ) does not depend on ν, it follows from Equation (19) that h νν − h 1ν = 0. Additionally, we use in Equation (32) the ansatz: to obtain a size distribution function n ν in the asymmetric Schulz-Zimm form [29]: where z represents the conventional polydispersity parameter, whereas C = y z+1 Γ(z+1) 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 σ = 1 √ z+1 , whereas the parameter y depends on the values < ν > and z only, i.e., y = (z + 1)/ < ν >.
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 ∆ν → 0 and α ν → α.

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 µ ν + ν 1/3 and the terms f or , f int , f mix represent the orientational, interaction, and mixing free energies, respectively. The expressions for these energies are provided below.

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 < L > P, where P is the filament persistence length. Indeed, in the monodisperse limit σ → 0, the correct limit f or → ln α − 1 is obtained for any value of the parameter γ [13]. For flexible filaments, where < L > P, 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 < L >= l m < ν >. Certainly, in the monodisperse limit, σ → 0 the correct limit f or → L P α 4 is obtained [11]. To calculate the orientational free energy for any persistence length, we combine these two asymptotic cases for f or using the following interpolating formula: The details on these calculations are provided in Appendix C.

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 z → ∞, the expression for g z (θ) 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 f int in the monodisperse case: with

Mixing Free Energy
The mixing free energy in Equation (43) can be written as: This expression is not well defined in the monodisperse limit (z → ∞) since the density distribution n ν → δ(ν− < ν >) and the mixing free energy diverges, rather than going to zero. To overcome this shortcoming, we extract the density ρ-independent term f mix (z → ∞) 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., f mix → 0 for σ → 0 (z → ∞). More details on these calculations are provided in Appendix E.

Free Energy in the Isotropic Phase
In the isotropic phase, there is no preferential direction for the filament orientations; thus, the orientational distribution function η ν (θ) = 1/(4π) 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 B I 2 is the second virial coefficient for polydisperse filaments in the isotropic phase Finally, substitution of Equation (A14) into Equation (63) generates the following expression for h in the isotropic phase where E 1 (x) is the elliptic integral given by Equation (A9).

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 µ ν = ∂F ∂N ν is the chemical potential of ν-sized filaments and µ 1 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 µ 1 , 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 G N = F N + Π ρ in Equation (69) yield In Equations (66) and (70), we use the formula Π = ρ 2 ∂ ∂ρ ( F N ) for practical calculations of the osmotic pressure.
Using Equation (4), we can also write the average filaments' length in Equation (70) as follows: where ρ A = ∑ ν≥1 νρn ν 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., ρ ≈ ρ G . 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 ρ A was accompanied by a change in the concentration of gelsolin ρ G 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 f mix and the term with µ (0) ν 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 f or 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, f ≡ f (c; α, γ). By minimizing f , we obtain the equilibrium parameter functions α eq (c) and γ eq (c) for each value of c. Thus, the expression for the nematic free energy at equilibrium becomes f eq (c) = f (c, α eq (c), γ eq (c)). The minimization is carried out using the downhill simplex method in multiple dimensions [31]. Finally, the coexisting dimensionless densities c I and c N are obtained from the coexistence conditions (66) and (68), whereas the corresponding actin densities ρ I A and ρ N A are obtained from Equations (71) and (73). In principle, the chemical potential of the non-depleted unimers (G-actin) µ 1 and the depleted actin concentration ρ A are related by Equation (70). We calculate the Helmholtz free energy in Equation (72) as a function of ρ A in the isotropic or nematic phase. Subsequently, we substitute the result into the r.h.s. of Equation (70) to obtain the relationship µ 1 = µ 1 (ρ A ) or, vice versa, ρ A = ρ A (µ 1 ). In the last case, the chemical potential of unimers µ 1 can be regarded as an input parameter, which is set by the reservoir condition.

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 D = 80Å, the monomer units' length l m = 27Å, the actin molar weight m A = 42 kDa, and the linear filament charge density λ = −4 e nm . We calculated the isotropic and nematic coexisting dimensionless densities for rigid monodisperse rods with length L = 1 µm for the electrolyte's concentration c e = 0.1 M to compare our results with those obtained by Borukhov et al. [2]. Our result was c I = 3.82 and c N = 5.10, whereas the corresponding one in [2] was c I = 4.25 and c N = 9.98 (see Figure 3 in [2]). The source of such discrepancy is associated with the cone approximation used in [2] for the angular distribution function.

Distribution Functions
In this section, we present the analysis on the distribution functions n ν and η ν (θ) in the I-N phase coexistence for polydisperse, semiflexible filaments. We used a persistence length P = 18 µm and an electrolyte concentration c e = 0.1 M. In Figure 1, we plot the Schulz length distribution function n ν for different average sizes < ν > in a weakly polydisperse system with normalized standard deviation σ = 0.5 (z = 3). 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´∞ 1 n ν dν = 1 can be split into the sum of two integrals´< ν> 1 n ν dν = 0.57 and´∞ <ν> n ν dν = 0.43 (z = 3), 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 n ν → δ(ν− < ν >).
In Figures 2 and 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 < ν >.  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. 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 ρ ν (θ)/ρ = n ν η ν (θ) on the filament size ν for several angles θ is displayed in Figure 6. 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 θ = 0. By increasing the angle θ, the amount of filaments decreases and is directed along the direction with angle θ.

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 ρ A for different average lengths < ν > and two values for the normalized standard deviation of Schulz distribution function: σ = 0 (monodisperse) and σ = 0.5 (weakly polydisperse). In these calculations, we used the values P = 18 µm and c e = 0.1 M. In the same figure, we plot the experimental data extracted from Figures 4 and 5 in [26].
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 P 2 = (3 cos 2 θ − 1)/2 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 S = 1. 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.
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 P = 7 µm and P = 18 µm, in physiological conditions (electrolyte concentration c e = 0.1 M).  . The isotropic-nematic phase diagram in terms of coexisting densities of actin ρ A for different average lengths of filaments < ν > and two values of persistence lengths P = 7 µm (dashed lines) and P = 18 µm (solid lines). The normalized standard deviation is σ = 0.5, and the electrolyte concentration c e = 0.1 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 P = 18 µm give the best fit against the experimental data. This result stems from the impact of the semiflexibility on the orientational free energy f or . It can be seen from Equation (49) that the increase of the persistence length P leads to the decrease of the orientational term f or , 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 D e f f on the concentration of monovalent ions c e .
Our results show an exponential decay of D e f f with increasing c e . 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 0.4 M, the effective diameter of filaments D e f f is smaller than the diameter of filaments D = 80Å, and the calculations may become inaccurate.
Furthermore, we show in Figure 11 the isotropic-nematic phase diagram for two values of electrolyte concentration c e = 0.01 M and c e = 0.1 M. It is seen that the increase in electrolyte concentration leads to the increase in both coexisting isotropic and nematic densities.
Clearly, the results for c e = 0.1 M give the best fit against the experimental data. This behavior is due to a substantial decrease of the effective diameter D e f f with increasing the electrolyte concentration c e (see Figure 10), as well as because the coexisting density ρ A 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 ρ A on the concentration of monovalent ions c e . Figure 11. The isotropic-nematic phase diagram in terms of coexisting densities of actin ρ A for different average lengths of filaments < ν > and two values of electrolyte concentrations c e = 0.01 M (dashed lines) and c e = 0.1 M (solid lines). The persistence length is P = 18 µm, and the normalized standard deviation σ = 0.5. The experimental data (triangles) were retrieved from [26]. The coexisting isotropic and nematic densities are plotted in black and red colors, respectively. It is seen that the increase of the concentration c e leads to the increase of both isotropic and nematic coexisting densities ρ A . On the other hand, for a fixed value of c e , the increase of the average size < ν > decreases both coexisting densities ρ A . Certainly, it is easier for long rods or filaments to form a nematic order rather than for shorter ones. Similarly, for longer filaments (< ν > = 3000) , the I-N phase transition occurs at lower densities ρ A compared to short ones (< ν > = 800) with the same concentration c e .
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 E 1 (t) 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 D e f f 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 D e f f 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 c e in Figure 13. For actin filaments, the parameter t decreases with the increase of the concentration c e . We also rescale the I-N phase diagram in terms of coexisting densities of actin ρ A for different values of the twisting parameter t in Figure 14. It is seen that both isotropic and nematic coexisting densities ρ A decrease with increasing twisting parameter t, whereas, for a fixed value t, the increase of the average size < ν > decreases both coexisting densities ρ A . Indeed, according to Figures 10 and 13, the increase of the twisting parameter t decreases the electrolyte concentration c e and, correspondingly, increases the effective diameter D e f f . Furthermore, our results plotted in Figure 14 show that stabilizing effects on the nematic state of larger effective diameters D e f f dominate destabilizing twisting effects with larger twisting parameters t [11].

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 f or and the n ν 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 n ν is an input parameter [17,19], actin filaments self-aggregate to produce a size distribution function n ν that satisfies the equilibrium condition given by Equation (27). Our results reveal that n ν 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 f or introduced in [13] to the case of a polydisperse system. Finally, we calculated the isotropicnematic 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 c e . We compared the obtained results with the corresponding experimental data presented in [26]. We found that the set of parameters P = 18 µm and c e = 0.1 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 ∆µ 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 filamentbinding 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.  Data Availability Statement: Some or all data, models, or code that support the findings of this study are available from the corresponding author upon reasonable request.

Conflicts of Interest:
The authors declare no conflict of interest.
Appendix A. Cluster Integral B ν 1 ν 2 (ω 1 , ω 2 ) in Equation (12) The volume integration in Equation (12) can be written as follows: where the x axis is chosen along the direction of the closest distance between rods, the z axis along the direction of first rod, and the y axis in the perpendicular direction to the z − x plane. L 1 , L 2 are the lengths of two rods and γ 12 the angle between the axis of two rods. Substitution of Equation (A1) into Equation (12) yields Integration in Equation (A2) on the y − z plane leads to Substitution of Equation (14) into Equation (A3) yields Substitution of Equation (15) into the integral in the r.h.s. of Equation (A4) yields the following expression: where The analytic solution of Equation (A5) is given by where γ E is Euler's constant: and E 1 (x) is the exponential integral: Next, we substitute Equations (A5)-(A7) into Equation (A4) to obtain We introduce the effective diameter D e f f to rewrite Equation (A10) as follows: where Finally, Equation (A11) can be written as where

Appendix B. Schulz Distribution Function in Equation (40)
Here, we provide some properties of Schulz distribution function n ν . We define the average size in the k-th power < ν k > as follows: Substitution of Equation (40) into Equation (A15) gives where Γ is the Gamma function: Equation (A16) yields the following expressions for k = 0, 1, 2: Substitution of Equations (A19) and (A20) in Equation (5) generates the following relation between the normalized standard deviation σ and parameter z:

Appendix C. Orientational Energy f or in Equation (49)
The orientational free energy f or in the nematic phase for monodisperse systems can be expressed as [13]: where L is the filament length and P the persistence length. Khokhlov  In the monodisperse case, i.e., for σ = 0, Equations (48), (A25), and (A26) reduce to those provided in [13]. It was also suggested in [13] to substitute (α − 1)N p instead of αN p in the monodisperse version of Equation (A26) to have reliable results for small α, even when N p ∼ 1.
where γ E = 0.5772 is Euler's constant, Γ(z + 1) = z! the gamma function, and H z = 1 + 1 2 + 1 3 + . . . + 1 z the z-th harmonic number. Substitution of Equation (A42) into the second term in the r.h.s. of Equation (A40) gives We use the following asymptotic formula for the harmonic number: to calculate the free energy f mix in the monodisperse limit. We substitute Equation (A44) into Equation (A43), and we use Equations (A19) and (A21), and Stirling's formula Γ(z + 1) → √ 2πz(z/e) z to finally obtain