Effects of electrostatic interactions on ligand dissociation kinetics

We study unbinding of multivalent cationic ligands from oppositely charged polymeric binding sites sparsely grafted on a flat neutral substrate. Our molecular dynamics (MD) simulations are suggested by single-molecule studies of protein-DNA interactions. We consider univalent salt concentrations spanning roughly a thousandfold range, together with various concentrations of excess ligands in solution. To reveal the ionic effects on unbinding kinetics of spontaneous and facilitated dissociation mechanisms, we treat electrostatic interactions both at a Debye-H\"{u}ckel (DH, or `implicit' ions, i.e., use of an electrostatic potential with a prescribed decay length) level, as well as by the more precise approach of considering all ionic species explicitly in the simulations. We find that the DH approach systematically overestimates unbinding rates, relative to the calculations where all ion pairs are present explicitly in solution, although many aspects of the two types of calculation are qualitatively similar. For facilitated dissociation (FD, acceleration of unbinding by free ligands in solution) explicit ion simulations lead to unbinding at lower free ligand concentrations. Our simulations predict a variety of FD regimes as a function of free ligand and ion concentrations; a particularly interesting regime is at intermediate concentrations of ligands where non-electrostatic binding strength controls FD. We conclude that explicit-ion electrostatic modeling is an essential component to quantitatively tackle problems in molecular ligand dissociation, including nucleic-acid-binding proteins.


I. INTRODUCTION
Electrostatic and non-electrostatic interactions between molecular ligands and their binding sites control many important aspects in biomolecular machinery from gene regulation to molecular recognition. Nonelectrostatic contributions (e.g., van der Waals) are usually attributed to interactions between charge-neutral groups, whereas effects due to structural charges (e.g., phosphate groups on nucleic acids) and solvated ionic species are the subjects of electrostatics. Cumulatively, non-electrostatic and electrostatic interactions determine the lifetime (i.e., inverse of unbinding rate) of a ligand on its binding site.
Experimentally, probing the role of non-electrostatic interactions in the unbinding process is possible, for instance, by testing different ligand-receptor pairs by varying one of the partners. However, changing either the ligand or the host molecule inevitably changes the electrostatic interactions as well, since the biomolecules involved, e.g., proteins, are usually amphiphilic and complex structures [1,2]. Manipulating solution strength is an alternative way of inferring binding thermodynamics. This is mainly because salt ions in solution impose an electrostatic screening length scale (i.e., the Debye length), which can be used as a ruler to probe related kinetic rates (the Debye length defines the volume, in which the electrostatic energy of concentration fluctuations is balanced by the thermal energy). Indeed, many workers have used salt concentration as a tool to study the role of electrostatic and non-elecrostatic interactions in the dissociation kinetics of biological ligands [3][4][5][6][7][8][9]. Extensive studies have shown that the affinity of nucleic acid binding proteins [8][9][10][11][12] and oligocations [13][14][15][16] decreases with increasing univalent salt concentration.
Solvated univalent salt ions in solution can affect the unbinding kinetics of a ligand in various ways. First, increasing the salt concentration weakens Coulomb interactions between a pair of charges separated by a distance larger than the Debye length. Similarly, the screened Coulomb interactions between structural charges on the binding site and those on the ligand can promote dissociation by lowering free energy barriers of binding. Second, changing the salt concentration can alter the ionic distributions and correlations, particularly near the charged macromolecules [17,18]. Hence, salt-induced competition between various entropic and enthalpic components can lead to ligand dissociation. Nevertheless, separation of the salt-related contributions from other effects is at least approximately possible via extrapolation of dissoci-ation rates to high salt limits [3,8].
Besides univalent salt ions, an excess amount of free ligands in solution can also facilitate dissociation of a bound ligand by decreasing the lifetime of the ligand on the binding site [19][20][21][22][23][24][25][26][27]. The free ligands can be identical to the initially bound ligand [21,22,[27][28][29][30]. However, different ligand molecules [27,30,31] or even nucleic acid fragments [20] can lead to the facilitated dissociation (FD). According to the proposed molecular mechanism for FD, a free ligand binds to an already occupied binding site and destabilizes the complex by forming a ternary complex (i.e., two ligands on the same site) [19,23,27]. This destabilized complex promotes shorter binding lifetimes for the ligands, while imposing an upper limit on the unbinding rates.
A recent study on FD has shown that free ligands with concentration ranges on the order of few hundred nano mole can significantly weaken the salt dependence of unbinding rates [27]. The ligand concentrations, at which strong deviations from ligand-free assays were observed, are protein concentrations found in cells. The fact that excess ligand and univalent salt can cause analogous effects on the unbinding kinetics suggests that the explicit nature of the ions and structural charges of ligands should be considered carefully.
Screened electrostatic potentials of the Debye-Huckel (DH) type have a mean-field nature, and provide a way to account for the decay of electrostatic forces between two charges separated by a distance larger than the Debye length in salt solutions. Although these mean-field potentials satisfy the Poisson-Boltzmann equation at long ranges, clearly they cannot account for the entropy or positional correlations of the screening ions. In addition, micro-dissociation events, whose time scales are likely comparable to the relaxation time of the ionic correlations, can be smeared out completely in mean fieldlevel treatments. Nevertheless, DH and other mean-field methods are attractive in studying charged biological, as well as synthetic systems, due to their advantage in computational and mathematical treatment [32][33][34][35][36][37].
To clearly understand the unbinding kinetics, the role of electrostatic interactions and the possible interplay between the ionic species and unbound ligands should be investigated. Given the variety of methods by which electrostatics can be considered, in this work we aim to answer specific questions including whether the way that electrostatic interactions are treated (i.e., by either modeling all ionic species explicitly or considering their effect via a mean-field DH approach) has any effect on the kinetics of unbinding. Furthermore, we inquire whether the electrostatic treatment influences the FD process, in which both the non-electrostatic binding energy and ionic contributions are expected to play a major role. Most importantly, we can compare the results obtained from the different electrostatic treatments to identify various electrostatic contributions to the unbinding mechanism. We investigate these issues by using Molecular Dynamics (MD) simulations, in which electrostatic interactions are considered carefully to distinguish explicit and meanfield level contributions of charged species. The simulations also allow us to manipulate the non-electrostatic pairwise interactions directly without significantly altering electrostatic interactions. Hence, we can probe the contributions of non-electrostatic binding strengths and ionic effects on the unbinding kinetics.
In our study, we monitor the unbinding events of single cationic ligands from short polymeric binding sites with opposite charges (Fig. 1). Many binding sites sparsely placed on a charge-neutral surface allow us to obtain accurate statistics from uncorrelated binding events. Similar setups are commonly used to experimentally investigate unbinding rates for nucleic acid binding proteins [38] and are also used in single-molecule studies of protein-DNA interactions [22,27]. In the simulations, we investigate two orders of magnitude of univalent salt concentration by considering two approaches that can allow us to distinguish the effect of the ionic component in the kinetic events. In the first approach, all ions (i.e., salt ions and solvated counterions of the ligands and binding sites) are modeled explicitly, and a pairwise Coulomb potential is calculated between all pairs of ions. We refer to these simulations as "explicit" throughout this work. In the second approach, the ions are removed from the simulation boxes, but a DH screened electrostatic potential implicitly takes the effect of the ions into account. We refer to these simulations as "implicit".
Our MD simulations and data analysis reveal that the implicit DH simulations systematically overestimate the unbinding rates of the ligands, and that this discrepancy is more dramatic near physiological salt conditions (c s ≈ 100 mM). The deviation between explicit and implicit treatments of ionic effects persists in the simulations where we have a "concentration quench" (i.e., at the initial step of simulations, all ligands are bound to their binding sites, and there is no ligand in solution), and also in the cases, where the FD effect is triggered by the excess amount of free ligands in solution. The finding that rebinding rates are nearly independent of electrostatic treatment (particularly at physiological salt conditions) indicate that the difference is related to local interactions and is not a simple consequence of the coarse-graining of the ions. Furthermore, aside from investigation of electrostatics, our simulations underline various dissociation regimes, in which the FD process at intermediate ligand concentrations depends on the non-electrostatic binding affinity between the ligands and the binding site. This dependency disappears at high ligand concentrations.
The paper is organized as follows. First, we describe the general simulation methodology and the two methods employed in the simulations to calculate electrostatic interactions. Next, the results obtained for spontaneous dissociation simulations are discussed by considering a multi-step dissociation model. In the next section, we focus on FD. We discuss our results and possible indications for biological systems together with future prospects in the Conclusion section. where all ionic species are treated explicitly. The red beads along the grafted semi-flexible chains attract bound and free proteins via an attractive short-range potential in addition to the electrostatic interactions. b) Top view of the actual simulation box with 10 × 10 binding sites. Ions and ligands are not shown for clarity. c) A single binding site is zoomed in. d) A single binding site used in the implicit-ion simulations, where the electrostatic interactions are calculated via Eq. 5. All beads in the systems except those forming the surface carry respective charges (see the text for details).

A. Details of MD simulations
In the simulation model, at least n 0 = 10 × 10 binding sites are sparsely grafted on an inert surface with an inter-site distance of d = 24 σ , where σ is the size of a unit bead (Fig. 1). The aqueous medium is modeled implicitly as a continuum (see below). The binding sites and ligands are modeled as coarse-grained Kremer-Grest (KG) bead-spring chains [39,40]. Each binding site is a linear semi-flexible polymer chain composed of N = 12 identical beads. A linear chain of p = 4 beads is placed onto each binding site to model initially bound ligand molecules (green beads in Fig. 1). For the FD simulations, a prescribed concentration of initially free ligands, c p , is added at random positions in the simulation boxes in addition to the initially bound ligands.
The bonding between the adjacent beads of the chains is taken care of by a nonlinear potential with finite ex- where the bond stiffness is k = 30 ǫ/σ 2 , the distance between adjacent beads is r = |r|, and the maximum bond length is r 0 = 1.5 σ [39]. The interaction strength ǫ was measured in the units of thermal energy k B T , where k B is the Boltzmann constant, and T is the absolute temperature.
The steric interactions between all beads are modeled by a truncated and shifted Lennard-Jones (LJ) potential, also know as WCA, where r c is the cutoff distance. A cut-off distance of r c = 2 1/6 σ is used with a shift factor v s = 1/4 to obtain good solvent conditions for the interactions between all beads unless otherwise noted with an interaction strength of 1 ǫ. Four of the N = 12 beads (the red beads in Fig. 1) interact with a strength of 2 ǫ and with r c = 2.5 σ and v s = 0 unless noted otherwise to mimic the specific binding sites for the ligands (e.g., protein specific binding sites along nucleic acid chains). A harmonic bending potential is introduced for the grafted N = 12 chains to account for the semi-flexible nature of the binding site (e.g., DNA) in the form of where the potential strength is k Θ = 30 ǫ/rad 2 , Θ is the angle formed by three adjacent beads, and Θ 0 is the reference angle. Θ 0 = π for all grafted chains except the beads connecting the chains to the surface. The grafted chains are kept at a right angle with the surface by setting Θ 0 = π/2 and k Θ = 90 ǫ/rad 2 for the grafted beads. All MD simulations are run with LAMMPS MD package [41] at constant volume V and T = 1.2 ǫ/k B . During the relaxation of initial configurations, each system is simulated for 10 5 MD steps by keeping the bound ligands on their binding sites by temporarily replacing 1 ǫ in Eq. 2 with 10 ǫ. The data production runs are carried out until a detailed balance is reached for the ligand un/binding events, which is between 10 6 − 10 8 MD steps. The simulations are run with a time step of ∆t = 0.005 τ , where the unit time scale in the simulations is τ = σ m/ǫ with a monomeric LJ mass of m = 1. The temperature is kept constant by a Langevin thermostat with a thermostat coefficient γ = 0.5 τ −1 . The volume of the total simulation box is set to L x × L y × L z = 232 × 232 × 58 σ 3 . The vertical size of the boxes (i.e.,L z = 58 σ) is higher than the effective effective Gouy-Chapman length of the surface due to the grafted charged chains (i.e., λ GC ≈ d 2 /2πℓ B N ≈ 10 σ).

B. Calculation of electrostatic interactions in the simulations
Each bead of the grafted chains is assigned a unit negative charge, whereas each ligand bead bears a positive unit charge regardless of the method used for the calculation of the electrostatics.
In the explicit case, for each charged bead, one oppositely charged counterion bead is added in the simulation box at a random position to obtain a charge-neutral systems even in the absence of salt. Fixed numbers of positively and negatively charged beads are added in the simulation boxes to model salt concentrations in the range of c s ≈ 5 − 1000 mM. The sizes of both counterions and salt monomers are taken to be 0.4 σ, which is smaller than the beads forming the binding sites and ligands. The ionic species interact with each other and with the other components via a shifted 9-6 LJ potential instead of the 12-6 potential given in Eq. 2 to account for the effect of (softer) hydration layers. The electrostatic interactions between two charged beads are accounted for by a pairwise Coulomb potential In Eq. 4, the Bjerrum length ℓ B quantifies the distance at which the Coulomb energy between two beads of size 1 σ at contact is equivalent to 1 k B T . In aqueous medium, ℓ B ≈ 0.7 nm. In the simulations, the dielectric constant is adjusted to obtain ℓ B = 1 σ. This adjustment corresponds to the charge density of Γ OM = 1/σ = 1/ℓ B as imposed by the Manning condensation effect [42]. The electrostatic cut-off distance in Eq. 4 is r e = 8 − 12σ, above which longer-range electrostatic interactions are set to zero. Particle-Particle-Particle Mesh (PPPM) Ewald solver with an error tolerance 10 −3 was also tested but no significant deviations were observed for the given salt concentrations (see SI Figure S1). Note that the maximum cut-off distance used here is half the distance between two binding sites on the surface. In the implicit case, all the counterion and salt beads are removed from the simulation boxes, and the electrostatic interactions between the remaining charged species (i.e., ligands and binding-sites) are calculated via a screened potential where κ −1 is the Debye screening length, which is defined as κ −1 = k B T / √ 8πℓ B c s in a solution of univalent salt. The values of κ −1 are varied to obtain effective salt concentrations in accord with the explicit-ion simulations. Note that the lowest salt concentration (i.e., c s = 5 mM) considered in this work corresponds to a Debye screening length of κ −1 ≈ 4 nm in real units.

C. Extraction of rates
In the simulations, the number of ligands remaining bound on their binding sites, n(t), is monitored as a function of the simulation time, t. If any bound ligand diffuses out of a spherical volume with a radius R c ≈ 4σ, centered around the center of mass of the binding site (red beads in Fig. 1), the ligand is tagged as unbound. If the ligand returns to the binding site, it is not counted as bound. To determine the off-rates both in the absence and presence of unbound ligands, the survival fraction data is fit by a single exponential where k off is defined as the inverse of the lifetime of the bound ligands on their binding sites (i.e., k off ≡ 1/τ off ). Error bars are calculated by averaging the results of multiple runs. In addition, n 0 in Eq. 6 is released as a fit parameter to enrich error statistics. In the fitting procedures, a weight function inversely proportional to the square of the data point is used. Error bars are not shown if they are smaller than the size of the corresponding data point. VMD is used for the visualizations [43].

A. Spontaneous dissociation
1. Off-rates for explicit salt are significantly less than those for implicit salt We first focus on the case, in which initially all binding sites are occupied by single ligands, and there is no free ligand in solution. This setup corresponds to a concentration quench, from which an initially high concentration of bound ligands desorbs from the binding sites to an initially zero-concentration bulk state. Similar setups have commonly been used in experiments studying unbinding rates for protein-DNA interactions [22,27,38].
In Fig. 2a, we show typical time traces for the normalized fraction, n(t)/n 0 , of bound ligands that are bound at time t, as a function of simulation time for both explicit (triangles in Fig. 2a) and implicit (open circles) cases at a salt concentration of c s ≈ 100 mM. As time progresses, chemical equilibrium between bound and unbound states is established, and all of the initially bound ligands leave their binding sites. For both explicit and implicit cases, the relaxation of the concentration quenches can be described well by single-exponential functions (dashed curves in Fig. 2a). However, the off-rate, k off , (i.e., the inverse dwelling time of the ligand on its binding site) obtained from the exponential fits for the explicit case is one order of magnitude slower than that for the implicit case. This indeed shows that the ligands remain on their binding sites for less time in the implicit simulations. Once the time axes are rescaled by the corresponding values of the off-rates, the data sets for the implicit and explicit cases are nearly indistinguishable (see the inset of Fig. 2a). This indeed confirms that both processes can be described by single-exponential decays independent of how we treat electrostatics.
To compare dependencies of the off-rates on univalent salt concentration (or on the Debye length), we performed similar titration simulations by varying the range of excess salt concentration between c s ≈ 5 − 1000 mM. Consistently, the simulations, in which we consider electrostatic interactions implicitly (open circles Fig. 2b) systematically exhibit higher unbinding rates compared to the explicit case (closed circles Fig. 2b) for the entire range of the salt concentration. However, in both cases, the off-rates have similar qualitative dependencies on the excess salt; at low salt, the off-rate increases gradually up to c s ≈ 100 mM, above which a saturation regime appears. In these saturation regimes, the salt has little or no effect on the off-rates. It is noteworthy that beyond c s > 100 , the electrostatic screening length of the solution is less than the size of a unit bead (i.e., κ −1 < 1 σ).
The plateau values of k off for the two cases are significantly different as seen in Fig. 2b.
To make a more systematic comparison between the off-rates obtained from the implicit and explicit cases and to quantify the observed regimes, we fit the data in Fig. 2b to a simplified version of a theoretical model. This model has been previously suggested for the unbinding rates of transcription factors (i.e., nucleic-acid binding proteins) from single binding sites along double stranded (ds)DNA [23,27]. The model assumes that unbinding of a multimeric ligand occurs in a multi-step process. Each step can have its own salt dependence. At the first step, the ligand partially unbinds. At the later steps, the partially dissociated ligand desorbs gradually into bulk solution. For the sake of simplicity and to minimize the number of free fit parameters, here we assume a two-step fashion and that only the first step has a salt dependence. This assumption is consistent with a previously reported analyses of the experimental data [27]. Thus, the off-rate can be expressed as the sum of two reaction times as where α, z, k sat are free fit parameters to be determined for the implicit and explicit data sets separately. While the first term in parentheses in Eq. 7 has a power-law dependence on the excess salt concentration, c s , with an exponent z, the second term in parenthesis accounts for a salt-independent (saturation) process. Fitting our data to Eq. 7 gives z exp = 1.35 ± 0.02 for the explicit case, which is slightly lower than that for the implicit case, z imp = 1.46 ± 0.04. The finding that z imp > z exp suggests a convergence of the off-rates as c s → 0 (or as κ −1 → ∞). However, we should note that even at vanishing excess salt concentrations, the counterions of the ligands and binding-sites are present in the explicit simulations, unlike the implicit case. The exponent z is often considered to correspond to the number of univalent salt ions replacing the ligand upon dissociation [3,5,7]. In our case, the explicit ion simulations surprisingly lead to a similar exponent with the implicit cases, for which neither ions nor their excluded volumes are in present. This suggests that in our simulations the exponent z may be rather a result of mixed effect of ionic entropy and electrostatic screening of the Coulomb interactions between the binding sites and the ligands. We will revisit this topic in the Conclusion section.
The saturation rates k sat obtained by fitting the offrates to Eq. 7 (Fig. 2b) differ roughly factor of 30 (i.e., k imp sat = 0.0058 τ −1 vs. k exp sat = 0.0002 τ −1 ). At high salt concentrations, the electrostatic interactions between the charges significantly weaken, thus, the difference in the saturation rates can be attributed to the contributions related to the degrees of freedom of the ions (e.g., ionic correlations and translational entropy), which are present only in the explicit-ion case. Ion-ion correlations near the binding site are stronger [44] and thus they promote the bound ligand in favor of weaker correlations in bulk. In addition, since translational entropy of the ions decreases near the binding site, this entropy component also favors a bound ligand to maximize entropy [18].
We note that while we see saturation behaviors for the off-rates in both cases above c s > 100 mM in Fig. 2b, the plateau regime in the explicit simulations appears at a slightly smaller threshold salt concentration (i.e., c * s = (k sat /α) 1/z ); the onset of the saturated regime in the implicit case manifests itself at a roughly factor of ∼ 2 more saltier solution. This implies that the effects beyond the screening of electrostatic forces play a role in the salt-dependent behavior of unbinding in Fig. 2b.
To further demonstrate the difference between the implicit and explicit cases, in the inset of Fig. 2b, we show the ratios of the fitted functions of the off-rates. The inset demonstrates that the difference between the implicit and explicit cases is not simply a constant and has rather a salt dependent profile. At c s > ∼ 100 mM, the difference increases gradually and reaches a plateau. However, the difference between the off-rates decreases almost to a factor of 5 as the salt concentration decreases to c s = 5 mM. This is indeed expected since at vanishing salt concentration, neither entropic nor enthalpic effects of the ionic atmosphere are present, and the dissociation kinetics is determined by contact energy between the ligand and binding site. The electrostatic component of this energy increases with decreasing salt concentrations (i.e., the Debye length increases).

On-rates only slightly differ for explicit and implicit salt
Our simulation trajectories also allow us to compute the ligand on-rate k on , allowing us check if the implicit electrostatics affects binding kinetics. We achieve this by considering the part of the trajectories, in which binding and unbinding events reach their quasi-steady states so that the on-rate can be expressed via a detailed balance expression as where n b and n f are the number of bound and free ligands at time t > 1/k off , respectively. In Fig. 3, we show the calculated on-rates as a function of the excess salt concentration. For comparison, we also add the fit curves representing the off-rates in Fig. 2b for implicit and ex- plicit cases. Unlike the off-rates, the on-rates do not show significant difference between the implicit and explicit treatments of the ions. According to Eq. 8, this requires that on average there should be more ligands bound in the explicit case. As we show in the inset of Fig. 3, the ratio n b /n f differ between the implicit and explicit cases; there are more bound ligands in the explicit case, but the corresponding off-rates are smaller compared to the implicit case. Therefore the change from explicit to implicit ions does not only change the unbinding kinetics, but also changes the equilibrium binding site occupation probability, and therefore the binding free energy. In Fig. 3, we have lower error bars particularly at high salt concentrations, for which the sampling is better due to the higher frequencies of the unbinding events (i.e., shorter lifetimes). For the same reason, the error bars in Fig. 3 at low salt concentrations are larger due to the computational cost of producing data when the electrostatic interactions are stronger. In experiments where protein unbinding from single binding sites was studied, a high concentration of free ligands (proteins) weakened the salt-dependence of the offrates [20,22,27]. In our simulations, we also investigated the effects of the electrostatic scheme on FD by introducing a fixed concentration of unbound ligands in the solution (see Fig. 1). In doing so, interplay between univalent ions and free cationic ligands (with valency +4) on the off-rates can be revealed. In Fig. 4, we show the off-rates of the ligands as a function of salt concentration for various free ligand concentrations, c p , for the cases where we model electrostatics implicitly and explicitly. Fig. 4 reveals that, regardless of the electrostatic treatment, if the number of free ligands is increased, the off-rates become less sensitive to the changes in the salt concentration. This result is in accord with previous experiments [27]. As c p is increased, there are more ligands competing for the binding sites, hence, the FD effects are more dominant. According to the data presented in Fig. 4, this increasing competition weakens the salt-dependent steps of the unbinding process significantly [27]. In Fig. 4a and b, as c p is increased, the slopes in logarithmic k off versus c s plots decrease and eventually at high enough c p values, respective plateau regimes appear. This behavior is qualitatively in agreement with the previous experiments [27]. However, in the explicit case (Fig. 4b), the plateau regime appears at smaller ligand concentrations with much slower off-rates compared to the implicit case (Fig. 4a); almost one order of magnitude more ligand concentration is needed in the implicit simulations to reach the salt-independent off-rates (for instance, compare open squares in Fig. 4a and pink circles in Fig. 4b). These results suggest that for the FD effect to take place, an increase in steric interactions on the binding site is essential [27,34]. In the explicit case, the presence of salt ions contribute this effect further by possibly increasing ion-binding site interactions [17], and thus, the plateau regime emerges at lower ligand concentrations compared to the implicit case.
A quantitative comparison of the off-rates obtained from the implicit and explicit simulations in Fig. 4 shows that at similar ligand concentrations, the off-rates for both cases are quite similar at low salt limits. For instance, at around c p ≈ 3 × 10 −4 M, the off-rates are near k off ≈ 10 −3 1/τ . This again suggests that once ligand concentration is high enough, they can replace the low concentration of salt ions near the binding sites and lead to similar kinetics rates for both cases.

Facilitated dissociation is slower for explicit ions than for implicit ions
The effect of electrostatic scheme on the FD can be seen more clearly in Fig. 5, where we show the offrates for implicit and explicit simulations as a function of ligand concentration at a salt concentration of c s ≈ 100 mM. At low ligand concentrations, the presence of explicit ions considerably decreases the off-rates by almost one order of magnitude, as seen in Fig. 5. As the free ligand concentration, c p , is increased, deviations from the spontaneous dissociation rates (indicated by vertical dashed lines in Fig. 5) in both explicit and implicit cases become more drastic. This can be seen more clearly in the inset of Fig. 5, where we show the same data sets but rescaled by the values of the spontaneous dissociation rates. The data in the inset demonstrates that in both cases, the deviations from the respective spontaneous-dissociation rates gradually increase. Also note that the deviations obey different slopes, thus, it is not possible to obtain a master curve for the two sets of off-rates.
In Fig. 5, as c p is further increased towards c p ≈ 10 −3 M, we observe that the difference between implicit and explicit treatment of electrostatics enters a decreasing trend. This behavior of the off-rates may be considered as result of a transition from dilute to semidilute solution with increasing ligand concentration. In semidilute regime, excluded volumes of ligands can overlap even in bulk. We can calculate the critical concentrations, c * p , for the ligands in various conditions. For our solution of ligand chains composed p = 4 beads, c * p on the scaling level is c * p ≈ p/(pσ) 3 ≈ 10 −2 M with σ = 0.7 for a strongly stretched conformation. In good-solvent conditions, c * p ≈ p/(p 3/5 σ) 3 ≈ 10 −1 M again with σ = 0.7 nm. Both of these critical concentrations are much larger than the concentration, at which the off-rates for implicit and explicit cases meet in Fig. 5 (i.e., c p ≈ 10 −3 M). Thus, the convergence is not due to crowding of ligand chains in bulk phase.
One explanation for proximity of the off-rates in Fig. 5 at around c p ≈ 10 −3 M for the implicit and explicit cases is that as the number of ligands competing for the same site increases, molecular overcrowding around the binding site by the ligands can deplete the ions from the binding site, and thus, lead to similar off-rates regardless of the electrostatic treatment. If we calculate the radial distribution functions between the site and the ligand at c p ≈ 10 −3 M (see SI Figure S2), indeed visually there is almost no difference between the two cases; the binding site is equally crowded by the ligands in both implicit and explicit cases at high ligand concentrations.
Interestingly, in Fig. 5, at c p > 10 −3 M, the difference between the off-rates obtained for the implicit and explicit cases tend to increase again. While the onset of the saturation is clearer in the explicit simulations (this result will be supported by additional data in the following section), for the implicit case, the saturation (if any) does not emerge clearly at the concentration range we study here. We note that at very high ligand concentrations, the ligands themselves can act like multivalent ions. Combined with explicit univalent ions, these ligands can yield saturation values lower than those in the implicit case, similar to that observed in Fig. 2b.
We underline that as the concentration of a polyelectrolyte solution is increased to c p ≫ c * p at high-salt conditions, the solution properties, such as effective viscosity, can change drastically [45]. Thus, at higher ligand and salt concentrations that we did not consider here due to their low relevance to biological systems, the effects of implicit and explicit treatments on the off-rates can be more complex.

Short-ranged non-electrostatic interactions modulate the rate of facilitated dissociation
In the previous subsections, we showed that use of explicit ions versus implicit ions modifies the unbinding kinetics by considering various salt and ligand concentrations. Those computations were done at a fixed value of non-electrostatic short-ranged interaction. Biological ligands (proteins) can have specific interactions with binding sites, and that specificity can be varied appreciably, e.g., with DNA sequence. In this subsection we vary the strength of the non-electrostatic interactions and focus on the cumulative effect of electrostatic and non-electrostatic interactions on unbinding kinetics dur- ing FD. In our simulations, we control the non-electrostatic binding energy between the ligands and binding sites by varying the strength of the attractive potential in Eq. 2. This attraction models enthalpic interactions between the ligands and the binding sites (e.g., between proteins and nucleic acid chains). In general, increasing the non-electrostatic attraction requires longer computational times to simulate complete titration kinetics because it shifts off-rates downwards without changing the salt dependence (see SI Figure S3). However, high non-electrostatic binding strengths also enhance the separation of bound and semi-bound states of the ligands. This allows clearer observation of FD [23,24].
In Fig. 6, we show the simulation results for the explicit electrostatic simulations only for various nonelectrostatic binding energies, since we have shown in the previous sections that implicit treatment cannot provide correct physical environment for unbinding. In the simulation, we vary the strength of the non-electrostatic binding energy per bead so that we can scan energy ranges around 10 k B T total, which is the typical molecular binding energies [10,34,46]. Our simulations reveal various regimes for unbinding events.
As can be seen in Fig. 6, at low ligand concentrations, below c p = 10 −5 M, the off-rates are weakly dependent on the concentration of free ligands in solution. In this regime, the time needed for a bound ligand to spontaneously dissociate is relatively little compared to that for a free ligand to bind and destabilize the complex at the binding site. Interestingly, at intermediate concentrations of excess ligands, there is a distinctive regime, whose slope strongly depends on the non-electrostatic in-teraction strength between the ligand and the binding site (Fig. 6); as the non-electrostatic binding strength is increased gradually from E 0 = 1.5 k B T to E 0 = 3.0 k B T , the c p -dependence of the off-rate becomes steeper in the concentration range c p ≈ 10 −5 − 10 −3 M. As the concentration approaches c p ≈ 10 −3 M, the off-rates for various strengths increase but tend to saturate at comparable values. A similar trend also emerges between the implicit and explicit simulations in Fig. 5 (recall that all data in Fig. 6 are obtained via explicit simulations). Analogous universal behavior was proposed in a computational study of neutral ligands for a broader range of concentrations [24]. We observe a convergence only at very high ligand concentrations in Fig. 6 for various binding strengths.
To quantify the binding isotherms in Fig. 6 obtained for various values of E 0 , we fit the data sets in Fig. 6 to an equation in the form of where A, B, D are the fit parameters to determine, and they correspond to low and high concentration limits in a few-state binding model (cf. [27]). The exponent m usually takes a value near unity for the familiar Langmuir binding isotherm in the absence of any cooperativity (i.e., binding events are not altered by already bound ligands) [23,27]. The dashed curves in Fig. 6 indeed are the fit functions obtained from Eq. 9 by settting m = 1 (see SI Table 1 for the fit parameters). For the binding strengths E 0 ≤ 2.5 k B T , Eq. 9 can describe both low and high concentration plateaus as well as the transition regimes successfully. However, for E 0 = 3.0 k B T , setting m = 1 in Eq. 9 does not lead to a reasonable fit. If we release m as a fit parameter in Eq. 9, the data set for E 0 = 3.0 k B T can also be fit to Eq. 9 but with a value of m ≈ 2, whereas other data sets can still be described by m ≈ 1 (see SI Table 2 for the fit parameters). The fit value of m > 1 for E 0 = 3.0 k B T in Fig. 6 suggests that high binding affinity can increase the binding of other ligands and increase cooperativity. With increasing binding affinity, the ligand bends the polymeric binding site stronger. Indeed, in the simulations, we observe this behavior qualitatively. Strong affinity towards the binding site, combined with local bending, can decrease rotational and translational mobility of the ligand on the binding site [47]. Thus, a newly attached ligand has a higher chance to bind to the site with minimum interference from the bound ligand.
Eq. 9 also provides valuable information for the low and high concentration limits, D/B and 1/A, respectively. While the high concentration saturation rate, 1/A, decreases with increasing E 0 , the variation among different strengths is not more than factor of 3 in Fig. 6. Consistently, the difference in 1/A's of E 0 = 3.0 k B T and E 0 = 2.5 k B T is rather insignificant (see SI Table  1  ues of D/B range between ≈ 10 −6 and 10 −1 1/τ (see SI Tables 1 and 2), suggesting exponentially decreasing off-rates with increasing E 0 at c s = 100 mM (i.e., D/B ∼ exp(2E 0 /k B T )) [23]. Interestingly, the effect of non-electrostatic interactions at the intermediate ligand concentrations is also evident in the implicit electrostatic simulations with the systematic overestimation of the off-rates (see SI Figure  S4). The similarities in the off-rates in the explicit and implicit as a function of solvated ligands confirm that steric interactions between two or more ligands on the binding site are essential for FD process.

IV. CONCLUSIONS
In this study, by using a generic polyelectrolyte model, we examined the effects of implicit (mean field) versus explicit treatments of ionic components on the dissociation kinetics of single ligands from polymeric binding sites. We studied the role of electrostatics in both spontaneous and facilitated disassociation mechanisms at various salt and excess ligand concentrations. Our simulations revealed that a mean-field level treatment of electrostatics overestimates the unbinding rates compared to explicit treatment of ionic species (e.g., by considering their net charge and steric interactions). The explicit ions provide translational entropy and ionic correlation contributions, both of which favor slower unbinding rates. Furthermore, the trend in unbinding rates versus the ligand concentration data does not allow us to obtain a master curve, in which the implicit and explicit simulation data could be overlapped. This suggest that the explicit ions change the underlying physics of the FD. We now discuss our major results in slightly more detail.

A. Spontaneous dissociation kinetics are significantly different in explicit versus implicit ion cases
Ionic effects are expected to play major roles in the unbinding kinetics of molecular ligands. First of all, with increasing salt concentration, the Debye length decreases. The ionic correlations in bulk are more correlated at length scales smaller than the Debye length. At scales larger than the Debye length, Coulomb interactions are screened. Secondly, entropic and enthalpic contributions of ionic atmosphere near the binding site and ligand can favor or inhibit dissociation, depending on the salt concentration. In our implicit simulations, only the screening effect of salt is available via manipulation of the Debye length in the pair potential given in Eq. 5. In our explicit simulations, where all the univalent ions are considered explicitly, charge-related and steric (excluded volume) effects of the ions are present. Hence, in addition to the screening effect, in the explicit case, the ions arrange themselves near the binding site and the ligand relative to the bulk by balancing corresponding entropy and enthalpic energy components. Below we discuss these components and their role in the unbinding processes along with our simulation results.

Salt dependence of unbinding is determined by electrostatic screening at low salt
In the implicit simulation, the binding affinity has two contributions; a non-electrostatic attraction characterized by E 0 , and an electrostatic Coulomb attraction between the oppositely charged binding site and the ligand. These two contributions favor a bound ligand. The nonelectrostatic binding energy has no salt dependence (we assume that the salt does not change chain conformations drastically). Contrarily, the electrostatic attraction becomes weaker with increasing salt concentration. As the electrostatic part of the binding affinity decreases with increasing salt concentration, thermal fluctuations promote (spontaneous) dissociation of the ligand at c s < ∼ 100 M. Thus, the slope of ≈ 1.4 that we observe in the off-rate versus salt concentration plots below physiological salt concentrations in Fig. 2b can only be due to the weakening of the electrostatic pairwise interactions (Eq. 5). At c s > 100 M, increasing salt concentration weakly change the electrostatic part of the attraction since κ −1 < 1σ. As a result, the non-electrostatic binding strength, E 0 , controls the unbinding. It is important to note that ion concentration per Debye volume scales as κ 3 ∼ c 3/2 s , which is somewhat close to the exponents observed in Fig. 2b.
At this point, we remind reader that on our polymeric binding sites, the 1D charge density is set to 1/ℓ B (i.e., condensed ions are not considered explicitly) for a proper comparison between explicit and implicit electrostatic treatments since the screening effect of "condensed" ions is not well defined [48]. We anticipate that adding "condensed" ions in the simulations would change the slope of the off-rates. In that case, the main contribution to the ion-release effect would be mainly due to "condensed" ions rather than relatively rapidly diffusing ions near the binding site, which we already include in our explicit-ion simulations. This suggests that in studying the kinetics of ligand unbinding, the contributions from condensed and diffuse ions should be considered separately to comprehend ionic effects in more detail.

Quantitatively lower unbinding rates in explicit simulations underline transnational entropy and correlations effects of ions
In the explicit simulations, the non-electrostatic binding affinity and electrostatic the Coulomb attraction also promote a bound ligand. The quantitatively lower offrates obtained in the explicit simulations relative to those in the implicit cases (Figs. 2 and 4) suggest that additional effects due to the explicit ions extends the lifetime of the bound ligand on the binding site.
On a naked binding site, an ion cloud interact with the binding site by forming transient ion pairs. Once a ligand binds to the site, the ionic cloud is distorted by the ligand. If the bound ligand compensates the loss ion pairs energetically, the corresponding energy change between bound and unbound states can be ignored [3]. If this is not the case (i.e., the loss ions pairs are not compensated upon binding), the unbound state can be favored by the interactions between the surrounding ion cloud and the binding site [18]. The strength of this effect intensifies with increasing salt concentration since more ions can interact with the binding site. Note that in the implicit simulations, there are no ions to compensate electrostatic energy loss once a ligand desorbs at low excess ligand concentrations.
In the explicit simulations, the ionic correlations near the binding site (and ligand) are "tighter" than in the bulk as a due to the electrostatic potential induced by the negatively charged binding site. The free energy associated with these ionic correlations favors a bound ligand [17], since dissociation of the bound ligand relocates more number of ions from bulk to nearby the binding site. This effect grows with increasing salt concentration since more ions must arrange themselves near the binding site once the ligand unbinds [18], Another effect that can slow down the off rates in the explicit simulations is the translational entropy of the ions. Ideally ions tend to diffuse away from the structural charges to increase their translational entropy. A bound ligand disperses a number of ions near the structural charges to bulk. Contrarily, upon unbinding, the ions need to replace the unbound ligand by sacrificing their translational entropy. Thus, the translational entropy gain of the ions in bulk works in the advantage of the bound ligand. This effect weakens as the salt concentration becomes more uniform throughout the solution [18].
Overall, the Coulomb attraction between the binding site and the ligand is present for both explicit and implicit cases and favors a bound ligand. In the explicit-ion simulations, both translational entropy and correlation contributions also favor a bound ligand, and they are absent in the implicit simulations. These two ionic contributions in the explicit case extend the lifetime of the ligands on the binding site, and thus, lead to much lower off rates. Interestingly, with increasing salt, the translational entropy contribution decreases whereas the contributions due to the correlations increases [18]. Moreover, in the explicit simulations, the interactions between the ion cloud and the structural charges can also favor or inhibit a bound state depending on the balance of the electrostatic energy upon unbinding. It is possible that these ion-related contributions cumulatively lead to the behavior observed for the ratios of the off rates in the inset of in Fig. 2b obtained for the implicit and explicit cases.

B. Interplay of ions and free ligands in facilitated dissociation
In the FD simulations, at low ligand concentrations the mechanism mentioned above for the spontaneous dissociation still controls the unbinding kinetics. At intermediate excess ligand concentrations, ideally a bound ligand is replaced by other ligand that has been previously in bulk. As a result, the energetic contributions we discuss above should not change as long as the ligands are identical. In this case, the ligand dissociation is determined by the repulsive and steric interactions between two (or more) ligands on the same binding site. As the non-electrostatic binding affinity, E 0 , is increased, it is harder for the bound ligands to kick one another out of the binding site (Fig. 6).
For both implicit and explicit cases, we obtain saltindependent regimes as the excess ligand concentrations is increased (Fig. 4). This result is consistent with the previous experiments, where the salt-dependence of the off-rate of the DNA binding proteins becomes weaker with increasing protein concentration in solution [27]. However, for the implicit case, one order of magnitude more ligand concentration is required to observe this regime (the off-rates in the implicit case are still faster compared to those in the explicit case). This suggests that the salt-independent regime requires a highly competitive binding site environment provided either by the ions or excess ligands competing to replace the bound ligand. In the explicit simulations, excluded volumes of the ions near the binding site enhance this effect by possibly promoting more partially bound ligands. These partially bound ligands are more prone to FD. For the implicit simulations, however, the excluded volumes of the ions are not present, and thus, a higher number of ligands per binding site is needed to observe the FD effect.
At high ligand concentrations, in Fig. 6 the off-rates exhibit a tendency towards a common saturation plateau, regardless of the strength of the non-electrostatic binding, E 0 . Such an effect requires that any contribution due to the non-electrostatic binding energy separating the bound and unbound states weakens or disappears. Such a scenario can arise if the local volume fraction of the monomers near the binding site reaches a value of unity (i.e., the limit of polymer melts). Under this condition,the two-body interactions between the monomers can be screened out (i.e., second virial coefficient approaches zero) [49]. In other words, the solvent quality changes at the local level and renders the nonelectrostatic binding strength irrelevant. Indeed, a comparison of the radial distribution functions between the binding sites and the ligands shows that the local ligand concentrations on the binding site for E 0 > 1.5 k B T approach each other at high enough values of c p (see SI Figure S5).
C. Explicit nature of ionic atmosphere is essential for modeling unbinding kinetics In the salt concentration range we consider (i.e., c s = 5 − 1000 mM), the distribution of explicitly modeled ions around polymeric binding site can lead to non-uniform screening effects. The small ions can destabilize the complex non-trivially and lead to local electrostatic screening. Furthermore, the bound cationic ligands alters the electrostatic potential of the binding site by bringing coions close to the binding site. Such effects cannot be described in the implicit electrostatic simulations, since the mean-field nature of Eq. 5 provides an electrostatic screening regardless of the conformation and relative position of the components (i.e., the Debye length is a function of total salt concentration).
In the on-rates we calculate via detailed balance (Fig. 3), we observe that the difference between the implicit and explicit cases is rather small. This suggest that the large difference in the off-rates reported here is not a simple consequence of procedure of coarse-graining ionic components. The inset in Fig. 2, where we show the ratio of the off-rates obtained via implicit and explicit simulations, also supports this argument, since this ratio is neither a simple rescaling constant nor a linear function of the salt concentration.

D. Implicit simulations can underestimate contact energies
In the implicit simulations, the contact energy between two beads (forming either the ligands or the binding site) at contact can be weaker due to meanfield nature of Eq. 5. This is because, by definition, Eq. 5 gives a rescaled contact (electrostatic) energy of ℓ B /σ exp(−κσ) < ℓ B /σ (recall that ℓ B /σ is the electrostatic contact energy rescaled by k B T between two beads of identical sizes of 1 σ at zero salt). The lower contact energy can also accelerate unbinding To test how the electrostatic contact energy changes the off-rates, we perform additional test simulations by increasing the contact energy in the implicit simulations. We achieve this by shifting the decay term in Yukawa potential (i.e., ℓ B /r exp(−κ(r − σ) = ℓ B /r) so that at r = σ, the implicit and explicit simulations give similar contact energies (here we ignore the effects of the ions surrounding the ligand and binding site). In those simulations, we obtain off-rates that are lower than those in the explicit-case (see SI Figure S6). While this test simulations can demonstrate the effect of the contact energy in unbinding, an arbitrary rescaling of the mean-field potential between the ligand and binding site is necessary to exactly match the results of the explicit-ion simulations.

E. Non-uniform dielectric constant can accelerate the binding
In our coarse-grained MD simulations, we use a continuum solvent model with a constant background dielectric constant. In a more realistic, orientation of water molecules near the binding sites, polar groups composing the ligand itself, and the salt concentration would influence the local dielectric constant and thus, ionic distributions [50,51]. Indeed, recently it has been shown that ionic mobility in salt solutions can be described more accurately with a non-uniform dielectric constant [52]. Similar effects can alter unbinding rates by pushing the ions from the binding site or increasing diffusion coefficient within the ionic atmosphere. Nonetheless, these and similar effects can enhance the importance of explicit treatment of ionic atmosphere in studies, where kinetics of biomolecules are considered.

F. Distribution of binding sites and chromatin in vivo
In the setup we consider, the binding sites arranged on a 2D surface. A 3D distribution of binding sites is expected in the case of protein specific binding sites along chromatin confined by nuclear envelope or bacterial cell membrane. Since, as we show here, the off-rates depend on the local interactions and ionic distributions around binding sites, we do not expect that the distribution of binding sites will change the results we present here as long as the distance between the binding sites is larger the Debye length. However, in a highly concentrated polyelectrolyte matrix, electrostatic potentials of constituting chains can broaden the ion distributions around neighboring chains [53]. In such cases, exchange rates can be a function of local polymer conformation [54].

G. Quantitative modeling of protein-dsDNA interactions
Our calculations have been motivated by consideration of experiments with proteins binding to dsDNA binding sites. One might ask how accurately our MD model can be used to physically describe the very high charge density of DNA. In our simulations, we set the size of each bead composing the binding sites and ligands to 0.7 nm. By setting the size of the beads to ≈ 1.4 nm and their charges to −2 without changing the size of the ions, a model mimicking binding sites along dsDNA is obtained. The simulations with those parameters do not lead to any drastic difference. At low salt concentrations, the offrates have a weaker dependence on salt (see SI Figure S7). In general, we should note that models like ours can be used to understand the fundamental physical principles of electrostatic effects. However, atomistic or less coarsegrained computational models are more suitable to study molecular aspects of DNA binding proteins [55,56]. We plan to use one of these models to further study unbinding kinetics in the near future.

H. Final remarks
Our results indicate that condensed ions can have stronger effects in determining the slope of dissociation constant in the low-salt limit since we obtain qualitatively similar behaviors for our unbinding rates in the implicit and explicit simulations, in which only the effects of uncondensed charges are introduced a priori. The contributions from condensed and diffuse ions should be considered separately to infer ionic effects in unbinding kinetics. Similarly, in our facilitated dissociation simulations, the implicit treatment led to higher unbinding rates for various excess ligand concentrations at all salt concentrations we considered. The salt concentration has a much weaker effect on the facilitated dissociation due to depletion of ions from binding sites. At high ligand concentrations, at which binding site is overcrowded by ligands, we observe a universal saturation regime of facilitated dissociation, regardless of non-electrostatic binding strength or electrostatic treatment. This leaves the possible role of nonidentical ligands in unbinding process as an open question for further investigations. Finally, our results also suggest cooperative binding of ligands for high non-electrostatic binding strengths, suggesting a general mechanism whereby any protein bound to a binding site can promote binding of a second protein near the same site, even without direct interactions of the adjacent proteins.