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.

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 pH along the spinning duct, suggesting an exquisitely controlled local rheology [34]. While lower pH 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.

FIG. 1.

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

FIG. 1.

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

Close modal

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,

(1)

required to induce FIC after a period time ts during which the system is subjected to the (experimentally controllable) transpose of the (local) velocity-gradient tensor κ=vT 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 ts render the specific work a control variable (Wσxyγ˙ts) 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 1% [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, ΔS. Hence, the nucleation rate increases as Jexp(ΔS/kB), where one may assume TΔSEW, with E 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, p and their lifetime, τs [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

(2)

where i is a chain segment at position Ri 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, ζi, the (friction-dependent) thermal forces, Fthermal,i, and the network forces, Fnetwork,i. 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, Fnetwork,i=Fintra,iFthermal,i. 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.

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 N monomers may be discretized using a number of nodes, Nnodes, 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 i is located at a spatial coordinate Ri relative to the center of mass of the chain. The strand between neighboring nodes i and i+1 has an end-to-end vector Q=Ri+1Ri and contains a fraction Δsi=Ns,i/(N+1) of all the monomers in the chain. At this level of coarse graining, the friction of each node is given by

(3)

with ζ0 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 Δsi0 at i=0 and at i=Nnodes [59].

FIG. 2.

(color online) The theory in Sec. II A applies to sticky entangled polymers that are parameterized using the locations of M 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 ζi that depends on the fraction of monomers of the chain, Δsi, that reside in each of the M+1 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.

FIG. 2.

(color online) The theory in Sec. II A applies to sticky entangled polymers that are parameterized using the locations of M 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 ζi that depends on the fraction of monomers of the chain, Δsi, that reside in each of the M+1 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.

Close modal

The equilibrium structure of the chain in quiescent conditions is determined by the end-to-end distance of the substrands, |Qi|=λb(ΔsiN)1/2, where the stretch ratio λ obeys the equilibrium distribution

(4)

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 N monomers with a mean stretch ratio of unity

(5)

where

(6)

approximately captures the anharmonicity of the spring force due to the finite extensibility of the substrand [60]. For the substrands i, 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 FintraFintra,i, NΔsiN, and λmaxΔsi1/2λmaxλmax,i. The direction of the force exerted by spring i on node i is Qi/|Qi|, while the direction of this force acted upon node i+1 is Qi/|Qi|. Hence, the net intramolecular force exerted on node i is

(7)

The thermal force is given by the equipartition theorem

(8)
(9)
(10)

with α,β=x,y,z being the Cartesian coordinates and kBT 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

(11)

where κ is the transpose of the velocity-gradient tensor, which in extension and shear is given by

(12)

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,

(13)

using which we recursively obtain the drift of the nodes as

(14)

The value of the first entry, R1/t, 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 i from chain A and sticker j 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 ζiA+ζjB, where ζiA and ζjB 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

(15)

Hence, the paired stickers i and j have an identical friction coefficient and experience the same net force Fintra,iA+Fintra,jB+Fthermal,iA (where Fthermal,iA=Fthermal,jB). Crucially to forced sticker dissociation, the net force that acts on the closed sticker pair is

(16)

which we assume, as in other cases of forces temporary unbinding, lowers the activation energy for sticker dissociation as

(17)

with Eact0 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] τs=ν1exp(Eact/kBT), 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 p/(1p)2=K0exp(0Fstic), with K0 being the equilibrium constant in the absence of any chain tension. Here, the free energy 0Fstic>0 captures the shift in detailed balance (i.e., the fraction of closed stickers decreases with an increase in chain tension), while Fstic 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 kopen=νexp([θ0FsticEact0]/kBT) and kclose=νK0exp([(1θ)0Fstic+Eact0]/kBT), respectively, where θ0, and where θ[0,1] 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 =0 and θ=1; 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 Δt a sticker is opened or closed with a probability (1exp[kopenΔt]) or (1exp[kcloseΔt]), 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.

FIG. 3.

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

FIG. 3.

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

Close modal

The dynamics of sticky polymers is complicated by the fact that a polymer with Zs stickers can be in 2Zs different states, as each individual sticker can be either opened or closed. An instructive simple case is a chain with Zs=2, 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],

(18)

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 kopen, and the closing rate is kclose. The time development of the probability distribution of the stretch ratio is described by [10]

(19)
(20)

with τR 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 ylnλ, so Pi/λ=(1/λ)Pi/lnλexp(y)Pi/y. Similarly, λPi/λ=Pi+Pi/y. Inserting this into the governing equations gives

(21)
(22)

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 ey0, which is equivalent to λ1.

In steady state, the left-hand side of the equation is zero and the equations can be cast in the form dP/dy=AP, with P=[P0,P1]T and A a constant 2 by 2 matrix. The solution of this system of first-ordinary differential equations is given by [10]

(23)
(24)

with c being the normalization constant (its value can in principle be determined by releasing the approximation ey0), and with the exponent of the power-law distribution given in terms of physical parameters by

(25)

[this is one of the eigenvalues of Eqs. (21) and (22); the other eigenvalue is 1 and is unphysical as a distribution of the form λ1 cannot be normalized.] The value of this stretching exponent diverges if the bare stretch transition at ε˙τR=1 is approached from small strain rates. However, because of the physics of the stickers, actual divergence already occurs at lower strain rates: At ε˙τR=(1p), the exponent becomes ν=1 and the stretch distribution can no longer be normalized. Depending on the sticker lifetime, at smaller strain rates the exponent may reach a value ν=2 if the “sticky Weissenberg number” (1p)ε˙τR reaches unity; here, the mean stretch diverges. While the mean stretch is finite for smaller strain rates, the variance of the stretch diverges for ν3, which happens if (1p)ε˙τR becomes larger than 1/2 [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,” λ(t) (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), Pi(t=0,λ)=ciδ(λλ(0) at λ(0)), 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, P0(t,λ)=c0δ(λλ(0)exp[ε˙t]) and retract to smaller stretches for the open state P1(t,λ)=c1δ(λλ(0)exp[(τR1ε˙)t]). This suggests that the “front,” λ(t), of any distribution with finite P0, shifts exponentially in time to higher values through λ(t)=λ(0)exp[ε˙t].

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 t0τSR where sufficient stickers have opened to facilitate chain relaxation and assume that the stretch distribution has reached a steady-state for small stretches λ<λ(t0) but is empty for larger stretch ratios. Here, λ(t0) 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 t>t0, further convergence of the stretch distribution takes place in the range of stretches λ(t0)<λ<λ(t), where the “front” of the distribution shifts to high stretch values as ln[λ(t)/λ(t0)]ε˙(tt0). Assuming that λ(t0)1, 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 t>t0 the λ<λ(t) portion of the stretch distribution becomes independent of time beyond t>t. The constancy of the distribution at λ(t0) 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 exp(sτ(λ))/s, 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

(26)
(27)

with ν being the “steady-state stretch exponent” in Eq. (25) and with

(28)

being the “dynamic stretch exponent,” which controls the growth of the front of the distribution as

(29)

In this equation, Wi=ε˙τR and Wisticky=ε˙τSR are the (extensional) Weissenberg numbers of the chain without and with stickers, respectively; within the two-state model, τSR=(1p)/kopen, see discussion under Eq. (25). Upon approaching the stretch transition Wisticky=1 where the mean stretch diverges, ν0 indicates “critical slowing down,” as the (late-stage) front of the distribution becomes immobile. For chains with strong stickers (1p)τsτR at the strain rate Wisticky=1/2 where the variance of the stretch diverges [see discussion under Eq. (25)], we find ν2, which indicates that the late-stage measure of the front is shifted from the early-stage measure for the outliers by a factor 2. We have also checked that the moving front is narrow for small strain rates Wisticky<1/2. 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 ΔrelpWiWisticky/(1Wisticky), where Δrel([P(λ,t)/Peq(λ,)]/lnλ)1/lnλ. As we shown in Subsection 2 of the  Appendix, typically this width is Δrel1, and the front of the distribution is narrow even close to the stretch transition.

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 M=4 to 64 (the beads are regularly along the backbone of the polymer, so Δsi=1/(M+1) for all i). Figure 4 shows that the choice of the number of beads has a negligible influence on the time evolution of the mean-square displacement, MSD, of the center of mass and is in all cases in agreement with the theoretical prediction

(30)

where the diffusivity, D, is for nonsticky polymers given by the bare Rouse diffusivity

(31)

Moreover, the inset of Fig. 4 shows that also the end-to-end-distance, Re, is distributed according to the physical equilibrium result of Eq. (4).

FIG. 4.

MSD of the center of mass of a nonsticky polymer against time (the main panel) and the time-averaged end-to-end length (Re) 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, M, per chain. The symbols and solid black curves represent the simulations and the theory, respectively.

FIG. 4.

MSD of the center of mass of a nonsticky polymer against time (the main panel) and the time-averaged end-to-end length (Re) 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, M, per chain. The symbols and solid black curves represent the simulations and the theory, respectively.

Close modal

For times shorter than the Rouse time of strands between stickers, i.e., for t<τR(Zs+1)2, 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 τSR, which is approximately given by [19]

(32)

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 10, i.e., when the center of mass of the chain has diffused 34 times its end-to-end distance.

FIG. 5.

Linear rheology of a sticky chain with Zs=10, p=0.9, τs=200τR within the rigid-network approximation (open symbols) and with this approximation released (closed symbols). (a) MSD of the center of mass against time. (b) Storage, G, and loss, G, modulus in units of G0 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 Zs=4 (dashed curves).

FIG. 5.

Linear rheology of a sticky chain with Zs=10, p=0.9, τs=200τR within the rigid-network approximation (open symbols) and with this approximation released (closed symbols). (a) MSD of the center of mass against time. (b) Storage, G, and loss, G, modulus in units of G0 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 Zs=4 (dashed curves).

Close modal

The elastic compliance not only affects the subdiffusive motion of the chain, but also the sticky-Rouse diffusivity DSR=DRτR/τSR 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 G and G against the frequency ω in Fig. 5(b). The data shown include nonassociating unentangled chains (Zs=0) and the unentangled sticky chains of Fig. 5(a); i.e., chains with Zs=10 stickers within the rigid-network approximation and with an elastically compliant network. The simulated data (symbols) were obtained from the relaxation modulus, G(t), through the multiple-tau-correlator algorithm discussed in [67]. To obtain the dynamic moduli G and G, we have used the finite-element approach from [68]. We have compared the data to the sticky-Rouse model (curves), which is given by

(33)

In this equation, the first summation captures the high-frequency bare Rouse modes (the number of Kuhn segments, N, truncates the highest frequencies) and the second summation captures the sticky-Rouse modes. The modulus G0 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 G,Gω1/2. Discrepancies, such as a roll-off of G 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 ω=1/τs. Within the rigid-network approximation (closed circles), the modulus of the plateau is G(ω)=G0Zs in agreement with the sticky-Rouse model in Eq. (33) for Zs=10. 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 Zs=4 (dashed curves). At lower frequencies ω<1/τs, 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 G is much narrower than in the theory. We have estimated that the shape of this peak is best described by Zs=4 within the rigid-network approximation and Zs=3 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 5,10,20 stickers with various sticker lifetimes, in particular, in the regime where the sticky-Rouse diffusivity scales with the sticker lifetime as DSR=DRτR/τSR1/τsZs2, see Eq. (32). Panel (b) shows that upon releasing the rigid-network approximation this scaling behavior persists, but rescaled with a prefactor 4. 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.

FIG. 6.

Sticky-Rouse diffusivity, DSR, against the sticker lifetime, τs for chains with Zs=2,5,10,20 stickers with p=0.9 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 DR and the bare Rouse time, τR.

FIG. 6.

Sticky-Rouse diffusivity, DSR, against the sticker lifetime, τs for chains with Zs=2,5,10,20 stickers with p=0.9 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 DR and the bare Rouse time, τR.

Close modal

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.

FIG. 7.

Representation of simulated chain conformations in extensional flow for ε˙τR=2 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 50Re, which is 65% of the fully extended chain.

FIG. 7.

Representation of simulated chain conformations in extensional flow for ε˙τR=2 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 50Re, which is 65% of the fully extended chain.

Close modal

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 (Zs=0) and sticky polymers below (Zs=2) and above (Zs=5) 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, Wisticky=ε˙τSR, with τSR being the sticky-Rouse time we expect divergent stretching (albeit that real divergence is obstructed by the maximum chain extensibility λmax=75). We have calculated the sticky-Rouse time as τSR=[DR/DSR]τR, with the ratio between the sticky and the bare diffusivity as presented above in Fig. 6. The relevant results are summarized in Table I.

FIG. 8.

Simulated steady-state stretch distributions of the end-to-end distance, Re, for various extension [(a), (c), and (e)] and shear [(b), (d), and (f)] rates for a linear unentangled, nonsticky (Zs=0), and sticky (Zs=2 and Zs=5) polymers. For these simulations, τSRτs=10τR (see Table I for all parameter values). The black curve represents the contour-length fluctuations in quiescent conditions, given by Eq. (4).

FIG. 8.

Simulated steady-state stretch distributions of the end-to-end distance, Re, for various extension [(a), (c), and (e)] and shear [(b), (d), and (f)] rates for a linear unentangled, nonsticky (Zs=0), and sticky (Zs=2 and Zs=5) polymers. For these simulations, τSRτs=10τR (see Table I for all parameter values). The black curve represents the contour-length fluctuations in quiescent conditions, given by Eq. (4).

Close modal
TABLE I.

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 modelDSR/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 modelDSR/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 (Zs=0), a broad stretch distribution with a cutoff set by λmax 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 0.9<ε˙τR<1.1 around the bare stretch transition, Wi=ε˙τR=1. Beyond the stretch transition, the stretch distribution is narrow and Gaussian and approaches λmax 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 ε˙τSR=1, with the sticky-Rouse time obtained from the sticky-Rouse diffusivity of Fig. 6 as τSR=τRDSR/DR.) 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, τs in Eq. (25), which models the simultaneous opening of all stickers, is replaced by τs/2, and results in

(34)

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.

FIG. 9.

(color online) Stretch exponent ν of the power-law tail of the stretch distribution Pλν for simulations of polymers with Zs=2 (circles) and 5 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 τs=10τR (see Table I for all physical parameter values). For ν>3 (horizontal line), the fluctuations in stretch diverge; this leads to a cutoff in the stretch distribution for chains with finite extensibility, see Fig. 8.

FIG. 9.

(color online) Stretch exponent ν of the power-law tail of the stretch distribution Pλν for simulations of polymers with Zs=2 (circles) and 5 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 τs=10τR (see Table I for all physical parameter values). For ν>3 (horizontal line), the fluctuations in stretch diverge; this leads to a cutoff in the stretch distribution for chains with finite extensibility, see Fig. 8.

Close modal

For the chains with two and five stickers and with a fraction p=0.9 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 ν3, this cutoff has a significant effect on the tail of the stretch distribution upon approaching ν=3. 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.

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, σxy/γ˙, and extensional stress (σyyσrr)/ε˙, with the stress tensor (in units of energy per molecule) given by

(35)

Focusing first on the results for nonsticky chains with a finite extensibility λmax=75 in Figs. 10(a) and 10(b), we reproduce the well-known qualitative features of their stress transient [57]. For small Weissenberg numbers, ε˙τR<1, γ˙τR<1 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 λmax.

FIG. 10.

(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 ε˙τs=0.5, which is below the stretch transition (at ε˙τs1, see Table I). (b) and (d) Transient stretch distribution of the end-to-end distance, Re, 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.

FIG. 10.

(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 ε˙τs=0.5, which is below the stretch transition (at ε˙τs1, see Table I). (b) and (d) Transient stretch distribution of the end-to-end distance, Re, 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.

Close modal

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 Zs=5. 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 ε˙τR>1, the sticky ones also show strain hardening for smaller strain rates ε˙τs>1. 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 0.3<ε˙τs<0.5, 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 (ε˙τs=2) 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 lnλε˙t until the steady state is reached after a time ε˙τlnλmax. This is in qualitative agreement with the predictions of the two-state model that we derived in Eq. (29) of Sec. II B.

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,” P2,K[0,1]. If P2,K1 (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, Ps, of chain segments of some length Δs[0,1] have stretched beyond a critical stretch ratio Lsλmax , where λmax=λmaxΔs is the maximum stretch of the chain segment and Ls[0,1] 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

(36)

where P(.) is the transient stretch distribution function and ts 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 ts 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 Δs=1/(Zs+1); in general, however, measures for Ps, Ls, and Δs will have to be determined through experimentation and (atomistic) MD simulations [15–18].

In this section, we will employ simulations with 50 chains of a fixed number of 11 beads (i.e., with 10 chain segments, giving Δs=1/10), and we will monitor the maximum stretch among the total of 500 chain segments (i.e., Ps=1/500). The time-evolution of the maximum stretch will enable us to screen how various values of Ls require a different processing time ts 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 P2,s[0,1], which is the largest eigenvalue of the nematic order tensor P2,s=(3uu1)/2, where u 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., P2,s>P2,K.) In Fig. 11, we have calculated the critical specific work, W, as given in Eq. (1), needed to achieve values of P2,s and Ls in the range from 0 to 1 for nonsticky (Zs=0) and sticky (Zs=5) chains for various shear and extensional rates.

FIG. 11.

(color online) Nematic order parameter, P2,s and characteristic stretch ratio, Ls, 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 Zs=5 with an elastically compliant network. All physical parameter values are given in Table I.

FIG. 11.

(color online) Nematic order parameter, P2,s and characteristic stretch ratio, Ls, 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 Zs=5 with an elastically compliant network. All physical parameter values are given in Table I.

Close modal

The top panels of this figure give the nematic order parameter, P2,s, and the measure for stretch fluctuations, Ls 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 ε˙τR=1). 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.

FIG. 12.

The critical time (left) and the specific critical work (right) against the sticky Weissenberg number, Wisticky=ε˙τSR,γ˙τSR, for various Ls and P2,s 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 Zs=5 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 Wisticky=1, where the “sticky” Weissenberg number is Wisticky10Wi=10ε˙τR, with Wi being the Weissenberg number of the nonsticky chain. This factor 10 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 Ls=0.6 and for Ls=0.8.

FIG. 12.

The critical time (left) and the specific critical work (right) against the sticky Weissenberg number, Wisticky=ε˙τSR,γ˙τSR, for various Ls and P2,s 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 Zs=5 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 Wisticky=1, where the “sticky” Weissenberg number is Wisticky10Wi=10ε˙τR, with Wi being the Weissenberg number of the nonsticky chain. This factor 10 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 Ls=0.6 and for Ls=0.8.

Close modal

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 ε˙τR0.1 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 Ps=1/500 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 Ls or P2,s criterion for critical nucleation, we are interested how the strain rate affects how much critical specific work, W, is needed, and at what time scale, ts this criterion is achieved. To investigate this, we focus on horizontal lines/cross sections of Fig. 12 (i.e., at fixed values 0.6 and 0.8 of both Ls and P2,s). For the data points along these lines, we plot the critical work, W, and the time scale, ts, in Fig. 12. The left panel shows that the time scale scales as tsWi1, 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, P2,s, or of stretch fluctuations, Ls, 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 Ls. We remark that for the system we studied, the stretch transition in the absence of stickers takes is located at Wisticky=ε˙τSR10 (because τSR10τR, see Table I). For smaller strain rates, Wisticky<10, we find there is a minimum in the specific critical work near the stretch transition Wisticky1. Indeed, while large stretches are achieved just below the stretch transition Wisticky<1 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,

(37)

with c being the constant, assuming that the open chains are in a relaxed state. Here, P0(λ,t) 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 W=0tsσxxε˙dt. To calculate W, we first will determine ts using the criterion

(38)

which, as before, implies a minimum concentration of chains with a stretch ratio of at least λs=Lsλmax,i. Second, we will need an expression for the time evolution of the probability density P0.

To obtain P0, 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 σxx. Therefore, we will only take into account the loss of strongly stretched chains by opening rate kopen and ignore the contribution of closing events by rate kclose. We will further use the initial condition P(λ,0)=δ(1λ), with δ(.) being the Dirac delta distribution to represent a narrow stretch distribution at time t=0. The dynamical equation in Eq. (21) then predicts that the Dirac delta distribution in time shifts to high stretch values along the λ axis, as

(39)

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 ts=lnλs/ε˙, which is in agreement with our simulated results displayed in Fig. 11. We insert this equation into the expression for the critical specific work, W=0tsσxxε˙dt, and find

(40)

where ε˙min 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 ε˙τs=1 (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 λs. 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.

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.

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.

The authors declare that they have no conflicts of interest.

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 Δt. Depending on the opening and closing rates, during this time span, Δt1Δt, 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 WT=Wa+Wd, with

(A1)

the sum of closing rates and

(A2)

the sum of dissociation rates, where kd,q differs for the different sticker pairs due to dispersity in chain tension. In these expressions, Nopen and Nclosed are the number of open and closed stickers, respectively; Nopen(Nopen1)/2 is the total number of possible associations, and the index q sums over all Nclosed/2pairs of closed stickers. Using this sum of rates, the time Δt2 at which the first opening or closing event occurs is

(A3)

with u(0,1] being the uniform random number (our code uses random numbers using the open-source SFMT library [71]). If Δt2 exceeds the time span Δt1, no opening or closing events occurs. However, if Δt2<Δt1, then a second random number [0,1] is drawn, and a closing event is selected with probability ka/WT, and a dissociation event q is selected with probability kd,q/WT. After updating the configurations of the stickers, the time span is updated to Δt1=Δt1+Δt2. The kMC scheme is terminated when Δt2>Δt1, 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 n, the time step for the next iteration is updated as

(A4)

where an error and tolerance are calculated for the change of each end-to-end vector Qi. We defined the error value for each change in Qi as error=|ΔQin|/Qmax, with Qmax set by λmax. For the tolerance value, we use scalar values tol and tol+ depending on whether |Qin| is smaller or larger than a certain cutoff set by λcutoff<λmax. Above the cutoff, we avoid numerical instabilities due to the singularity at λmax by using

(A5)

For continuity of the derivative, α=4c2/(34c2+c4), with c=λcutoff/λmax; for a cutoff λcutoff=0.9λmax even this smooth potential is steep (α8) and in practice we use a softer potential (α=4).

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 λ>λ1, which is

(A6)
(A7)

where Pi~(s,y)L{Pi(t,y)} is the Laplace transform of Pi for i=0,1 [hence, we have used the standard Laplace transform of the time derivative L{Pi/t}=sP~i(s,y)Pi(0,y)]. We will obtain the early- and late-stage solutions by using different initial conditions Pi(0,y) at t=0 and boundary conditions that we will discuss below.

Focusing first on the early-stage limit, we consider a narrow distribution P(λ,0)=δ(1λ) of chain segments between closed stickers at time t=0, 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

(A8)

with P~0(λ,s) being the Laplace transform of P0(λ,s). The solution is of the standard form P~0exp(sτ), 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 t=t the distribution is at steady state for small stretches λλ, while the large-stretch tail of the distribution is unoccupied. Hence, at t=t, the distribution is given by

(A9)
(A10)

where ylnλ and Θ is the Heaviside step function. The prefactor,

(A11)

normalizes the distribution. We now set λ to a large value, so cc, and at late times t>t 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

(A12)
(A13)

with ν(s) and ν+(s) being the eigenvalues given by

(A14)

and where the coefficients, ci±, follow from the boundary condition at y=y.

At late times, i.e., for small s, we have ν(s)νeq(s/ε˙)ν+(1/2)(s/ε˙)2ν, where νeq is given by Eq. (25), and where

(A15)
(A16)

are both positive, provided that the sticky Weissenberg number is sufficiently small, WistickyWi/(1p)<1 [10], where Wi=ε˙τR is the Weissenberg number of the chain without stickers.

From the boundary condition, we find that the coefficients must be of the form ci±1/s. As the “+” solution leads to a non-normalizable solution, however, ci+=0, and the solution is

(A17)
(A18)

Finally, after taking the inverse Laplace transform, we have

(A19)
(A20)

Hence, the exponentially extending front of the distribution is located at the stretch ratio

(A21)

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 f(t)=P(y,t)/Peq(y) with again y=lnλ, and P and Peq are the transient and steady-state stretch distributions, respectively. A narrow front that reaches y at time τ and reaches a steady state at time τ+Δ may be approximated by

(A22)

The Laplace transform of this function is

(A23)

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

(A24)

From the linear term, we find τ+Δ/2=(ν/ε˙)lny [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

(A25)

The relative width, compared to the location of the front (τ+Δ/2), is

(A26)

The relative width calculated in the time domain also represents the relative width of the (logarithmic) stretch distribution,

(A27)

Upon approaching the strain rate where the mean stretch diverges, i.e., at Wisticky=1, the relative width of the front diverges as Δrel24pWiWisticky/(1Wisticky). In this equation, the bare Weissenberg number is Wi=Wisticky(1p)τR/τs. Hence, if the sticker lifetime is τs=10τR and the fraction of closed stickers is p=0.9 (as in our simulations), then significant broadening of the front only happens very close to the stretch transition, Wisticky>0.99. This verifies that our approximation of a narrow front is indeed accurate.

3. Power-law regression

To determine the sticky-Rouse diffusivity, DSR, from the MSD of the center of mass,

(A28)

as a function of time t, and the stretch exponent, ν, from the probability distribution,

(A29)

as a function of the stretch ratio, λ, we write both equations in the form

(A30)

and perform common linear regression. However, because both power-laws represent asymptotic behavior for large x, there is also a cutoff value, xcutoff, above which they apply. We determine the cutoff by minimizing

(A31)

with respect to a, b, and i0 (note that xi0=xcutoff); σi is the uncertainty on the simulated y data. Here, we set b=1 fixed and the number of free parameters Npar=1 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 b as a free fitting parameter and set Npar=2.

1.
Graham
,
R. S.
, and
P. D.
Olmsted
, “
Coarse-grained simulations of flow-induced nucleation in semicrystalline polymers
,”
Phys. Rev. Lett.
103
,
115702
(
2009
).
2.
Troise
,
E. M.
,
H. J. M.
Caelers
, and
G. W. M.
Peters
, “
Full characterization of multiphase, multimorphological kinetics in flow-induced crystallization of IPP at elevated pressure
,”
Macromolecules
50
,
3868
3882
(
2017
).
3.
Nicholson
,
D. A.
, and
G. C.
Rutledge
, “
An assessment of models for flow-enhanced nucleation in an n-alkane melt by molecular simulation
,”
J. Rheol.
63
,
465
475
(
2019
).
4.
Moghadamand
,
S.
,
I. S.
Dalal
, and
R. G.
Larson
, “
Unraveling dynamics of entangled polymers in strong extensional flows
,”
Macromolecules
52
,
1296
1307
(
2019
).
5.
Graham
,
R. S.
, “
Understanding flow-induced crystallization in polymers: A perspective on the role of molecular simulations
,”
J. Rheol.
63
,
203
214
(
2019
).
6.
Read
,
D. J.
,
C.
McIlroy
,
C.
Das
,
O. G.
Harlen
, and
R. S.
Graham
, “
Polystrand model of flow-induced nucleation in polymers
,”
Phys. Rev. Lett.
124
,
147802
(
2020
).
7.
Holland
,
C.
,
F.
Vollrath
,
A. J.
Ryan
, and
O. O.
Mykhaylyk
, “
Silk and synthetic polymers: Reconciling 100 degrees of separation
,”
Adv. Mater.
24
,
105
109
(
2012
).
8.
Laity
,
P. R.
,
E.
Baldwin
, and
C.
Holland
, “
Changes in silk feedstock rheology during cocoon construction: The role of calcium and potassium ions
,”
Macromol. Biosci.
19
,
1800188
(
2018
).
9.
Schaefer
,
C.
,
P. R.
Laity
,
C.
Holland
, and
T. C. B.
McLeish
, “
Silk protein solution: A natural example of sticky reptation
,”
Macromolecules
53
,
2669
2676
(
2020
).
10.
Schaefer
,
C.
, and
T. C. B.
McLeish
, “
Power-law stretching of associating polymers in steady-state extensional flow
,”
Phys. Rev. Lett.
126
,
057801
(
2021
).
11.
Schaefer
,
C.
,
P. R.
Laity
,
C.
Holland
, and
T. C. B.
McLeish
, “
Stretching of Bombyx mori silk protein in flow
,”
Molecules
26
,
1663
(
2021
).
12.
Asakura
,
T.
, “
Structure of silk I (Bombyx mori silk fibroin before spinning)-type II β-turn, not α-helix-
,”
Molecules
26
,
3706
(
2021
).
13.
Dunderdale
,
G. J.
,
S. J.
Davidson
,
A. J.
Ryan
, and
O. O.
Mykhaylyk
, “
Flow-induced crystallisation of polymers from aqueous solution
,”
Nat. Commun.
11
,
3372
(
2020
).
14.
Laity
,
P. R.
, and
C.
Holland
, “
Seeking solvation: Exploring the role of protein hydration in silk gelation
,”
Molecules
25
,
551
(
2020
).
15.
Donets
,
S.
, and
J.-U.
Sommer
, “
Molecular dynamics simulations of strain-induced phase transition of poly(ethylene oxide) in water
,”
J. Phys. Chem. B
122
,
392
397
(
2018
).
16.
Donets
,
S.
,
O.
Guskova
, and
J.-U.
Sommer
, “
Flow-induced formation of thin PEO fibres in water and their stability after the strain release
,”
J. Phys. Chem. B
124
,
9224
9229
(
2020
).
17.
Donets
,
S.
,
O.
Guskova
, and
J.-U.
Sommer
, “
Searching for aquamelt behavior among silklike biomimetics during fibrillation under flow
,”
J. Phys. Chem. B
125
,
3238
3250
(
2021
).
18.
Mkandawire
,
W. D.
, and
S. T.
Milner
, “
Simulated osmotic equation of state for poly(ethylene oxide) solutions predicts tension-induced phase separation
,”
Macromolecules
54
,
3613
3619
(
2021
).
19.
Leibler
,
L.
,
M.
Rubinstein
, and
R. H.
Colby
, “
Dynamics of reversible networks
,”
Macromolecules
24
,
4701
4707
(
1991
).
20.
Colby
,
R. H.
,
X.
Zheng
,
M. H.
Rafailovich
,
J.
Sokolov
,
D. G.
Peiffer
,
S. A.
Schwarz
,
Y.
Strzhemechny
, and
D.
Nguyen
, “
Dynamics of lightly sulfonated polystyrene ionomers
,”
Phys. Rev. Lett.
81
,
3876
3879
(
1998
).
21.
Hackelbusch
,
S.
,
T.
Rossow
,
P.
van Assenbergh
, and
S.
Seiffert
, “
Chain dynamics in supramolecular polymer networks
,”
Macromolecules
46
,
6273
6286
(
2013
).
22.
Chen
,
Q.
,
Z.
Zhang
, and
R. H.
Colby
, “
Viscoelasticity of entangled random polystyrene ionomers
,”
J. Rheol.
60
,
1031
1040
(
2016
).
23.
Tomkovic
,
T.
, and
S. G.
Hatzikiriakos
, “
Nonlinear rheology of poly(ethylene-co-methacrylic acid) ionomers
,”
J. Rheol.
62
,
1319
1329
(
2018
).
24.
Zhang
,
Z.
,
Q.
Chen
, and
R. H.
Colby
, “
Dynamics of associative polymers
,”
Soft Matter
14
,
2961
2977
(
2018
).
25.
Hinton
,
Z. R.
,
A.
Shabbir
, and
N. J.
Alvarez
, “
Dynamics of supramolecular self-healing recovery in extension
,”
Macromolecules
52
,
2231
2242
(
2019
).
26.
Golkaram
,
M.
, and
K.
Loos
, “
A critical approach to polymer dynamics in supramolecular polymers
,”
Macromolecules
52
,
9427
9444
(
2019
).
27.
Zuliki
,
M.
,
S.
Zhang
,
K.
Nyamajaro
,
T.
Tomkovic
, and
S. G.
Hatzikiriakos
, “
Rheology of sodium and zinc ionomers: Effects of neutralization and valency
,”
Phys. Fluids
32
,
023104
(
2020
).
28.
Liu
,
S.
,
X.
Cao
,
C.
Huang
,
R. A.
Weiss
,
Z.
Zhang
, and
Q.
Chen
, “
Brittle-to-ductile transition of sulfonated polystyrene ionomers
,”
ACS Macro Lett.
10
,
503
509
(
2021
).
29.
Witten
,
T.
, “
Associating polymers and shear thickening
,”
J. Phys.
49
,
1055
1063
(
1988
).
30.
Jin
,
H.-J.
, and
D. L.
Kaplan
, “
Mechanism of silk processing in insects and spiders
,”
Nature
424
,
1057
1061
(
2003
).
31.
Tanaka
,
F.
,
R.
Takeda
, and
K.
Tsurusaki
, “
Critical shear rate for gelation in aqueous solutions of associating polymers under shear flows
,”
J. Phys. Soc. Jpn.
87
,
074801
(
2018
).
32.
Séréro
,
Y.
,
V.
Jacobsen
, and
J.-F.
Berret
, “
Evidence of nonlinear chain stretching in the rheology of transient networks
,”
Macromolecules
33
,
1841
1847
(
2000
).
33.
Tripathi
,
A.
,
K. C.
Tam
, and
G. H.
McKinley
, “
Rheology and dynamics of associative polymers in shear and extension: Theory and experiments
,”
Macromolecules
39
,
1981
1999
(
2006
).
34.
Koeppel
,
A.
,
N.
Stehling
,
C.
Rodenburg
, and
C.
Holland
, “
Spinning beta silks requires both pH activation and extensional stress
,”
Adv. Funct. Mater.
31
,
2103295
(
2021
).
35.
Janeschitz-Kriegl
,
H.
, “
How to understand nucleation in crystallizing polymer melts under real processing conditions
,”
Colloid. Polym. Sci.
281
,
1157
1171
(
2003
).
36.
Mykhaylyk
,
O. O.
,
P.
Chambon
,
R. S.
Graham
,
J.
Patrick
,
A.
Fairclough
,
P. D.
Olmsted
, and
A. J.
Ryan
, “
The specific work of flow as a criterion for orientation in polymer crystallization
,”
Macromolecules
41
,
1901
1904
(
2008
).
37.
Seo
,
J.
,
A. M.
Gohn
,
R. P.
Schaake
,
D.
Parisi
,
A. M.
Rhoades
, and
R. H.
Colby
, “
Shear flow-induced crystallization of poly(ether ether ketone)
,”
Macromolecules
53
,
3472
3481
(
2020
).
38.
Likhtman
,
A. E.
, and
T. C. B.
McLeish
, “
Quantitative theory for linear dynamics of linear entangled polymers
,”
Macromolecules
35
,
6332
6343
(
2002
).
39.
Graham
,
R. S.
,
A. E.
Likhtman
,
T. C. B.
McLeish
, and
S. T.
Milner
, “
Microscopic theory of linear, entangled polymer chains under rapid deformation including chain stretch and convective constraint release
,”
J. Rheol.
47
,
1171
1200
(
2003
).
40.
Sefiddashti
,
M. H. N.
,
B. J.
Edwards
, and
B.
Khomami
, “
Individual chain dynamics of a polyethylene melt undergoing steady shear flow
,”
J. Rheol.
59
,
119
153
(
2015
).
41.
Mohagheghi
,
M.
, and
B.
Khomami
, “
Elucidating the flow-microstructure coupling in the entangled polymer melts. Part I: Single chain dynamics in shear flow
,”
J. Rheol.
60
,
849
859
(
2016
).
42.
Mohagheghi
,
M.
, and
B.
Khomami
, “
Elucidating the flow-microstructure coupling in the entangled polymer melts. Part II: Molecular mechanism of shear banding
,”
J. Rheol.
60
,
861
872
(
2016
).
43.
Huang
,
Q.
,
J.
Ahn
,
D.
Parisi
,
T.
Chang
,
O.
Hassager
,
S.
Panyukov
,
M.
Rubinstein
, and
D.
Vlassopoulos
, “
Unexpected stretching of entangled ring macromolecules
,”
Phys. Rev. Lett.
122
,
208001
(
2019
).
44.
O’Connor
,
T. C.
,
T.
Ge
,
M.
Rubinstein
, and
G. S.
Grest
, “
Topological linking drives anomalous thickening of ring polymers in weak extensional flows
,”
Phys. Rev. Lett.
124
,
027801
(
2020
).
45.
Bonato
,
A.
,
D.
Marenduzzo
, and
D.
Michieletto
, “
Simplifying topological entanglements by entropic competition of slip-links
,”
Phys. Rev. Res.
3
,
043070
(
2021
).
46.
Khalatur
,
P. G.
, and
D. A.
Mologin
, “
Rheologlcal properties of self-associating polymer systems: Nonequilibrium molecular dynamics simulation
,”
J. Mol. Liq.
91
,
205
217
(
2001
).
47.
Mohottalalage
,
S.
,
M.
Senanayake
,
T.
O’Connor
,
G.
Grest
, and
D.
Perahia
, “
Nonlinear elongation flows effects on aggregation in associating polymer melts
,” in APS March Meeting 2021 (APS, 2021), Vol. 66.
48.
Shivokhin
,
M. E.
,
T.
Narita
,
L.
Talini
,
A.
Habicht
,
S.
Seiffert
,
T.
Indei
, and
J. D.
Schieber
, “
Interplay of entanglement and association effects on the dynamics of semidilute solutions of multisticker polymer chains
,”
J. Rheol.
61
,
1231
1241
(
2017
).
49.
Boudara
,
V. A. H.
, and
D. J.
Read
, “
Stochastic and preaveraged nonlinear rheology models for entangled telechelic star polymers
,”
J. Rheol.
61
,
339
362
(
2017
).
50.
Cui
,
G.
,
V. A. H.
Boudara
,
Q.
Huang
,
G. P.
Baeza
,
A. J.
Wilson
,
O.
Hassager
,
D. J.
Read
, and
J.
Mattsson
, “
Linear shear and nonlinear extensional rheology of unentangled supramolecular side-chain polymers
,”
J. Rheol.
62
,
1155
1174
(
2018
).
51.
Hua
,
C. C.
, and
J. D.
Schieber
, “
Segment connectivity, chain-length breathing, segmental stretch, and constraint release in reptation models. I. Theory and single-step strain predictions
,”
J. Chem. Phys.
109
,
10018
10027
(
1998
).
52.
Shanbhag
,
S.
,
R.
Larson
,
J.-I.
Takimoto
, and
M.
Doi
, “
Deviations from dynamic dilution in the terminal relaxation of star polymers
,”
Phys. Rev. Lett.
87
,
195502
(
2001
).
53.
Doi
,
M.
, and
J.-I.
Takimoto
, “
Molecular modelling of entanglement
,”
Philos. Trans. R. Soc.
361
,
641
652
(
2003
).
54.
Schieber
,
J. D.
,
J.
Neergaard
, and
S.
Gupta
, “
A full-chain, temporary network model with sliplinks, chain-length fluctuations, chain connectivity and chain stretching
,”
J. Rheol.
47
,
213
233
(
2003
).
55.
Likhtman
,
A. E.
, “
Single-chain slip-link model of entangled polymers: Simultaneous description of neutron spin-echo, rheology, and diffusion
,”
Macromolecules
38
,
6128
6139
(
2005
).
56.
Andreev
,
M.
, and
G. C.
Rutledge
, “
A slip-link model for rheology of entangled polymer melts with crystallization
,”
J. Rheol.
64
,
213
222
(
2020
).
57.
Dealy
,
J. M.
,
D. J.
Read
, and
R. G.
Larson
,
Structure and Rheology of Molten Polymers
(
Hanser
,
Munich
,
2018
).
58.
Doi
,
M.
, and
S. F.
Edwards
,
The Theory of Polymer Dynamics
(
Clarendon
,
Oxford
,
1986
).
59.
Andreev
,
M.
,
R. N.
Khaliullin
,
R. J. A.
Steenbakkers
, and
J. D.
Schieber
, “
Approximations of the discrete slip-link model and their effect on nonlinear rheology predictions
,”
J. Rheol.
57
,
535
557
(
2013
).
60.
Cohen
,
A.
, “
A padé approximant to the inverse Langevin function
,”
Rheol. Acta
30
,
270
273
(
1991
).
61.
Rubinstein
,
M.
, and
A. N.
Semenov
, “
Thermoreversible gelation in solutions of associating polymers. 2. Linear dynamics
,”
Macromolecules
31
,
1386
1397
(
1998
).
62.
Raffaelli
,
C.
,
A.
Bose
,
C. H. M. P.
Vrusch
,
S.
Ciarella
,
T.
Davris
,
N. B.
Tito
,
A. V.
Lyulin
,
W. G.
Ellenbroek
, and
C.
Storm
,
Rheology, Rupture, Reinforcement and Reversibility: Computational Approaches for Dynamic Network Materials
(
Springer
,
Cham
,
2020
).
63.
Vinu
,
R.
, and
L. J.
Broadbelt
, “
Unraveling reaction pathways and specifying reaction kinetics for complex systems
,”
Annu. Rev. Chem. Biomol. Eng.
3
,
29
54
(
2012
).
64.
Lukkien
,
J. J.
,
J. P. L.
Segers
,
P. A. J.
Hilbers
,
R. J.
Gelten
, and
A. P. J.
Jansen
, “
Efficient Monte Carlo methods for the simulation of catalytic surface reactions
,”
Phys. Rev. E
58
,
2598
2610
(
1998
).
65.
Binder
,
K.
, and
D.
Heermann
,
Monte Carlo Simulation in Statistical Physics
, 5th ed. (
Springer-Verlag
,
Berlin
,
2019
).
66.
Jansen
,
A. P. J.
,
An Introduction to Kinetic Monte Carlo Simulations of Surface Reactions
, 1st ed. (
Springer-Verlag
,
Berlin
,
2012
).
67.
Ramírez
,
J.
,
S. K.
Sukumaran
,
B.
Vorselaars
, and
A. E.
Likhtman
, “
Efficient on the fly calculation of time correlation functions in computer simulations
,”
J. Chem. Phys.
133
,
154103
(
2010
).
68.
Evans
,
R. M. L.
,
M.
Tassieri
,
D.
Auhl
, and
T. A.
Waigh
, “
Direct conversion of rheological compliance measurements into storage and loss moduli
,”
Phys. Rev. E
80
,
012501
(
2009
).
69.
de Greef
,
T. F. A.
,
M. M. J.
Smulders
,
M.
Wolffs
,
A. P. H. J.
Schenning
,
R. P.
Sijbesma
, and
E. W.
Meijer
, “
Supramolecular polymerization
,”
Chem. Rev.
109
,
5687
5754
(
2009
).
70.
Koeppel
,
A.
,
P. R.
Laity
, and
C.
Holland
, “
The influence of metal ions on native silk rheology
,”
Acta Biomater.
117
,
204
212
(
2020
).
71.
Saito
,
M.
, and
M.
Matsumoto
, “Simd-oriented fast mersenne twister: A 128-bit pseudorandom number generator,” in Monte Carlo and Quasi-Monte Carlo Methods 2006, edited by A. Keller, S. Heinrich, and H. Niederreiter (Springer, Berlin, 2008), pp. 607–622.