Branched polycations such as poly-(amidoamine) (PAMAM) dendrimers have been extensively investigated for applications including drug delivery and gene transfection, where phosphate-buffered solutions are commonly employed to maintain pH. Experimental and simulation studies have reported phosphate accumulation around PAMAM and attributed it, at least in part, to phosphate-specific interactions. Here, we investigate whether such specific interactions are necessary to explain phosphate association. Combining constant-pH molecular simulations, potentiometric titrations, and Site Binding theory, we study the coupled ionization of G2 PAMAM dendrimers and phosphate ions in aqueous solution. By employing a coarse-grained model that intentionally excludes phosphate-specific short-range interactions, we isolate the effects of electrostatics and charge regulation. The simulations quantitatively reproduce the ionization behavior of PAMAM and phosphate and reveal a strongly asymmetric charge-regulation response: phosphate has only a minor effect on PAMAM ionization, whereas PAMAM significantly enhances phosphate ionization in its vicinity. This coupling drives phosphate accumulation around PAMAM, with maximum adsorption at circumneutral pH where both species are highly charged. At phosphate concentrations comparable to those used in phosphate-buffered saline, the adsorbed phosphates reduce the effective charge of the PAMAM-phosphate complex. The qualitative agreement with previous experimental observations demonstrates that phosphate accumulation does not necessarily require phosphate-specific interactions but can emerge as a generic consequence of electrostatic attraction and coupled ionization equilibria. More broadly, the results highlight that multivalent buffer ions can act as active regulators of macromolecular electrostatics rather than merely as pH-control agents or passive screening species.

Heparin-mimicking polyelectrolytes act as anticoagulants in vivo due to their highly negative charge and high charge density. Unlike heparin, these polymers can be structurally tailored to tune chain length, sulfation degree and binding affinity, thereby modulating their therapeutic action. Here, we present two such heparin-mimicking polyelectrolytes, derived from natural rubber as an inexpensive and readily available raw material. These polymers, sodium poly((sulfamate-carboxylate) isoprene) (PISC) and poly((amino-carboxylate) isoprene (PACIS), contain sulfonate groups, as well as weakly acidic and basic groups on adjacent carbons, which enable us to tune their net charge in the biological pH range, while maintaining solubility and anticoagulation activity. Moreover, their activity may be further optimized by adjusting chain length. To test this tunability, we prepared a set of PISC- and PACIS-based copolymers varying in chain length and dispersity by post-polymerization modification of cis-1,4-polyisoprene. We characterized their structure and functionalization degree using various methods, further analyzing with potentiometric titration and computer simulations how their net charge varies with pH and affects aggregation. Despite aggregation, all PISC- and PACIS-based copolymers displayed heparin-mimicking activity, measured as activated partial thromboplastin time, while remaining biocompatible. In this context, the longest-chain copolymers performed best, extending aPTT by more than threefold in comparison with the control at 3.8 mu g/mL, while maintaining biocompatibility up to 100 mu g/mL across various cell lines. These findings demonstrate that natural rubber-derived polyelectrolytes provide a promising alternative to heparin for antithrombotic therapy.

The charge of peptides and proteins, determined by the effective pKA of their residues, is one of the key factors determining their function. Molecular interactions modulate the effective pKAs, shifting them relative to the bare pKA of isolated functional groups. Computational costs of predicting the effective pKAs range from CPU days for all-atom simulations to fractions of a second for the most simplified approaches. Herein, we examined to what extent the higher computational cost is balanced by higher accuracy. As a model system, we used short peptide sequences systematically combining acidic and basic residues, previously studied experimentally and by all-atom simulations. We supplemented the existing results with our own all-atom simulations, a coarse-grained approach, and two fast approximate approaches: FPTS and pepKalc. Our hypothesis that the predictions from all-atom approaches are more accurate than the approximate ones was confirmed only partly. The pKA shifts from most approaches agreed qualitatively but differed quantitatively. Two out of the four all-atom approaches matched the experiments significantly better than the others, while the remaining two were comparable or worse than the approximate ones. Our key finding is that all-atom simulations can be more accurate; however, their accuracy depends not only on the chosen force field, but also on other technical details that are often overlooked. Furthermore, comparing pKA shifts predicted by various simulations requires accounting for different pKA inputs. We show how to shift them to a common reference, enabling for the first time a rigorous and unified comparison across diverse pKA prediction methodologies.

Abstract Polypeptides with ionizable side chains are weak polyelectrolytes. When multiple ionizable groups are present along the polymer chain, their acid–base equilibria couple through local electrostatic fields. This coupling shifts the effective pKa of each group relative to its intrinsic value. As a result, the ionization behavior systematically deviates from the ideal Henderson–Hasselbalch prediction. Coarse-grained (CG) bead-spring models are widely used to study this pH-dependent ionization of polypeptides, but the accuracy of their predictions depends on the structural details present in the model. To investigate this, we compare two CG models: a one-bead and a two-bead representation against potentiometric titration and steady-state fluorescence spectroscopy data for two random copolypeptides: a polyacid poly(Glu,Tyr) and a polyampholyte poly(Lys,Tyr). Both models qualitatively reproduce the experimentally observed deviations from ideal ionization behavior. However, the two-bead model, which explicitly separates backbone and side-chain positions, gives better quantitative agreement in most cases, while the one-bead model consistently overestimates the effect of charge regulation. Some features of the experimental curves were not reproduced by our model, probably because it neglected nonelectrostatic interactions, which support α-helix formation at low degrees of ionization. These findings highlight the importance of geometric detail in CG models. The use of two independent experimental techniques provides a strong benchmark for the validation of the simulation model.

We explore the possible reversible formation of hydrogels through electrostatic interactions between four-arm star-shaped block copolymers, consisting of a polyethylene glycol (PEG) inner block and either an anionic polystyrene sulfonate [PEG(27)-b-PSS108](4) or a zwitterionic polybetaine [PEG(27)-b-PCBMAAm(110)](4) as outer block. The combination of both can induce attractive or repulsive electrostatic interactions depending on the solution pH value and ionic strength. The polymers were synthesized using controlled atom transfer radical polymerization (ATRP). The charge of [PEG(27)-b-PCBMAAm(110)](4) was further investigated by potentiometric titration and zeta potential measurements. Using oscillatory shear rheology, we demonstrated the required conditions for hydrogel formation. Stable hydrogel formation is observed within a wide pH range (6.8 -9.5), corresponding to the protonation states of the carboxylic acid groups that facilitate electrostatic interactions. We also showed how the hydrogel stability is influenced by parameters like block copolymer concentration and ionic strength. Coarse-grained simulations provided molecular-scale insights, revealing charge regulation effects and the energetic favorability of electrostatic complexation up to high pH values. Overall, our results demonstrated the key design principles, as the polyelectrolyte length, ionic strength, and charge regulation effects, for the formation of partially reversible hydrogels, triggered by changes in the solution pH. Furthermore, we showed that understanding the desired conditions for hydrogel formation requires a combination of experimental characterization with modeling approaches.

Poly(amidoamine) (PAMAM) dendrimers are promising candidates for nucleic acid delivery; however, biocompatibility and transfection efficiency remain a challenge. Here, we investigated how the composition of short peptide tails conjugated to generation 2 PAMAM (G2) dendrimers influence DNA association and condensation across a range of pH values. Using a combination of potentiometric titrations, DNA precipitation assays, and coarse-grained molecular simulations with charge regulation, we show that the ionization of G2 dendrimers is strongly affected by both pH and proximity to DNA. Although charge regulation enhances dendrimer protonation and strengthens DNA association at low pH, DNA condensation by unmodified G2 remains largely insensitive to pH within the studied range. In contrast, conjugation of a single peptide tail introduces a pronounced pH dependence to DNA condensation. Histidine-containing conjugates exhibit the strongest response, with condensation efficiency decreasing markedly as the pH increases. Simulations reveal that the interaction strength between conjugates and DNA depends on both peptide composition and pH and that histidine-containing peptide tails become nearly neutral at physiological pH, contributing little to DNA binding. While single-conjugate simulations explain the trends in DNA association, they do not fully account for the observed condensation behavior, highlighting the importance of collective effects involving multiple conjugates. Overall, peptide conjugation transforms G2 PAMAM dendrimers from relatively pH-insensitive DNA condensing agents into pH-responsive DNA-binding systems. These findings provide molecular-level insight into the interplay between charge regulation, peptide composition, and DNA condensation.

The constant-pH Monte Carlo method is a popular algorithm to study acid-base equilibria in coarse-grained simulations of charge regulating soft matter systems including weak polyelectrolytes and proteins. However, the method suffers from systematic errors in simulations with explicit ions, which lead to a symmetry-breaking between chemically equivalent implementations of the acid-base equilibrium. Here, we show that this artifact of the algorithm can be corrected a-posteriori by simply shifting the pH-scale. We present two analytical methods as well as a numerical method using Widom insertion to obtain the correction. By numerically investigating various sample systems, we assess the range of validity of the analytical approaches and show that the Widom approach always leads to consistent results, even when the analytical approaches fail. Overall, we provide practical guidelines on how to use constant-pH simulations to avoid systematic errors, including cases where special care is required, such as polyampholytes and proteins.

When using dialysis ultra- or diafiltration to purify protein solutions, a dialysis buffer in the permeate is employed to set the pH in the protein solution. Failure to achieve the target pH may cause undesired precipitation of the valuable product. However, the pH in the permeate differs from that in the retentate, which contains the proteins. Experimental optimization of the process conditions is time-consuming and expensive, while accurate theoretical predictions still pose a major challenge. Current models of dialysis account for the Donnan equilibrium, acid-base properties, and ion-protein interactions, but they neglect the patchy distribution of ionizable groups on the proteins and its impact on the solution properties. Here, we present a simple computational model of a colloidal particle with weakly acidic sites on the surface, organized in patches. This minimalistic model allows systematic variation of the relevant parameters, while simultaneously demonstrating the essential physics governing the acid-base equilibria in protein solutions. Using molecular simulations in the Grand-Reaction ensemble, we demonstrate that interactions between ionizable sites significantly affect the nanoparticle charge and thereby contribute to pH difference between the permeate and retentate. We show that the significance of this contribution increases if the ionizable sites are located on a smaller patch. Protein solutions are governed by the same physics as our simple model. In this context, our results show that models which aim to quantitatively predict the pH in protein solutions during dialysis need to account for the patchy distribution of ionizable sites on the protein surface.

Uptake of proteins and ampholytic solutes into polyelectrolyte brushes underlies some biological processes and also applications in sensing or biomedicine. Especially uptake on the "wrong" side of the isoelectric point (pI) remains puzzling, with charge regulation and solute patchiness proposed as possible mechanisms. Using a hierarchy of approximations, coarse-grained molecular simulations, self-consistent mean-field, and a simple phenomenological model, we investigated the uptake of model ampholytic solutes into polyanionic brushes across varying pH, salt concentrations, pKa values, and peptide sequences. In a narrow pH range on the wrong side of pI, charge regulation enables uptake of the ampholytes by inducing charge inversion so that they become positively charged in the brush despite being negatively charged in the bulk. This charge inversion can be calculated from the pH difference between the brush and the bulk, which is related to the Donnan potential. It is strongest for ampholytes with small differences between acidic and basic pKa values and decreases with increasing salt. Our phenomenological model reproduces the universal effect of charge regulation promoting ampholyte uptake into brushes but fails to be quantitative. The mean field model is close to explicit simulations for alternating sequences, but fails to describe the effect of charge patchiness, which is only captured by explicit simulations. Thus, our phenomenological framework offers a practical rule of thumb for estimating uptake from experimentally accessible parameters without sophisticated calculations. Deviations from this rule of thumb for complex ampholytes, such as proteins or peptides with patterned charge sequences, are captured only by explicit simulations.

Particle-based coarse-grained simulations bridge the scales between models with atomistic details and a continuum description. They are typically applied to systems ranging from the nanometer scale to the micron scale. Examples for these systems are polymers, colloids or biological cells. In this chapter we will first provide an overview of simulation techniques suitable for coarse-grained models, in particular molecular dynamics and Monte Carlo schemes. Afterwards the simulation package ESPResSo is shortly described, together with simulation models from different fields that illustrate its new capabilities. In particular, we highlight simulations of weak polyelectrolytes, hydrogels, soft magnetic materials and chemical reactions subject to a flow field.

We present the Python-based Molecule Builder for ESPResSo (pyMBE), an open source software application to design custom coarse-grained (CG) models, as well as pre-defined models of polyelectrolytes, peptides, and globular proteins in the Extensible Simulation Package for Research on Soft Matter (ESPResSo). The Python interface of ESPResSo offers a flexible framework, capable of building custom CG models from scratch. As a downside, building CG models from scratch is prone to mistakes, especially for newcomers in the field of CG modeling, or for molecules with complex architectures. The pyMBE module builds CG models in ESPResSo using a hierarchical bottom-up approach, providing a robust tool to automate the setup of CG models and helping new users prevent common mistakes. ESPResSo features the constant pH (cpH) and grand-reaction (G-RxMC) methods, which have been designed to study chemical reaction equilibria in macromolecular systems with many reactive species. However, setting up these methods for systems, which contain several types of reactive groups, is an error-prone task, especially for beginners. The pyMBE module enables the automatic setup of cpH and G-RxMC simulations in ESPResSo, lowering the barrier for newcomers and opening the door to investigate complex systems not studied with these methods yet. To demonstrate some of the applications of pyMBE, we showcase several case studies where we successfully reproduce previously published simulations of charge-regulating peptides and globular proteins in bulk solution and weak polyelectrolytes in dialysis. The pyMBE module is publicly available as a GitHub repository (https://github.com/pyMBE-dev/pyMBE), which includes its source code and various sample and test scripts, including the ones that we used to generate the data presented in this article.

Electrostatic interactions between charged macromolecules are ubiquitous in bio- logical systems and they are important also in materials design. Attraction between oppositely charged molecules is often interpreted as if the molecules had a fixed charge, which is not affected by their interaction. Less commonly, charge regulation is invoked to interpret such interactions, i.e., a change of the charge state in response to a change of the local environment. Although some theoretical and simulation studies suggest that charge regulation plays an important role in intermolecular interations, experi- mental evidence supporting such view is very scarce. In the current study, we used a model system, composed of a long polyanion interacting with cationic oligolysines, containing up to 8 lysine residues. We showed using both simulations and experiments that while these lysines are only weakly charged in the absence of the polyanion, they charge up and condense on the polycations if the pH is close to the pKa of the lysine side chains. We show that the lysines coexist in two distinct populations within the same solution: 1. practically non-ionized and free in solution; 2. highly ionized and condensed on the polyanion. Using this model system, we demonstrate under what conditions charge regulation plays a significant role in the interactions of oppositely charged macromolecules and generalize our findings beyond the specific system used here.

Mixing of oppositely charged polyelectrolytes can result in phase separation into a polymer-poor supernatant and a polymer-rich polyelectrolyte complex (PEC). We present a new coarse-grained model for the Grand-reaction method that enables us to determine the composition of the coexisting phases in a broad range of pH and salt concentrations. We validate the model by comparing it to recent simulations and experimental studies, as well as our own experiments on poly(acrylic acid)/poly(allylamine hydrochloride) complexes. The simulations using our model predict that monovalent ions partition approximately equally between both phases, whereas divalent ones accumulate in the PEC phase. On a semiquantitative level, these results agree with our own experiments, as well as with other experiments and simulations in the literature. In the sequel, we use the model to study the partitioning of a weak diprotic acid at various pH values of the supernatant. Our results show that the ionization of the acid is enhanced in the PEC phase, resulting in its preferential accumulation in this phase, which monotonically increases with the pH. Currently, this effect is still waiting to be confirmed experimentally. We explore how the model parameters (particle size, charge density, permittivity, and solvent quality) affect the measured partition coefficients, showing that fine-tuning of these parameters can make the agreement with the experiments almost quantitative. Nevertheless, our results show that charge regulation in multivalent solutes can potentially be exploited in engineering the partitioning of charged molecules in PEC-based systems at various pH values.

The constant-pH ensemble method is a popular algorithm to study acid-base equilibria in charge regulating soft matter systems including weak polyelectrolytes and proteins. However, the method suffers from systematic errors in simulations with explicit ions, which lead to a symmetry-breaking between chemically equivalent implementations of the acid-base equilibrium. Here, we show that this artifact of the algorithm can be corrected a-posteriori by simply shifting the pH-scale. We present two analytical methods as well as a numerical method using Widom insertion to obtain the correction. By numerically investigating various sample systems, we assess the range of validity of the analytical approaches and show that the Widom approach always leads to consistent results, even when the analytical approaches fail. Overall, we provide practical guidelines on how to use constant-pH simulations to avoid systematic errors, including cases where special care is required, such as polyampholytes and proteins.

We present the explicit bonding Reaction ensemble Monte Carlo (eb-RxMC) method, designed to sample reversible bonding reactions in macromolecular systems in thermodynamic equilibrium. Our eb-RxMC method is based on the reaction ensemble method; however, its implementation differs from the latter by the representation of the reaction. In the eb-RxMC implementation, we are adding or deleting bonds between existing particles, instead of inserting or deleting particles with different chemical identities. This new implementation makes the eb-RxMC method suitable for simulating the formation of reversible linkages between macromolecules, which would not be feasible with the original implementation. To enable coupling of our eb-RxMC algorithm with molecular dynamics algorithm for the sampling of the configuration space, we biased the sampling of reactions only within a certain inclusion radius. We validated our algorithm using a set of ideally behaving systems undergoing dimerization and polycondensation reactions, for which analytical results are available. For dimerization reactions with various equilibrium constants and initial compositions, the degree of conversion measured in our simulations perfectly matched the reference values given by the analytical equations. We also showed that this agreement is not affected by the arbitrary choice of the inclusion radius or the stiffness of the harmonic bond potential. Next, we showed that our simulations can correctly match the analytical results for the distribution of the degree of polymerization and end-to-end distance of ideal chains in polycondensation reactions. Altogether, we demonstrated that our eb-RxMC simulations correctly sample both reaction and configuration spaces of these reference systems, opening the door to future simulations of more complex interacting macromolecular systems.

Levin and Bakhshandeh suggested in their comment that (1), we stated in our recent review that pH-pK(A) is a universal parameter for titrating systems, that (2), we omitted to mention in our review the broken symmetry of the constant pH algorithm, and that (3), a constant pH simulation must include a grand-canonical exchange of ions with the reservoir. As a reply to (1), we point out that Levin and Bakhshandeh misquoted and hence invalidated our original statement. We therefore explain in detail under which circumstances pH-pK(A) can be a universal parameter, and also demonstrate why their numerical example is not in contradiction to our statement. Moreover, the fact that pH-pK(A) is not a universal parameter for titrating systems is well known in the pertinent literature. Regarding (2), we admit that the symmetry-breaking feature of the constant pH algorithm has escaped our attention at the time of writing the review. We added some clarifying remarks to this behavior. Concerning (3), we point out that the grand-canonical coupling and the resultant Donnan potential are not features of single-phase systems, but are essential for two-phase systems, as was shown in a recent paper by some of us, see J. Landsgesell et al., Macromolecules, 2020, 53, 3007-3020.

Recent experiments on weak polyelectrolyte brushes found marked shifts in the effective pK_{a} that are linear in the logarithm of the salt concentration. Comparing explicit-particle simulations with mean-field calculations we show that for high grafting densities the salt concentration effect can be explained using the ideal Donnan theory, but for low grafting densities the full shift is due to a combination of the Donnan effect and the polyelectrolyte effect. The latter originates from electrostatic correlations that are neglected in the Donnan picture and that are only approximately included in the mean-field theory. Moreover, we demonstrate that the magnitude of the polyelectrolyte effect is almost invariant with respect to salt concentration but depends on the grafting density of the brush. This invariance is due to a complex cancellation of multiple effects. Based on our results, we show how the experimentally determined pK_{a} shifts may be used to infer the grafting density of brushes, a parameter that is difficult to measure directly.

To study macroscopic systems with coarse grained simulations one typically simulates a micro- scopic part of this macroscopic system. By reducing the size of the simulated system one introduces finite size effects. In this work we study the finite-size effects in the reaction ensemble, which is used to simulate reactive system. We calculate the finite-size effects in a non-interacting systems by explicitly calculating the partition function. This approach provides high precision data at low computational costs. For a grand canonical insertion/deletion of a pair of particles our results reproduces previously published results, validating our approach. Further, we show that a sim- ple isomerization reaction is not affected by finite size effects. For a decomposition reaction we show that previous estimates were overestimating the finite-size effects, and one can simulate much smaller systems while avoiding the finite-size effects. For previously studied acid-base equilibria the finite-size effects are only relevant at extreme conditions. The tool we provide allows to a priori estimate the finite-size effects and find the limits of the applicability of the reaction ensemble.

We synthesized three different polyzwitterions-poly(N,N-diallyl glutamate) (PDAGA), poly(dehydroalanine) (PDha), and poly(2-(imidazol-1-yl)acrylic acid) (PImAA)-and investigated how their ionization states respond to changes in solution pH. We used molecular simulations to determine how the net charge per monomer and the ionization states of individual acidic and basic groups differ from the ideal (Henderson-Hasselbalch) behavior. To complement the theoretical predictions, we performed potentiometric titrations and zeta-potential measurements of all studied polyzwitterions. By comparing these experiments with theoretical predictions, we could show that molecular simulations can predict and explain the origin of the differences between the effective and bare pK(a) values of individual titratable groups. Furthermore, we have shown that it is not possible to obtain these effective pK(a) values directly from the equivalence point recognition criterion (ERC), commonly used in potentiometric titrations. However, the effective pK(a) values can be reliably obtained by calculating the net charge per monomer from the potentiometric titration curves and validating these results against theoretical predictions. The approach we propose works reliably for polyzwitterions in which the ionization response is dominated by electrostatic interactions, such as PDAGA or PDha; however, it fails if other specific interactions contribute significantly, such as in the case of PImAA.

We developed a new method for coarse-grained simulations of acid-base equilibria in a system coupled to a reservoir at a given pH and concentration of added salt, that we term the Grand-reaction method. More generally, it can be used for simulations of any reactive system coupled to a reservoir of a known composition. Conceptually, it can be regarded as an extension of the reaction ensemble, combining explicit simulations of reactions within the system and Grand-canonical exchange of particles with the reservoir. To demonstrate its strength, we applied our method to a solution of weak polyelectrolytes in equilibrium with a reservoir. Our results show that the ionization and swelling of a weak polyelectrolyte are affected by the Donnan effect due to the partitioning of ions and by the polyelectrolyte effect due to electrostatic repulsion along the chain. Both effects lead to a similar shift in ionization and swelling as a function of pH, albeit for different physical reasons. By comparison with published results, we showed that neglecting one or the other effect may lead to erroneous predictions or misinterpretations of results. In contrast, the Grand-reaction method accounts for both effects on the results and allows us to quantify them. Finally, we outline possible extensions and generalizations of the method and provide a set of guidelines for its safe application by a broad community of users.