Self-assembling multiblock amphiphiles:molecular design, supramolecular structure, and mechanical properties

Self-assembling multiblock amphiphiles:
molecular design, supramolecular structure, and mechanical properties

Hamed Mortazavi and Cornelis Storm Department of Applied Physics and Institute for Complex Molecular Systems, Eindhoven University of Technology, P.O. Box 513, 5600 MB Eindhoven, The Netherlands
July 5, 2019

We perform off-lattice, canonical ensemble molecular dynamics simulations of the self-assembly of long segmented copolymers consisting of alternating, tunably attractive and hydrophobic binder domains, connected by hydrophilic linker chains whose length may be separately controlled. In such systems, the molecular design of the molecule directly determines the balance between energetic and entropic tendencies. We determine the structural phase diagram of this system, which shows collapsed states (dominated by the attractive linkers’ energies), swollen states (dominated by the random coil linkers’ entropies) as well as intermediate network hydrogel phases, where the long molecules exhibit partial collapse to a single molecule network state. We present an analysis of the connectivity and spatial structure of this network phase, and relate its basic topology to mechanical properties, using a modified rubber elasticity model. We find that it is possible to optimize the mechanical performance by an appropriate choice of molecular design, which may point the way to novel synthetics that make optimal mechanical use of constituent polymers.

82.20.Wt, 82.35.Jk, 81.16.Fg

I Introduction

Biological materials consist mostly of networks of filamentous fibers. A prime example is the cytoskeleton, which - among its many other functions - governs mechanical stability, the response to stress, and the motility of the cell. Its mechanical properties, particularly in the non-linear regime are uniquely linked to its molecular architecture: an open meshworks of interconnected semi-flexible fibers Bausch and Kroy (2006); Broedersz and MacKintosh (2014). Mimicking these mechanical properties in a synthetic material would offer important new possibilities both for cell research and health technologies such as drug delivery or cell growth, as well as for novel bio-insipred and biobased performance materials. In this paper, we investigate a self-assembly route to recreating the topology of such materials in synthetic associative polymers, and ask the question how these should be designed in order to spontaneously produce networks that present optimal mechanical performance. To this end we study, numerically, a self-assembling system of block copolymers that forms an effective cross-linked network.

Figure 1: We simulate long copolymers consisting of alternating binder (yellow) and linker (blue) blocks. The binder blocks consist of repeating PEG units, the linker blocks are bisurea motifs. In synthesis, the length of both these blocks is separately tuneable. In our simulations, we coarse grain these molecules into beads representing many atoms at once. The molecular design is fixed by specifying values for , and .

The self assembling system we focus on here is a multiblock amphiphile, consisting of many regular repeats of alternating hydrophilic and hydrophobic blocks. Previous work, which we will briefly review in a moment, has identified such molecules as quite promising candidates for precise regulation of self-organized architecture, but so far experimental realizations have largely been lacking. Polysaccharides with hydrophobic substitutions - methylcellulose, in particular - possess the required alternatingly hydrophilic/hydrophobic backbone structure, but in more or less randomly distributed blocks Sarkar (1979). Recently, however, a novel molecule was synthesized and studied Pawar et al. (2012) that does allow precise control over the spatial arrangement of the alternating blocks, and in fact also over the attractiveness of the individual hydrophobic blocks. Its basic building blocks (shown in Fig. 1) are macromolecules consisting of alternating binder (attractive) and linker blocks, which are, respectively, bis-urea hydrophobic motifs (essentially identical to those discussed in Koenigs et al. (2014)), and polyethylene glycol (PEG) hydrophilic chains Pawar et al. (2012). The objective of this work is to provide an initial survey, based upon Molecular Dynamics simulation and elementary rubber elastic theory, of the expected microstructure of the self-assembled states of this molecule, as well as the mechanical properties of its gel/network phase.

In our simulations, we coarse grain many repeat units into beads that interact via standard hard-core repulsive and - for the binder blocks - attractive potentials.

In aqeous solvent the amphiphilic design of these molecules leads to aggregation and causes the long molecules to partially collapse into clusters that are connected via one or more links. Clearly, if the attraction is too weak the chains will adopt random coil configurations, while if the attraction is too large the molecules will collapse on themselves forming compact globular aggregates. In between, as we will show, there is an intermediate regime where the system spontaneously aggregates into a network phase of connected clusters linked to each other. Here we present simulations in which we control the linker length and the interaction potential between the binders to tune the relative importance of energy and entropy in this self-assembly process, and characterize the resultant network in terms of cluster size and connectivity. We summarize these findings in a phase diagram, identifying the experimentally most promising regime for recovering biopolymer network-like architectures in self-assembling copolymers. We then analyze the mechanical properties of these networks, first in the context of a rubber-elastic theory for networks of variable functionality, then in direct numerical simulations. We are not the first to consider models of self assembling multiblock copolymers: Rooted in extensive work on the assembly and rheology of amphiphilic block copolymers (see, e.g., Alexandridis and Lindman (2000) and references therein) more recent work has also explicitly considered the multiblock design. Early scaling theory by Halperin Halperin (1992) anticipated the emergence of flowerlike micellar structures, either single or multiple such structures connected to each other in the manner of pearls on a string, depending on the quality of the solvent. Subsequent efforts Zhang et al. (2003); Kawata et al. (2007) refined the conditions for formation and stability of the flowerlike micelles. Later numerical work Tanaka and Koga (2000); Tanaka (2002) confirmed that indeed such structures may form in dilute systems. In Hugouvieux et al. (2009, 2011), the Monte Carlo (MC) technique on a lattice model was applied on the system to establish an equilibrium phase diagram, which features characteristic phases of isloated micelles, connected micelles, but also laminar and tubular phases. Similar results were reported in Tanaka and Koga (2000); Cooke and Williams (2003) in yet more MC simulations. More recently, the first dynamical simulations (Brownian Dynamics, BD) were presented, varying systematically the solvent quality. Again, the familiar structural phases were reported Zhang et al. (2013). The picture that emerges is a very rich self-assembling system, which - provided one has some synthetic control over the moelcular architecture - provides a direct control over mesoscopic organization of a hydrogel. A question that has received very little attention, so far, has been the implications of such structure for mechanical performance. Our ultimate ambition is to simulate dynamical rheology for these networks, in conjunction with the self-assembly. We focus on those structural aspects of the hydrogel phase that govern mechanical response. This response is analyzed first in the context of simple rubber elastic models, and then directly, using a computational implementation of the oscillatory rheology protocol on the networks presented here. While - aside from the polysaccharides mentioned above - there are few experimental realizations to compare to, the molecular design we coarse grain here was, in fact, studied in previous experimental work Pawar et al. (2012), which revealed a tendency to form strong hydrogels, whose rheological properties make them uniquely suited as injectable substrates for drug delivery. However, despite the detailed AFM and SAXS analysis of the self-assembled hydrogels, the question of the precise molecular structure of the gels remained unanswered. In this paper, we use molecular dynamical simulations to address this question, which enables us to determine a relation between molecular design, gel structure, and mechanical properties that breaks ground for further rational design in these materials.

This paper is organized as follows: Section II presents the experimental model system, the modeling environment and approach we have adopted, and the relevant interaction potential we employ. In Section III the simulation protocol in coarse-graining the model system is presented. In Section IV, we present data on the phase behavior resulting from the self-assembly process. Section V analyzes the topology of the network phase and presents results on the size and distribution of the connectivity of clusters. In the remainder of the paper we investigate the mechanical response of these networks, in Section VI, we recall and apply a modified rubber elastic model, in Section VII we measure the rheological response of such networks numerically, using a computational bulk rheology approach. This section is followed by qualitative result section. We then relate these findings to experimental observations. We present our conclusions for the relation between molecular design, supramolecular network structure and the mechanical properties in our system.

Ii Model system and coarse grained MD setup

As indicated in Fig. 1, the multiblock copolymer we study here is a polymer consisting of repeating diblock units. Its structure may be denoted as [[PEG][bisurea]]; , and are the repeat numbers of the PEG molecule, the bisurea motif, and the resultant diblock unit respectively. In typical experimental settings, Pawar et al. (2012). The PEG chain is hydrophilic and highly flexible which results in a tendency to swell in water due to entropy maximization. The bis-urea block, in contrast, is stiff (nonswellable) and hydrophobic, resulting in a generic tendency to aggregate in water. This tendency is further enhanced in by the ability of the ureas to hydrogen bond, leading to an increased stability of their aggregates. Thus, the ultimate behavior of this long copolymer is determined by conflicting tendencies imparted by binders and linkers; depending on the dominating feature the system will lean towards collapsed (binder dominated) states or swollen (linker dominated) states. Interestingly, the geometry of the molecule may be used to directy affect this balance between energy and entropy. The longer the linker chain (), the larger the entropic contribution; the longer the bisurea block () the stronger the hydrophobicity/hydrogen bond-driven tendency to collapse.

Figure 2: Graph of the Lennard-Jones (LJ) potentials between binder beads (left). Linker beads, interacting with either other linker beads or binder beads do so with a LJ potential truncated at to approximate hard-core repulsion between beads(right).

In our simulations, we are initially interested in the generic features of this class of systems. We therefore model a copolymer chain as a sequence of permanently attached beads of fixed (and equal) radius . There are two distinct types of beads, and they are distinguished by their mutual interactions (see Fig. 2). All beads repel strongly at distances shorter than , to enforce hard-core repulsion and non-crossing of the copolymer (we address this in more detail below). Binders mutually attract at distances larger than , whereas linkers experience no interactions, either with other linkers or binders, beyond their hard-core radius. As we are initially interested in open, meshworked phases we shall use as our two main control parameters in this system the number of binder beads , the number of linker beads , and the strength of the LJ attraction as measured by the LJ well depth . As the experimental molecules have either or , the physical size of the binder region does not vary much, which is why we choose to fix , reducing the number of free parameters. Simulated coarse grained molecules are denoted as BL, so that, for instance denotes a molecule whose repeating diblock motif consists of 3 binder beads and 24 linker beads. For all simulations reported here, we use 500 repeats of this diblock, so that each molecule in each simulation contains 1500 binder beads.

Figure 3: Evolution of the self-assembly of a network during a typical MD run in the network phase. Starting from a random arrangement, the macromolecule gradually collapses on itself to form clusters of a well defined size, connected by one or more linker chains. Note, that these images are still of single (albeit very large) molecules.

As was previously observed in lattice MC and BD simulations Hugouvieux et al. (2009, 2011); Zhang et al. (2013), the equilibrium phase diagram of such systems is already very rich. At low values of , or very large linker regions, the entropy dominates and swollen phases are expected. For very high values of the LJ attraction, or for short linkers, energy dominates and complete collapse (similar to a polymer in poor solvent) is observed. In the intermediate regime, chains of micelles are reported, as are nonspherical aggregates. The objective of the present paper is to establish whether these structures are also present in off-lattice MD simulations. As will show, most do indeed feature in our simulations, although a large region of the dynamical phase space is occupied by meshlike structures which are best termed flowerlike micellar networks, which appear to be dynamically arrested intermediates that are unable to fully equilibrate to a single collapsed state, see Fig. 3. As these flowerlike micellar networks most closely resemble the biopolymer gels whose topology we aim to recreate, we characterize the connectivity structure and relate it to mechanical response using modified rubber elastic theory.

Our simulations are carried out using the LAMMPS molecular dynamics package Plimpton (1995). We represent the long copolymers using a bead-spring model, under fixed boundary conditions, with two distinct types of beads to represent binder and linker segments. Adjacent beads interact via a harmonic bond potential of the form , and a harmonic bending potential , is applied to each set of three neighboring binder beads (where is the distance between the centers of mass of pairs of beads, is the equilibrium bond length, the size (diameter) of the beads and the angle formed by the two bonds that connect the middle bead to its two adjacent beads). and are the stretching and bending stiffnesses, which are both fixed at fairly high values: and , to enforce inextensibility of the backbone, and inflexibility of the binder region during the simulations. The bending stiffness is set to zero for linker beads, so that the linkers are represented as flexible polymers of molecular weight proportional to . All binder beads that are not first or second nearest neighbors (i.e., those that are part of different binder domains) interact through a full Lennard-Jones (LJ) potential, whose attractive strength measured by :


The interactions involving linker beads (i.e, linker-linker and binder-linker) are purely repulsive, and are modeled with a truncated LJ potential, (see Fig. 2)


In the latter equation, measures the strength of the repulsion, which in our simulations is set equal to ; is still the center-to-center distance between the beads. The repulsive potential is used to ensure that the polymer obeys self-avoidance: We geometrically ensure that the chain cannot cross itself by setting the bead-bead distance to 1.2, and we truncate the LJ potential at . The repulsive LJ potential takes care of the self avoidance of all beads, and to prohibit self-crossing of the polymer we choose an equilibrium bond length of 1.2. This leaves a gap with an equilibrium size of 0.2 between neighbouring beads, and as one may see in Fig. 4, the bond length never fluctuates to lengths sufficient to leave a space through which another segment of polymer may pass.

Figure 4: Histogram of bond length distribution during one run. The spacing between two beads is never sufficiently large to allow other beads through, and hence the chain cannot self cross during the dynamics.

To verify that indeed, this reproduces correct polymer scaling, we turn off the attractive interactions (setting , as well as the bond-bending term) and record the radius of gyration as a function of the length of the polymer. As one may see in Fig. 6, we recover the correct self-avoiding scaling with a Flory exponent of about 0.6. This scaling is unaffected by the actual size of beads. The repulsive LJ potential 2 exists between all beads, its value affects the effective radius of a bead. If we choose it to be too small, the beads become ’soft’ and may overlap slightly. As one may see from Fig. 6, the exponent rises quickly from (the ideal chain value) to its self-avoiding value upon for increasing values of . We fix at 0.01 because at this value, we first see (within our numerical accuracy) the correct exponent of .

Iii Simulation Protocol

In our LAMMPS simulations, each simulation is repeated from random initial conditions for different systems, which are then averaged over. Our protocol is as follows: We begin each simulation by an equilibration process starting from a random initial geometry. To capture the effects of different attractive strengths (representing the different cohesive energies, which rises with in the molecule) and different linker lengths, we simulate twelve types of coarse grained molecules: B3L3, B3L4, B3L5, B3L6, B3L9, B3L12, B3L18, B3L24, B3L45, B3L75, B3L120, B3L150, and each of these is simulated for twelve values of . After an initial randomization, we turn on the attractive interactions and watch the system evolve.

Figure 5: scaling for contour lengths of and beads. The attractive potentials are switched off. The Flory exponent for beads of diameters of and is in agreement with theoretical prediction.
Figure 6: Scaling exponent for different . We have chosen in all our simulations which as this is the value where correct real chain scaling is first observed.

Iv Phase Diagram

What happens next obviously depends on the strength of the attractive interactions, and the geometry of the molecule. For low values of , much like swollen polymers in a good solvent, we see random coil configurations, determined by the entropy of the linkers. Upon increasing , the chains start to self-assemble and form clusters. We see the size of clusters grow with increasing but their growth rate is limited by the entropy of the linkers. For short linkers, collapse to a single cluster at high values is observed. For even longer linkers, however, the average cluster size remains limited since the entropic penalty for packing the linkers into the corona of the aggregate rises strongly with their length. The intermediate regime is where we see a network of flowerlike micellar cores (small aggregates of binders), linked to other such clusters by one or more linker chains - see Fig. 7. This phase diagram is in good agreement with a similar result reported in earlier lattice Monte Carlo simulations, though we obviously cannot reproduce the lattice-induced layered aggregate phase reported there.

A second, more striking difference is that the network phase we see appears to be more prevalent in our MD simulations than in the earlier MC work Hugouvieux et al. (2011, 2009); Zhang et al. (2013). We attribute this to kinetic trapping of the structures: while truly long relaxation might recover collapsed phases as true thermal equilibria, on the timescales that we are able to access in our MD sims the system tends to arrest in long lived, but possibly metastable, stationary states. In the following, we characterize the topology in this network phase.

Figure 7: The phase diagram of the self-assembled states of the macromolecule is divided into three distinct phases: large single clusters (violet), single molecule networks (green) and random coils (purple). The network phase is also loosely divided into three sub-phases: One with a broad distribution of cluster sizes (), ane with much tighter distribution of cluster sizes (), and one with a combination of clusters and isolated single beads ().

V Network structure analysis

Clearly, the region in the phase diagram that most closely resembles that of a crosslinked biopolymer mesh is the network phase. As is known from extensive previous work Wilhelm and Frey (2003); Gardel et al. (2004, 2003); Storm et al. (2005); Huisman et al. (2007); Broedersz and MacKintosh (2014); Bausch and Kroy (2006), the topology of such a network (connectivity, crosslinker density) greatly influence the mechanical response. We will quantify the connectivity and cluster size of the networks that form in the intermediate regime, as these parameters are expected to directly translate into effective parameters for the crosslinked network: We identify each cluster with a network node, whose cohesive energy is determined by the number of binders it contains, and we count the number of connecting links between such network nodes to determine the distribution of functionalities in the network.

Figure 8: Frequency distribution of the step length for B3L7, at . The rapid decay justifies the approximation we will make, which is that all linkers are one step linkers.

Our method to identify the clusters is as follows: After equilibration, we export the MD trajectories of all beads. To identify individual clusters, we sort all binders that are within from each other in a group. is the mean value of the longest distance between two center of mass (the heart-to-heart distance between the central beads of two aligned binder blocks) and the shortest such distance (i.e., the heart-to-heart distance between the central beads of two parallel binder blocks). The average cluster size is then an average over the distribution of number of binders in each group. This computation is averaged over five separate runs. We also calculate the functionality of each cluster, defined as the number of linker chains that bridge to other clusters emanating from one particular cluster. In our simulations we count all of the chains that connect two separate clusters together. However, the length of these bridges can be one or more times the linker length, and the length of the linker determines its effective elasticity. To keep track of this, we also record the step length of bridges, denoting the number of linkers that connect one cluster to the next. We find, however, that the distribution of the step length decays rapidly with , suggesting that the initial approximation that most bridges are one-step is justified - see Fig. 8.

As Fig. 10 shows, the cluster size increases with increasing attractive LJ strength . Basically, we are increasing the importance of energetic over entropic effects. The upper bound of 500 on the cluster size is determined by our choice of system - as stated we have 500 total binder motifs in each simulation (i.e., 1500 binder beads). A system that reaches a cluster size of 500 is therefore fully collapsed. For very low values of , the cluster size is 1 which means no binder is connected to any other binder, and the system is swollen. The cluster size measure is thus able to distinguish between the extremes of swollen random coil phase and collapsed state. For intermediate values of , the system settles into the network phase (see also Fig. 10) and exhibits a finite cluster size, below 500 but significantly greater than one.

The other way in which we can alter the balance between energy and entropy is by increasing the linker length. The longer the highly flexible linker chains become, the more configurational entropy they possess. This is indeed borne out by Fig. 10, which shows that for short linkers, we generally see full collapse but that for larger values of , the system swells and - at a rate that depends on - reverts to the fully random swollen state.

Clearly, the combination of and determines the size of the clusters, as well as the connectivity between them. We may now ask what the projected mechanical properties of the resultant network phase will be. These too will depend on the structure of the network: in a dilute system, such as the ones we consider here, both the random coil state and the collapsed state are expected to have poor mechanical performance, if not complete lack of rigidity - below the overlap concentration, distinct molecules will not entangle to any significant degree in the random coil phase, and most certainly will not do so in the collapsed state. In the intermediate, network phase different molecules will perticipate in a single, connected network of flowerlike micelles. In the following, we estimate the mechanical modulus of such a network based on the structural features we have just discussed.

Figure 9: Average cluster size vs . Networks with large tend to form bigger clusters. This tendency is more pronounced in molecules with shorter linkers.
Figure 10: Average cluster size vs linker length. Increasing the proportion of linker beads increases the contribution of the entropy of the linker to the overall free energy, favoring smaller clusters.
Figure 9: Average cluster size vs . Networks with large tend to form bigger clusters. This tendency is more pronounced in molecules with shorter linkers.

Vi the mechanical response

In the simplest approximation, we shall consider the effective network that emerges in the regime of intermediate and as a rubber with Guassian linkers connecting crosslinkers of varying functionality. In previous work, Yeo et al. Yeo et al. (1981) established that the shear modulus of such a network may be computed according to a modified classical rubber model (see, for instance, Rubinstein and Colby (2007)) , as


In this equation, is the number concentration of elastically effective chains, and g is the so-called front factor, which is to be determined according to the functionality distribution. In principle, counts both the contributions from physical entanglements () and crosslinkers (), but since our network - outside of the network phase and at sufficiently low concentrations - does not possess entanglements that contribute to the modulus we will set . The front factor accounts for changes in density due to the contraction of the network upon crosslinking, as well as for the role of functionality. The combined result, obtained first by Yeo and based upon earlier work by Duiser and Staverman Duiser and Staverman (1965), Graessly Graessley (1975) and Tobolsky Tobolsky et al. (1961) then yields, for a network of average functionality


in this equation is the ratio of the mean squared end-to-end distance of the polymer chains in the crosslinked network to that of the same chains in the uncrosslinked state. Becasue, as we have seen, primarily unit step length bridges occur we will count only those as elastically effective in the following. Using Eq. 4, we may now estimate the shear modulus of the various network states in our simulations, using the functionality distribution to compute , and counting the number density of clusters to determine . Sample results are presented in Fig. 11. The general trends we observe are consitent with what the phase diagram also suggests: As a function of both linker length and , there exists and intermediate regime where cluster sizes are sufficiently large to permit high functionalities (many opportunities for connecting to other clusters), but are not yet so large that the number of potential partner clusters becomes limiting and the clusters become single collapsed entities. For shorter linker lengths (Fig. 11) the dropoff in modulus at higher values of is very pronounced, as the clusters quickly become fully isolated - for larger linker lengths (Fig. 11), the wider reach of every cluster (even when it is already fairly large) allows the system to retain some connectivity even at larger cluster sizes. Similar figures may be drawn for the dependence on the linker length, and in Fig. 12 we collect these into a modulus diagram.

This diagram illustrates what we feel is a crucial point about these networks: more is not always better for increased rigidity. An optimum in modulus exists at intermediate as well as , where the balance between connectivity opportunities and functionality is optimal for overall mecahnical response. It would be interesting to explore whether indeed such an optimal regime exists.

Figure 11: Rubber modulus (according to Eq. 4) vs for B3L12. At large , the high degree of association of the cluster suppresses the potential for connecting to other clusters - large clusters spaced further apart. Thus, the modulus of the network drops off at higher . The solid line is a guide to the eye.
Figure 12: Overview of computed moduli for all molecular architectures studied here. Mechanical performance that is optimal in the sense that it provides the highest modulus arising due to a favorable balance between cluster size and cluster functionality is obtained for B3L75 and .

Finally we note that, obviously, fancier models are available to predict the mechanical properties of networks such as those we consider. Once fully formed, the structure of the hydrogel is no different from that of a telechelic gel, which was shown to be well-described by the classical theories of Flory Flory (1953) and Stockmayer Stockmayer (1944), and whose kinetics of aggregation Wilson et al. (2011) and de-aggregation may be used to assess even their visco-elastic response Skrzeszewska et al. (2010a, b).

Vii Computational rheology

In order to verify the accuracy of rubber elastic models, we now turn to a direct, quantitative characterization of the viscoelastic response of the self-assembled network. Experimentally, establishing these properties is often challenging. There are various rheological methods, including micro- and macroscopic probing methods that can characterize the mechanical behavior of polymer networks. The macroscopic method is called bulk rheology. By this method one can relate the shear deformation (i.e., the strain ), to the shear response (the stress ) in a sample. We copy this protocol and implement it directly in our simulations. Fig. [13] shows a simulation box under oscillatory shear. For purely elastic response this relation is , where is called the shear modulus. For viscoelastic materials, that exhibit viscous responses as well as elastic responses, is usually decomposed into a real and an imaginary part:

Figure 13: In our oscillatory shear simulations the box is deformed (right) in the -plane, and executes a periodic oscillation characterized by a frequency and an amplitude . The network contained inside will accomodate this dynamic deformation differently depending on its ability - or inability - to undergo structural relaxations on the timescale of the applied strain.

In this equation, is the storage modulus, that characterizes the elastic response of the material. The loss modulus quantifies the viscous response. For a purely Hookean solid, is constant, and . For a Newtonian fluid, conversely, and with the dynamic viscosity. In oscillatory rheology, a time-dependent strain of the form


where , amplitude, is the maximum deformation applied in each cycle. The time lag between two sinusoidal signals determines the viscoelastic moduli of the system:


In the above equations, is the phase lag between the stress and strain signals. The phase lag for a viscoelastic system is shown in Fig. [14].

Figure 14: The phase lag between sinusoidal signal of deformation and response determines the elastic and viscous modulus. In this figure an example of is shown. This is an example of an ideal viscoelastic response .

In the case of linear response, this is the entire story - stress and strain always have the same frequency as no higher harmonics are generated. The only degree of freedom that linear response allows is the phase lag , but in linear response the moduli themselves cannot in any way be functions of time, or of the strain amplitude, themselves. Since we will be interested in the non-linear (finite strain) response as well, it is useful to expand the stress response in the higher harmonics, expressing it as a Fourier series


Where is storage modulus, a measure of elastic energy stored in the material. Provided the applied strain is sinusoidal, indicates the magnitude of the harmonic of the stress response. The loss modulus is correlated with the viscous dissipated response and measures the out-of-phase component of the stress. If one has the full signals and , all the moduli can be determined from the complex coefficients of the discrete Fourier transform of the discrete stress time series


where is the period of the applied strain times the sampling frequency. Using the relationships between the Fourier coefficients and , along with knowledge that is the complex conjugate of , each harmonic modulus can then be determined as


With and being the imaginary and real parts of the complex Fourier coefficients respectively. To assess whether or not a particular system, exposed to a strain with period and amplitude can be considered in a regime of linear response, we typically monitor the ratio of the first to the third harmonic modulus (the second is zero, by symmetry).

vii.1 Model and simulation protocol

As stated in section II, in all of simulations so far we have studied a single, very long chain consisting of 500 repeats of the binder+linker motif. In the experiments, chains are typically much shorter. We expect there to be little difference between studying one long chain, or chopping this long chain up into many smaller chains. Now, we verify this explicitly: we compare the results of one repeating chain of binder-linker groups with those obtained for individual chains of binder-linker groups each (this is the typical repeat number for the experimentally used Poly(8kU4U) molecule). Both systems are allowed to self-assemble into a crosslinked network. Each binder is made up of three spherical beads and every linker is made up of beads, and we let vary corresponding to the length of the linker from to beads. To ensure a fixed bond length in our simulations, we choose a large value for the strength of the bond potential . The hydrophobic chains in the experimental system are quite short, roughly , which makes them act like a stiff rod (Pawar et al., 2012). To simulate this we choose a large value for the strength of bending potential . We also compare the results of simulations for harmonic bond potential with FENE bond potential (Bird et al., 1977) by which the nearest-neighbour beads along the chain interact through an anharmonic, finitely extensible, non-linear, elastic potential given by


The parameter values of and prevent entanglement and overlapping of chains (Baschnagel and Varnik, 2005). The FENE potential diverges logarithmically as , providing a finite distance between chain beads. It is used in coarse grained simulations to encode the finite extensibility of real polymers.

According to the above equation, particles closer than interact via a repulsive Lennard-Jones potential whereas beyond the interaction is zero to ensure real chain scaling - See also section II. We study the effect of temperature and density on the mechanical properties of the system and compare them with the experimental results. To simulate experiments at different concentrations, we keep the number of the beads in simulation constant, but change the size of the simulation box to decrease or increase the concentration. To do so, we perform NPT simulations: here, we dial in the pressure to control the box volume. We note, that this may induce other effects besides those of concentration alone, as we necessarily have to go to elevated pressures to realize higher concentrations. We have not studied these prestress effects separately. We choose three different initial pressures to work with: , and (molecular dynamics units111The molecular dynamic units and their conversion to real units are discussed in Appendix A). Our procedure for the equilibration and determining the pressure and volume of the system is as follows:

Figure 15: Left: the system as it is started off, from a random initial geometery in a big box. Right: the pressurized system, self-assembled in an NPT run

All simulations start from a random geometry in an enormously large box. In the first stage of the simulation, we perform a NVT simulation, that is coupled to a Nosé -Hoover thermostat, in a large box to thermostat the temperature of the system to a desired temperature (figure 15). In the next stage, we shrink the simulation box to a very high density via a NPT ensemble simulation while also self-assembly of the block copolymers takes place. Once the temperature and the total energy of the system become stationary, we perform a series of NPT runs to expand the box until the pressure is just zero. This is the volume at which the chains precisely fit inside the box, but do not exert any outward (or inward) pressures on it. We take this volume to be the proper volume of a particular self-assembled configuration, and use it to determine the density. This density, , corresponds - though likely not identical - to the overlap concentration in experimental systems which is:


where is degree of polymerization and is radius of gyration of polymer (Daoud et al., 1975). Since the number of beads in our simulations do not vary, we can say that . This run is followed by more NPT runs which the simulation box is compressed to provide higher density systems.

The rheology simulations are performed with periodic boundary conditions using the NVT ensemble, which serves to ensure that the pressure of the system changes corresponding to the applied deformation - in LAMMPS the stress must be determined from the pressure tensor. We solve the SLLOD equations of motion which were proven to be equivalent to Newton’s equations of motion for shear (Evans and Morriss, 1984). In this method, instead of boundary driven deformation , where shearing deformation is induced to the particles by the motion of the boundaries, a velocity gradient is generated to move all the particles proportional to the deformation of the boundaries. These equations are coupled to a Nosé -Hoover chain thermostat in a velocity Verlet formulation. The oscillatory shear strain - according to Eq. (6) - is applied to the simulation box in the -plane, for various amplitudes and oscillation periods. To obtain the shear stress response, we compute the component of the symmetric (Cauchy) stress tensor, , for every bead in the simulation box and sum over all atoms every simulation time-steps. In the following sections we investigate the feasibility of performing rheology with the LAMMPS molecular dynamics package (Plimpton, 1995).

vii.2 Computing stress and temperature control

The oscillatory shear simulation imparts a continuously changing deformation to the simulation box. As a result, for affine deformation, each atom (bead) in the simulation box can be thought of as being forced to drift at a given velocity. For example, if the box is being sheared in the -plane the atoms at the bottom of the box (low ) have a smaller velocity in the -direction than those atoms at top of the box (at high ). LAMMPS subtracts this spurious position-dependent drifting velocity from each atom while shearing.

To obtain the viscoelastic modulus of a molecular system using Eq. (8), one needs to compute the stress response of each particle inside the system to external deformation. The Cauchy stress tensor, a by tensor, completely defines the state of stress at any particular point of a material structure in a given deformed state. Written out in components, the Cauchy stress tensor has the following form:

Figure 16: Stress fluctuations. Here we show the stress response over one period of the strain, for two different maximal shear values. The stress response is drawn in red every cycles, and the blue line is the average of all the cycles. (a) The shear rate is slow compared to fluctuation timescales, so the stress response is slightly different in every cycle. (b) For higher deformation, larger than , the system becomes more rigid and fluctuations are suppressed.

Thus, is the force in the -direction on the surface whose normal is in the -direction. The diagonal components of the stress tensor are force per area in all three dimensions, and contain the hydrostatic pressure as well as any normal stresses that develop inside the material. In our simulations, the complete stress tensor for atom , multiplied by the volume that the atom occupies, in our simulations is computed as follows:

Figure 17: The stress signal for , and , reconstructed from only the first harmonic of the full Fourier expansion fits very well to the numerical stress data, evidencing that we am in a regime of linear response.

The first term is a kinetic energy contribution for atom . The second term is a pairwise energy contribution where runs over the neighbors of atom which are all the atoms that are within the cut-off range of the potential, and are the position of the two atoms involved in the pairwise interaction, and and are the forces on the atoms resulting from the pairwise interaction. The third term is a bond contribution of similar form for the bonds which atom is part of. There is also a similar term for the angle interactions that atom is involved in. The so-called virial stress tensor is similar to the Cauchy stress tensor, and has all the above terms except for the contribution from the kinetic energy. In our systems we find no significant difference in the moduli calculated by these two different methods of computing the stress tensors, because the actual velocities remain low and the energy is dominated by the bonded and non-bonded contributions. The comparison is presented in Fig. [20]. As discussed in section VII.1, we convert the stress signal to dynamic moduli using discrete Fourier function. In the linear regime only the first Fourier harmonic of the stress signal has contribution to the modulus. To check that the calculated modulus is in the linear regime we reconstruct the stress signal from only the first Fourier harmonic using the following equation


where is the component of Eq. (11). The reconstructed stress, obtained from the above formula, fits very well to the stress data from simulations. This is shown in the Fig. [17].

Figure 18: Strain dependent modulus for B3L45, , . The rise of modulus in the high shear rate region is due to bond stretching and hard core potential of the beads during the simulation. This regime is unphysical - at such high strain rates the implicit solvent interactions break down and our simulations can no longer capture the actual deformations of the material. We discard points after the minumum, and interpret the rapid dropoff as yielding behavior for our system.

In molecular dynamics simulation of oscillatory shear, care should be taken in choosing the right deformation rate. If the box deformation rate is larger than the time-scale in which the beads interact and fluctuate then the position and interaction of the beads are the same in all of the oscillatory cycles (see Fig. [16]) and the response is similar to the response of a crystalline materials which has significantly higher modulus than soft materials. Obviously increasing the shear rate eventually would lead in stretching the bonds more than equilibrium bond distance which breaks the bonds and stops the simulation. In the Fig. [18] is shown that the mechanical response of the system increases for which is towards suppressing the fluctuations and solid-like structures. This rise in modulus eventually terminates at where the simulation stops because of excessive bond stretching.

In MD simulations a thermostat must be used as a means of controlling the temperature of particles. Typically a target temperature is specified by the user, and the thermostat attempts to equilibrate the system to the requested . To compute the temperature, the kinetic energy is divided by the Boltzmann constant and , number of degrees of freedom present in the system:


Since the kinetic energy is a function of particle velocity, there is often a need to distinguish between a particle’s streaming velocity which occurs due to group motion of aggregated particles and its thermal velocity due to thermal fluctuation. The sum of the two is the particle’s total velocity. Using the Nosé -Hoover equations (Nosé, 1984) we thermostat the translational velocity of particles and subtract a velocity bias that is the result of deforming the simulation box. Hence the dependence of the computed stress on temperature is only the contribution of thermal velocity which is due to the fluctuations of the particles.

Figure 19: Simulation snapshots of (a) single long chain with 500 repeating binder+linker blocks and (b) split chain, where the separate chains are distinguished by different colors.(c) Strain-dependent simulation results for two different systems of single chain (solid line) and split chain (dashed line). Storage modulus is shown in blue color and loss modulus is shown in red color. As may be seen here, the results are identical for the split and the long system.

Viii Results

Now, we present the results of a series of oscillatory rheology simulations of the self-assembling multiblock copolymer system. First, we study the effect of changing the molecular architecture from a single long chain to several short chains: Strain-controlled simulations results at different densities are shown in Fig. [19]. To justify the validity of our coarse-grained model to simulate a block copolymer, we compare the results of two different system, a one repeating chain that is self-assembled to a flowerlike network split chains which are chains of hydrophobic-hydrophilic block each. The system sizes in both cases are equal to particles. Fig. [19] shows equivalent results for both systems.

Second, we demonstrate the practical equivalence (in our simulation settings) for the virial and Cauchy stress approaches (Zhou, 2003; Zimmerman et al., 2010). Here we compare modulus of one system using Cauchy and Virial stresses. As shown in Fig. [20], the difference of two stresses in our model is small.

Figure 20: Strain-dependent modulus of a of BL crosslinked network at , . The complex modulus was calculated from the two different stress tensors: Cauchy and Virial.
(a) Harmonic spring bond,
(b) FENE bond,
Figure 21: Strain-dependent comparison between a bead-spring polymer network of harmonic bonds vs a bead spring network of FENE bonds for BL, and . The storage modulus is shown in solid symbols, the loss modulus is shown in open symbols.

Finally, we examine the effects of the details of the bond potential. As discussed in section VIII, FENE bonds are often used as finitely extensible bonds in coarse grained MD simulations to better reflect the finite backbone length of polymers. To compare the effect of harmonic bonds vs. FENE bonds in the shear deformation simulations we compare the modulus for two systems containing two above bonds in Fig. [21]. Again, the results are completely similar: This may be understood from the fact that in the regime where we measure, the non-linear stretching regime is never engaged and the FENE bonds act as linear springs.

The effect of concentration is also shown in Fig. [22]. The overlap concentration - which is the concentration at which the single molecule network precisely fits within the box - is obtained by relaxing the network volume until the pressure first reaches zero. Higher concentrations are shown in units of this overlap concentration.

(a) ,
Figure 22: Dynamic moduli and , measured in strain-controlled settings at and at various concentrations. The storage modulus is shown in solid symbols, the loss modulus is shown in open symbols.
Figure 23: Dynamic moduli and , measured in strain-controlled settings at and at various temperatures. is shown by solid symbols and is indicated by open symbols.
Figure 24: The linear modulus as a function of temperature and density for BL at . (a) The slope of power-law fit is .

The temperature of the simulation also influences the modulus of the system. Using Eq. (VII.2), we can change the temperature in the simulation. It affects the mechanics in two ways: first, the self assembly process is affected as elevated temperatures favor high entropy states more than high energy states, which skews the architectures towards smaller cluster sizes. This tends to lower . Secondly, the mechanical response due to the linkers is itself temperature dependent: entropically elastic effects scale with temperature. This latter effect would raise the storage modulus as a function of temperature, as is the case in rubber elasticity. In the Fig. [23] the strain dependent shear modulus for at different temperatures is shown. The modulus is calculated based on the stress data from simulations that include all the contributions of the stress tensor. The result shows a linear regime up to , after which drops off rapidly. Linear modulus as a function of temperature and density is shown in a logarithmic plot in Fig. [24]. The slope of the power law fit to the temperature dependant linear modulus from numerical results is , that is close to theoretical prediction for spherical flower-like micelles, i.e. 1.95.

The yield point for all of the systems in Fig. [23] is where the starts to drop, at . That this is the yield point is also evidenced by the concurrent and characteristic rise in . To investigate the effect of the temperature on the modulus we plot the elastic modulus for all temperatures as a function of strain in 25. The experimental measurements at large strains gives a quantitative comparison of the contribution of the elastic response at different temperatures(Fahimi, 2014). A decrease in the magnitude of the modulus with the temperature is observed which is in good agreement with the simulation results we find here Fig. [25]. The fact that the overall effect of temperature in these networks is to lower the modulus suggests that, apparently, the effect of temperature on the aggregate size is dominant over the single-chain elastic contribution. This suggests that rubber elasticity theories should be modified to include also the temperature dependent changes in connectivity to fully capture the mechanical response of the system (and that, therefore, the model presented in section II may not offer a complete description.)

Figure 25: Dynamic moduli , measured in strain-controlled settings at and at various temperatures. Different colors corresponds to different temperatures as indicated in the legend.

Ix Conclusion

Our MD simulations suggest that the motif of using repeating hydrophilic and hydrophobic blocks gives rise to interesting self-assembling phase behavior. In a large regime of phase space (defined by a combination of and ) we observe, in NVT MD simulations, network phases where small clusters of the hydrophobic blocks are connected by (typically) single hydrophilic linker chains to other such clusters, giving rise to a ”single molecule networks”. The connectivity of this network is generally determined by a combination of the size of an aggregate, which sets the number of outgoing linkers and thus the potential for bridges, and the number of aggregates which sets the number of potential partners for such bridging. At intermediate values of and , the combination between ”supply and demand” of linkers is optimal which should result, in the highest moduli for the network material. The multiblock amphiphile system thus allows direct control over its supramolecular arrangement through molecular design, and thus to mechanical properties. This suggests much richer applications for these systems than what has been established thus far, and in particular makes them suitable candidates for exploring future biomimetic mechanical performance.

We thank Hans Wyss, Rint Sijbesma, Gajanan Pawar, Zahra Fahimi, Marcel Koenigs, Christian Moerland, Thijs Michels, Alexey Lyulin and Arlette Baljon for enlightening discussions and for sharing experimental details and results. Funding from the Eindhoven University of Technology’s High Potential programme is gratefully acknowledged.


  • Bausch and Kroy (2006) A. Bausch and K. Kroy, Nat. Phys. 2, 231 (2006).
  • Broedersz and MacKintosh (2014) C. Broedersz and F. MacKintosh, Rev. Mod. Phys. 86, 995 (2014).
  • Sarkar (1979) N. Sarkar, J. Appl. Polym. Sci. 24, 1073 (1979).
  • Pawar et al. (2012) G. M. Pawar, M. Koenigs, Z. Fahimi, M. Cox, I. K. Voets, H. M. Wyss, and R. P. Sijbesma, Biomacromolecules 13, 3966 (2012).
  • Koenigs et al. (2014) M. Koenigs, A. Pal, H. Mortazavi, G. Pawar, C. Storm, and R. Sijbesma, Macromolecules 47, 2712 (2014).
  • Alexandridis and Lindman (2000) P. Alexandridis and B. Lindman, Amphiphilic block copolymers: Self assembly and applications (Elsevier, Amsterdam, 2000).
  • Halperin (1992) A. Halperin, Macromolecules 24, 1418 (1992).
  • Zhang et al. (2003) G. Zhang, F. Winnik, and C. Wu, Phys. Rev. Lett. 90, 035506 (2003).
  • Kawata et al. (2007) T. Kawata, A. Hashidzume, and T. Sato, Macromolecules 40, 1147 (2007).
  • Tanaka and Koga (2000) F. Tanaka and T. Koga, Comp. Theor. Polym. Sci. 10, 259 (2000).
  • Tanaka (2002) F. Tanaka, J. Non-Cryst. Solids 307-310, 688 (2002).
  • Hugouvieux et al. (2009) V. Hugouvieux, M. A. V. Axelos, and M. Kolb, Macromolecules 42, 392 (2009).
  • Hugouvieux et al. (2011) V. Hugouvieux, M. A. V. Axelos, and M. Kolb, Soft Matter 7, 2580 (2011).
  • Cooke and Williams (2003) I. R. Cooke and D. R. M. Williams, Macromolecules 36, 2149 (2003).
  • Zhang et al. (2013) J. Zhang, Z.-Y. Lu, and Z.-Y. Sun, Soft Matter 9, 1947 (2013).
  • Plimpton (1995) S. Plimpton, Journal of computational physics 117, 1 (1995).
  • Wilhelm and Frey (2003) J. Wilhelm and E. Frey, Phys. Rev. Lett. 91, 108103 (2003).
  • Gardel et al. (2004) M. L. Gardel, M. T. Valentine, J. C. Crocker, A. R. Bausch, and D. A. Weitz, Science 304, 1301 (2004).
  • Gardel et al. (2003) M. L. Gardel, M. T. Valentine, J. C. Crocker, A. R. Bausch, and D. A. Weitz, Phys. Rev. Lett. 91, 158302 (2003).
  • Storm et al. (2005) C. Storm, J. J. Pastore, F. C. MacKintosh, T. C. Lubensky, and P. A. Janmey, Nature 435, 191 (2005).
  • Huisman et al. (2007) E. M. Huisman, T. van Dillen, P. R. Onck, and E. van der Giessen, Phys. Rev. Lett. 99, 208103 (2007).
  • Yeo et al. (1981) J. Yeo, L. Sperling, and D. Thomas, J. Appl. Pol. Sci. 26, 3977 (1981).
  • Rubinstein and Colby (2007) M. Rubinstein and R. H. Colby, Polymer Physics (Oxford University Press, 2007).
  • Duiser and Staverman (1965) J. Duiser and A. Staverman, in Physics of non-crystalline solids, edited by J. Prins (North-Holland, Amsterdam, 1965).
  • Graessley (1975) W. W. Graessley, Macromolecules 8, 186 (1975).
  • Tobolsky et al. (1961) A. Tobolsky, D. Carlson, and N. Indicator, J. Pol. Sci. 54, 175 (1961).
  • Flory (1953) P. J. Flory, Principles of polymer chemistry (Cornell University Press, 1953).
  • Stockmayer (1944) W. Stockmayer, J. Chem. Phys 12, 125 (1944).
  • Wilson et al. (2011) M. Wilson, A. Rabinovitch, and A. R. C. Baljon, Phys. Rev. E 84, 061801 (2011).
  • Skrzeszewska et al. (2010a) P. J. Skrzeszewska, F. A. de Wolf, M. A. Cohen Stuart, and J. van der Gucht, Soft Matter 6, 416 (2010a).
  • Skrzeszewska et al. (2010b) P. J. Skrzeszewska, J. Sprakel, F. A. de Wolf, R. Fokkink, M. A. Cohen Stuart, and J. van der Gucht, Macromolecules 43, 3542 (2010b).
  • Bird et al. (1977) R. B. Bird, R. C. Armstrong, O. Hassager, and C. F. Curtiss, Dynamics of polymeric liquids, vol. 1 (Wiley New York, 1977).
  • Baschnagel and Varnik (2005) J. Baschnagel and F. Varnik, Journal of Physics: Condensed Matter 17, R851 (2005).
  • Daoud et al. (1975) M. Daoud, J. Cotton, B. Farnoux, G. Jannink, G. Sarma, H. Benoit, C. Duplessix, C. Picot, and P. De Gennes, Macromolecules 8, 804 (1975).
  • Evans and Morriss (1984) D. J. Evans and G. Morriss, Physical Review A 30, 1528 (1984).
  • Nosé (1984) S. Nosé, The Journal of Chemical Physics 81, 511 (1984).
  • Zhou (2003) M. Zhou, Proceedings of the Royal Society of London. Series A: Mathematical, Physical and Engineering Sciences 459, 2347 (2003).
  • Zimmerman et al. (2010) J. A. Zimmerman, R. E. Jones, and J. A. Templeton, Journal of Computational Physics 229, 2364 (2010).
  • Fahimi (2014) Z. Fahimi, Ph.D. thesis, Eindhoven University of Technology (2014), an optional note.
Comments 0
Request Comment
You are adding the first comment!
How to quickly get a good reply:
  • Give credit where it’s due by listing out the positive aspects of a paper before getting into which changes should be made.
  • Be specific in your critique, and provide supporting evidence with appropriate references to substantiate general statements.
  • Your comment should inspire ideas to flow and help the author improves the paper.

The better we are at sharing our knowledge with each other, the faster we move forward.
The feedback must be of minimum 40 characters and the title a minimum of 5 characters
Add comment
Loading ...
This is a comment super asjknd jkasnjk adsnkj
The feedback must be of minumum 40 characters
The feedback must be of minumum 40 characters

You are asking your first question!
How to quickly get a good answer:
  • Keep your question short and to the point
  • Check for grammar or spelling errors.
  • Phrase it like a question
Test description