Self-consistent field description of polyelectrolyte-grafted colloidal actuators

Abstract

We present a theoretical description of actuators in prototype artificial muscle tissue by means of a self-consistent (mean field) lattice computational scheme. The actuators are composed of pH-responsive polyelectrolytes grafted at both ends between plate- or rod-like colloidal particles and immersed in an aqueous solution. We build on a model developed for grafted rods [A. Ianiro, J. A. Berrocal, R. Tuinier, M. Mayer and C. Weder, J. Chem. Phys., 2023, 158, 14901], but we specifically include weakly acidic monomers to incorporate the effects of a pH variation to trigger the expansion of the material. As a first toy model, we consider strong polyelectrolyte chains: for both plate- and rod-like colloidal particles, we obtain pressure differences of the order of tens of MPa; sufficient to generate volume variation. During actuation, the system expands and contracts by approximately one third of the polymer contour length and about 100kT of work per polymer chain is performed. Secondly, we show that for weakly charged polyelectrolyte chains the salt concentration can be used to tune the actuation window by multiple pH units, which is important to obtain a biocompatible range of pH values. In this scenario, we find smaller actuation pressures than for strong polyelectrolytes, but still sufficient expansion and contraction for practical purposes.

---

1 Introduction

We present a theoretical description of a prototype colloidal actuator or artificial muscle. Artificial muscles are actuating materials that can expand, contract, or change conformation to generate displacement when subjected to external input.¹ Because of this ability to change their shape or volume, artificial muscles are useful for a wide variety of applications in biomedical engineering and soft robotics. With this purpose, a great variety of materials have been studied: in most cases, these are soft and biocompatible materials, and the nature of the external input which triggers actuation varies greatly, depending on the material and the principle of the actuation. Possible triggers are temperature, pH, light, or electric and magnetic fields.^2–4 Examples of soft and biocompatible materials include molecular motors,^5,6 shape-memory polymers,⁷ or electroactive polymers.^7,8 Shape-memory polymers contain some segment in the chain that can reversibly change shape when crossing some typical transition temperature, while electroactive polymers can undergo deformations in response to an electrically induced transport of ions and/or solvent, or in response to Coulomb forces, which have been employed in tissues which mimic biological muscles, as well as in robotics.^7,8

Liquid crystal elastomers are another example of possible constituents of artificial muscles.^9,10 In addition to the elastic properties of the polymers, they can (reversibly) form nematic phases upon a temperature change.⁹ Elastomers can be used as part of wearable devices with prosthetic or rehabilitation purposes, or in components of bio-mimetic robots, which can manipulate objects.^11,12

Hydrogels (composed of water-swollen cross-linked polymers) are also a good candidate for artificial muscles, due to their softness and their ability to change their conformation, taking up or releasing water, by responding to external stimuli such as a light, temperature, or pH.^13–15 Differently from the swelling of hydrogels, the model actuator system described in this paper is based on the idea of generating actuation using linear stimuli-responsive polymers grafted to colloidal particles.^13–17

Linear stimuli-responsive polymer chains, which can swell or contract when changing the conditions of the surroundings, have, for example, been used to control the diffusion of ions through nano-pores,¹⁸ to make surface properties responsive to an electrical potential,¹⁹ but also to induce displacement of nanoparticles patched to the free end of the polymer layer.^20–22 Santer and Rühe²² made use of a layer of polymers tethered to a hard surface at one end (see Section 2.2). These polymer chains consisted of either di- or tri-block copolymers or of a mixture of linear chains of two different kinds A and B, and block-copolymers composed of A and B blocks. The authors induced a different topography on the polymer layer by switching between two solvents. When using a good solvent for both blocks, the surface was a more flat substrate; when using a selective solvent (good solvent for only one of the blocks and poor for the other block), the polymers reorganized such that the more soluble block would be more exposed to the solvent, making the top of the layers more irregular. As a consequence, a particle patched at the free end of the polymer layer moved when the topography switched.

In contrast, in our present work, the stimuli-responsive properties of the polymer chain are used to generate a directional displacement and perform work. We envision a tissue in which the actuation units are composed of linear stimuli-responsive polymers, grafted at both ends to two hard anisotropic particles. Through an external stimulus (in this case, pH variation), the polymer chains expand or contract, taking up or expelling solvent, and generating a directed displacement. If the actuating units within the tissue are properly aligned, the anisotropy of the colloids could provide directionality to the expansion and contraction of the full tissue too.

We explore two geometries, with rod-like and plate-like colloidal particles. Fig. 1 shows a schematic representation of both these systems. The doubly-grafted chains are modelled as weak acidic polyelectrolytes, where the degree of dissociation of the acidic groups along the backbone of the polymer depends on the pH of the solution and can be influenced by other factors such as the salinity of the solution.²³ An example of a weakly charged polyelectrolyte is poly(methacrylic acid) (PMAA), which contains carboxylic acid groups that are partially deprotonated at neutral pH, but have a lower degree of dissociation at lower pH values, causing the polymer to become more hydrophobic in a more acidic solution.^23,24 As a “toy model”, we also consider the case of strong anionic polyelectrolytes (also known as quenched polyelectrolytes, e.g., partially sulfonated polystyrene).²³

Fig. 1 Schematic representation of the cylindrical system with rod-like colloidal particles (at the top, cross section) and planar system with plate-like colloidal particles (at the bottom).

For the colloidal components of the actuator system, the particles should have an anisotropic shape to obtain directional actuation, and a relatively large surface area to allow grafting of the polyelectrolyte chains. We therefore focused on rod-shaped colloids and colloidal platelets. In the former case, we envision, for instance, cellulose nanocrystals (CNCs),^25–27 or silica rods,^28,29 which can form different liquid crystal phase states. In the latter case, it would be possible to use clay materials such as platelets of montmorillonite.^30–32

Theoretical descriptions of grafted polyelectrolytes have been previously reported for singly end-grafted chains, where the polymers are modelled as strong (or quenched) polyelectrolytes, carrying a fixed amount of charge^33–36 or as weak (or annealing) polyelectrolytes, having a variable degree of dissociation, depending on (local) pH.³⁷ These studies provide results for the structure as well as for the equilibrium thickness of the polymer layer at the surface by making use of self-consistent field (SCF) theory calculations.

We apply SCF theory to quantify the performance of a single actuating unit of this prototype artificial muscle, in order to study its actuation properties and possibly to optimize them. We start from a model similar to the one presented in ref. 38, but we specifically consider pH-responsive polyelectrolyte chains and incorporate a more realistic description of the deprotonation of the acidic groups on these chains. For the weak polyelectrolyte case, we explicitly model the dissociating character of the weakly charged monomers of the polyelectrolytes. This allows us to account for pH variations as a trigger for the expansion mechanism by introducing a local degree of dissociation, depending on the local pH. We also analyse a simpler toy model of linear strong polyelectrolytes. We model these polymer chains as carrying a fixed amount of monomers with charge −e, and we model the deprotonation process by varying the number of these negatively charged monomers, to study the actuation properties as a function of the amount of charge in the polymer chains.

Furthermore, we extend the approach described in ref. 38 to colloidal platelets, and we compare the models of weak and strong polyelectrolytes in planar and cylindrical geometries and study the effects of different parameters (such as salt concentration) on the actuation properties. Section 2 summarizes the theoretical background on polyelectrolyte chains grafted to a hard surface. Subsequently, the main results and actuation properties predicted by our model are described and discussed in Section 3. In the case of weakly acidic polyelectrolytes results obtained with these computational method are compared to previous experimental work.^17,39 Finally, some comments are made about the limitations of the model and what could be improved.

2 Theoretical approaches

2.1 Self-consistent field theory lattice computations

In this section, the modeling of the actuation system using self-consistent field theory is described. First, Scheutjens–Fleer SCF (SF-SCF) theory is discussed in Section 2.1.1. Section 2.1.2 focuses on the application of SF-SCF in this work, and Section 2.1.3 describes the parameters used to predict the actuation properties and the way they are computed. Section 2.2 summarizes some scaling results of a polyelectrolyte layer grafted to a hard surface.

2.1.1 Scheutjens–Fleer approach. The properties of the system are determined by means of numerical SCF lattice computations as developed by Scheutjens and Fleer.⁴⁰ Here we provide a brief summary of this model; a more detailed explanation of SF-SCF can be found in ref. 40–45.

The system is modelled on a lattice, and the polymers are represented as linear chains whose segments occupy one lattice site each. Other components of the system (such as solvent or salt ions) are also described in terms of beads that occupy one lattice site. We use the complete filling constraint, which means all lattice sites are occupied.

The equilibrium configuration is computed using a mean-field approach: the system is described using the volume fraction profiles for each component as a function of the space coordinates φ(x), where x is a point in the 3D (lattice) space. The interactions between different components and the constraints imposed by the boundary conditions result in a potential field u(x)/kT. In the following, it is assumed that the concentrations are constant in two dimensions (for example, x and y), so they depend only on the z-direction.

The density profiles φ_i(z) (where i is an index running over the different components of the system) directly depend on the potentials, since the latter determine the probability for a single monomer or component to be in a specific layer. In fact, in this case the statistical weight corresponds to the Boltzmann probability e^−u(z)/kT.

When considering a monomer number s along a chain, the probability of having that monomer in a specific layer z also depends on the position of the previous monomers in the chain. Therefore, in a polymer chain, this probability is expressed through a propagator (or statistical weight) G(z, s), which is a function of the spatial coordinate z and of the position of the monomer s along the chain: G(z, s) = G(z)〈G(z, s − 1)〉. The latter term 〈G(z, s − 1)〉 is an average over the current layer and the two adjacent ones, which takes into account the number of nearest neighbours that the (s − 1)th segment has in layer z − 1, z, and z + 1. In general, if segment s − 1 has a fraction λ₀ of nearest neighbours in the same layer z and a fraction λ₊ and λ₋ in the two adjacent layers z − 1 and z + 1, 〈G(z, s − 1)〉 will be:

〈G(z, s − 1)〉 = λ₋G(z − 1, s − 1) + λ₀G(z, s − 1) + λ₊G(z + 1, s − 1). (1)

For example, in a simple cubic planar lattice each segment has 4 nearest neighbours in the same layer, 1 in the layer before and 1 in the following layer:

If the lattice is not planar, this average has to be corrected for the number of lattice sites in each layer. In a cylindrical configuration, for example, the number of lattice sites in each layer increases along the radial direction. As a consequence, the amount of contact a lattice site has with neighbouring sites also increases along the radial direction (here for convenience also indicated by the letter z). The factors λ have to be multiplied by the contact surface area between adjacent layers per unit length of the cylinder. The number of sites in a cylindrical concentric layer z is proportional to L(z) = π(z² − (z − 1)²), so this surface area in layer z is S(z) = 2πz. Therefore:

λ^cyl₀(z) = 1 − λ^cyl₊(z) − λ^cyl₋(z). (3c)

Once G(z, 1) is known from the potential u, it is possible to compute G(z, s) for any segment s along the chain. From that, the volume fraction of any interior monomer s in layer z can be computed by considering any segment as the end point of two walks: the first, starting at s = 1, has s segments, while the second, starting at s = N, has N − s + 1 segments. Therefore, the volume fraction of monomer s in layer z is proportional to the product of the two distributions of the two (partial) chains (composition law):

However, if the chains are grafted at one end, the first segment must be in the first layer, so the statistical weight of segment s = 1 is:
The segment statistical weight and the volume fraction profile of the grafted chains are therefore:In the case of a linear chain grafted on both sides, the same reasoning can be applied also for the last segment of the chain s = N. The volume fraction in layer z can then be obtained by summing over all φ(z, s): .

Potential profiles u(z) can also be expressed as functions of the volume fraction profile, depending on the type of interactions in the system. For example, in the system considered in this work, the potential of a component i has three contributions. Using z as the spatial coordinate as before, the potential profile for component i reads:

where in eqn (6) j runs over all other components of the system. e is the elementary charge, q is the valence of the component i and ψ(z) is the electrostatic potential. The first term is a Lagrange parameter which is introduced to meet the filling constraint, implying: in each layer z. The second Flory–Huggins term describes the excluded volume and non-electrostatic interactions, summarized in the interaction parameter χ between component i and component j. φ^b_j is the bulk concentration of component j, which is 0 in the case of the grafted polyelectrolytes. The last term describes the electrostatic interactions. A more detailed description of how electrostatic interactions are included in the SCF framework can be found in ref. 36, 37, 45 and 46.

The equilibrium configuration is found by minimization of the free energy (or another thermodynamic potential, depending on the ensemble description adopted), which is expressed as a function of these concentrations and potential profiles. The minimization constraints with the addition of the filling constraint result in a system of coupled self-consistent equations (as concentrations depend on the potentials and vice versa), which is solved numerically using a Newton-like algorithm. For this, we use the software package sfbox developed at Wageningen University. A new open source version of the code is available on Github.⁴⁷

2.1.2 Description of the system. We use the SCF computational framework to describe the actuating properties of a single unit of artificial muscle tissue, which is composed of a colloidal particle (either a platelet or a rod) and doubly end-grafted polyelectrolyte chains in an aqueous solution containing salt and acid/base. A schematic representation of the systems is shown in Fig. 2. An actuation unit is modelled on either a cylindrical lattice (resembling a rod grafted with polyelectrolyte chains) or a planar lattice (resembling a platelet grafted with polyelectrolyte chains). In the first case, a hard cylindrical surface is placed before the innermost layer, at the centre of the lattice: periodic boundary conditions are applied along the cylinder axis, while mirror boundary conditions are applied in the radial direction, to mimic the presence of other grafted colloidal particles around the single unit that is considered. The lattice layers are concentric cylindrical shells surrounding the surface. In the case of planar geometry, a three-dimensional cubic lattice is considered, with a planar surface at one end representing the platelet and a mirror boundary condition at the opposite end. Lattice layers are slabs parallel to the surface. Concentrations of all components are assumed constant within one lattice layer, so that variations are only possible in the radial direction in the cylindrical system and in one spatial direction (perpendicular to the surface) in the planar one.

Fig. 2 Schematic representation of the cylindrical (a) and flat (b) systems. The black dashed line represents the mirror boundary condition, while the surface is indicated in purple. The lines denote the boundary between the layers. Polymers are represented as chains of alternating light blue beads (carboxylic acid groups) and dark blue beads (rest of the monomer), while green and orange beads stand for salt anions and cations. The remaining lattice sites in the system are filled with solvent. The lattice sites themselves and solvent are not drawn for clarity.

Polymer chains are end-grafted to the surface and to the mirror at the end of the system. We compare two different approaches to describe the electrostatic properties of the polyelectrolyte chains. In the first toy model, a certain fraction of the polymer beads has a fixed valence −1, and we study the actuation properties of the material as a function of the fraction of charged monomers. In the second approach, these chains are described as weak polyelectrolytes. To be consistent with the approach of previous work,³⁸ the chain is composed of alternating weakly acidic monomers HA/A⁻ and neutral beads. The weakly acidic monomers undergo an acid–base reaction: HA + H₂O ⇌ H₃O⁺ + A⁻. In each calculation, the acidic groups' pK_a and the bulk concentration of [H₃O⁺] are fixed parameters, while the degree of dissociation can vary throughout the system. The bulk concentration of H₃O⁺ is varied in the different calculations. By applying pH = −log([H₃O⁺]), it is thus possible to estimate the actuation properties as a function of the pH of the system. The pK_a (of dilute, separate monomers) is set to 4.5 (typical value for carboxylic acid groups).⁴⁸ A different pK_a would lead to a shift in the pH values necessary to generate the actuation. The pK_a value based on molar concentrations has to be converted into a ‘lattice’ value, as described in the appendix. For a more detailed description of how the dissociation degree is calculated through the whole system, we refer to ref. 46 and 49–52. Ref. 49 presents an extension of the original SCF theory of Scheutjens and Fleer,⁴⁰ describing the case where monomers along a polymer chain can undergo a reaction changing their state, such as an acid–base reaction. The code we used to solve the equation is however based on the description in ref. 46.

The lattice constant is chosen as b = 0.3 nm.³⁸ The input parameters of the system are the grafting density of the chains σ (number of chains per unit area), the degree of polymerisation N_P, the Flory–Huggins parameter for the interaction between the water and the polymer segments χ and the molar salt concentration c_salt.

In order to keep the same grafting density σ of polymers on the colloidal surface the same for both geometries, the number of polymer chains N_chains in the system has to be set accordingly. Therefore, in the flat geometry, N_chains = b²σ and in the cylindrical geometry, N_chains = 2πrb²σ (where r is the radius of the cylinder in units of number of layers).

To set the salt concentration, the concentration of either cations or anions is fixed to c_salt, whereas the other salt species is computed to meet the constraint that the whole system is neutral.

2.1.3 Computation of actuation properties. Since the number of grafted polymers in the system is fixed and the salt concentration is determined by its fixed chemical potential, the calculations are performed in the semi-grand canonical ensemble. The equilibrium configuration can be found by minimizing the partial open free energy or semi-grand potential , where F is the free energy, μ_j the chemical potential of component j, N_j the number of molecules of the non-grafted species j.

The partial open free energy is computed for systems with varying numbers of lattice layers, in order to study its dependency on the system size. The number of lattice layers is varied from a maximum, corresponding to the configuration in which the polymer chains are completely extended, to a minimum, where the chains are completely compressed and no solvent is present in the system. To quantify the extension of the system, a normalized system strain is defined as δ = Λ/(Λ_max − Λ_min), where Λ indicates the number of lattice layers.

From the semi-grand potential, the disjoining pressure as a function of system strain can be derived by Π = −∂Ω/∂V|_μjNiT, where Ω is the partial open free energy and V the volume of the system, μ_j is the chemical potential of the non grafted species (index j runs over the non grafted species), N_i is the number of grafted components (index i runs over the grafted components of the system), and T is the temperature. In the lattice system, this is computed by taking the ratio between the change in Ω and the change in volume when removing one layer from the system. The equilibrium strain for a system at rest is found via Π(δ) = 0.

The process to compute the actuation properties is exemplified in Fig. 3. The pressure computation is performed between two different conditions (generally indicated by A and B in the next paragraph), such as different values of charge or pH.

Fig. 3 Example of the pressure curves and work performed during a cycle computation for polyelectrolyte chains grafted at a planar surface. The pressure is presented as a function of the system strain parameter , where Λ indicates the number of lattice layers. System parameters: N_P = 100, χ = 0.4, c_s = 0.001 mol L⁻¹, σ = 1 nm⁻².

Starting from the equilibrium state at Π^B(δ_B) = 0, a sudden increase in pH will cause deprotonation and expansion of the system. The actuation pressure upon deprotonation is computed as the difference between the two pressure curves at fixed strain ΔΠ_dep = Π^A(δ_B) − Π^B(δ_B). The system then relaxes by reversibly expanding to the new equilibrium strain under the new conditions, after which the pH is changed back again to compute the protonation pressure: ΔΠ_pro = Π^B(δ_A) − Π^A(δ_A) (which has a negative value). The system finally relaxes back to its initial conformation. One important assumption is that the pH change—and the resulting pressure jump—happens instantly, i.e., the system remains at the strain of the initial state while the pH change occurs, and only changes in volume afterwards.

In addition to the actuation pressure jumps, the maximum (reversible) work that can be performed by the system in the expansion–contraction cycle is a useful quantity to quantify the actuation properties. This work is computed by integrating the difference between the two pressure curves over the volume (from one equilibrium value δ_B to the other δ_A), and is normalized over the number of polymer chains in the system: . Here, the work performed by the chains is defined as a positive quantity. The work is indicated in Fig. 3 as the shaded area.

2.2 Analytical approximations

The equilibrium thickness of polyelectrolyte chains tethered to a surface can also be described by means of analytical theories. The extension of these chains from the surface is influenced by the grafting density, which determines the distance between the grafted groups. In the so-called mushroom regime, the chains are sparsely grafted so that they are far enough from each other to not interact; each chain will just coil up along the surface. However, if the grafting density of the chains is high enough, the coils will stretch to minimize contact between segments of different chains. Such a system of densely grafted polymer chains is called a polymer brush, and it can be compared to the polymer layer in our system, where the chains are grafted at both ends. In this section, two analytical models for strong and weak polyelectrolytes are described and compared.

An analytical description for a polyelectrolyte brush with a fixed amount of charge of the polymer chains at a single solid surface is provided in ref. 33–36. According to this framework, the polymer concentration is approximated as being constant as a function of the distance from a solid surface. Consequently, the number density ρ of polymer segments in the grafted polymer layer and the thickness of the grafted chain H (brush thickness) are related by the expression: ρ = Nσ/H, where σ is the anchor density. This approximation is defined as the ‘box model’. In Section 3.2, we compare our results with predictions from this model to investigate the (qualitative) difference between single- and double-grafted systems.

Within this framework, the brush structure depends on the system conditions, and it is possible to identify different regions of the system parameters where the brush thickness shows different trends.³⁶ The most important parameter in determining the transition from one regime to another is the amount of salt present in the solution. This is expressed by the volume fraction of the salt φ_s, which is the volume occupied by the salt molecules over the total volume of the system. In a lattice scheme, this is related to the salt concentration c_s (in mol L⁻¹) by: φ_s = c_s·10³N_Avb³, where N_A is the Avogadro number, b the lattice size in m³ and the factor 10³ is necessary to convert m³ into L. The chosen value of φ_s determines whether the polymer–solvent interactions dominate over the electrostatics and the osmotic pressure of the counterions in the system. More specifically, at low salt concentrations, the Debye screening length is high, and deprotonation of acidic groups on the chain would lead to strongly unfavourable electrostatic repulsion. Therefore, the amount of charge on the polyelectrolyte is still negligible, and the brush can be treated as neutral (neutral brush regime). In this case, the scaling relation for neutral brushes derived using the box model can be applied, and the brush thickness is given by H ∼ N(νσ)^1/3 (with ν the excluded volume parameter).^53,54

When the salt concentration increases, the osmotic brush regime is reached: the salt screens the interaction between the weakly acidic groups, which start to dissociate, leading to an increase in the charge fraction in the brush. In this regime, the resulting counterions are still trapped inside the brush, causing it to swell under the osmotic pressure; the relation between H and the charge density can be expressed as:

H ∼ Nm^−1/2, (7)

where m is the distance between two charged groups inside the brush. Finally, in the salted brush regime, the excess of salt screens the electrostatic repulsion between charged groups on the chains, leading to a decrease in the brush thickness. The equilibrium thickness then scales as H ∼ σ^1/3m^−2/3φ_s^−1/3.

A similar division in brush regimes can also be used to predict the equilibrium properties of a weakly charged single-grafted polyelectrolyte brush.³⁷ In this case, the brush equilibrium height H is expressed as a function of the degree of dissociation α rather than as a function of the fraction of charged monomers along the chain f = 1/m. In the neutral regime both α and f are 0. With increasing salt concentration, the osmotic brush regime is reached and the degree of dissociation increases following a sigmoid curve given by:

with φ_H3O⁺ the bulk volume fraction of H₃O⁺ and σ the grafting density³⁷). When the salt concentration is further increased, the brush enters the salted brush regime and α ≈ 1. The boundaries of these different regimes are determined by the grafting density σ, the pK_a of the acidic groups and the pH of the system (via the bulk concentration of H₃O⁺). Fig. 4 shows the degree of dissociation α as a function of the salt volume fraction φ_s for different pH values for two different systems for pK_a values 4.5 (solid curves) and 7.0 (dashed curves). The left graph was obtained using parameters which describe our system: grafting density σ = 1 nm⁻², N_p = 200, χ = 0.4, pK_a = 4.5 (solid curves) or pK_a = 7.0 (dashed). For comparison, the dashed curves in the right figure refer to polymer chains with pK_a = 7.0 and grafting density σ = 0.0056 nm⁻²: these trends agree with those of Israëls et al.³⁷

Fig. 4 Salt concentration dependence of the dissociation degree of the weakly acidic groups on the polymers α for different pH values, in a planar configuration. Left: N = 200, σ = 1 nm⁻², χ = 0.4, lattice constant 0.3 nm, pK_a = 4.5 (solid curves, closed symbols) or pK_a = 7.0 (dashed curves, open symbols). Right: N = 500, σ = 0.0056 nm⁻², χ = 0, lattice constant 0.6 nm, pK_a = 4.5 (solid curves, closed symbols) or pK_a = 7.0 (dashed curves, open symbols).

Furthermore, scaling relations can also be derived for the brush equilibrium thickness H in the different regimes, as in the case of the strong polyelectrolyte. The osmotic brush height H scales as , while in the salted brush regime the equilibrium thickness scales as H ∼ (σν_eff)^1/3, so the brush can then be described as a neutral brush with a higher ‘effective’ excluded volume parameter ν_eff = ν + α/φ_s.

3 Results and discussion

This section is divided into three parts: firstly, results for the toy model of strong polyelectrolytes are presented. Secondly, a comparison is drawn between our weak polyelectrolyte model with double-grafted chains and previous studies on weak polyelectrolytes,³⁷ and finally, the last subsection is dedicated to our results for the actuation properties using weak polyelectrolytes.

3.1 Toy model: strong polyelectrolytes

We will first discuss polyelectrolytes with a fixed amount of charge. Here, a negative charge −e is assigned to a fraction f = 1/m of the polymer beads, which is varied. The system strain, the actuation pressure between neutral and partially charged, and the work performed in one expansion–contraction cycle are computed as functions of the fraction of charged segments along the polymer chain.

In Fig. 5, results for the equilibrium strain δ_eq, disjoining pressure ΔΠ and work per polymer chain W_chain are plotted for different values of salt concentrations. The disjoining pressure and work show a roughly parabolic increase with the charge fraction, both in the flat and cylindrical case. In general, it can be noticed that the amount of salt plays only a modest role, in agreement with previous results.³⁸ Consequently, actuation is still possible at rather high salt concentrations.

Fig. 5 SCF results of the equilibrium strain, the disjoining pressure and the work performed per polymer chain in an expansion and contraction cycle, for the system with strongly charged polyelectrolytes, as a function of the fraction of charged groups on the polymer chain. For example, a charge fraction of 0.1 means that 1 in 10 segments carries a charge −e. In the first graph, the horizontal dashed line corresponds to the equilibrium system strain in the neutral case. Open symbols: cylindrical. Filled symbols: flat. Parameters: N = 200, σ = 1 nm⁻², χ = 0.4, lattice constant 0.3 nm.

In both flat and cylindrical configurations, the system strain δ_eq follows the same trend of brush height for single-grafted chains in the osmotic brush regime (see eqn (7)), indicated by the black dashed curve. Most of the parameter combinations considered in this study fall into the osmotic brush regime, except for very low values of f. The actuation pressures during expansion (ΔΠ > 0) have significantly large absolute values, due to the counterions, compared to the contraction (ΔΠ < 0), which is mainly driven by entropy (neutral brush going back to its equilibrium position). In the flat case, the pressure increase is noticeably higher than in the cylindrical one. Finally, the work depends similarly on the salt concentration but attains higher values in the cylindrical case. We believe that this is due to the fact that for the cylindrical configuration the integration is carried over a much larger volume than for the planar one.

The screening effect of the polyelectrolytes in the system only starts to affect the actuation performance to a small extent at 1 mol L⁻¹ salt concentration in the flat configuration, where δ_eq and the work per chain are noticeably lower. The salt concentration has a small influence on the disjoining pressure. A similar dependence on the salt concentration can be observed in the cylindrical configuration, although the effect of the ionic strength is larger: the system strain is considerably lower already at a salt concentration of 0.1 M. The greater effect of salt on the cylindrical case is probably due to the fact that, due to the lattice shape, the available volume increases towards the radial ends of the system, and thus the chains are relatively more diluted between the cylindrical layers at fixed grafting density.

3.2 Comparison with singly end-grafted weak polyelectrolytes

To justify our approach for weak polyelectrolytes described in Section 2.1.2, we compare the equilibrium strain δ in the flat configuration to the equilibrium thickness H of a singly end-grafted weak polyelectrolyte brush (Section 2.2) as derived using SCF by ref. 37.

We perform this comparison in the flat geometry using two different parameters: the equilibrium strain δ_eq and the root-mean-square brush thickness, H_rms. Here, H_rms is defined as

where z is the direction normal to the surface. The root-mean-square thickness was computed for a singly and a doubly-grafted brush. In the latter case, we computed the equilibrium strain δ_eq and then used the corresponding concentration profile to calculate H_rms in the equilibrium dimension of the system.

The results of our comparison are shown in Fig. 6, illustrating that the equilibrium strain δ_eq behaves similarly as H_rms for both singly and doubly-grafted brushes. This indicates that the loss in entropy due to the grafting at both ends of the brush does not significantly affect the brush equilibrium thickness. Moreover, it follows that δ_eq and H_rms depend on φ_s and pH in a similar fashion. In the doubly end-grafted case, we roughly find that δ_eq ≈ 1.7H_rms. Therefore, we conclude that the analytical models of Section 2 for singly grafted brushes can also be applied to our system.

Fig. 6 Comparison between the system equilibrium dimensions in our computations with weakly charged polyelectrolytes (a) and the equilibrium brush thickness for a singly end-grafted weak polyelectrolyte brush, reprinted with permission from ref. 37 (b). The left picture shows the equilibrium dimension of the system δ_eq in units of layers. The centre picture shows the H_rms for the singly end-grafted and doubly end-grafted weak polyelectrolyte brush. In this picture, the continuous lines and the dots represents the H_rms for singly end-grafted chains, the dashed lines and the squares represent the H_rms for doubly end-grafted chains. System parameters: N = 500, σ = 0.0056 nm⁻², χ = 0, lattice constant 0.6 nm, pK_a = 7.0.

3.3 pH-Responsive actuators: weak polyelectrolytes

In this section, results for weak polyelectrolytes are presented.

Results for the cylindrical and flat configurations are shown in Fig. 7. For each salt concentration, all the relevant quantities (equilibrium strain, disjoining pressures in expansion and contraction, and work) are computed for a switch between the current pH and the lowest pH considered at that salt concentration. In both the cylindrical and flat case δ_eq shows a sigmoid-like dependence on the pH of the system: at low pH, most of the carboxylic acid groups are protonated and the system is contracted; with increasing pH, the weakly acidic groups deprotonate and the system expands due to the repulsion between charges and the pressure of the counterions, which are still mostly inside the brush. A plateau is reached when the fraction of dissociated groups does not increase any more with increasing pH.

Fig. 7 SCF results of the system equilibrium strain, the disjoining pressure and the work performed per polymer chain in an expansion and contraction cycle, for the system with weakly charged polyelectrolytes, as a function of pH. In the strain graph, the horizontal dashed line indicates the equilibrium strain at the lowest pH possible at each salt concentration. Open symbols: cylindrical. Filled symbols: flat. Parameters: N = 200, σ = 1 nm⁻², χ = 0.4, lattice constant 0.3 nm, pK_a = 4.5.

The case of pH-responsive polyelectrolytes shows that the actuation pH range depends sensitively on the salt concentration c_s, shifting to higher pH values as c_s decreases. This is because at low salt concentrations, the dissociation is hindered by the proximity of other groups: electrostatic repulsion suppresses the (further) release of counterions from the polyelectrolytes. As a result, the pH actuation range is shifted to a higher pH.

Although the flat and cylindrical cases are qualitatively similar, there are some quantitative differences. In the flat configuration, for all salt concentrations, all three parameters show a plateau at approximately the same value. This is not the case for the cylindrical configuration, where especially the 0.1 M case differs significantly. We believe that this can be explained by the fact that in a cylindrical geometry the polymers are more dilute further from the cylinder surface; as a consequence, the screening effect of the salt is higher.

Furthermore, when comparing the results obtained with planar and cylindrical geometries in both strong and weak polyelectrolytes models, the flat configuration gives larger disjoining pressure than the cylindrical configuration for both the deprotonation ΔΠ_dep and protonation ΔΠ_pro steps. Especially ΔΠ_pro is markedly lower in the cylindrical case. However, the cylindrical configurations do have a larger stroke. We believe that the ‘radial dilution’ of the charges in the cylindrical systems leads to a lower actuation pressure than for the flat configuration. At the same time, as the polymer chains are more dilute, we think these give less ‘entropic resistance’ against swelling. In turn, these effects combined cause smaller internal pressure jumps ΔΠ and larger strokes in the cylindrical configuration. The work performed in the cylindrical case is almost double the planar one as the volume expansion scales in a quadratic way with the radial expansion, which compensates for the smaller pressures in the cylindrical system.

Our results for the brush extension (for singly-grafted polyelectrolytes) and equilibrium strain (for doubly-grafted polyelectrolytes) for the weakly acidic case are in line with previous experimental results with grafted weak polyelectrolytes.^39,55 For example, Parnell et al.³⁹ measured the brush height of poly(methacrylic acid) grafted to a substrate when changing pH of the solution in the range from 7 to 11. They observed an increase in the brush height of roughly 3–4 times, similar to our results in Fig. 7. Furthermore, their results indicate that high salt concentrations lead to a less pronounced swelling of the brush, as we also obtained (even though only for the cylindrical case) for the higher salt concentration 0.1 M in Fig. 7.

On the other hand, comparing the work per cycle that we obtain with experimental results is somewhat more difficult due to the lack of comparable experimental systems. Howse et al.,¹⁷ however, presented a prototype of a synthetic muscle based on block copolymers with hydrophobic end blocks of poly(methyl methacrylate) and a central block of poly(methacrylic acid), which swells and contracts upon changing pH. Each cycle in their system takes roughly half an hour, and the peak power is reported to be 20 mW kg⁻¹, translating into a peak power of roughly 5kT per hour per polymer chain, although it should be noted that the average power is much lower. Nevertheless, assuming that each cycle in our system also takes half an hour (performing 50kT of work per chain), a very rough estimation of the power of our system would be 100kT per hour per polymer chain. There are probably many factors contributing to this large difference, one for instance being that in our highly idealized system pH jumps are assumed to occur instantaneously, while in ref. 17 they occur on a much slower timescale. Nevertheless, the comparison allows our results to be placed in some experimental context.

Finally, both the fixed-charge ‘toy model’ and the pH-responsive model predict actuation properties of about the same order of magnitude, suggesting that the fixed charge approximation can serve as a useful first approximation to estimate the actuation properties of these systems.

4 Concluding remarks

We have presented a theoretical model for an actuator unit of muscle-like tissue, which consists of pH-sensitive polyelectrolyte doubly end-grafted between two plate- or rod-like colloidal particles. In a previous model,³⁸ the polyelectrolytes were modeled as chains with a fixed amount of charge. Here, a new model is introduced where the local degree of dissociation is also a variable and only the bulk pH is imposed as an external parameter to the system. This is also compared with a ‘toy model’ for strong polyelectrolytes, where the fraction of charged monomers along the chain is fixed. Our results show that it is possible to generate actuation using weakly acidic polymer chains. This could be used as a prototype for the development of actuating materials, for example, for prosthetic devices. The two cases of strong and weak polyelectrolytes show results for the functionality of the actuator of the same order of magnitude. The strong polyelectrolyte model with fixed charges could be used for a rough estimate of the actuation properties of a given system, provided the degree of dissociation is known a priori. On the other hand, the weak polyelectrolyte model can predict pH range where the actuation is possible, which, we show, depends strongly on salt concentration and the pK_a. We show that it is possible to tune the system parameters to attain the desired biocompatible pH range, or to get the desired actuation properties. For the weak polyelectrolyte case, we also checked that the equilibrium system strain for doubly-grafted polymer chains follows the same trends as the H_rms for a singly-grafted brush.

In the future, some limitations of the model could be further examined. First of all, the values of the disjoining pressure are considerably higher than common human muscles.⁵⁶ This is probably due to the fact that the lattice model employed for the calculations does not include any chemical detail and assumes that the sizes of all species (water, ions and monomer units) are the same, in this case 0.3 nm. This value was chosen to roughly reproduce the dimensions of a water molecule. However, a single unit of the polymer (whether it is a monomer or some other representative unit, for example based on the Kuhn length) can be much higher. Chain stiffness is not taken into account: the polyelectrolytes are modelled as fully flexible chains, where part of the segments can carry a negative charge. Furthermore, the flexibility is neither affected by the charge nor by the salt concentration, and we are not considering any dynamics in the system: we assume a sudden jump in pressure, followed by a reversible expansion (or contraction) of the system. Finally, note that we presented a model for a single actuating unit of the artificial muscle tissue. In order to study the full system, the alignment between the rods or platelets is of major importance, as well as the correlation between the motion of the different blocks, as the expansions and contractions should be coordinated. Still, our model does allow us to explore how changing parameters (such as salt concentration, pH, chain length) affect the actuation properties and can help to pave the way to the experimental realization of such actuators.

Author contributions

Eleonora Foschino: formal analysis, software, visualization, writing – original draft. Irene E. Hulsen: formal analysis, software, writing – review and editing. Alessandro Ianiro: conceptualization, methodology, funding acquisition, writing – review and editing. Remco Tuinier: conceptualization, funding acquisition, methodology, supervision, project administration, writing – review and editing. Mark Vis: conceptualization, supervision, methodology, project administration, writing – review and editing.

Conflicts of interest

There are no conflicts to declare.

Data availability

Data for this article, including data files and analysis scripts, are available at 4TU.ResearchData at https://doi.org/10.4121/be927bb6-5c0b-4abc-a5e9-410fde1df7cf.

Appendices

Appendix A

To perform SCF lattice computations, the equilibrium constants K_a for acid–base reactions (as well as the pK_a values) should be expressed in terms of dimensionless (lattice) volume fractions, and not in terms of molar concentrations. The equilibrium constant K^M_a based on molar concentrations is related to the ‘lattice’ equilibrium constant based on volume fractions K^latt_a via:
In this expression, the constant c is given by 10³b³N_Av, where N_Av is Avogadro's number, b the lattice size in m³ and the factor 1 × 10³ is necessary to convert from mol L⁻¹ to mol m⁻³. Using b = 3 × 10⁻¹⁰ m, c has a value of approximately 0.016294. For the pK_a, the value of pK^latt_a can be calculated as:
In a similar way, the equilibrium constant of the autoionization of water K^M_w (given by K^M_w = [H₃O⁺][OH⁻]) is converted in ‘lattice’ equilibrium constant K^latt_w by:

K^latt_w = φ_H3O⁺φ_OH⁻ = c[H₃O⁺]c[OH⁻] = c²[H₃O⁺][OH⁻] = c²K^M_w.

pK^latt_w then becomes:

Appendix B

The value of the pK_a of the weakly acidic chains that we mentioned in this work is the one of dilute chains. However, if the monomer density increases, we expect the pK_a to shift to higher values because of the presence of other acidic groups along the chain, which oppose dissociation.

Fig. 8 shows the degree of dissociation of a brush grafted to a flat surface as a function of varying grafting density σ. Every other monomer A of the polymer chain can dissociate according to the reaction HA + H₂O ⇌ H₃O⁺ + A⁻, with a pK_a of 4.5. Therefore, we have the same system as in the planar configuration case presented in Section 3.3, but using singly end-grafted chains. The degree of polymerization is kept constant at N_p = 200 and the salt concentration is c_s = 0.001 M, a relatively low value so that the effect of the increasing grafting density is not ‘hidden’ by a large Debye screening constant. As can be seen from the graph, the ‘effective pK_a’ approximately corresponds to the dilute pK_a in the dilute case (low grafting density), but it shifts to higher values when the grafting density, and therefore the monomer density in the brush, is increased.

Fig. 8 Dissociation degree as a function of the grafting density of the polyelectrolyte chains. N_p = 200, c_s = 0.001 M, χ = 0.4.

Acknowledgements

This work was supported by the EIC Pathfinder Open “INTEGRATE” Grant Agreement No. 101046333. The authors thank all members of the “INTEGRATE” consortium for the comments and discussions about the work presented in this paper. The authors thank professor Frans Leermakers for the advice and fruitful discussions and for providing the code for the SCF computations, and professor Jens-Uwe Sommer for useful suggestions and discussions. E. F. thanks A. F. den Ouden for the advice and useful comments.

Notes and references

1. C. Greco, P. Kotak, L. Pagnotta and C. Lamuta, Int. Mater. Rev., 2022, 67, 575–619.
2. M. Li, A. Pal, A. Aghakhani, A. Pena-Francesch and M. Sitti, Nat. Rev. Mater., 2021, 7, 235–249.
3. I. Apsite, S. Salehi and L. Ionov, Chem. Rev., 2022, 122, 1349–1415.
4. L. M. de Espinosa, W. Meesorn, D. Moatsou and C. Weder, Chem. Rev., 2017, 117, 12851–12892.
5. J. Chen, F. K.-C. Leung, M. C. A. Stuart, T. Kajitani, T. Fukushima, E. V. D. Giessen and B. L. Feringa, Nat. Chem., 2018, 10, 132–138.
6. T. Nitta, Y. Wang, Z. Du, K. Morishima and Y. Hiratsuka, Nat. Mater., 2021, 20, 1149–1155.
7. P. Ariano, D. Accardo, M. Lombardi, S. Bocchini, L. Draghi, L. D. Nardo and P. Fino, J. Appl. Biomater. Funct. Mater., 2015, 13, 1–9.
8. F. Carpi, R. Kornbluh, P. Sommer-Larsen and G. Alici, Bioinspiration Biomimetics, 2011, 6, 045006.
9. S. Li, H. Bai, Z. Liu, X. Zhang, C. Huang, L. W. Wiesner, M. Silberstein and R. F. Shepherd, Sci. Adv., 2021, 7, eabg3677.
10. Liquid Crystal Elastomers: Materials and Applications, ed. W. H. de Jeu, Springer Verlag, 2012.
11. G. Agarwal, N. Besuchet, B. Audergon and J. Paik, Sci. Rep., 2016, 6, 34224.
12. W. B. Han, G.-J. Ko, K.-G. Lee, D. Kim, J. H. Lee, S. M. Yang, D.-J. Kim, J.-W. Shin, T.-M. Jang, S. Han, H. Zhou, H. Kang, J. H. Lim, K. Rajaram, H. Cheng, Y.-D. Park, S. H. Kim and S.-W. Hwang, Nat. Commun., 2023, 14, 2263.
13. M. C. Koetting, J. T. Peters, S. D. Steichen and N. A. Peppas, Mater. Sci. Eng., R, 2015, 93, 1–49.
14. S. Maeda, Y. Hara, T. Sakai, R. Yoshida and S. Hashimoto, Adv. Mater., 2007, 19, 3480–3484.
15. R. Yoshida, Biophysics, 2012, 8, 163–172.
16. P. D. Topham, J. R. Howse, C. J. Crook, S. P. Armes, R. A. L. Jones and A. J. Ryan, Macromolecules, 2007, 40, 4393–4395.
17. J. R. Howse, P. Topham, C. J. Crook, A. J. Gleeson, W. Bras, R. A. L. Jones and A. J. Ryan, Nano Lett., 2006, 6, 73–77.
18. G. W. de Groot, M. G. Santonicola, K. Sugihara, T. Zambelli, E. Reimhult, J. Vörös and G. J. Vancso, ACS Appl. Mater. Interfaces, 2013, 5, 1400–1407.
19. J. Lahann, S. Mitragotri, T.-N. Tran, H. Kaido, J. Sundaram, I. S. Choi, S. Hoffer, G. A. Somorjai and R. Langer, Science, 2003, 299, 371–374.
20. S. A. Prokhorova, A. Kopyshev, A. Ramakrishnan, H. Zhang and J. Rühe, Nanotechnology, 2003, 14, 1098.
21. S. Santer, A. Kopyshev, J. Donges, H.-K. Yang and J. Rühe, Adv. Mater., 2006, 18, 2359–2362.
22. S. Santer and J. Rühe, Polymer, 2004, 45, 8279–8297.
23. M. Geoghegan, Soft Matter, 2022, 18, 2500–2511.
24. W.-X. Hong, V. Y. Shevtsov and Y.-T. Shieh, J. Polym. Res., 2022, 29, 345.
25. C. Endes, S. Camarero-Espinosa, S. Mueller, E. J. Foster, A. Petri-Fink, B. Rothen-Rutishauser, C. Weder and M. J. D. Clift, J. Nanobiotechnol., 2016, 14, 78.
26. J. C. Natterodt, A. Petri-Fink, C. Weder and J. O. Zoppe, Chimia, 2017, 71, 376.
27. I. A. Sacui, R. C. Nieuwendaal, D. J. Burnett, S. J. Stranick, M. Jorfi, C. Weder, E. J. Foster, R. T. Olsson and J. W. Gilman, ACS Appl. Mater. Interfaces, 2014, 6, 6127–6138.
28. B. Sun, G. Zhou and H. Zhang, Prog. Solid State Chem., 2016, 44, 1–19.
29. A. Kuijk, D. V. Byelov, A. V. Petukhov, A. van Blaaderen and A. Imhof, Faraday Discuss., 2012, 159, 181–199.
30. N. Yaghmaeiyan, M. Mirzaei and R. Delghavi, Results Chem., 2022, 4, 100549.
31. S. Morimune-Moriya, M. Kotera and T. Nishino, Polymer, 2022, 256, 125202.
32. J. Chanra, E. Budianto and B. Soegijono, IOP Conf. Ser.: Mater. Sci. Eng., 2019, 509, 012057.
33. P. Pincus, Macromolecules, 1991, 24, 2912–2919.
34. O. V. Borisov, T. M. Birshtein and E. B. Zhulina, J. Phys. II, 1991, 1, 521–526.
35. E. B. Zhulina, O. V. Borisov and T. M. Birshtein, J. Phys. II, 1992, 2, 63–74.
36. R. Israëls, F. A. M. Leermakers, G. J. Fleer and E. B. Zhulina, Macromolecules, 1994, 27, 3249–3261.
37. R. Israëls, F. A. M. Leermakers and G. J. Fleer, Macromolecules, 1994, 27, 3087–3093.
38. A. Ianiro, J. A. Berrocal, R. Tuinier, M. Mayer and C. Weder, J. Chem. Phys., 2023, 158, 14901.
39. A. J. Parnell, S. J. Martin, C. C. Dang, M. Geoghegan, R. A. Jones, C. J. Crook, J. R. Howse and A. J. Ryan, Polymer, 2009, 50, 1005–1014.
40. J. M. H. M. Scheutjens and G. J. Fleer, J. Phys. Chem., 1978, 83, 1619–1635.
41. F. A. M. Leermakers and J. M. H. M. Scheutjens, J. Phys. Chem., 1989, 93, 7417–7426.
42. G. J. Fleer, Adv. Colloid Interface Sci., 2010, 159, 99–116.
43. C. M. Wijmans and E. B. Zhulina, Macromolecules, 1993, 26, 7214–7224.
44. C. M. Wijmans, J. M. H. M. Scheutjens and E. B. Zhulina, Macromolecules, 1992, 25, 2657–2665.
45. G. J. Fleer, M. A. Cohen Stuart, J. M. H. M. Scheutjens, T. Cosgrove and B. Vincent, Polymers at Interfaces, Springer, Dordrecht, 1993, p. 496.
46. J. van Male, Self-consistent-field theory for chain molecules: extensions, computational aspects, and applications, Wageningen Universiteit, Wageningen, 2003.
47. F. A. M. Leermakers, R. Varadharajan, D. Emmery and A. Kazakov, Namics – MD-SCF Hybrid simulation tool, 2023, https://github.com/leermakers/Namics.
48. L. F. W. Vleugels, S. Ricois, I. K. Voets and R. Tuinier, Food Hydrocolloids, 2018, 81, 273–283.
49. P. Linse and M. Bjoerling, Macromolecules, 1991, 24, 6700–6711.
50. E. B. Zhulina and O. V. Borisov, Langmuir, 2011, 27, 10615–10633.
51. A. A. Mercurieva, T. M. Birshtein, E. B. Zhulina, J. V. M. P. Iakovlev and F. A. M. Leermakers, Macromolecules, 2002, 35, 4739–4752.
52. O. Rud, T. Richter, O. Borisov, C. Holm and P. Košovan, Soft Matter, 2017, 13, 3264–3274.
53. S. Alexander, J. Phys., 1977, 38, 983–987.
54. A. Naji, C. Seidel and R. R. Netz, Surface-Initiated Polymerization II, Springer-Verlag, 2006, vol. 198, pp. 149–183.
55. M. Geoghegan, L. Ruiz-Pérez, C. C. Dang, A. J. Parnell, S. J. Martin, J. R. Howse, R. A. L. Jones, R. Golestanian, P. D. Topham, C. J. Crook, A. J. Ryan, D. S. Sivia, J. R. P. Webster and A. Menelle, Soft Matter, 2006, 2, 1076–1080.
56. R. K. Josephson, Annu. Rev. Physiol., 1993, 55, 527–546.

---

Footnote

† Present address: Laboratory of Self-Organizing Soft Matter, Department of Chemical Engineering and Chemistry & Institute for Complex Molecular Systems (ICMS), Eindhoven University of Technology, P.O. Box 513, 5600 MB, Eindhoven, The Netherlands.

---