The crowded cellular environment can affect biomolecular binding energetics, with specific effects depending on the properties of the binding partners and the local environment. Often, crowding effects on binding are studied on particular complexes, which provide system-specific insights but may not provide comprehensive trends or a generalized framework to better understand how crowding affects energetics involved in molecular recognition. Here, we use theoretical, idealized molecules whose physical properties can be systematically varied along with samplings of crowder placements to understand how electrostatic binding energetics are altered through crowding and how these effects depend on the charge distribution, shape, and size of the binding partners or crowders. We focus on electrostatic binding energetics using a continuum electrostatic framework to understand effects due to depletion of a polar, aqueous solvent in a crowded environment. We find that crowding effects can depend predictably on a system’s charge distribution, with coupling between the crowder size and the geometry of the partners’ binding interface in determining crowder effects. We also explore the effect of crowder charge on binding interactions as a function of the monopoles of the system components. Finally, we find that modeling crowding via a lowered solvent dielectric constant cannot account for certain electrostatic crowding effects due to the finite size, shape, or placement of system components. This study, which comprehensively examines solvent depletion effects due to crowding, complements work focusing on other crowding aspects to help build a holistic understanding of environmental impacts on molecular recognition.

## I. INTRODUCTION

The cellular environment is crowded, heterogeneous, and dynamic, containing a variety of biomolecules and metabolites that occupy up to 40% of the cellular volume.^{1,2} Macromolecular crowding can impact protein folding,^{3} binding,^{4} and kinetics,^{5} with such effects having broader implications, such as impacting the formation of fibrils linked to neurodegenerative diseases.^{6–8} Crowding effects were originally thought to be mediated primarily through “volume exclusion,” whereby crowders reduce the available volume for systems of interest.^{9–11} In such models, it is entropically favorable for systems to assume compact states.^{12–15} More recently, enthalpic interactions between systems and crowders have been studied. These “soft” interactions, sometimes referred to as the “quinary” structure, can counteract the volume exclusion effect in certain cases, demonstrating the importance of accounting for the chemical behavior of crowders.^{16–22}

Crowding effects on solvent properties can also be important.^{23–25} First, crowding causes altered water mobility, as hydration shell water differs in its properties from bulk water.^{26–28} Second, and a focus of this paper, crowding can cause water depletion, leading to changes in electrostatic screening and desolvation energetics upon binding.^{29}

As many entropic and enthalpic mechanisms can underlie how crowding can affect systems, crowding effects in a particular scenario may depend on physical and chemical properties of the crowders, such as charge, shape, flexibility, and size.^{30–39} Additionally, the physicochemical properties of the system being crowded may also affect the ways in which it is impacted by crowding.^{40–48} There is thus utility in exploring how the characteristics of the biomolecular system of interest, the crowders in its environment, and the solvent could lead to differences in crowding effects.

Experiments have probed various contributions of crowding effects, and both “hard core” repulsions and “soft” chemical interactions continue to be modeled.^{49} For example, phase diagrams can predict the relative enthalpic and entropic contributions to crowder-induced changes in a process free energy as a function of simple interaction potentials between system components.^{50} Models, from scaled particle theory to those built on the depletion forces of the Asakura and Oosawa model,^{51} have brought the field to consider how trade-offs between excluded volume and other crowding effects depend on system properties. Nevertheless, there is also value in focusing more exclusively on the details of these “other” crowding effects through a detailed exploration of how their contributions depend on system properties.

Here, we aim to provide a more comprehensive understanding of some of these “other” crowding effects—namely, electrostatic and solvation effects—with an exclusive focus on binding energetics, which can be affected to differing extents by crowding.^{52–58} Our goal is to systematically characterize solvent-mediated electrostatic binding energetics in crowded systems as a function of physicochemical properties of the binding partners and surrounding crowders. Rather than providing system-specific insights into these effects as previous studies have done,^{21,29,45,59,60} we aim here to provide broader, generalizable insight, intuition, and understanding through employing theoretical model systems. Although idealized molecules clearly deviate from chemical reality, they allow us to separate physical properties from chemical constraints for greater control, better sampling, and efficient investigation. Similar to previous work that explored electrostatic specificity and affinity in a dilute (uncrowded) solution,^{61} the idealized models here are built from overlapping, spherical “atoms” and can be easily manipulated to vary physical properties. By using abstracted models to evaluate the solvent-mediated electrostatic effects of crowding as a function of the charge distribution, size, and shape of the binding partners and crowders (Fig. 1), we aim to create a framework to more holistically understand electrostatic interactions in the presence of crowded, solvated environments.

## II. METHODS

### A. Idealized model systems

Idealized models were constructed by overlaying “atoms” of radius 2 Å, with a distance of 1.5 Å separating adjacent centers, to create an initial rectilinear grid. Each atom was then offset in a random direction and magnitude of up to 0.5 Å to break the symmetry of the grid. See Figs. 2(a), 4(a), and 6(a) for images of sample complexes. The sizes and shapes for each model complex varied by the physical property isolated (details are given in Table SI in the supplementary material). The total complex volume was kept roughly constant within each physical property for maximal control. Dummy atoms were placed at identical corner points within models directly compared such that the grid resolution (grids/Å) for continuum electrostatic calculations was constant (see Table SII in the supplementary material for details). Coordinates for the maximal possible area of interfacial atoms were kept identical across compared models (Table SII).

Each model atom could hypothetically take on any point charge value. In most cases, to sample a wide variety of charge distributions, four interfacial charges on each binding partner were varied from −1e to 1e in increments of 0.5e (used for unweighted averages; see Sec. II E) or 1e (used for Boltzmann-weighted averages), yielding either 625 or 81 different charge distributions for each partner [see bottom left panel of Fig. 2(a) for an image showing a typical set of four interfacial atoms considered before their random offset to break the symmetry]. In some specified cases, additional atoms were also charged by randomly selecting atoms throughout each partner biased by probabilities based on the density of buried and solvent-exposed charged amino acids in a dataset of protein structures.^{62}

### B. Crowding models

#### 1. Explicit crowders

Crowders of a predefined shape were randomly placed around the model complex in a fictitious box. To define the box, dummy atoms were placed 50 Å in each dimension from the most extreme *x*-, *y*-, and *z*-coordinate of each model complex (see Table SII). Unless otherwise stated, the crowders occupied ∼30% of the box volume not occupied by the ligand/receptor/complex itself. Generally, spherical crowders of varied radii randomly sampled between 5 and 20 Å were used. When comparing spherical and rod-like crowder shapes, each spherical crowder had a radius of 8 Å; rod-shaped crowders consisted of five overlaying spherical atoms of radii 5.19 Å, each separated by a distance of 5 Å, leading to an overall length of 30.37 Å. Volumes of spherical and rod-shaped crowders were ∼2145 and 2141 Å^{3}, respectively, differing by <0.5%. When spherical crowders were charged, 6 points within each crowder, each 1 Å away from maximal or minimal *x*, *y*, and *z* extents of the sphere, were charged to +1/6 e or −1/6 e. Unless otherwise noted, bound and unbound states were crowded separately and independently for each trial and 200 such trials were carried out for each system.

#### 2. Implicit crowding

As an implicit model of crowding to account for water depletion in an average way,^{27,59,63} the solvent’s dielectric constant was lowered to 57, approximately a 70/30 weighted average of the pure aqueous solvent and macromolecular dielectric constants of 80 and 4, respectively.

### C. Calculation of electrostatic binding free energies

To compute the electrostatic binding free energy between partners with a given, fixed charged distribution, the Linearized Poisson–Boltzmann Equation (LPBE) was numerically solved for the solvated bound and unbound states, with charges aligned in each state such that artificial grid energies canceled out when taking differences. To enumerate binding free energies, desolvation penalties, and solvent-screened interactions across a combinatorial range of charge distributions for a given trial, each charge center on the ligand or receptor to be varied was charged to +1, in turn, and the potential was solved for the bound and unbound states via the LPBE and used to construct unit potential matrices, described in detail in previous work.^{64–66} The electrostatic binding free energy can be calculated for any charge distribution via the equation

where *q*_{L} and *q*_{R} are vectors whose elements are the atom-centered point charges on the ligand and receptor, respectively, and *L*, *R*, and *C* are unit potential difference matrices. Specifically, *L _{ij}* or

*R*is one-half the potential difference between the bound and unbound states at charge center

_{ij}*j*when charge center

*i*is set to +1e.

*C*is the potential at charge center

_{ij}*j*on the receptor in the bound state when charge center

*i*on the ligand is set to +1e.

The first term in Eq. (1) quantifies the ligand desolvation penalty (LDP), or the energetic cost of replacing a receptor-shaped region of the high dielectric solvent with a low-dielectric cavity. Physically, it is the free energy required to remove the necessary aqueous solvent from the ligand interface upon binding. The second term quantifies the analogous process for the receptor, the receptor desolvation penalty (RDP). The final term quantifies the bound-state solvent-screened interaction between the partners (INT).

To calculate free energies of binding in the presence of charged crowders, a thermodynamic cycle was used, in which the total binding free energy was calculated as

which accounts for the desolvation of each crowded unbound state ($\Delta Gdesolv,unbound\u25e6$), the removal of Coulombic interactions between each unbound partner and its crowders ($\Delta Gcoul,unbound,crowders\u25e6$), the Coulombic interaction between the (now uncrowded) partners as they come together ($\Delta Gcoul,partners\u25e6$), the Coulombic interaction between the complex and its crowders upon their introduction ($\Delta Gcoul,bound,crowders\u25e6$), and the solvation of the bound, crowded state ($\Delta Gsolv,bound\u25e6$). In order to correct for the fact that the number of crowders is not kept constant between the unbound and bound states, the last two terms represent the subtraction of a fictitious cycle in which the difference in solvation energies is calculated for only the crowders (in the same positions they are in for each bound and unbound state). This method does not explicitly account for direct interactions between the crowders as a function of their different placements in the bound and unbound states. Doing so may require keeping the total number of crowders in bound and unbound states constant by modeling both partners separated by a great distance in crowded unbound states; such length scales are prohibitive for the current study. For Coulombic calculations, a dielectric constant of 4 was used. See supplementary material, Fig. S1 for a visual diagram of the cycle.

### D. Continuum electrostatic calculations

To calculate binding and solvation free energies or to obtain the unit potentials for each charge to generate the relevant matrix elements, the Linearized Poisson–Boltzmann Equation was solved using the DelPhi^{67} software package (version 8.4). To generate molecular surfaces, a probe radius of 1.4 Å was used. The dielectric constants of model binding partners and crowders were set to 4 with the solvent dielectric constant set to 80 unless lowered to implicitly model crowding effects. A 401 × 401 × 401 grid was used for the numerical solution of the LPBE with a focusing procedure in which the longest system dimension occupied first 23% and then 92% of the grid. This yielded grid densities between 2.2 and 3.0 grids/Å (systems that were directly compared had identical dummy atom placement to ensure identical grid spacing). The ionic strength was set to 0.145 M, the temperature to 310 K, and the Stern layer to 2 Å. Energetic components were computed from matrices using MATLAB (The Mathworks, Inc., Natick, MA).

### E. Averaging, Boltzmann weighting, and statistical analyses

Bound and unbound state potential matrices were computed for 200 trials for each system. In general, we paired the (independent, unless otherwise noted) bound state with the unbound states and the unbound states with each other within each trial. To provide insight into sampling effects, we averaged over trials in two ways using MATLAB (The Mathworks, Inc., Natick, MA). First, we performed simple unweighted averages, with standard errors showing the statistical variation in average energetic components due to random crowder placement. *p*-values were estimated by calculating standard errors between differences of samples. When differences were paired, the standard error was computed on the pairwise difference. Otherwise, the variance in the mean of the sum was computed by summing the variances in the mean for each sample. A result is considered significant to *p* < 0.05 if the sample difference is at least twice the standard error of their average difference.

We also performed Boltzmann-weighted averages over each set of trials using a temperature of 310 K. When crowders were uncharged, the bound and unbound total electrostatic grid energies for each set of crowder placements were used to weight that trial’s contribution for each component (LDP, RDP, INT) toward the binding free energy. We assumed that the electrostatic free energy roughly equals electrostatic energy and used separate unit potential matrices for the bound and (paired) unbound states for each trial, ultimately taking the difference of their separately weighted averages across all trials. To preserve the additivity of LDP, RDP, and INT toward the overall binding free energy, the additional entropic difference due to the energy spectrum created by 200 trials was neglected. The robustness of Boltzmann-weighted averages to the arbitrary pairing of unbound states was evaluated through re-calculating averages using partition functions of unpaired unbound states, and the robustness to finite sampling was estimated by calculating standard errors via bootstrapping for a subset of our results; quantitative averages for uncharged crowders were generally robust to within ∼0.1 kcal/mol or less with bootstrapped standard errors also generally ∼0.1 kcal/mol or less, and qualitative results were unaffected for the considered subset.

When crowders were charged, each bound state energy used for Boltzmann weighting was calculated as

using the above definitions for Eq. (2). The unbound state energy was

For charged crowders, where only one complex charge distribution was considered at a time rather than being combinatorially varied, the standard errors of the Boltzmann-weighted averages were estimated through computing the standard deviation over 200 bootstrapped samples; for each sample, 200 crowder placements for each bound and unbound state were sampled with replacement from the 200 unique trials independently, except for when the same crowder placements were considered in the bound and unbound states. Here, we present our observed weighted average along with the standard errors computed from the bootstrapped samples. The observed averages in certain cases noted in the results were not quantitatively robust (different by >1 kcal/mol) to sampling, as seen via bootstrapping, or to pairing of unbound states, suggesting the need for more sampling in future studies requiring quantitative predictions about charged crowder effects. In these cases, the qualitative ranges suggested by these resampling methods are noted in Sec. III.

### D. Figure generation

All model complex and crowding visualizations were generated with Visual Molecular Dynamics (VMD).^{68} Plots were generated with MATLAB (The Mathworks, Inc.).

## III. RESULTS

### A. Electrostatic effects of crowding due to physical properties of the binding partners

We first consider crowding effects on electrostatic binding free energies as a function of physical properties of the binding partners themselves—their charge distributions, shape, and relative size. For each system considered, uncharged, spherical crowders of varying radii were randomly placed at a volume density of 30%, and bound and unbound states were crowded separately for each of 200 trials.

#### 1. Electrostatic crowding effects depend predictably on the charge distributions of the binding partners

To systematically assess how charge distributions of binding partners impact crowding effects on their electrostatic binding free energies, the charges of four indicated interfacial atoms on each binding partner shown in Fig. 2(a) were varied combinatorially in both uncrowded and crowded environments.

We first discuss results from weighting each of the 200 trials of random crowder placements equally before presenting results using Boltzmann weighting. We note that equally weighting the crowded trials may not realistically model “true” crowder spatial distributions, but these results allow us to highlight and build intuition for phenomena that can become obscured or complicated upon Boltzmann weighting. Here, each atom’s charges were varied from −1e to +1e in 0.5e increments, yielding 625 charge distributions per partner. Figures 2(b)–2(e) show the average change due to crowding in the LDP, RDP, INT, and total electrostatic binding free energy as a function of each respective uncrowded quantity. On average, crowding lowers the LDP and RDP in a manner linearly proportional to the original LDP or RDP [Figs. 2(b) and 2(c), *p* < 0.05 for all nonzero charge distributions]. The linear relationship is expected because desolvation penalties are quadratic functions of charge in a linear-response model, so a doubling of all charges should quadruple both a desolvation penalty and a change in desolvation penalty due to crowding. These results generalize previous work that predicted lower desolvation penalties for the barnase/barstar system in a crowded environment.^{29}

Lower desolvation penalties may result from “predesolvation”^{29,59}—the idea that crowders in the unbound state may partially occupy volume ultimately taken up by the binding partner, pre-depleting the interfacial solvent in the unbound state and reducing the amount displaced by the partner upon binding. Indeed, an implicit model of crowding, in which a lower solvent dielectric constant (ε_{solvent} = 57) captures ∼30% water depletion in a spatially averaged sense, also shows lowered desolvation penalties [pink points in Figs. 2(b) and 2(c)], as it accounts for lowered solute–solvent interactions.

Using either explicit crowders or a lowered solvent dielectric constant, solvent-screened interactions were amplified in a manner linearly proportional to their uncrowded interactions [Fig. 2(d)], again generalizing a previous system-specific observation.^{29} Because the solvent screens interactions, its depletion via crowding would descreen—or amplify—those interactions.

Interestingly, because the desolvation penalties are lowered on average but interactions can either be increased or decreased, crowding effects on the overall electrostatic binding free energies produce a more intriguing curve [Fig. 2(e)]. Partners that bind favorably in uncrowded environments have more favorable desolvation penalties and interactions upon crowding, and so their binding free energies are robustly lowered (*p* < 0.05 for >99% of such charge distributions). Unfavorable electrostatic binding free energies are lowered to a lesser extent because of the trade-off between lowered desolvation penalties and amplified repulsions. The same qualitative trend is seen with a lowered solvent dielectric.

Whether interaction or desolvation dominates in this trade-off for unfavorably interacting partners can depend on the locations of the relevant charges. Here, we considered interfacial charges that were buried upon binding; desolvation effects dominated and overall binding free energies were lowered upon crowding for all charge distributions. Preliminary work shows that upon adding non-interfacial charges, interaction effects can dominate;^{69} this idea is also seen in Fig. S2, in which 39 non-interfacial “background” charges on each partner were set to +0.25e, as described in Sec. II, to model unfavorable longer-range interactions. When crowding was modeled as a lowered solvent dielectric and each partner’s interfacial charges were varied as before, many complexes with unfavorable uncrowded binding free energies became more unfavorable (Fig. S2D) due to dominating unfavorable interactions (Fig. S2C) despite less unfavorable desolvation penalties (Figs. S2A and S2B).

Thus far, we equally weighted each crowded trial to build intuition for how energetics depend on crowder positions. Here, to begin to provide insight into what happens when accounting for the unequal sampling of such positions, we show average crowding effects that were Boltzmann-weighted by electrostatics, i.e., assuming the limiting case of electrostatics and polar solvation dominating energetics. With this weighting, desolvation penalties were still lowered, on average, for most charge distribution pairs [Figs. 3(a) and 3(b)], but effects are highly muted, yielding only a total of tenths of kcal/mol at most when combining LDP and RDP effects (for a given ligand/receptor charge distribution, there are now multiple values for changes in LDP/RDP because bound state weights depend on the partner’s charge distribution). This result makes sense because unbound state trials that are more “predesolvated” have higher electrostatic energy than their more solvated counterparts and will be weighted less. Nevertheless, this result shows that predesolvation plays a smaller role in cases where crowder placement is dominated by electrostatics. Additionally, the differing probabilities across crowder placements are not adequately captured by modeling crowding implicitly via a lowered solvent dielectric constant either, and thus, the results in Figs. 3(a) and 3(b) differ from the both the implicit and unweighted explicit results in Figs. 2(b) and 2(c).

Boltzmann-weighted interaction effects are similar to those predicted using unweighted averages [Fig. 3(c)], although less destabilization upon crowding is seen for unfavorably interacting partners, perhaps because unstable bound state configurations with the most desolvation—and therefore the greatest amplification of repulsive interactions—would be weighted less. Because predesolvation effects are smaller, the effects on interactions dominate and the overall binding free energy is slightly increased for complexes with unfavorable interactions [Fig. 3(d)], which is the opposite of what was seen using unweighted averages. This result again highlights the sensitivity of crowding effects to the sampling of crowder placements.

In summary, we see that solvent depletion due to crowding lowers the average electrostatic binding free energies of complexes that interact favorably and could either lower or increase binding free energies of complexes that interact unfavorably, with predesolvation effects showing greater dependence on sampling than interaction effects. The observed change in relative binding affinities between differently charged partners upon crowding implies that crowding could modulate binding specificity. Additionally, because changes in ionic strength often affect interactions more drastically than they affect desolvation penalties,^{61} the charge-dependent interplay between changes in desolvation and interaction due to crowding could be rationally modulated by varying the ionic strength. Finally, electrostatic binding free energies are often calculated to be unfavorable—sometimes larger than 50 or even 100 kcal/mol—between macromolecules, such as proteins.^{70} In such cases, binding is driven by other nonpolar factors. As a result, the electrostatic effects of crowding for complexes with unfavorable electrostatic binding energetics may be important to consider, and it is in these cases where the potential competing interplay between interaction and desolvation components is more relevant.

#### 2. The geometry of the binding interface can affect changes in desolvation energetics upon binding due to crowding

The effect of crowding on electrostatic binding free energies was assessed for two different interfacial shapes, “flat” and “cavity,” as shown in Fig. 4(a). For maximal control, the ligand was kept identical in both cases, and the two complex volumes were within ∼1% of each other (see Table SI). Using either unweighted averages or Boltzmann weighting, Figs. 4(b) and 5(a) show that the different interface shapes do not greatly impact crowding effects on LDPs (although the cavity-bound ligand not surprisingly pays a higher desolvation penalty when highly charged). However, there are differences between the interfaces for changes in RDPs. When crowding is modeled explicitly, the receptor in the cavity system does not show lowered RDPs, on average, due to crowding when using unweighted averages [Fig. 4(c), blue error bars], unlike that in the flat system (*p* < 0.05 for all charge distributions). This trend is also qualitatively seen when assuming Boltzmann-weighted averages, although, again, the quantitative differences are much smaller [Fig. 5(b)]. Predesolvation of the receptor in the cavity case is more difficult as a crowder must partially occupy the deep pocket in the unbound state, which, due to spatial constraints, is less likely than predesolvation for a flat interface. Indeed, an implicit model that cannot account for the finite size of crowders shows the “typical” RDP lowering upon crowding [Fig. 4(c), cyan points].

To demonstrate the relationship between the finite size of crowders and cavity desolvation, we repeated the unweighted average calculations using a fixed crowder size of either 8 Å in radius (too large to fit within the unbound “cavity” receptor pocket) or 4 Å in radius (a size that can more easily enter the pocket). Because of the large number of crowders in the 4 Å case, we crowded to a total volume of ∼20% rather than ∼30% in both cases. Figures S3C and S3D show that while predesolvation does not appreciably occur for the cavity receptor with 8 Å crowders, it does occur to some extent with 4 Å crowders, as there is now a lowering of RDP upon crowding. Interestingly, smaller crowders are also able to better descreen interactions in the bound state for the flat interface (Figs. S3E and S3F), likely because such crowders can more closely approach interfacial charges in the bound state. In the cavity case, interfacial charges are more buried and inaccessible to the solvent in the bound state. This shows that the relationship between the crowder size and interfacial geometry can have subtle but potentially consequential effects on both overall electrostatic binding energetics and the interplay between various contributions toward these overall energetics. This interplay could be capitalized upon to modulate molecular recognition in different environments.

In the explicit crowding model, the flat INT was slightly more enhanced by crowding than the cavity-interface INT [Figs. 4(d) and 5(c)], perhaps due to the ability of crowders to more closely approach and descreen interacting charges in the bound state when the interface is flat; with a lower solvent dielectric that cannot account for finite crowder size, one may expect this effect to be more accentuated, as is observed. Figures 4(e) and 5(d) show that while not identical, there are no consistent, significant differences across all charge distributions in the total overall binding free energies between the two interfacial geometries. Nevertheless, these data show the manner in which crowding can change the subtle interplay between desolvation and interaction depending on the geometry at the binding interface and the charge distributions of the partners. Importantly, an implicit crowding model makes qualitatively different predictions in certain cases, as it cannot account for the finite-size effects of actual crowders.

As we have highlighted certain physical phenomena via unweighted averages and shown also how physically relevant sampling alters their relative contributions, we exclusively present the more realistic Boltzmann-weighted averages for the remainder of the properties considered, with a subset of unweighted results included as supplementary material.

#### 3. Crowding effects on electrostatic binding free energetics appear to depend modestly at best on relative sizes of the binding partners

To assess electrostatic crowding effects for partners of different sizes, the relative sizes of partners were systematically changed, keeping local interfacial geometries and the total complex volume the same (to within <1%). Figure 6(a) shows the three systems considered, ranging from one in which the partners were of the same size to one in which the so-called ligand was much smaller than the receptor. The Boltzmann-weighted changes in LDP, RDP, INT, and overall electrostatic binding free energies are shown as a function of their respective uncrowded values for explicitly modeled crowders [Figs. 6(b)–6(e); for unweighted averages see Fig. S4 in the supplementary material] as well as the implicit model using a solvent dielectric constant of 57 (Fig. S5 in supplementary material). With explicit crowders, there are no robust differences in crowding effects between partners of differing relative size. The implicit model predicts a greater strengthening of interactions for the most asymmetrically sized partners (Fig. S5). This differing result may arise because the most asymmetric complex has its interfacial charges most solvent exposed such that lowering the solvent dielectric would more drastically descreen bound state interactions than finite-sized crowders that may not be able approach the interface as robustly.

### B. Electrostatic effects of crowding due to physical properties of the crowders

In this section, we consider how the crowder charge, shape, and size may affect electrostatic binding interactions within our model system.

#### 1. Crowder charge can affect binding affinity

To assess the effect of crowder charges on binding free energies, four scenarios were studied, shown schematically in Fig. 7(a). In each case, the binding partner geometries and the set of 200 crowder placements are the same as those used to generate the results shown in Fig. 2. Monopoles on each partner were generated by charging the four interfacial atoms on each partner as well as certain surface and buried atoms, as described in Sec. II. Each atom was charged to +0.1e to generate an overall monopole of +4.3e in the case where partners both had positive monopoles. For opposite monopoles, each charged atom on the positive partner was charged to +0.2 to create asymmetry between the two crowded systems. In this case, the total monopoles on the two partners were −4.3e and +8.6e. Each crowder was charged to either +1e or −1e, as described in Sec. II.

Figure 7(b) shows the Boltzmann-weighted average ΔΔ*G*° due to the charge distributions on the crowders (i.e., the binding free energy difference when the crowders were charged vs when they were uncharged, assuming identical placements in each case). Here, the limiting case of Boltzmann weighting dominated by electrostatics is more appropriate. Upon charging crowders, the measured average binding free energy was lowered in three cases, with the result significant to *p* < 0.05 in the “neg++” and “pos++” cases. The “neg−++” scenario had a higher measured average binding free energy upon charging the crowders, although with our current sampling, it was not statistically significant to *p* < 0.05. The Boltzmann-weighted average neg−++ and pos++ values were least robust to bootstrapping methods noted in Sec. II, and therefore, more sampling will be needed in future work to provide quantitative predictions on specific systems. Nevertheless, results ranged from suggesting no significant difference to a significant weakening in the neg−++ case and from suggesting a strengthening that was not statistically significant to the one that was significant in the pos++ case (data not shown), and our following analyses will assume these more conservative ranges.

Charged crowders can affect binding free energies in at least two ways. First, interactions between each partner and the crowders are strengthened upon binding relative to the unbound state due to less solvent available to screen each bound state partner. This effect would be expected to make binding stronger in the case of negatively charged crowders and weaker with positively charged crowders in each of our cases because the complex bears an overall positive monopole. Second, binding can be affected due to differing interactions between unbound and bound state crowders due to differences in crowder placement. In each unbound state, crowders can occupy locations ultimately occupied by the binding partner. Here, because of the overall positive monopole of the complex in all cases, positively charged crowders may be expected to destabilize the unbound state relative to the bound state when compared to identical uncharged crowders and therefore favor binding more than negatively charged crowders, which would interact favorably in the unbound state. Therefore, the two effects should counteract one another.

We can separate out these two effects to see which dominates by calculating the change in binding free energy due to charging the crowders when the crowder placement is kept the same in the bound and unbound states. In such a case, only the first effect noted above would be seen. The results of this case are shown in Fig. S6. Negatively charged crowders strengthen binding, as expected. While not statistically significant at our current sampling, our measured average suggests destabilization due to positively charged crowders. The reason that positively charged crowders do not destabilize as much as negatively charged crowders stabilize may be because the highly unfavorable bound state configurations with positively charged crowders are weighted far less [indeed, using unweighted averages, the destabilization was statistically significant to *p* < 0.05 (data not shown)]. Nevertheless, the data overall are reasonably consistent with the expected prediction based only on the first effect.

The data in Fig. 7(b) do not qualitatively match those in Fig. S6 for multiple cases, suggesting that the second effect—the differing interactions between the partners and crowders in the bound and unbound states—is important and can dominate. The predicted lack in gain of binding affinity in the “neg−++” case may result from highly weighted stabilized unbound states in which negative crowders favorably interact with the highly positive binding partner. Likewise, the observed gain of binding affinity in the “pos++” case may result from the destabilizing crowder-unbound state repulsive interactions.

Previous work that, in part, used identical crowder placements in the bound and unbound states—and therefore accounted for only the first of the above two effects in a controlled manner—also showed a lowering of binding free energy when crowders had the opposite monopole as the bound state and a raising of binding free energy when crowders had the same monopole as the bound state.^{59} Here, we show the importance of considering both effects explicitly. Probing the causes and consequences of charged crowder configurations is important ongoing work; indeed, experimental work has shown that a redistribution of charged protein crowders, partially in response to the system charge, can create a “reaction field” that affects the electric fields of surface side chains.^{39}

#### 2. Crowder shape has minimal effects on binding free energies in our model system

To probe the effect of crowder shape on electrostatic binding free energies, the model complexes in Fig. 2(a) were surrounded by either spherical or rod-like crowders of nearly identical volumes, randomly placed to a complex and crowder total volume fraction of ∼20% [Fig. 8(a)]. Figures 8(b)–8(e) show the Boltzmann-weighted average changes in LDP, RDP, INT, and total electrostatic binding free energies due to crowding with either spherical or rod-like crowders. The results are nearly identical for both crowder shapes. Unweighted averages also yielded no statistically significant differences between the two crowder shapes (Fig. S7).

Earlier work predicted a small but statistically significant effect on peptide-DNA binding due to the crowder shape.^{59} However, that model assumed identical crowder placements in bound and unbound states and highly charged binding partners with charged atoms throughout, not only at the interface. To understand which model difference accounted for the difference in predictions, we considered two additional systems here: (a) one in which the binding partners were more robustly charged throughout, with monopoles of −4.3e and +8.6e (as described in the charged crowders subsection) and (b) one using the charged binding partners in (a) but also using the same crowder placements in the bound and unbound states of a given trial. In case (a), again no significant differences were seen between the spherical and rod-shaped crowders using unweighted averaging over the trials, as was done for the placements in the previous study (Fig. S8A, *p* > 0.05). However, in case (b), which most closely resembles the previous study, rod-shaped crowders favored binding more than the spherical ones (Fig. S8B, *p* < 0.05). This suggests that the predicted effect of crowder shape can depend on whether differing crowder placements between bound and unbound states are explicitly modeled.

#### 3. Smaller crowders have greater impact on electrostatic binding free energetics

Finally, to study the effect of crowder size on electrostatic binding energetics, the binding partners in Fig. 4(a), left complex (“flat” binding interface), were crowded, in turn, to a volume fraction of ∼20% with crowders of a fixed radius of either 4, 8, or 16 Å. Boltzmann-weighted average crowding effects on desolvation, interaction, and overall binding free energy were greatest with the 4 Å crowders and smallest with the 16 Å crowders (Fig. 9; see Fig. S9 for unweighted averages) This result generalizes a previous one for the barnase–barstar system^{29} and is likely due to the ability of smaller spheres to more easily approach the partners and effectively “predesolvate” unbound states and descreen bound state interactions. Although Boltzmann-weighted predesolvation effects tend to be modest due to small weights for higher-energy predesolvated unbound states, the combined effects of predesolvation in the 4 Å case approach 1 kcal/mol for certain systems [Figs. 9(a) and 9(b)].

Modeling crowding implicitly using a lower solvent dielectric constant again predicted lowered desolvation penalties and strengthened interactions (Fig. S9, pink), but such a model cannot capture effects dependent on the crowder size.

## IV. DISCUSSION

Our goal was to provide a comprehensive, generalizable framework to better conceptualize how the effects of crowding on electrostatic binding energetics depend on the physical properties (charge distributions, shape, and size) of the interacting molecules and the crowders. We predicted that crowding not only can impact electrostatic binding energetics but also differentially does so as a function of certain system physical properties, which can have implications for binding specificity, or the relative binding energetics between a molecule and its potential partners. We also showed that in certain cases, modeling crowding implicitly via a lowered solvent dielectric constant cannot account for differences observed with a more explicit model nor can it capture consequences arising from the differing likelihoods of possible crowder placements. Figure S10 shows a summary diagram depicting the qualitative predicted effects of uncharged crowders on desolvation penalties and interactions using explicit or implicit crowding models.

As we focused on the electrostatic component of binding here, it is useful to view our observations in the context of other studies that looked at coupling between system properties and overall crowding effects. For example, through experiments on two differently shaped homodimers, each comprised of a particular variant of the GB1 protein, it was shown that effects of crowding can depend on the shapes of the binding partners, with the more compact shape stabilized to a greater extent by “inert crowders;”^{44} this result was supported by theory that focused on “excluded volume” repulsions. By combining models accounting for both excluded volume and electrostatic effects of solvent depletion, a greater set of experimental observations can be holistically addressed.

The effects of crowder charge on binding energetics have also been experimentally studied using the GB1 dimer variant in the above study as well as a second study.^{45} In the latter study, protein crowders bearing the same monopole as each GB1 monomer in the homodimer were found to favor dimerization, while crowders bearing the opposite monopole as the monomers were found to destabilize the dimer under pH conditions that further amplified their opposite monopoles. The “pos++” and “neg++” scenarios shown in Fig. 7(a) can serve as a very general model for these experiments. To make a more direct comparison to this experiment, Fig. S11 shows the difference in electrostatic binding free energies between crowded and uncrowded states for the scenarios shown in Fig. 7(a) [as opposed to showing just the difference due to adding charges on the crowders, shown in Fig. 7(b)]. The standard errors are too large to predict a difference between the two scenarios with the current sampling. Nevertheless, our data are consistent with the observation that the crowders with the same monopole as the binding partners can have a more stabilizing effect on binding than those with an opposite monopole. Of course, a more specific model is required to quantitatively treat this particular experimental system, and no experimental system can perfectly control for crowder monopole without also affecting other characteristics; this analysis shows the value of pursuing both computational studies that enable us to conceptually tease apart different contributions to binding and experiments that demonstrate what happens when such factors must couple.

The role of crowder shape in biological systems has been studied previously. For example, a model incorporating a rod-like—and not spherical—crowder shape can, in part, explain the experimentally observed elongated conformation of the urea-unfolded protein, apoazurin, in the presence of dextran-20 crowders.^{71} Here, we did not incorporate conformational flexibility into our idealized models, and such studies would be interesting future work. In terms of binding free energies, previous computational work predicts that increasingly elongated crowders yield greater increases in binding affinity in the barnase/barstar system.^{72} Here, we did not see clear effects due to the crowder shape, suggesting that volume exclusion—and not electrostatic effects due to solvent depletion—may be a dominant mechanism by which the crowder shape can influence binding.

Previous work has also supported the idea that the crowder size can affect biomolecular processes, with smaller crowders potentially having greater effects, in agreement with conclusions in the current work. However, such trends were observed via hard-sphere excluded volume effects^{48} or altered hydration shell dynamics.^{23} Understanding the interplay between these effects and electrostatic consequences of solvent depletion is important future work. Additionally, the specific heterogeneous size distribution of crowders affects molecular recognition.^{13}

To comprehensively compute binding free energies, we relied on an implicit solvent model and focused on water depletion effects. As crowding is known to alter water mobility, however, it is useful to complement this work with explicit solvent simulations. Alternatively, one can implicitly account for altered water mobility through modeling the remaining solvent surrounding the explicit crowders via a lowered solvent dielectric constant, as constrained water molecules in crowded hydration shells will exhibit a reduced polarization response.^{27,73} Preliminary calculations on the same system shown in Fig. 2 but with a lowered solvent dielectric constant of 60 surrounding explicit crowders and binding partners show a greater lowering of desolvation penalties and strengthening of interactions (Fig. S12). This result suggests that both water depletion and lowered water mobility can combine to have more robust effects in crowded environments.

Our focused study on the interplay between various electrostatic impacts of crowding due to water depletion should be interpreted in conjunction with studies focused on other consequences of crowding, such as excluded volume effects, nonpolar enthalpic effects, or solvent entropy changes via the hydrophobic effect. Additionally, we sampled crowder placements either randomly or via Boltzmann weighting based solely on electrostatics, neither of which may reflect the “true” distribution of sampled states. Sampling based on MD simulations or algorithms using a holistic energy evaluation is important for more quantitative predictions.

Our exploration of a model system generated several experimentally testable hypotheses. Building off work noted above,^{45} one can continue to experimentally assess the effects of crowders bearing different monopoles on binding free energetics between partners with particular charge distributions. To help isolate electrostatic effects of crowding as a function of either the shape or size of the partners or the crowders, one can look at how ionic strength effects on binding change upon crowding in systems differing in these properties. The results from this theoretical study create a generalized framework that can be sampled, explored, and further developed through experiment to create an ongoing cycle and productive dialogue about environmental effects on molecular recognition.

## SUPPLEMENTARY MATERIAL

Supporting tables (SI and SII) and supporting figures (Figs. S1–S12) are found in the supplementary material.

## ACKNOWLEDGMENTS

This work was supported by the National Science Foundation (Grant No. NSF MCB1615313). M.L.R. was the Whitehead Chair in Critical Thought at Wellesley College during initial phases of this project. The authors thank Artemis Metaxa-Kakavouli, Helena Qi, and Ying-Yi Zhang for helpful software, Jaydeep Bardhan and Bruce Tidor for helpful discussions, and Cassandra Pattanayak for helpful discussions regarding statistical analyses.

## DATA AVAILABILITY

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

## AUTHORS’ CONTRIBUTIONS

M.L.R. oversaw the study. R.K. and M.L.R. jointly planned the study, carried out all analyses, and wrote the manuscript.

## DEDICATION

This paper is dedicated to Professor Nancy H. Kolodny who for decades catalyzed women’s excitement about physical chemistry at Wellesley College and has served as a role model for finding joy and purpose in life through both one’s career and family.