## Abstract

Nanoparticle (NP) size and charge play key roles in bioconjugation chemistry, imaging and drug delivery. Although the electrophoretic mobility and hydrodynamic size are routinely measured, interpreting these data can be extremely difficult. Here, the challenge is addressed via an electrokinetic model for spheres bearing a soft amphoteric corona, the charge of which is regulated by a multi-component electrolyte. The model is applied to NPs with a metallic core to which are grafted poly(ethylene glycol) chains with either weak acid or amphiprotic end groups. The results elucidate the separate roles of electrolyte pH and ionic strength on the electrophoretic mobility and diffusion coefficient. In this study, the forces were evaluated directly, rather than from the Stokeslet velocity disturbances. While the second-order convergence was demonstrated by both methods, the direct approach, which uses only the inner part of the global solution, furnished superior accuracy and robustness. This may benefit future attempts to model the dielectric and electroacoustic properties of these complex nanoparticulates.

## 1. Introduction

Free-solution electrophoresis and gel electrophoresis are widely adopted to characterize nanoparticles (NPs) [1–5]. Despite their ubiquity, the electrophoretic mobility and its accompanying *ζ*-potential are challenging to interpret because they reflect the distribution of polymer segments and charge, which may be coupled to the composition of the electrolyte in which the particles are dispersed [6].

The challenges are exemplified by noting that the mobility of an impenetrable sphere with a prescribed surface charge density (no corona) depends on the scaled *ζ*-potential *ζe*/(*k*_{B}*T*) and the scaled radius *κa* (*k*_{B}*T*/*e* is the ratio of the thermal energy to the fundamental charge, and *κ* is the reciprocal Debye length). When *κa*∼1, which is often the case for NPs (as compared with their larger counterparts with size

Adding an inhomogeneous, charge-regulating polyelectrolyte coating significantly compounds the challenges, motivating recent initiatives to engage molecular-modelling methodologies, e.g. to interpret mobility reversal of polyelectrolyte-grafted NPs [9]. Despite an obvious potential for continuum breakdown at such small scales, continuum electrokinetic models—such as the one advanced in this paper—have demonstrated remarkable consistency with their discrete-particle/molecular counterparts [9,10]. In addition to providing physical insight, the computational expediency of continuum models furnishes powerful quantitative interpretations of experiments, e.g. NP gel electrophoresis [11], dielectric spectroscopy of colloidal phase transitions [12], and the electrokinetic sonic amplitude of NP-doped hydrogels [13].

This work is developed using—as a guiding example—the general features of the prototypical functional NP to which bio-conjugation is undertaken [14], namely an impenetrable spherical core (e.g. Au) that bears surface charge (e.g. from adsorbed electrolyte ions, such as Cl^{−}) and a soft corona of grafted polymer (e.g. PEG) with peripheral end groups (e.g. a carboxymethyl or azido moiety). Note that the model is readily adapted (by varying its scalar parameters) to a variety of soft spherical colloidal particulates, including microgels (vanishing core) and biological cells. These have unique electrokinetic signatures when the electrophoretic mobility is correlated against electrolyte composition.

This study advances the MPEK software package [15], which has been available (by request, without charge) to the research community for a decade. Several enhancements to expedite its application to dispersions with multi-component electrolytes and particles with non-uniform, charge-regulating polymer coatings are addressed. These unite other initiatives in the literature to address charge-regulating surfaces and polyelectrolyte coatings, and will pave the way to better models of dielectric, electroacoustic and radio-frequency-heating properties [16].

## 2. Theory

The fundamental equations—comprising the Poisson, ion-conservation, and fluid momentum and mass conservation equations—to describe the steady-state dynamics for a particle that may be translating with velocity ** V** are

*ψ*is the electrostatic potential,

**and**

*u**p*are the fluid velocity and pressure, and

*n*

_{j}are the mobile-ion concentrations with each ion having its distinct charge

*z*

_{j}

*e*and mobility

*D*

_{j}/(

*k*

_{B}

*T*), where

*D*

_{j}is a diffusion coefficient,

*k*

_{B}

*T*is the thermal energy and

*e*is the fundamental charge. Furthermore,

*ρ*

_{m}and

*ρ*

_{f}denote the mobile (in the electrolyte) and immobile (on the skeleton) charge densities. Physical properties include the solvent/electrolyte and particle dielectric constants,

*ϵ*

_{0}

*ϵ*

_{s}and

*ϵ*

_{0}

*ϵ*, and the solvent/electrolyte viscosity,

*η*. The Brinkman screening length ℓ is related to a number density of Stokes resistance centres

*n*

_{s}, each of which exerts a drag force on the fluid −6

*πηF*

_{s}

*a*

_{s}(

**−**

*u***), where**

*v**F*

_{s}

*a*

_{s}is the hydrodynamic radius of a centre with physical radius

*a*

_{s}, and

*F*

_{s}is a dimensionless function of the resistance-centre volume fraction

^{−2}(

*r*)=6

*πF*

_{s}(

*r*)

*a*

_{s}

*n*

_{s}(

*r*).

### (a) Equilibrium solution

The equilibrium state, identified with the superscript 0, corresponds to a stationary particle, corona and electrolyte in the absence of any external forcing. The equations above reduce to:

### (b) Perturbations

With the application of an external driving force, which is taken to be a vector directed along the *z*-axis, all the variables are expressed as the sum of the foregoing equilibrium solutions and their perturbations. As the model will be truncated to linear order in these perturbations, linearity and symmetry considerations demand solutions that have the following forms:
** X**=

*X*

*e*_{z}denotes the vector (directed along the

*z*-axis) that drives the perturbations: either the applied electric field vector

**=**

*E**E*

*e*_{z}or particle velocity

**=**

*V**V*

*e*_{z}. Note that

**=0 is the far-field pressure gradient, and**

*P***=**

*U**U*

*e*_{z}is the far-field fluid velocity; moreover,

**=**

*r**r*

*e*_{r}is the position vector with

*e*_{r}the radial unit vector.

In the linear approximation, the perturbations satisfy

Taking the curl of the momentum equation (to eliminate the pressure) and writing the perturbation equations in terms of the radial perturbation functions (hatted variables) furnishes the linear ordinary differential system
*r* denote radial differentiation, and

### (c) Charge-regulation closure

The charge-regulation model here reduces to the single dissociation reaction in the original MPEK package. It is similar to the two-site (acid–base) model adopted by Yeh *et al.* [19] for spherical polyelectrolytes and is the same as the model adopted by Chen *et al.* [20] for the surface charge of pH-regulated spheres. Here, however, the model will be perturbed from the equilibrium state, thereby capturing the secondary charge-density perturbations that arise from the disturbances induced by the external forcing (diffusion, electro-osmosis or electrophoresis).

The immobile charge associated with the *j*th mobile ion species results from dissociation/association equilibria of the form
_{2} denote complexes with one and two ions of species Y. Moreover, subscripts *j* on symbols X, Y, *n**_{f}, *K*_{1} and *K*_{2} have been temporarily discarded for notational convenience, and the concentrations (mol l^{−1}) are denoted using square brackets. From this model, the immobile charge density associated with the *j*th mobile ion can be written as
*j* are now affixed so that the variables are consistent with the conservation equations above, i.e.
*n**_{f,j}, *K*_{1,j} and *K*_{2,j} (and bulk electrolyte composition) are prescribed. In this work, *n**_{f,1} (for the amphiprotic moiety associated with the H^{+} ion) is a radial function of position with the other *n**_{f,j}=0.

### (d) Force evaluation

The forces have customarily been evaluated from the strength of the far-field Stokeslet velocity disturbances in the so-called *U*- and *E*-problems of [7] for the bare-sphere electrophoresis problem. However, there are several benefits to be gained by directly evaluating the forces. Firstly, the total force can be broken down into separate contributions, each of which has a clear physical origin; and, secondly, a direct evaluation will be demonstrated to exhibit much improved computational accuracy.

Evaluating the Stokeslet velocity disturbance requires the correct asymptotic decay of all the perturbed variables to be captured at the boundary of a finite computational domain. Authors who have adopted finite-element packages to solve similar problems have invariably set these disturbances to zero at a large, but finite, distance from the particle contained within and none have quantified the errors of such an approximation.

The following derivation of the force is inspired by López-Garcia *et al.* [21], who used it to ascertain the forces on soft spheres with uniform coatings. Here it is generalized for the non-uniform layers. Accordingly, the force on the particle with polymer coating (assumed to be rigid) is
** X**, furnishing

*r*) confirm that

*P*=0 because ℓ

^{−2}and

*r*=

*a*

_{c}. Thus, substituting the pressure disturbance from equation (2.1) gives

*r*=

*a*

_{c}has been applied,

*Ω*and its radial derivative

*Ω*

_{r}. In this work, the permeability ℓ

^{2}is not a constant, so the radial integrals must be evaluated numerically. This is accomplished on the non-uniform finite-difference mesh using the trapezoidal rule. These discretization errors were negligible compared with those from the second-order discretization of the differential equations. Note that the integrands decay rapidly, because of the rapid decays of ℓ

^{−2}(

*r*),

*ψ*

^{0}(

*r*) beyond the polymer coating.

### (e) Electrophoretic mobility

The electrophoretic mobility *V*/*E* comes from the particle equation of motion

In practice, the equations for the *V* -problem above—in which the particle translates at velocity ** V**=

*V*

*e*_{z}in a quiescent electrolyte with

*E*=0—are solved as a

*U*-problem, i.e. with

*V*=

*E*=0 and

*F*^{V}=−

*F*^{U=V}, so the electrophoretic mobility

*F*

^{E}and

*F*

^{U=V}denote the (signed) magnitude of these forces when

*E*=1 and

*U*=1, respectively.

## 3. Results

The electrokinetic model is now applied to explore how the ionic strength and pH of the bulk electrolyte impact the electrophoretic mobility of Au NPs bearing terminally anchored 5 kDa poly(ethylene glycol) (PEG) ligands with either weak acid or amphiprotic end groups. In both cases, any coupling of the corona structure to its electrostatic state is completely neglected. Here, the number densities of the PEG segments (Stokes resistance centres) and end groups are specified using the Gaussian functions
*N*_{a} is the aggregation number for PEG ligands with *N*_{s} the number of segments in a chain. Note that the number of end groups *N*_{c} is equal to *N*_{a}, because each grafted PEG chain bears one carboxymethyl (CM) (or other) end group. For convenience, *N*_{s} was set so that the product *n*_{s}(*r*)6*πa*_{s}=ℓ^{−2}(*r*) furnished an acceptable corona permeability (and thus NP hydrodynamic radius) with *F*_{s}≈1.

With a prescribed aggregation number (e.g. available from thermogravimetric analysis), there are three geometrical model parameters: (i) the radial extent of the coating, *δ*_{1}; (ii) the nominal radial position of the end groups, *L*_{2}∼*δ*_{1}; and (iii) the nominal width of the peripheral layer of end groups, *δ*_{2}. Ideally, these geometrical parameters should be prescribed using independent knowledge of the PEG-chain size and conformation, or from the hydrodynamic size, which may be available from diffusion experiments, such as dynamic light scattering, fluorescence correlation spectroscopy or sedimentation. In this study, they are set to physically acceptable constants that also furnished reasonable predictions and interpretations of experimental data [22].

The model parameters adopted for the NP core, corona and electrolyte are summarized in tables 1 and 2. Note that there are four electrolyte ions: the concentrations of H^{+} and OH^{−} are set according to the prescribed pH, with the concentration of Na^{+} chosen to vary the bulk ionic strength, and the concentration of Cl^{−} calculated to ensure an electroneutral bulk electrolyte. This corresponds to mixing NaCl, NaOH and HCl to form an electrolyte with a prescribed pH and ionic strength or conductivity. There is clearly a lower limit to the ionic strength at each pH, and this increases as the bulk pH deviates from 7.

### (a) Weak-acid/dissociation-only charge regulation

Figures 1–4 present results for the CM end groups when they are modelled as weak acids with pK_{1,1}=4.85 and *ζ*-potential in mV, i.e. *ζ*=*Mη*/(*ϵ*_{0}*ϵ*_{s})—is plotted versus the ionic strength in figure 1 for several pH values in the range 4.5–7. Elsewhere, the mobility is reported as the dimensionless quantity (Hückel scaling)
*ζ*∼−40 mV), and (ii) PEGylated Au NPs *without charged end groups* are dispersed in similar electrolytes (*ζ*∼−5 mV) [11].

The theory captures two important qualitative features of the experimental data: firstly that mobilities are highest at low ionic strength, when much of the corona countercharge—and possibly some of the core countercharge—resides beyond the corona (in the free electrolyte); and secondly that mobilities approach those of particles with zero charge on the corona as the corona charge vanishes (at low pH). The theory is consistent with experiments only to the extent that the mobilities at higher ionic strengths correspond well at the pH values of the TAE and PBS buffers used in the experiments (6.1 for RO water to 8.2 for 1×TAE buffer). Thus, the theory does not explain the lower mobilities measured at lower ionic strengths, unless the pH were to decrease with decreasing ionic strength.

It is plausible that the buffering capacity of the electrolytes would be exhausted at low ionic strength, so that the pH decreases in response to dissolved CO_{2} and the CM counterions. However, this would require the mobilities measured in RO water to be much lower, perhaps corresponding to −*ζ*∼20–30 mV rather than the measured *ζ*∼−40 mV. To address these uncertainties, future experimental studies should report dispersion conductivities *and* pH, and correlate these with the concentrations of the NPs and added electrolytes, particularly at low ionic strengths.

The mobilities in figure 1 are replotted in figure 2 versus the electrolyte conductivity. Here, the solid lines are based on the conductivity of the bulk electrolyte at infinite NP dilution. The dashed lines, however, have corrections from electrical and ion-concentration polarization, and the green/dashed-dotted lines include additional terms from added counterions and non-specific adsorption [12,23]. These corrections are proportional to the NP concentration, so the calculations were undertaken with an NP volume fraction (based on the core) *ϕ*_{p}=10^{−4}, which corresponds to a molar concentration ∼2 μmol l^{−1}, a characteristic particle separation ∼90 nm, or a Au core mass concentration ∼2 *g* *l*^{−1}. As demonstrated in figure 2, a small but finite NP volume fraction can have a significant impact on the dispersion conductivity when the particles are dispersed at low ionic strength/conductivity [8].

Figure 3 shows how the hydrodynamic size (typically measured using light scattering) depends on the ionic strength. These data are furnished by the *V* -problem that contributes to the mobility calculations reported in figures 1 and 2. Note that the hydrodynamic size is reported as a dimensionless drag coefficient *F*, which is equivalent to the hydrodynamic radius *a*_{h} divided by the NP core radius *a*_{c}, and the ionic strength is reported as the scaled reciprocal Debye length, *κa*_{c}. Similar to the mobilities in figure 1, the lower limits of the curves correspond to the lowest amount of added salt (1 μmol l^{−1} at the corresponding pH).

The limiting value of the hydrodynamic radius at high ionic strength (*κa*_{c}≫1) is the value corresponding to an uncharged corona. This shows that the corona has a hydrodynamic thickness ≈(3.2−1)*a*_{c}=2.7*a*_{c}≈6 nm, so the apparent NP radius (core plus corona) is ≈9 nm. At lower ionic strengths, when *κa*_{c}∼0.1, the hydrodynamic size is somewhat larger because of electroviscous effects. These arise from polarization of the electrical charge, which imparts a flow-induced electrical force that enhances the apparent hydrodynamic drag. This clearly increases with the increasing corona charge that accompanies an increase in the electrolyte pH.

Note that the changing apparent hydrodynamic size in figure 3 does not reflect a change in the structure or dimensions of the corona (table 1). Hill *et al.* [22] recently demonstrated that a corona thickness that decreases with increasing electrolyte ionic strength, albeit with a fixed corona charge, is necessary to quantitatively interpret the experimental data in figure 1. Thus, the calculations undertaken here (with charge regulation, but fixed corona dimensions) support their hypothesis that electrostatic interactions—among the CM end groups and, perhaps, involving the core surface charge—modulate corona structure. Such interactions were demonstrated to be particularly effective at very low ionic strengths, when the Debye length becomes comparable to the unperturbed layer thickness. At these ionic strengths, there is a delicate balance of hydrodynamic drag, charge polarization and electro-osmosis (all of which are enhanced by the corona expansion), so the mobilities are particularly sensitive to the corona structure.

Representative streamlines of the flows (in the laboratory frame) for the *V* - and *E*-problems, and for electrophoresis, are shown in figure 4. The *V* -problem demonstrates how the corona entrains fluid, imparting a hydrodynamic radius *a*_{h}≈3.2*a*_{c} that is modestly greater than the characteristic radius *a*_{c}+*δ*_{1}≈2.48*a*_{c}. The *E*-problem identifies electrolyte being electro-osmotically pumped through the charged corona periphery, and the electrophoresis problem demonstrates that this electro-osmotic pumping combined with the particle translation yields a toroidal recirculation, where the centre is located close to the particle. For a negatively charged corona, the electro-osmotic flow in the *E*-problem is directed from left to right in figure 4*b* (with *E*>0), and the particle undergoing electrophoresis in figure 4*c* is translating from right to left, so the recirculating flow in the top (bottom) half of figure 4*c* is clockwise (counterclockwise).

### (b) Amphiprotic charge regulation

The results above were undertaken with the assumption that the CM end groups are weak acids. Here, the end groups are amphiprotic, so that they may also attain a net positive charge by binding two H^{+} ions. The dissociation constants *K*_{1,1} and *K*_{1,2} are set equal (to 4.85), as indicated in table 2, thus producing an isoelectric point at pH=4.85. Note that, because the corona is not uniform, the local pH is a non-monotonic function of radial position. Nevertheless, the electrophoretic mobilities shown in figure 5 indicate that the NPs behave as if they have a point of zero charge when the bulk electrolyte pH≈4.85. As expected, the results with pH=7 are the same as in figure 1, because at this pH all the groups are fully dissociated, bearing the maximum possible negative charge. Interestingly, the theoretical calculations with pH=6 superficially capture the experimental trends in figure 1, almost bridging the transition from buffer-containing electrolytes to ‘pure’ water. Note, however, that the CM moieties (experimental data in figure 1) are not amphiprotic.

Representative streamlines of flows (in the laboratory frame) for the *V* - and *E*-problems, and for electrophoresis, are shown in figure 6. While the streamlines for the *V* - and *E*-problems are similar to those in figure 4, the streamlines for electrophoresis now reveal two counter-rotating toroidal flows. Here, the particle should be considered to translate from right to left (the mobility is negative) when the electric field is directed from left to right, while the electro-osmotic flow in the corona periphery is from left to right (the countercharge is positive). Clearly, the strength of the electro-osmotic flow in the *E*-problem can have a profound influence on the fluid dynamics during electrophoresis. This is exemplified by the sequence of electrophoresis streamlines shown in figure 7 at several successively increasing ionic strengths, all with the same bulk pH=4.5. The particles in this figure should be considered to translate from left to right when the electric field is directed from left to right, as the mobilities are positive at the prevailing pH. At low ionic strength, the counterions reside mostly outside the corona, so the electro-osmotic flow is very weak. Here, the particle motion reflects the net bare Coulomb electrical force balancing the hydrodynamic drag. Increasing the ionic strength brings the counterions into the corona periphery, strengthening the role of the electo-osmotic flow, and increasing the viscous dissipation. At the highest ionic strength (figure 7*d*), the electro-osmotic flow in the corona periphery is from right to left (negative countercharge), and so the outermost circulation is in the opposite direction to those at lower ionic strengths (figure 7*a*–*c*), even though the particle still translates in the same direction; this is possible because of the two counter-rotating vortexes.

### (c) Computational performance

For a given adaptive-mesh algorithm, the computational accuracy depends on (i) the prescribed number of grid points *N* and (ii) the radial extent of the computational domain *R* (located at *r*=*a*_{c}+*R*). Note that the far-field boundary conditions at *r*=*a*_{c}+*R* permit every perturbation variable to realize its far-field asymptotic decay [23]. Such behaviour occurs in the electrically neutral fluid, beyond the corona and diffuse charge layer, i.e. where *r*=*a*_{c}+*R*). Here, we examine the computational merits of directly evaluating the forces via the surface and volume integrations in equation (2.2).

Figure 8*a* shows the drag coefficient *F* (*a*(i)) and dimensionless mobility *M** (*a*(ii)) versus *κR* with *N*=6400, which is smaller than the *N*=8000 grid points used for the production results above, but still large enough to achieve high accuracy with the largest values of *κR*. The figures compare the Stokeslet (open symbols) and direct (filled symbols) force evaluations. These clearly show that directly evaluating the forces using equation (2.2) achieves superior accuracy, thus affirming that the inner solution is somewhat tolerant of errors in the outer region. Figure 8*b* shows the relative errors in *F* (*b*(i)) and *M** (*b*(ii)) when increasing *N* with fixed *κR*. Here, the ordinates |*ϵ*| are the absolute value of the error relative to the value furnished by Richardson extrapolation of the data with the two largest values of *N*, using an interpolation of the form *f* (either *F* or *M**) when *c* is a constant.

As expected from the second-order finite-difference approximations used to discretize the model equations, both force-evaluation methods exhibit approximately the second-order convergence. The direct-force evaluation, however, produces notably weaker fluctuations—particularly in the drag force. Note that the Richardson extrapolations of both force-evaluation methods furnish the same values of *F* and *M** to five and four significant figures (*F*≈3.2823 and *M**≈−2.212, dashed lines in figure 8*a*(i,ii)), indicating that the radial integrations required to evaluate equation (2.2) contribute a negligible error. Similar results (shown in appendix C) were achieved with a much thinner diffuse layer, *κa*_{c}≈3.96.

All the production calculations presented above were undertaken with *N*=8000 grid points and *κR*=200(1+*κa*+*κδ*_{1}). Note that *R* is considerably larger than recommended in the literature to solve comparable problems using finite-element software: these have generally required only that the decaying perturbed variables vanish on an outer boundary. While the errors will surely be larger than those achieved in this study, the relative magnitudes may be—on account of the analysis in figure 8—approximately less than 10%, and thus imperceptible to the eye. Notable exceptions will be addressed elsewhere.

## 4. Conclusion

An electrokinetic model for a wide variety of amphoteric charge-regulating soft spheres with complex core–shell structure has been developed. The computational solution, which efficiently addresses the inherent stiffness arising from several disparate length scales, accounts for important electrokinetic processes, such as nonlinear electrostatics, charge regulation, and ion-concentration and charge perturbations.

The model superficially captured trends in the experimental data of Hill *et al.* [22] for soft PEGylated NPs bearing a corona with CM (weak acid) end groups. However, the mobility-ionic strength relationship for these NPs seems more likely to reflect a CM-PEG layer that changes its structure/extent with the bulk electrolyte ionic strength, as an acid-dissociation-only model did not adequately capture all the available experimental trends over the full range of ionic strengths. Therefore, this study of charge regulation supports the recent interpretation of Hill *et al.* [22], whereby the CM groups were taken to be fully dissociated and the corona structure was varied with the electrolyte ionic strength.

Particular attention was given to the computational accuracy, and the merits of directly evaluating the forces. The second-order convergence was demonstrated, and a direct force evaluation was found to significantly improve the computational accuracy and robustness, as compared with the customary Stokeslet analysis. Interestingly, a computational domain extending several hundred particle radii from the particle surface was required to achieve four-significant-figure accuracy in the electrophoretic mobilities and drag coefficient, even with an adaptive grid and boundary conditions that permit all the perturbation variables to decay according to their respective asymptotic decay rates.

## Data accessibility

This work has no accompanying data.

## Authors' contributions

The author performed all the work.

## Competing interests

The author has no competing interests.

## Funding

This work was supported by an NSERC Discovery Grant.

## Acknowledgements

The author thanks two anonymous referees for carefully reading the manuscript and offering helpful suggestions.

## Appendix A. Vorticity

The vorticity is

## Appendix B. Differential identities

The following identities are used to derive results in the main text (** I** is the identity tensor):

- Received July 29, 2015.
- Accepted October 27, 2015.

- © 2015 The Author(s)