Phase-field collagen fibrils: Coupling chirality and density modulations

To describe the interaction between longitudinal density modulations along collagen fibrils (the D band) with a radial twist field of molecular orientation (double twist), we couple a phase-field crystal (PFC) with liquid-crystalline free energies to obtain a hybrid model of equilibrium collagen fibril structure. We numerically compute the resulting axial and radial structure. We find two distinct fibrillar phases, L and C, with a coexistence line that ends in an Ising-like critical point. We propose that the coexistence between these phases can explain the bimodal distribution of fibril radii that has been widely reported within tendon tissues. Tensile strain applied to our model fibrils straightens the average fibrillar twist and flattens the D-band modulation. Our PFC approach should apply directly to other longitudinally modulated chiral filaments, such as fibrin and intermediate filaments.

implies a D-band period d ∝ cos ψ (r) within an individual fibril [2], whereas only a single D-band period is observed in experiment. One previous model of both D band and double twist therefore assumes a radially constant twist [19], though this is energetically unfavorable for the double-twist field [13]. A second model has an approximately constant twist gradient [9], which is energetically preferable for the double twist but implies local stretching or compression of the D band [2]. These two models of the radial twist ψ (r) represent the opposite limits of a stiff or soft D band, respectively. How can we explore the coupling of the D band and double twist more generally?
While atomistic molecular dynamics (MD) simulations of the D band, e.g., Ref. [20], can explore axial strain [21], they currently ignore radial twist. This is because a continually varying ψ (r) precludes unit cells with small numbers of molecules. Fortunately, coarse-grained phase-field crystal (PFC) approaches for addressing periodic modulations of crystalline materials [22] are compatible with coarse-grained models of the radial twist. PFC models impose density modulations with coarse-grained fields, and so do not require the inefficiently short atomic timescales of MD. PFC models allow us to quantitatively explore both the mechanical and structural properties of the fibril for arbitrary values of the D-band stiffness.
Phase-field crystal fibril. PFC theory adds terms to the coarse-grained free energy to generate a periodic structure. PFC is particularly simple in one dimension, such as along the collagen fibril. From the Hodge-Petruska model, molecules pack along their long axis with a periodd = 67 nm in the absence of molecular twist. We can write the PFC contribution to the free-energy per unit volume as where we integrate over a cylindrical fibril of radiusR and lengthL, andφ is the amplitude of modulations due to the D band.∇ is the gradient operator in the direction parallel to the local molecular orientation. The first integral of Eq. (1) is minimized when the modulations have the same local periodicity asd , and˜ characterizes the D-band stiffness. For χ 2 > 0, the second integral is minimized whenφ 2 is nonzero, which determines a preferred nonzero D-band amplitude.ω characterizes the energetics of D-band formation. We can further simplify these PFC contributions. Under a single-mode approximation for the D band, we takẽ φ(z) =δ cos(ηz), whereη is the observed D-band wave number. Furthermore, we work within the ansatz that the molecular orientation is determined by the twist field, n = sin ψ (r)φ + cos ψ (r)ẑ [11,13]. Because the local orientation is not along the fibril axis, this couples the D-band periodicity to ψ (r) through the gradient in the frame parallel to n,∇ = cos ψ (r)∂/∂z. We then obtain a simpler expression, To this we add the volume-averaged Frank free-energy density for the double-twist director field (see Ref. [13]), whereK 22 ,K 33 , andk 24 are the usual Frank elastic constants (twist, bend, and saddle splay, respectively). We also add the average surface energy per unit volume due to free interfaces, given byẼ surf = 2γ /R, whereγ is the surface tension. For the remainder of this Rapid Communication, variables without a tilde will be dimensionless. We do this by measuring energies in units ofK 22q 2 , measuring density in units of χ , measuring radius in units of 1/q, and the D-band wave number in units of 1/d . Combining E pfc , E Frank , and E surf , we obtain the total average free-energy density of the fibril as a function of radius R, D-band modulation amplitude δ, D-band modulation period 2π/η, and twist angle field ψ (r), There are five dimensionless parameters that control the behavior of our system. K 33 ≡K 33 /K 22 and k 24 ≡k 24 /K 22 characterize the bend and saddle-splay elastic constants. Consistent with our previous work, we will fix K 33 = 30 [13,23]. γ ≡γ /(K 22q ) controls the surface tension. ≡ 2˜ χ 2 /(3K 22q 2d 4 ) controls the coupling strength between the D band and the molecular twist and can be related to the Young's modulus at zero twist. ω ≡ 2ωχ 4 /(3K 22q 2 ) controls the strength of the D-band double-well potential, which cannot be any larger than the polymerization energy of collagen fibrils [24].
We minimize E tot with respect to ψ (r) for chosen initial values of R, η, and δ using a numerical implementation of the corresponding Euler-Lagrange equations [11], and obtain ψ 0 (r). We then minimize E 0 ≡ E [R, η, δ; ψ 0 (r)] with respect to R to obtain a global (thermodynamic) minimization of E tot . Our numerical minimization routines are available online via GitHub [25]. The default parameter values, which apply unless otherwise stated, are γ = 0.04, k 24 = 0.5,q = 4 μm −1 , = 600, and ω = 20, and are discussed below.
Coexistence. In Fig. 1(a) we show the phase diagram of the PFC fibril model in the γ , k 24 plane. The D-band modulation leads to a coexistence line (indicated by the thick black line) between qualitatively distinct fibril phases, which ends at a critical point at γ c 0.035, k c 24 0.26. This coexistence line, along which the fibril radius and twist field changes discontinuously, is not observed without the PFC terms [13].
For the default parametrization, indicated by a white star in Fig. 1(a), we show the two coexisting twist-field solutions which minimize E tot in Fig. 1(b). For smaller R (blue line), the equilibrium twist field ψ (r) has an approximately constant gradient and we call it the "linear twist" phase (L). For larger R (orange line), ψ (r) has a large region of approximately constant twist and we call it the "constant twist" phase (C). From a molecular twist perspective, the C phase can be viewed as a "core-shell" structure with the core showing a strong linear twist gradient while the shell shows a constant twist, albeit with a narrow region of additional twist gradient at the surface.
Qualitatively, fibrils are energetically stabilized by surface twist through k 24 as well as by the D-band amplitude δ. Larger values of surface tension γ drive larger R to reduce the surface area per unit volume. These larger radii also lead to increasingly large D-band elastic energy through , which can then be reduced by a region of constant twist.
The L and C twist fields are visualized schematically in Figs. 1(c) and 1(d), respectively. The inset circles within Fig. 1(a) indicate the axial D-band strain in a cross section of unstretched fibrils corresponding to the blue (L) and orange (C) curves of Fig. 1 For the L phase, most of the fibril is strained ranging from the center under axial compression of 0.3% to the surface under tension of 0.1%. For the C phase, only the center and surface are significantly strained.
Our coarse-grained free-energy approach does not include any fluctuations, so we expect mean-field critical exponents [26] to describe the discontinuities across the coexistence line sufficiently close to the critical point. In Fig. 2, we show the difference in radii of the coexisting linear and constant twist phases, R * L and R * C , respectively, versus the distance from the critical point t ≡ k 24 /k c 24 − 1. We find an Ising-like mean-field critical exponent β = 1/2, as indicated by the dashed black line. In the inset, we show the splitting of R into R * L and R * C near the critical point. Near t 3 we see that the ratio R * C /R * L can be as large as 100.
Remarkably, the coexistence of widely different fibril radii has been reported [27] (see also Refs. [28,29]) within tendon samples. The ratio of radii increases with age, and approaches 5 for older tendons [27]. We propose that tendon fibrils are in coexistence within our equilibrium model. Indeed, our default parameters are chosen to approximately recapitulate the tendon fibril properties. γ = 0.04 and k 24 = 0.5 on the coexistence line are chosen to recover R * C /R * L 4. Taking q = 4 μm −1 is consistent with in vitro studies of chiral nematic collagen solutions [30], and leads to radiiR C ≈ 106 nm andR L ≈ 27 nm.
Elastic properties. We can probe the tensile response of our model fibril to an applied strain . To do this, we axially strain the fibril and D band while conserving volume by imposing η = η eq /(1 + ) and R = R eq / √ 1 + , where η eq and R eq are the unstrained values. We then minimize E tot with respect to δ and ψ (r) at these strained values of η and R. Figure 3 shows the resulting stressσ = dE/d [ Fig. 3(a)], D-band amplitude δ [ Fig. 3(b)], and volume average twist ψ (r) [Fig. 3(c)] versus strain . As increases, the fibrils stiffen until they reach a maximum stress, after which the fibrils are unstable. This stiffening is delayed in strain and Untwisting [decreasing ψ , as seen in Fig. 3(b)] and flattening of the D-band modulation (as seen in Fig. 3) occur as strain increases. The untwisting and flattening with strain result from the coupling between the D band and the radial twist in our model.
By choosing = 600 we effectively determineỸ high ≈ 100 MPa to be comparable with the maximal slope ofỸ 30 MPa observed in noncrosslinked fibrils [31], usingK 22 = 6 pN [13]. However, the observation of a significant low-slope regime is also observed experimentally [31] with a compara-bleỸ low 1-5 MPa. Only ω is relatively unconstrained, but that is because (data not shown) our results are qualitatively independent of it. We take ω = 20.
Discussion. The phase coexistence (Fig. 1) that emerges from the interaction between the axial D-band modulation and the radial double twist results in strikingly different fibril radii coexisting in thermodynamic equilibrium (Fig. 2) and provides a natural explanation for experimentally observed bimodal radius distributions of tendon fibrils [27][28][29].
Coupling the D band and twist fields also leads us to predict that elastic strain straightens the twist field [ Fig. 3(b)] and flattens the D-band amplitude [ Fig. 3(c)]. Qualitatively similar strain straightening has been observed in recent synchrotron x-ray scattering studies of corneal fibrils [32]. Similarly, recent creep studies of modestly strained tendons completely degrade the D band after several hours [33], which is also consistent with our results.
We have found that coexisting phases have distinct twist fields ψ (r), with the smaller radii fibrils qualitatively similar to constant gradient models discussed by Raspanti [2], while the larger fibrils are qualitatively similar to constant-tilt models proposed for corneal fibrils [10,19,34]. We find that the equilibrium axial strain [ Fig. 1(a) insets] within these fibrils is quite distinct in local variation, though not in magnitude. Both L and C fibrils are under considerable compression at the center and tension at the surface.
We also find that the twist field has a dramatic impact on the elastic response of fibrils, as seen when comparing the L-and C-phase fibrils in Fig. 3(a). This contrast may not be observed for in vivo fibrils where covalent crosslinking greatly stiffens fibrils [35]. Nevertheless, carefully controlled assembly conditions in vitro can achieve a wide range of fibril radii and twist fields [13]. While coexistence has not been directly reported for in vitro assembly, there is a striking tenfold increase in fibril radius over a narrow range of precursor concentration from 75 to 100 mg/ml [3]. This is consistent with the discontinuity observed along coexistence in Fig. 2 and we would therefore expect significant differences in the stressstrain curves as well. Further exploration of in vitro assembled fibrils may also offer a way of exploring the critical point of Fig. 1(a), and of characterizing fluctuation effects there.
We can obtain a general relationship between the observed twist-field and D-band periodicity by minimizing the first term in Eq. (2) with respect to d/d || = 2π/η, where this equation holds for all parameter values. Approximating the twist field as constant, we obtain d C /d || = cos ψ (R). Approximating the twist field as having a constant gradient, i.e., ψ (r) = ar, we obtain d L /d || = 1 − ψ (R) 2 /4 + 19ψ (R) 4 /288 + O(ψ (R) 6 ). While we have not discussed corneal fibrils so far, ψ (R) 0.3, while d cornea /d || = 64 nm/67 nm 0.955 [2]. This indicates that corneal fibrils are in the C phase, confirming results from a previous electron tomography study [36]. While we have described the equilibrium fibril structure that minimizes E tot , we can expand around that minimum and consider equilibrium fluctuations. In particular, consider D-band periods d = 2π (1 + u)/η eq , where 2π/η eq is the equilibrium period. We have that E (u) − E (0) ≈ 1/2Ỹ u 2 . Multiplying by the volume in a single D-band period, and using the Boltzmann distribution we obtain P(u) ∝ With k B T ≈ 4.114 × 10 −3 pN μm,R, andỸ , we obtain σ C σ L 0.001. Atomic force microscopy (AFM) studies of uncrosslinked fibrils have reported a narrow, approximately Gaussian distribution of D-band spacings with a fractional width σ d/d ≈ 1% [18], approximately tenfold larger than our prediction. While this indicates that the coupling we have used between the D band and twist field is plausible, it suggests that additional physics is needed to describe the observed D-band variability. One possibility is that the reported longitudinal variability of mechanical properties along collagen fibrils [37] may lead to an increased D-band variability as well. We have seen that using a one-dimensional phase-field crystal (PFC) approach to couple longitudinal D-band modulations with the radial twist field leads us to predict two distinct structural phases for collagen fibrils: linear twist (L) and constant twist (C). We find a phase coexistence between L and C that could describe the bimodal radius distribution of tendon collagen fibrils observed in vivo. It should be possible to explore the critical point with in vitro fibril assembly systems.
This PFC approach to the coupling of density modulations with orientation fields could also be applied more generally in chiral materials [38,39]. It should also directly apply to a number of chiral self-assembling fibrillar systems that exhibit longitudinal modulations, with only the parametrization to be determined. These include fibrin [40,41], keratin filaments [42], and nuclear lamin paracrystals [43].