---
abstract: 'When a highly charged globular macromolecule, such as a dendritic polyelectrolyte or charged nanogel, is immersed into a physiological electrolyte solution, monovalent and divalent counterions from the solution bind to the macromolecule in a certain ratio and thereby almost completely electroneutralize it. For charged macromolecules in biological media, the number ratio of bound mono- versus divalent ions is decisive for the desired function. A theoretical prediction of such a sorption ratio is challenging because of the competition of electrostatic (valency), ion-specific, and binding saturation effects. Here, we devise and discuss a few approximate models to predict such an equilibrium sorption ratio by extending and combining established electrostatic binding theories such as Donnan, Langmuir, Manning as well as Poisson–Boltzmann approaches, to systematically study the competitive uptake of mono- and divalent counterions by the macromolecule. We compare and fit our models to coarse-grained (implicit-solvent) computer simulation data of the globular polyelectrolyte dendritic polyglycerol sulfate (dPGS) in salt solutions of mixed valencies. The dPGS has high potential to serve in macromolecular carrier applications in biological systems and at the same time constitutes a good model system for a highly charged macromolecule. We finally use the simulation-informed models to extrapolate and predict electrostatic features such as the effective charge as a function of the divalent ion concentration for a wide range of dPGS generations (sizes).'
author:
- Rohit Nikam
- Xiao Xu
- Matej Kanduč
- Joachim Dzubiella
title: 'Competitive sorption of mono- versus divalent ions by highly charged globular macromolecules'
---

\[intro\] Introduction
======================

Polyelectrolytes in polar solvents such as water are important and ubiquitous in biological as well as in synthetic matter. [@Muthukumar2017; @Katchalsky1964; @Rubinstein2012; @Boroudjerdi2005; @Forster1995; @Dobrynin2005; @Liu2003] In these systems, electrostatic interactions, regulated by free ions and water, play a dominant role in shaping the structural and electrostatic characteristics of the polyelectrolyte, and the subsequent function of the system. [@Muthukumar2017; @Rubinstein2012; @Boroudjerdi2005] The electrostatic attraction between the isolated polyelectrolyte molecule and the oppositely charged counterions in the solution leads to strong counterion condensation on the molecule. This significantly modifies its interaction with other charged molecules (*e.g.*, proteins, DNA, etc.) and its electric properties such as the electrophoretic mobility in an external electric field. [@Boroudjerdi2005; @Forster1995; @Liu2003] Therefore, understanding counterion condensation is of utmost importance in order to understand the properties of polyelectrolytes and their implications in the biological and synthetic environments. [@Chremos2016; @Boroudjerdi2005] Condensation effectively leads to neutralizing an equivalent amount of the structural charge ${Z_\mathrm{d}}$ of the macromolecule. [@alexander1984charge; @belloni1998ionic] Hence, the charged substrate plus its confined counterions may be considered as a single entity with an effective (or renormalized) charge ${Z_\mathrm{eff}}$, which is significantly lower than the bare structural charge ${Z_\mathrm{d}}$. One can then identify the difference ${Z_\mathrm{d}}-{Z_\mathrm{eff}}$ as the amount of counterions condensed in the surface region. [@Bocquet2002]

The phenomenon of counterion condensation and the effect of ionic strength on the configurational properties of different types of polyelectrolyte molecules such as chains, [@Forster1992; @dobrynin1995; @Wenner2002; @Raspaud1998; @Dobrynin2005; @Liu2002; @Liu2003; @Muthukumar2004; @Chremos2016] brushes, [@ruhe2004polyelectrolyte; @pincus1991colloid; @borisov1991collapse; @Zhulina1995; @Zhulina2000] or polyelectrolyte nanogels [@nanogel1; @nanogel2; @nanogel3; @arturo2017; @arturo2018] have been studied extensively in the past. Through the knowledge of the distribution of the salt ions around the polyelectrolyte, *e.g.*, measured in terms of the radial distribution function in simulations and experiments, it is possible to derive important properties such as charge–charge correlation, osmotic compressibility and shear viscosity of the system. [@forster1995polyelectrolytes] Muthukumar, in his extensive and comprehensive review of the experimental, theoretical and simulation based research done on polyelectrolyte chains, described the effect of salt concentration, valency of counterions, chain length and polyelectrolyte concentration on counterion condensation. [@Muthukumar2017; @manning2012poisson] Besides the properties of a single isolated polyelectrolyte molecule, the ionic strength of the solution also influences the interaction of polyelectrolytes with other entities, such as adsorption on substrates, [@VandeSteeg1992; @Dahlgren1993; @Netz1999; @Hariharan1998; @Gittins2001; @caruso2000hollow] formation of ultra-thin polyelectrolyte multilayer membranes, [@Decher1992; @Decher1997; @Ladam2000; @McAloney2001; @Dubas2001] the structure and solubility of polyelectrolyte complexes [@hugerth1997effect; @rusu2003formation; @winkler2002complex; @Kudlay2004a; @Mende2002] or coacervates. [@spruijt2010binodal; @Gucht2011; @biesheuvel2004electrostatic; @Perry2014]

As an emerging class of functional polyelectrolytes, polyelectrolyte nanogels [@nanogel1; @nanogel2; @nanogel3; @arturo2017; @arturo2018] and dendritic or hyperbranched polyelectrolytes [@JensDernedde2010; @Khandare2012; @Groeger2013; @Maysinger2015; @Reimann2015] have attracted considerable interest in the scientific community in the last years due to their multifaceted bioapplications, such as biological imaging, drug delivery and tissue engineering. [@Leereview; @Ballauff2004; @Tian2013] In particular, the hyperbranched or dendritic polyglycerol sulfate molecules (hPGS or dPGS, respectively) are found to possess strong anti-inflammatory properties,[@Maysinger2015; @Reimann2015] act as a transport vehicle for drugs towards tumor cells,[@Sousa-Herves2015; @Groeger2013; @Vonnemann2014] and can be used as imaging agents for the diagnosis of rheumatoid arthritis. [@Vonnemann2014] This wide variety of applications, thus, have proven them to be high potential candidates for the use in medical treatments. [@Khandare2012] Hence, the understanding of dPGS interaction with the *in vivo* environment becomes important. The highly symmetric dendritic topology, terminated with monovalent negatively charged sulfate groups, makes dPGS also an excellent representative model in the class of highly charged globular polyelectrolytes. [@xu2017charged; @nikam2018charge] Because of the charged terminal groups, dPGS mainly interacts through electrostatics, rendering counterion condensation and subsequent charge renormalization effects to become substantial for function.

There have been past efforts to investigate the counterion condensation and to define the effective charge as a result of the charge renormalization on charged hard-sphere colloids. [@Ohshima1982; @Zimm1983; @alexander1984charge; @Belloni1984; @belloni1998ionic; @Ramanath1988; @Manning2007; @Bocquet2002; @Gillespie2014] However, the characterization of open-structure nanogel particles or dendrites such dPGS, which in part are penetrable to ions and a surface is not well defined, remains challenging. [@Ohshima2008] Recently, Xu *et al.* implemented a simple but accurate scheme to define and determine the effective surface potential and its location for dPGS, by mapping potentials obtained from simulations to the Debye–Hückel potential in the far-field regime. [@xu2017charged] This scheme is widely known as the Alexander prescription. [@alexander1984charge; @Trizac2002; @bocquet2002effective; @Levin2004] Based on this criterion, a systematic electrostatic characterization of dPGS has been performed via coarse-grained [@xu2017charged] and all-atom [@nikam2018charge] simulations by defining the number of condensed (bound) ions. It was then established that the strong binding of dPGS to lysozyme – an abundant protein in the human biological environment – and a sequential formation of a protein corona around dPGS in the presence of NaCl salt solution, is dominantly governed by the entropic gain due to the release of a few Na$^+$ counterions during binding. [@xu:biomacro] Proteins typically bind strongly to the macromolecular surface, thereby forming a protein ‘corona’, a dense shell of proteins that can entirely coat the macromolecule. [@Owens2006; @Cedervall2007; @Lindman2007; @Monopoli2012; @wang2013biomolecular; @LoGiudice2016; @Boselli2017]

Considering the medicinal applications of dPGS, it is important to study its interactions with [*divalent metal cations*]{}, *viz.* magnesium(II) and calcium(II) ions, which are key constituents of the human blood serum. Mg$^{2+}$ is essential for the stabilization of proteins, polysaccharides, lipids and DNA/RNA molecules, while Ca$^{2+}$ is critical for bone formation and plays a key role in signal transduction. [@Friesen2019a; @da2001biological] Human serum blood contains approximately $0.75-0.95$mM Mg$^{2+}$ ions, $1-4$mM Ca$^{2+}$ ions and around $150$mM NaCl salt in a dissociated form. [@meyers2004encyclopedia; @kretsinger2013encyclopedia] Thus, upon the administration of dPGS into the human biological environment, it is imperative for the competitive adsorption between the divalent (Mg$^{2+}$/Ca$^{2+}$) and monovalent (Na$^+$) ions to establish on the dPGS molecule, which can change the effective charge, and subsequently the interaction properties of dPGS with other charged entities such as proteins. This microscopic mechanism has a potential to significantly alter the attributes of protein corona around dPGS, [@XiaoCPS] thus, the biological immune response to the dPGS–protein corona complex, its metabolic fate, and the function of such a complex in biomedical and biotechnological applications. The competitive ion binding can be observed also in a wide variety of the biological and industrial ion-exchange processes such as the alkaline-earth/alkali-metal ion-exchange onto polyelectrolytes, [@Pochard1999] desalination of saline water to produce potable water, [@Birnhack2019] demineralization of whey, acid and alkali recovery from waste acid [@kobuchi1986application] and alkali solutions [@sata1993new] by diffusion dialysis, [@Sata2002] etc.

Interactions of multivalent ions with polyelectrolyte solutions have been theoretically studied in the past, in terms of their thermodynamic properties, [@Kuhn1999] ionic and potential distributions, [@Gavryushov1997] accurate calculation of the effective charge, [@DosSantos2010] and the effect on the interaction between polyelectrolyte macromolecules. [@Arenzon1999; @Naji2004; @Kanduc2010; @Rudi2016] In this paper, the focus is to theoretically analyze the competitive sorption of [*mono- versus divalent counterions*]{} by highly charged spherical dPGS-like polyelectrolytes with the help of mean-field continuum and discrete binding site models, informed by coarse-grained computer simulations of dPGS of various generations. The theoretical models are generally formulated for globular charged macromolecules and include ion-specific effects in a parametric way and can thus be straightforwardly modified or adapted to other charged globules, where mono-/divalent ion-exchange plays a role. In particular, we begin with the simple Donnan model, modified for ion-specific uptake, assuming that the electrostatic potential and the ionic concentrations are constant within the macromolecule phase and the bulk phase. [@Basser1993; @arturo2016; @Ahualli2014] Despite being simple, still, for the mixed case of monovalent and divalent ions the resultant composition is a non-trivial outcome. We continue with the mean-field Poisson–Boltzmann (PB) model, widely used in colloidal science and electrochemistry, [@Rubinstein2012; @israelachvili2011intermolecular; @adamson1967physical; @verwey1947theory; @Borukhov2000] and with the limitations well known and discussed, in particular the neglect of electrostatic and steric correlations, [@eigen1954; @kralj1996; @Cuvillier1997; @Borukhov1997; @DosSantos2010] or ion-specific sorption effects. [@Kalcher2010; @Kalcher2010a; @Chudoba2018; @arturo2014; @Yaakov2009a; @Lima2008; @Koelsch2007; @okur2017beyond; @LoNostro2012; @Schwierz2013] The PB model has also been implemented to address the problem of competitive counterion binding in a mixed salt for the cases of linear polyelectrolytes such as DNA [@Burak2004; @Chen2002; @Rouzina1997; @Misra1994; @Paulsen1987] and planar geometries. [@Rouzina1994] We also devise a two-state approximation model for an ion condensation around a charged globule.

![image](images/fig1.pdf){width="14cm" height="14cm"}

The two-state approach was firstly used in the Oosawa–Manning model [@Manning1969; @Oosawa] for the counterion condensation around polyelectrolyte chains, according to which, counterions in a solution can be classified into two categories: ‘free’ counterions, which are able to explore the whole solution volume $V$ and the ‘condensed’ (or ‘bound’) counterions, which are localized within a small volume surrounding the polyelectrolyte macromolecule. An equivalent model for an impenetrable sphere with a surface charge was developed by Manning, where the number of condensed counterions on the macromolecule per bare unit surface charge is obtained by a free energy minimization, pointing to the competition between the electrostatic binding of counterions to the macromolecule and their dissociation entropy. [@Manning2007] We extend this model by introducing a discrete binding site model by considering the finite configurational volume of the ion in the condensed state and that the macromolecule has a finite number of charged binding sites by adopting the mixing entropy from the works of McGhee and von Hippel. [@Mcghee1974] Ion-binding models in the same spirit have been developed in the past to describe the ionization equilibrium of [*linear*]{} polyelectrolytes in monovalent salt [@flory1953molecular; @Raphael1990; @Muthukumar2004], multivalent salt, [@Friedman1984] and in mixtures of mono- and divalent salts. [@Kundagrami2008] All our models are compared to molecular simulations and used to study systematically the key electrostatic features of a highly charged globule, such as the effect of competitive adsorption on the variation of the number of condensed monovalent and divalent counterions, effective charge, and its variation with divalent ion concentration.

Coarse-grained computer simulations {#cg_sim}
===================================

Simulation methods, force fields, and systems {#sec:cg_ff}
---------------------------------------------

The coarse-grained (CG) monomer-resolved models of the dPGS macromolecule have been developed previously [@xu2017charged] and maintain the essential dPGS structural and electrostatic features with affordable computing expense. In brief, the dPGS branching units (–) and inner core (–) (both of which are a part of the glycerol chemical group, respectively), and the terminal sulfate groups (–) are individually represented by the CG segments of specific type. The gross number of the CG segments is equal to the dendrimer polymerization $N_g =  3 \times 2^{n+1} - 2$ of generation index $n$. Only the terminal segments are charged with $-1e$ (where $e$ is the elementary charge), leading to the dPGS bare valency $|Z_n| = 3 \times 2^{n+1}$. The CG segments are connected by bonded and angular potentials both in harmonic form. In the previous work [@xu2017charged] we only studied monovalent ions. Here we extend it to study the competitive uptake of mono- and divalent ions for generations 2 and 4. The bare charge valencies of the G$_2$-dPGS are thus ${Z_\mathrm{d}}= Z_{n=2} = -24$ and ${Z_\mathrm{d}}= Z_{n=4} = -96$. Snapshots are shown in Fig. \[snap\].

The non-bonded interactions between CG beads are described by the Lennard-Jones (LJ) potential together with the Lorentz–Berthelot mixing rules. In particular, the energy parameter $\epsilon_{\rm LJ} = 0.1\,{k_\mathrm{B}T}$ and the diameter $\sigma_{\rm LJ} = 0.4$ nm are set identical for all ions (mono- and divalent) and thus any ion-specific effects are not explicitly included. In our simulations we place the dPGS in the center of a periodically repeated cubic box with a volume of $V$ (side-length of $L = 30$ nm). The solvent is implicitly assumed as a dielectric continuum with a dielectric constant $\epsilon_{\rm w} = 78$. The CG simulations employ the stochastic dynamics (SD) integrator in Gromacs $4.5.5$ as in our previous work. [@xu2017charged]

![image](images/fig2.pdf){width="12.5cm" height="12.5cm"}

![image](images/fig3.pdf){width="13cm" height="13cm"}

All simulations are performed in the canonical ensemble. The divalent cations (DCs), monovalent cations (MCs) and monovalent anions in the system are referred to with subscripts $++,\; +$ and $-$, respectively. The dPGS is accompanied by the corresponding number of monovalent counterions ${N_\mathrm{s}}$ ($24$ for G$_2$-dPGS and $96$ for G$_4$-dPGS) electrically neutralizing the macromolecule and having the same chemical identity as the MCs of the salt. The number of salt ions $i$ ($i = +\!+,+,-$) is denoted as $n_i$, while the corresponding total salt concentrations are denoted as ${c^0_{i}}=n_i/V$. Bulk concentrations are defined as ${c^{\mathrm{b}}_{i}}=(n_i - {N^\mathrm{b}_{i}})/(V - {v_\mathrm{eff}})$ (for $i=+\!+,-$) and ${c^{\mathrm{b}}_{{\scalebox{.8}{$\scriptscriptstyle +$}}}}=(n_{{\scalebox{.8}{$\scriptscriptstyle +$}}} +{N_\mathrm{s}}- {N^\mathrm{b}_{{\scalebox{.8}{$\scriptscriptstyle +$}}}})/(V - {v_\mathrm{eff}})$, where ${v_\mathrm{eff}}= 4\pi {r_\mathrm{eff}}^3/3$ is the volume enclosed by the effective radius ${r_\mathrm{eff}}$ of dPGS and ${N^\mathrm{b}_{i}}$ is the number of ions $i$ condensed (bound) on the dPGS. The definitions of both ${r_\mathrm{eff}}$ and ${N^\mathrm{b}_{i}}$ are discussed in Section \[sim\_analysis\].

The simulations are performed at the total DC concentrations ${c^0_{{\scalebox{.8}{$\scriptscriptstyle ++$}}}}$ of 0.98, 2.95, 3.75, 9.96 and 14.94mM. G$_2$-dPGS simulation snapshots for different ${c^0_{{\scalebox{.8}{$\scriptscriptstyle ++$}}}}$ values are shown in Fig. \[snap\](a)-(c), while the whole simulation box is displayed in Fig. \[snap\](d). The MC concentration ${c^0_{{\scalebox{.8}{$\scriptscriptstyle +$}}}}$ is fixed to $150.37$mM and the monovalent anion concentration is adjusted in a way to ensure electroneutrality in the simulation box. The bulk ionic strength $I = \frac{1}{2}\sum_i z^2_i {c^{\mathrm{b}}_{i}}$ ($i = +,+\!+,-$ with the charge valency $z_i$) ranges from $150.5$mM to $195$mM. The corresponding Debye screening length $\kappa^{-1} = \left(8\pi {l_\mathrm{B}}I \right)^{-1/2}$ (where ${l_\mathrm{B}}$ is the Bjerrum length) ranges from $0.8$nm (${c^{\mathrm{b}}_{{\scalebox{.8}{$\scriptscriptstyle ++$}}}}=0$ and ${c^{\mathrm{b}}_{{\scalebox{.8}{$\scriptscriptstyle +$}}}}=150.5$mM) to $0.7$nm (${c^{\mathrm{b}}_{{\scalebox{.8}{$\scriptscriptstyle ++$}}}}=14.94$ and ${c^{\mathrm{b}}_{{\scalebox{.8}{$\scriptscriptstyle +$}}}}=150.5$mM). As a reference, we also perform CG simulations in the limit of only monovalent salt, with total concentrations ${c^0_{{\scalebox{.8}{$\scriptscriptstyle +$}}}}$ of 10.02, 25.06 and 150.37mM.

Simulation results: radial density distributions {#sim_results}
------------------------------------------------

The dPGS structure and its response to the addition of the DCs, is examined by the density distribution of the terminal sulfate beads ${c_\mathrm{s}}(r)$ as a function of the distance $r$ from the center-of-mass (COM) of the dPGS, for different DC concentrations ${c^0_{{\scalebox{.8}{$\scriptscriptstyle ++$}}}}$, as shown in Fig. \[rho\_sulfate\]. Interestingly, the presence of DCs does not lead to a notable change in the dPGS structure. Instead, the ${c_\mathrm{s}}(r)$ profiles in the operated range of ${c^0_{{\scalebox{.8}{$\scriptscriptstyle ++$}}}}$ and for both G$_2$-dPGS and G$_4$-dPGS are reasonably coincident. Fig. \[rho\_sulfate\](a) shows that for G$_2$-dPGS, a single-peak distribution is found, indicating that most of the sulfate beads reside on the molecular surface. However, in Fig. \[rho\_sulfate\](b), a bimodal distribution is seen for G$_4$-dPGS with a small peak at $r \simeq 0.6$nm. This backfolding phenomenon, contributing to a dense-core arrangement due to the dense macromolecular shell, [@Ballauff2004] is also found in our previous works [@xu2017charged; @nikam2018charge] and has been detected for other terminally charged CG dendrimer models. [@Huismann2010; @Huismann2010B; @Klos2010; @Klos2011] After the major peak, ${c_\mathrm{s}}(r)$ gradually subsides to zero. The location where the charge density ${c_\mathrm{s}}(r)$ falls to $150$mM, which we set as the physiological NaCl concentration, is defined as the bare (intrinsic) radius of dPGS ${r_\mathrm{d}}$, [^1] shown as vertical dashed blue lines in Fig. \[rho\_sulfate\]. The ${r_\mathrm{d}}$ values for G$_2$-dPGS and G$_4$-dPGS are obtained as $1.40$nm and $2.11$nm, respectively. Fig. \[rho\_sulfate\](b) also shows that a slight shift in the location of the major peak and an enrichment of the lower peak appears as ${c^0_{{\scalebox{.8}{$\scriptscriptstyle ++$}}}}$ increases, indicating a slow shrinking of the dPGS molecule due to the condensation of DCs (see Fig. \[density\_ions\]).

Figs. \[density\_ions\](a) and (b) show the cation density distributions $c_i(r)$ ($i =+,+\!+$) for G$_2$-dPGS and G$_4$-dPGS, respectively. Let us focus first on G$_2$ in Fig. \[density\_ions\](a). The MC distribution $c_{{\scalebox{.8}{$\scriptscriptstyle +$}}}(r)$ shows a high accumulation of counterions close to the sulfate groups, with a global maximum at distances $r \sim 1.2$nm, slightly larger than the sulfate peak (peaking roughly at $\sim 1$nm). This means that the most strongly bound ‘condensed’ MCs are thus distributed more on the surface layers of the dPGS. At larger distances, $r \sim 2$nm, a Debye–Hückel like decay is observed. Adding more DCs, the MC distribution gradually diminishes, as expected from the exchange of MCs with DCs within the dPGS. However, interestingly, the DC distribution peaks at distances distinctively smaller than the location of the sulfate peak, roughly $0.5-0.6$nm shifted towards the dPGS center away from the peak of the MC distribution. This more interior binding might be attributed to different binding mechanisms between DCs and sulfate, *e.g.*, bridging of two sulfate groups by one DC, which might be sterically favored closer to the dPGS core. These subtle structural effects may have important consequences in the context of the counterion-release mechanism driving the dPGS–protein binding, [@xu:biomacro] which should be interesting for future studies. The ion profiles for G$_4$ shown in Fig. \[density\_ions\](b) show qualitatively the same behavior but are broader and double-peaked because of the significant sulfate backfolding as previously presented in Fig. \[rho\_sulfate\](b).

It is worth noting that simulations of DCs in general are more challenging than for MCs only. DC are more heavily hydrated than MCs (*e.g.*, Mg$^{2+}$ and Na$^+$ ions), [@Stokes1948; @Marcus2006] therefore future studies should scrutinize the ionic size used in the implicit solvent. Furthermore, quantum mechanical charge transfer effects as a result of the ion-induced powerful electronic polarization of the surrounding media, [@Yao2015] which are much more prevalent in the case of DCs [@Pavlov1998; @Kohagen2014] than MCs, may also be subsumed in ionic sizes in the implicit water. These model details may subtly change the density profiles shown in Figs. \[density\_ions\](a) and (b). However, the effects on total competitive uptake should be relatively minor as they are dominantly driven by valency and electrostatic correlations, and size effects are typically of second order importance.

Using the density distributions of the charged entities shown above, the electrostatic properties of dPGS can be studied in the presence of the mixture of DCs and MCs. The analysis methods described in Section \[sim\_analysis\] are used to define the effective radius ${r_\mathrm{eff}}$, charge valency ${Z_\mathrm{eff}}$ and potential ${\phi_\mathrm{eff}}$ of dPGS.

Electrostatic properties of dPGS {#sim_analysis}
--------------------------------

**dPGS effective radius** The first step to study the ion condensation behavior is to adopt a characteristic distance ${r_\mathrm{eff}}$ to distinguish a bound ion from an unbound one. A practical method in that respect has been summarized in our previous work. [@xu2017charged] In short, we first consider the dPGS radial electrostatic potential profile $\phi$ (scaled by ${k_\mathrm{B}T}/e$), through the framework of the Poisson’s equation $$\nabla^2 \phi = -4\pi {l_\mathrm{B}}\sum_{i} z_i c_i(r) \qquad  i= s, ++, +, -
\label{possion}$$ Here, $c_i(r)$ refers to the radial number density profiles with respect to the distance to the dPGS-COM $r$ for all charged species in the CG simulation, namely, sulfates ($s$), DCs, MCs, and monovalent anions. For all ionic species, $c_i(r)$ reaches the bulk number density ${c^{\mathrm{b}}_{i}}$ in the far-field. The simulation results for the profiles are shown in Figs. \[rho\_sulfate\] and \[density\_ions\]. The Poisson’s equation is numerically integrated twice to obtain $\phi(r)$, which is then compared with the dimensionless Debye–Hückel potential $\phi^{}_\mathrm{DH}$, given by [@alexander1984charge; @xu2017charged; @nikam2018charge] $$\phi^{}_\mathrm{DH}(r) = {Z_\mathrm{eff}}{l_\mathrm{B}}\frac{\mathrm{e}^{\kappa {r_\mathrm{eff}}}}{1+\kappa {r_\mathrm{eff}}} \frac{\mathrm{e}^{-\kappa r}}{r}.
\label{DH_eq}$$ $\phi^{}_{\rm DH}$ is applicable to a charged sphere with radius ${r_\mathrm{eff}}$ and valency ${Z_\mathrm{eff}}$. It approaches to $\phi$ only after the distance $r^*$ where non-linear effects, including the correlation and condensation of ions, subside. Thus, $r^* = {r_\mathrm{eff}}$ is eligible to serve as the dPGS effective radius to define the bound ions. The effective surface potential of dPGS obtained from simulations is then defined as ${\phi_\mathrm{eff}}= \phi({r_\mathrm{eff}})$, which is shown in Table \[cg\_table\]. Comparing Eq.  to the radial electrostatic potentials from the simulations, [@xu2017charged] the value of ${r_\mathrm{eff}}$ for dPGS in the simulations for G$_2$ and G$_4$ was found to be $1.65$nm and $2.40$nm, respectively, under the operated concentration range in the mixture of DCs with MCs as well as in the monovalent limit, as shown in Table \[cg\_table\]. These values are different than the ones obtained in our previous work, [@xu2017charged] which operates at ${c^0_{{\scalebox{.8}{$\scriptscriptstyle +$}}}}= 10$mM, unlike the current work where ${c^0_{{\scalebox{.8}{$\scriptscriptstyle +$}}}}= 150.37$mM. The newly obtained ${r_\mathrm{eff}}$ values in this work are then used as an input for the MMvH model, as discussed in Section \[mvh\], to describe the competitive sorption. It is thus implicitly assumed that ${r_\mathrm{eff}}$ does not depend on the sorption of DCs, within the operated range of DC concentrations. The same prescription will be used to define ${r_\mathrm{eff}}$ (denoted as ${r_\mathrm{eff}}^\mathrm{PB}$) from the solutions of the PPB model, as discussed in Section \[sec:ppbl\]. The results for ${r_\mathrm{eff}}^\mathrm{PB}$ are also shown in Table \[cg\_table\].

  ---------------------------------------------------- ------------------ -------------------- -------------------- ----------------------- ------------------ -------------------- -------------------- ----------------------- -------------------------------- -------------------------------- ----------------------------------- -------------------------------- -------------------------------- -----------------------------------
                                                                                                                                                                                                                                                                                                                                                                                                         
                                                                                                                                                                                                                                                                                                                                                                                                         
   ${c^0_{{\scalebox{.8}{$\scriptscriptstyle ++$}}}}$   ${r_\mathrm{d}}$   ${r_\mathrm{eff}}$   ${Z_\mathrm{eff}}$   ${\phi_\mathrm{eff}}$   ${r_\mathrm{d}}$   ${r_\mathrm{eff}}$   ${Z_\mathrm{eff}}$   ${\phi_\mathrm{eff}}$   ${r_\mathrm{eff}}^\mathrm{PB}$   ${Z_\mathrm{eff}}^\mathrm{PB}$   ${\phi_\mathrm{eff}}^\mathrm{PB}$   ${r_\mathrm{eff}}^\mathrm{PB}$   ${Z_\mathrm{eff}}^\mathrm{PB}$   ${\phi_\mathrm{eff}}^\mathrm{PB}$
                        ([mM]{})                            ([nm]{})            ([nm]{})                                                         ([nm]{})            ([nm]{})                                                                ([nm]{})                                                                                              ([nm]{})                                              
                         $0.00$                                                                      $-10.09$               $-1.26$                                                       $-20.04$               $-1.27$                                                      $-11.72$                          $ -1.32$                                                            $-23.60$                           $-1.56$
                         $0.98$                                                                      $-8.85$                $-1.15$                                                       $-17.75$               $-1.14$                                                      $-9.79$                            $-1.12$                                                            $-20.03$                           $-1.38$
                         $2.95$                                                                      $-7.40$                $-0.98$                                                       $-14.21$               $-0.93$                                                      $-8.89$                            $-0.88$                                                            $-15.54$                           $-1.05$
                         $3.75$                                                                      $-6.84$                $-0.85$                                                       $-12.25$               $-0.77$                                                      $-8.29$                            $-0.83$                                                            $-14.34$                           $-0.97$
                         $9.96$                                                                      $-6.33$                $-0.75$                                                       $-10.11$               $-0.62$                                                      $-7.03$                            $-0.57$                                                            $-8.86$                            $-0.60$
                        $14.94$                                                                      $-5.86$                $-0.68$                                                       $-9.65$                $-0.55$                                                      $-6.36$                            $-0.46$                                                            $-6.13$                            $-0.44$
  ---------------------------------------------------- ------------------ -------------------- -------------------- ----------------------- ------------------ -------------------- -------------------- ----------------------- -------------------------------- -------------------------------- ----------------------------------- -------------------------------- -------------------------------- -----------------------------------

**Number of bound ions and effective charge** The cumulative number of ions of species $i$ as a function of the distance $r$ from the COM of dPGS is calculated as $$N_{\mathrm{acc},i}(r) = \int^{r}_{0} c_i(r') 4\pi {r'}^2 \mathrm{d}r' \quad  \quad i = ++,+,-.
\label{nbi}$$ Summing up the contribution of all charged species, the cumulative charge valency of the system as a function of the distance $r$ reads $$Z_\mathrm{acc}(r) = {Z_\mathrm{d}}(r) + 2N_{\mathrm{acc},{\scalebox{.8}{$\scriptscriptstyle ++$}}}(r) + N_{\mathrm{acc},{\scalebox{.8}{$\scriptscriptstyle +$}}}(r) - N_{\mathrm{acc},{\scalebox{.8}{$\scriptscriptstyle -$}}}(r),
\label{neff}$$ where ${Z_\mathrm{d}}(r)$ denotes the spatial distribution of bare charge valency of the dPGS, obtained from the simulation. With that, the number of the bound ions and the effective charge valency of the dPGS follow as ${N^\mathrm{b}_{i}}= N_{\mathrm{acc},i}({r_\mathrm{eff}})$ and ${Z_\mathrm{eff}}= Z_{\mathrm{acc}}({r_\mathrm{eff}})$, respectively. The values are shown in Table \[cg\_table\]. ${Z_\mathrm{eff}}$ and ${\phi_\mathrm{eff}}$ exhibit a strong decrease in the magnitude with higher ${c^0_{{\scalebox{.8}{$\scriptscriptstyle ++$}}}}$, indicating an enhanced dPGS charge renormalization.

Theoretical models
==================

Basic model {#basic}
-----------

In our theoretical models, the macromolecule is represented as a perfect sphere with the bare radius ${r_\mathrm{d}}$, the bare charge valency ${Z_\mathrm{d}}$, the effective radius ${r_\mathrm{eff}}$ and the effective charge valency ${Z_\mathrm{eff}}$, enclosed in a spherical domain of radius $R$ and volume $V$, as shown in Fig. \[fig:theoretical-model\]. The total number of charged monomers in the macromolecule is ${N_\mathrm{s}}$, each of which is negatively charged with a charge valency ${z_\mathrm{s}}$.

![\[fig:theoretical-model\] Schematic of a theoretical model representing the system shown in Fig. \[snap\](d). The computational cell domain (blue) is assumed to be spherical with the same volume as that of the simulation box, $V$, and with a uniform dielectric constant of water $\epsilon_\mathrm{w}=78$. dPGS is assumed to be a perfect sphere (orange) at the center of the domain. The dPGS bare and effective charge valencies are ${Z_\mathrm{d}}$ and ${Z_\mathrm{eff}}$, respectively. ${r_\mathrm{d}}$ is the bare radius of dPGS, while ${r_\mathrm{eff}}$, the effective radius, representing the distance separating the electric double layer regime ($r > {r_\mathrm{eff}}$) from the non-linear counterion ’condensation’ regime ($r < {r_\mathrm{eff}}$).](images/fig4.pdf){width="5.5cm" height="5.5cm"}

All ionic species and the macromolecule are assumed to be in an aqueous bath with an implicitly modeled solvent, having a uniform dielectric constant $\epsilon_\mathrm{w}=78$ at a temperature $T=298$K.

\[dm\] The Donnan model (DM)
----------------------------

The arguably simplest model for competitive uptake is the Donnan model. The Donnan equilibrium assumes two strictly electroneutral and mutually exclusive regions, *i.e.*, the macromolecule region with the Donnan radius set to be the bare radius ${r_\mathrm{d}}$ taken from simulations (*i.e.*, with a bare macromolecular volume ${v_\mathrm{d}}=4\pi {r_\mathrm{d}}^3 /3$) and total homogeneously distributed bare charge of valency ${Z_\mathrm{d}}= {z_\mathrm{s}}{N_\mathrm{s}}$ with a concentration ${c_\mathrm{s}}= {N_\mathrm{s}}/{v_\mathrm{d}}$ of charged groups of the macromolecule, and the bulk region outside the molecule with a bulk ion concentration ${c^{\mathrm{b}}_{i}}$ ($i=+,+\!+,-$). Charge neutralization of the macromolecule by the counterions leads to the Donnan potential, which is a potential having a constant non-zero value in the macromolecule region. The potential in the bulk region is set to zero. The equilibrium distribution (partitioning) of ions among the regions results in the concentrations of ionic species $i$ as ${c^\mathrm{m}_{i}}$ and ${c^{\mathrm{b}}_{i}}$ in the macromolecule and bulk regions, respectively. These concentrations are related via the partition coefficient ${\mathcal{K}_{i}}$, given by $${\mathcal{K}_{i}}= \frac{{c^\mathrm{m}_{i}}}{{c^{\mathrm{b}}_{i}}} \quad \quad i=+\!+,+,-
\label{ki-don1}$$ Neglecting ion–ion correlations, an approximate expression for ${\mathcal{K}_{i}}$ can be obtained using the condition that the equilibrium electrochemical potential of ion $i$ is equal in both the macromolecule and bulk regions, implying that $${\mathrm{ln}}\,{c^{\mathrm{b}}_{i}}=  z_i{\phi^{}_\mathrm{D}}+ {\mathrm{ln}}\,{c^\mathrm{m}_{i}}+ \beta {\Delta \mu_{\mathrm{int},\,i}}\label{phase_eqm}$$ where ${\phi^{}_\mathrm{D}}$ is the dimensionless Donnan potential (scaled by ${k_\mathrm{B}T}/e$) in the macromolecule region and $\beta^{-1} ={k_\mathrm{B}T}$ is the thermal energy. With ${\Delta \mu_{\mathrm{int},\,i}}$ we account for additional non-electrostatic effects that can drive adsorption, *e.g.*, dispersion and hydrophobic forces in the net ion–macromolecule interaction, and is termed the ion-specific binding chemical potential of the condensed ion. The inclusion of ${\Delta \mu_{\mathrm{int},\,i}}$ has been considered in previous work, for example, as a term reflecting the steric ion–ion packing effects in a Donnan model for ion binding by polyelectrolytes or charged hydrogels. [@Chudoba2018; @arturo2014; @Ahualli2014]

Eq.  with the help of Eq.  then leads to $${\mathcal{K}_{i}}= \frac{{c^\mathrm{m}_{i}}}{{c^{\mathrm{b}}_{i}}}  = \mathrm{e}^{-\beta {\Delta \mu_{\mathrm{bind},\,i}}}= \mathrm{e}^{-\beta {\Delta \mu_{\mathrm{int},\,i}}} \mathrm{e}^{-z_i {\phi^{}_\mathrm{D}}}
\label{ki_don2}$$ where ${\Delta \mu_{\mathrm{bind},\,i}}$ is the total transfer chemical potential for ion $i$ from the bulk to the macromolecule region. This allows us to define the intrinsic partition ratio for ionic species $i$ as $${\mathcal{K}_{\mathrm{int},\,{i}}}= \mathrm{e}^{-\beta {\Delta \mu_{\mathrm{int},\,i}}} \quad \quad i=+,+\!+
\label{kiz}$$ and the Donnan partition ratio as a contribution from pure electrostatic interaction between the ion and the macromolecule environment as $${\mathcal{K}_{\mathrm{el},\,{i}}}= \mathrm{e}^{-z_i {\phi^{}_\mathrm{D}}} \quad \quad i=+,+\!+
\label{kid}$$ The electrostatic component of total binding chemical potential of a counterion $i$ is then defined as $\beta {\Delta \mu_{\mathrm{el},\,i}}= -{\mathrm{ln}}\,{\mathcal{K}_{\mathrm{el},\,{i}}}= z_i{\phi^{}_\mathrm{D}}$. Eq.  can then be conveniently shortened as $${\mathcal{K}_{i}}= {\mathcal{K}_{\mathrm{el},\,{i}}}\, {\mathcal{K}_{\mathrm{int},\,{i}}}\label{ki-don3}$$ where ${\mathcal{K}_{i}}$ is shown as a composition of intrinsic and electrostatic effects. The signature assumption behind the Donnan model is the electroneutrality in the macromolecule region expressed as $${z_\mathrm{s}}{c_\mathrm{s}}+ \sum_i z_i {c^{\mathrm{b}}_{i}}\,{\mathcal{K}_{i}}= 0   
\label{neutral}$$ Solving Eq.  for ${\phi^{}_\mathrm{D}}$ enables us to evaluate the net partition coefficient ${\mathcal{K}_{i}}$. Eq.  has no closed solution for multivalent ions, but it exists for the case of only monovalent ions in the system, ($i=\pm$) and is given as [@cemil2] $${\phi^{}_\mathrm{D}}= -{\mathrm{ln}}\left(-\frac{ \sqrt{1 + \chi_{{\scalebox{.8}{$\scriptscriptstyle +$}}} \chi_{{\scalebox{.8}{$\scriptscriptstyle -$}}}} + 1 }{\chi_{{\scalebox{.8}{$\scriptscriptstyle +$}}}} \right)
\label{donnan-mono-expr}$$ where $\chi^{}_i = 2{\mathcal{K}_{\mathrm{int},\,{i}}}\, {c^{\mathrm{b}}_{i}}/{z_\mathrm{s}}{c_\mathrm{s}}$. Note that $\chi^{}_i<0$, since the valency of charged groups ${z_\mathrm{s}}$ is negative. Using Eqs. , and , for the monovalent-only case, the number of ions of species $i(=\pm)$ partitioned into the macromolecule region is then given as $$N^\mathrm{b}_{\pm} = c^0_{\pm} {v_\mathrm{d}}\, \mathcal{K}_{\mathrm{int}, \pm}  \left( -\frac{ \sqrt{1 + \chi_{{\scalebox{.8}{$\scriptscriptstyle +$}}} \chi_{{\scalebox{.8}{$\scriptscriptstyle -$}}}} + 1 }{\chi_{{\scalebox{.8}{$\scriptscriptstyle +$}}}} \,\right)^{\pm 1}  
\label{Nb-don}$$ To evaluate the competition between MCs and DCs in the Donnan model we evaluate Eqs.  and numerically, *cf.* section \[numerical\].

Because of the electroneutrality assumption, the Donnan prediction for the amount of counterion sorption by the macromolecule in the monovalent-only case is given by ${N^\mathrm{b}_{{\scalebox{.8}{$\scriptscriptstyle +$}}}}= |Z| + N^\mathrm{b}_{{\scalebox{.8}{$\scriptscriptstyle -$}}}$. For highly charged macromolecules, *i.e.*, $\chi_i \rightarrow 0$, Eq.  this trivially yields ${N^\mathrm{b}_{{\scalebox{.8}{$\scriptscriptstyle +$}}}}\simeq |Z|$. For the competitive sorption case, however, the result is non-trivial and can give a useful orientation with little effort. The Donnan model should become more quantitative for large dPGS generations, i.e., large size and/or high salt concentrations, so that $\kappa {r_\mathrm{d}}\gg 1$, for which the electroneutrality assumption is then very well justified.

\[sec:ppbl\] Ion-specific Penetrable Poisson–Boltzmann (PPB) model
------------------------------------------------------------------

We now put forward a penetrable PB (PPB) model in which the charge profiles can be resolved in $r$, the radial distance from the macromolecular center. Since our charged macromolecules we have in mind (dPGS above and similar) are polymer-based with open structures and typically internally smeared out charge distributions, we opted (as in the Donnan model) for a penetrable model instead of a PB model for surface adsorption as typically used in studies of colloidal charge renormalization. [@Wall1957; @Ohshima1982; @Ohshima2008] Based on the parametrization described in the basic model (Section \[basic\]), we assume the macromolecule as a perfect penetrable sphere with a charge valency ${Z_\mathrm{d}}= {z_\mathrm{s}}{N_\mathrm{s}}$ and radius ${r_\mathrm{d}}$, as shown in Fig. \[fig:theoretical-model\]. ${r_\mathrm{d}}$ is taken from the dPGS internal charge distribution obtained from simulations, *cf.* Section \[sim\_results\] and Fig. \[rho\_sulfate\]. The charged monomers of the macromolecule, thus, have a uniform number distribution ${c_\mathrm{s}}= {N_\mathrm{s}}/{v_\mathrm{d}}$ (where ${v_\mathrm{d}}= 4\pi {r_\mathrm{d}}^3/3$) within the volume ${v_\mathrm{d}}$. ${c_\mathrm{s}}$ is applicable only within the macromolecule domain, *i.e.*, ${c_\mathrm{s}}(r) = {c_\mathrm{s}}\left(1-H(r-{r_\mathrm{d}})\right)$, where $H(r)$ is the Heaviside-step function. As an improvement to the standard PB model, here we also consider a contribution of the intrinsic non-electrostatic ion-specific interaction ${\Delta \mu_{\mathrm{int},\,i}}$ between the ion and the macromolecule, [@Kalcher2010; @Kalcher2010a] analogous to Eq.  in the Donnan model above. Assuming the electrostatic potential far away from the macromolecule, $\phi\left(r \rightarrow R \right)=0$, we first balance the chemical potential for each ion, between the bulk regime far from the macromolecule and the regime at the finite distance $r$ from the center of the macromolecule $${\mathrm{ln}}\,{c^{\mathrm{b}}_{i}}= z_i\phi(r) + {\mathrm{ln}}\,c_{i}(r) + \beta {\Delta \mu_{\mathrm{int},\,i}}(r), 
\label{pb_chem_pot_bal}$$ which is similar to Eq. , but in a distance-resolved manner. ${\Delta \mu_{\mathrm{int},\,i}}$ is considered on a local level, *i.e.*, ${\Delta \mu_{\mathrm{int},\,i}}(r) = {\Delta \mu_{\mathrm{int},\,i}}\left(1-H(r-{r_\mathrm{d}})\right)$. The Boltzmann ansatz then becomes $$c_i(r) = {c^{\mathrm{b}}_{i}}\,\mathrm{e}^{-z_i\phi(r) - \beta {\Delta \mu_{\mathrm{int},\,i}}(r) }
\label{pb_ansatz0}$$ The distance-resolved electrostatic potential can be calculated from Eq.  together with the Poisson’s equation as $$\nabla^2 \phi(r) = -4\pi {l_\mathrm{B}}\left(\sum_{i} z_i c_i(r) + {z_\mathrm{s}}{c_\mathrm{s}}(r) \right) \quad \quad i = +\!+,+,-
\label{pb-pot}$$ which establishes the PPB model including ion-specific binding effects. The boundary conditions used are $({{\mathop{}\!\mathrm{d}}}\phi/{{\mathop{}\!\mathrm{d}}}r)\left(r \rightarrow 0 \right) = 0$ and $({{\mathop{}\!\mathrm{d}}}\phi/{{\mathop{}\!\mathrm{d}}}r)\left(r \rightarrow R \right)= 0$.

An effective radius for dPGS is calculated independently for this model (labeled ${r_\mathrm{eff}}^\mathrm{PB}$) using the Alexander prescription [@alexander1984charge; @Trizac2002; @bocquet2002effective; @Levin2004] on the obtained potential $\phi$, the same recipe used to calculate ${r_\mathrm{eff}}$ from simulations, *cf.* Section \[sim\_analysis\]. The values of ${r_\mathrm{eff}}^\mathrm{PB}$ for G$_2$-dPGS and G$_4$-dPGS are obtained as $1.42$nm and $2.36$nm, respectively, under the operated range of ${c^0_{{\scalebox{.8}{$\scriptscriptstyle ++$}}}}$ and at ${c^0_{{\scalebox{.8}{$\scriptscriptstyle +$}}}}=150.37$mM. The ${r_\mathrm{eff}}^\mathrm{PB}$ values are thus found to be close to those obtained from the simulations, as shown in Table \[cg\_table\]. The effective surface potential of the macromolecule is then given by ${\phi_\mathrm{eff}}^\mathrm{PB} = \phi({r_\mathrm{eff}}^\mathrm{PB})$. The number of bound ions of species $i$ within ${r_\mathrm{eff}}$, is then given by $${N^\mathrm{b}_{i}}= \int_0^{{r_\mathrm{eff}}} c_i(r)\, 4 \pi r^2 {\mathop{}\!\mathrm{d}}r \quad  \quad  i=+,+\!+
\label{Nib}$$ The corresponding effective charge valency ${Z_\mathrm{eff}}^\mathrm{PB}$ is calculated using Eq. . The PPB equations are solved numerically, *cf.* Section \[numerical\].

\[mvh\] Manning–McGhee–von Hippel binding model (MMvH)
------------------------------------------------------

In this section, we introduce a model based on a discrete two-state (condensed or free) perspective for the counterions, built to capture the essential physics of polyelectrolyte–ion binding in an accurate but minimalistic fashion. The model is an extension of ideas by Manning, [@Manning2007] in which ion-condensation on charged spherical surfaces was described on a mean-field free energy level as a competition between the charging (Born) self-energy of the macromolecule in salt solution and the entropy cost of binding for one-component counterions. Here, we extend this model to the case of mixtures of MCs and DCs, including binding saturation for a fixed number of binding sites like in Langmuir isotherms. The extension of the latter to binary binding of one or two binding sites by mono- or divalent solutes, respectively, was put forward buy McGhee and von Hippel. [@Mcghee1974] Therefore, we name the model Manning–McGhee–von Hippel binding model (MMvH).

Following Manning, [@Manning2007] we treat the macromolecule as an impenetrable sphere of radius ${r_\mathrm{eff}}$ and charge valency ${Z_\mathrm{d}}= z_s{N_\mathrm{s}}$ taken from simulations, and extend the Manning’s model into a discrete binding site model, where the ${N_\mathrm{s}}$ charged monomers act as a finite collection of discrete binding sites for both the MCs and DCs. For the case of the DCs, two adjacent charged monomers can collectively act as a single binding site for a DC. The resulting combinatorial ways to arrange the bound MCs and DCs lead to mixing entropies worked out by McGhee and von Hippel. [@Mcghee1974] Pertaining to the canonical ensemble, we fix the total number of salt ions $n_i$, the corresponding concentrations ${c^0_{i}}$ ($i= +\!+,+,-$), the number of monovalent counterions ${N_\mathrm{s}}$ to the macromolecule, the total number of binding sites on the macromolecule and the total domain volume $V$. The coions in this model simply serve the function of maintaining electroneutrality in the total domain and their explicit adsorption is neglected.

A counterion $i$ ($=+,+\!+$) is assumed to bind to the macromolecule and to occupy $f_i$ consecutive (spatially adjacent) charged terminal groups of the macromolecule. We designate $f_{{\scalebox{.8}{$\scriptscriptstyle +$}}}=1$ and $f_{{\scalebox{.8}{$\scriptscriptstyle ++$}}}=2$ for MCs and DCs, respectively, implying that, in a bound state, one MC occupies only one charged terminal group, while one DC occupies two consecutive charged terminal groups, owing to the fact that each terminal group has a charge valency ${z_\mathrm{s}}= -1$. Consider at a given state, ${N^\mathrm{b}_{{\scalebox{.8}{$\scriptscriptstyle +$}}}}$ MCs and ${N^\mathrm{b}_{{\scalebox{.8}{$\scriptscriptstyle ++$}}}}$ DCs are bound to the macromolecule. The *binding density*, *i.e.*, the number of bound counterions per charged terminal group is then ${N^\mathrm{b}_{{\scalebox{.8}{$\scriptscriptstyle +$}}}}/{N_\mathrm{s}}$ and ${N^\mathrm{b}_{{\scalebox{.8}{$\scriptscriptstyle ++$}}}}/{N_\mathrm{s}}$ for MCs and DCs, respectively. By multiplying with $f_i$, we then define the fraction of the binding sites occupied by the counterions, *i.e.*, coverages ${\Theta_{{\scalebox{.8}{$\scriptscriptstyle +$}}}}= f_{{\scalebox{.8}{$\scriptscriptstyle +$}}} {N^\mathrm{b}_{{\scalebox{.8}{$\scriptscriptstyle +$}}}}/ {N_\mathrm{s}}= {N^\mathrm{b}_{{\scalebox{.8}{$\scriptscriptstyle +$}}}}/ {N_\mathrm{s}}$ and ${\Theta_{{\scalebox{.8}{$\scriptscriptstyle ++$}}}}= f_{{\scalebox{.8}{$\scriptscriptstyle ++$}}} {N^\mathrm{b}_{{\scalebox{.8}{$\scriptscriptstyle ++$}}}}/{N_\mathrm{s}}= 2{N^\mathrm{b}_{{\scalebox{.8}{$\scriptscriptstyle ++$}}}}/{N_\mathrm{s}}$. Thus, the total number of binding sites on the macromolecule available for MCs, is $N_{{\scalebox{.8}{$\scriptscriptstyle +$}}} = {N_\mathrm{s}}/f_{{\scalebox{.8}{$\scriptscriptstyle +$}}} = {N_\mathrm{s}}$, and those available for DCs, is $N_{{\scalebox{.8}{$\scriptscriptstyle ++$}}} = {N_\mathrm{s}}/f_{{\scalebox{.8}{$\scriptscriptstyle ++$}}} = {N_\mathrm{s}}/2$. The effective charge valency of the macromolecule is then ${Z_\mathrm{eff}}= -{N_\mathrm{s}}+ {N^\mathrm{b}_{{\scalebox{.8}{$\scriptscriptstyle +$}}}}+ 2{N^\mathrm{b}_{{\scalebox{.8}{$\scriptscriptstyle ++$}}}}= {-{N_\mathrm{s}}(1-{\Theta_{{\scalebox{.8}{$\scriptscriptstyle +$}}}}-{\Theta_{{\scalebox{.8}{$\scriptscriptstyle ++$}}}})}$. The total Helmholtz free energy ${\mathcal{F}_\mathrm{tot}}$ depends on the coverages ${\Theta_{{\scalebox{.8}{$\scriptscriptstyle +$}}}}$ and ${\Theta_{{\scalebox{.8}{$\scriptscriptstyle ++$}}}}$ and the ionic concentrations ${c^0_{i}}$. The coverages can then be obtained by minimizing ${\mathcal{F}_\mathrm{tot}}$ simultaneously with respect to ${\Theta_{{\scalebox{.8}{$\scriptscriptstyle +$}}}}$ and ${\Theta_{{\scalebox{.8}{$\scriptscriptstyle ++$}}}}$. The total Helmholtz free energy ${\mathcal{F}_\mathrm{tot}}$ is given by the expression $${\mathcal{F}_\mathrm{tot}}= {\mathcal{F}_\mathrm{el}}+ {\mathcal{F}_\mathrm{tr}}+ {\mathcal{F}_\mathrm{mix}}+ {\mathcal{F}_\mathrm{int}}\label{eq:mvh_ftot}$$ where the four additive contributions, ${\mathcal{F}_\mathrm{el}}$, ${\mathcal{F}_\mathrm{tr}}$, ${\mathcal{F}_\mathrm{mix}}$ and ${\mathcal{F}_\mathrm{int}}$ are defined respectively as (i) electrostatic (Born) self-energy of charge renormalized macromolecule, (ii) ideal gas entropy of free ions in the bulk regime, (iii) mixing entropy of the condensed counterions in the macromolecule, and (iv) the non-electrostatic ion-specific binding free energy between the condensed counterion and the corresponding binding site on the macromolecule.

The Born charging self-energy of the macromolecule immersed in an electrolyte solution associated with the Debye screening length $\kappa^{-1}$, refers to the work required to charge the macromolecule from its electroneutral to a certain charged state. Following Manning, such a charged state is associated with the effective charge ${Z_\mathrm{eff}}e$, corresponding to the sum of the intrinsic bare charge of the macromolecule ${Z_\mathrm{d}}$ and its captive, neutralizing counterions. [@Manning2007] Thus, the expression for the Born charging free energy of the macromolecule (or the self energy of the charge renormalized macromolecule) per monovalent binding site is thus expressed as $$\beta {\mathcal{F}_\mathrm{el}}= \frac{{Z_\mathrm{eff}}^2 {l_\mathrm{B}}}{2 {N_\mathrm{s}}{r_\mathrm{eff}}(1+\kappa {r_\mathrm{eff}})} = \frac{\zeta}{2} (1-{\Theta_{{\scalebox{.8}{$\scriptscriptstyle +$}}}}-{\Theta_{{\scalebox{.8}{$\scriptscriptstyle ++$}}}})^2
\label{fel1}$$ where $\zeta/2$ is the Born free energy per monovalent binding site in the absence of counterion condensation, and $\zeta$ is given for surface charging by[@mcquarrie2000statistical] $$\zeta =  \frac{{N_\mathrm{s}}{l_\mathrm{B}}}{{r_\mathrm{eff}}(1 + \kappa{r_\mathrm{eff}})}
\label{zeta}$$

Considering the effective volume of dPGS ${v_\mathrm{eff}}$ to be very small compared to the total volume $V$ (${v_\mathrm{eff}}\ll V$), the bulk concentrations of MCs and DCs are given by $$\begin{aligned}
\begin{split}
&{c^{\mathrm{b}}_{{\scalebox{.8}{$\scriptscriptstyle +$}}}}= {c^0_{{\scalebox{.8}{$\scriptscriptstyle +$}}}}+ \frac{{N_{{\scalebox{.8}{$\scriptscriptstyle +$}}}}(1-{\Theta_{{\scalebox{.8}{$\scriptscriptstyle +$}}}})}{V} \\
&{c^{\mathrm{b}}_{{\scalebox{.8}{$\scriptscriptstyle ++$}}}}= {c^0_{{\scalebox{.8}{$\scriptscriptstyle ++$}}}}- \frac{{N_{{\scalebox{.8}{$\scriptscriptstyle ++$}}}}{\Theta_{{\scalebox{.8}{$\scriptscriptstyle ++$}}}}}{V} 
\end{split}
\label{eq:bulk_conc_mvh}\end{aligned}$$ owing to the depletion of the ions in the bulk due to partitioning. ${c^{\mathrm{b}}_{{\scalebox{.8}{$\scriptscriptstyle +$}}}}$ above is calculated considering the monovalent counterions remaining in the solution, in the salt-free limit. We assume that no anions are bound to the macromolecule binding sites, hence their bulk concentration is assumed to be the same as their salt concentration, *i.e.*, ${c^{\mathrm{b}}_{{\scalebox{.8}{$\scriptscriptstyle -$}}}}= {c^0_{{\scalebox{.8}{$\scriptscriptstyle -$}}}}$.

The ideal gas free energy of free cations in the bulk, normalized by the number of monovalent binding sites ${N_\mathrm{s}}$, is given as\
$$\begin{aligned}
\begin{split}
&\beta {\mathcal{F}_\mathrm{tr}}= -\frac{S_\mathrm{id}}{{N_\mathrm{s}}{k_\mathrm{B}}} = \sum_{i=+,++} \left(\frac{n_i - {N^\mathrm{b}_{i}}}{{N_\mathrm{s}}}\right)\left({\mathrm{ln}}\,{c^{\mathrm{b}}_{i}}\Lambda^3_i -1\right) \\
&= \sum_{i=+,++} \left(\frac{n_i - N_i {\Theta_{i}}}{{N_\mathrm{s}}}\right)\left[{\mathrm{ln}}\left({c^0_{i}}\Lambda^3_i - \frac{N_i {\Theta_{i}}\Lambda_i^3}{V}\right) -1\right]
\end{split}\end{aligned}$$ where $\Lambda_i$ and $n_i$ are the thermal (de Broglie) wavelength and the number of salt ions $i$.

The bound DCs and MCs can occupy the binding sites on the macromolecule in different proportions, and can distribute among the occupied sites in multiple ways at a certain bound coverages ${\Theta_{{\scalebox{.8}{$\scriptscriptstyle +$}}}}$ and ${\Theta_{{\scalebox{.8}{$\scriptscriptstyle ++$}}}}$. We exert constraints to such possibilities of binding compositions and configurations, such that, ($i$) one bound DC can only bind to two adjacent monovalent binding sites, ($ii$) all non-overlapping configurations between the bound ions are possible, ($iii$) there are no designated binding sites for DCs, and ($iv$) the position of the bound DC can be shifted by a single adjacent monovalent binding site. The number of possible combinatorial binding arrangements under these constraints, adopted from the work by McGhee and von Hippel, [@Mcghee1974] is given by $$W =  \frac{ \gamma^{{N^\mathrm{b}_{{\scalebox{.8}{$\scriptscriptstyle +$}}}}}_{{\scalebox{.8}{$\scriptscriptstyle +$}}} \gamma^{{N^\mathrm{b}_{{\scalebox{.8}{$\scriptscriptstyle ++$}}}}}_{{\scalebox{.8}{$\scriptscriptstyle ++$}}} ({N_\mathrm{s}}- {N^\mathrm{b}_{{\scalebox{.8}{$\scriptscriptstyle ++$}}}})!}{{N^\mathrm{b}_{{\scalebox{.8}{$\scriptscriptstyle +$}}}}!{N^\mathrm{b}_{{\scalebox{.8}{$\scriptscriptstyle ++$}}}}!({N_\mathrm{s}}- 2{N^\mathrm{b}_{{\scalebox{.8}{$\scriptscriptstyle ++$}}}}- {N^\mathrm{b}_{{\scalebox{.8}{$\scriptscriptstyle +$}}}})!}
\label{w}$$ where we define $\gamma_i={v^0_{i}}/\Lambda^3_i$ in terms of the effective configurational volume ${v^0_{i}}$ in the bound state. [@XiaoCPS] ${v^0_{i}}$ takes into account the rotational and vibrational degrees of freedom of a bound counterion $i$. We now define the free energy associated with the partition function $W$, normalized by the number of monovalent binding sites ${N_\mathrm{s}}$, as the free energy of mixing of the bound ions per binding site, $$\begin{aligned}
\begin{split}
&\beta {\mathcal{F}_\mathrm{mix}}= -\frac{S_\mathrm{mix}}{{N_\mathrm{s}}k_{\rm B}}= -\frac{1}{{N_\mathrm{s}}}{\mathrm{ln}}\, W\\
 & \simeq {\Theta_{{\scalebox{.8}{$\scriptscriptstyle +$}}}}{\mathrm{ln}}{\Theta_{{\scalebox{.8}{$\scriptscriptstyle +$}}}}+ \frac{{\Theta_{{\scalebox{.8}{$\scriptscriptstyle ++$}}}}}{2}{\mathrm{ln}}\frac{{\Theta_{{\scalebox{.8}{$\scriptscriptstyle ++$}}}}}{2} -\left(1-\frac{{\Theta_{{\scalebox{.8}{$\scriptscriptstyle ++$}}}}}{2}\right){\mathrm{ln}}\left(1-\frac{{\Theta_{{\scalebox{.8}{$\scriptscriptstyle ++$}}}}}{2}\right) \\
&+  (1 - {\Theta_{{\scalebox{.8}{$\scriptscriptstyle +$}}}}- {\Theta_{{\scalebox{.8}{$\scriptscriptstyle ++$}}}}){\mathrm{ln}}(1 - {\Theta_{{\scalebox{.8}{$\scriptscriptstyle +$}}}}- {\Theta_{{\scalebox{.8}{$\scriptscriptstyle ++$}}}}) \\
&- {\Theta_{{\scalebox{.8}{$\scriptscriptstyle +$}}}}{\mathrm{ln}}\frac{{v^0_{{\scalebox{.8}{$\scriptscriptstyle +$}}}}}{\Lambda_+^3} - \frac{{\Theta_{{\scalebox{.8}{$\scriptscriptstyle ++$}}}}}{2} {\mathrm{ln}}\frac{{v^0_{{\scalebox{.8}{$\scriptscriptstyle ++$}}}}}{\Lambda_{++}^3}
\end{split}
\label{fmx}\end{aligned}$$ where the Stirling approximation has been used for the logarithm of the factorials. This description of condensed counterion entropy is different than the ion-binding models proposed in previous works for linear polyelectrolytes [@flory1953molecular; @Raphael1990; @Muthukumar2004] in terms of the localization of counterions within volume ${v^0_{i}}$.

We express this intrinsic interaction ${\mathcal{F}_\mathrm{int}}$ by the intrinsic binding chemical potential ${\Delta \mu_{\mathrm{int},\,i}}$ of each bound ion $i$. The sum of such interactions for all bound ions, normalized by the total number of monovalent binding sites gives $$\begin{aligned}
\begin{split}
\beta {\mathcal{F}_\mathrm{int}}&= \frac{1}{{N_\mathrm{s}}}\left({N^\mathrm{b}_{{\scalebox{.8}{$\scriptscriptstyle +$}}}}\beta {\Delta \mu_{\mathrm{int},\,{{\scalebox{.8}{$\scriptscriptstyle +$}}}}}+ {N^\mathrm{b}_{{\scalebox{.8}{$\scriptscriptstyle ++$}}}}\beta {\Delta \mu_{\mathrm{int},\,{{\scalebox{.8}{$\scriptscriptstyle ++$}}}}}\right) \\
&= {\Theta_{{\scalebox{.8}{$\scriptscriptstyle +$}}}}\beta {\Delta \mu_{\mathrm{int},\,{{\scalebox{.8}{$\scriptscriptstyle +$}}}}}+ \frac{{\Theta_{{\scalebox{.8}{$\scriptscriptstyle ++$}}}}}{2}\beta {\Delta \mu_{\mathrm{int},\,{{\scalebox{.8}{$\scriptscriptstyle ++$}}}}}\end{split}\end{aligned}$$

The equilibrium coverages $\Theta_i$ are then obtained by the minimization condition $$\frac{\partial}{\partial {\Theta_{i}}} {\mathcal{F}_\mathrm{tot}}\overset{!}{=} 0 \qquad \quad i=+,+\!+$$ This leads to the relation $${\Delta \mu_{\mathrm{tr},\,i}}+ {\Delta \mu_{\mathrm{el},\,i}}+ {\Delta \mu_{\mathrm{mix},\,i}}+ {\Delta \mu_{\mathrm{int},\,i}}= 0 \quad  \quad i=+,+\!+
\label{k_and_g_all}$$ where ${\Delta \mu_{\mathrm{tr},\,i}}$ denotes the translational entropy change associated with one ion $i$ when it transfers from the bulk environment to the bound state in the macromolecule. ${\Delta \mu_{\mathrm{el},\,i}}$ is the *electrostatic binding chemical potential* and ${\Delta \mu_{\mathrm{mix},\,i}}$ is the *mixing chemical potential*. Eq. , similar to the PPB (Eq. ) and DM (Eq. ) models, indicates the counterion chemical potential components contributing to its condensation on the macromolecule. The expressions for the constituent chemical potential contributions in Eq.  are given by $$\begin{aligned}
\begin{split}
&\beta {\Delta \mu_{\mathrm{tr},\,i}}=-{\mathrm{ln}}\,{c^{\mathrm{b}}_{i}}{v^0_{i}}\qquad \qquad \qquad  \quad i=+,+\!+\\
&\beta {\Delta \mu_{\mathrm{el},\,i}}=-z_i\zeta(1-{\Theta_{{\scalebox{.8}{$\scriptscriptstyle +$}}}}-{\Theta_{{\scalebox{.8}{$\scriptscriptstyle ++$}}}}) \quad  \quad i=+,+\!+\\
&\beta {\Delta \mu_{\mathrm{mix},\,i}}=
\begin{cases}
\begin{aligned}
 {\mathrm{ln}}\frac{{\Theta_{{\scalebox{.8}{$\scriptscriptstyle ++$}}}}\left(2-{\Theta_{{\scalebox{.8}{$\scriptscriptstyle ++$}}}}\right)}{4(1-{\Theta_{{\scalebox{.8}{$\scriptscriptstyle +$}}}}-{\Theta_{{\scalebox{.8}{$\scriptscriptstyle ++$}}}})^2} \quad \quad i=++\\[1.1ex]
 {\mathrm{ln}}\frac{{\Theta_{{\scalebox{.8}{$\scriptscriptstyle +$}}}}}{(1-{\Theta_{{\scalebox{.8}{$\scriptscriptstyle +$}}}}-{\Theta_{{\scalebox{.8}{$\scriptscriptstyle ++$}}}})} \quad  \quad i=+
\end{aligned}
\end{cases} 
\end{split}
\label{all_free_energies}\end{aligned}$$

Using Eqs.  and leads to the final form of the MMvH model, given by $$K_{{\scalebox{.8}{$\scriptscriptstyle ++$}}} = v^0_{{\scalebox{.8}{$\scriptscriptstyle ++$}}} {\mathcal{K}_{\mathrm{int},\,{{\scalebox{.8}{$\scriptscriptstyle ++$}}}}}\mathrm{e}^{2\zeta \left( 1-{\Theta_{{\scalebox{.8}{$\scriptscriptstyle +$}}}}- {\Theta_{{\scalebox{.8}{$\scriptscriptstyle ++$}}}}\right)} = \frac{{\Theta_{{\scalebox{.8}{$\scriptscriptstyle ++$}}}}(2 - {\Theta_{{\scalebox{.8}{$\scriptscriptstyle ++$}}}})}{4 {c^{\mathrm{b}}_{{\scalebox{.8}{$\scriptscriptstyle ++$}}}}{(1 - {\Theta_{{\scalebox{.8}{$\scriptscriptstyle +$}}}}- {\Theta_{{\scalebox{.8}{$\scriptscriptstyle ++$}}}})}^2} 
\label{Kb}$$ $$K_{{\scalebox{.8}{$\scriptscriptstyle +$}}} = v^0_{{\scalebox{.8}{$\scriptscriptstyle +$}}} {\mathcal{K}_{\mathrm{int},\,{{\scalebox{.8}{$\scriptscriptstyle +$}}}}}\mathrm{e}^{\zeta \left(1-{\Theta_{{\scalebox{.8}{$\scriptscriptstyle +$}}}}- {\Theta_{{\scalebox{.8}{$\scriptscriptstyle ++$}}}}\right)} = \frac{ {\Theta_{{\scalebox{.8}{$\scriptscriptstyle +$}}}}}{{c^{\mathrm{b}}_{{\scalebox{.8}{$\scriptscriptstyle +$}}}}(1 - {\Theta_{{\scalebox{.8}{$\scriptscriptstyle +$}}}}- {\Theta_{{\scalebox{.8}{$\scriptscriptstyle ++$}}}})} 
\label{Ka}$$ where $K_i$ are the *equilibrium binding constant* associated with the binding of ion $i$ to its corresponding binding site on the macromolecule. The relationship between $K_i$, the total binding chemical potential ${\Delta \mu_{\mathrm{bind},\,i}}$ and the total partition ratio ${\mathcal{K}_{i}}$ is given as $$\beta {\Delta \mu_{\mathrm{bind},\,i}}= -{\mathrm{ln}}\,\frac{K_i}{{v^0_{i}}} = -{\mathrm{ln}}\,{\mathcal{K}_{i}}\quad  \quad i=+,+\!+
\label{kbind}$$ Or in other words, referring back to Eq. , $${\mathcal{K}_{i}}= {\mathcal{K}_{\mathrm{int},\,{i}}}\,{\mathcal{K}_{\mathrm{el},\,{i}}}= {\mathcal{K}_{\mathrm{int},\,{i}}}\, \mathrm{e}^{z_i \zeta \left(1-{\Theta_{{\scalebox{.8}{$\scriptscriptstyle +$}}}}- {\Theta_{{\scalebox{.8}{$\scriptscriptstyle ++$}}}}\right)}$$ where the electrostatic contribution of the total partition ratio is defined as $${\mathcal{K}_{\mathrm{el},\,{i}}}= \mathrm{e}^{-\beta {\Delta \mu_{\mathrm{el},\,i}}} = \mathrm{e}^{z_i\zeta(1-{\Theta_{{\scalebox{.8}{$\scriptscriptstyle +$}}}}-{\Theta_{{\scalebox{.8}{$\scriptscriptstyle ++$}}}})} \quad  \quad i=+,+\!+
\label{kiel-mvh}$$ From Eq. , for a given magnitude of $K_i$, the absolute magnitude of ${\Delta \mu_{\mathrm{bind},\,i}}$ depends on ${v^0_{i}}$, which we calculate from our simulations and predict respective values of ${\Delta \mu_{\mathrm{bind},\,i}}$.

Finally, we consider the limit of the MMvH model for vanishing DCs (MCs only). Without DCs, we have $$\begin{aligned}
\begin{split}
&\beta {\Delta \mu_\mathrm{tr}}=-{\mathrm{ln}}\,{c^{\mathrm{b}}_{{\scalebox{.8}{$\scriptscriptstyle +$}}}}v^0_{{\scalebox{.8}{$\scriptscriptstyle +$}}}  \\
&\beta {\Delta \mu_\mathrm{el}}=-\zeta(1-{\Theta_{{\scalebox{.8}{$\scriptscriptstyle +$}}}}) \\
&\beta {\Delta \mu_\mathrm{mix}}= {\mathrm{ln}}\,\frac{{\Theta_{{\scalebox{.8}{$\scriptscriptstyle +$}}}}}{(1 - {\Theta_{{\scalebox{.8}{$\scriptscriptstyle +$}}}})}
\label{g_langmuir}
\end{split}\end{aligned}$$ Combining Eqs.  and leads to the “Manning–Langmuir" (ML) model $$K_{{\scalebox{.8}{$\scriptscriptstyle +$}}} = v^0_{{\scalebox{.8}{$\scriptscriptstyle +$}}} {\mathcal{K}_{\mathrm{int},\,{{\scalebox{.8}{$\scriptscriptstyle +$}}}}}\mathrm{e}^{\zeta (1 - {\Theta_{{\scalebox{.8}{$\scriptscriptstyle +$}}}})} = \frac{{\Theta_{{\scalebox{.8}{$\scriptscriptstyle +$}}}}}{{c^{\mathrm{b}}_{{\scalebox{.8}{$\scriptscriptstyle +$}}}}(1 - {\Theta_{{\scalebox{.8}{$\scriptscriptstyle +$}}}})} 
\label{manning-langmuir}$$ The McGhee–von Hippel combinatorics here reduces to the standard one-component Langmuir picture, *i.e.*, the right-hand-side of Eq.  reflects the Langmuir isotherm. The standard Langmuir model is thus extended to include charging free energies by ion condensation (charge renormalization) and ion-specific binding. From another perspective, it extends the Manning model for the counterion condensation on spheres [@Manning2007; @Gillespie2014] to include ion-specific effects as well as the saturation of binding sites in terms of the translation entropy of the condensed ions.

Numerical evaluation {#numerical}
--------------------

The PPB model, with the assumption of the uniform intrinsic macromolecular volume charge distribution ${c_\mathrm{s}}(r)e$ and with the knowledge of the bare radius ${r_\mathrm{d}}$ of the macromolecule inherited from simulations, generates the distance-resolved number density profiles of charged species, similar to Fig. \[density\_ions\]. Hence, it performs the same analysis as that for simulations (*cf.* Section \[sim\_analysis\]), to calculate the effective radius ${r_\mathrm{eff}}$ and other electrostatic properties of the macromolecule, such as ${Z_\mathrm{eff}}$, ${\phi_\mathrm{eff}}$, etc. The DM model also assumes uniform ${c_\mathrm{s}}(r)e$ and requires the knowledge of the electroneutrality radius, which is taken as ${r_\mathrm{d}}$ from simulations as an input parameter, similar to the PPB model. The MMvH (ML) model, on the other hand, assumes the macromolecule as a hard sphere with a uniform surface charge distribution. The effective radius of the hard sphere ${r_\mathrm{eff}}$ is taken from simulations as an input parameter. The results from the DM, PPB and MMvH (ML) models and simulations are compared in terms of the coverages ${\Theta_{i}}$ ($i=+\!+, +$), which are defined as ${\Theta_{i}}={N^\mathrm{b}_{i}}/{N_{i}}$, where ${N^\mathrm{b}_{i}}$ is the number of condensed counterions $i$ and ${N_{i}}$ is the corresponding number of binding sites available on dPGS, defined in Section \[mvh\]. Since the PPB model deals with a volume sorption, while the DM model deals with the ion partitioning between two electroneutral phases, “coverage" ${\Theta_{i}}$ in these cases are interpreted as a load or an extent of neutralization of dPGS. For the DM, PPB and MMvH (ML) models, the intrinsic partition coefficients ${\mathcal{K}_{\mathrm{int},\,{i}}}$ for both ions ($i=+\!+,+$) are unknowns and taken as fitting parameters in order to match the coverages from the simulations, which are described in Section \[sec:cg\_ff\]. Regarding the PPB and DM models, we make a further assumption that intrinsic non-electrostatic ion–binding site interaction for the MCs is identical to that for the monovalent anions, *i.e.*, ${\mathcal{K}_{\mathrm{int},\,{{\scalebox{.8}{$\scriptscriptstyle +$}}}}}= {\mathcal{K}_{\mathrm{int},\,{{\scalebox{.8}{$\scriptscriptstyle -$}}}}}$.

Mathematically, the PPB model represents a boundary-value problem having a second order differential equation (Eq. ) non-linear in the electrostatic potential paired with the boundary conditions, while the MMvH model (Eqs.  and ) represents two non-linear simultaneous equations in coverages ${\Theta_{{\scalebox{.8}{$\scriptscriptstyle +$}}}}$ and ${\Theta_{{\scalebox{.8}{$\scriptscriptstyle ++$}}}}$. Both PPB and MMvH models are evaluated self consistently for the potential and coverages, respectively. To solve Eq. , we employ the `solve_bvp` function in the SciPy library (version 1.3.1) from Python (version 3.7.4), which solves a boundary-value problem for a system of ordinary differential equations using the fourth order collocation algorithm. [@kierzenka2001bvp] The bulk concentration ${c^{\mathrm{b}}_{i}}$ is obtained using the law of conservation of mass, in an iterative manner. Eqs.  and are solved using `fsolve` function from the SciPy library, which is also used to evaluate the DM model (Eq. ) representing the single non-linear equation in the Donnan potential ${\phi^{}_\mathrm{D}}$.

The effective configurational volume ${v^0_{i}}$ of bound counterions, used in the MMvH model is assumed to be equal for both counterions, *i.e.*, $v_{{\scalebox{.8}{$\scriptscriptstyle ++$}}} = v_{{\scalebox{.8}{$\scriptscriptstyle +$}}} = v_0$. It is worth considering that the volume $v_0$ depends on the precise nature of the bound state and it is infeasible to have its knowledge in experiments due to unknown microscopic details, although it can be computed using simulations. [@xu:biomacro; @Yu2015] According to the convention in experiments, the standard volume is defined as $v_0 = 1\,\mathrm{M}^{-1} \simeq 1.6$nm$^{3}$, corresponding to the standard concentration $c^\mathrm{std} = 1$M. [@atkins; @Gilson2007; @General2010a] In this case, the total binding chemical potential ${\Delta \mu_{\mathrm{bind},\,i}}$ can be referred to as the standard binding energy $\Delta G^0$. [@Gilson2007; @General2010a]

Results and Discussion
======================

Monovalent limit: theoretical comparison and best fit to simulations
--------------------------------------------------------------------

![image](images/fig5.pdf){width="12cm" height="12cm"}

Considering the monovalent limit as reference case, we now start with the application of aforementioned theoretical binding models. Fig. \[monocg\](a) shows the predictions of the PPB and ML (monovalent-only limit of MMvH) models for the variation of the binding coverage of MCs, ${\Theta_{{\scalebox{.8}{$\scriptscriptstyle +$}}}}$, as a function of the MC concentration, ${c^0_{{\scalebox{.8}{$\scriptscriptstyle +$}}}}$. It can be observed that ${\Theta_{{\scalebox{.8}{$\scriptscriptstyle +$}}}}$ increases sharply for a small increase in ${c^0_{{\scalebox{.8}{$\scriptscriptstyle +$}}}}$ from $0$ to $\sim 10$mM, while it increases gradually for larger ${c^0_{{\scalebox{.8}{$\scriptscriptstyle +$}}}}$. This is attributed to the combined contribution of the electrostatics and an entropy of a bound counterion, facilitating condensation. In the low ${c^0_{{\scalebox{.8}{$\scriptscriptstyle +$}}}}$ regime, the bare charge of G$_2$-dPGS is weakly renormalized, and some of the dPGS binding sites are unoccupied. This leaves a high propensity of condensation for new incoming counterions. This can be conveniently explained via the ML model. Referring to Eq. , the increase in the condensation of MCs at the limit of low ${c^0_{{\scalebox{.8}{$\scriptscriptstyle +$}}}}$, $\lim_{{c^0_{{\scalebox{.8}{$\scriptscriptstyle +$}}}}\to 0} {{\mathop{}\!\mathrm{d}}}{\Theta_{{\scalebox{.8}{$\scriptscriptstyle +$}}}}/{{\mathop{}\!\mathrm{d}}}{c^0_{{\scalebox{.8}{$\scriptscriptstyle +$}}}}$ is directly proportional to the total binding constant $K_{{\scalebox{.8}{$\scriptscriptstyle +$}}}$, while at high ${c^0_{{\scalebox{.8}{$\scriptscriptstyle +$}}}}$, $\lim_{{c^0_{{\scalebox{.8}{$\scriptscriptstyle +$}}}}\to \infty} {{\mathop{}\!\mathrm{d}}}{\Theta_{{\scalebox{.8}{$\scriptscriptstyle +$}}}}/{{\mathop{}\!\mathrm{d}}}{c^0_{{\scalebox{.8}{$\scriptscriptstyle +$}}}}= 0$. This implies that at low ${c^0_{{\scalebox{.8}{$\scriptscriptstyle +$}}}}$, the resultant low coverage ${\Theta_{{\scalebox{.8}{$\scriptscriptstyle +$}}}}$ leads to a high electrostatic driving force for condensation as well as entropy of a bound counterion, thus a high amount of condensation. On the other hand, at high ${c^0_{{\scalebox{.8}{$\scriptscriptstyle +$}}}}$, the macromolecule charge is almost entirely renormalized and most of the binding sites are occupied, resulting in hardly any increase in condensation.

Comparing the coverage profiles from PPB and ML models that neglect ion-specific effects, *i.e.*, with ${\Delta \mu_{\mathrm{int},\,{{\scalebox{.8}{$\scriptscriptstyle +$}}}}}=0$ (dotted curves), we find that the PPB coverage values are close to the ML values in the low ${c^0_{{\scalebox{.8}{$\scriptscriptstyle +$}}}}$ regime, however, attain higher values than the ML counterpart at high ${c^0_{{\scalebox{.8}{$\scriptscriptstyle +$}}}}$. This is attributed to the effects of discrete binding sites incorporated in the ML model, in the form of the configurational volume $v_0$ (here, we used $v_0=1.04$ M$^{-1}$ obtained from our previous simulations [@xu:biomacro]). The PPB model, on the other hand, assumes the condensed ions as point charges, leaving no entropic penalty for new incoming counterions as they condense on the binding sites. Another reason is that the PPB model also incorporates, to some extent, the non-linear effects in the electrostatic interactions, which are not considered in the DH-level Born energy used in the ML model. Both models, however, underestimate the simulations if we do not include corrections via ${\Delta \mu_{\mathrm{int},\,{{\scalebox{.8}{$\scriptscriptstyle +$}}}}}$. The reason is likely the approximative treatments of the electrostatic energy in both models, PPB and ML, which are mean-field and do not include the discrete nature of the charged binding sites and the complex spatial charge correlations inside the macromolecule. The DM model, in addition to these assumptions, takes the macroscopic view of macromolecule and bulk phases in a segregated form. The model then predicts the ion partitioning while imposing electroneutralities of phases. In that respect, for highly charged macromolecules like dPGS, the DM model predicts ${N^\mathrm{b}_{{\scalebox{.8}{$\scriptscriptstyle +$}}}}\simeq {N_\mathrm{s}}$, implying ${\Theta_{{\scalebox{.8}{$\scriptscriptstyle +$}}}}\simeq 1$. This plot is not shown, since it does not provide a useful insight for us in the context of counterion condensation. The case of salt concentration ${c^0_{{\scalebox{.8}{$\scriptscriptstyle +$}}}}= 0$ is referred to as the counterion-only case, and gives ${\Theta_{{\scalebox{.8}{$\scriptscriptstyle +$}}}}\sim 0.28$ for the PPB model. Note that ${\Theta_{{\scalebox{.8}{$\scriptscriptstyle +$}}}}$ in this limit is system specific, since the size of the simulation box/computational domain determines the counterion concentration and subsequently the coverage. The coverage ${\Theta_{{\scalebox{.8}{$\scriptscriptstyle +$}}}}$ in the ML model in this limit is undefined, since the electrostatic binding energy of MCs depends on the screening length $\kappa^{-1}$, which is undefined in this model in the absence of the salt.

In the next step, ${\Theta_{{\scalebox{.8}{$\scriptscriptstyle +$}}}}$ values for PPB and ML models are fitted (bold curves in Fig. \[monocg\](a)) to the simulation results for G$_2$-dPGS in the monovalent limit by allowing ion-specific effects in the counterion–macromolecule binding, *i.e.*, ${\Delta \mu_{\mathrm{int},\,{+}}}$ as a fitting parameter. The values of ${\Delta \mu_{\mathrm{int},\,{+}}}$ are found to be $-0.45\,{k_\mathrm{B}T}$ and $-1.81\,{k_\mathrm{B}T}$ for PPB and ML models, respectively. Recall that the simulations have not really included ion-specific effects in terms of specific hydration phenomena, etc., still, they include excluded-volume, dispersion attraction, and importantly, all electrostatic charge–charge correlations, not captured in the mean-field theories. Hence, the ion-specific fitting parameters can be viewed in general as correction factors, including all ionic contributions that are beyond the mean-field treatment of the PPB and ML models. The larger fitting parameter for ML than PPB (in the absolute value) may indicate the higher level of approximations in the ML model. Having the models now informed using the benchmark data from simulations, they can be utilized to predict the binding at other ion concentrations.

Fig. \[monocg\](b) shows the numerical fitting of ${\Theta_{{\scalebox{.8}{$\scriptscriptstyle +$}}}}$ values (dashed curves) to those obtained from G$_4$-dPGS simulations. The values of ${\Delta \mu_{\mathrm{int},\,{+}}}$ as a fitting parameter are $-0.56\,{k_\mathrm{B}T}$ and $-1.85\,{k_\mathrm{B}T}$ for PPB and ML models, respectively, which are close to those obtained for G$_2$-dPGS, within the error difference of $\sim 0.1\,{k_\mathrm{B}T}$. The ML model fits better to both G$_2$-dPGS and G$_4$-dPGS CG results than the PPB model at large ${c^0_{{\scalebox{.8}{$\scriptscriptstyle +$}}}}$, which may indicate that the dPGS charge in the simulations acts more as finite binding sites, as assumed in the ML model.

Divalent case: theoretical comparison and best fit to simulations
-----------------------------------------------------------------

![image](images/fig6.pdf){width="14cm" height="14cm"}

We now aspire to use the obtained ${\Delta \mu_{\mathrm{int},\,{{\scalebox{.8}{$\scriptscriptstyle +$}}}}}$ to inform the MMvH and PPB models with the help of the reference data obtained from simulations, in order to capture the competitive ion binding in a mixture of MCs and DCs. The models fitted to the benchmark data can then be used to predict the binding coverages ${\Theta_{{\scalebox{.8}{$\scriptscriptstyle ++$}}}}$ and ${\Theta_{{\scalebox{.8}{$\scriptscriptstyle +$}}}}$ for different dPGS generations and salt concentrations. In practice, we perform the numerical fitting of ${\Theta_{{\scalebox{.8}{$\scriptscriptstyle ++$}}}}$ and ${\Theta_{{\scalebox{.8}{$\scriptscriptstyle +$}}}}$ obtained from the MMvH and PPB models to those from simulations, by fixing ${\Delta \mu_{\mathrm{int},\,{{\scalebox{.8}{$\scriptscriptstyle +$}}}}}$ for MCs obtained from the monovalent-only case, and then subsequently fitting ${\Delta \mu_{\mathrm{int},\,{{\scalebox{.8}{$\scriptscriptstyle ++$}}}}}$ for DCs. The values of ${\Delta \mu_{\mathrm{int},\,{{\scalebox{.8}{$\scriptscriptstyle +$}}}}}$ for MCs obtained from the monovalent limit are, for a given binding model (ML or PPB), found to be approximately independent of the dPGS generation (with $\sim 0.1\,{k_\mathrm{B}T}$ as margin of error). Therefore, ${\Delta \mu_{\mathrm{int},\,{{\scalebox{.8}{$\scriptscriptstyle +$}}}}}$ is averaged over generations (G$_2$ and G$_4$), as shown in Table \[mvh\_ppb\_kint\_table\]. Fig. \[divalent\] depicts the behavior of MMvH, PPB and the DM model in terms of the binding coverages ${\Theta_{i}}$, in a mixture of DCs and MCs. The MMvH model uses the effective configurational volumes $v_0=1.04$ M$^{-1}$ and $0.57$ M$^{-1}$ for G$_2$-dPGS and G$_4$-dPGS, respectively, as obtained from our previous simulations. [@xu:biomacro] At low DC concentration, *i.e.* in the monovalent limit, MCs act as the only counterions to the macromolecule, resulting in the highest MC coverage ${\Theta_{{\scalebox{.8}{$\scriptscriptstyle +$}}}}$. In this limit at ${c^0_{{\scalebox{.8}{$\scriptscriptstyle +$}}}}= 150.37$mM, both MMvH and PPB models show ${\Theta_{{\scalebox{.8}{$\scriptscriptstyle +$}}}}\simeq 0.57$ for G$_2$-dPGS, and ${\Theta_{{\scalebox{.8}{$\scriptscriptstyle +$}}}}\simeq 0.8$ for G$_4$-dPGS. As ${c^0_{{\scalebox{.8}{$\scriptscriptstyle ++$}}}}$ increases, more DCs bind to the macromolecule and more of the previously bound MCs get released into the bulk. Table \[cg\_table\] shows the resultant effective charge valency ${Z_\mathrm{eff}}^\mathrm{PB}$ and potential ${\phi_\mathrm{eff}}^\mathrm{PB}$ of G$_2$-dPGS and G$_4$-dPGS evaluated by the PPB model. Quantitatively consistent with the ${Z_\mathrm{eff}}$ and ${\phi_\mathrm{eff}}$ obtained from simulations, ${Z_\mathrm{eff}}^\mathrm{PB}$ and ${\phi_\mathrm{eff}}^\mathrm{PB}$ show a strong decrease in magnitude with a higher ${c^0_{{\scalebox{.8}{$\scriptscriptstyle ++$}}}}$, depicting higher dPGS charge renormalization.

[cdd]{} &\
& ++ & +\
& 5.13 & 3.37\
& -1.87 & -0.50\
($v_0$ [CG]{}) & -2.85 & -1.83\
($v_0$ [Std.]{}) & -2.86 & -1.44\

![image](images/fig7.pdf){width="11.8cm" height="11.8cm"}

Corresponding to the fitting of binding coverages ${\Theta_{i}}$ on G$_2$-dPGS and G$_4$-dPGS binding sites, as shown in Fig. \[divalent\], the resulting ${\Delta \mu_{\mathrm{int},\,{{\scalebox{.8}{$\scriptscriptstyle ++$}}}}}$ values are calculated as $-2.73\,{k_\mathrm{B}T}$ (G$_2$) and $-2.98\,{k_\mathrm{B}T}$ (G$_4$) for the MMvH model, whereas $-1.77\,{k_\mathrm{B}T}$ (G$_2$) and $-1.98\,{k_\mathrm{B}T}$ (G$_4$) for the PPB model. Table \[mvh\_ppb\_kint\_table\] shows the values of ${\Delta \mu_{\mathrm{int},\,{{\scalebox{.8}{$\scriptscriptstyle ++$}}}}}$ averaged over G$_2$-dPGS and G$_4$-dPGS cases. It can be observed that both ${\Delta \mu_{\mathrm{int},\,{{\scalebox{.8}{$\scriptscriptstyle ++$}}}}}$ and ${\Delta \mu_{\mathrm{int},\,{{\scalebox{.8}{$\scriptscriptstyle +$}}}}}$ values from the MMvH model exceed (in magnitude) those from the PPB model across the whole ${c^0_{{\scalebox{.8}{$\scriptscriptstyle ++$}}}}\sim 0-25$mM range. This can again be attributed to higher approximations in the electrostatic partition coefficient designed in the MMvH model, based on the Debye–Hückel charging free energy, as compared to that from the PPB model, incorporating non-linear effects in the electrostatic potential near the macromolecule vicinity. The standard intrinsic chemical potentials ${\Delta \mu_{\mathrm{int},\,i}}^0$ after fitting the MMvH model ${\Theta_{i}}$ with those from simulations are also given in Table \[mvh\_ppb\_kint\_table\].

Unlike the other models, we simultaneously fit both ${\Delta \mu_{\mathrm{int},\,{{\scalebox{.8}{$\scriptscriptstyle +$}}}}}$ and ${\Delta \mu_{\mathrm{int},\,{{\scalebox{.8}{$\scriptscriptstyle ++$}}}}}$ to perform numerical fitting of ${\Theta_{{\scalebox{.8}{$\scriptscriptstyle +$}}}}$ and ${\Theta_{{\scalebox{.8}{$\scriptscriptstyle ++$}}}}$ obtained from the DM model with the simulation data. As shown in Table \[mvh\_ppb\_kint\_table\], the values of ${\Delta \mu_{\mathrm{int},\,{{\scalebox{.8}{$\scriptscriptstyle +$}}}}}$ and ${\Delta \mu_{\mathrm{int},\,{{\scalebox{.8}{$\scriptscriptstyle ++$}}}}}$ for the model turn out large and positive compared with those from other models, since the DM model tries to neutralize the entire dPGS charge via the electroneutrality condition in the dPGS phase. The DM fits for ${\Theta_{{\scalebox{.8}{$\scriptscriptstyle +$}}}}$ differ to an extent with those from simulations, while those for ${\Theta_{{\scalebox{.8}{$\scriptscriptstyle ++$}}}}$ are found to be reasonably good. The DM provides much better fits for ${\Theta_{{\scalebox{.8}{$\scriptscriptstyle +$}}}}$ in the case of G$_4$-dPGS as compared to G$_2$-dPGS. This is attributed to the bigger size of G$_4$-dPGS, which better satisfies the criterion $\kappa {r_\mathrm{d}}\gg 1$, under which the DM electroneutrality condition holds comparatively well.

Having established the model frameworks by informing ${\Delta \mu_{\mathrm{int},\,i}}$ by fitting the coverages ${\Theta_{i}}$ to those from simulations and averaging the values of obtained ${\Delta \mu_{\mathrm{int},\,i}}$ over generations (See Table \[mvh\_ppb\_kint\_table\]), we finally utilize their predictive ability to explore the electrostatic characterization of dPGS for different generations and salt concentrations. As an example, Fig. \[prediction\](d) shows the MMvH model predictions for the binding coverages ${\Theta_{{\scalebox{.8}{$\scriptscriptstyle +$}}}}$ and ${\Theta_{{\scalebox{.8}{$\scriptscriptstyle ++$}}}}$ for the case of a competitive ion binding on G$_6$-dPGS, similar to Fig. \[divalent\] on G$_2$-dPGS and G$_4$-dPGS. We also study the effective charge valency ${Z_\mathrm{eff}}$ of dPGS along with the composition of condensed ions on the molecule. Figs. \[prediction\](a) and \[prediction\](b) show the variation of the effective charge valency ${Z_\mathrm{eff}}$ of G$_2$-dPGS and its normalized form ${Z_\mathrm{eff}}/{Z_\mathrm{d}}$, respectively, as a function of the DC concentration ${c^0_{{\scalebox{.8}{$\scriptscriptstyle ++$}}}}$, as predicted by the MMvH model. It can be clearly seen from Fig. \[prediction\](a) that the introduction of DCs leads to a net charge renormalization of dPGS, which further decreases its ${Z_\mathrm{eff}}$. The inset shows that, with reference to the monovalent limit, the dPGS effective charge is $30-35\%$ further renormalized upon introducing DCs in the range of $1-4$mM, which is close to the physiological concentration range for calcium(II) ions. Fig. \[prediction\](b) shows that the fraction of the bare dPGS charge that gets renormalized increases with the dPGS generation. The inset shows the variation for ${c^0_{{\scalebox{.8}{$\scriptscriptstyle ++$}}}}$ varying from $0$mM to $10$mM. The rate of dPGS charge renormalization with respect to ${c^0_{{\scalebox{.8}{$\scriptscriptstyle ++$}}}}$ is the highest at the low ${c^0_{{\scalebox{.8}{$\scriptscriptstyle ++$}}}}$ regime and subsides as ${c^0_{{\scalebox{.8}{$\scriptscriptstyle ++$}}}}$ increases, since the charge renormalized dPGS results in lower electrostatic binding chemical potential ${\Delta \mu_{\mathrm{el},\,i}}$. The reduced amount of renormalization is not attributed to the ion packing, which is evident from Fig. \[prediction\](c) showing the total number of condensed ions (including both DCs and MCs) per dPGS sulfate group. As ${c^0_{{\scalebox{.8}{$\scriptscriptstyle ++$}}}}$ increases, the total number of condensed ions decreases, indicating that the ion packing effects diminish as ${c^0_{{\scalebox{.8}{$\scriptscriptstyle ++$}}}}$ increases. The decrease in the amount of renormalization thus predominantly has electrostatic origin. Fig. \[prediction\](a) shows that $80-90\%$ of the dPGS bare charge is renormalized as ${c^0_{{\scalebox{.8}{$\scriptscriptstyle ++$}}}}$ increases from $0-100$mM, however, the total number of condensed counterions effectively decreases, according to Fig. \[prediction\](c). This in effect would significantly hamper the binding affinity of protein with dPGS. It has been well established through our previous works that the dPGS–protein complexation is dominantly influenced by the release of a few MCs that were highly confined due to strong charge renormalization. [@xu:biomacro] The introduction of DCs, however, decreases the confinement of these condensed counterions, thus less counterions to be released during dPGS-protein binding. In addition, the strongly charge renormalized dPGS leads to lower electrostatic contribution to its overall binding affinity with the protein or any other multivalent ligand.

Conclusion
==========

In this paper, we have addressed the biologically and industrially relevant problem of the competitive sorption of mono- and divalent counterions into a highly charged globular polyelectrolyte, with direct comparison to CG simulations of the dendritic macromolecule dPGS. Beyond simple Donnan and ion-specific penetrable PB models, we introduced a two-state discrete binding site model (MMvH) applicable for heterogeneous ligand systems (counterions with mixed valencies/stoichiometries). The broad classification of surrounding counterions as “bound" and “free" gives the MMvH model a computationally unique advantage over the PPB model, which involves the calculation of the distance-resolved counterion density profiles. The fitting results with simulations highlight the key differences in the MMvH and PPB models. Although being on a mean-field level, the PPB model incorporates non-linear electrostatic effects, which become more prominent near the surface of dPGS, delivering a relatively accurate picture of the dPGS–counterion electrostatic binding affinity, compared to the MMvH model, which approximates dPGS–counterion electrostatic interaction on a linearized PB (DH) level by absorbing these non-linear electrostatic effects into the effective charge valency ${Z_\mathrm{eff}}$ of dPGS. On the contrary, the MMvH model provides more accurate values of the extent of counterion adsorption $\Theta$ at high concentrations (*i.e.*, in the binding site saturation regime) than the PPB model. The reason is that the MMvH model assumes discrete binding sites, whereas the PPB model treats dPGS charge as continuum and allows an unlimited uptake of counterions, which is not realistic.

Future extensions of the MMvH model could include an extra level of competition between adsorbed ions explicitly, namely through a non-linear term in Eq.  (of the type used in the regular solution theory or the Flory–Huggins approximation in polymer theories) that describes the interaction between two adsorbed ions in proximal positions (sites). The effects of this generalization in a different context can be found in a study on ion induced lamellar-lamellar phase transition in charged surfactant systems. [@Harries2006] In general, this type of competition results in non-continuous adsorption equilibria and could be interesting in the present context.

The simplest presented model, the Donnan model (DM) extended for ion-specific effects, is also useful for a quick, qualitative prediction of the adsorption ratio. Per construction it should become more accurate for large globules and/or large salt concentrations (for which the globule size becomes larger than the DH screening length), where the electroneutrality condition is better justified.

The models presented in this work can be used to accurately extrapolate and predict the competitive ionic sorption in experiments for a wide range of salt concentrations and salt compositions. They can be also easily generalized to more ionic components and valencies. The electroneutrality radius required for the DM model and the intrinsic macromolecular charge distribution required for the PPB model as an input parameter (in the form of the bare radius ${r_\mathrm{d}}$), are taken from simulations. However, they can also be derived by measuring the form factors from, *e.g.*, neutron scattering. [@Boris1996; @Berndt2006] The MMvH (ML) model requires the effective radius ${r_\mathrm{eff}}$ of the macromolecule as an input parameter, which besides simulations, can also be derived from independent experiments such as electrophoresis and fitting structure factors (of non-dilute colloidal suspensions) by DLVO interactions. [@hunter; @israelachvili2011intermolecular] As we showed, ${r_\mathrm{eff}}$ can also be obtained using PB models and related theories provided the intrinsic macromolecular charge distribution is available.

The authors are indebted to Matthias Ballauff for insightful discussions. R.N. thanks Jacek Walkowiak for fruitful discussion. This project has received funding from the European Research Council (ERC) under the European Union’s Horizon 2020 research and innovation programme (grant agreement No. 646659). X. X. acknowledges the National Science Foundation of China (21903045) and China Postdoctoral Science Foundation (2019M661842) for financial support. M.K. acknowledges the financial support from the Slovenian Research Agency (research core funding No. P1-0055).

Data availability {#data-availability .unnumbered}
=================

The data that support the findings of this study are available from the corresponding author upon reasonable request.

[100]{}

M. Muthukumar, Macromolecules [**50**]{}, 9528 (2017).

A. Katchalsky, Biophys. J. [**4**]{}, 9 (1964).

M. Rubinstein and G. A. Papoian, Soft Matter [**8**]{}, 9265 (2012).

H. Boroudjerdi et al., Phys. Rep. [**416**]{}, 129 (2005).

S. Förster and M. Schmidt, Adv. Polym. Sci. [**120**]{} (1995).

A. Dobrynin and M. Rubinstein, Prog. Polym. Sci. [**30**]{}, 1049 (2005).

S. Liu, K. Ghosh, and M. Muthukumar, J. Chem. Phys. [**119**]{}, 1813 (2003).

A. Chremos and J. F. Douglas, Soft Matter [**12**]{}, 2932 (2016).

S. Alexander et al., J. Chem. Phys. [**80**]{}, 5776 (1984).

L. Belloni, Colloids Surf., A [**140**]{}, 227 (1998).

L. Bocquet, E. Trizac, and M. Aubouy, J. Chem. Phys. [**117**]{}, 8138 (2002).

S. F[ö]{}rster, M. Schmidt, and M. Antonietti, J. Phys. Chem. [**96**]{}, 4008 (1992).

A. Dobrynin, R. H. Colby, and M. Rubinstein, Macromolecules [**28**]{}, 1859 (1995).

J. R. Wenner, M. C. Williams, I. Rouzina, and V. A. Bloomfield, Biophys. J. [**82**]{}, 3160 (2002).

E. Raspaud, M. [Olvera de la Cruz]{}, J.-L. Sikorav, and F. Livolant, Biophys. J. [**74**]{}, 381 (1998).

S. Liu and M. Muthukumar, J. Chem. Phys. [**116**]{}, 9975 (2002).

M. Muthukumar, J. Chem. Phys. [**120**]{}, 9343 (2004).

J. R[ü]{}he et al., Polyelectrolyte brushes, in [*Polyelectrolytes with Defined Molecular Architecture I*]{}, pages 79–150, Springer, 2004.

P. Pincus, Macromolecules [**24**]{}, 2912 (1991).

O. Borisov, T. Birshtein, and E. Zhulina, J. Phys. II (Paris) [**1**]{}, 521 (1991).

E. B. Zhulina, T. M. Birshtein, and O. V. Borisov, Macromolecules [**28**]{}, 1491 (1995).

E. B. Zhulina, J. [Klein Wolterink]{}, and O. V. Borisov, Macromolecules [**33**]{}, 4945 (2000).

P. K. Jha, J. W. Zwanikken, F. A. Detcheverry, J. J. de Pablo, and M. Olvera de la Cruz, Soft Matter [**7**]{}, 5965 (2011).

J. Hua, M. K. Mitra, and M. Muthukumar, The Journal of Chemical Physics [**136**]{}, 134901 (2012).

J. Landsgesell, D. Sean, P. Kreissl, K. Szuttor, and C. Holm, Phys. Rev. Lett. [**122**]{}, 208002 (2019).

I. Adroher-Ben[í]{}tez et al., Phys. Chem. Chem. Phys. [**19**]{}, 6838 (2017).

M. Quesada-P[é]{}rez, J. A. Maroto-Centeno, A. Mart[í]{}n-Molina, and A. Moncho-Jord[á]{}, Phys. Rev. E [**97**]{}, 042608 (2018).

S. F[ö]{}rster and M. Schmidt, Polyelectrolytes in solution, in [*Physical Properties of Polymers*]{}, volume 120, Springer, Berlin, Heidelberg, 1995.

G. S. Manning, Soft Matter [**8**]{}, 9334 (2012).

H. G. M. [Van de Steeg]{}, M. A. [Cohen Stuart]{}, A. [De Keizer]{}, and B. H. Bijsterbosch, Langmuir [**8**]{}, 2538 (1992).

M. A. G. Dahlgren et al., J. Phys. Chem. [**97**]{}, 11769 (1993).

R. R. Netz and J. F. Joanny, Macromolecules [**32**]{}, 9026 (1999).

R. Hariharan, C. Biver, J. Mays, and W. B. Russel, Macromolecules [**31**]{}, 7506 (1998).

D. I. Gittins and F. Caruso, J. Phys. Chem. B [**105**]{}, 6846 (2001).

F. Caruso, Chem. Eur. J. [**6**]{}, 413 (2000).

G. Decher, J. D. Hong, and J. Schmitt, Thin Solid Films [**210**]{}, 831 (1992).

G. Decher, Science [**277**]{}, 1232 (1997).

G. Ladam et al., Langmuir [**16**]{}, 1249 (2000).

R. A. McAloney, M. Sinyor, V. Dudnik, and M. C. Goh, Langmuir [**17**]{}, 6655 (2001).

S. T. Dubas and J. B. Schlenoff, Langmuir [**17**]{}, 7725 (2001).

A. Hugerth, N. Caram-Lelham, and L.-O. Sundel[ö]{}f, Carbohydr. Polym. [**34**]{}, 149 (1997).

L. Rusu-Balaita, J. Desbrieres, and M. Rinaudo, Polym. Bull. [**50**]{}, 91 (2003).

R. G. Winkler, M. O. Steinhauser, and P. Reineker, Phys. Rev. E [**66**]{}, 021802 (2002).

A. Kudlay and M. [Olvera de la Cruz]{}, J. Chem. Phys. [**120**]{}, 404 (2004).

M. Mende, G. Petzold, and H.-M. Buchhammer, Colloid Polym. Sci. [**280**]{}, 342 (2002).

E. Spruijt, A. H. Westphal, J. W. Borst, M. A. Cohen Stuart, and J. van der Gucht, Macromolecules [**43**]{}, 6476 (2010).

J. van der Gucht, E. Spruijt, M. Lemmers, and M. A. [Cohen Stuart]{}, J. Colloid Interface Sci. [**361**]{}, 407 (2011).

P. M. Biesheuvel and M. A. Cohen Stuart, Langmuir [**20**]{}, 2785 (2004).

S. Perry, Y. Li, D. Priftis, L. Leon, and M. Tirrell, Polymers (Basel). [**6**]{}, 1756 (2014).

J. Dernedde et al., Proc. Natl. Acad. Sci. U. S. A. [**107**]{}, 19679 (2010).

J. Khandare, M. Calder[ó]{}n, N. M. Dagiaa, and R. Haag, Chem. Soc. Rev. [**41**]{}, 2824 (2012).

D. Gr[ö]{}ger et al., Bioconjug. Chem. [**24**]{}, 1507 (2013).

D. Maysinger et al., Biomacromolecules [**16**]{}, 3073 (2015).

S. Reimann et al., Adv. Healthc. Mater. [**4**]{}, 2154 (2015).

C. C. Lee, J. A. MacKay, J. M. Fr[é]{}chet, and F. C. Szoka, Nat. Biotechnol. [**23**]{}, 1517 (2005).

M. Ballauff and C. N. Likos, Angew. Chem., Int. Ed. [**43**]{}, 2998 (2004).

W. Tian and Y. Ma, Chem. Soc. Rev. [**42**]{}, 705 (2013).

A. Sousa-Herves et al., Nanoscale [**7**]{}, 3923 (2015).

J. Vonnemann et al., Theranostics [**4**]{}, 629 (2014).

X. Xu, Q. Ran, R. Haag, M. Ballauff, and J. Dzubiella, Macromolecules [**50**]{}, 4759 (2017).

R. Nikam, X. Xu, M. Ballauff, M. Kandu[č]{}, and J. Dzubiella, Soft Matter [**14**]{}, 4300 (2018).

H. Ohshima, T. W. Healy, and L. R. White, J. Colloid Interface Sci. [**90**]{}, 17 (1982).

B. H. Zimm and M. L. Bret, J. Biomol. Struct. Dyn. [**1**]{}, 461 (1983).

L. Belloni, M. Drifford, and P. Turq, Chem. Phys. [**83**]{}, 147 (1984).

G. V. Ramanathan, J. Chem. Phys. [**88**]{}, 3887 (1988).

G. S. Manning, J. Phys. Chem. B [**111**]{}, 8554 (2007).

D. A. J. Gillespie et al., Soft Matter [**10**]{}, 566 (2014).

H. Ohshima, J. Colloid Interface Sci. [**323**]{}, 92 (2008).

E. Trizac, L. Bocquet, and M. Aubouy, Phys. Rev. Lett. [**89**]{}, 248301 (2002).

L. Bocquet, E. Trizac, and M. Aubouy, J. Chem. Phys. [**117**]{}, 8138 (2002).

A. Diehl and Y. Levin, J. Chem. Phys. [**121**]{}, 12100 (2004).

X. Xu et al., Biomacromolecules [**19**]{}, 409 (2018).

D. E. Owens and N. A. Peppas, Int. J. Pharm. [**307**]{}, 93 (2006).

T. Cedervall et al., Proc. Natl. Acad. Sci. [**104**]{}, 2050 (2007).

S. Lindman et al., Nano Lett. [**7**]{}, 914 (2007).

M. P. Monopoli, C. [Å]{}berg, A. Salvati, and K. A. Dawson, Nat. Nanotechnol. [**7**]{}, 779 (2012).

F. Wang et al., Nanomedicine: Nanotechnology, Biology and Medicine [**9**]{}, 1159 (2013).

M. C. [Lo Giudice]{}, L. M. Herda, E. Polo, and K. A. Dawson, Nat. Commun. [**7**]{}, 13475 (2016).

L. Boselli, E. Polo, V. Castagnola, and K. A. Dawson, Angew. Chemie Int. Ed. [**56**]{}, 4215 (2017).

S. Friesen, G. Hefter, and R. Buchner, J. Phys. Chem. B [**123**]{}, 891 (2019).

J. F. Da Silva and R. J. P. Williams, , Oxford University Press, 2001.

R. A. Meyers and R. A. Meyers, , Wiley Online Library, 2004.

R. H. Kretsinger, V. N. Uversky, and E. A. Permyakov, , Springer, 2013.

X. Xu and J. Dzubiella, Colloid Polym. Sci. , 17 (2020).

I. Pochard, A. Foissy, and P. Couchot, Colloid Polym. Sci. [**277**]{}, 818 (1999).

L. Birnhack, O. Keller, S. C. Tang, N. Fridman-Bishop, and O. Lahav, Sep. Purif. Technol. [**223**]{}, 24 (2019).

Y. Kobuchi, H. Motomura, Y. Noma, and F. Hanada, J. Membr. Sci. [**27**]{}, 173 (1986).

T. Sata, New applications of ion exchange membranes, in [*Macromolecules 1992*]{}, page 451, VSP Utrecht, 1993.

T. Sata, T. Sata, and W. Yang, J. Memb. Sci. [**206**]{}, 31 (2002).

P. S. Kuhn, Y. Levin, and M. C. Barbosa, Phys. A Stat. Mech. its Appl. [**266**]{}, 413 (1999).

S. Gavryushov and P. Zielenkiewicz, J. Phys. Chem. B [**101**]{}, 792 (1997).

A. P. dos Santos, A. Diehl, and Y. Levin, J. Chem. Phys. [**132**]{}, 104105 (2010).

J. Arenzon, J. Stilck, and Y. Levin, Eur. Phys. J. B [**12**]{}, 79 (1999).

A. Naji, A. Arnold, C. Holm, and R. R. Netz, Europhys. Lett. [**67**]{}, 130 (2004).

M. Kandu[č]{}, A. Naji, and R. Podgornik, J. Chem. Phys. [**132**]{}, 224703 (2010).

N. Ad[ž]{}i[ć]{} and R. Podgornik, J. Chem. Phys. [**144**]{}, 214901 (2016).

P. J. Basser and A. J. Grodzinsky, Biophys. Chem. [**46**]{}, 57 (1993).

A. Moncho-Jord[á]{} and J. Dzubiella, Phys. Chem. Chem. Phys. [**18**]{}, 5372 (2016).

S. Ahualli, A. Mart[í]{}n-Molina, and M. Quesada-P[é]{}rez, Phys. Chem. Chem. Phys. [**16**]{}, 25483 (2014).

J. N. Israelachvili, , Academic Press London; San Diego, 2nd ed. edition, 1991.

A. W. Adamson and A. P. Gast, , volume 15, Interscience New York, 1967.

E. J. W. Verwey, J. Phys. Chem. [**51**]{}, 631 (1947).

I. Borukhov, D. Andelman, and H. Orland, Electrochim. Acta [**46**]{}, 221 (2000).

M. Eigen and E. Wicke, J. Phys. Chem. [**58**]{}, 702 (1954).

V. Kralj-Iglič and A. Iglič, J. Phys. II [**6**]{}, 477 (1996).

N. Cuvillier et al., Prog. Colloid Polym. Sci. [**105**]{}, 118 (1997).

I. Borukhov, D. Andelman, and H. Orland, Phys. Rev. Lett. [**79**]{}, 435 (1997).

I. Kalcher, J. C. F. Schulz, and J. Dzubiella, Phys. Rev. Lett. [**104**]{}, 097802 (2010).

I. Kalcher, J. C. F. Schulz, and J. Dzubiella, J. Chem. Phys. [**133**]{}, 164511 (2010).

R. Chudoba, J. Heyda, and J. Dzubiella, Soft Matter [**14**]{}, 9631 (2018).

A. Moncho-Jord[á]{} and I. Adroher-Ben[í]{}tez, Soft Matter [**10**]{}, 5810 (2014).

D. Ben-Yaakov, D. Andelman, D. Harries, and R. Podgornik, J. Phys. Condens. Matter [**21**]{} (2009).

E. R. A. Lima et al., J. Phys. Chem. B [**112**]{}, 1580 (2008).

P. Koelsch et al., Colloids Surfaces A Physicochem. Eng. Asp. [**303**]{}, 110 (2007).

H. I. Okur et al., J. Phys. Chem. B [**121**]{}, 1997 (2017).

P. [Lo Nostro]{} and B. W. Ninham, Chem. Rev. [**112**]{}, 2286 (2012).

N. Schwierz, D. Horinek, and R. R. Netz, Langmuir [**29**]{}, 2602 (2013).

Y. Burak, G. Ariel, and D. Andelman, Curr. Opin. Colloid Interface Sci. [**9**]{}, 53 (2004).

S.-w. W. Chen and B. Honig, J. Phys. Chem. B [**101**]{}, 9113 (1997).

I. Rouzina and V. A. Bloomfield, Biophys. Chem. [**64**]{}, 139 (1997).

V. K. Misra, K. A. Sharp, R. A. Friedman, and B. Honig, J. Mol. Biol. [**238**]{}, 245 (1994).

M. D. Paulsen, B. Richey, C. F. Anderson, and M. Record, Chem. Phys. Lett. [**139**]{}, 448 (1987).

I. Rouzina and V. A. Bloomfield, J. Phys. Chem. [**100**]{}, 4292 (1996).

G. S. Manning, J. Chem. Phys. [**51**]{}, 924 (1969).

F. Oosawa, , Marcel Dekker, New York, 1971.

J. D. McGhee and P. H. von Hippel, J. Mol. Biol. [**86**]{}, 469 (1974).

P. Flory, J. Chem. Phys. [**21**]{}, 162 (1953).

E. Raphael and J. F. Joanny, EPL [**13**]{}, 623 (1990).

R. A. G. Friedman and G. S. Manning, Biopolymers [**23**]{}, 2671 (1984).

A. Kundagrami and M. Muthukumar, J. Chem. Phys. [**128**]{}, 244901 (2008).

S. Hui[ß]{}mann, A. Wynveen, C. N. Likos, and R. Blaak, J. Phys.: Condens. Matter [**22**]{}, 232101 (2010).

S. Hui[ß]{}mann, A. Wynveen, C. N. Likos, and R. Blaak, J. Mat. Chem. [**20**]{}, 10486 (2010).

J. S. K[ł]{}os and J.-U. Sommer, Macromolecules [**43**]{}, 4418 (2010).

J. S. K[ł]{}os and J.-U. Sommer, J. Chem. Phys. [**134**]{}, 204 (2011).

$r_\protect \mathrm {d}$ in our previous works is defined as the location of the major peak of the sulfate density distribution. .

R. H. Stokes and R. A. Robinson, J. Am. Chem. Soc. [**70**]{}, 1870 (1948).

Y. Marcus and G. Hefter, Chem. Rev. [**106**]{}, 4585 (2006).

Y. Yao, M. L. Berkowitz, and Y. Kanai, J. Chem. Phys. [**143**]{}, 241101 (2015).

M. Pavlov, P. E. M. Siegbahn, and M. Sandstr[ö]{}m, J. Phys. Chem. A [**102**]{}, 219 (1998).

M. Kohagen, P. E. Mason, and P. Jungwirth, J. Phys. Chem. B [**118**]{}, 7902 (2014).

C. Yigit, N. Welsch, M. Ballauff, and J. Dzubiella, Langmuir [**28**]{}, 14373 (2012).

F. T. Wall and J. Berkowitz, J. Chem. Phys. [**26**]{}, 114 (1957).

D. McQuarrie, , University Science Books, 2000.

J. Kierzenka and L. F. Shampine, ACM Transactions on Mathematical Software [**27**]{}, 299 (2001).

S. Yu et al., Soft Matter [**11**]{}, 4630 (2015).

P. W. Atkins and J. De Paula, , W. H. Freeman and Company, 2010.

M. Gilson and H.-X. Zhou, Annu. Rev. Biophys. [**36**]{}, 21 (2007).

I. J. General, J. Chem. Theory Comput. [**6**]{}, 2520 (2010).

D. Harries, R. Podgornik, V. A. Parsegian, E. Mar-Or, and D. Andelman, J. Chem. Phys. [**124**]{}, 224702 (2006).

D. Boris and M. Rubinstein, Macromolecules [**29**]{}, 7251 (1996).

I. Berndt, J. S. Pedersen, and W. Richtering, Angew. Chemie Int. Ed. [**45**]{}, 1737 (2006).

R. J. Hunter, , Oxford university press, 2001.

[^1]: ${r_\mathrm{d}}$ in our previous works is defined as the location of the major peak of the sulfate density distribution. [@xu2017charged; @nikam2018charge]
