Electrostatically assisted macroion association

A model system of highly asymmetric polyelectrolyte with directional short-range attractive interactions was studied by canonical Monte Carlo computer simulations. Comparison of MC data with previously published theoretical results shows good agreement. For moderate values of binding energies, which matches those of molecular docking, a dynamic equilibrium between free and dimerized macroions is observed. Fraction of dimerized macroions depends on macroion concentration, binding energy magnitude, and on the valency of small counterions. Divalent counterions induce an effective attraction between macroions and enhance dimerization. This effect is most notable at low to moderate macroion concentrations.


Introduction
Charged colloidal systems are of great scientific and technological importance, as they include systems such as micellar solutions, metal-oxide nanoparticle suspensions, and protein solutions [1]. These systems have long been subject to extensive experimental and theoretical research. One of the most widely used theoretical approaches used to study the properties of charged colloidal systems is the classical Derjaguin-Landau-Verwey-Overbeek (DLVO) theory [2,3]. In contrast to this theory, which predicts purely repulsive forces between similarly-charged colloidal particles and was used to explain stability of colloidal systems, several studies presented evidence of effective attractive forces between macroions in these systems. Khan et. al. [4] reported an unusual phase behaviour of surfactant systems with divalent counterions. Kjellander and Marčelja observed an attractive force between two mica surfaces in presence of divalent ions [5]. Ise observed void formation in salt-free latex suspensions [6]. In addition to experimental work, computer simulations were also extensively used to address these observations [7,8]. Later, it has been confirmed that asymmetric polyelectrolytes undergo phase separation [9,10].
Associating systems are particularly important, as they mimic protein association [11][12][13][14], enzymeligand binding, and molecular docking. Many proteins are prone to dimerization or oligomerization [15,16]; this process often increases protein stability or grants them additional functionality. In these systems, directional short-ranged interactions are at work along with other intermolecular forces. Structural and thermodynamic properties of liquids are successfully described with theories based on the Ornstein-Zernike integral equation. Associative hypernetted chain approach, developed by Kalyuzhnyi and coworkers [17,18], was further extended with Wertheim's [19,20] two-density theory for associating liquids to yield a powerful tool for studies of dimerizing charged colloidal systems [21].
Theoretical description of associating fluids has seen additional progress over the last decade in a form of the resummed thermodynamic perturbation theory for associating fluids with multiple bonding sites [22]. Patchy colloids can also be treated as associating systems and have been studied with integral equation theories and computer simulations [23,24].
The aim of this work is to estimate a contribution of long-range electrostatic interactions to macroion association and to provide a comparison with already published theoretical data [21]. J. Reščič

Model and method
The salt-free polyelectrolyte solution model consists of spherical macroions with radius M = 15 Å and spherical counterions with radius C = 2 Å. Macroion charge M is equal to −15 0 while the counterion charge C is 1 0 or 2 0 .
Total potential energy of the model system is assumed to be pairwise additive and consists of three contributions: (i) hard-sphere repulsion, (ii) electrostatic interactions, and (iii) associative interactions.
with HS = < HS ( ) and and are radii of particles and , respectively. The is the corresponding centre-to-centre separation.
Electrostatic term has the following form: where is the valency of particle , 0 is the dielectric permittivity of vacuum and is the relative permittivity of water which is equal to 78.4 at a room temperature.
Association is possible only between macroions, as they possess one associative site per molecule. The origin of an associative site is located on the surface of a macroion and the site is spherical in shape. The associative part of the intermolecular potential is given by: Δ is the square well width and is chosen to be equal to 0.264 M , which allows only for formation of dimers. Since one of the goals of this research is to provide a comparison with the existing theoretical study, the model parameters are the same as in the work of Kalyuzhnyi and Vlachy [21]. The model is schematically represented in figure 1. Linse thoroughly investigated salt-free colloidal suspensions and proposed a set of three parameters which describe these systems [8]. These are macroion-counterion charge ratio = | M / C |, macroion volume fraction Φ M = 4π 3 M /3, and electrostatic coupling parameter Γ = 2 C / M , with = 2 0 /(4π 0 ) being the Bjerrum length and = 1/( B ). Large values of Γ indicate phase transition of a colloidal system.

Monte Carlo simulations
The model system was used with the canonical ensemble Monte Carlo computer simulation following the standard Metropolis sampling scheme. A simulation box was a cube with dimensions which correspond to a given macroion concentration. The number of macroions M was 64 while the number of counterions was C = M × | M / C | to fulfill electroneutrality. Macroions were initially placed on a 3D grid while counterions were inserted randomly. For a production run, at least 1 million moves per particle were performed. Electrostatic interactions were calculated using Ewald summation. Along with the excess internal energy, radial distribution functions, and cluster size distributions were calculated. Moreover, the evidence of bonds (which macroions are bonded) was maintained and updated at every macroion move. A fraction of bonded macroions was calculated from the collected bond histogram. All simulations were performed at a temperature of 298 K. Several square well depth values were studied: assoc = 5, 8, 10, 20, 25, 33, and 43. Each assoc value can represent a different enzyme since each one has a unique binding energy for a matching ligand. Macroion concentration varied from 0.004 M to 0.05 M. Highly charged macroions are tightly surrounded by counterions, resulting in a very high 33502-2 rejection rate of macroion single-particle trial moves. To overcome this issue, cluster moves were applied. Typically, displacement parameters are adjusted to the values which yields 50% acceptance rate. This is not always the best strategy to achieve a fast convergence and some techniques have been proposed for choosing optimal displacement parameters [25]. For highly charged systems, a total r.m.s. displacement for each particle type is monitored and should be of the size of a simulation box or larger. Associated particles should be allowed to explore various bonded configurations [13]. All these approaches are included in the MOLSIM simulation package [26] which was used to perform simulations.

Results and discussion
The same model was treated via the AHNC/Wertheim formalism [21]. At first, several simulations have been carried out to compare theoretical results with MC data, which is displayed in figure 2. For a system without associative potential the agreement is very good. It must be noted that the values of associative parameter assoc of 20 B or higher result in almost unbreakable bond -once a dimer is formed, both macroions forming it stay bonded. Energies are therefore simply shifted to more negative values by the same amount for all concentrations studied at each value of the assoc . Concentrations lower than 0.004 M were not studied due to a slow convergence.
Realistic enzyme/ligand binding energies are smaller in magnitude and are in a range of 10 B . Values of assoc = 5, 8, and 10 were used throughout the rest of the study. Excess internal energies are presented in figure 3 and cover polyelectrolyte concentrations in the range from 0.004 M to 0.05 M. Weaker association potential allows for dynamic creation, destruction, and rearrangement of bonds between macroions. Excess internal energies are at lower polyelectrolyte concentrations similar to each other. At larger concentrations, more macroions are bonded which is reflected in lower energies than a system without associative part of the intermolecular potential. Panel a) of figure 3 applies to monovalent counterions. Divalent counterions induce a weak effective attraction between macroions, which results in a lower excess internal energy of the system. This is presented in the panel b) of figure 3.
Fraction of bonded macroions was calculated in all cases. Figure 4 presents results for three squarewell depth values. Larger assoc yields a larger fraction of dimers for all concentrations studied. At  energy which is evident for the smallest assoc = 5. Increasing the counterion valency from 1 to 2 has a profound effect. It is known that increasing electrostatic coupling destabilizes a colloidal suspension.
The Γ parameter, which is a measure for the electrostatic coupling, is for the model system with monovalent counterions equal to 0.476 and rises fourfold to 1.9 when divalent ions are present. Weak to moderate effective attractive forces start to appear among macroions. This in turn increases the probability  of overlapping the associative regions on macroions. A similar model system without associative potential undergoes the gas-liquid phase separation and has a critical point at Γ value of 2.6 [9]. Presence of divalent counterions leads to a larger fraction of bonded macroions regardless of the assoc value. Both electrostatics and association lower the energy of the system, which is displayed in figure 3. This effect is even more pronounced at low to moderate macroion concentration, where the increase in the fraction of dimerized macroions is nearly 10-fold. Each panel contains information about dimerization dynamics for two values of assoc and for both monovalent and divalent counterions. In all cases, systems with divalent counterions contain a larger fraction of dimerized macroions, while this fraction is the smallest in systems with monovalent counterions and smaller assoc . Larger value of assoc increases longevity of bonds. Systems with divalent counterions also reach the equilibrium fraction of bonded macroions faster than their monovalent counterparts.
To better illustrate the effect of increased electrostatic coupling on macroion association when switching from monovalent to divalent counterions, macroion-macroion radial distributions ( ) are plotted in figure 6. The upper panel presents ( ) in systems with monovalent counterions and different assoc . Only at the strongest association energy of 10 B , the ( ) differs qualitatively from the rest, as a small peak appears at the macroion-macroion distance of 36 Å.
Association becomes much more pronounced in systems with divalent counterions. Slight effective macroion-macroion attraction is deduced from the ( ) for a system without the associative potential as the broad peak around 56 Å disappears and a small peak appears at 38 Å. Electrostatic forces bring macroions close enough so that short range associative potential can take dimerization over. Smaller average macroion-macroion distance increases the probability of dimerization. The peak at 36 Å in ( ) is much larger and is also visible at smaller values of assoc .

Conclusions
A salt-free colloidal solution model was studied using the canonical Monte Carlo simulations. Each highly charged macroion possesses one site with short-range associative potential. Such a model is suitable to study the protein dimerization in aqueous solutions with focus on the effect of electrostatic forces on macroion-macroion dimerization. First, we compared the excess internal energies obtained from simulations with theoretical results available in the literature, finding a good agreement between the two approaches. Structures of the solutions in a form of radial distribution functions were also very similar. Since divalent counterions induce a weak effective attraction among macroions, it was of particular interest to evaluate the effect of long-range electrostatic forces on macroion dimerization. Much larger fractions of dimerized macroions were found even at low macroion concentration in presence of divalent counterions, showing the effect of increased electrostatic coupling.

Acknowledgement
Author acknowledges Slovenian Research Agency for financial support through the research program P1-0201.