Next Article in Journal
On Quantum Collapse as a Basis for the Second Law of Thermodynamics
Next Article in Special Issue
Fluctuation-Driven Transport in Biological Nanopores. A 3D Poisson–Nernst–Planck Study
Previous Article in Journal
“Over-Learning” Phenomenon of Wavelet Neural Networks in Remote Sensing Image Classifications with Different Entropy Error Functions
Previous Article in Special Issue
Taxis of Artificial Swimmers in a Spatio-Temporally Modulated Activation Medium
 
 
Font Type:
Arial Georgia Verdana
Font Size:
Aa Aa Aa
Line Spacing:
Column Width:
Background:
Article

Brownian Dynamics Computational Model of Protein Diffusion in Crowded Media with Dextran Macromolecules as Obstacles

1
Department of Material Science and Physical Chemistry, Barcelona University, 08028 Barcelona, Spain
2
Institute of Theoretical and Computational Chemistry (IQTC), Barcelona University, 08028 Barcelona, Spain
3
Department of Chemistry, University of Lleida (UdL), 25003 Lleida, Spain
*
Authors to whom correspondence should be addressed.
Entropy 2017, 19(3), 105; https://doi.org/10.3390/e19030105
Submission received: 22 December 2016 / Revised: 1 March 2017 / Accepted: 3 March 2017 / Published: 9 March 2017
(This article belongs to the Special Issue Nonequilibrium Phenomena in Confined Systems)

Abstract

:
The high concentration of macromolecules (i.e., macromolecular crowding) in cellular environments leads to large quantitative effects on the dynamic and equilibrium biological properties. These effects have been experimentally studied using inert macromolecules to mimic a realistic cellular medium. In this paper, two different experimental in vitro systems of diffusing proteins which use dextran macromolecules as obstacles are computationally analyzed. A new model for dextran macromolecules based on effective radii accounting for macromolecular compression induced by crowding is proposed. The obtained results for the diffusion coefficient and the anomalous diffusion exponent exhibit good qualitative and generally good quantitative agreement with experiments. Volume fraction and hydrodynamic interactions are found to be crucial to describe the diffusion coefficient decrease in crowded media. However, no significant influence of the hydrodynamic interactions in the anomalous diffusion exponent is found.

1. Introduction

The study of reaction and diffusion processes in biological media has been a challenging topic of recent research. Although single macromolecules are present in low concentrations, there is a high total concentration of macromolecules (such as proteins, polysaccharides, etc.) in cellular environments. In this context, macromolecular crowding can be defined as “macromolecular cosolutes that are nominally inert with respect to the reaction of interest” [1]. In general, cell cytosol presents an occupied volume fraction of 20%–30%, which means an approximate macromolecule concentration of 200–300 g/L. Moreover, macromolecular crowding is also relevant outside the cells [2] (e.g., blood plasma has a non-negligible 80 g/L protein concentration).
Since in vitro experiments are usually carried out at low concentrations (1–10 g/L), alternative approaches are necessary in order to evaluate the effect of macromolecular crowding in the thermodynamic and kinetic properties of the system. At the experimental level, high concentrations of crowding agents (usually dextran or ficoll macromolecules), which are considered to interact only by means of steric non-specific interactions, have been used to mimic the in vivo environment (usually called in vivo-like media) [3,4]. The experimental studies of macromolecule diffusion are mainly based on two experimental techniques: Fluorescence Correlation Spectroscopy (FCS) [5] and Fluorescence Recovery After Photobleaching (FRAP) [6,7].
FCS relies on the detection and temporal analysis of the fluorescence signal emitted from a small confocal detection volume. In these studies, fluorescent markers are specifically bound to the tracer molecule allowing a wide range of possible tracers (e.g., proteins, polymers, metal-complexes, etc.). The fluctuations in fluorescence intensity in a small region of the sample are recorded and used to calculate the temporal auto-correlation function, which allows determining the diffusion coefficient and the anomalous diffusion exponent of the tracer.
On the other hand, in FRAP, a small volume of the sample is lighted with a laser beam. As a result, the molecules in the lighted region become bleached and no fluorescence is exhibited. The diffusion coefficient and the anomalous diffusion exponent of the fluorescent molecule can be estimated by means of the velocity of fluorescence recovery in the bleached region.
In addition, computational studies have been performed in order to understand the effect of macromolecular crowding [8]. These studies use different approaches such as on-lattice Monte Carlo simulations [9] or off-lattice Brownian Dynamics (BD) simulations [10,11]. Recently, macromolecular crowding effect on diffusion has been modelled by means of atomistic models of cytoplasm using Molecular Dynamics [12].
Macromolecular crowding is crucial to describe macromolecular diffusion in biological media. Under these conditions, the well-known Einstein–Smoluchowski equation is no longer valid and three different difusion regimes are observed, as shown in Figure 1 [9,13]. At short times, the particles of the system have not collided yet and the diffusion coefficient ( D short ) remains constant. Generally, D short is not equal to the dilute solution diffusion coefficient because the motion of the particles is slowed down due to Hydrodynamic Interactions (HI) with the other particles. As a result of the inter-particle collisions, the diffusion coefficient decays until it reaches a new stationary state at long times ( D long ). In the intermediate regime, also known as anomalous regime, the mean square displacement of the particles ( < r 2 > ) becomes non linear over time and can be expressed as:
< r 2 > = ( 2 d ) Γ t α .
This is formally equivalent to having a diffusion coefficient that is not constant over time [14]:
D ( t ) = Γ t α 1 2 d ,
where α is the so-called anomalous diffusion exponent, d is the topological dimension of the system and Γ is a generalized transport coefficient (also known as anomalous diffusion coefficient).
The effect of HI is also a key factor in macromolecular dynamics [15]. These interactions emerge from the fact that the motion of different particles becomes correlated by means of solvent interactions. Since in BD the solvent is not explicitly included in the simulations, HI need to be included in the algorithm in an effective way. In general, two different approaches have been used to take into account the HI: those which rely on calculating the diffusion tensor [16] using the Rotne–Prager–Yamakaya (RPY) method [17,18], and those which benefit from the Tokuyama model [19,20,21].
The RPY method assumes a far-field approximation involving pairs of equal sized particles. The calculation of the diffusion tensor and its factorization is a computationally expensive procedure as it scales with the number of particles N as N 3 . Several approximations can be implemented in order to speed up the results of BD simulations [22]. These procedures are mainly based on avoiding the Cholesky decomposition of the diffusion tensor since it is the bottleneck of the simulation [23,24]. Recently, several efforts have been made to improve the RPY approach such as generalization to different-sized particles [25] or inclusion of many-body and near-field HI by means of the Durlofsky–Brady–Bossis method [26,27]. Moreover, several studies have shown that the inclusion of the short ranged HI or lubrication forces are crucial for describing macromolecular diffusion in crowded media since long ranged HI become screened [28,29].
Unlike the RPY approach, the Tokuyama model is a mean-field approximation for equal-sized soft core spheres that mimic the self-diffusion of biomolecules in solution. This model provides a better description of short range HI than the conventional RPY diffusion tensor. This procedure has started to be used in BD simulations in crowded media [30,31] because it is computationally cheaper than the RPY approach since it allows for introducing the HI contributions without calculating the diffusion tensor.
In the present paper, macromolecule diffusion in crowded media is studied by using BD simulations. Two different experimental systems have been modelled using a new dextran model involving an effective radius which allows accounting for their steric compression in crowded media. The HI are included using the Tokuyama model and their role in diffusion is discussed. Finally, the effect of non-specific interactions, obstacle size and HI on the diffusion coefficient and the anomalous diffusion exponent is analyzed.

2. Methodology and System Parametrization

2.1. BD Simulations with HI

The large amount of solvent molecules in the macromolecule solution makes all-atom Molecular Dynamics simulations computationally very expensive. In this context, Langevin equation [32] provides a suitable procedure since the solvent is implicitly included by adding a stochastic force in the classic Newton equations of motion which accounts for the collisions of the Brownian particles with the solvent. In BD time scales, it is more convenient to apply the over-damped limit [33] for which the equation of motion reads [34]:
r ( t + Δ t ) = r ( t ) D Δ t k B T V ( r , t ) + 2 D Δ t ξ ( t ) ,
where ξ is a vector of 3N Gaussian random numbers with zero mean and unit variance.
Since we wish to model highly concentrated macromolecular solutions, we have chosen a coarse-grained approach where each macromolecule is modelled as a single sphere using a proper effective radius. In order to avoid overlapping, a harmonic pairwise repulsion is applied which acts when the distance between two macromolecules is smaller than the sum of their radii:
V ij ( r i , r j ) = 1 2 k ( d ij R i j ) 2 , d i j < R i j , 0 , d i j R i j ,
where R i j is the sum of the radius of the interacting particles i and j, d i j is the distance between i and j particles and k is a parameter accounting for the stiffness of the potential [35]. The cubic simulation box is used as outlined in Figure 2. The length of the simulation box is adapted to the dextran size, which is taken as 18 nm, 38 nm and 77 nm for D5, D50 and D400 dextran molecules, respectively (see the Table in Section 3). This allowed the number of particles to range between 50 and 200 for all the dextran sizes. Periodic Boundary Conditions (PBC) are applied in all directions. The Mean Square Displacement of the particles is calculated after a thermalisation time of 10 ns to avoid the effect of unrealistic particle overlapping due to the random initial configuration. All the performed simulations are 1000 ns long with a time step of 0.1 ns. The system temperature is 298.15 K. For each studied system, a minimum of 400 different BD realizations are averaged over. D long values have been calculated by fitting the Einstein–Smoluchowski equation to the simulated mean square displacement ( < r 2 > ) in the long-time regime.
HI are included in the BD equation of motion using the Tokuyama model [19]. In this approach, a Fokker–Planck equation for the single-particle distribution function is proposed which is coupled with the Navier–Stokes equation. This equation is analytically solved at short times considering a stationary configuration of equal-sized Brownian particles. Correlation effects and direct interactions are neglected. As a result, the diffusion coefficient of the particle at short times ( D short ) is found to be a function of the volume fraction ϕ , which has been calculated using [5,36,37,38]
ϕ = 4 π ( R T 3 + N O R O 3 ) 3 L 3 ,
where R T is the radius of the tracer protein, and R O is the radius of the dextran obstacles, N O is the number of dextran obstacles and L is the side of the simulation box. Following the Tokuyama approach, the short time diffusion coefficient is related to the volume fraction by the equation
D short ( ϕ ) = D 0 [ 1 + H ( ϕ ) ] ,
where D 0 is the diffusion coefficient at dilute solution and H ( ϕ ) is the contribution due to the HI valid in the short-time regime. H ( ϕ ) is given by
H ( ϕ ) = 2 b 2 1 b c 1 + 2 c b c ( 2 + c ) ( 1 + c ) ( 1 b + c ) ,
where b = 9 8 ϕ and c = 11 16 ϕ . Recently, Tokuyama et al. [20] have obtained good agreement for the diffusion coefficient at long times ( D long ) between Equation (6) and Molecular Dynamics with explicit solvent. It is also worth mentioning that Tokuyama [21] has proposed an analytical expression of the D long in terms of the volume fraction:
D long ( ϕ ) = D short ( ϕ ) 1 + κ D short ( ϕ ) D 0 ϕ ϕ c 1 ϕ ϕ c 2 ,
where κ and ϕ c are parameters settled to κ = 2 . 0 and ϕ c = 1 . 09 . In this work, simulations with and without HI are compared. In the simulations without HI, the diffusion coefficient D in Equation (3) corresponds to the diffusion coefficient in dilute solution D 0 . HI are included by means of the Tokuyama method [20] replacing D in Equation (3) by D short calculated using Equations (6) and (7).
Tokuyama equations for D short and D long are derived assuming equal-sized spheres. Since our systems contain two different-sized particles, the tracer protein and the dextran obstacles, the use of Equations (6) and (7) for D short involves an approximation. D long , however, is calculated using direct simulated mean square displacements at long times, so that the difference size effect is included in a straightforward way via the inter-particle maximum approach R i j in Equation (4). The resulting D long values are then compared to those predicted by Equation (8) in Section 2.3.

2.2. Effect of the Interaction Potential

The effect of the interaction potential stiffness (Equation (4)) has been also analysed to ensure that the selected constant ( k = 50 , 000 J·mol 1 ·nm 2 ) was enough to prevent particles to overlap. With this aim, preliminary computations were performed which are shown in Figure 3. No significant difference was observed between the three proposed values for k ( k = 1 × 10 5 , k = 5 × 10 4 and k = 2 . 5 × 10 4 J·mol 1 ·nm 2 ), which means that the potential was rigid enough to avoid overlapping.

2.3. Effect of Long Range Hydrodynamic Interactions

The effect of the HI in macromolecular diffusion in crowded media has been evaluated by selecting a simple system of equal sized spheres at different volume fractions. We have performed calculations comparing the conventional RPY tensor [17,18] and Tokuyama method [19,20,21].
Long range HI experience a slow decay with inter-particle distances, a well-known effect in electrostatic interactions. This fact could lead to size-dependent results and to unphysical non-definite positive diffusion tensors [39]. Therefore, in principle, procedures like Ewald sums [40] should be implemented. In order to avoid this problem, we performed simulations at different system sizes. We found that for N > 50, the finite size box effect becomes rather small [41]. It was also checked that the diffusion tensor was always definite positive at every time step. This agrees with previous studies [28,29] indicating that long range HI become considerably screened at high volume fractions.
The obtained values for D long without HI, with HI described by RPY diffusion tensor, and with HI using the Tokuyama method are shown in Figure 4. We can observe that the simulations without HI and with HI using the RPY diffusion tensor provide almost identical D long values. This fact supports previously reported results [28,29] highlighting the long range HI screening in crowded media. Conversely, significant differences arise in implementing the Tokuyama method, which account for short range HI (usually called lubrication forces). A marked decay of the diffusion coefficient is observed, in agreement with the relevance of lubrication forces reported in previous works. It is also important to highlight the almost perfect agreement, for a system of equal-sized particles, between the Tokuyama D long analytical prediction (Equation (8)) and BD simulations using the Tokuyama model, as Tokuyama previously reported [20].

2.4. Effect of the Difference in Size between Tracer and Obstacles

Tokuyama analytical expressions were derived assuming the same size for all of the particles. However, in our experimental systems, the tracer particle (a protein) and the obstacles (dextran macromolecules) present different sizes. In order to check the accuracy of the Tokuyama prediction for D long , we have performed simulations that implement the difference in size by a suitable modification of the minimum inter-particle distance R i j in the interaction potential (Equation (4)). The Tokuyama effective diffusion coefficient D short is still used to describe the short range HI forces since, to our knowledge, no expression for D short is available yet for different sized particles.
D long resulting from BD simulations and the ones calculated using the Tokuyama analytical expression (Equation (8)) versus volume fraction are depicted in Figure 5. It can be observed that when the tracer (with radius R = 2 . 33 nm) and obstacle ( R = 2 . 9 nm) present similar sizes, very good agreement between Tokuyama equations and simulations is found, as shown in Figure 5a. In contrast, if the tracer particle is far bigger ( R = 4 . 90 nm) than the obstacles ( R = 1 . 2 nm), important differences arise. As shown in Figure 5b, D long values predicted using the Tokuyama analytic method clearly overestimate the ones obtained by BD computations, with differences up to 60% between both approaches. In conclusion, Tokuyama analytic equations provide a reliable prediction for D long only for systems of similarly-sized spheres so that BD computations are unavoidable otherwise.

2.5. Dextran Model

Here, we assume that each macromolecule is a single sphere. In similar coarse-grained models [12,31], the radius of the particles is calculated using their dilute solution diffusion coefficient via the Stokes–Einstein equation. In these approaches, the dilute solution diffusion coefficient is estimated using HYDROPRO software (Version 10) [42]. HYDROPRO needs as input the atomic coordinates of the macromolecule, which is usually obtained from crystallographic data or NMR spectra. HYDROPRO uses this atomistic structure to compute its hydrodynamic properties such as the diffusion coefficient in dilute solution. Unfortunately, this information is not available for dextran macromolecules yet.
However, experimental hydrodynamic radius ( R H ) obtained using quasi-elastic light scattering [43], and radius of gyration ( R G ) obtained from capillary viscometry [44] are available. We have fitted R H and R G versus molecular weight data to the power law:
R x = K · M w γ ,
where R x can be either R G or R H . K and γ are the best fitted parameters reported in Table 1. Unexpectedly, we have found that these experimental radii lead to high volume fractions at experimental concentrations (50–300 g/L). This behaviour increases dramatically with dextran size, resulting in unphysical volume fractions larger than one for the biggest sizes. This can be explained taking into account the very complex branched structure of dextran macromolecules [45,46]. As dextran macromolecular size increases, the molecular structure becomes more branched. Moreover, the structure of dextran is flexible allowing steric compression in concentrated solution.
An alternative to calculating the volume fraction at high dextran concentrations is the use of an average dextran specific volume (ν = 0.625 cm 3 /g) [47]. Considering the dextran molecules as compact spheres with a radius R c :
R c = 3 ν 4 π N A 3 M W 1 3 = K · M w γ ,
where N A is the Avogadro number, γ = 1 3 and K = 3 ν 4 π N A 3 = 0 . 063 nm · Da 1 3 . However, previous studies using Equation (10) showed bad agreement with the experimental diffusion coefficients [5,6,7] obtained using FRAP and FCS. More specifically, the obtained results exhibit a slower decay of the diffusion coefficient with dextran concentration than that experimentally observed. This was caused by the neglect of solvatation effects resulting in volume fractions that were too low.
In summary, a suitable effective radius ( R eff ) lying in between the hydrodynamic and the compact radius is necessary for dextran macromolecules at high concentration. Such an effective radius has been determined by performing BD simulations scanning the physically meaningful values for K and γ parameters. In these simulations, K ranged from 0.043 to 0.063 nm·Da γ while γ ranged from 1 3 to 0.445. Among the performed calculations, the best agreement between experimental [5,6,7] and computed D long was obtained with γ = 0 . 387 and K = 0.045 nm·Da γ .
It is also worth noting that the chosen parameters for R eff are closer to those of R c than to those of R H (Figure 6). This means that steric compression is crucial to properly describe dextran size in crowded media. Moreover, the large difference between R H and R c at high dextran molecular weight reveals the importance of the steric compression as dextran size increases.

3. Results and Discussion

Two different experimental systems have been selected and computationally modeled in the present work. The first one is the FCS study of Streptavidin protein ( R H = 4 . 90 nm ) diffusion at different concentrations of six distinct dextran macromolecules performed by Banks et al. [5]. The second one uses FRAP to record α-Chymiotrypsin protein ( R H = 2 . 33 nm ) diffusion [6,7] at different concentrations of three different-sized dextran macromolecules acting as crowding agents. We have chosen the dextran molecules D5, D50 and D400 as obstacles since they have been used in both experimental studies, allowing better comparison between them. The characteristic parameters of the used dextran macromolecules are shown in Table 2. In FCS, the fluorescence fluctuations of a sample of tracer proteins are associated with the tracer auto-correlation function [5]. This allows for experimentally measuring the anomalous exponent α and a characteristic residence time τ D . τ D is defined as the time which the tracer particle needs to traverse the characteristic length of the detection volume. A similar approach is followed in the analysis of the fluorescence recovery measured by FRAP [6,7]. α and τ D are used to calculate the effective diffusion coefficient ( D eff ) at a time equal to τ D . The experimentally found τ D , around 1–3 ms [5,6], is much longer than the computational times at which D long is calculated. As a consequence, the computed values of D long should reasonably correspond to the experimentally obtained D eff .

3.1. Long Time Diffusion Coefficient

As mentioned above, diffusion in crowded media has three different temporal regimes over time (Figure 1). In order to compare computations and experiments, it is necessary to calculate the diffusion coefficient in the long time regime ( D long ). In the selected experimental systems, the diffusion of a tracer protein is studied at different concentrations of three different-sized dextran macromolecules (D5, D50 and D400), which act as inert obstacles. Two different approaches have been applied in the simulations. In the first one, the HI are not considered and thus only a pairwise repulsive harmonic potential is acting. In the second one, HI are taken into account using the Tokuyama method.
Figure 7 and Figure 8 show the effect of macromolecular crowding in the diffusion of Streptavidin and α-Chymiotrypsin, respectively. D long is normalized using the diffusion coefficient at infinite dilution ( D 0 ) of the particle. The simulations including HI using the Tokuyama model exhibit, in general, better qualitative and quantitative agreement with the experimental data for all the dextran sizes studied. As a consequence, the inclusion of HI is crucial to describe the reduction of D long experimentally observed. This is in concordance with the importance of the short range HI in crowded media previously reported in the literature [29]. It is also worth mentioning that D long always decays as obstacle concentration increases, as a result of the increase of particle collisions, which hinder protein diffusion.
BD simulations show a clear effect of the tracer particle size, which is different in both systems. Streptavidin shows a faster decay with dextran concentration than α-Chymiotrypsin. This is logical since Streptavidin (with R = 4 . 9 nm) is quite larger than α-Chymiotrypsin ( R = 2 . 33 nm). As a result, the loss of mobility due to an increase in the number of collisions of Streptavidin is larger than in the case of α-Chymiotrypsin. On the other hand, BD simulations show a small contribution of the dextran size in D long for the same dextran concentration. This is probably due to the limitations of Tokuyama approach, which assumes equal sizes for the particles in the derivation of the equations. As a consequence, Tokuyama equations only account for the dependence on the volume fraction, but not on the obstacle size. This is the case of α-Chymiotrypsin, but not that of Streptavidin, for which D long clearly decays faster as the obstacle size increases. This fact points out the interest of future extensions of the Tokuyama approach to systems of different-sized particles.

3.2. Anomalous Diffusion Exponent

The anomalous diffusion exponents α for the two studied diffusing proteins are depicted in Figure 9 and Figure 10. In concordance with the obtained results for D long , the α exponent also decays with dextran concentration. A good quantitative prediction of the experimental data is obtained for both proteins except for the smallest dextran size.
In general, the obtained results with and without HI are very similar indicating that, unlike D long , the α exponent is not significantly influenced by HI. Since the anomalous diffusion exponent is related to the transition rate between D short and D long regimes, this means that although HI slow down D short and D long , they do not affect the rate of the transition. One possible explanation is that the collision rate of particles, the driving mechanism of the transition, does not change substantially with the inclusion of HI in the BD simulations of crowded media. For Streptavidin, it is also observed that the anomalous exponent decreases with dextran size. However, the non-monotonous dependence of the anomalous exponent for α-Chymiotrypsin is not clear, given the high experimental error for this parameter. Probably, more experiments should be necessary in order to fully clarify this point.

4. Conclusions

A Brownian Dynamics computational model for two experimental systems of diffusing proteins (Streptavidin and α-Chymiotrypsin) in crowded media has been proposed. In both cases, dextran macromolecules have been used as crowding agents. Dextran macromolecules have been modelized as spheres with effective radii accounting for macromolecular compression. The obtained long time diffusion coefficients ( D long ) are, in general, close to the ones experimentally observed. In the analysis of D long , the inclusion of the HI ends up being fundamental to explain the experimental decay. This means that steric effects and HI are the two main contributions to the decrease of D long in crowded media. However, no significant influence of HI in the anomalous diffusion exponent is detected. Although the qualitative behaviour is properly modelled using the Tokuyama approach, in some cases, the obtained results by simulation are still far from the experimental ones. This could be mainly due to two different reasons.
The first one is that the Tokuyama model assumes equal-sized spheres and our systems are heterogeneous. Therefore, this method is less accurate when the difference in size between the tracer protein and obstacles particles increases. A generalization of the Tokuyama model for different spheres is not available yet and it would be a promising procedure to improve our results.
The second factor that can be responsible for the difference between simulation and experiments is the procedure used to model dextran molecules. This method does not take into account that the degree of compression of the dextran macromolecules increases as their concentration increases. This suggests the possibility to improve the dextran description by changing the harmonic potential for a new potential able to allow dextran radius to evolve from its hydrodynamic radius to its compact radius as the concentration increases [48]. This new potential could also account for the entanglement of the polymer branches, allowing some overlap between the macromolecules. This effect could be also relevant in the diffusion properties of macromolecules in crowded media.

Acknowledgments

We acknowledge the financial support from: the Spanish Ministry of Science and Innovation (project CTM2012-39183 and CTM2016-78798-C2-1-P) and Generalitat de Catalunya (Grants 2014SGR1017, 2014SGR1132 and XrQTC). Sergio Madurga and Francesc Mas acknowledge the funding of the project 8SEWP-HORIZON 2020 (692146).

Author Contributions

Pablo M. Blanco, Josep Lluís Garcés, Sergio Madurga and Francesc Mas designed the BD code and the BD simulations. Pablo M. Blanco and Mireia Via developed and tested the computational code and performed the BD simulations. Pablo M. Blanco, Sergio Madurga and Francesc Mas discussed and proposed the new dextran macromolecular model. Pablo M. Blanco, Josep Lluís Garcés, Sergio Madurga and Francesc Mas collaborated to write the present paper. All authors have read and approved the final manuscript.

Conflicts of Interest

The authors declare no conflict of interest. The founding sponsors had no role in the design of the study; in the collection, analyses, or interpretation of data; in the writing of the manuscript, and in the decision to publish the results.

Abbreviations

The following abbreviations are used in this manuscript:
BDBrownian Dynamics
FCSFlourescence Correlation Spectroscopy
FRAPFluorescence Recovery After Photobleaching
HIHydrodynamic Interactions
RPYRotne–Prager–Yamakawa

References

  1. 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. [Google Scholar] [CrossRef] [PubMed]
  2. Ellis, R.J. Macromolecular crowding: Obvious but underappreciated. TRENDS BioChem. Sci. 2001, 26, 597–604. [Google Scholar] [CrossRef]
  3. Pastor, I.; Pitulice, L.; Balcells, C.; Vilaseca, E.; Madurga, S.; Isvoran, A.; Mas, F. Effect of crowding by Dextrans in enzymatic reactions. Biophys. Chem. 2014, 185, 8–13. [Google Scholar] [CrossRef] [PubMed]
  4. Balcells, C.; Pastor, I.; Pitulice, L.; Hernández, C.; Via, M.; Garcés, J.L.; Madurga, S.; Vilaseca, E.; Isvoran, A.; Cascante, M.; et al. Macromolecular Crowding upon in-vivo-Like Enzyme-Kinetics: Effect of Enzyme-Obstacle Size Ratio. New J. Chem. 2015, 24, 3–16. [Google Scholar]
  5. Banks, D.S.; Fradin, C. Anomalous Diffusion of Proteins Due to Molecular Crowding. Biophys. J. 2005, 89, 2960–2971. [Google Scholar] [CrossRef] [PubMed]
  6. Pastor, I.; Vilaseca, E.; Madurga, S.; Garcés, J.L.; Cascante, M.; Mas, F. Diffusion of α-Chymiotrypsin in Solution-Crowded Media. A Flourescence Recovery After Photobleaching Study. J. Phys. Chem. 2010, 114, 4028–4034. [Google Scholar] [CrossRef] [PubMed]
  7. Vilaseca, E.; Pastor, I.; Isvoran, A.; Madurga, S.; Garcés, J.L.; Mas, F. Diffusion in macromolecular crowded media: Monte Carlo simulation of obstructed diffusion vs. FRAP experiments. Theor. Chem. Acc. 2011, 128, 795–805. [Google Scholar] [CrossRef]
  8. Elcock, A.H. Models of macromolecular crowding effects and the need of quantitative comparisons with experiment. Curr. Opin. Struct. Biol. 2010, 20, 196–206. [Google Scholar] [CrossRef] [PubMed]
  9. Vilaseca, E.; Isvoran, A.; Madurga, S.; Pastor, I.; Garceés, J.L.; Mas, F. New insights into diffusion in 3D crowded media by Monte Carlo simulations: Effect of size, mobility and spatial distribution of obstacles. Phys. Chem. Chem. Phys. 2011, 13, 7396–7407. [Google Scholar] [CrossRef] [PubMed]
  10. Hasnain, S.; McClendon, C.L.; Hsu, M.T.; Jacobson, M.P.; Bandyopadhyay, P. A New Coarse-Grained Model for E. coli Cytoplasm: Accurate Calculation of the Diffusion Coefficient of Proteins and Observation of Anomalous Diffusion. PLoS ONE 2014, 9. [Google Scholar] [CrossRef] [PubMed]
  11. Gomez, D.; Klumpp, S. Biochemical reactions in crowded environments: Revisiting the effects of volume exclusion with simulations. Front. Phys. 2015. [Google Scholar] [CrossRef]
  12. Yu, I.; Mori, T.; Ando, T.; Harada, R.; Jung, J.; Sugita, Y.; Feig, M. Biomolecular interactions modulate macromolecular structure and dynamics in atomistic model of a bacterial cytoplasm. eLife 2016, 5, e19274. [Google Scholar] [CrossRef] [PubMed]
  13. Saxton, M.J. Anomalous Diffusion Due to Obstacles: A Monte Carlo Study. Biophys. J. 1994, 66, 394–401. [Google Scholar] [CrossRef]
  14. Bauchaud, J.P.; Georges, A. Anomalous diffusion in disordered media: Statistical mechanics, model and physical application. Phys. Rep. 1990, 185, 127–293. [Google Scholar] [CrossRef]
  15. Elcock, A.H. Molecule-Centered Method for Accelerating the Calculation of Hydrodynamic Interactions in Brownian Dynamics Simulations Containing Many Flexible Biomolecules. J. Chem. Theory Comput. 2013, 9, 3224–3239. [Google Scholar] [CrossRef] [PubMed]
  16. Ermak, D.L.; McCammon, J.A. Brownian dynamics with hydrodynamic interactions. J. Chem. Phys. 1978, 69, 1352–1360. [Google Scholar] [CrossRef]
  17. Rotne, J.; Prager, S. Variational treatment of hydrodynamic interaction in polymers. J. Chem. Phys. 1969, 50, 4831–4837. [Google Scholar] [CrossRef]
  18. Yamakawa, H.J. Transport properties of polymer chains in dilute solution: Hydrodynamic interaction. Chem. Phys. 1970, 53, 436–443. [Google Scholar] [CrossRef]
  19. Tokuyama, M.; Oppenheim, I. Dynamics of hard-sphere suspensions. Phys. Rev. E 1994, 50, R16. [Google Scholar] [CrossRef]
  20. Tokuyama, M.; Moriki, T.; Kimura, Y. Self-diffusion of biomolecules in solution. Phys. Rev. E 2011, 83, 081402. [Google Scholar] [CrossRef] [PubMed]
  21. Tokuyama, M. Mean-field theory of glass transitions. Physica A 2006, 364, 23–62. [Google Scholar] [CrossRef]
  22. Rodríguez, R.; Hernández, J.G.; García, J. Comparison of Brownian dynamics algorithms with hydrodynamic interactions. J. Chem. Phys. 2011, 135, 084116. [Google Scholar]
  23. Fixman, M. Implicit algorithm for Brownian Dynamics of Polymers. Macromolecules 1986, 19, 1195–1204. [Google Scholar] [CrossRef]
  24. Geyer, T.; Winter, U. An O(N2) approximation for hydrodynamic interactions in Brownian dynamics simulations. J. Chem. Phys. 2009, 130, 114905. [Google Scholar] [CrossRef] [PubMed]
  25. Zuk, P.J.; Wajnryb, E.; Mizerski, K.A.; Szymczak, P. Rotne–Prager–Yamakawa approximation for different-sized particles in application to macromolecular bead model. J. Fluid. Mech. 2014, 741, R5. [Google Scholar] [CrossRef]
  26. Durlofsky, L.; Brady, J.F.; Bossis, G. Dynamic simulation of hydrodynamically interacting particles. J. Fluid Mech. 1987, 180, 21–49. [Google Scholar] [CrossRef]
  27. Brady, J.F.; Bossis, G. Stokesian Dynamics. Annu. Rev. Fluid Mech. 1988, 20, 111–157. [Google Scholar] [CrossRef]
  28. Phillips, R.J.; Brady, J.F.; Bossis, G. Hydrodynamic transport properties of hard-sphere dispersions. I. Suspensions of freely mobile particles. Phys. Fluids 1988, 31, 3462–3472. [Google Scholar] [CrossRef]
  29. Ando, T.; Chow, E.; Skolnick, J. Dynamic simulation of concentrated macromolecular solutions with screened long-range hydrodynamic interactions: Algorithm and limitations. J. Chem. Phys. 2013, 139, 121922. [Google Scholar] [CrossRef] [PubMed]
  30. Sun, J.; Weinstein, H. Toward realistic modeling of dynamic processes in cell signaling: Quantification of macromolecular crowding effects. J. Chem. Phys. 2007, 127, 155105. [Google Scholar] [CrossRef] [PubMed]
  31. Mereghetti, P.; Wade, R.C. Atomic Detail Brownian Dynamics Simulations of Concentrated Protein Solutions with a Mean Field Treatment of Hydrodynamic Interactions. J. Chem. Phys. 2012, 116, 8523–8533. [Google Scholar] [CrossRef] [PubMed]
  32. Gillespie, D.T. The multivariate Langevin and Fokker–Planck equations. Am. J. Phys. 1996, 64, 1246–1257. [Google Scholar] [CrossRef]
  33. Mori, H. Transport, Collective Motion and Brownian Motion. Prog. Theor. Phys. 1965, 33, 423–455. [Google Scholar] [CrossRef]
  34. Dhont, J.K.G. Brownian motion of non-interacting particles. In An Introduction to Dynamics of Colloids; Elsevier: Amsterdam, The Netherlands, 1996; Chapter 2. [Google Scholar]
  35. Schöneberg, J.; Noé, F. ReaDDy—A Software for Particle-Based Reaction-Diffusion Dynamics in Crowded Cellular Environments. PLoS ONE 2013, 8, e74261. [Google Scholar] [CrossRef] [PubMed]
  36. Roosen-Runge, F.; Hennig, M.; Zhang, F.; Jacobs, R.M.; Sztucki, M.; Schober, H.; Seydel, T.; Schreiber, F. Protein self-diffusion in crowded solutions. Proc. Natl. Acad. Sci. USA 2011, 108, 11815–11820. [Google Scholar] [CrossRef] [PubMed]
  37. Qin, S.; Zhou, H.X. Atomistic modeling of macromolecular crowding predicts modest increases in protein folding and binding stability. Biophys. J. 2009, 97, 12–19. [Google Scholar] [CrossRef] [PubMed]
  38. Sastry, S.; Truskett, T.M.; Debenedetti, P.G.; Torquato, S.; Stillinger, F.H. Free volume in the hard sphere liquid. Mol. Phys. 1998, 95, 289–297. [Google Scholar] [CrossRef]
  39. Brady, J.F.; Phillips, R.J.; Lester, J.C.; Bossis, G. Dynamic simulation of hydrodynamically interacting suspensions. J. Fluid Mech. 1988, 195, 257–280. [Google Scholar] [CrossRef]
  40. Beenakker, C.W.J. Ewald sum of the Rotne–Prager tensor. J. Chem. Phys. 1986, 85, 1581–1582. [Google Scholar] [CrossRef]
  41. Ladd, A.J. Hydrodynamic transport coefficients of random dispersions of hard spheres. J. Chem. Phys. 1990, 93, 3484–3494. [Google Scholar] [CrossRef]
  42. Ortega, A.; Amorós, D.; de la Torre, J.G. Prediction of hydrodynamic and other solution properties of rigid proteins from atomic-and residue-level models. Biophys. J. 2011, 101, 892–898. [Google Scholar] [CrossRef] [PubMed]
  43. Armstrong, J.K.; Wenby, B.; Meiselman, H.J.; Fisher, T.C. The Hydrodynamic Radii of Macromolecules and Their Effect on Red Blood Cell Aggregation. Biophys. J. 2004, 87, 4259–4270. [Google Scholar] [CrossRef] [PubMed]
  44. Fundueanu, G.; Nastruzzi, C.; Carpov, A.; Desbrieres, J.; Rinaudo, M. Physico-chemical characterization of Ca-alginate microparticles produced with different methods. Biomaterials 1999, 20, 1427–1435. [Google Scholar] [CrossRef]
  45. Sloan, J.W.; Alexander, B.H.; Lohmar, R.L.; Wolff, I.A.; Rist, C.E. Determination of dextran structure by periodate oxidation techniques. J. Am. Chem. Soc. 1954, 76, 4429–4434. [Google Scholar] [CrossRef]
  46. Sabatie, J.; Choplin, L.; Moan, M.; Doublier, J.L.; Paul, F.; Monsan, P. The effect of synthesis temperature on the structure of dextran NRRL B 512F. Carbohydr. Polym. 1988, 9, 87–101. [Google Scholar] [CrossRef]
  47. Schaefer, D.W. A unified model for the structure of polymers in semidilute solution. Polymer 1984, 25, 387–394. [Google Scholar] [CrossRef]
  48. Vilaseca, P.; Franzese, G. Softness dependence of the anomalies for the continuous shouldered well potential. J. Chem. Phys. 2010, 133, 084507. [Google Scholar] [CrossRef] [PubMed]
Figure 1. Time evolution of the diffusion coefficient obtained in a system with volume fraction of ϕ = 0.17. The red line is the result of averaging 400 BD simulations of a tracer particle diffusing among 68 obstacles with radius 4.9 and 1.2 nm respectively. The simulation time is 250 ns and the length of the simulation box is 18 nm.
Figure 1. Time evolution of the diffusion coefficient obtained in a system with volume fraction of ϕ = 0.17. The red line is the result of averaging 400 BD simulations of a tracer particle diffusing among 68 obstacles with radius 4.9 and 1.2 nm respectively. The simulation time is 250 ns and the length of the simulation box is 18 nm.
Entropy 19 00105 g001
Figure 2. Snapshot of one of the performed dynamics. A protein (in red) diffuses among dextran molecules (in yellow) which act as crowding agents.
Figure 2. Snapshot of one of the performed dynamics. A protein (in red) diffuses among dextran molecules (in yellow) which act as crowding agents.
Entropy 19 00105 g002
Figure 3. Decrease of the normalised diffusion coefficient at long times ( D long ) with the volume fraction at three different values of the stiffness constant k of the interaction potential. BD simulation perfomed: (a) without HI; (b) with HI using the Tokuyama model. The results show that the potential is rigid enough to avoid overlapping. Continuous lines are only to guide the lecturer.
Figure 3. Decrease of the normalised diffusion coefficient at long times ( D long ) with the volume fraction at three different values of the stiffness constant k of the interaction potential. BD simulation perfomed: (a) without HI; (b) with HI using the Tokuyama model. The results show that the potential is rigid enough to avoid overlapping. Continuous lines are only to guide the lecturer.
Entropy 19 00105 g003
Figure 4. Comparison of the diffusion coefficient calculated using BD without HI (red triangles), BD with HI using the conventional RPY diffusion tensor (blue squares), BD with HI using the Tokuyama model (cyan circles) and Tokuyama analytical expression (Equation (8), orange discontinuous line). The results show the important role of short range HI (lubrication forces) in the diffusion reduction of the diffusion coefficient, while the effect of the of long range HI is very low. It is also worth noting the good agreement between the Tokuyama analytical prediction of D long and the BD simulation with HI using Tokuyama model for D short . The lines are to guide the lecturer.
Figure 4. Comparison of the diffusion coefficient calculated using BD without HI (red triangles), BD with HI using the conventional RPY diffusion tensor (blue squares), BD with HI using the Tokuyama model (cyan circles) and Tokuyama analytical expression (Equation (8), orange discontinuous line). The results show the important role of short range HI (lubrication forces) in the diffusion reduction of the diffusion coefficient, while the effect of the of long range HI is very low. It is also worth noting the good agreement between the Tokuyama analytical prediction of D long and the BD simulation with HI using Tokuyama model for D short . The lines are to guide the lecturer.
Entropy 19 00105 g004
Figure 5. Comparison between the normalized D long obtained using Tokuyama analytic equations (Equations (6)–(8)) and the ones resulting from BD simulations at different volume fractions. In (a), the tracer particle (with radius R = 2 . 33 nm) and the obstacles ( R = 2 . 9 nm) present a similar size. D long values obtained by Tokuyama method (orange circles) and the ones obtained by BD simulations (red squares) clearly agree. In contrast, in (b), the tracer particle ( R = 4 . 90 nm) is far bigger than the obstacles ( R = 1 . 2 nm). In this case, the Tokuyama method (blue circles) clearly overestimates D long values obtained by BD simulations (cyan squares).
Figure 5. Comparison between the normalized D long obtained using Tokuyama analytic equations (Equations (6)–(8)) and the ones resulting from BD simulations at different volume fractions. In (a), the tracer particle (with radius R = 2 . 33 nm) and the obstacles ( R = 2 . 9 nm) present a similar size. D long values obtained by Tokuyama method (orange circles) and the ones obtained by BD simulations (red squares) clearly agree. In contrast, in (b), the tracer particle ( R = 4 . 90 nm) is far bigger than the obstacles ( R = 1 . 2 nm). In this case, the Tokuyama method (blue circles) clearly overestimates D long values obtained by BD simulations (cyan squares).
Entropy 19 00105 g005
Figure 6. Radius of gyration ( R G , red circles) taken from [44] and hydrodynamic radius ( R H , green triangles) taken from [43] versus molecular weight. The fittings corresponding to the power law (Equation (9)) are represented with red and green lines. The compact radius ( R c , blue) and the effective radius ( R eff , purple line) versus molecular weight are also depicted. The effective radius for the chosen dextran obstacles (purple squares) in our computations is also plotted. The lines are plotted only to guide the lecturer.
Figure 6. Radius of gyration ( R G , red circles) taken from [44] and hydrodynamic radius ( R H , green triangles) taken from [43] versus molecular weight. The fittings corresponding to the power law (Equation (9)) are represented with red and green lines. The compact radius ( R c , blue) and the effective radius ( R eff , purple line) versus molecular weight are also depicted. The effective radius for the chosen dextran obstacles (purple squares) in our computations is also plotted. The lines are plotted only to guide the lecturer.
Entropy 19 00105 g006
Figure 7. Decay of D long corresponding to Streptavidin protein versus dextran concentration. Three different sizes of dextran obstacles: (a) D5; (b) D50; and (c) D400 have been used. The obtained results including the HI with Tokuyama model exhibit better agreement in general with the experimental data showing the relevance of HI in the D long values in crowded media. The lines are only for guiding the lecturer.
Figure 7. Decay of D long corresponding to Streptavidin protein versus dextran concentration. Three different sizes of dextran obstacles: (a) D5; (b) D50; and (c) D400 have been used. The obtained results including the HI with Tokuyama model exhibit better agreement in general with the experimental data showing the relevance of HI in the D long values in crowded media. The lines are only for guiding the lecturer.
Entropy 19 00105 g007
Figure 8. α-Chymiotrypsin protein D long decay with dextran concentration. Three different sizes of dextran obstacles: D5 (a); D50 (b); and D400 (c) have been used. The obtained results including HI with Tokuyama model exhibit in general better agreement with the experimental data showing the relevance of HI in the D long values in crowded media. Continuous lines are only to guide the lecturer.
Figure 8. α-Chymiotrypsin protein D long decay with dextran concentration. Three different sizes of dextran obstacles: D5 (a); D50 (b); and D400 (c) have been used. The obtained results including HI with Tokuyama model exhibit in general better agreement with the experimental data showing the relevance of HI in the D long values in crowded media. Continuous lines are only to guide the lecturer.
Entropy 19 00105 g008
Figure 9. Anomalous diffusion exponent α of Streptavidin versus dextran concentration for three different sizes of dextran obstacles: (a) D5; (b) D50; and (c) D400. The obtained results exhibit good quantitative agreement with the experimental data except for the smallest dextran size D5. In general, no significant influence of HI in the α exponent is obtained. The lines are to guide the lecturer.
Figure 9. Anomalous diffusion exponent α of Streptavidin versus dextran concentration for three different sizes of dextran obstacles: (a) D5; (b) D50; and (c) D400. The obtained results exhibit good quantitative agreement with the experimental data except for the smallest dextran size D5. In general, no significant influence of HI in the α exponent is obtained. The lines are to guide the lecturer.
Entropy 19 00105 g009
Figure 10. α-Chymiotrypsin anomalous diffusion exponent α versus dextran concentration for three different sizes of dextran obstacles: (a) D5; (b) D50; and (c) D400. Computations show good quantitative agreement with experiments except for the smallest dextran size D5. The small difference between the results obtained with and without HI reveals HI are not important in the α exponent value. The lines are only for guiding the lecturer.
Figure 10. α-Chymiotrypsin anomalous diffusion exponent α versus dextran concentration for three different sizes of dextran obstacles: (a) D5; (b) D50; and (c) D400. Computations show good quantitative agreement with experiments except for the smallest dextran size D5. The small difference between the results obtained with and without HI reveals HI are not important in the α exponent value. The lines are only for guiding the lecturer.
Entropy 19 00105 g010
Table 1. Best fitted parameters of power law (9) to the experimental radius of gyration and the hydrodynamic radius versus molecular weight.
Table 1. Best fitted parameters of power law (9) to the experimental radius of gyration and the hydrodynamic radius versus molecular weight.
R G R H R c R eff
K (nm·Da γ )0.018 ± 0.0010.043 ± 0.0020.0630.045
γ0.544 ± 0.0050.445 ± 0.004 1 3 0.387
Table 2. Characteristic parameters of the dextran molecules chosen as obstacles. M W : averaged molecular mass in weight; R G : radius of gyration [5,6,7]; R H : hydrodynamic radius [43]; R c : compact spheres radius; R eff : the chosen effective radius. They have been obtained using power law (9) with the parameters shown in Table 1.
Table 2. Characteristic parameters of the dextran molecules chosen as obstacles. M W : averaged molecular mass in weight; R G : radius of gyration [5,6,7]; R H : hydrodynamic radius [43]; R c : compact spheres radius; R eff : the chosen effective radius. They have been obtained using power law (9) with the parameters shown in Table 1.
Dextran M W ( kDa ) R G ( nm ) R H ( nm ) R c ( nm ) R eff ( nm )
D55.21.71.91.11.2
D5048.65.85.22.32.9
D400409.81713.54.76.7

Share and Cite

MDPI and ACS Style

Blanco, P.M.; Via, M.; Garcés, J.L.; Madurga, S.; Mas, F. Brownian Dynamics Computational Model of Protein Diffusion in Crowded Media with Dextran Macromolecules as Obstacles. Entropy 2017, 19, 105. https://doi.org/10.3390/e19030105

AMA Style

Blanco PM, Via M, Garcés JL, Madurga S, Mas F. Brownian Dynamics Computational Model of Protein Diffusion in Crowded Media with Dextran Macromolecules as Obstacles. Entropy. 2017; 19(3):105. https://doi.org/10.3390/e19030105

Chicago/Turabian Style

Blanco, Pablo M., Mireia Via, Josep Lluís Garcés, Sergio Madurga, and Francesc Mas. 2017. "Brownian Dynamics Computational Model of Protein Diffusion in Crowded Media with Dextran Macromolecules as Obstacles" Entropy 19, no. 3: 105. https://doi.org/10.3390/e19030105

Note that from the first issue of 2016, this journal uses article numbers instead of page numbers. See further details here.

Article Metrics

Back to TopTop