Silk is a semidilute solution of randomly coiled associating polypeptide chains that crystallize following the stretch-induced disruption, in the strong extensional flow of extrusion, of the solvation shell around their amino acids. We propose that natural silk spinning exploits both the exponentially broad stretch distribution generated by associating polymers in extensional flow and the criterion of a critical concentration of sufficiently stretched chains to nucleate flow-induced crystallization. To investigate the specific-energy input needed to reach this criterion in start-up flow, we have coupled a model for the Brownian dynamics of a bead-spring-type chain, whose beads represent coarse-grained Gaussian chain segments, to the stochastic, strain-dependent binding and unbinding of their associations. We have interpreted the simulations with the aid of analytic calculations on simpler, tractable models with the same essential physical features. Our simulations indicate that the associations hamper chain alignment in the initial slow flow, but, on the other hand, facilitate chain stretching at low specific work at later, high rates. We identify a minimum in the critical specific work at a strain rate just above the stretch transition (i.e., where the mean stretch diverges), which we explain in terms of analytical solutions of a two-state master equation. We further discuss how the silkworm appears to exploit the chemical tunability of the associations to optimize chain alignment and stretching in different locations along the spinning duct: this delicate mechanism also highlights the potential biomimetic industrial benefits of chemically tunable processing of synthetic association polymers.
I. INTRODUCTION
The manufacturing of both natural and artificial polymer-based fibers relies on flow-induced crystallization in nonlinear rheological conditions [1–6]. The energy input required by this process may be significantly reduced in natural silk-spinning, though the mechanism by which this efficiency is achieved has been far from clear [7]. There is evidence, however, that locally tailored macromolecular interactions are involved [8–11]: The silk protein, of which the conformation in solution closely resembles a random coil [12], self-assembles in flow in aqueous conditions under energy requirements orders of magnitude lower than its synthetic counterparts [7]. It has been hypothesized that flow-induced stretching of the chain disrupts a solvation layer and, in turn, enables crystallization to commence [7,13,14]. This mechanism was supported by molecular dynamics (MD) simulations [15–18] and was employed to induce crystallization of synthetic poly-ethylene oxide by flow at similarly low energetic requirements as silk, however, at much higher molecular weight and/or strain rates [13]. The low-energy mechanism for natural silk-spinning, therefore, remains to be identified. Clues may be present in the subtle electrostatically modified rheo-physics of associating polymers [19–28].
We previously found, in collaboration with Laity and Holland, that the silk protein exhibits calcium bridges that act as intermolecular reversible cross-links [8,9]. Such associations, sometimes referred to as “stickers” that can be in a bound/closed or unbound/open state [19], shift the alignment-to-stretch transition to smaller strain rates by replacing the usual Rouse relaxation dynamics for “sticky-Rouse” relaxation [19–28]. Inspired by these observations, we envision a mechanism of flow-induced crystallization (FIC) where the reversible network is initially equilibrated (in Stark contrast to the typical mechanism for the solgel transition of associating polymers, where shear flow breaks metastable intramolecular associations and facilitates the formation of an intermolecular network [29–31]). In our case, strong flow stretches the “bridging” strands between the stickers [32,33]. This stretch, in turn, aligns the strands at the scale of the Kuhn segments (in which water-soluble systems may disrupt the solvation layer [7,13]), so nucleating crystals as structural elements within (silk) fibers. It will turn out that such a picture contains within it a mechanism for the super-efficiency of natural silk-spinning through a surprisingly strong heterogeneity in the chain stretch distribution.
While this mechanism seems plausible, it is not evident how this process may be controlled and/or optimized by the number of stickers per chain and by their lifetime. Intriguingly, however, it has been observed that the Bombyx mori silkworm tunes the sticker lifetime and, hence, the (non)linear rheology, before and during spinning through local chemical control variables. Prior to pupation, i.e., when the silkworm is not required to spin a cocoon, the silk is stored in the gland at high viscosity using long sticker lifetimes [8,9]. When pupation commences, potassium cations are added to decrease the sticker lifetime and reduce the viscosity [8,9].
We first hypothesize, as schematically indicated in Fig. 1, that the decrease in the sticker lifetime decreases the specific work needed to align the chains in the direction of the flow field well upstream from the spinnerette. The group of Holland also discovered that the structural features of the silk fiber are significantly enhanced through a gradient in the H along the spinning duct, suggesting an exquisitely controlled local rheology [34]. While lower H may induce partial folding of the protein [12], it is also expected to enhance the lifetime of the stickers. Crucially, inspired by our previous finding that broad conformational distributions emerge due to the stochastic nature of binding and unbinding stickers [10,11], we, therefore, hypothesize second that crystallization may be initiated by reaching a critical concentration of highly stretched chain segments. This would require significantly less energy input than for stretching the entire population of chain segments.
Schematic representation of the gland of a B. mori silkworm, where extensional flow initially aligns at low deformation rates and subsequently stretches the intrinsically disordered silk protein at higher rates. Within the interpretation of a sticky-reptation model, we hypothesize that the experimentally observed gradients in the H serve to control the lifetime of intermolecular cross-links locally within the process, which, in turn, minimizes the energetic requirements to deform the network and induce fiber crystallization.
Schematic representation of the gland of a B. mori silkworm, where extensional flow initially aligns at low deformation rates and subsequently stretches the intrinsically disordered silk protein at higher rates. Within the interpretation of a sticky-reptation model, we hypothesize that the experimentally observed gradients in the H serve to control the lifetime of intermolecular cross-links locally within the process, which, in turn, minimizes the energetic requirements to deform the network and induce fiber crystallization.
To theoretically investigate this hypothesis, we focus our attention on the flow-induced preparation of the conditions for crystallization (rather than crystallization itself). We are, in particular, interested in the specific critical work,
required to induce FIC after a period time during which the system is subjected to the (experimentally controllable) transpose of the (local) velocity-gradient tensor and the (local) stress response . The integral is taken in the (local) Lagrangian co-moving frame of a fluid element. In experimental works (see [35–37] and citations therein), the shear rate and duration render the specific work a control variable () that controls the number of nuclei generated in the system. As the efficiency to converse the energy input into nucleation events is rather limited (estimated [37]), it is worth investigating how the energy loss may be reduced, e.g., by making use of intermolecular associations.
Clearly, the formation of nuclei must be controlled by the underlying molecular conformations. In particular, strong flow aligns and stretches chain segments and thereby reduce the entropic penalty of crystallization, . Hence, the nucleation rate increases as , where one may assume , with being the efficiency of reducing the chain entropy by flow [37]. We argue that in strongly heterogeneous systems, i.e., with broad conformational distributions, it is not the ensemble-averaged change in the entropy [38,39], but the local changes, i.e., the outliers in the tail of the distribution [10,11], that will determine the nucleation rates. Hence, by controlling the distribution, the overall efficiency of the process may be optimized. These broad distributions emerge for sticky polymers in strong flow, which we focus on in the present work; other examples are nonsticky linear polymers in shear [40–42] and ring polymers in extension [43–45].
To calculate the specific work in Eq. (1), one could in principle determine the stress response, , to the transpose of the velocity-gradient tensor, , using MD simulations, where the interactions between stickers are modeled using attractive potentials. At this level of computational detail, sticker dissociation may occur following attempts to escape the attractive potential through molecular vibrations [46,47]. These MD simulations are, however, computationally very demanding, as the dissociation events are quite rare. However, because of this rarity of events, the local equilibration of the chains enables a much simpler description of the chain dynamics in terms of the fraction of closed stickers, and their lifetime, [19]. In a coarse-grained picture, this sticker lifetime is an elementary rather than an emergent time scale. This allows a description of the problem in terms of the dynamics of a single chain in a crowded environment [10,11,48–50], an approach similar to the modeling of entangled polymers through slip-link and slip-spring models [45,48,51–56], where the generation and destruction of entanglements are modeled as elementary processes.
While there is no unique way of formulating a coarse-grained single-chain model [57], all variants of bead-spring, slip-link, and slip-spring models can be written in the general form
where is a chain segment at position that is thermally equilibrated at the relevant time scales [58]. We will refer to this chain segment as a “node” of an elastic network, which may represent a nonsticky segment of a chain (a purely frictional “bead”), a segment with a reversible association (a “sticker”), or it may be an entangled segment (a “slip-link” or a “slip-spring”). Which of these representations is invoked manifests itself in the definition of the friction coefficient, , the (friction-dependent) thermal forces, , and the network forces, . For instance, in classes of models where nodes move affinely with the flow field, the network force exactly cancels the sum of the (conformation-dependent) intramolecular force and the thermal force, . This “rigid-network approximation” is tacitly invoked in the slip-link model by Hua and Schieber [51] and in our recently published model for sticky polymers in a rigid network [10,11]. Within Likhtman’s slip-spring model, the slip-spring may diffuse within a potential energy landscape that represents the elastic compliance of the entangled network [55]. In the present work, we will account for the compliance experienced by the stickers in a reversible network.
In the following, in Sec. II A, we present the usual intramolecular, thermal, and drag forces that act on single chains. To capture how the stickers modify the intermolecular forces (i.e., the “elastic compliance” of the surrounding network) and the segmental drag, we present a nonspatially explicit multichain approach. In Sec. II B, we present a two-state master equation that generates analytical predictions of the impact of sticker opening and closing on both the steady-state and transient stretch distributions of the chains, which enables us to interpret our simulated data in Sec. III. By first mapping the results in the linear flow regime to the analytic sticky-reptation (SR) model, in Sec. III A we discuss how the stochastic nature of sticker opening and closing and the elastic compliance affects the linear rheological data. Then, in Sec. III B, we show how a broad steady-state distribution of chain conformations emerges in strongly nonlinear flows of shear and extension. By simulating the transient emergence of these distributions in start-up flow in Sec. III C, we show that the stickers initially hamper the collective alignments of the chains in mildly nonlinear aligning flows, but facilitates the emergence of stretched outliers. In Sec. III D, we discuss how these outliers may reduce the critical specific work for FIC. In the discussion and conclusions of Sec. IV, we use our findings to interpret the experimental observations of silk spinning and argue that the chemical tuning of associations is indeed a promising mechanism to control the FIC of artificial materials.
II. MODEL AND THEORY
A. Brownian dynamics of sticky polymers in flow
In this section, we will present a coarse-grained description of associating polymers, where the dynamics of sticker opening and closing will depend on the number of open and closed stickers in a nonspatially explicit collection of chains. Any linear polymer that consists of monomers may be discretized using a number of nodes, , see Fig. 2. We use the wording “node” to emphasize that the node may not just represent a traditional, frictional bead of a bead-spring model, but may also represent a sticker that can be in an open or closed state or a slip-link or slip-spring (which, unlike traditional beads, may fluctuate in numbers). Each node is located at a spatial coordinate relative to the center of mass of the chain. The strand between neighboring nodes and has an end-to-end vector and contains a fraction of all the monomers in the chain. At this level of coarse graining, the friction of each node is given by
with being the monomeric friction. The assumption that the dangling chain ends are relaxed may be released by explicitly modeling the position of the chain ends and setting at and at [59].
(color online) The theory in Sec. II A applies to sticky entangled polymers that are parameterized using the locations of nodes. Each node may be a bead (green disk), a sliplink/entanglement (blue ellipses), a closed sticker (orange disk), or an open sticker (orange circles). All nodes are assigned a friction that depends on the fraction of monomers of the chain, , that reside in each of the substrands, see Eq. (3). In general, the number of beads and entanglements may fluctuate during a simulation. In the present work, we focus on the physics of the stickers and fix the number of beads and do not include any entanglements.
(color online) The theory in Sec. II A applies to sticky entangled polymers that are parameterized using the locations of nodes. Each node may be a bead (green disk), a sliplink/entanglement (blue ellipses), a closed sticker (orange disk), or an open sticker (orange circles). All nodes are assigned a friction that depends on the fraction of monomers of the chain, , that reside in each of the substrands, see Eq. (3). In general, the number of beads and entanglements may fluctuate during a simulation. In the present work, we focus on the physics of the stickers and fix the number of beads and do not include any entanglements.
The equilibrium structure of the chain in quiescent conditions is determined by the end-to-end distance of the substrands, , where the stretch ratio obeys the equilibrium distribution
This distribution emerges as a consequence of the intramolecular and thermal forces in Eq. (2).
In order to derive the intramolecular spring forces, we consider the spring force of the entire chain of monomers with a mean stretch ratio of unity
where
approximately captures the anharmonicity of the spring force due to the finite extensibility of the substrand [60]. For the substrands , the harmonic spring force is larger than that of the full chain and the maximum stretch ratio is smaller. This is captured by the renormalization , , and . The direction of the force exerted by spring on node is , while the direction of this force acted upon node is . Hence, the net intramolecular force exerted on node is
The thermal force is given by the equipartition theorem
with being the Cartesian coordinates and being the thermal energy.
The force acted upon the nodes by flow is, provided that our coordinate system moves with the flow field, given by
where is the transpose of the velocity-gradient tensor, which in extension and shear is given by
respectively. As the coordinate system moves with the flow field, the spatial quantities of physical interest to calculate are the deformation of the individual substrands,
using which we recursively obtain the drift of the nodes as
The value of the first entry, , is adjusted to fix the center of mass of the chain (this assumes that the center of mass moves affinely with the flow field).
The dynamics of the chain conformation depends on the state of the stickers through the network force, which, in turn, depends on the dynamics of sticker opening and closing and so, finally, on the chain conformation itself. In particular, when chain segments are highly stretched, the network forces may cause the stickers to dissociate. To obtain these forces, we simulate multiple chains and track the collection of open and closed stickers. When sticker from chain A and sticker from chain B are closed to form a pair, the friction coefficient, the thermal force, and the network force are modified until the sticker pair opens again. The friction coefficient of both nodes becomes , where and are given by Eq. (3), and the thermal forces are given by the equipartition theorem Eq. (10) as before, but with this modified friction coefficient. The network forces are now given by
Hence, the paired stickers and have an identical friction coefficient and experience the same net force (where ). Crucially to forced sticker dissociation, the net force that acts on the closed sticker pair is
which we assume, as in other cases of forces temporary unbinding, lowers the activation energy for sticker dissociation as
with being the activation energy in quiescent conditions and being the typical length scale associated with sticker dissociation [11]. We remark that the (apparent) activation energy obtained from experiments using the Arrhenius-type equation [24] , for the sticker lifetime with an attempt frequency, may be much larger than this activation energy for dissociation. This is due to fast sticker recombination processes [9,61] or due to the mixing of various mechanisms of sticker opening and closing, such as bondswapping [11,62].
For now, we assume a well-defined pairwise association-dissociation reaction whose equilibrium condition is described by the detailed balance , with being the equilibrium constant in the absence of any chain tension. Here, the free energy captures the shift in detailed balance (i.e., the fraction of closed stickers decreases with an increase in chain tension), while in Eq. (17) modifies the rate by which the equilibrium is reached. Indeed, in terms of transition state theory, we may write the opening and closing rates as and , respectively, where , and where is the so-called Brønsted–Evans–Polanyi coefficient [63]. While its value may be determined using experiments or atomistic simulations, we know that must be larger than zero in order to capture strain-induced sticker dissociation [29–33]. We argue that the rheological physics of a reversible polymer network does not necessitate exact knowledge of : When a sticker opens, it may freely diffuse and find conditions to bind to another sticker that is not subject to the influence of strongly stretched chain segments: association will typically take place in conditions where the activation barrier is equal to that in quiescent conditions. Indeed, in our simulations, we find that the mean fraction of open stickers in conditions of strong flow remains similar to the fraction in quiescent conditions, despite noticable acceleration of sticker dissociation.
These arguments have enabled us to conveniently set and ; the latter avoids the need for on-the-fly calculations of association rates during our simulation. We have implemented the opening and closing of stickers using a kinetic Monte Carlo (kMC; also known as a Discrete Event Simulation) scheme, where after a time interval a sticker is opened or closed with a probability or , respectively. In our simulation algorithm, shown in Fig. 3 and discussed in detail in Subsection 1 of the Appendix, we take time steps during which the chain conformations are approximately fixed and for which the time-independent (but conformation-dependent) rates of sticker opening and closing are calculated. The dynamics of the stickers is simulated during the time step using a kMC scheme. This essentially creates and destroys constraints in a similar way as in the slip-link model [54], but where the constraints physically represent closed stickers instead of entanglements (hence, our approach may be generalized using appropriate kMC algorithms [64–66] to go beyond the unentangled chains with pairwise association and dissociation of stickers focused on in the present work, and also capture entanglements, stickers that dimerize through bondswapping, and stickers that may assemble into larger aggregates). After this step of “constraint-dynamics,” the Brownian dynamics are solved, the conformations are updated, and the next time step is commenced.
Flow chart of the algorithm to simulate the conformational dynamics of sticky polymers and the dynamics of sticker association and dissociation (detailed discussion: see Subsection 1 of the Appendix).
Flow chart of the algorithm to simulate the conformational dynamics of sticky polymers and the dynamics of sticker association and dissociation (detailed discussion: see Subsection 1 of the Appendix).
B. Approximate theory in transient extensional flow: Two-state model
The dynamics of sticky polymers is complicated by the fact that a polymer with stickers can be in different states, as each individual sticker can be either opened or closed. An instructive simple case is a chain with , as the chain is either completely free to relax when either of the stickers is open (state 1) or can only be extended by flow when both stickers are closed (state 0). Hence, we can accurately distinguish between an extension state where the polymer is unable to relax and a relaxation state where the polymer is able to relax. Using this “two-state” description, we previously discovered that stickers give rise to enormous stretch fluctuations in extensional flow below the strain rate at which the mean stretch diverges, i.e., below the “stretch transition,” which are described by the steady-state power-law stretching distribution [10],
It turned out that this two-state prediction, which is exact for chains with two stickers, also described the steady-state stretch distribution for chains with multiple stickers. In the present work, we recapitulate our previous analysis of the steady-state situation and extend it for transient start-up flow. In all of this analysis, we will consider a single relaxation mode of the polymer at time scales beyond the relaxation time of the surrounding network; hence, we invoke the rigid-network approximation in this entire section.
The starting point is to consider a chain in two states where the chain is either unable to retract (state 0) or is free to retract (state 1). The opening rate is , and the closing rate is . The time development of the probability distribution of the stretch ratio is described by [10]
with being the bare Rouse time of the chain without stickers. In this equation, we have neglected the high-frequency relaxation modes of the polymer as well as the (potentially much slower) relaxation of the surrounding network; the latter is justified in view that the network rapidly stiffens with an increasing strain. To approximate this equation analytically, we first make the substitution , so . Similarly, . Inserting this into the governing equations gives
The nonlinear contributions can then be omitted by considering the limit of large stretches where their contribution to the distribution is exponentially small, i.e., we approximate , which is equivalent to .
In steady state, the left-hand side of the equation is zero and the equations can be cast in the form , with and a constant by matrix. The solution of this system of first-ordinary differential equations is given by [10]
with being the normalization constant (its value can in principle be determined by releasing the approximation ), and with the exponent of the power-law distribution given in terms of physical parameters by
[this is one of the eigenvalues of Eqs. (21) and (22); the other eigenvalue is and is unphysical as a distribution of the form cannot be normalized.] The value of this stretching exponent diverges if the bare stretch transition at is approached from small strain rates. However, because of the physics of the stickers, actual divergence already occurs at lower strain rates: At , the exponent becomes and the stretch distribution can no longer be normalized. Depending on the sticker lifetime, at smaller strain rates the exponent may reach a value if the “sticky Weissenberg number” reaches unity; here, the mean stretch diverges. While the mean stretch is finite for smaller strain rates, the variance of the stretch diverges for , which happens if becomes larger than [10], at which point (considerably slower than the bare stretch transition) we expect a long tail of very high stretched chains to develop in the distribution.
This analytic approach can be extended to predict the transient dynamics of the distribution in start-up flow. As we will show, the late-stage dynamics in which the tail of the distribution “fills up” is independent of the initial conditions. In those late stages, the distribution reaches a steady state for stretches below a certain “front,” (above which the distribution function has a value of zero), which shifts to high stretch values over time. The precise number of chains with a certain stretch also depends on the width of this moving front. We assess analytical predictions on the front position and width using the two-state model using solutions in an early- and late-stage regime, where the time scale is, respectively, much shorter and much larger than the sticker lifetime. While the long-time regime will slow down the progression of the front due to sticker opening, in the early-stage regime we will obtain an upper limit of the rate by which the front moves.
In the early-stage regime, we approximate the stretch distribution using the Dirac-delta distribution (justified by the very wide long-time distribution), at ), from which it can be easily seen that the distributions shift initially, when pure advection dominates over sticker dynamics, to higher stretches for the closed state, and retract to smaller stretches for the open state . This suggests that the “front,” , of any distribution with finite , shifts exponentially in time to higher values through .
To develop an analytic approximation for the long-time limiting behavior of the sticky polymers in start-up flow, we consider some point in time where sufficient stickers have opened to facilitate chain relaxation and assume that the stretch distribution has reached a steady-state for small stretches but is empty for larger stretch ratios. Here, can be thought of as the establishment of the “front” of the stretch distribution at later times moving to higher stretches. In the following, we will show that the ansatz of this moving front is indeed a good approximation for the tail of the transient stretch distribution and that for later times , further convergence of the stretch distribution takes place in the range of stretches , where the “front” of the distribution shifts to high stretch values as . Assuming that , the steady-state portion of the distribution is negligibly affected by the loss of small-stretch contributions to the tail of the distribution [see discussion around Eq. (A11) in Subsection 2 of the Appendix], and for any time the portion of the stretch distribution becomes independent of time beyond . The constancy of the distribution at provides a fixed-boundary condition. Hence, this problem essentially models the dynamical response to a unit step and lends itself to an analysis through a Laplace transform to give a solution for the distribution at each stretch ratio of the form , which is the Laplace transform of a time-dependent function that becomes nonzero at the time . The inverse function is then the trajectory of the “front” of the distribution. In Subsection 2 of the Appendix, we detail the Laplace transform of Eqs. (21) and (22) with the boundary condition in this long-time regime, which as a solution gives
with being the “steady-state stretch exponent” in Eq. (25) and with
being the “dynamic stretch exponent,” which controls the growth of the front of the distribution as
In this equation, and are the (extensional) Weissenberg numbers of the chain without and with stickers, respectively; within the two-state model, , see discussion under Eq. (25). Upon approaching the stretch transition where the mean stretch diverges, indicates “critical slowing down,” as the (late-stage) front of the distribution becomes immobile. For chains with strong stickers at the strain rate where the variance of the stretch diverges [see discussion under Eq. (25)], we find , which indicates that the late-stage measure of the front is shifted from the early-stage measure for the outliers by a factor . We have also checked that the moving front is narrow for small strain rates . In Subsection 2 of the Appendix, we provide more analytical analysis of the two-state model to estimate the width of the front (relative to its extent) as , where . As we shown in Subsection 2 of the Appendix, typically this width is , and the front of the distribution is narrow even close to the stretch transition.
III. RESULTS
A. Linear dynamics
We have verified the physics of our model in the linear viscoelastic regime by first simulating nonsticky chains of fixed length but a varying number of beads from to (the beads are regularly along the backbone of the polymer, so for all ). Figure 4 shows that the choice of the number of beads has a negligible influence on the time evolution of the mean-square displacement, , of the center of mass and is in all cases in agreement with the theoretical prediction
where the diffusivity, , is for nonsticky polymers given by the bare Rouse diffusivity
Moreover, the inset of Fig. 4 shows that also the end-to-end-distance, , is distributed according to the physical equilibrium result of Eq. (4).
of the center of mass of a nonsticky polymer against time (the main panel) and the time-averaged end-to-end length () distribution (the inset). The number of real monomers per chain is fixed, while the level of coarse-graining is varied through varying the number of beads, , per chain. The symbols and solid black curves represent the simulations and the theory, respectively.
of the center of mass of a nonsticky polymer against time (the main panel) and the time-averaged end-to-end length () distribution (the inset). The number of real monomers per chain is fixed, while the level of coarse-graining is varied through varying the number of beads, , per chain. The symbols and solid black curves represent the simulations and the theory, respectively.
For times shorter than the Rouse time of strands between stickers, i.e., for , the dynamics of a sticky polymer are governed by the same Rouse diffusion as nonsticky chains, see Fig. 5(a). For later times than that, the motion of the polymer is subdiffusive until the sticky-Rouse time , which is approximately given by [19]
Focusing on the crossover from early-stage bare Rouse diffusion to subdiffusive motion, one would expect this crossover to occur at the point in time where the substrands between stickers have just relaxed and where further relaxation requires sticker dissociation. Indeed, we find this is the case within the rigid-network approximation. However, for the elastically compliant network, the closed stickers themselves are able to diffuse. The friction experienced by the closed sticker depends on the level of deformation of the surrounding network, which is initially small. As the sticker diffuses further, a larger portion of the surrounding network is deformed and the contribution of “next-neighbor” stickers starts to contribute to the friction. Clearly, the increase of the friction increases rapidly beyond a certain characteristic distance. It is unknown what this distance might be, but it is likely to be strongly dependent on the topology of the network. The plateau value in Fig. 5(a) shows that for our simulations this happens to occur when the MSD of the center of mass of chain is approximately , i.e., when the center of mass of the chain has diffused – times its end-to-end distance.
Linear rheology of a sticky chain with , , within the rigid-network approximation (open symbols) and with this approximation released (closed symbols). (a) of the center of mass against time. (b) Storage, , and loss, , modulus in units of against the frequency, , plotted for the chain in (a) as well as for an nonsticky chain (triangles). There is fair agreement with the analytical sticky-Rouse model in Eq. (33) (solid curves) for the sticky chain within the rigid-network approximation for the nonsticky chain. For the sticky chain with an elastically compliant network, the plateau modulus decreases to that of the theory with (dashed curves).
Linear rheology of a sticky chain with , , within the rigid-network approximation (open symbols) and with this approximation released (closed symbols). (a) of the center of mass against time. (b) Storage, , and loss, , modulus in units of against the frequency, , plotted for the chain in (a) as well as for an nonsticky chain (triangles). There is fair agreement with the analytical sticky-Rouse model in Eq. (33) (solid curves) for the sticky chain within the rigid-network approximation for the nonsticky chain. For the sticky chain with an elastically compliant network, the plateau modulus decreases to that of the theory with (dashed curves).
The elastic compliance not only affects the subdiffusive motion of the chain, but also the sticky-Rouse diffusivity at times beyond the sticky-Rouse time. While the analytical expression for the sticky-Rouse diffusivity accurately describes our simulations within the rigid-network approximation, we find that it overestimates the diffusivity of chains in an elastically compliant network. We have investigated the consequence of this to the interpretation of linear viscoelastic data, which are often used experimentally to estimate the number of associations per chain, by calculating the dynamic moduli and against the frequency in Fig. 5(b). The data shown include nonassociating unentangled chains () and the unentangled sticky chains of Fig. 5(a); i.e., chains with stickers within the rigid-network approximation and with an elastically compliant network. The simulated data (symbols) were obtained from the relaxation modulus, , through the multiple-tau-correlator algorithm discussed in [67]. To obtain the dynamic moduli and , we have used the finite-element approach from [68]. We have compared the data to the sticky-Rouse model (curves), which is given by
In this equation, the first summation captures the high-frequency bare Rouse modes (the number of Kuhn segments, , truncates the highest frequencies) and the second summation captures the sticky-Rouse modes. The modulus is proportional to the number density of monomers and to the thermal energy.
Figure 5(b) shows dominance of bare Rouse relaxation at high frequencies, where all moduli will approach (in principle) the scaling relation . Discrepancies, such as a roll-off of at high frequencies, emerge due to the finite number of modes/beads that are included in the simulations. At decreasing frequencies, the moduli of the nonsticky chains (triangles) decrease rapidly, while the moduli of the sticky chains reach a plateau value that ranges down to . Within the rigid-network approximation (closed circles), the modulus of the plateau is in agreement with the sticky-Rouse model in Eq. (33) for . However, if the network is elastically compliant (open circles), the plateau value decreases and is better described if the theory would be adjusted with an apparent number of stickers (dashed curves). At lower frequencies , the moduli rapidly decrease. In the simulations, the moduli decrease much more rapidly than in the theory, as also noted earlier in [50]. We find that this terminal relaxation time (we remind the reader that this relaxation time is for unentangled chains entirely determined by sticker relaxation, i.e., not by SR [22,24]) is even further reduced for the chain in an elastically compliant network. Consequently, the peak of the dynamic modulus is much narrower than in the theory. We have estimated that the shape of this peak is best described by within the rigid-network approximation and for the compliant network. This clearly indicates that analysis of the dynamic modulus peak in rheological data (which is required when high frequencies are experimentally inaccessible [9]) provides an underestimate of the actual number of stickers per chain.
To obtain a wider view of the impact of the elastic compliance on the dynamics of chains with various numbers of stickers and sticker lifetimes, we have calculated the diffusivities of various chains within the rigid-network approximation and with a compliant network in Fig. 6. Panel (a) shows that the predictions of [19] describe our simulations well within the rigid-network approximation for chains with stickers with various sticker lifetimes, in particular, in the regime where the sticky-Rouse diffusivity scales with the sticker lifetime as , see Eq. (32). Panel (b) shows that upon releasing the rigid-network approximation this scaling behavior persists, but rescaled with a prefactor . While this scaling regime is reached for the chains with more than five stickers (i.e., above the percolation threshold for network formation), this is not the case for the chains with two stickers. Within the rigid-network approximation, this originates from the fact that at sticker lifetimes a plateau is reached where the chains with all stickers open dominate the dynamics. Without the rigid-network approximation, the chains cluster into linear “supramolecular” dimers, trimers, etc. through an exponentially decaying cluster-size distribution [69], which implies a distribution of diffusivities that strongly differs from that predicted by the sticky-Rouse model. Hence, while our simulation approach accounts for the elastic compliance of the percolating network, it also captures the contributions of cluster diffusion near and below the percolation threshold for network formation.
Sticky-Rouse diffusivity, , against the sticker lifetime, for chains with stickers with within a rigid network (a) and a compliant one (b). The symbols are our simulation results, and the curves represent the sticky-Rouse model in [19]. The units are given in terms of the bare Rouse diffusivity and the bare Rouse time, .
Sticky-Rouse diffusivity, , against the sticker lifetime, for chains with stickers with within a rigid network (a) and a compliant one (b). The symbols are our simulation results, and the curves represent the sticky-Rouse model in [19]. The units are given in terms of the bare Rouse diffusivity and the bare Rouse time, .
B. Nonlinear dynamics: Steady state
Ordinary Gaussian polymer melts and solutions of narrow molecular-weight distribution exhibit broad conformational distributions in shear flow due to dynamic stretching, tumbling, and recoiling of the chains [40–42]. In extensional flow, however, such chains do not tumble and recoil, and their stretch distributions are narrow, see Fig. 7(a). Perhaps surprisingly, by incorporating stickers into the chain these stretch distributions become much wider, see Fig. 7(b). This figure shows that the sticky chains exhibit an enormous dispersity in the chain stretch, as well as occasional hairpin conformations [Fig. 7(b)]. These are caused by the stochastic binding and unbinding of stickers, where the network forces may occasionally act in the opposite direction of the drag forces exerted by flow.
Representation of simulated chain conformations in extensional flow for for nonsticky (a) and sticky (b) polymers. While the variations in stretch are narrow for nonsticky polymers, these variations are broad for the sticky polymers: when a sticker in a retracting chain segment binds to a neighboring chain segment, this may disrupt the neighboring chain. The scale bar represents approximately a length , which is of the fully extended chain.
Representation of simulated chain conformations in extensional flow for for nonsticky (a) and sticky (b) polymers. While the variations in stretch are narrow for nonsticky polymers, these variations are broad for the sticky polymers: when a sticker in a retracting chain segment binds to a neighboring chain segment, this may disrupt the neighboring chain. The scale bar represents approximately a length , which is of the fully extended chain.
To go beyond these qualitative observations, we have quantified this phenomenon using steady-state stretch distributions of polymers at various extension and shear rates in Fig. 8. We have selected nonsticky polymers () and sticky polymers below () and above () the percolation threshold for network formation: the chains with only two stickers may only assemble into high-molecular weight chains, while chains with five stickers may branch into percolating networks. We have modeled the physics of the stickers using the same description as in our previous work on chains that are prealigned in the flow field [11]. We have summarized the associated parametrization in the caption of Table I. In extensional flow, above the sticky Weissenberg number, , with being the sticky-Rouse time we expect divergent stretching (albeit that real divergence is obstructed by the maximum chain extensibility ). We have calculated the sticky-Rouse time as , with the ratio between the sticky and the bare diffusivity as presented above in Fig. 6. The relevant results are summarized in Table I.
Simulated steady-state stretch distributions of the end-to-end distance, , for various extension [(a), (c), and (e)] and shear [(b), (d), and (f)] rates for a linear unentangled, nonsticky (), and sticky ( and ) polymers. For these simulations, (see Table I for all parameter values). The black curve represents the contour-length fluctuations in quiescent conditions, given by Eq. (4).
Simulated steady-state stretch distributions of the end-to-end distance, , for various extension [(a), (c), and (e)] and shear [(b), (d), and (f)] rates for a linear unentangled, nonsticky (), and sticky ( and ) polymers. For these simulations, (see Table I for all parameter values). The black curve represents the contour-length fluctuations in quiescent conditions, given by Eq. (4).
In our simulations of sticky polymers in nonlinear flow conditions, we use as parameters p = 0.9 as the fraction of closed stickers (in quiescent conditions), a sticker lifetime τs = 10τR, an activation energy Eact = 8kBT, and a sticker dissociation length of ℓ = 1 nm. The maximum extension ratio of the chain is λmax = 75. The intramolecular forces in Eq. (5) are calculated by assuming a total number of N = 5525 Kuhn segments and a Kuhn length of b = 0.4 nm. As we focus on chains with Zs = 2 and 5 stickers, we here tabulate the ratio between the bare Rouse and sticky-Rouse diffusivities, [DR/DSR], and relaxation times, [τR/τSR]. The diffusivities were determined in Fig. 6, and the sticky-Rouse time is calculated as τSR = [DR/DSR]τR [19].
Polymer model . | DSR/DR . | τSR/τR . |
---|---|---|
Zs = 2 (rigid) | 0.0949 ± 0.0002 | 10.54 ± 0.02 |
Zs = 5 (rigid) | 0.021 56 ± 0.000 04 | 46.38 ± 0.09 |
Zs = 2 (compliant) | 0.4331 ± 0.001 | 2.309 ± 0.005 |
Zs = 5 (compliant) | 0.1050 ± 0.0002 | 9.52 ± 0.02 |
Polymer model . | DSR/DR . | τSR/τR . |
---|---|---|
Zs = 2 (rigid) | 0.0949 ± 0.0002 | 10.54 ± 0.02 |
Zs = 5 (rigid) | 0.021 56 ± 0.000 04 | 46.38 ± 0.09 |
Zs = 2 (compliant) | 0.4331 ± 0.001 | 2.309 ± 0.005 |
Zs = 5 (compliant) | 0.1050 ± 0.0002 | 9.52 ± 0.02 |
Equation (4) shows that in all cases the equilibrium stretch distribution for zero-flow conditions (black curve) is approached for small strain rates. For nonsticky chains (), a broad stretch distribution with a cutoff set by emerges in shear due to the dynamic stretching, tumbling, and recollapsing of the chains. In extensional flow, the distribution broadens only within a narrow range of strain rates around the bare stretch transition, . Beyond the stretch transition, the stretch distribution is narrow and Gaussian and approaches with an increasing strain rate. This behavior qualitatively changes upon incorporating stickers.
Figure 8 shows that the steady-state stretch distributions in shear are similar to those of the nonsticky chains, while in extensional flow the distributions of sticky polymers are remarkably distinct from the nonsticky ones. In contrast to the nonsticky polymers, the sticky polymers show broad stretch distributions in steady-state extensional flow over a broad range of flow rates. We have observed this behavior previously in simulations where the chains were prealigned in the flow field and where we invoked the rigid-network approximation [11]. Our current simulations show that this phenomenon persists when these approximations are released, but also show a dynamic coexistence of stretched chains, relaxed coils, and hairpins. Interestingly, there is a qualitative similarity between the distributions of the chains with two or five stickers, despite the fact that these are below and above the percolation threshold for network formation, respectively. This indicates that the enormous reduction of the chain retraction rate due to the stickers does not necessitate network formation: the formation of high-molecular weight assemblies suffices.
We also find that the large fluctuations in stretch below the formal stretch transition carry over from case of two stickers per chain to multiple stickers [11]. (The stretch transition is defined at the condition , with the sticky-Rouse time obtained from the sticky-Rouse diffusivity of Fig. 6 as .) In particular, we find that for small strain rates and large stretch ratios the stretch distribution has a power-law tail [see Eq. (18)] of which the width is set by a -dependent stretch exponent (see Sec. II B). We have determined the stretch exponent from the distributions of the chains with two and five stickers (we discuss the numerical method in Subsection 3 of the Appendix) in extensional flow with and without the rigid-network approximation and finite extensibility and plot these against the strain rate in Fig. 9. As anticipated, we have been able to map the stretch exponent of the chain with two stickers onto the analytical result in Eq. (25). To achieve that, it has to be taken into account that the open state of the chain can be achieved by opening either of the stickers; hence, in Eq. (25), which models the simultaneous opening of all stickers, is replaced by , and results in
For chains with multiple stickers, no such analytic theory is yet available; however, we do find a qualitative agreement of the increasing power-law exponent with an increasing strain rate.
(color online) Stretch exponent of the power-law tail of the stretch distribution for simulations of polymers with (circles) and stickers (squares), within the rigid-network approximation (closed symbols) and using elastic compliance and finite chain extensibility (open symbols). The solid curve is given by the two-state model in Eq. (34) with (see Table I for all physical parameter values). For (horizontal line), the fluctuations in stretch diverge; this leads to a cutoff in the stretch distribution for chains with finite extensibility, see Fig. 8.
(color online) Stretch exponent of the power-law tail of the stretch distribution for simulations of polymers with (circles) and stickers (squares), within the rigid-network approximation (closed symbols) and using elastic compliance and finite chain extensibility (open symbols). The solid curve is given by the two-state model in Eq. (34) with (see Table I for all physical parameter values). For (horizontal line), the fluctuations in stretch diverge; this leads to a cutoff in the stretch distribution for chains with finite extensibility, see Fig. 8.
For the chains with two and five stickers and with a fraction of closed stickers, we also simulated the stretch distributions while including finite extensibility and an elastically compliant network. Finite extensibility implies that there is a cutoff of the power-law tail, which becomes apparent with increasing (less negative) . Since the fluctuations in diverge for , this cutoff has a significant effect on the tail of the stretch distribution upon approaching . Figure 9 does confirm a broadening power-law stretch distribution for the chains in a compliant network, but shifted to higher strain rates, as expected from the faster sticky-diffusion rates from Fig. 4.
C. Nonlinear dynamics: Transients
In our pursuit to understand the FIC of associating polymers such as the silk protein, we are interested in capturing the macroscopically observable stresses in start-up flow, and to interpret crystallization rates in terms of the chain conformations that underlie these stresses. To address these challenges, in this section, we will present the time-dependent rate-normalized transient shear stress, , and extensional stress , with the stress tensor (in units of energy per molecule) given by
Focusing first on the results for nonsticky chains with a finite extensibility in Figs. 10(a) and 10(b), we reproduce the well-known qualitative features of their stress transient [57]. For small Weissenberg numbers, , the polymers are able to relax, while for large strain rates there is an overshoot in shear flow, which is related to the onset of tumbling and recollapsing of stretched chains, and in extensional flow there is a sharp increase in the extensional stress until a plateau due to the finite extensibility of the chains is reached. Because of the thermal fluctuations and dispersity in the initial chain conformations, Fig. 10(b) shows broadening of the stretch distribution at early times. At late times, when all chains are aligned (at the level of the beads), a sharp peak emerges at high stretches near the maximum extensibility .
(a) and (c) Simulated rate-normalized transient extensional and shear stresses averaged over 50 polymers for the nonsticky (a) and the sticky (c) case. The sticky polymer exhibits strong fluctuations for , which is below the stretch transition (at , see Table I). (b) and (d) Transient stretch distribution of the end-to-end distance, , in extensional flow for the nonsticky (b) and sticky (d) chain at selected strain rates. The error bars in (d) represent half of the standard error of the mean. All physical parameter values are given in Table I.
(a) and (c) Simulated rate-normalized transient extensional and shear stresses averaged over 50 polymers for the nonsticky (a) and the sticky (c) case. The sticky polymer exhibits strong fluctuations for , which is below the stretch transition (at , see Table I). (b) and (d) Transient stretch distribution of the end-to-end distance, , in extensional flow for the nonsticky (b) and sticky (d) chain at selected strain rates. The error bars in (d) represent half of the standard error of the mean. All physical parameter values are given in Table I.
This sharp peak in the stretch distribution is a fingerprint for nonsticky linear polymers in extensional flow and will not be visible for the sticky polymers, as we will now show for . We plot the resulting start-up stresses and stretch distributions in Figs. 10(c) and 10(d).
Qualitatively, we find similar shear and extensional viscosities as in the nonsticky case, although there is now no distinctive overshoot in shear flow. In extensional flow, the stresses at long time scales have shifted to higher values because of the contribution by the reversible cross-links. Further, while nonsticky polymers show strain hardening only for , the sticky ones also show strain hardening for smaller strain rates . For strain rates smaller than that we identify large fluctuations in the transient extensional stress, which are caused by temporary exponential stretching of chain segments between closed stickers that rapidly retract to a near-relaxed state when the stickers open [10]. For strain rates , these fluctuations fill up a power-law distribution whose stretch exponent is depicted in Fig. 9. For higher rates, the finite extensibility causes a truncation of this power law tail.
The dynamics by which the stretch distributions evolve in extensional flow above the stretch transition () is shown in Fig. 10(d). At early times, the stretch distribution closely resembles the equilibrium distribution of Eq. (4). As time proceeds, the distribution broadens exponentially with time as until the steady state is reached after a time . This is in qualitative agreement with the predictions of the two-state model that we derived in Eq. (29) of Sec. II B.
D. Critical specific work
Now that we have captured how stickers lead to broad stretch distributions, we will investigate how these distributions affect the critical work for FIC. The usual predictor for FIC is the “Kuhn segment nematic order parameter,” . If (see, e.g., [3]), virtually all chains are aligned at the level of the Kuhn segments, i.e., they are completely extended/stretched in the direction of the flow field. However, in this case of high chain heterogeneity, we expect this average measure to be a poor descriptor. We know that the critical nuclei will be dominated by the small fraction of highly stretched chains and that it is the oriented segments in these chains only that promote crystallization. To model this extremum-dominated physics, therefore, we will assume that FIC may commence when a critical fraction, , of chain segments of some length have stretched beyond a critical stretch ratio , where is the maximum stretch of the chain segment and a parameter that may be viewed as proxy for chain stretch at the Kuhn length of this extremely stretched chain fraction. Hence, the criterion for FIC may within our interpretation be formulated as
where is the transient stretch distribution function and is the time into the process of startup flow at which the criterion is satisfied. Essentially, this criterion provides a prediction for the time required to form the first nuclei and, hence, the time should not be confused with the fixed time in FIC experiments [35–37] during which a different number of nuclei may form depending on the strain rate. A comparison to those experiments would require knowledge of the physical relationship between the nucleation rate and the conformational distribution; here, we have proposed a hypothetical condition that is likely to correlate to a fixed nucleation rate. For associating polymers, a natural measure for the length of flow-crystallizable chain segments is ; in general, however, measures for , , and will have to be determined through experimentation and (atomistic) MD simulations [15–18].
In this section, we will employ simulations with chains of a fixed number of beads (i.e., with chain segments, giving ), and we will monitor the maximum stretch among the total of 500 chain segments (i.e., ). The time-evolution of the maximum stretch will enable us to screen how various values of require a different processing time and a different input of specific energy. We obtain statistics on this relationship by averaging our results over five simulations with different initialization “seeds” of the random-number generator. We will discuss the implications of the criterion in Eq. (36) by comparing it to a measure of the (mean-field-type) nematic order parameter. At our level of coarse graining, the highest resolution of nematic chain alignment is captured using the nematic order parameter , which is the largest eigenvalue of the nematic order tensor , where is the unit vector tangential to the backbone of the chain. (We remark that this nematic order parameter is an overestimate of the Kuhn segment nematic order, i.e., .) In Fig. 11, we have calculated the critical specific work, , as given in Eq. (1), needed to achieve values of and in the range from 0 to 1 for nonsticky () and sticky () chains for various shear and extensional rates.
(color online) Nematic order parameter, and characteristic stretch ratio, , against the specific work (see the main text) for sticky (closed symbols) and nonsticky (open symbols) polymers in shear (left) and extensional (right) flow. The symbols are obtained from simulations with various strain rates for a chain with with an elastically compliant network. All physical parameter values are given in Table I.
(color online) Nematic order parameter, and characteristic stretch ratio, , against the specific work (see the main text) for sticky (closed symbols) and nonsticky (open symbols) polymers in shear (left) and extensional (right) flow. The symbols are obtained from simulations with various strain rates for a chain with with an elastically compliant network. All physical parameter values are given in Table I.
The top panels of this figure give the nematic order parameter, , and the measure for stretch fluctuations, against the critical specific work. For large values of the critical work, both measures converge, which suggests that both measures can interchangeably used as predictors for FIC for nonsticky chains. We notice that the critical work in shear (left) and extensional flow (right) shows similar trends well above the stretch transition (the stretch transition of the bare chain is ). Just above this transition, the critical work required is relatively large. This implies a monotonically decreasing critical work with an increasing strain rate, which is due to the suppression of energy dissipation by recoiling of the chains (we discuss this in more detail in Fig. 12). This is in contrast to the typical behavior in experiments on nonassociating polymers (e.g., the FIC of HDPE [7]), where the critical work increases with an increasing strain rate. We argue this discrepancy occurs because we here consider unentangled rather than entangled chains. Finally, the top panels of Fig. 11 confirm the expected behavior where the nematic order parameter (closed symbols) is typically larger than the stretching parameter (open symbols): with an increasing specific work, the chains first align and then stretch.
The critical time (left) and the specific critical work (right) against the sticky Weissenberg number, , for various and criteria for the critical condition. The open symbols were calculated in shear and the closed ones in extensional flow. The values are obtained for a chain with with an elastically compliant network. It is useful to interpret the strain rates in relation to the stretch transition for the sticky chains in extension at , where the “sticky” Weissenberg number is , with being the Weissenberg number of the nonsticky chain. This factor is nonuniversal and depends on the number and lifetime of stickers, see Table I for all physical parameter values. The solid curves are given by Eq. (40) for and for .
The critical time (left) and the specific critical work (right) against the sticky Weissenberg number, , for various and criteria for the critical condition. The open symbols were calculated in shear and the closed ones in extensional flow. The values are obtained for a chain with with an elastically compliant network. It is useful to interpret the strain rates in relation to the stretch transition for the sticky chains in extension at , where the “sticky” Weissenberg number is , with being the Weissenberg number of the nonsticky chain. This factor is nonuniversal and depends on the number and lifetime of stickers, see Table I for all physical parameter values. The solid curves are given by Eq. (40) for and for .
This behavior is crucially altered for the sticky polymers, as shown in the bottom panels of Fig. 11. We find that the alignment of the chains requires more critical work in both shear (left) and extensional flow (right), which is due to the fact that the full alignment of the chains requires the opening of intermolecular associations. On the other hand, the stretching of chain segments can take place before global chain alignment. (Note that the stretch transition is for this system, see Table I.) The stretching parameter (open symbols) follows a sharp sigmoidal dependence against the critical work and rapidly outgrows the alignment parameter (closed symbols). This is possible because the stretching parameter provides information about a fraction of chains in the tail of the distribution, while the alignment parameter provides information about the mean properties. This supports our hypothesis that FIC may be achieved at a small critical specific work by exploiting the stochastic nature of associating polymers.
Given either a or criterion for critical nucleation, we are interested how the strain rate affects how much critical specific work, , is needed, and at what time scale, this criterion is achieved. To investigate this, we focus on horizontal lines/cross sections of Fig. 12 (i.e., at fixed values and of both and ). For the data points along these lines, we plot the critical work, , and the time scale, , in Fig. 12. The left panel shows that the time scale scales as , as one may expect and discuss in more detail below. Below the stretch transition, this dependence becomes stronger: under these conditions many chain stretches are attempted, but fail due to sticker opening and lead to energy dissipation through chain retraction. This crossover between two regimes qualitatively agrees with that found in Fig. 2 of the work by Holland et al. on silk [7]; more dedicated research is needed to investigate this observation.
The right panel of Fig. 12 shows the critical specific work needed to achieve a certain degree of alignment, , or of stretch fluctuations, , in shear (open symbols) and extensional flow (closed symbols), against the sticky Weissenberg number. Evidently, a high degree of overall alignment/nematic order requires much larger specific work than a small fraction of large stretch fluctuations does, as discussed in Fig. 11. Having in mind our overarching proposition that crystallization may occur in response to stretch fluctuations, we now focus on the measure for . We remark that for the system we studied, the stretch transition in the absence of stickers takes is located at (because , see Table I). For smaller strain rates, , we find there is a minimum in the specific critical work near the stretch transition . Indeed, while large stretches are achieved just below the stretch transition due to long power-law tails in the stretch distribution [10], many attempt fluctuations are needed before the required stretch value is achieved. Due to the energy dissipation of such unsuccessful attempts, the specific critical work increases with a decrease in strain rates. Above the minimum, the specific work increases and eventually reaches a plateau.
We explain the increase in the critical specific work with an increasing strain rate in terms of the two-state model that we introduced in Sec. II. We argue that the stress is dominated by the contributions of stretched chains in the closed state,
with being the constant, assuming that the open chains are in a relaxed state. Here, is the stretch distribution of the closed chains, of which we will discuss the dynamics below. We will then calculate the critical specific work as . To calculate , we first will determine using the criterion
which, as before, implies a minimum concentration of chains with a stretch ratio of at least . Second, we will need an expression for the time evolution of the probability density .
To obtain , we will assume that all chains that have (temporarily) opened are sufficiently relaxed compared to the most stretched chains to have a negligible contribution to the overall stress . Therefore, we will only take into account the loss of strongly stretched chains by opening rate and ignore the contribution of closing events by rate . We will further use the initial condition , with being the Dirac delta distribution to represent a narrow stretch distribution at time . The dynamical equation in Eq. (21) then predicts that the Dirac delta distribution in time shifts to high stretch values along the axis, as
with an amplitude that decreases in time due to sticker opening (we present the derivation in the first two paragraphs of Subsection 2 of the Appendix).
Equation (39) shows that the critical stretch and the critical time are related by , which is in agreement with our simulated results displayed in Fig. 11. We insert this equation into the expression for the critical specific work, , and find
where is the minimum strain rate for which the criterion in Eq. (38) is obeyed. This function is plotted in Fig. 12(b). It diverges at (this divergence is not followed by the simulation data, because stochastic closing events that generate new bound chain segments), reaches a minimum, and then monotonically increases toward a plateau value. Physically, this plateau value represents the case where the entire distribution of chains is stretched to reach the critical stretch value . In this case, the concentration of stretched segments far exceeds the critical concentration, and more energy has been put into the system then needed. By decreasing the strain rate, an increasing number of stickers are able to open and the stress is relaxed, in turn decreasing the critical specific work to achieve the critical condition in Eq. (38). This supports our proposition that the stochastic nature of the binding and unbinding of associations enables molecularly engineer associating polymers to undergo FIC at low energetic costs. In particular, we have shown, using simulations and an approximate theory in Eq. (40), that there is an optimum strain rate at which the critical work for critical stretch is minimized.
IV. DISCUSSION AND CONCLUSIONS
This work has shown that the transient evolution of the chain-stretch distribution of associating “sticky” polymers in shear, and especially extensional, flow possesses an extremely rich structure. The theoretical and numerical investigations reported here were driven by the observation that the silk protein (i) undergoes efficient, chemically tunable, FIC and (ii) can be modeled as an associating/sticky polymer. Our findings have implications for the interpretation of silk-spinning data, as well as to the development of novel associating polymers and the computational modeling tools (we introduced a “sticky” sliplink model and an analytical two-state master equation, which may be transferable to also address the peculiar dynamics of ring polymer in flow [43–45]).
Regarding silk rheology, we have theoretically confirmed our hypothesis that the stickers between chains may reduce the critical specific work to induce FIC under reasonable assumptions for critical crystallization criteria. In our approach, we have adopted the view that FIC may commence when a sufficient concentration of chains is aligned at the level of the Kuhn segments. However, in contrast to the ensemble-averaged approach where the Kuhn segmental nematic order parameter is measured as a predictor for FIC, we have assumed that a critical concentration of strongly stretched chain segments in the tail of the distribution is a sufficient condition for crystallization. Indeed, by comparing a measure for the stretch fluctuations to the (ensemble-average) nematic order parameter, we have found that the stickers hamper initial chain alignment (chain alignment is slowed down by the need for stickers to dissociate), while the segmental stretch is facilitated by the closed stickers. Importantly, our analysis revealed that the incorporation of stickers enables a significant reduction in the input of specific work needed to achieve large stretch fluctuations and, consequently, may reduce the energy requirements for FIC.
Focusing on our finding that chain alignment at low, nonstretching, flow rates requires less specific work in the absence of stickers (and presumably for low sticker lifetimes) than with stickers, while the stretching of the chains at high rates is helped by long sticker lifetimes, we speculate that control over both the structural aspects of the final material and over the specific work needed is possible through time- or position-dependent sticker lifetimes. We argue this can be achieved through external chemical control. Indeed, during its larval life cycle, the silkworm stably stores its silk solution at a high viscosity, but just prior to silk spinning it lowers the viscosity through an increase in the potassium concentration through a decreasing lifetime of calcium bridges (stickers) [8,9]. This we can now interpret as a mechanism to ease chain alignment in flow. Intriguingly, downstream the spinning duct the acidity increases [34], which we expect to increase the stability and hence the lifetime of the calcium bridges, and hence enhance local chain stretching, see Fig. 1, which may, in turn, disrupts the solvation layer of the protein and induce efficient crystallization [7,13,15–18].
While this seems a compelling mechanism for efficient FIC, it is not yet clear how this process may be optimized. The experimental accessibility of these and other questions has come in reach owing to recent advances in controlling the content of metal cations in silk feedstock [70]. In the case of B. mori silk, we identified a regular spacing of the negative charges along the backbone of the chain, with strands of approximately 500 uncharged amino acids between; the length of these sticker strands is of the order of the entanglement molecular weight [9]. The regularity of the spacing and the coincidental similarity between the number of stickers and entanglements suggests some degree of evolutionary optimization. The functionality of ordered- versus random copolymers is of high importance from a synthetic polymer chemistry point of view and needs to be addressed using simulations that include both associations and entanglements.
We conclude that our modeling approach leaves us well prepared to investigate the ways in which the evolution of silk-producing organisms may have exploited the potential optimal strategies for efficient fiber processing. The next piece of physics to add to this account of the rheology of polymers with temporary associations, not only for modeling silk proteins but also general associating polymers, concerns the interaction between entanglements and associations in strong flow. We anticipate that this will further enrich the ongoing debate in polymer physics on the physics of entanglement generation and destruction (i.e., “entanglement stripping”) in nonlinear rheo-physics, as well as continue the account of how silk-forming organisms point to novel rheo-physics of flow-induced phase-transformations.
ACKNOWLEDGMENTS
This research was funded by the Engineering and Physical Sciences Research Council (Grant No. EPSRC EP/N031431/1). C.S. is grateful for the computational support from the University of York high-performance computing service (the Viking Cluster) and thanks Daniel J. Read for his suggestion to pair closed stickers of different chains. C.S. and T.C.B.M. thank Pete Laity and Chris Holland for useful discussions on the mechanism of silk self-assembly.
AUTHOR DECLARATIONS
Conflict of Interest
The authors declare that they have no conflicts of interest.
APPENDIX
1. Algorithm
Because of the large distribution of chain stretch in the conditions we are interested in, there is also a large distribution of opening rates; in our previous work, we used small time steps in which the chain conformation was updated, and each closed pair had a sufficiently small opening probability. Here, we significantly improve this algorithm by enabling much larger time steps between conformational updates, and during which the stickers may open and close many times, see Fig. 3.
In our algorithm, we update the chain conformation using the Brownian dynamics equation from Sec. II A using a time span . Depending on the opening and closing rates, during this time span, , the sticker configuration may be updated many times or not at all according to a kMC scheme [64–66]. In every kMC step, the rate at which any opening or closing event may occur is calculated as , with
the sum of closing rates and
the sum of dissociation rates, where differs for the different sticker pairs due to dispersity in chain tension. In these expressions, and are the number of open and closed stickers, respectively; is the total number of possible associations, and the index sums over all pairs of closed stickers. Using this sum of rates, the time at which the first opening or closing event occurs is
with being the uniform random number (our code uses random numbers using the open-source SFMT library [71]). If exceeds the time span , no opening or closing events occurs. However, if , then a second random number is drawn, and a closing event is selected with probability , and a dissociation event is selected with probability . After updating the configurations of the stickers, the time span is updated to . The kMC scheme is terminated when , following which the chain conformation is updated.
While in the linear rheological conditions we solve the dynamics using a fixed time step, in strong flow we implemented an adaptive time step to handle the large and fast fluctuations in stretch that emerge in some parameter regimes of the system. In every iteration , the time step for the next iteration is updated as
where an error and tolerance are calculated for the change of each end-to-end vector . We defined the error value for each change in as , with set by . For the tolerance value, we use scalar values and depending on whether is smaller or larger than a certain cutoff set by . Above the cutoff, we avoid numerical instabilities due to the singularity at by using
For continuity of the derivative, , with ; for a cutoff even this smooth potential is steep () and in practice we use a softer potential ().
2. Asymptotic limits of the two-state model
The two-state master equation in Eqs. (21) and (22) has analytical solutions for early times where advection dominates over the sticker dynamics, and for late times where the sticker dynamics is fast compared to the rate by which the deep tail of the stretch distribution fills up. We obtain these analytical solutions in both cases using the Laplace transform of Eqs. (21) and (22) in the limit of large stretches , which is
where is the Laplace transform of for [hence, we have used the standard Laplace transform of the time derivative ]. We will obtain the early- and late-stage solutions by using different initial conditions at and boundary conditions that we will discuss below.
Focusing first on the early-stage limit, we consider a narrow distribution of chain segments between closed stickers at time , with being the Dirac delta distribution. For early times, these segments stretch exponentially with time until the stickers open. To inspect how these segments evolve, we insert the initial conditions into Eq. (A6), which gives
with being the Laplace transform of . The solution is of the standard form , which after inverse Laplace transform gives Eq. (39) in the main text.
To solve Eqs. (21) and (22) in the long-time limit, we make the useful approximation that at an intermediate time the distribution is at steady state for small stretches , while the large-stretch tail of the distribution is unoccupied. Hence, at , the distribution is given by
where and is the Heaviside step function. The prefactor,
normalizes the distribution. We now set to a large value, so , and at late times the filling of the tail of the distribution (for ) occurs with a negligible effect on the distribution at small stretches.
Of which the solution is of the form
with and being the eigenvalues given by
and where the coefficients, , follow from the boundary condition at .
At late times, i.e., for small , we have , where is given by Eq. (25), and where
are both positive, provided that the sticky Weissenberg number is sufficiently small, [10], where is the Weissenberg number of the chain without stickers.
From the boundary condition, we find that the coefficients must be of the form . As the “+” solution leads to a non-normalizable solution, however, , and the solution is
Finally, after taking the inverse Laplace transform, we have
Hence, the exponentially extending front of the distribution is located at the stretch ratio
We have checked the validity of our interpretation of a narrow moving-front by also calculating the width of this front. To do this, we consider the relaxation function with again , and and are the transient and steady-state stretch distributions, respectively. A narrow front that reaches at time and reaches a steady state at time may be approximated by
The Laplace transform of this function is
We compare this result to the solution of the two-state model in Eq. (A14) through a second-order Taylor expansion of the exponential terms
From the linear term, we find [as expected from Eq. (29)]. After substitution into the second term and elimination of this variable, we find the width of the front to be
The relative width, compared to the location of the front (), is
The relative width calculated in the time domain also represents the relative width of the (logarithmic) stretch distribution,
Upon approaching the strain rate where the mean stretch diverges, i.e., at , the relative width of the front diverges as . In this equation, the bare Weissenberg number is . Hence, if the sticker lifetime is and the fraction of closed stickers is (as in our simulations), then significant broadening of the front only happens very close to the stretch transition, . This verifies that our approximation of a narrow front is indeed accurate.
3. Power-law regression
To determine the sticky-Rouse diffusivity, , from the MSD of the center of mass,
as a function of time , and the stretch exponent, , from the probability distribution,
as a function of the stretch ratio, , we write both equations in the form
and perform common linear regression. However, because both power-laws represent asymptotic behavior for large , there is also a cutoff value, , above which they apply. We determine the cutoff by minimizing
with respect to , , and (note that ); is the uncertainty on the simulated data. Here, we set fixed and the number of free parameters for extracting the sticky-Rouse diffusivity from the MSD data. To determine the stretch exponent () from the stretch distributions, we use the same approach but leave as a free fitting parameter and set .