Theoretical rheo-physics of silk: Intermolecular associations reduce the critical specific work for flow-induced crystallisation

Silk is a semi-dilute solution of randomly coiled associating polypeptide chains that crystallise 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 crystallisation. To investigate the specific-energy input needed to reach this criterion in start-up flow, we have coupled a model for the coarse-grained Brownian dynamics of the chain to the stochastic, strain-dependent binding and unbinding of their associations. 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 optimise chain alignment and stretching in different locations along the spinning duct: this delicate mechanism also highlights the potential biomimetic industrial benefits of chemically tunable processing of synthetic association polymers.


I. INTRODUCTION
The manufacturing of both natural and artificial polymer-based fibres relies on flowinduced crystallisation in non-linear rheological conditions [1][2][3][4][5].The energy input required by this process may be significantly reduced through tailored macromolecular interactions, as exemplified by natural silk [6][7][8][9]: This protein, of which the conformation closely resembles a random coil [10], self-assembles in flow in aqueous conditions under energy requirements orders of magnitude lower than its synthetic counterparts [6].It has been hypothesised that flow-induced stretching of the chain disrupts a solvation layer and in turn enables crystallisation to commence.This mechanism was employed to induce crystallisation of synthetic poly-ethylene oxide by flow at similarly low energetic requirements as silk, however, at much higher molecular weight and/or strain rates [11].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 [12][13][14][15][16][17][18][19][20][21].
We previously found in collaboration with Laity and Holland that the silk protein exhibits intermolecular reversible cross-links [7,22].While these associations shift the alignment-tostretch transition to smaller strain rates by replacing the usual Rouse relaxation dynamics for 'sticky Rouse' relaxation [12][13][14][15][16][17][18][19][20][21], this is not the full story, for the Bombyx mori silkworm manages also to generate the opposite effect during pupation: when the silkworm starts spinning the material it actually chemically reduces this relaxation time through the addition of potassium cations [7,22].Intriguingly, the group of Holland discovered that the structural features of the silk fibre is significantly enhanced through a gradient in the pH along the spinning duct, suggesting an exquisitely controlled local rheology [23].While lower pH may induce partial folding of the protein [10], it is also expected to enhance the lifetime of the intermolecular associations, or 'stickers'.As schematically indicated in Fig. 1, we hypothesise first that the initial reduction of the viscosity decreases the specific work needed to align the chains in the direction of the flow field, while the subsequent increase of the sticker lifetime downstream promotes the stretching of chain segments, which in turn induces crystallisation.
Crucially, inspired by our previous finding that broad conformational distributions emerge due to the stochastic nature of binding and unbinding stickers [8,9], we therefore secondly hypothesise that crystallisation may be initiated by reaching a critical concentration of highly stretched chain segments.This would require significantly less energy input than for 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 hypothesise that the experimentally observed gradients in the pH serve to control the lifetime of intermolecular cross-links locally within the process, which in turn minimises the energetic requirements to deform the network and induce fibre crystallisation.
stretching the entire population of chain segments.
To theoretically investigate this hypothesis, we seek to simulate the conformational response of associating (i.e, 'sticky') polymers in various non-linear flow conditions, and calculate the specific critical work needed to achieve both alignment (quantified using a nematic order parameter) and finite chain stretch in a critical fraction at time t s (this integral is taken in the Lagrangian comoving frame of a fluid element).In Eq. ( 1), κ is the strain-rate tensor, and σ the stress response.Because the stretch distribution is extremely broad for associating polymers in strong flow, it is required to capture how the entire distribution of conformational properties (rather than only the ensemble averages [24,25]) evolves with time.There are now several such classes of polymer rheology known, for which mean conformations are very poor estimates of the distribution.Other examples are linear polymers in shear [26][27][28] and ring polymers in extension [29][30][31][32].
The strain-rate tensor in Eq. ( 1) may in principle be predicted using molecular dynamics (MD) simulations, where the interactions between stickers are modelled using attractive potentials.At this level of computational detail, sticker dissociation may occur following attempts to escape the attractive potential through molecular vibrations [31].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 [12]: In a coarse-grained picture, this sticker lifetime is an elementary rather than an emergent timescale.This allows a description of the problem in terms of the dynamics of a single chain in a crowded environment [8,9,[33][34][35], an approach similar to the modelling of entangled polymers through slip-link and slip-spring models [32,33,[36][37][38][39][40][41], where the generation and destruction of entanglements are modelled as elementary processes.
While there is no unique way of formulating a coarse-grained single-chain model [42], all variants of bead-spring, slip-link and slip-spring models can be written in the general form where i is a chain segment at position R i that is thermally equilibrated at the relevant time scales [43].We will refer to this chain segment as a 'node' of an elastic network, which may represent a non-sticky segment of a chain (a regular 'bead'), a segment with a reversible association (a 'sticker'), or it may be an entangled segment (a 'slip-link' or a 'slipspring').Which of these representations is invoked manifests itself in the definition of the friction coefficient, ζ i , the (friction-dependent) thermal forces, F thermal,i , and the network forces, F network,i .For instance, in classes of models where nodes move affinely with the flow field, the network force exactly cancels the sum of the intramolecular force set by the chain conformation, F network,i = −F intra,i − F thermal,i .This 'rigid-network approximation' is tacitly invoked in the slip-link model by Hua and Schieber [39] and in our recently published model for sticky-polymers in a rigid network [8,9]).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 [40].In the present work, we will account for the compliance experienced by the stickers in a reversible network.
In the following, in Section II A we develop 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 non-spatially-explicit multi-chain approach.In Section II B, we present a two-state master equation that generates analytical predictions of the impact of sticker

A. Brownian dynamics of Sticky Polymers in Flow
In this section we will present a coarse-grained description of associating polymers, where the dynamics of sticker opening and closing will depend on the number of open and closed stickers in a non-spatially-explicit collection of chains.Any linear polymer that consists of N monomers may be discretised using a number of nodes, N nodes , see Fig.At this level of coarse-graining, the friction of the nodes is given by with ζ 0 the monomeric friction.The assumption that the dangling chain ends are relaxed may be released by explicitly modelling the position of the chain ends and setting ∆s i ≡ 0 at i = 0 and at i = N nodes [44].
The equilibrium structure of the chain in quiescent conditions is determined by the endto-end distance of the substrands, , where the stretch ratio λ obeys the equilibrium distribution This distribution emerges as a consequence of the intramolecular and thermal forces in Eq. (2).
In order to derive the intramolecular spring forces, we consider the spring force of the entire chain of N monomers with a mean stretch ratio of unity where approximately captures the anharmonicity of the spring force due to the finite extensibility of the substrand [45].For the substrands i the terms in this equation have to be renormalised as F intra → F intra,i , and N → ∆s i N and λ max → ∆s 1/2 i λ max ≡ λ max,i .(a short strand has a higher harmonic spring force and a smaller maximum stretch ratio than a long chain.)The direction of the force exerted by spring i on node i is Q i /|Q i |, while the direction of this force acted upon node i + 1 is −Q i /|Q i |.Hence, the net intramolecular force exerted on node i is The thermal force is given by the equipartition theorem F thermal,i,α (t)F thermal,i,β (t ) = 0, for α = β (9) with α, β = x, y, z the Cartesian coordinates and k B T the thermal energy.
The force acted upon the nodes by flow is, provided that our coordinate system moves with the flow field, given by where κ is the strain-rate tensor, which in extension and shear is given by respectively.As the coordinate system moves with the flow field, the spatial quantities of physical interest to calculate are the deformation of the individual substrands using which we recursively obtain the drift of the nodes as The value of the first entry, ∂R 1 /∂t is adjusted to fix the centre of mass of the chain (this assumes that the centre 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.3), and the thermal forces are given by the equipartition theorem Eq. ( 10) as before, but with this modified friction coefficient.The network forces are now given by Hence, the paired stickers i and j have an identical friction coefficient and experience the same net force F A intra,i + F B intra,j + F A thermal,i (where F A thermal,i = F B thermal,j ).Crucially to forced sticker dissociation, the net force that acts on the closed sticker pair is which we assume, as in other cases of forces temporary unbinding, lowers the activation energy for sticker dissociation as E act = E 0 act − F stic with E 0 act the activation energy in quiescent conditions and the typical length scale associated with sticker dissociation [9].
We remark that the (apparent) activation energy obtained from experiments using the Arrhenius-type equation [17] τ s = ν exp(E act /k B T ), 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 [7,46] or due to the mixing of various mechanisms of sticker opening and closing, such as bondswapping [9].For now, we assume a well-defined pairwise association-mechanism, for which the rate of sticker opening is k open = τ −1 s and the rate of s (1−p)/p with p the equilibrium fraction of closed stickers [9,12].
In our non-equilibrium simulations, the actual number of open stickers fluctuates and k close is rescaled accordingly.We implement 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 ( respectively.In the Appendix we discuss the algorithm that couples a kMC scheme for sticker opening and closing to the Brownian dynamics of the chain conformation.This coupling is independent of the assumed mechanism of pairwise association-dissociation, and may in principle be used to also describe pairwise bondswapping and the formation of larger sticker aggregates if appropriate kMC algorithms are applied [47][48][49].

B. Approximate theory in transient extensional flow: Two-state model
The dynamics of sticky polymers is complicated by the fact that a polymer with Z s stickers can be in 2 Zs different states, as each individual sticker can be either open or closed.
An instructive simple case is a chain with Z s = 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).We have previously shown that this simple 'sticky dumbbell' or 'two-state' model possesses the double advantage of capturing the essential qualitative physics of the more complex cases and offering some analytical tractability if the chain stretch remains below the finite extensibility limit.In particular, this theory predicts steadystate power-law stretching distributions in extensional flow below the stretching transition of chains with various numbers of stickers [8], Before embarking on simulations with the full multi-sticker and finitely-extensible chain model described in the last section, therefore, we will anticipate the phenomenology with the simple two-state model.In particular, we extend our previous analysis to describe the dynamics by which the two-state stretch distribution evolves in start-up flow.
The starting point is to consider a chain in two states where the stickers are either closed (state 0) or open (state 1).The opening rate is κ − and the closing rate is κ + .The time development of the probability distribution of the stretch ratio is described by [8] ∂P To approximate this equation analytically, we reduce it to a system of first-order ordinary differential equations by first introducing the variable y ≡ ln λ, which using the chain rule The non-linear contributions can then be omitted by considering the limit of large stretches, i.e., we approximate e −y ≈ 0, which is equivalent to λ 1.
In steady state, the left-hand side of the equation is zero, and the solution is given by[8] with c a normalisation constant (its value can in principle be determined by releasing the approximation e −y ≈ 0), and with the exponent of the power-law distribution given in terms of physical parameters by (this is one of the eigenvalues of Eq. ( 20) and Eq. ( 21); the other eigenvalue is −1 and is unphysical as a distribution of the form λ −1 cannot be normalised.)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 at (1 − p) ετ R = 1, the exponent becomes ν = −1 and the stretch distribution can no longer be normalised.Depending on the sticker lifetime, at smaller strain rates the exponent may reach a value ν = −2 if the 'sticky Weissenberg number' (1 − p) ε 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 (1 − p) ε becomes larger than 1/2 [8], 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 equilibrates 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 a the Dirac-delta distribution (justified by the very wide long-time distribution), at λ * (0), from which it can be easily seen that the distributions shift to higher stretches for the closed state, P 0 (t, λ) = c 0 δ(λ − λ * (0) exp [ εt]) and retract to smaller stretches for the open state ).This suggests that the 'front', λ * (t), of any distribution with finite P 0 , shifts exponentially in time to higher values through To develop an analytic approximation for the long-time limiting approach to the steadystate, the stretch distribution is equilibrated below the stretch ratio of the front λ * (t 0 ) at time t 0 , while the tail of the distribution is empty for larger stretch ratios.Assuming that λ * (t 0 ) 1, the equilibrated part of the distribution is virtually independent of time t; this provides a fixed-boundary condition at λ * (t 0 ).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 non-zero at the time τ (λ).The inverse function λ(τ ) is then the trajectory of the 'front' of the distribution.In the Appendix, we detail the Laplace transform of Eqs.(20-21) with the boundary condition in this long-time regime, which as a solution gives with ν the 'equilibrium stretch exponent' in Eq. ( 24) and with the 'dynamic stretch exponent', which controls the growth of the front of the distribution as Upon approaching the stretch transition Wi sticky = 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 (1 − p)τ s τ R at the strain rate Wi sticky = 1/2 where the variance of the stretch diverges (see discussion under Eq.( 24)), we find ν ≈ 1/3, which indicates that the late-stage measure of the front is shifted from the early-stage measure for the outliers by a factor 3. We have also checked that the moving front is narrow for small strain rates Wi sticky < 1/2.In the Appendix, we provide more analytical analysis of the two-state model to estimate the width of the front (relative to its extent) as ∆ rel ∝ pWi 0 Wi sticky /(1 − sticky), where ∆ ≈ (∂[P (λ, t)/P eq (λ, ∞)]/∂ ln λ) −1 / ln λ.As we show in the Appendix, typically this width is ∆ rel 1, and the front of the distribution is narrow even close to the stretch transition.

A. Linear dynamics
We have verified the physics of our model in the linear viscoelastic regime by first simulating non-sticky 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 ∆s i = 1/(M + 1) for all i).
Fig. 3 shows that the choice of the number of beads has a negligible influence on the time evolution of the mean-square displacement, MSD, and is in all cases in agreement with the theoretical prediction where the diffusivity, D, is for non-sticky polymers given by the bare Rouse diffusivity Moreover, the inset of Fig. 3 shows that also the end-to-end-distance is distributed according to the physical equilibrium result of Eq. ( 4).For times shorter than the Rouse time of strands between stickers, i.e., for the dynamics of a sticky polymer are governed by the same Rouse diffusion as non-sticky chains, see Fig. 5(a).For later times than that, the motion of the polymer temporarily arrests within the rigid-network approximation (i.e., were the closed stickers are unable to diffuse, see Introduction), while it is continuous in the case the network is elastically compliant.As the mean-square displacement increases, a larger substrand of the chain itself is engaged in the motion (requiring as in normal Rouse-chain diffusion longer Rouse wavelengths to relax), while in the case of sticky polymers also a larger portion of the surrounding network is engaged in the motion (which increases the friction).Consequently, the MSD will reach also a plateau in the compliant-network case, accounting for inter-chain interactions explicitly,  There appears a non-trivial mapping of equivalent parameters between this compliantnetwork model and analytic Rouse-chain results.For low frequencies, the dynamic moduli obtained within the rigid-network approximation correspond to the sticky-Rouse model albeit with Z s = 3 in the analytical model, while the stochastic model has Z s = 10.A similar discrepancy was previously discussed in Ref. [35].More importantly, we see an even stronger deviation if the network is compliant instead of rigid: The plateau modulus decreases and the (apparent) plateau in G narrows down.This is in agreement with Fig. 5(a), which suggests that the apparent (rheological) number of stickers becomes smaller than the actual number of stickers counted from in the polymer sequence.Indeed, the number of stickers of the silk protein determined from its linear rheology was smaller by a similar factor than the number of stickers counted from the sequence of amino acids [7].Fig. 5(a) shows that the analytical model of Ref. [12] describes our simulations well 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 see Ref. [12].Fig. 5(b) shows that upon releasing the rigid-network approximation this scaling behaviour persists, but rescaled with a prefactor ≈ 4. While this scaling regime is reached for the chains with more than 5 stickers (i.e., above the percolation threshold for network formation), this is not the case for the chains with 2 stickers.Within the rigidnetwork 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 rigidnetwork approximation, the chains cluster into linear 'supramolecular' dimers, trimers, etc.
through an exponentially decaying cluster-size distribution [52], 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.

B. Non-Linear Dynamics: Steady State
While ordinary Gaussian polymer melts and solutions of narrow molecular-weight distribution exhibit correspondingly narrow stretch distributions in steady extensional flow, in shear their conformational distributions become wide due to dynamic stretching, tumbling and recoiling of the chains.In this section, we investigate what these distributions look like when stickers are implemented, where we take finite extensibility of the chains and elastic compliance of the surrounding network into account.We illustrate a striking example of how the stickers affect the steady-state chain conformations in Fig. 6.While intuitively stickers help to stretch the chains, we find that above the bare stretch transition ετ R = 2, where non-sticky chains are all stretched, the sticky chains exhibit an enormous dispersity in the chain stretch, as well as occasional hairpin conformations.These are cause by the stochastic binding and unbinding of stickers, where the network forces may occasionally act in the opposite direction of the flow field.This mechanism rensembles the tumbling of non-sticky polymers in shear flow, but has an entirely different origin [26][27][28].
To go beyond these qualitative observations, we have calculated the stretch distributions of these chains with Z s = 0 and Z s = 5 at various extension and flow rates in Fig. 7.In order to enable a quantitative comparison with the analytic results on the two-state model, we have also included results for a chain with Z s = 2.
Eq. (4) shows that in all cases the equilibrium stretch distribution for zero-flow conditions (black curve) is approached for small strain rates.For non-sticky chains (Z s = 0), a broad stretch distribution with a cutoff set by λ max emerges in shear due to the dynamic stretching, tumbling and re-collapsing of the chains.In extensional flow, the distribution broadens only within a narrow range of strain rates 0.9 < ετ R < 1.1 around the stretch transition.Beyond the stretch transition, the stretch distribution is narrow and Gaussian and approaches λ max with an increasing strain rate.This behaviour qualitatively changes upon incorporating stickers.
Polymers with 2 or 5 stickers form supramolecular branches below and above the percolation threshold for network formation, respectively.Although the steady-state stretch distributions in shear are similar to those of the non-sticky chains, in extensional flow the distributions of sticky polymers are remarkably different.In contrast to the non-sticky polymers, the sticky polymers show broad stretch distributions in steady-state extensional flow over a broad range of flow rates.We have observed this behaviour previously in simulations where the chains were pre-aligned in the flow-field and were we invoked the rigid-network approximation [9].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.
We also find that the large fluctuations in stretch below the formal stretch transition carry over [9].(Thestretch transition is defined at the condition τ SR ε = 1, with the sticky Rouse time obtained from the sticky-Rouse diffusivity of Fig. 5 as we find that for small strain rates and large stretch ratios λ the stretch distribution has a  (a,c,e) and shear (b,d,f) rates for a linear unentangled, non-sticky (Z s = 0) and sticky (Z s = 2 and Z s = 5) polymers with a maximum stretch ratio λ max = 75.The black curve represents the quiescent contour-length fluctuations, given by Eq. ( 4).
power-law tail (see Eq. ( 17)) of which the width is set by a ε-dependent stretch exponent ν (see Section II B).We have determined the stretch exponent from the distributions of the chains with 2 and 5 stickers (we discuss the numerical method in the Appendix) in extensional flow with and without the rigid-network approximation and finite extensibility, and plot these against the strain rate in Fig. 8.As anticipated, we were able to map the stretch exponent of the chain with two stickers onto the analytical result in Eq. ( 24).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. ( 24) is replaced by τ s /2, which gives We further confirmed this agreement by also simulating a chain with two stickers and p = 0.5, and comparing it to this theory.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.Two-State Model Z s = 2, p = 0.9 (compliant) Z s = 2, p = 0.9 (rigid) Z s = 5, p = 0.9 (compliant) Z s = 5, p = 0.9 (rigid) FIG. 8. Stretch exponent ν of the power-law tail of the stretch P ∝ λ ν for simulations of polymers with Z s = 2 (blue symbols) and 5 stickers (red symbols) with p = 0.9, within the rigidnetwork approximation (closed symbols) and using elastic compliance and finite chain extensibility (open symbols).The solid curve is given by the sticky dumbbell model in Eq. (33).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. 7.
For the chains with 2 and 5 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.Fig. 8 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. 3.

C. Non-Linear Dynamics: Transients
In our pursuit to understand the flow-induced crystallisation of associating polymers such as the silk protein, we are interested in capturing the macroscopically observable stresses in start-up flow, and to interpret crystallisation rates in terms of the chain conformations that underlie these stresses.To address these challenges, in this section we will present the timedependent rate-normalised transient shear stress, σ xy / γ, and extensional stress (σ yy −σ rr )/ ε, with the stress tensor (in units of energy per molecule) given by extension Z s = 0; ετ R = 2 Focussing first on the results for non-sticky chains with a finite extensibility λ max = 75 in Fig. 9, we reproduce the well-known qualitative features of their stress transient [42]: For small Weissenberg numbers, ετ R < 1, ετ γ < 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 re-collapsing of stretched chains, and in extensional flow there is a sharp increase in the extensional stress until a plateau due to the finite extensbility of the chains is reached.Because of the thermal fluctuations and dispersity in the initial chain conformations, Fig.
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 .
This sharp peak in the stretch distribution is a fingerprint for non-sticky linear polymers in extensional flow, and will not be visible for the sticky polymers, as we we will now show for Z s = 5.We have modelled the physics of the stickers using the same description as in our previous work on chains that are pre-aligned in the flow field [9]: the fraction of closed stickers is p = 0.9, the sticker lifetime is τ s = 10τ R , and the activation energy is E act = 8k B T .
Sticker dissociation is described using a dissociation length of = 1 nm; the intramolecular force to obtain the change in activation energy is obtained by assuming a chain of N = 5525 monomers and Kuhn length 0.4 nm).We plot the resulting start-up stresses and stretch distributions in Fig. 10.Qualitatively, we find similar shear and extensional viscosities as in the non-sticky 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 non-sticky 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, as expected [8].For strain rates 0.3 < ετ s < 0.5 these fluctuations fill up a power-law distribution whose stretch exponent is depicted in Fig. 8.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(b).At early times, the stretch distribution is equilibrated according to Eq. ( 4).As time proceeds, a 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. ( 28) of Section II B.

D. Critical specific work
Now that we have captured how stickers lead to broad stretch distributions, we will investigate how these distributions affect the critical work for flow-induced crystallisation (FIC).The usual predictor for FIC is the Kuhn segment nematic order parameter, P 2,K ∈ [0, 1].If P 2,K → 1 (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.We argue that FIC is achieved above a critical fraction of chain segments, P s , that are stretched above some critical stretch λ s,i ≡ L s λ max,i , with λ max i = λ max √ ∆ i the maximum stretch ratio of a segment (note that in this definition we assume ∆ i to be constant for all chain segments i, see Section II A), and L s ∈ [0, 1] a critical stretch parameter, which may be viewed as a proxy for chain stretch at the Kuhn length.This 2-parameter family of critical criteria is formally captured using Lsλ max,i P (λ, t s )dλ ≥ P s .
In our simulations, we have monitored the maximum stretch of a single chain segment in an ensemble of 50 chains that each have 10 chain segments (i.e., ∆s i = 1/10 and P s = 1/500; we averaged the results over 5 simulations with different random number seeds).At this level of coarse graining, the highest resolution of nematic chain alignment is captured using the nematic order parameter P 2,s ∈ [0, 1], which is the largest eigenvalue of the nematic order tensor P 2,s = (3 uu − 1)/2, where u is the unit vector tangential to the backbone of the chain (note that this is an overestimate of the Kuhn segment nematic order, i.e., P 2,s > P 2,K ).In Fig. 11, we have calculated the critical specific work, W , as given in Eq. ( 1), needed to achieve values of P 2,s and L s in the range from 0 to 1 for non-sticky (Z s = 0) and sticky (Z s = 5) chains for various shear and extensional rates.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.This is in contrast to the typical behaviour in experiments on non-associating polymers (e.g., the flow-induced crystallisation of HDPE [6]), 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 behaviour where the nematic order parameter (red) is typically larger than the stretching parameter (blue): with an increasing specific work the chains first align and then stretch.
This behaviour 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 both in 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.The stretching parameter (blue) follows a sharp sigmoidal dependence against the critical work, and rapidly outgrows the alignment parameter (red).This supports out hypothesis that flow-induced crystallisation may be achieved a small critical specific work by exploiting the stochastic nature of associating polymers.
We have investigated how the critical specific work, W , and the timescale, t s , at which the critical specific work is reached depend on the strain rate in Fig. 12.The left panel shows that the timescale scales as t s ∝ Wi −1 , 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 Figure 2 of the work by Holland et al. on silk [6]; more dedicated research is needed to investigate this observation.
The right panel of Fig. 12 shows that the critical specific work shifts to higher values for an increasing strain rate, which is in agreement with experimental observations on the flow-induced crystallisation of synthetic polymers and silk [6].
We explain the increase of the critical specific work with an increasing strain rate in terms of the two-state model that we introduced in the Theory section.We argue that the stress crystallisation and (ii) can be modelled 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 modelling 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 [29][30][31][32]).
Regarding silk rheology, we have confirmed our hypothesis that the stickers between chains may reduce the critical specific work to induce flow-induced crystallisation under reasonable assumptions for critical crystallisation criteria.While intuitively this is achieved by the presence of stickers shifting the stretch transition to smaller strain rates (which is of course highly relevant to the processing conditions both in animals and in the industry), the actual reason for creating a minimum in critical work is more subtle: The stickers fix a network at time scales below the sticker dissociation time, and in fact hinder the alignment of the chains in the direction of the flow field.We found, using the measure L s , that owing to the stochastic nature of the binding and unbinding of the stickers a significant portion of the chain segments can be stretched to a large extend, while the global alignment, measured by the nematic order parameter P 2,s , remains small.Assuming that crystallisation commences above a critical concentration of stretched chain segments, rather than when the entire population is stretched, we found a significant reduction of the critical specific work is possible.
Focussing on our finding that chain alignment at low, non-stretching, 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 argue 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 of the potassium concentration through a decreasing lifetime of calcium bridges (stickers) [7,22].This, as we can now interpret as a mechanism to ease chain alignment in flow.Intriguingly, downstream the spinning duct the acidity increases [23], 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 crystallisation [6,11].or not at all according to a kinetic Monte Carlo (kMC) scheme [47][48][49].In every kMC step, the rate at which any opening or closing event may occur is calculated as W T = W a + W d , with the sum of closing rates and with u ∈ (0, 1] a uniform random number (our code uses random numbers using the opensource SFMT library [53]).If ∆t 2 exceeds the time span ∆t 1 , no opening or closing events From the boundary condition, we find that the coefficients must be of the form c ± i ∝ 1/s.As the '+' solution leads to a non-normalisable solution, however, c + i = 0, and the solution is P0 (s, λ) = c s (λ/λ * ) ν − −(s/ ε)ν − 1 2 (s/ ε) 2 ν +O(s 3 ) , (56) Finally, after taking the inverse Laplace transform, we have P 0 (t, λ) = c λ λ * (0) νeq Θ(ν ln λ/λ * − εt) (58) Hence, the exponentially extending front of the distribution is located at the stretch ratio We have checked the validity of our interpretation of a narrow moving-front by also calculating the width of this front.To do this, we consider the relaxation function f (t) = P (y, t)/P eq (y) with again y = ln λ, and P and P eq the transient and steady-state stretch distributions, respectively.A narrow front that reaches y at time τ and equilibrates at time τ + ∆ may be approximated by 0, for t < τ (t − τ )/∆, for τ ≤ t < τ + ∆ 1, for, t ≥ τ + ∆. (61) The Laplace transform of this function is We compare this result to the solution of the two-state model in Eq. ( 53) through a secondorder Taylor expansion of the exponential terms

FIG. 1 .
FIG. 1. Schematic representation of the gland of a B. mori silkworm, where extensional flow opening and closing on both the steady-state and transient stretch distributions of the chains, which enables us to interpret our simulated data in Section III.By first mapping the results in the linear flow regime to the analytic sticky-reptation (SR) model, in Section III A we discuss how the stochastic nature of sticker opening and closing and the elastic compliance affects the linear rheological data.Then, in Section III B we show how a broad steadystate distribution of chain conformations emerges in strongly non-linear flows of shear and extension.By simulating the transient emergence of these distributions in start-up flow in Section III C, we show that the stickers initially hamper the collective alignments of the chains in mildly non-linear aligning flows, but facilitates the emergence of stretched outliers.In Section III D we discuss how these outliers may reduce the critical specific work for flow-induced crystallisation.In the conclusions in Section 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 flow-induced crystallisation of artificial materials.

FIG. 3 .
FIG.3.Mean-square displacement, MSD, of the centre of mass of a non-sticky polymer against time (main panel) and the time-averaged end-to-end length (R e ) distribution (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.

= 2 FIG. 4 .
FIG. 4. Linear rheology of a sticky chain with Z s = 10, p = 0.9, τ s = 200τ R within the rigidnetwork approximation (open symbols) and with this approximation released (closed symbols).(a) Mean-square displacement MSD against time.(b) Loss, G , and storage, G , moduli of the chain against the frequency, ω.The moduli are scaled with respect to the plateau modulus G 0 s within the rigid-network approximation.The curves are the SR model of Eq. (31) with Z s = 2 and 3 as the apparent number of stickers per chain.

Z s = 2 Z s = 5 Z s = 10 Z s = 20 FIG. 5 .
FIG. 5.Sticky Rouse diffusivity, D SR , against the sticker lifetime, τ s for chains with Z s = 2, 5, 10, 20 stickers within a rigid network (a) and a compliant one (b).The symbols are our simulation results, and the curves represents the sticky Rouse model in Ref.[12].The units are given in terms of the bare Rouse diffusivity D R and the bare Rouse time, τ R .

FIG. 6 .
FIG. 6. Representation of simulated chain conformations in extensional flow for ετ R = 2 for non-sticky (a) and sticky (b) polymers.While the variations in stretch are narrow for non-sticky polymers, these variations are broad for the sticky polymers: when a sticker in a retracting chain segment binds to a neighbouring chain segment, this may disrupt the neighbouring chain.The scale bar represents approximately a length 50R e , which is 65% of the fully extended chain.

FIG. 9 .
FIG. 9. (a) Simulated rate-normalised transient extensional and shear stresses of non-sticky polymers with a maximum stretch ratio λ max = 75 (averaged over 50 polymers).(b) Transient stretch distribution in extensional flow for the chain in (a) at strain rate ˙ τ R = 2.
FIG. 10.(a) Rate-normalised transient extensional and shear stresses of sticky polymers with a maximum stretch ratio λ max = 75 (averaged over 50 polymers).The sticky polymers have Z s = 5 stickers with lifetime τ s = 10τ R and closed fraction p = 0.9.Below the stretch transition large fluctuations are visible.(b) Transient stretch distribution in extensional flow for the chain in (a) at strain rate ˙ τ s = 2; the error bars represent half of the standard error of the mean.

FIG. 11 .
FIG.11.Nematic order parameter, P 2,s and characteristic stretch ratio, L s , against the specific work (see main text) for sticky (red) and non-sticky (blue) polymers in shear (left) and extensional (right) flow.The symbols are obtained from simulations with various strain rates; the curves serve as a guide to the eye.

FIG. 13 .
FIG.13.Flow chart of the algorithm so simulate the conformational dynamics of sticky polymers and the dynamics of sticker association and dissociation.
the sum of dissociation rates, where k d,q differs for the different sticker pairs due to dispersity in chain tension.In these expressions, N open and N closed are the number of open and closed stickers, respectively; N open (N open − 1)/2 is the total number of possible associations, and the index q sums over all N closed /2 pairs of closed stickers.Using this sum of rates, the time ∆t 2 at which the first opening or closing event occurs is ∆t 2 = − 1 W T ln(u),

2
+ ∆τ ) (ν / ε2 ) ln y s 2 .(63)determine the cutoff by minimising χ 2 (a, b, i 0 ) ≡ 1 N data + 1 − i 0 − N par a, b and i 0 (note that x i 0 = x cutoff ); σ i is the uncertainty on the simulated y data.Here, we set b = 1 fixed and the number of free parameters N par = 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 N par = 2.