Nanotribology of biopolymer brushes in aqueous solution using dissipative particle dynamics simulations: an application to PEG covered liposomes in a theta solvent

Abstract

We undertake the investigation of sheared polymer chains grafted onto flat surfaces to model liposomes covered with polyethylene glycol brushes as a case study for the mechanisms of efficient drug delivery in biologically relevant situations, for example, as carriers for topical treatments of illnesses in the human vasculature. For these applications, specific rheological properties are required, such as low viscosity at high shear rates, to improve the transport of the liposomes. Therefore, extensive non-equilibrium, coarse-grained dissipative particle dynamics simulations of polymer brushes of various lengths and shear rates are performed to obtain the average viscosity and the friction coefficient of the system as functions of the shear rate and polymerization degree under theta-solvent conditions, and we find that the brushes experience considerable shear thinning at large shear rates. The viscosity (η) and the friction coefficient (μ) are shown to obey the scaling laws η ∼[small gamma, Greek, dot above]^−0.31 and μ ∼ [small gamma, Greek, dot above]^0.69 at high shear rates ([small gamma, Greek, dot above]) in a theta solvent, irrespective of the degree of polymerization of brushes. These results confirm recent scaling predictions and reproduce very well trends in measurements of the viscosity at a high shear rate ([small gamma, Greek, dot above]) of red blood cells in a liposome containing medium.

---

I. Introduction

Assemblies of long-chain macromolecules which are grafted onto surfaces at sufficiently high densities are called polymer brushes. Some of their most important current applications are found in the stability of colloidal dispersions,^1,2 in the process of enhanced oil recovery,³ in the design and applications of stimuli-responsive materials,⁴ in the characterization of the mechanical properties of cancerous human cells,⁵ and in the emerging field of nano-patterning, for example.⁶ They display some general fundamental properties, and scaling relations have been derived for them under circumstances such as high grafting density or negligible chain–chain interaction.^7–9 When polymer brushes are sheared a plethora of phenomena that are intrinsic to the non-equilibrium nature of the shearing appear. For example, experiments on fluids compressed between plates under the influence of external shear have shown that the viscosity of the fluid is reduced when polymer chains are grafted onto each plate.¹⁰ Biological systems that share similarities with polymer brushes occur in articular cartilage surfaces¹¹ and synovial joints,¹² in glycocalyx filaments that coat the human circulatory system,¹³ and in glycosylated cell surfaces and liposomes, which can be used as carriers for drug delivery.^14–16 Despite these studies, the experimental understanding of the molecular mechanisms in various environments is still incomplete, because, among other reasons, parameters such as the thickness of the polymer brushes on the sheared surfaces are still difficult to determine.⁶ These difficulties arise partly from the self-assembly processes which are typically used for the modification of surfaces and depend on chemical reactions that do not necessarily yield a well-defined brush length.⁶ In this regard, computer simulations have come into play with an increasingly significant role, and there are now various studies investigating the role of the polymer brush length, the spacing between the plates, the properties of the solvent and the influence of the electrostatic interactions, the fluid's shear thinning as the shear rate is increased, and the effect of changing the polymer grafting density on the plates, to name but a few.^17–25

Most of the studies have been carried out for polymers under good solvent conditions, with very few exceptions.^20,22–24 However, it is known that some biopolymers found in an aqueous environment are in the borderline between good solvent and theta solvent conditions.^22,26 An important example is that of polyethylene glycol (PEG) polymer brushes on liposomes, which are used as carriers under flow for efficient drug delivery.^27,28 PEG is known to have a decreasing solubility as its molecular weight increases,²⁹ and for molecular weights in the range of 2000 to 10 000 g mol⁻¹ in aqueous solution it has been reported to be slightly below its theta temperature, when placed at the human body temperature.³⁰ A detailed understanding of drug delivery mechanisms where PEG is involved, which is of paramount importance for optimal pharmaceutical design, is still far from complete.³¹ Therefore, our purpose here is to study the viscosity and friction of PEG brushes on surfaces that experience an external shear force under theta solvent conditions, for the first time. To reach scales comparable to those of typical non-equilibrium experiments we have carried out particle-based, mesoscopic molecular dynamics simulations, using the method known as dissipative particle dynamics (DPD),³² which has successfully predicted properties of polymer brushes both in equilibrium and non-equilibrium situations.³³

This work reports explicit-solvent DPD simulations of the rheological properties of biopolymer (PEG) brushes grafted onto parallel flat surfaces (liposomes). Including the solvent explicitly in the simulations is crucial to reproduce experiments measuring friction between polymer brushes.^18,19 We want to model the behavior of tethered proteins and biopolymers in environments such as vessels, where the separation between plates is essentially constant but we allow for variations in the degree of polymerization, at a fixed grafting density and substrate separation. In particular, we would like to determine the optimal conditions under which PEG brushes on liposomes promote the flow of the latter as mechanisms for efficient drug delivery. We obtain the viscosity and the friction coefficient as functions of the shear rate and degree of polymerization. This is the first report of scaling behavior for the viscosity and the friction coefficient at high shear rates of polymer brushes in a theta solvent using soft DPD interactions. Understanding of these phenomena is useful for the interpretation of several recent experiments in biological and other mesoscopic systems of current academic and industrial interest.³⁴

In Section II, called “Models and methods”, we provide details about the DPD interaction model, the simulation algorithm and the systems studied. In Section III the main results are presented, accompanied by their discussion. The final section is devoted to our major conclusions.

II. Models and methods

We have performed DPD simulations of linear grafted polymers in the canonical ensemble (fixed number of particles, volume and temperature). The DPD model is by now well known; therefore we shall be brief here. Three forces make up the basic DPD structure: a conservative force (F⃑^C_ij) that accounts for the local pressure and is proportional to an interaction constant, a_ij; a dissipative force (F⃑^D_ij), which represents the viscosity arising from collisions between particles, proportional to the relative velocity of the particles and to a constant, γ; and a random force (F⃑^R_ij) that models the Brownian motion of the particles with an intensity given by the constant σ, all acting between any two particles, i and j. The spatial dependence of these forces is usually chosen as shown in eqn (1)–(3).

F⃑^C_ij = a_ij(1 − r_ij/R_c)ê_ij (1)

F⃑^D_ij = −γ(1 − r_ij/R_c)²[ê_ij × v⃑_ij]ê_ij (2)

F⃑^R_ij = σ(1 − r_ij/R_c)ê_ijξ_ij (3)

In eqn (1)–(3), r_ij = r_i − r_j represents the relative position vector between particles i and j, ê_ij is the unit vector in the direction of r_ij, while v_ij = v_i − v_j is the relative velocity between particles i and j, with r_i and v_i being the position and velocity of particle i, respectively. The random variable ξ_ij is generated between 0 and 1 with a Gaussian distribution, zero mean and unit variance; R_c is the cut off distance, beyond which all forces are zero. The DPD beads are all of the same size, with radius R_c, which is set equal to 1. The constants in the dissipation and random forces are not independent and satisfy the relation:³⁵

which, in effect, fixes the temperature; here, k_B is Boltzmann's constant. This natural thermostat that arises from the balance between the dissipative and random forces is a defining feature of the DPD model.³³ Also, the short-range nature of the DPD forces and their linearly decaying spatial dependence allow the use of relatively large time steps when integrating the equation of motion. The DPD particles are representations of sections of fluid rather than physical particles, and can group several atoms or molecules, making the DPD method an attractive alternative to study systems at the mesoscopic level. Note also that the forces in eqn (1)–(3) obey Newton's third law, which means momentum is conserved locally, and globally, which in turn preserves any hydrodynamic modes present in the fluid. This is a feature of fundamental importance when studying non-equilibrium properties of fluids, since loss of information about these modes can lead to a different phase from that obtained when they are fully accounted for.³⁶ The DPD interaction model has been used successfully to predict equilibrium properties of polymer melts,³⁷ surfactants in solution,³⁸ and colloidal stability,³⁹ among others. For further reading, see ref. 33.

All our simulations are performed in reduced units (marked with asterisks); distances are reduced with the cutoff radius, R_c, which for a coarse-graining degree equal to three water molecules per DPD particle is equal to R_c = 6.46 Å,³⁸ hence r = r*R_c. The time step δt is reduced with δt = (mR²_c/k_BT)^1/2δt*, where m is the mass of a DPD particle, while energy is reduced with k_BT. All other quantities can be reduced through combinations of these relations, for example, the viscosity: η = η*(k_BTδt/R³_c). Using the mass of 3 water molecules per DPD particle, at room temperature, one obtains δt ≈ (6.3 × 10⁻¹² s)δt* and η ≈ (9.64 cP)η*. All our simulations are performed for (k_BT)* = m* = R^*_c = 1, δt* = 0.03. In addition to the fluid monomeric particles, our system contains soft parallel surfaces as models for biological membranes and polymer chains attached to them, forming brushes. The polymers are built as linear chains of DPD beads joined by freely rotating harmonic springs.⁴⁰ The spring constant is chosen as κ = κ*(k_BT/R²_c), where κ* = 100 and the spring's equilibrium position is r_eq = r^*_eqR_c, with r^*_eq = 0.7 in all cases.⁴¹ In our simulations we have fixed the distance between the plates, which are placed perpendicular to the z-direction, at a distance of D* = 7. Periodic boundary conditions are applied on the xy-plane but not in the z-direction, to reinforce the confinement. The equation of motion is solved using the velocity Verlet algorithm, adapted to DPD.⁴² The parameters of the dissipative and random force intensities are, respectively, γ = 4.5 and σ = 3.0, so that k_BT* = 1 (see eqn (4)). The value of the conservative force intensity (see eqn (1)) was set to a_ij = 78.0 for all cases, namely for the interaction between particles of the same type (solvent–solvent, monomer–monomer) and for particles of different types (solvent–monomer). The choice a_ij = 78.0 is obtained when one uses a coarse graining degree that groups 3 water molecules in a single DPD bead.³⁸ Since all particle–particle interactions are equal, the polymers are in a theta solvent, as is the case for PEG-covered liposomes in the human circulatory system.³⁰ It has been shown that this model for PEG on colloidal surfaces can successfully predict brush scaling laws,²⁸ as well as adsorption isotherms on metallic oxides and disjoining pressure profiles.³⁹ The soft membranes on which these biopolymers are tethered are modeled as linearly decaying forces that act on the ends of the simulation box (in the z-direction), given by the following expression:³⁹

In eqn (5), F_w(z_i) symbolizes the force exerted by the effective wall on particle i, whose position component in the z-direction is z_i. The intensity of the force is given by the constant a_w, while z_C is the cutoff distance also set equal to 1, which defines the reach of this force, i.e., F_w(z_i) = 0 for z_i > z_C. Eqn (5) represents a soft membrane because the maximum repulsion between the beads and the surface is a_w when z_i = 0; harder walls can be used for other DPD applications.³⁹ It is known that liposomes can be up to about 3 orders of magnitude larger than their brush thickness, and the solvent,^2,43 therefore we consider planar walls only. To fix one end of the polymers on these substrates we chose their interaction equal to a_w = 70.0 (reduced energy over length DPD units), while the rest of the polymer beads and the solvent interact with the walls with an intensity given by a_w = 140.0. With this choice of parameters, the polymer end experiences a less repulsive force toward the walls than the rest of the polymer (and solvent) beads, and is therefore adsorbed on them. The tethered ends of the polymers remain free to move on the xy-plane, of course, as occurs also for biopolymers interacting with membranes. Although this model for membranes is soft, it remains impenetrable to polymer and solvent molecules, and we assume moreover that the drug molecules have already been incorporated into the structure.

To carry out non-equilibrium simulations of these systems, we apply a constant shear velocity to the polymer ends adsorbed on the walls, of equal magnitude but opposite sign for beads on different plates. This is equivalent to moving the plates in opposite directions under the influence of a fixed external force that keeps the plates moving with a constant speed, known as the Couette flow,^20–24 see Fig. 1. The units of velocity are given by [v₀] = v^*₀[R_c/δt] ≈ 1.01 × 10²v^*₀ m s⁻¹, using the values of R_c and δt given previously. Some of the solvent monomers are seen to penetrate the PEG brushes and carry them along in their flow to produce lubrication, while there is no interpenetration of the brushes. As pointed out before, our simulations are performed in the canonical ensemble, where particle density and temperature are kept constant; previous studies have shown that completely equivalent results are obtained if one performs the simulations in the grand canonical ensemble, at fixed chemical potential, volume and temperature.¹⁸ To keep the simulation results invariant to the particular choice of a_ij interaction parameters, we fix the global dimensionless density to 3 (ref. 38) in all the simulations reported here; this means that when the degree of polymerization of the PEG brushes (N) is increased, the number of solvent molecules is reduced.

Fig. 1 Snapshot of two linear PEG brushes made up of N = 7 beads. The lateral dimensions of the simulation box are L^*_x = 35, L^*_y = 35, and the spacing between the plates is L^*_z = D* = 7. The grafting density is equal to Γ* = 0.3. A constant shear velocity of magnitude v^*₀ = 1.0 is applied to the tethered beads of each polymer (in dark blue); the rest of the polymer beads are shown in yellow. The solvent monomers are shown in red. The dark blue beads seen above and below the blue layers are the tethered ends of chains that were adsorbed closer to the surfaces when shear was applied, due to bead–bead collisions. Notice how the solvent penetrates the polymer brushes all the way up to the surfaces for this grafting density. See text for details.

The grafting density is defined as Γ* = N_p/A*, where N_p is the number of polymer chains tethered on the surface and A* = L^*_xL^*_y is its area. Our simulations are performed at grafting densities that lie within the brush regime (Γ^*_θ ∼ 1/N, for theta-solvent conditions, where N is the degree of polymerization; for polymer chains in a good solvent, Γ^*_good ∼ 1/N^2(0.588), which means that the brush regime is quickly reached for large chains in a good solvent, see ref. 28), where the “stealth” properties of the PEG brushes work best.⁴⁴ These properties are the increased lifetime in the bloodstream and the kinetic stability of nanoparticles, such as liposomes, due to the protective layer formed by the brushes.³¹ We obtain the viscosity (η) and the friction coefficient (μ) of the fluid through the relations.²⁰

In eqn (6) and (7) above, 〈F_x([small gamma, Greek, dot above])〉 and 〈F_z([small gamma, Greek, dot above])〉 are the mean forces that the particles on the surfaces experience along the flow direction and perpendicularly to it, respectively; the brackets indicate time average over all particles. Eqn (6) can be understood as the local definition of viscosity, applied to the entire sample. Considering a liquid confined between two planar walls, a linear flow can be generated by moving for example the top wall at a constant velocity. This motion requires that a steady force be applied on the top wall and an equal in magnitude but opposite in direction force on the bottom wall, see Fig. 1. In linear flow, the velocity gradient (shear rate, [small gamma, Greek, dot above]) is constant throughout the liquid. The shear rate [small gamma, Greek, dot above] is defined as 2v^*₀/D*, where v^*₀ is the flow velocity exerted on the grafted monomers, and D* is the surface separation (see Fig. 1). To obtain first a local definition of viscosity, a small cubic volume in the fluid (fluid element) can be considered. This fluid element undergoes strain, being deformed from a cubic to a parallelepiped shape, with increasing deformation in time. The amount of strain is the added lateral displacement of top and bottom planes of the fluid volume divided by its height. This strain increases at a constant rate in time for a linear flow, such that the strain rate, i.e. the strain per unit time, is the velocity gradient in the limit of an infinitesimal fluid element. The strain rate of deformation is therefore the shear rate, [small gamma, Greek, dot above]. On the other hand, a force acts on the top and bottom planes of the fluid element, with identical magnitude and opposite direction. The center of mass of the fluid element suffers no acceleration and therefore the total external force on it must vanish. These forces must be proportional to the area of the top and bottom planes of the fluid element. Also, the force-per-unit-area or stress, σ, must be uniform from top to bottom in the fluid element since the fluid is not accelerating. For a linear flow, this is true not only for a fluid element but for the whole fluid: every part of the fluid is under the same stress, which produces everywhere the same velocity gradient ([small gamma, Greek, dot above]). The definition of viscosity assumes that the shear rate ([small gamma, Greek, dot above]) produced in a fluid element is proportional to the shear stress (σ) exerted on it. The viscosity is defined as the constant of proportionality between them: σ = η[small gamma, Greek, dot above], or η = σ/[small gamma, Greek, dot above]. The total shear stress on the sample (σ) can be determined from the average total force on the brush heads divided by the wall area (F_x/A) while the shear rate ([small gamma, Greek, dot above]) is extracted from the slope of the linear fit of the average velocity profile. The viscosity is then expressed by eqn (6). As for the friction coefficient, shown in eqn (7), its calculation follows directly from the analogy with the friction coefficient between solid surfaces, namely F_x = μF_z, where F_z is the force acting perpendicularly to the surface, and F_x is the force acting parallel to it, which is responsible for the shear stress (σ = F_x/A).²⁵ Our results were obtained from averages of simulations of up to 4 × 10³ blocks, with each block composed of 2 × 10⁴ time steps, using the first 2 × 10⁵ time steps for equilibration and the rest for the production phase; when properly dimensionalized this represents a time observation window of 0.12 ms. Typical density profiles of the solvent and polymer monomers are shown in Fig. 2 for two polymerization degrees, where it is clear that the solvent penetrates the PEG brushes and reaches the surfaces, at that given grafting density.

Fig. 2 Density profiles of the PEG brushes (dashed lines, in blue for N = 7 and in green for N = 14), and the solvent (full lines, in black for N = 7 and in red for N = 14). The polymer monomers order near the walls, which leads to the layering effect, represented here by peaks near the walls, while there is a relatively free (bulk like) fluid made up of solvent molecules at the center of the simulation box. For the case shown in this figure, L^*_x = L^*_y = 7, D* = 7, v^*₀ = 1.0 and Γ* = 1.0. Note there is very little interpenetration between the brushes. The brush layering is rather strong for the N = 7 brush because the chains are relatively short, while for N = 14 the brush layering is considerably reduced.

The profiles of the opposing polymer brushes show very little interpenetration and an almost bulk like concentration of solvent monomers at the center of the simulation box (see Fig. 2), which reduces the viscosity of the system, because there are more momentum – carrying solvent molecules in the gap between the surfaces than if there were longer brushes. The total density (solvent plus polymer monomers) of the system was fixed at three; hence at the center of the simulation box one finds mostly solvent particles. The structuring of the PEG brush monomers, represented by the maxima in Fig. 2, coincides with that of the solvent monomers, indicating that even for a relatively large grafting density as that shown in the figure, the solvent is able to reach the membrane. The maxima of both types of profiles (the solvent's and the brushes') occur at the same positions because the interactions between them are the same (recall a_ij = a_ii = 78.0 for all i and j), as it is befitting for theta-solvent conditions. Theta solvent conditions are obtained whenever the interactions between polymer chains are the same as the interactions between the solvent and the polymers.² A strong layering is expected when the chains are relatively short, as for N = 7 in Fig. 2, because in such a case monomers are more easily arranged than when long chains form the brushes.^17–21 In Fig. 3 we present the velocity profile for PEG brushes of polymerization degree N = 14 and grafting density Γ* = 0.30 under a shear velocity equal to v^*₀ = 0.1. Clearly, there appears a linear gradient in the center of the channel formed between liposomes, from which the shear rate [small gamma, Greek, dot above] can be obtained, through the relation [small gamma, Greek, dot above] = 2v^*₀/D*, as pointed out before. The inset of Fig. 3 shows the temperature profile of the complex fluid in the z-direction, which indicates the brushes are at the temperature fixed by the thermostat (k_BT* = 1), although at large shear rates one expects the brush to “heat up” somewhat,²⁰ as a consequence of the increased dissipation rate.

Fig. 3 Velocity profile at the center of the channel between liposomes for PEG brushes with N = 14 and Γ* = 0.30, using v^*₀ = 0.1. The box size is L^*_x = L^*_y = 7, D* = 7. The inset shows the temperature profile (T^*_x) in the z-direction. All quantities are reported in reduced units.

Most simulations were performed for a brush grafting density Γ* = 0.30, except where indicated otherwise, because for this value the average distance between grafted heads on the surface is smaller than their radius of gyration in a theta solvent, for all values of N we studied.²⁸ By doing so one makes sure that the grafted polymer chains are in the brush, rather than the mushroom regime,⁴⁵ which is important for the situation we are interested in modeling. The polymerization degree was varied in the range N = 1 up to N = 25, so that scaling laws could be extracted from the data. Finally, the shear rate was chosen to vary from [small gamma, Greek, dot above] = 0 to [small gamma, Greek, dot above] = 0.30, except where indicated otherwise, so that we could compare our results with those available in the literature.^17–21 Molecular dynamics simulations like ours solve the equation of motion essentially exactly,³⁶ although the use of soft forces in DPD (see, for example, eqn (1)) might limit its applicability to situations where atomistic details are not important, as is the case for the nanotribology studies that are the focus of this work.

III. Results and discussion

One would like to simulate systems with a large number of particles so that macroscopic properties can be reproduced, while simultaneously keeping such numbers manageable to obtain results in a timely fashion. Therefore, we have carried out simulations to quantify the extent of finite size effects in the viscosity of the polymer brushes, if any. Fig. 4 shows the mean viscosity of the PEG brushes and the solvent, calculated using eqn (6), as a function of the lateral size of the membranes (L^*_x = L^*_y), at fixed separation between them (D* = 7), grafting density (Γ* = 0.30) and shear rate ([small gamma, Greek, dot above] = 0.28). The results in Fig. 4 cover a range that goes from 10³ particles up to 10⁵ particles while the relative change in the viscosity amounts to less than 1%, showing that finite size effects are not important in the calculation of the viscosity, which allows us to make correct predictions using relatively small systems.

Fig. 4 Finite size effects in the viscosity, η*. The symbol L^*_x represents the size of the simulation box in the x-direction, which is equal to that in the y-direction, L^*_y. In all cases, L^*_z = D* = 7, the shear rate is [small gamma, Greek, dot above] = 0.28, and the grafting density is Γ* = 0.30. All quantities are reported in reduced units. Lines are only guides for the eye.

Similar results have been obtained in equilibrium DPD simulations⁴⁶ which proved that finite size effects in the interfacial tension between two model fluids were as small as those found in the viscosity and shown in Fig. 4, when the appropriate ensemble was used. This behavior is to be attributed to the soft, short-range repulsive nature of the DPD forces (see eqn (1)–(3)), although these effects may become non-negligible if a long-range interaction is included, such as the electrostatic one. In the remainder of this work we present results obtained using only the smallest simulation box in Fig. 4.

Let us now proceed to discuss the results obtained for the mean viscosity of a complex fluid made up of PEG brushes and an aqueous solvent, as a function of the applied shear rate on the surfaces, for different degrees of polymerization (N), as seen in Fig. 5. The simulations were performed with a cubic box of lateral size L* = 7 and for Γ* = 0.30. For very small values of [small gamma, Greek, dot above] (less than 10⁻⁴) the calculation of η* requires very long simulations due to poor signal-to-noise ratios,^20,24 which makes it very difficult to obtain accurately the zero-shear viscosity, η₀, from simulations. However, a good estimate of it can be obtained from the extrapolation of the rheology profiles in Fig. 5 as [small gamma, Greek, dot above] approaches zero. For the shorter PEG brushes (N = 1 and N = 3) the fluid displays essentially Newtonian behavior, whereas for the case with N = 7 the fluid shows shear-thinning. One observes also that as N is increased so is the extrapolated η₀ value, which is a consequence of the combined effects of an increasing brush thickness and a reduced number of solvent particles (to keep the global density ρ* = 3). Increasing N is equivalent to reducing the separation between the membranes, D*, when the grafting density is kept constant, and yields increasing values of η₀.¹⁹ Hence, we have opted to vary N rather than D* in this work. The simulations shown in Fig. 5 were carried out for relatively short polymers with the purpose of establishing the minimum polymerization degree at which the fluid started to show shear thinning, which, from Fig. 5, is N = 7. For polymerization degrees smaller than 7, it is difficult for the polymer chains to interact with the solvent monomers enough to reduce the viscosity as the shear rate is increased, hence the fluid remains approximately Newtonian (N = 1 and N = 3 in Fig. 5).

Fig. 5 Viscosity for a fluid made of polymer brushes of various degrees of polymerization (N) and monomeric solvent particles. Notice how the shorter brushes show virtually no shear-thinning, and behave almost like Newtonian fluids. In all cases, L^*_x = L^*_y = D* = 7, and the grafting density is Γ* = 0.30. All quantities are reported in reduced units. Lines are only guides for the eye.

It has been reported that increasing the polymerization degree of PEG brushes on liposomes increases their circulation longevity,²⁹ which would be a desirable aspect to improve drug delivery mechanisms. However, as its degree of polymerization grows, PEG's solubility in aqueous environment is reduced,³⁰ hence there would be a competition between these trends that can be investigated with DPD simulations, varying N. Also, it is important to find out if the viscosity behaves qualitatively the same as a function of shear rate, regardless of the value of N, so that general conclusions are obtained that can be used to interpret a variety of experiments.²⁵ Before presenting our results for the effective mean viscosity of systems where PEG brushes have larger polymerization degrees we comment on the behavior of the mean shear stress on the membranes along the direction of the flow as a function of the shear, 〈F_x([small gamma, Greek, dot above])〉. Following Galuschko and coworkers,⁴⁷ we find a critical shear rate for our systems, [small gamma, Greek, dot above]*, as the value of the shear rate when the behavior of 〈F_x([small gamma, Greek, dot above])〉 changes from linear to sub linear, which depends on N, being [small gamma, Greek, dot above]* ∼ 8.6 × 10⁻³ for N = 14 and [small gamma, Greek, dot above]* ∼ 2.9 × 10⁻³ for N = 25. Finding [small gamma, Greek, dot above]* is tantamount to finding the Weissenberg number, We, since We = [small gamma, Greek, dot above]/[small gamma, Greek, dot above]*, and it signals when shear-thinning behavior starts to set in (We ≥ 1). Fig. 6 shows the dependence of the normalized shear stress, u ≡ 〈F_x([small gamma, Greek, dot above])〉/〈F_x([small gamma, Greek, dot above]*)〉, on We for fluids with brushes with polymerization degrees equal to N = 14, 20 and 25. The data for the three brushes are seen to collapse reasonably well, especially at low We; this is the so-called linear regime, u ∼ We, (dashed blue line in Fig. 6). For larger values of the Weissenberg number, we find that the power law u ∼ We^κ is obeyed with κ = 0.69 (solid black line in Fig. 6) for our systems under theta-solvent conditions. Galuschko et al.⁴⁷ predicted κ = 9/13 ≈ 0.69 using scaling arguments for what they call “dry” polymer brushes, which are brushes where the hydrodynamic and excluded volume interactions are screened out, for fairly concentrated solutions and polymer melts, where the corresponding Flory exponent is ν = 0.5.⁴⁸ Our simulations predict a κ value in excellent agreement with scaling arguments⁴⁷ because under theta-solvent conditions and at the concentrations we modeled, the brushes are almost indistinguishable from the solvent and act as a melt.

Fig. 6 Normalized mean shear stress (u) on the surfaces in the direction of the flow as a function of the Weissenberg number, We, for PEG brushes of three different polymerization degrees, N. In all cases, L^*_x = L^*_y = D* = 7, and the grafting density is Γ* = 0.30. The dashed blue line represents the linear regime, u ∼ We (We < 1), and the solid black line the power law u ∼ We^κ, with κ = 0.69. See text for details.

Recent molecular dynamics simulations⁴⁹ of strongly compressed polyelectrolyte brushes under shear at melt concentrations find that the shear force 〈F_x([small gamma, Greek, dot above])〉 that appears in eqn (6) and (7) behaves at high shear rates ([small gamma, Greek, dot above]) as 〈F_x([small gamma, Greek, dot above])〉 ∼ [small gamma, Greek, dot above]^0.69, in agreement with our predictions in Fig. 6, and with the scaling arguments of Galuschko et al.⁴⁷ Although there are explicit electrostatic interactions in the work of Spirin and Kreer,⁴⁹ they argue that their brush and solvent system behaves as a melt of neutral polymers at high We because of the immobilization of the counterions. Therefore, the scaling exponent κ = 0.69 we obtained in Fig. 6 is very robust.

The effective mean viscosity data of these systems can also be collapsed as a function of We if one defines s = η/η₀,⁴⁷ for We > 1, as seen in Fig. 7(a). The value of the zero shear viscosity η₀ was obtained from extrapolation, using the same procedure as the one described in the discussion of Fig. 5. The data are scattered at low values of the applied shear rate because the signal-to-noise level in those cases is rather poor and one is required to perform very long simulations to get a good estimate of the extrapolated, zero-shear viscosity, as has been reported by others.^20,24 All three cases experience considerable shear thinning as the Weissenberg number is increased, and in this regime the data follow the power law s ∼ We^ζ, with ζ = −0.31, see the black line in Fig. 7(a). To our knowledge, this is the first time such scaling law has been obtained in a theta solvent.

Fig. 7 (a) Normalized effective mean viscosity s = η/η₀ as a function of the Weissenberg number, We, for fluids with PEG brushes of three different polymerization degrees, N. In all cases, L^*_x = L^*_y = D* = 7, and the grafting density is Γ* = 0.30. The line represents the power law, s ∼ We^ζ, with ζ = −0.31. (b) Viscosity of a fluid made up of red blood cells (RBC) dispersed in a 40 kDa liposome–dextran medium, as a function of the shear rate.⁵¹ The blue line represents the power law η ∼ [small gamma, Greek, dot above]^−0.31. See text for details.

In ref. 47 the value ζ = −0.43 was obtained for linear Lennard-Jones brushes in a good solvent. A similar value (ζ = −0.42) was found in simulations of bottle-brush polyelectrolytes in a good solvent.¹⁹ Hence, the polymer architecture appears to be of secondary importance, compared with the quality of the solvent or the polymer concentration. Galuschko and collaborators⁴⁷ proposed also the following relation between κ and ζ:

κ − ζ = 1 (8)

which is clearly fulfilled in our theta solvent simulations. Therefore the relation between the scaling exponents is the same regardless of the solvent quality, and is more fundamental than the architecture of the polymer chains or the particular values of the exponents themselves, separately, as is known to be the case for equilibrium scaling exponents.⁵⁰

In Fig. 7(b) we have included the high shear rate values of the viscosity of human red blood cells (RBC) in a liposomal suspension,⁵¹ which shows evidently shear-thinning. This non-Newtonian behavior is interpreted to be the result of the dissociation of flocculated liposomes at the higher values of the shear rate, an interpretation that is consistent with our predictions in Fig. 7(a). The blue line in Fig. 7(b) represents the power law we have obtained for PEG-grafted surfaces in a theta solvent, which models accurately the results for RBC at high values of the shear rate, providing confirmation for our predictions. Experiments carried out for gels used as vehicles for drug delivering liposomes,⁵² such as hydroxyethylcellulose (HEC) and a mixture of HEC and an acrylic acid-based polymer (carbopol 974), yield values of the viscosity at high shear rates very close to that found in Fig. 7(a), namely −ζ = 0.30 ± 0.01 for HEC and −ζ = 0.310 ± 0.009 for the mixture. The fact that the exponent in these experiments is found to be very close to our prediction suggests that it is of paramount importance to take into account the properties of the solvent (theta conditions), as well as those of the brush (PEG), and lends additional support to these calculations.

Adding polymer brushes to sheared surfaces not only modifies the viscosity of the system but it has been shown to reduce its friction coefficient¹⁰ also. Since this is a parameter that can be directly compared with experiments, we have obtained the friction coefficient of PEG brushes in an explicitly included theta solvent. Using eqn (7) one calculates the mean friction coefficient between the two sliding surfaces. The data in Fig. 8(a) show that the friction coefficient μ at a constant PEG grafting density increases with the shear rate, in agreement with trends found by others;^19,47,53 it increases also with N because the overlap between the brushes also grows. Moreover, the μ values for PEG brushes in a theta solvent are somewhat larger than those obtained in simulations of linear brushes in a good solvent;¹⁸ this is of course expected because the brush thickness is reduced when it is immersed in a theta solvent, and the contacts between the polymers grow when going from a good to a theta solvent.² Experiments by Klein and collaborators⁵⁴ on the friction of polymer brushes in a theta solvent compared to a good solvent show that the friction coefficient can be larger in the former by up to three orders of magnitude than in the latter. The reason relies on the fact that under theta solvent conditions the brushes can interdigitate and the solvent can penetrate the brushes (see Fig. 2), reducing the freely flowing solvent in the center of the channel, thereby increasing the shear force and the friction. Klein and coworkers⁵⁴ obtain values for μ which are close to our predictions, shown in Fig. 8(a). When electrostatic interactions are included the values for μ obtained from computer simulations turn out to be slightly larger and closer to 1 for good solvent conditions at high values of the shear rate,^19,53 due to the added osmotic pressure of the ions. However, in biological systems like those modeled here, the quality of the solvent plays a key role in providing a mechanism that promotes a lower viscosity environment.

Fig. 8 (a) Friction coefficient for polymer brushes of various degrees of polymerization (N), obtained using eqn (7). In all cases, L^*_x = L^*_y = D* = 7, and the grafting density is Γ* = 0.30. Lines are only guides for the eye. (b) The friction coefficient data shown in (a), normalized by its value at [small gamma, Greek, dot above]* for the largest 3 brushes, μ/μ*, as a function of the Weissenberg number, We = [small gamma, Greek, dot above]/[small gamma, Greek, dot above]*. The black line represents the fit to the power law μ/μ* ∼ We^κ, with κ = 0.69. The inset shows the friction coefficient for the polymer with the smallest degree of polymerization (N = 7, pink rhombi), with the blue line representing a linear fit, μ ∼ [small gamma, Greek, dot above].

In Fig. 8(b) we show the same data as in Fig. 8(a) on logarithmic scales to emphasize the power law behavior. For the shortest brush (see the inset of Fig. 8(b)), made up of chains with polymerization degree equal to N = 7, there is a linear relation between the friction coefficient and the shear rate, μ ∼ [small gamma, Greek, dot above], for the entire range of shear rate values used in the simulations (blue line). To see sublinear behavior here one would have to impose larger shear. For the other polymerization degrees shown (N = 14 to 25) the main panel in Fig. 8(b) shows that the data collapse rather well on a single curve when the friction coefficient (μ) is normalized by its value (μ*) at the start of the shear-thinning regime ([small gamma, Greek, dot above]*). The black line is the fit to the power law μ/μ* ∼ We^κ, with κ = 0.69. This is the same exponent as the one obtained from the normalized shear stress data, in Fig. 6, indicating that the force 〈F_z([small gamma, Greek, dot above])〉 (see eqn (7)) is essentially independent of shear once the brushes have reached a certain length. At the same time, the fact that the exponent is κ = 0.69 at We > 1 in all three cases (N = 14, 20 and 25) shows that 〈F_x([small gamma, Greek, dot above])〉 (see eqn (7)) is responsible for the behavior of the friction coefficient in the high shear rate regime. The same behavior is found for Lennard-Jones brushes,⁴⁷ namely the friction coefficient dependence on We scales with the same exponent as that of the shear stress, although in such cases κ = 0.57 because that system is under good-solvent conditions. In ref. 55 we have included a table listing all the data shown in Fig. 8(a) and (b).

IV. Conclusions

We have shown that complex fluids with polymer brushes under theta solvent conditions display rheological characteristics that differ in detail from their good solvent counterparts, but that obey nevertheless the same general scaling properties, in particular for the shear stress and the viscosity. The exponent for the shear stress as a function of We obtained from our simulations, κ = 0.69, is in excellent agreement with the value predicted using scaling arguments⁴⁷ for dense polymer brushes where excluded volume interactions are screened out, as in our model. The viscosity we obtained scales with the shear rate with an exponent equal to ζ = −0.31, which reproduces remarkably well measurements of the viscosity of red blood cells dispersed in a liposome carrying aqueous fluid at a high shear rate,⁵¹ as well as other experiments.⁵² The friction coefficient data as a function of shear rate for different polymerization degrees of the chains making up the brushes were found to collapse on a universal curve whose behavior at We > 1 follows the same power law as the shear force, i.e., μ ∼ We^κ, with κ = 0.69, in agreement with trends found in simulations under good-solvent conditions,⁴⁷ albeit with a value of κ particular to those conditions. Ours are the first simulations, to the best of our knowledge, of the scaling of viscosity and the friction coefficient for systems under theta-solvent conditions.

It is argued that our simulations are useful for understanding the behavior of biopolymer brushes coating drug-carrying liposomes in an aqueous environment that acts as a theta solvent, under shear. Increasing the thickness of the brush through the degree of polymerization or the shear velocity on the surfaces is shown to raise the friction coefficient, which is a deleterious effect for the transport properties of these liposomes. However, when liposomes are covered by PEG brushes in a theta solvent at high shear rates, their flowing characteristics make them optimal carriers for drug delivery, as the polymer brushes imprint them with efficient injectable characteristics (low viscosity at high shear rates) while at the same time providing them with thermodynamically stable (“stealth”) mechanisms. These simulations have the additional advantage of including hydrodynamic interactions, as well as the solvent explicitly, and being mesoscopic they reproduce length and time scales comparable with those of environments of biological interest. Therefore, we believe our work should be useful in the improved design of drug-carriers, the rheological characterization of sheared brushes, and in establishing general scaling laws for non-equilibrium polymer brushes.

Acknowledgements

The authors would like to thank ABACUS, CONACyT grant EDOMEX-2011-C01-165873, for funding. AGG thanks M. A. Balderas Altamirano, C. Carmín, E. de la Cruz, J. P. López Neria and E. Pérez for instructive discussions.

References

1. D. H. Napper, Polymeric Stabilization of Colloidal Dispersions, Academic Press, London, 1983.
2. J. N. Israelachvili, Intermolecular and Surfaces Forces, Academic Press, New York, 3rd edn, 2011.
3. A. Mollaei and B. Maini, J. Can. Pet. Technol., 2010, 49, 65.
4. M. A. Cohen Stuart, W. T. S. Huck, J. Genzer, M. Müller, C. Ober, M. Stamm, G. B. Sukhorukov, I. Szleifer, V. V. Tsukruk, M. Urban, F. Winnik, S. Zauscher, I. Luzinov and S. Minko, Nat. Mater., 2010, 9, 101.
5. S. Iyer, R. M. Gaikwad, V. Subba-Rao, C. D. Woodworth and I. Sokolov, Nat. Nanotechnol., 2009, 4, 389.
6. Polymer Brushes, ed. R. Advincula, W. J. Brittain, K. C. Caster and J. Rühe, Wiley-VCH, 2004.
7. P. G. de Gennes, Macromolecules, 1980, 13, 1069.
8. S. Alexander, J. Phys., 1977, 38, 983.
9. S. T. Milner, T. A. Witten and M. E. Cates, Europhys. Lett., 1988, 5, 413.
10. J. Klein, E. Kumacheva, D. Mahalu, D. Perahia and L. J. Fetters, Nature, 1994, 370, 634.
11. L. Han, D. Dean, P. Mao, C. Ortiz and A. J. Grodzinsky, Biophys. J., 2007, 92, 1384.
12. J. Klein, Science, 2009, 323, 47.
13. S. Weinbaum, J. M. Tarbell and E. R. Damiano, Annu. Rev. Biomed. Eng., 2007, 9, 121.
14. A. Samal, Y. Sultana and M. Aqil, Curr. Drug Delivery, 2007, 4, 297.
15. Proteoglycans: Structure, Biology and Molecular Interactions, ed. R. V. Iozzo, CRC Press, New York, 2000.
16. P. Y. Lai and C. Y. Lai, Phys. Rev. E: Stat. Phys., Plasmas, Fluids, Relat. Interdiscip. Top., 1996, 54, 6958.
17. T. Kreer, M. H. Müser, K. Binder and J. Klein, Langmuir, 2001, 17, 7804.
18. F. Goujon, P. Malfreyt and D. J. Tildesley, Mol. Phys., 2005, 103, 2675.
19. J. M. Carrillo, W. M. Brown and A. V. Dobrynin, Macromolecules, 2012, 45, 8880.
20. C. Pastorino, K. Binder, T. Kreer and M. Müller, J. Chem. Phys., 2006, 124, 064902.
21. F. Goujon, P. Malfreyt and D. J. Tildesley, Soft Matter, 2010, 6, 3472.
22. M. Deng, X. Li, H. Liang, B. Caswell and G. E. Karniadakis, J. Fluid Mech., 2012, 711, 192.
23. D. Irfachsyad, D. J. Tildesley and P. Malfreyt, Phys. Chem. Chem. Phys., 2002, 4, 3008.
24. C. Pastorino, K. Binder and M. Müller, Macromolecules, 2009, 42, 401.
25. K. Binder and A. Milchev, J. Polym. Sci., Part B: Polym. Phys., 2012, 50, 1515.
26. F. Alarcón, E. Pérez and A. Gama Goicochea, Eur. Biophys. J., 2013, 42, 661.
27. Stealth Liposomes, ed. D. D. Lasic and F. Martin, CRC Press, Florida, 1995.
28. P. L. Hansen, J. A. Cohen, R. Podgornik and V. A. Parsegian, Biophys. J., 2003, 84, 350.
29. C. Allen, N. Dos Santos, R. Gallagher, G. N. C. Chiu, Y. Shu, W. M. Li, S. A. Johnstone, A. S. Janoff, L. D. Mayer, M. S. Webb and M. B. Bally, Biosci. Rep., 2002, 22, 2.
30. C. Ö. Dinç, G. Kibarer and A. Güner, J. Appl. Polym. Sci., 2010, 117, 1100.
31. A. Bunker, Phys. Procedia, 2012, 34, 24.
32. P. J. Hoogerbrugge and J. M. V. A. Koelman, Europhys. Lett., 1992, 19, 155.
33. T. Murtola, A. Bunker, I. Vattulainen, M. Deserno and M. Karttunen, Phys. Chem. Chem. Phys., 2009, 11, 1869.
34. R. Suresh, S. N. Borkar, V. A. Sawant, V. S. Shende and S. K. Dimble, International Journal Pharmaceutical Drugs Sciences Technology, 2010, 3, 901.
35. P. Español and P. Warren, Europhys. Lett., 1995, 30, 191.
36. D. Frenkel and B. Smit, Understanding Molecular Simulation, Academic, New York, 2002.
37. R. D. Groot, J. Chem. Phys., 2003, 118, 11265.
38. R. D. Groot and P. B. Warren, J. Chem. Phys., 1997, 107, 4423.
39. A. Gama Goicochea, Langmuir, 2007, 23, 11656.
40. M. Murat and G. S. Grest, Phys. Rev. Lett., 1989, 63, 1074.
41. A. Gama Goicochea, M. Romero-Bastida and R. López-Rendón, Mol. Phys., 2007, 105, 2375.
42. I. Vattulainen, M. Karttunen, G. Besold and J. M. Polson, J. Chem. Phys., 2002, 116, 3967.
43. D. D. Verma, S. Verma, G. Blume and A. Fahr, Int. J. Pharm., 2003, 258, 141.
44. D. Needham, K. Hristova, T. J. McIntosh, M. Dewhirst, N. Wu and D. D. Lasic, J. Liposome Res., 1992, 2, 411.
45. M. Doi and S. F. Edwards, The Theory of Polymer Dynamics, Oxford University Press, 1988.
46. M. E. Velázquez, A. Gama Goicochea, M. González-Melchor, M. Neria and J. Alejandre, J. Chem. Phys., 2006, 124, 084104.
47. A. Galuschko, L. Spirin, T. Kreer, A. Johner, C. Pastorino, J. Wittmer and J. Baschnagel, Langmuir, 2010, 26, 6418.
48. P. G. de Gennes, Scaling Concepts in Polymer Physics, Cornell University Press, New York, 1979.
49. L. Spirin and T. Kreer, ACS Macro Lett., 2013, 2, 63.
50. N. Goldenfeld, Lectures on Phase Transitions and the Renormalization Group, Addison Wesley Publishing Company, 1992.
51. H. Sakai, A. Sato, N. Okuda, S. Takeoka, N. Maeda and E. Tsuchida, Am. J. Physiol., 2009, 297, 583.
52. S. Mourtas, S. Fotopoulou, S. Duraj, V. Sfika, C. Tsakiroglou and S. G. Antimisiaris, J. Colloid Interface Sci., 2008, 317, 611.
53. J. M. Carrillo, D. Russano and A. Dobrynin, Langmuir, 2011, 27, 14599.
54. J. Klein, E. Kumacheva, D. Perahia, D. Mahalu and S. Warburg, Faraday Discuss., 1994, 98, 173.
55. See ESI.†.

---

Footnote

† Electronic supplementary information (ESI) available. See DOI: 10.1039/c3sm52486h

---