Quantifying the Influence of the Crowded Cytoplasm on Small

Jun 21, 2016 - We find that in typical, moderately crowded cell cytoplasm (ϕ ≈ 0.8), ... delineate the extent that cytoplasmic crowders influence s...
1 downloads 0 Views 3MB Size
Download PDF
Article pubs.acs.org/JPCB

Quantifying the Influence of the Crowded Cytoplasm on Small Molecule Diffusion Peter M. Kekenes-Huskey,* Caitlin E. Scott, and Selcuk Atalay* Department of Chemistry, University of Kentucky, Lexington, Kentucky 40506, United States S Supporting Information *

ABSTRACT: Cytosolic crowding can influence the thermodynamics and kinetics of in vivo chemical reactions. Most significantly, proteins and nucleic acid crowders reduce the accessible volume fraction, ϕ, available to a diffusing substrate, thereby reducing its effective diffusion rate, Deff, relative to its rate in bulk solution. However, Deff can be further hindered or even enhanced, when long-range crowder/diffuser interactions are significant. To probe these effects, we numerically estimated Deff values for small, charged molecules in representative, cytosolic protein lattices up to 0.1 × 0.1 × 0.1 μm3 in volume via the homogenized Smoluchowski electro-diffusion equation. We further validated our predictions against Deff estimates from ϕ-dependent analytical relationships, such as the Maxwell−Garnett (MG) bound, as well as explicit solutions of the time-dependent electrodiffusion equation. We find that in typical, moderately crowded cell cytoplasm (ϕ ≈ 0.8), Deff is primarily determined by ϕ; in other words, diverse protein shapes and heterogeneous distributions only modestly impact Deff. However, electrostatic interactions between diffusers and crowders, particularly at low electrolyte ionic strengths, can substantially modulate Deff. These findings help delineate the extent that cytoplasmic crowders influence small molecule diffusion, which ultimately may shape the efficiency and timing of intracellular signaling pathways. More generally, the quantitative agreement between computationally expensive solutions of the time-dependent electro-diffusion equation and its comparatively cheaper homogenized form suggest that the latter is a broadly effective model for diffusion in wide-ranging, crowded biological media.

■

INTRODUCTION Intracellular biochemical reactions commonly rely on the diffusion of charged, small molecules1 between regions where substrates are stored or synthesized to where they are ultimately utilized. These intracellular reactions can be strongly influenced by effective diffusion rates of their substrates,2 which may be depressed by up to an order of magnitude relative to unrestricted diffusion in bulk solution.3,4 In part, protein-, carbohydrate-, and nucleic acid−based “crowders” restrict the intracellular volume accessible to diffusing substrates and thereby reduce their rates of transport (reviewed here5). However, additional factors beyond volume exclusion can influence the mobility of molecular diffusers. These factors include the distribution and activity of proteins or charged surfaces of organelles that selectively bind substrates, as well as the substrate’s shape and charge.3,6−8 Improved quantitative predictions that delineate the relative contributions of these factors on shaping biomolecule diffusion is a challenging, but necessary, endeavor for understanding biochemical reactions and signaling in vitro and in vivo.9−11 A wide range of experimental and simulation approaches have provided considerable insight into the effective diffusion rates of small molecules involved in intracellular signaling (see reviews3,12). Methods such as fluorescence or raster image correlation spectroscopy,13,14 and nuclear magnetic resonance15−17 have yielded refined estimates for diffusion rates of small molecules such as calcium, magnesium, and AMP in © XXXX American Chemical Society

crowded intracellular media. A common theme emerging from these studies is that small molecule diffusion rates are substantially smaller than measurements in crowder-free solutions and are frequently anisotropic, owing to the organization of intracellular structures, such as actin filament lattices. In part, these reductions in diffusion rates can be attributed to a reduced free volume fraction,3 although the nature of substrate/crowder interactions can substantially strongly modulate the effective rate.6 However, given the considerable variations in the sizes, molecular composition, and electrostatic properties of cytoplasmic crowders, isolating the contributions of each factor via experimental means is challenging.9 In this regard, computational modeling is a strong complement to experimental techniques for describing the molecular underpinnings of experimentally observed diffusion rates and their influence on intracellular signaling.9,18 A variety of computational approaches for modeling diffusion in crowded domains has been reported in the literature.18 The most common of which are based on explicit, particle-based representations of the diffuser or both diffuser and crowders, lattice-based models that assume discrete, finite-length hops for Special Issue: J. Andrew McCammon Festschrift Received: April 16, 2016 Revised: June 18, 2016

A

DOI: 10.1021/acs.jpcb.6b03887 J. Phys. Chem. B XXXX, XXX, XXX−XXX

Article

The Journal of Physical Chemistry B

Figure 1. Configurations of domains used in this study. (a) Diffusion of charged particle (red spheres) in the absence of crowders (black spheres). Steady-state concentration gradient is indicated by higher concentration of diffusers (red background) at origin x = 0 and lower concentrations (blue background) at x = L. (b) Same as panel a, but with neutrally charged diffuser (black spheres) (c) Same as panel b, except crowders are assigned repulsive electrostatic potentials.

microscopic scale, which ultimately yields an effective diffusion tensor that is appropriate for transport on the macroscopic scale. We and others have shown that homogenized PDE models applied to ionic and nucleotide diffusion in structurally detailed, intracellular environments are in excellent agreement with experimentally measured trends in hindered diffusion rates and anisotropy.8,42−48 To our knowledge, applications of homogenization theory to biomolecular systems has been limited to a single or a small number of crowder types, thus the challenge remains to understand transport in highly heterogeneous systems with variable protein sizes, positions, and charge distributions.5 In this study, we therefore consider the diffusion of a small, anionic molecule (such as adenosine monophosphate (AMP)) in an immobile, submicrometer scale lattice of common cytoplasmic proteins (Figure 1), for which we precisely control the crowder distribution and charge density, as well as the electrolyte concentration. Specifically, we used realistic threedimensional geometries of the cytoplasm, as well as twodimensional (2D) and three-dimensional (3D) geometries comprising charged, spherical inclusions to unravel molecular determinants of small molecule diffusion rates. Our modeling results for the realistic and 2D/3D simplified geometries are consistent with prior studies,43,49,50 in that decreasing accessible volume progressively hinders diffusion, while long-range interactions can further slow or accelerate diffusion (repulsive versus attractive, respectively). Given that analogous trends are exhibited in both the 3D and 2D geometries, the majority of our simulations are done using the latter geometries, namely for rigorously examining the effects of ionic strength, surface charge potential and heterogeneous crowder shapes and distributions. Overall, our simulations offer quantitative insight into the effective diffusion rates of small biomolecules in heterogeneous cytoplasmic fractions, thereby providing important guidelines for electrokinetic transport in confined51−53 or crowded media.

the diffuser, or spatially continuous differential descriptions of the diffuser and its diffusion domain. These systems are prominently modeled as either deterministic or stochastic processes,18 or combinations thereof.19 Among these include particle-based approaches such as molecular dynamics20−22 and Brownian dynamics6,7,23−25 that explicitly represent the geometry of the diffuser and its crowders. Their computational complexity, however, generally limit their applications to nanometer spatial and microsecond temporal scales that are insufficient for describing subcellular transport. On the other hand, particle-based methods that represent diffusers as points, including Smoldyn26 and MCell27 scale well to longer time scales and to greater spatial extents appropriate for subcellular phenomena. Point representations, however, clearly sacrifice details of crowder shape and crowder/diffuser interactions that are known to influence transport and reaction kinetics.18,24,28−31 As a compromise between these two extremes, partial differential equation (PDE)-based continuous diffusion models enable the preservation of atomistic details of crowder shape and composition, as well as nonbound interactions (like electrostatic forces32) with an implicitly represented diffuser.33−35 Such PDE approaches are most appropriate for a continuum of diffusing species, such as ions and nucleotides, that are fairly concentrated, significantly smaller than typical cellular crowders (100 kDa19), and do not substantially influence the local electrostatic potential of the medium. In earlier works, we have applied such continuum models to quantify the influence of molecular crowding around enzymes and electrostatic interactions on biochemical reaction rates and efficiency.29,30,36,37 Nevertheless, the computational expense grows exponentially with increasing number of explicitly represented crowders, which renders these approaches impractical for the nanometer length scales relevant to subcellular transport phenomena. To reduce this computational burden, multiscale theories such as “homogenization” can facilitate the extrapolation of molecular-scale information onto transport processes occurring over submicrometer and longer length scales.38−41 In essence, homogenization theory posits that a continuum transport process occurs on two decoupled spatial scales: a microscopic scale within which the system is in steady-state, and a multifold larger macroscopic scale. This leads to a modified PDE (described in the Methods section) that is solved on the

■

METHODS Structural Models of the Crowded Cytoplasm. The model cytoplasm considered in this study was based on Brownian dynamics (BD) trajectory snapshots of common cytosolic proteins and nucleic acids in E. coli by McGuffee et al.24 (see Figure 2). We assume lattice proteins are immobile relative to the rapidly diffusing small molecule, given the large B

DOI: 10.1021/acs.jpcb.6b03887 J. Phys. Chem. B XXXX, XXX, XXX−XXX

Article

The Journal of Physical Chemistry B

additionally randomize the protein size and position in accordance with Figure S2, based on a Monte Carlo protocol (see Supporting Information, the Monte Carlo Protocol for Generating Crowder Ensembles section). We additionally assumed crowders have boundary potentials of ψ0 = ±19.2 [mV], based on ζ potential measurements ranging from 0.5 to 0.75 kBT for bovine serum albumin in 150 mM potassium chloride (KCl) and approximately neutral pH.60 These parameters are summarized in Table S1. Diffusion Models. Our numerical approach is based on a time-dependent diffusion model governing the diffuser concentration c(x, t) on a domain Ωm, ∂c = −∇·D∇c ∂t

Figure 2. (a) 808 Å × 808 Å × 808 Å snapshot of crowded bacterial cytoplasm simulation from McGuffee et al.24 (b) Prototypical representative volume element (RVE) used for homogenization represented in color.

on

Ωm

(1)

where D is the small molecule diffusion constant. To reflect the influence of an electrostatic potential of mean force acting on the diffuser, we formulate the concentration function as a Smoluchowski equation such that

difference in diffusion constants (e.g., diffusion constant for AMP is at least 10-fold greater than the median for E. coli cytoplasmic proteins54). Our numerical approach (see Supporting Information) uses three-dimensional finite element meshes based on these BD data, which were constructed by (1) exporting their solvent-accessible molecular surfaces from VMD,55 (2) meshlab56 refinement to remove defects, (3) tetrahedralization via tetgen,57 (4) conversion to.xml format for compatibility with the PDE solver FEniCS58 via DOLFINCONVERT. Details of this procedure are available in the Supporting Information (Volumetric Mesh Generation for Cytoplasm Model section). The final result was a finite element mesh of the representative volume element (RVE) representing the aqueous phase between the crowders, of which an example is shown in Figure 4a. We additionally created lattices of cylindrical (2D) or spherical (3D) primitives via GMSH57 to afford more control over testing the influence of crowder shape, size, and surface charge on diffusion rates. The size distribution were chosen to best approximate those of the original cytoplasm model, based on our characterization of the crowder radius of gyration (Rg) (see Figure S2), for which the average Rg was approximately 30 Å. The appropriate accessible volume fraction (accessible to the diffusing species), ϕ, was determined based on simulation data from McGuffee et al.24 that used 1008 proteins within an 808 × 808 × 808 Å3 box, or ϕ = 0.78. We note that our definition of ϕ as the accessible volume fraction is consistent with59 and our earlier publications,8,45 although in other contexts the volume fraction is based on the inaccessible fraction (e.g., 1 − ϕ). On the basis of a perfect simple cubic arrangement of Rg = 30 Å crowders with a center of mass distance, dcom, of 80.6 Å would yield ϕ = 0.78, which is consistent with crowding estimates in E. coli.4 Given the computational expense in evaluating large numbers of three-dimensional finite element meshes, we additionally defined two-dimensional reference systems using dcom = 80.6 Å (based on surface distances for spherical crowders) and dcom = 113.0 Å for volume fractions of ϕ = 0.57 and ϕ = 0.78, respectively. We consider both volume fractions in order to examine the potential influences of intercrowder distances and crowder volume fractions analogous to the 3D system. The intercrowder distance, H, is related to the distance between centers of mass (e.g., dCOM = H + 2Rg); while dcom may be a more familiar metric, H is better suited for assessing the range of electrostatic interactions between charged protein surfaces. Where indicated in the text, we

∂c = −∇·De−βqψ ∇(e βqψ c) ∂t = −∇·D̃ ∇c ̃

on Ωm

(2)

on Ωm

(3)

where q is the diffuser charge and ψ is the electrostatic potential. β ≔ 1/kBT where kB is the Boltzmann constant, and T is the temperature. The second line reflects the Slotboom form using D̃ (y) ≔ D(y) e−βqψ(y) and c̃(y) ≔ eβqψ(y)c(y). The evaluation of the electrostatic potential arising from charged crowders is described in the Supporting Information, Poisson− Boltzmann Model for the Electrostatic Potential. Direct numerical simulation of eq 2 can yield effective diffusion parameters, as we describe in Sec. Estimating effective diffusion rate from continuum diffusion simulations; in practice, however, this entails evaluating the PDE over a considerable number of time intervals, which can be prohibitively expensive. Therefore, we propose solving a ‘homogenized’ steady-state problem that directly yields effective diffusion constants, albeit at a computational cost that can be at least an order of magnitude smaller. The two-scale homogenization approach39,61 applied to eq 2 assumes that Ωm has two length scales of interest, x and y ≔ x/ϵ where the constant ϵ is assumed to be small compared to the dimensions of Ωm, and that Ωm is periodic in y where y corresponds to a nanometer length-scale, whereas x corresponds to submicrometer or longer. The y-scale periodic unit, otherwise known as the RVE, is denoted by Ω, within which the accessible region is defined as Ωϵ. This relationship is illustrated in Figure 2b, where the RVE Ω (colored) is a subset of the macroscopic domain Ωm (gray). We further assume that D is y-periodic, such that D(x, y) can be written as Dϵ(y). Before continuing, we briefly discuss the implications of this assumption. Namely, assuming “microscopic” periodicity allows us to represent the influence of interacting crowders on diffusion implicitly through an effective diffusion coefficient, Deff. In practice, time-dependent and spatially dependent simulations of intracellular signaling in cell-shaped geometries could then proceed using the Deff in eq 2, for which the Deff term effectively “coarse-grains” the impact of crowders on subcellular diffusion. Two-scale homogenization proceeds by rewriting eq 2 using the expansion of c in powers of ϵ (i.e., c = ∑i ciϵi) and the revised gradient operator (i.e., ∇ = ∂x + ϵ−1∂y). In previous C

DOI: 10.1021/acs.jpcb.6b03887 J. Phys. Chem. B XXXX, XXX, XXX−XXX

Article

The Journal of Physical Chemistry B

Figure 3. Two-dimensional (a) and three-dimensional (b) simulations domains used for validating continuum and homogenized Smoluchowski equation (HSE) models of hindered diffusion. (c) Comparison of 2D (red) and 3D (black) effective diffusion coefficient for HSE (solid symbol) and continuum model (open symbol) as a function of accessible volume fraction with appropriate analytical bounds for noninteracting crowders (Maxwell-Garnett (MG) for cylinders and Hashin-Shtrikman (HS) for spheres) are shown with dashed and solid lines, respectively.

works,8,45 we demonstrated that the resulting equation yields the homogeneous y-scale steady state problem: find χ such that ⎛ ⎛ ⎞ ∂χ ⎞ ∇·⎜Dijϵ⎜⎜δjk + k ⎟⎟⎟ = 0, ⎜ ∂yj ⎠⎟⎠ ⎝ ⎝

on Ωϵ

⎛ ⎛ ⎞⎞ ⎜D ϵ⎜δ + ∂χk ⎟⎟ ·n ̂ = 0, ⎜ ij ⎜ jk ∂yj ⎟⎠⎟⎠ ⎝ ⎝

on Γϵ :=∂Ωϵ

⎛ ⎛ ⎞⎞ ⎜D ϵ⎜ ∂χk ⎟⎟ ·n ̂ = 0, ⎜ ij ⎜ ∂y ⎟⎟ ⎝ ⎝ j ⎠⎠

dimensional crowder domains in Figure 3a,b by (1) approximating the effective diffusion constant using semianalytical bounds and (2) explicitly solving the time-dependent diffusion eq (eq 1) in the presence of crowders to obtain an effective diffusion constant. Figure 3c depicts the normalized effective diffusion coefficients of a neutral diffuser predicted via HSE (solid symbols) for two- and three-dimensional crowder domains with accessible volume fractions ranging from ϕ = 0.5−1.0. For ϕ →1.0, Deff approaches the (normalized) bulk rate at small occluded volume fractions. As the free volume fraction decreases toward ϕ = 0.5 with increasing crowder density, Deff monotonically decreases to roughly 30% and 40% of its bulk value for spherical and cylindrical inclusions, respectively. We further note that Deff values predicted for cylindrical inclusions tend to be slightly smaller (up to approximately 25%) than those for spherical crowders; therefore our simulations using two-dimensional geometries will very modestly overestimate the effects of crowders on small molecular diffusion in analogous three-dimensional domains. Nevertheless, near typical intracellular volume fractions (ϕ ≈ 0.78), Deff values for both spherical and cylindrical inclusions are comparable at roughly approximately 60−70% of the bulk diffusion rates. These HSE predictions are consistent with estimates based on two semianalytical bounds for uniformly distributed inclusions, the Hashin−Shtrikman (HS) upper bound63 for spherical crowders (solid line),

on Γ\Γϵ := ∂Ω\∂Ωϵ (4)

where Γϵ corresponds to a molecular boundary and Γ\Γϵ is typically the RVE boundary, while Dϵ is the small molecule diffusion rate within the RVE. Equation 4 is similar to eq 2, but includes the Kronecker delta, δjk, which essentially couples the microscopic and macroscopic coordinate systems. Given solutions for χ, the effective diffusion tensor (applicable to eq 2 on the macroscopic domain) is determined via Dij =

1 |Ω|

⎛

∫Ω Dijϵ(y)⎜⎜δjk + ϵ

⎝

∂χk ⎞ ⎟ dy ∂yj ⎟⎠

(5)

where Ω is the RVE volume. In this study, we report normalized self-diffusion constants along the x-direction; for D example, Deff≡ Dϵ xx(y) , except where otherwise noted. The reader

DHS =

xx

2ϕ (3 − ϕ)

(6)

and the MG bound for cylindrical particles (dashed line),64

may consult the references8,42,61,62 for additional details of the derivation. The PDEs defined in this study were numerically evaluated via the finite element method using FEniCS,58 for which a piece-wise linear basis and the default direct linear solver were used. Details of the numerical procedure are identical to protocols outlined in refs 8 and 45. All code written in support of this publication is publicly available at https:// bitbucket.org/pkhlab/pkh-lab-analyses. Simulation input files and generated data are available upon request.

DMG =

ϕ (2 − ϕ)

(7)

where ϕ is the accessible volume fraction. Lastly, we find comparable Deff predictions based on explicit time-dependent simulations of the continuum diffusion equation (open symbols, see Supporting Information section Deff from Continuum Diffusion Simulations for more information and error analyses). Given that Deff estimates for two- and threedimensional lattices show similar dependencies on accessible volume fraction and are comparable in magnitude, subsequent analyses will primarily use 2D geometries for reasons of computational efficiency, except where otherwise indicated. Influence of Protein Shape and Distribution on Effective Diffusion Rates for Neutral Diffusers. After validating our HSE model against primitive, uniformly distributed crowders, we next determined Deff based on three-dimensional snapshots

■

RESULTS AND DISCUSSION Diffusion of a Neutral Diffuser in the Crowded Cytoplasm. Model Validation for Two- and Three-Dimensional Domains. The homogenized Smoluchowski equation (HSE) enables the prediction of effective diffusion parameters based on the crowder configurations in nano- to micro-scale domains. We first demonstrate the accuracy of our implementation for a neutral diffuser in the two- and threeD

DOI: 10.1021/acs.jpcb.6b03887 J. Phys. Chem. B XXXX, XXX, XXX−XXX

Article

The Journal of Physical Chemistry B taken from published BD simulation data.24 We note that the HSE model permits a computationally inexpensive evaluation of Deff for a complex, crowder-laden cytoplasmic domain (several minutes on 16 processors), and to our knowledge, is the first application of homogenization to a structurally detailed, angstrom-resolution cytosol fraction. In Figure 4a,

as overlapping and self-intersecting vertices. Moreover, our Monte Carlo approach permits more exhaustive sampling of crowder locations and spatial configurations than would be possible from the original Brownian dynamics trajectory data made available by McGuffee et al.24 In Figure 5, we report the

Figure 4. (a) Three-dimensional solution of homogenized diffusion equation for an uncharged diffuser (see eq 4) based on the angstromresolution cytoplasm data in Figure 2. (b) HSE-predicted effective diffusion constants along the x (red), y (yellow), and z (black) directions as a function of RVE size (100, 200, and 300 Å) and the quantitative comparison of them with the analytical Hashin− Shtrikman (HS) upper bound (blue).

Figure 5. Predicted effective diffusion rates, Deff, for a neutral diffuser in crowded, two-dimensional RVEs with volume fraction ϕ = 0.78, as a function of domain size (cell length). Geometries are derived from the crowder distributions in Figure 2.

mean Deff values as a function of RVE size for configurations with randomized positions and radii (RandPos/RandRad) as well as randomized positions and constant (Rg = 30 Å) radii (RandPos/ConstRad). We found that Deff ≈ 0.65 across RVE lengths, regardless of whether randomized or fixed radii were used. As expected for larger sample sizes, the standard error in the mean Deff decreased with increasing RVE size. These data therefore indicate that for neutral diffusers, (1) reliable Deff estimates may be obtained using relatively small RVE lengths (300 nm) provided sufficiently large sample sizes and (2) Deff estimates are apparently insensitive to crowder position and size variations at cytosolic volume fractions (ψ ≈ 0.78). In other words, fine-resolution details of cytosolic globular proteins’ three-dimensional shape are likely unnecessary for predicting their influence on diffusion rates of small molecules. Therefore, we expect that the deviations in Deff as a function of direction and cell size relative to the HS bound in Figure 4 were due to small sample sizes, as opposed to variations in shape. We clearly anticipate exceptions to this trend if the protein positions or shapes are strongly anisotropic, such as for myosin/actin filaments in contractile cells,67 especially if the proteins strongly interact with the diffuser (examined in the next section). Diffusion of a Charged Molecule in the Crowded Cytoplasm. Effects of Protein Charge on Electric Potential. The electrostatic potential arising from surface-exposed polar amino acids is known to influence protein−protein and protein−ligand association thermodynamics and kinetics.68−70 Therefore, it can be expected that the electrostatic interactions of a small, charged diffuser with densely packed, immobile proteins comprising the cytoplasm would influence their effective diffusion rate, although the contributions of (1) typical protein surface potentials and (2) protein distribution in large submicrometer cytoplasmic fractions has not been wellexplored. Before investigating these contributions, we first predict the electrostatic potential from distributions of charged 2D crowders by numerical solution of the Poisson−Boltzmann (PB) equation. Further details on the determination of the electrostatic potential and assignment of crowder surface potential can be found in the Supporting Information section Estimating Effective Diffusion Rate from Continuum Diffusion Simulations.

we present the three-dimensional solution for χ (see eq 4); generally speaking, the reddish and bluish tones represent regions where diffusion is more strongly hindered relative to the grayish regions. In Figure 4b we report the effective diffusion constants along the x, y, and z axes for RVE sizes of 100−300 Å. The diffusion tensors are nearly isotropic, with ranges from 0.40 and 0.55; such a result is not surprising, given that the cytoplasmic proteins are roughly spherical in shape and have no preferential, anisotropic arrangement. Furthermore, our predictions are consistent with estimates indicating that translational diffusion rates are within an order of magnitude of bulk diffusion rates for diffuser sizes < 500 kDa3,4 and closely match the Deff = 0.51 predicted by the HS bound. The Deff predictions from these data have no obvious dependence on the RVE size, which we attribute to small sample sizes (one case for each cell length). To better understand the variance of effective diffusivity with respect to crowder positions and sizes, we characterized the Rg (see Figure S2) of crowders reflected in the McGuffee cytoplasm data. The skewed distribution had an average Rg of approximately 30 Å, with prominent peaks at 15 and 30 Å, although some radii were as much as 80 Å in size. Given the sensitivity of the effective diffusion tensor to the free volume fraction and the difficulty in constructing atomistic-resolution three-dimensional meshes of complex topologies,65 we generated randomized, 2D crowder distributions (11 per cell length ranging from 100 to 2000 Å). Specifically, each randomized distribution was chosen such that their radii of gyration well-approximated the statistics derived from the BD data shown in Figure S2. Their randomized positions were generated using a Monte Carlo approach (see Supporting Information section Monte Carlo Protocol for Generating Crowder Ensembles) that ensured crowders were nonoverlapping and within the RVE boundaries. We note that the accessible volume fraction of the original three-dimensional geometry could be modulated by shrinking or enlarging the solvent probe radius for generation of the Connolly solventaccessible surface.66 However, there is considerable burden of manually revising commonly occurring meshing defects65 such E

DOI: 10.1021/acs.jpcb.6b03887 J. Phys. Chem. B XXXX, XXX, XXX−XXX

Article

The Journal of Physical Chemistry B

Figure 6. Electric potential in (a) two dimensions and (b) one dimension between negatively charged crowders for various background KCl concentrations (5, 20, 150, and 500 mM, corresponding to a κH of 1.06, 2.14, 5.86, and 10.69, respectively), given a crowder edge spacing (H) of 5.2 nm and ϵ = −0.75 kBT.

expected to compete with the bulk electrolyte and potentially decrease the Debye length. To better approximate the diverse electrostatic properties of proteins comprising the cytoplasm, we also present electric potentials corresponding to an electro-neutral system (nearly equal numbers of positively and negatively charged crowders) for alternating (Figure S4) and randomly distributed crowders (Figure S5). For alternating charged crowders, ψ = 0 midway between crowders (e.g., d = H/2). Despite the inverse symmetry of the electric potential across the midpoint, however, an asymmetric distribution of diffusers is expected, as c(r = 0, zieψ = 0.75kT) ≈ 0.47 and c(r = 0,zieψ = 0.75kT) ≈ 2.12 for diffusers with zi = 1 and zi = −1, respectively. For randomly distributed charged crowders (Figure S5), we observe localized electropositive (red) and electronegative (blue) regions that we expect will give rise to strongly heterogeneous diffuser distributions. Effective Diffusion Rates As a Function of Charged Crowder Distribution. To investigate the influence of charged crowders on small molecule diffusion, we evaluated the HSE subject to the PB solutions from the previous section. Similar to the neutral diffuser case in the previous section, we predicted Deff values for charged cylindrical and spherical crowders; however, since we were unaware of effective diffusion rate relationships that account for both electrostatic interaction strength and volume fractions, we validated our HSE predictions against the time-dependent Smoluchowski model (eq 2). Namely, in Figure 7 we compute the normalized effective diffusion coefficient for charged diffusers (z = −1, 0, +1) using a lattice of cylindrical crowders with varying volume fractions (ϕ = 0.5−1.0). For these cases, we assumed that the crowders have uniform, negatively charged surface potentials of −19.2 mV and are immersed in a background 150 mM KCl solution (λD = 7.8 Å). Similar to our results assuming a neutral diffuser, the normalized Deff decreases with reduced accessible volume fraction regardless of charge, although the rates of decrease are charge-dependent. While the effective diffusion rates of the

In Figure 6, we predict the electric potentials arising from a uniform distribution (ϕ = 0.78) of crowders with −19.2 mV surface potentials subject to ionic strengths ranging from 5 to 500 mM. The 2D surface plot electric potential distributions are shown, as well as the spatial distribution of electric potential between two neighboring particles in one dimension (see Figure 6b). The ionic strength is reported as the dimensionless quantity, κH, for which κ is the inverse Debye length, and H is the edge-to-edge average distance between proteins. The magnitude of the electrostatic potential is clearly maximal at the protein surface (ψ0 = −19.2 [mV]) and minimal midway between neighboring proteins. As would be expected for strong electrolyte solutions, the decay of the potential as a function of distance from the crowder increases with increasing ionic strength. At κH = 10.69 corresponding to a concentration of 500 mM, the electric potential rapidly decays to zero within 4.3 Å of the crowders, thus the influence of electrostatic interactions are expected to be negligible within most of the accessible diffusion volume. However, as κH decreases, significant overlap of potentials from neighboring crowders is observed throughout the entire domain. To understand the impact of this potential on the distribution of charged diffusers, we assume local equilibrium inside the RVE, which is commonly done to ensure validity of the constitutive model at the microscale despite nonequilibrium conditions at the macroscale.40 In this regime, the Boltzmann relation c(x,zieψ) = c0 exp[−βzieψ] may be used to describe the diffuser concentration, c, in the diffuse layer (DL) relative to the bulk solution (c0). Hence, electronegative DLs are predicted to exclude negatively charged diffusers, which is analogous to increasing the crowders’ effective radii; in complement, attractive forces would localize positively charged diffusers toward the charger, thereby reducing the crowders’ effective size. This interpretation holds for the weakly attractive electrostatic interaction energies assumed in this study (less than −1 to −2 kBT) and low diffuser concentrations relative to the bulk electrolyte; beyond this regime, the diffuser would be F

DOI: 10.1021/acs.jpcb.6b03887 J. Phys. Chem. B XXXX, XXX, XXX−XXX

Article

The Journal of Physical Chemistry B

these calculations a physiologically reasonable volume fraction (ϕ = 0.78) and a bulk electrolyte with λD = 7.8 Å. The predicted Deff values from these three approaches are in quantitative agreement, with the most hindered diffusion rates (0.54) reported at ψ = −25.6 mV, corresponding to an electrostatic interaction energy of ϵ ≡ qψ = 1kBT, given z = −1 for q ≡ ze; on the other hand, the fastest rates (0.81) were found for ϵ → −1kBT. In the absence of electrostatic interactions, ϵ = 0kBT, and we recovered the MG bound (Deff = 0.64). It is interesting to note that log(Deff) is nearly linear with respect to the electrostatic interaction energy as shown in D(ψ ) Figure 8b forlog D ′ = βγ ϵ. We found the greatest agree-

Figure 7. Predicted normalized diffusion coefficients based on negatively charged crowders with volume fraction ϕ, similar to the configuration in Figure 6. Coefficients are computed using the HSE (open symbols), continuum (closed symbol-dashed line), and Maxwell-Garnett (MG) bound (cross-symbol). Diffuser charges of z = −1 (circle), z = 0 (square), and z = +1 (triangle) are given, assuming CKCl = 150 mM and ψ0 = −19.2 mV (0.75 kBT for q = 1) for the crowders. The black dotted-line corresponds to the center of mass distance between crowders. We present an analogous plot for spherical inclusions in Figure S6.

ψ =0

ment at low ionic strengths (see dashed-dot line for CKCl = 1 mM), for which the spatial decay in the electrostatic potential energy, ϵ(r), and diffuser probability density, CAMP(r), were mostly linear (see Figure S7). Given this correspondence for modest electrostatic interaction energies, these data suggest a simple correction to the MG bound of the form.

neutral diffusers (z = 0, black squares in Figure 7) follow the MG analytical bound (red crosses), negatively charged diffusers exhibit slower diffusion (blue dots). In contrast, for attractive interactions (red triangles), the normalized diffusion coefficients exceed those of the uncharged and negatively charged diffusers. Moreover, the HSE results (hollow symbols) and the time-dependent continuum models (line-symbols) are in quantitative agreement, which validates our use of the HSE in modeling electrokinetic phenomena, albeit at a fraction of the computational expense. We attribute hindered diffusion under repulsive interactions to the diffuser’s exclusion from the crowders’ DLs, as illustrated in Figure 6. In contrast, the diffusional acceleration by weakly attractive interactions observed here and in related experiments7,23,29,30,43,71−73 has been rationalized as the “smoothing out” of the chemical potential energy surface roughened by crowder-induced excluded-volume regions.7 We additionally examined the dependence of chargedmediated Deff values on the electrostatic interaction strength. In Figure 8a, Deff values are reported from solutions based on the time-dependent diffusion equation, in which ψ is determined from the nonlinear (solid line) and linearized (dashed-line) PB equation, and based on the HSE (squares symbols) using the linearized PB equation. We assumed for

Deff =

ϕ e βγϵ 2−ϕ

(8)

where γ may be determined by a linear fit to the data in Figure 8b. This expression is unlikely to hold for strongly attractive interactions (−qψ0 ≫ (2−3)kBT), which would result in adsorption of the diffuser and a corresponding reduction in diffusion rates.7 This topic is further discussed in the Limitations section. In Figure 8a we compared the HSE model (red square) with continuum simulations (linearized and nonlinearized PB equations) for different surface potentials when λ = 7.85 Å, which corresponds to a 150 mM background electrolyte concentration and ϕ = 0.78. Strong agreement for the predictions from the HSE model and the time-dependent models are observed. We additionally present in Figure 8a effective diffusion rates that were calculated using the HSE model for electro-neutral (nearly equal numbers of positivelyand negatively charged crowders) crowder distributions as representations of a typical crowded cytoplasm. For these calculations, we assumed alternating positive/negative crowders (diamond) as well as randomly charged (circle) or randomly charged and randomly distributed crowders (triangle symbols). For alternating positive and negatively charged crowders, Deff

Figure 8. (Left) Normalized effective diffusion coefficient of the diffuser as a function of electrostatic potential energy, ϵ = qψ[kBT], where we have assumed crowder surface potentials of |ψ| < 25.6 [mV] and z = −1 for q ≡ ze. Results based on the linearized and nonlinear PB equations are compared with HSE when λ = 7.85 Å, which corresponds to 150 mM salt concentration and ϕ = 0.78. (Right) Logarithm of Deff from left figure is shown, in addition to a linear fit to Deff values predicted at 1 mM ionic strength, based on eq 8. G

DOI: 10.1021/acs.jpcb.6b03887 J. Phys. Chem. B XXXX, XXX, XXX−XXX

Article

The Journal of Physical Chemistry B

Figure 9. Normalized effective diffusion coefficient (via HSE) for (a) repulsive and (b) attractive electrostatic interactions as a function of κH and varying crowder surface distances (H). Results for H = 2.57, 5.2, 9, or 14 nm are denoted with crosses, circles, diamond, triangles, and squares, respectively. Evaluation of MG bound for each volume fraction is represented by a pink star. Solid lines represent empirical fits based on “effective” volume fractions from eq 9 (see Figure S9).

strengths as shown in Figure 9a,b. The γ parameter scales approximately linearly with surface potential independent of accessible volume fraction as illustrated in Figure S9. The revised analytical bound (line), eq 9, agrees well with the diffusion coefficient calculated using HSE for moderate to high κH when γ is 0.45 and 0.66 for repulsive and attractive interactions, respectively, and it is independent of the accessible volume fraction. The small deviation at ϕ = 0.717 when κH < 2.5 can be attributed to strong DL overlap, for which eq 8 is a more appropriate bound. Limitations. Our HSE simulations of small molecule diffusion through matrices of weakly charged cytosolic proteins offer insight into the influence of protein shape, size, charge, and density on transport within the cell. The fundamental assumption for homogenization theory is that there exists an underlying periodic structure. In biological media, this assumption can only be satisfied approximately, as the distribution of proteins is not strictly periodic. Nevertheless, it is apparent from this study and others45,75 that homogenization-based predictions of effective diffusion tensors are in excellent agreement with experimental measurements. We attribute this success to the fact that the position of each inclusion does not strongly deviate from the “typical” unit cell, or similarly, the statistical distribution within a given unit cell is consistent across all cells. In fact, this is the basis for a homogenization scheme valid for nonperiodic materials,76 wherein it is assumed that periodic spherical inclusions are remapped to nonperiodic positions that remain strictly confined within the unit cell. Additionally, for examination of systems with periodicity-breaking defects, such as the presence of intracellular organelles in a cytoplasmic protein lattice, periodic correctors for localized perturbations may be appropriate.77,78 There are a number of additional model improvements that may be considered, particularly for describing highly charged systems. For the crowder surface potentials considered in this study (|ψ| ≤ 25 mV), the linearized PB was sufficient to describe the diffuse layer about charged crowders as well as their overlapping DLs. For more highly charged diffusers or crowders (|ψ| ≫ 25. mV), a numerical solution of the nonlinear PB will be necessary to adequately model the electric potential.79 More importantly, stronger attractive electrostatic interactions will introduce several complications that did not arise in our study. As illustrated by the Boltzmann equation, ρ = ρbulk exp(−βqψ), attractive interaction energies as low as 3 kcal/mol are sufficient to increase the local diffuser

appears to modestly increase (approximately 5%) for increasing |ϵ|. This may arise due to the stronger influence of attractive interactions on diffusion relative to repulsive interactions of the same amplitude. This effect, however, is largely lost when the crowder positions or distribution of positive versus negative crowders are randomized; for these cases, the Deff coefficients are nearly uniform for interaction strengths of |ϵ| ≤ 1kBT and comparable to predictions for an uncharged diffuser. To understand the dependence of Deff on the bulk ionic strength for uniformly charged crowders illustrated in the previous section, we present HSE simulations for crowders in electrolyte solutions ranging from dilute to concentrated. For weak electrostatic interactions occurring predominantly for thin DLs (i.e., κH ≫ 10), Deff approaches the neutral diffuser limit (star symbol). As κH → 0, the electrostatic interactions between crowders and the diffusion medium are shielded to a lesser extent; as a result, the hindrance of diffusion due to repulsive crowder/diffuser interactions is intensified, particularly as the accessible volume fraction is reduced. Similarly, the enhancement of diffusion due to attractive interactions is more apparent at lower κH values. These trends are consistent with the ionic strength dependence of diffusion rates reported in experiment.71,72 In sharp contrast, for the randomized ensembles of negatively and positively charged crowders introduced in the previous section, Deff is almost entirely independent of κH, except at very low κH where minor enhancement is observed. Together, these data suggest that for randomly distributed cytoplasmic crowders, the excluded volume effect overwhelmingly determines the effective diffusion rate of small diffusing molecules; however, for regions with a disproportionate share of electropositive or electronegative crowders, as might be expected for compacted DNA in the nucleus,74 effective diffusion rates may significantly deviate from estimates based on excluded volume alone. For the latter scenario, we propose a Deff relationship based on an effective crowder radius determined by the Debye length, λ. For this approximation, we use the functional form of eq 7, but assume ϕ is dependent on λ as given by ⎛ π (R ± γλD)2 N ⎞ ⎟ ϕ = ⎜1 − L2 ⎝ ⎠

(9)

where +γ and −γ are corrections used for repulsive and attractive interactions, respectively. In contrast to eq 8, which is most accurate at low ionic strength (e.g., long Debye lengths), this approximation performs well at moderate to high ionic H

DOI: 10.1021/acs.jpcb.6b03887 J. Phys. Chem. B XXXX, XXX, XXX−XXX

Article

The Journal of Physical Chemistry B

typical intracellular accessible volume fractions and crowder interaction strengths are sufficient to reduce effective diffusion rates within an order of magnitude, which is in line with conventional estimates.3 Surprisingly, partitioning of likecharged crowders into subcellular domains can strongly perturb diffusion in localized regions of the cell. This codependence of Deff on ϕ and electrostatic interactions can be empirically modeled by slight modifications of the MG bound for neutral crowders (see eq 8 and eq 9). In our opinion, this raises the possibility that nature may exploit such perturbations to control the efficiency and kinetics of reactions in small, subcellular compartments, as is postulated for intracellular nucleotide-, Na+- and Ca2+-dependent “nanodomain” signaling.62,86,87 Moreover, given the strong agreement between Deff estimates from the time-dependent and homogenized electro-diffusion equations, the adaptability of the homogenized Smoluchowski equation to arbitrary potentials of mean force (see eq 2) could facilitate examination of how additional prominent physical phenomena, including adsorption and buffering,41 shape intracellular communication in crowded, biological systems.

concentration by 2 orders of magnitude relative to bulk. For the 1 mM AMP concentration assumed in our time-dependent modeling, the local AMP concentration near crowders could approach that of the bulk KCl and would likely significantly enhance electrostatic shielding. This contribution would mandate including both the diffuser concentration and bulk electrolyte (KCl) in the determination of the electrostatic potential; Poisson−Nernst−Planck (PNP) approaches,80 whereby the diffusion and PB equations are solved iteratively and concurrently, are routinely used to describe such phenomena. Here, protocols to homogenize the PNP would be a natural extension of our current work.48 It may also be advantageous to consider the finite sizes of electrolyte and diffuser particles, which are well-known to attenuate apparent concentrations of counterions along charged surfaces compared to theories that assume point charge electrolytes.81−83 Along these lines, significant surface interactions would likely give rise to adsorption phenomena that may require the consideration of adsorption isotherms when defining the model boundary conditions.84 Consideration of additional surface interactions, such as hydrophobic interactions and colloidal aggregation forces, may be advantageous for larger diffusing species, and in many cases could be accommodated by augmenting the qψ term of eq 2 with Lennard-Jones or DLVO potentials,79 respectively. We assumed that the diffuser was substantially smaller than the crowders, which justified our use of immobile protein lattices. Inherent in this assumption is that small molecules diffuse much faster than the crowders and thereby rapidly achieve a local steady state before significant displacement of the crowders. A recent Brownian dynamics study from Putzel et al.7 examining particle diffusion among larger crowders was also based on this assumption; however, in their Supplement, they also reported effective diffusion rates for systems composed either entirely or partially of mobile crowders that were 20% to 40% higher than the fixed system. Hence, we anticipate fairly modest increases in diffusivity under the assumption of mobile crowders. However, for larger diffusers than the small biomolecules considered here, hydrodynamic effects and subdiffusive behavior would likely demand a mobile crowder treatment.24,25,85

■

ASSOCIATED CONTENT

S Supporting Information *

The Supporting Information is available free of charge on the ACS Publications website at DOI: 10.1021/acs.jpcb.6b03887. Details on the methods including a table of simulation parameters and descriptions of Monte Carlo protocol, volumetric mesh generation, Poisson−Boltzmann model, and an effective diffusion rate estimation, results pertaining to the effective diffusion rate, as well as nine figures (PDF).

■

AUTHOR INFORMATION

Corresponding Authors

*E-mail: [email protected], (859) 257-4741. *E-mail: [email protected], (859) 218-5406. Notes

The authors declare no competing financial interest. All code written in support of this publication is publicly available at https://bitbucket.org/pkhlab/pkh-lab-analyses. Simulation input files and generated data are available upon request.

■

CONCLUSIONS We investigated the effective diffusion rates for small, charged biomolecules within crowded, submicrometer-scale cytoplasm domains as functions of crowder shape, size, distribution, electrostatic interaction energies and ionic strength. Our HSE approach recapitulated observations that neutral diffusers exhibit hindered diffusion as the accessible volume fraction is decreased due to crowding. However, our studies further indicate that in the absence of electrostatic interactions, crowder shape, size, and distribution variation have minor effects on the effective diffusion coefficient at physiological free volume fractions (ϕ ≈ 0.8). Our method also demonstrates that electrostatic interactions between crowders and diffusers substantially alter Deffs in manners highly dependent on whether interactions are attractive versus repulsive, the surface potential amplitude, and solvent ionic strength. Generally, media presenting exclusively repulsive interactions tend to hinder diffusion, whereas enhanced diffusion is observed for weakly attractive interactions. In the more likely crowded configuration consisting of randomly distributed and randomly charged crowders, our homogenization results indicate that

■

ACKNOWLEDGMENTS PKH is indebted to Andy McCammon for his mentorship in all things that diffuse and more importantly, for his foresight to apply these insights to wide-ranging problems in biophysics and physiology. His rigor and creativity remain an inspiration as PKH launches his career. This work used the Extreme Science and Engineering Discovery Environment (XSEDE), which is supported by National Science Foundation Grant No. ACI1053575.88 We are grateful to Adrian Elcock and his lab for sharing their cytoplasm simulation trajectories from ref 24 as well as helpful suggestions.

■

REFERENCES

(1) Olson, T. M.; Terzic, A. Human K-ATP channelopathies: Diseases of metabolic homeostasis. Pfluegers Arch. 2010, 460, 295− 306. (2) KIm, J. S.; Yethiraj, A. Effect of macromolecular crowding on reaction rates: A computational and theoretical study. Biophys. J. 2009, 96, 1333−1340. I

DOI: 10.1021/acs.jpcb.6b03887 J. Phys. Chem. B XXXX, XXX, XXX−XXX

Article

The Journal of Physical Chemistry B

(26) Andrews, S. S.; Addy, N. J.; Brent, R.; Arkin, A. P. Detailed simulations of cell biology with Smoldyn 2.1. PLoS Comput. Biol. 2010, 6, e1000705. (27) Kerr, R. A.; Bartol, T. M.; Kaminsky, B.; Dittrich, M.; Chang, J.C. J.; Baden, S. B.; Sejnowski, T. J.; Stiles, J. R. Fast Monte Carlo simulation methods for biological reaction-diffusion systems in solution and on surfaces. SIAM J. Sci. Comput. 2008, 30, 3126−3149. (28) Długosz, M.; Trylska, J. Diffusion in crowded biological environments: Applications of Brownian dynamics. BMC Biophys. 2011, 4, 3. (29) Eun, C.; Kekenes-Huskey, P. M.; Metzger, V. T.; McCammon, J. A. A model study of sequential enzyme reactions and electrostatic channeling. J. Chem. Phys. 2014, 140, 105101. (30) Metzger, V. T.; Eun, C.; Kekenes-Huskey, P. M.; Huber, G.; McCammon, J. A. Electrostatic channeling in P. falciparum DHFR-TS: Brownian dynamics and Smoluchowski modeling. Biophys. J. 2014, 107, 2394−2402. (31) Kekenes-Huskey, P. M.; Eun, C.; McCammon, A. Enzyme localization, crowding, and buffers collectively modulate diffusioninfluenced signal transduction: Insights from continuum diffusion modeling. J. Chem. Phys. 2015, 143, 094103. (32) Baker, N. A.; Sept, D.; Joseph, S.; Holst, M. J.; McCammon, J. A. Electrostatics of nanosystems: Application to microtubules and the ribosome. Proc. Natl. Acad. Sci. U. S. A. 2001, 98, 10037−10041. (33) Cheng, Y.; Suen, J. K.; Zhang, D.; Bond, S. D.; Zhang, Y.; Song, Y.; Baker, N. A.; Bajaj, C. L.; Holst, M. J.; McCammon, J. A. Finite element analysis of the time-dependent Smoluchowski equation for acetylcholinesterase reaction rate calculations. Biophys. J. 2007, 92, 3397−3406. (34) Song, Y.; Zhang, Y.; Shen, T.; Bajaj, C. L.; McCammon, J. A.; Baker, N. A. Finite element solution of the steady-state Smoluchowski equation for rate constant calculations. Biophys. J. 2004, 86, 2017− 2029. (35) Zhou, H.; Szabo, A. Theory and simulation of stochasticallygated diffusion-influenced reactions. J. Phys. Chem. 1996, 100, 2597− 2604. (36) Eun, C.; Kekenes-Huskey, P. M.; McCammon, J. A. Influence of neighboring reactive particles on diffusion-limited reactions. J. Chem. Phys. 2013, 139, 044117−044117. (37) Kekenes-Huskey, P. M.; Gillette, A.; Hake, J.; McCammon, J. A. Finite-element estimation of protein-ligand association rates with postencounter effects: Applications to calcium binding in troponin C and SERCA. Comput. Sci. Discovery 2012, 5, 014015. (38) Allaire, G. Introduction to Homogenization Theory; Ecole Polytechnique: 2012 (39) Auriault, J. L.; Boutin, C.; Geindreau, C. Homogenization of Coupled Phenomena in Heterogenous Media; Wiley-ISTE: London, U.K., 2010. (40) Schmuck, M.; Berg, P. Effective macroscopic equations for species transport and reactions in porous catalyst layers. J. Electrochem. Soc. 2014, 161, E3323−E3327. (41) Schmuck, M.; Pavliotis, G. A.; Kalliadasis, S. Effective macroscopic interfacial transport equations in strongly heterogeneous environments for general homogeneous free energies. Appl. Math. Lett. 2014, 35, 12−17. (42) Goel, P.; Sneyd, J.; Friedman, A. Homogenization of the cell cytoplasm: The calcium bidomain equations. Multiscale Model. Simul. 2006, 5, 1045−1062. (43) Bourbatache, K.; Millet, O.; Aït-Mokhtar, A. Ionic transfer in charged porous media. Periodic homogenization and parametric study on 2D microstructures. Int. J. Heat Mass Transfer 2012, 55, 5979− 5991. (44) Chen, K. C.; Nicholson, C. Changes in brain cell shape create residual extracellular space volume and explain tortuosity behavior during osmotic challenge. Proc. Natl. Acad. Sci. U. S. A. 2000, 97, 8306−8311. (45) Kekenes-Huskey, P. M.; Liao, T.; Gillette, A. K.; Hake, J. E.; Zhang, Y.; Michailova, A. P.; Mcculloch, A. D.; McCammon, J. A.

(3) Dix, J. A.; Verkman, A. S. Crowding effects on diffusion in solutions and cells. Annu. Rev. Biophys. 2008, 37, 247−263. (4) Ellis, R. J. Macromolecular crowding: An important but neglected aspect of the intracellular environment. Curr. Opin. Struct. Biol. 2001, 11, 114−119. (5) Zhou, H.-X.; Rivas, G.; Minton, A. P. Macromolecular crowding and confinement: Biochemical, biophysical, and potential physiological consequences. Annu. Rev. Biophys. 2008, 37, 375−397. (6) Zhang, X.; Hansing, J.; Netz, R. R.; DeRouchey, J. E. Particle transport through hydrogels is charge asymmetric. Biophys. J. 2015, 108, 530−539. (7) Putzel, G. G.; Tagliazucchi, M.; Szleifer, I. Nonmonotonic diffusion of particles among larger attractive crowding spheres. Phys. Rev. Lett. 2014, 113, 138302. (8) Kekenes-Huskey, P.; Gillette, A. K.; McCammon, J. Predicting the influence of long-range molecular interactions on macroscopicscale diffusion by homogenization of the Smoluchowski equation. J. Chem. Phys. 2014, 140, 174106. (9) Elcock, A. H. Models of macromolecular crowding effects and the need for quantitative comparisons with experiment. Curr. Opin. Struct. Biol. 2010, 20, 196−206. (10) Schoffelen, S.; van Hest, J. Multi-enzyme systems: Bringing enzymes together in vitro. Soft Matter 2012, 8, 1736−1746. (11) Chen, A. H.; Silver, P. A. Designing biological compartmentalization. Trends Cell Biol. 2012, 22, 662−670. (12) Hake, J.; Kekenes-Huskey, P. M.; Mcculloch, A. D. Computational modeling of subcellular transport and signaling. Curr. Opin. Struct. Biol. 2014, 25, 92−97. (13) Bacia, K.; Kim, S. A.; Schwille, P. Fluorescence cross-correlation spectroscopy in living cells. Nat. Methods 2006, 3, 83−89. (14) Vendelin, M.; Birkedal, R. Anisotropic diffusion of fluorescently labeled ATP in rat cardiomyocytes determined by raster image correlation spectroscopy. Am. J. Physiol., Cell Physiol. 2008, 295, C1302−C1315. (15) Kinsey, S. T.; Locke, B. R.; Penke, B.; Moerland, T. S. Diffusional anisotropy is induced by subcellular barriers in skeletal muscle. NMR Biomed. 1999, 12, 1−7. (16) Kinsey, S. T.; Locke, B. R.; Dillaman, R. M. Molecules in motion: Influences of diffusion on metabolic structure and function in skeletal muscle. J. Exp. Biol. 2011, 214, 263−274. (17) Kinsey, S. T.; Moerland, T. S. Metabolite diffusion in giant muscle fibers of the spiny lobster Panulirus argus. J. Exp. Biol. 2002, 205, 3377−3386. (18) Takahashi, K.; Arjunan, S. N. V.; Tomita, M. Space in systems biology of signaling pathways - towards intracellular molecular crowding in silico. FEBS Lett. 2005, 579, 1783−1788. (19) Ridgway, D.; Broderick, G.; Lopez-Campistrous, A.; et al. Coarse-grained molecular simulation of diffusion and reaction kinetics in a crowded virtual cytoplasm. Biophys. J. 2008, 94, 3748−3759. (20) Hummer, G. Position-dependent diffusion coefficients and free energies from Bayesian analysis of equilibrium and replica molecular dynamics simulations. New J. Phys. 2005, 7, 34. (21) Lee, H.; Venable, R. M.; MacKerell, A. D., Jr.; Pastor, R. W. Molecular dynamics studies of polyethylene oxide and polyethylene glycol: Hydrodynamic radius and shape anisotropy. Biophys. J. 2008, 95, 1590−1599. (22) Thomas, A.; Elcock, A. Direct measurement of the kinetics and thermodynamics of association of hydrophobic molecules from molecular dynamics simulations. J. Phys. Chem. Lett. 2011, 2, 19−24. (23) Jardat, M.; Hribar-Lee, B.; Dahirel, V.; Vlachy, V. Self-diffusion and activity coefficients of ions in charged disordered media. J. Chem. Phys. 2012, 137, 114507. (24) McGuffee, S. R.; Elcock, A. H. Diffusion, crowding and protein stability in a dynamic molecular model of the bacterial cytoplasm. PLoS Comput. Biol. 2010, 6, e1000694. (25) Ando, T. T.; Skolnick, J. J. Crowding and hydrodynamic interactions likely dominate in vivo macromolecular motion. Proc. Natl. Acad. Sci. U. S. A. 2010, 107, 18457−18462. J

DOI: 10.1021/acs.jpcb.6b03887 J. Phys. Chem. B XXXX, XXX, XXX−XXX

Article

The Journal of Physical Chemistry B Molecular and subcellular-scale modeling of nucleotide diffusion in the cardiac myofilament lattice. Biophys. J. 2013, 105, 2130−2140. (46) Millet, O.; Aît-Mokhtar, A.; Amiri, O. Periodic homogenization of Nernst Planck equation. Application to the measurement of chloride macroscopic diffusibility. 8me Congrs de Mcanique 2007, 1−3. (47) Muntean, A.; Chalupecky, V. Homogenization Method and Multiscale Modeling; Eindhoven University of Technology: 2012 (48) Schmuck, M.; Bazant, M. Z. Homogenization of the poissonnernst-planck equations for ion transport in charged porous media. SIAM J. Appl. Math. 2015, 75, 1369−1401. (49) Moyne, C.; Murad, M. A. A two-scale model for coupled electro-chemo-mechanical phenomena and Onsager’s reciprocity relations in expansive clays: I homogenization analysis. Transp. Porous Media 2006, 62, 333−380. (50) Bourbatache, K.; Millet, O.; At-Mokhtar, A. Multi-scale Periodic Homogenization of Ionic Transfer in Cementitious Materials; In Advances in Bifurcation and Degradation in Geomaterials, Vol. 11; Bonelli, S., Dascalu, C., Nicot, F., Eds.; Springer Netherlands: Dordrecht, 2011. (51) Daiguji, H.; Yang, P.; Majumdar, A. Ion transport in nanofluidic channels. Nano Lett. 2004, 4, 137−142. (52) Yeh, L.-H.; Zhang, M.; Qian, S. Ion transport in a pH-regulated nanopore. Anal. Chem. 2013, 85, 7527−7534. (53) Atalay, S.; Yeh, L.-H.; Qian, S. Proton enhancement in an extended nanochannel. Langmuir 2014, 30, 13116−13120. (54) Kalwarczyk, T.; Tabaka, M.; Holyst, R. Biologistics−diffusion coefficients for complete proteome of Escherichia coli. Bioinformatics 2012, 28, 2971−2978. (55) Humphrey, W.; Dalke, A.; Schulten, K. VMD: Visual molecular dynamics. J. Mol. Graphics 1996, 14, 33−38. (56) Cignoni, P.; Callieri, M.; Corsini, M. Meshlab: An open-source mesh processing tool; Eurographics Italian Chapter Conference; Scarano, V., De Chiara, R., Erra, U., Eds.; 2008. (57) Geuzaine, C.; Remacle, J.-F. Gmsh: A 3-D finite element mesh generator with built-in pre- and post-processing facilities. Int. J. Numer. Meth. Eng. 2009, 79, 1309−1331. (58) Logg, A.; Wells, G. N.; Hake, J. DOLFIN: A C++/Python finite element library. In Automated Solution of Differential Equations by the Finite Element Method; Logg, A., Mardal, K.-A., Wells, G., Eds.; Springer: Berlin, Heidelberg, 2012. (59) Auriault, J.-L.; Boutin, C.; Geindreau, C. Numerical and Analytical Estimates for the Effective Diffusion Coefficient. In Homogenization of Coupled Phenomena in Heterogenous Media; WileyISTE: London, U.K., 2010. (60) Salgn, S.; Salgn, U.; Bahadr, S. Zeta potentials and isoelectric points of biomolecules: The effects of ion types and ionic strengths. Int. J. Electrochem. Sci. 2012, 7, 12404−12414. (61) Bensoussan, A.; Lions, J. L.; Papanicolaou, G. Asymptotic Analysis for Periodic Structures; American Mathematical Soc.: 2011. (62) Higgins, E. R.; Goel, P.; Puglisi, J. L.; Bers, D. M.; Cannell, M.; Sneyd, J. Modelling calcium microdomains using homogenisation. J. Theor. Biol. 2007, 247, 623−644. (63) Hashin, Z.; Shtrikman, S. Conductivity of polycrystals. Phys. Rev. 1963, 130, 129−133. (64) Garnett, J. Colours in metal glasses, in metallic films, and in metallic solutions. II. Philos. Trans. R. Soc., A 1905, 203, 385−420. (65) Young, P. G.; Beresford-West, T. B. H.; Coward, S. R. L.; Notarberardino, B.; Walker, B.; Abdul-Aziz, A. An efficient approach to converting three-dimensional image data into highly accurate computational models. Philos. Trans. R. Soc., A 2008, 366, 3155−3173. (66) Connolly, M. Solvent-accessible surfaces of proteins and nucleic acids. Science 1983, 221, 709−713. (67) Shorten, P. R.; McMahon, C. D.; Soboleva, T. K. Insulin transport within skeletal muscle transverse tubule networks. Biophys. J. 2007, 93, 3001−3007. (68) Sheinerman, F. B.; Norel, R.; Honig, B. Electrostatic aspects of protein-protein interactions. Curr. Opin. Struct. Biol. 2000, 10, 153− 159.

(69) Majhi, P. R.; Ganta, R. R.; Vanam, R. P.; Seyrek, E.; Giger, K.; Dubin, P. L. Electrostatically driven protein aggregation: βlactoglobulin at low ionic strength. Langmuir 2006, 22, 9150−9159. (70) Gilson, M. K.; Zhou, H.-X. Calculation of protein-ligand binding affinities. Annu. Rev. Biophys. Biomol. Struct. 2007, 36, 21−42. (71) Duan, C.; Majumdar, A. Anomalous ion transport in 2-nm hydrophilic nanochannels. Nat. Nanotechnol. 2010, 5, 848−852. (72) Smith, J. J.; Zharov, I. Ion transport in sulfonated nanoporous colloidal films. Langmuir 2008, 24, 2650−2654. (73) Moyne, C.; Murad, M. A. A two-scale model for coupled electro-chemo-mechanical phenomena and Onsager’s reciprocity relations in expansive clays: II Computational validation. Transp. Porous Media 2006, 63, 13−56. (74) Kuhlman, T. E.; Cox, E. C. Gene location and dna density determine transcription factor distributions in Escherichia coli. Mol. Syst. Biol. 2012, 8, 610. (75) Shorten, P. R.; Sneyd, J. A mathematical analysis of obstructed diffusion within skeletal muscle. Biophys. J. 2009, 96, 4764−4778. (76) Shkoller, S. An approximate homogenization scheme for nonperiodic materials. Comput. Math. Appl. 1997, 33, 15−34. (77) Blanc, X.; Le Bris, C.; Lions, P. L. A Possible Homogenization Approach for the Numerical Simulation of Periodic Microstructures with Defects. Milan J. Math. 2012, 80, 351−367. (78) Le Bris, C. Homogenization theory and multiscale numerical approaches for disordered media: Some recent contributions. ESAIM: ProcS 2014, 45, 18−31. (79) Israelachvili, J. N. Intermolecular and Surface Forces; Academic Press: San Diego, C.A., 2011. (80) Eisenberg, B. Mass action in ionic solutions. Chem. Phys. Lett. 2011, 511, 1−7. (81) Boschitsch, A. H.; Danilov, P. V. Formulation of a new and simple nonuniform size-modified Poisson-Boltzmann description. J. Comput. Chem. 2012, 33, 1152−1164. (82) Wang, N.; Zhou, S.; Kekenes-Huskey, P. M.; Li, B.; McCammon, J. A. Poisson-Boltzmann versus size-modified PoissonBoltzmann electrostatics applied to lipid bilayers. J. Phys. Chem. B 2014, 118, 14827−14832. (83) Li, B. Continuum electrostatics for ionic solutions with nonuniform ionic sizes. Nonlinearity 2009, 22, 811−833. (84) Valiullin, R.; Kortunov, P.; Karger, J.; Timoshenko, V. Concentration-dependent self-diffusion of liquids in nanopores: A nuclear magnetic resonance study. J. Chem. Phys. 2004, 120, 11804. (85) Straube, R.; Ridgway, D. Investigating the effects of molecular crowding on Ca2+ diffusion using a particle-based simulation model. Chaos 2009, 19, 037110. (86) Alekseev, A. E.; Reyes, S.; Selivanov, V. A.; Dzeja, P. P.; Terzic, A. Compartmentation of membrane processes and nucleotide dynamics in diffusion-restricted cardiac cell microenvironment. J. Mol. Cell. Cardiol. 2012, 52, 401−409. (87) Aronsen, J. M.; Swift, F.; Sejersted, O. M. Cardiac sodium transport and excitation-contraction coupling. J. Mol. Cell. Cardiol. 2013, 61, 11−19. (88) Towns, J.; Cockerill, T.; Dahan, M.; Foster, I.; Gaither, K.; Grimshaw, A.; Hazlewood, V.; Lathrop, S.; Lifka, D.; Peterson, G. D.; et al. XSEDE: Accelerating scientific discovery. Comput. Sci. Eng. 2014, 16, 62−74.

K

DOI: 10.1021/acs.jpcb.6b03887 J. Phys. Chem. B XXXX, XXX, XXX−XXX