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

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^{\mathrm{short}}$) remains constant. Generally, $D^{\mathrm{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^{\mathrm{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>=\left(2d\right)\mathrm{\Gamma }{t}^{\alpha }.$$

(1)

This is formally equivalent to having a diffusion coefficient that is not constant over time [14]:

$$D\left(t\right)=\frac{\mathrm{\Gamma }{t}^{\alpha -1}}{2d},$$

(2)

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]:

$$\mathbf{r}(\mathit{t}+\mathrm{\Delta }\mathit{t})=\mathbf{r}\left(\mathit{t}\right)-\frac{\mathit{D}\mathrm{\Delta }\mathit{t}}{k_B\mathit{T}}\nabla \mathbf{V}(\mathbf{r},\mathit{t})+\sqrt{2\mathit{D}\mathrm{\Delta }\mathit{t}}\mathit{\xi }\left(\mathit{t}\right),$$

(3)

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:

$${\mathbf{V}}_{\mathit{ij}}({\mathbf{r}}_{\mathit{i}},{\mathbf{r}}_{\mathit{j}})=\left\{\begin{array}{cc}\frac{1}{2}\mathit{k}{(d_{\mathit{ij}}-R_{ij})}^2,\hfill & d_{ij}<R_{ij},\hfill \\ 0,\hfill & d_{ij}\ge R_{ij},\hfill \end{array}\right.$$

(4)

where $R_{ij}$ is the sum of the radius of the interacting particles i and j, $d_{ij}$ 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^{\mathrm{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^{\mathrm{short}}$) is found to be a function of the volume fraction $\varphi$, which has been calculated using [5,36,37,38]

$$\varphi =\frac{4\pi (R_{\mathrm{T}}^3+N_{\mathrm{O}}{R}_{\mathrm{O}}^3)}{3L^3},$$

(5)

where $R_{\mathrm{T}}$ is the radius of the tracer protein, and $R_{\mathrm{O}}$ is the radius of the dextran obstacles, $N_{\mathrm{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^{\mathrm{short}}\left(\varphi \right)=\frac{D_0}{[1+H(\varphi \left)\right]},$$

(6)

where $D_0$ is the diffusion coefficient at dilute solution and $H\left(\varphi \right)$ is the contribution due to the HI valid in the short-time regime. $H\left(\varphi \right)$ is given by

$$H\left(\varphi \right)=\frac{2b^2}{1-b}-\frac{c}{1+2c}-\frac{bc(2+c)}{(1+c)(1-b+c)},$$

(7)

where $b=\sqrt{\frac{9}{8}\varphi }$ and $c=\frac{11}{16}\varphi$. Recently, Tokuyama et al. [20] have obtained good agreement for the diffusion coefficient at long times ($D^{\mathrm{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^{\mathrm{long}}$ in terms of the volume fraction:

$$D^{\mathrm{long}}\left(\varphi \right)=\frac{D^{\mathrm{short}}\left(\varphi \right)}{\left[1+\kappa \frac{D^{\mathrm{short}}\left(\varphi \right)}{D_0}\left({\displaystyle \frac{\varphi }{{\varphi }_c}}\right){\left(1-{\displaystyle \frac{\varphi }{{\varphi }_c}}\right)}^{-2}\right]},$$

(8)

where κ and ${\varphi }_c$ are parameters settled to $\kappa =2.0$ and ${\varphi }_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^{\mathrm{short}}$ calculated using Equations (6) and (7).

Tokuyama equations for $D^{\mathrm{short}}$ and $D^{\mathrm{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^{\mathrm{short}}$ involves an approximation. $D^{\mathrm{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_{ij}$ in Equation (4). The resulting $D^{\mathrm{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\times {10}^5$, $k=5\times {10}^4$ and $k=2.5\times {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^{\mathrm{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^{\mathrm{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^{\mathrm{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^{\mathrm{long}}$, we have performed simulations that implement the difference in size by a suitable modification of the minimum inter-particle distance $R_{ij}$ in the interaction potential (Equation (4)). The Tokuyama effective diffusion coefficient $D^{\mathrm{short}}$ is still used to describe the short range HI forces since, to our knowledge, no expression for $D^{\mathrm{short}}$ is available yet for different sized particles.

$D^{\mathrm{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^{\mathrm{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^{\mathrm{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_{\mathrm{H}}$) obtained using quasi-elastic light scattering [43], and radius of gyration ($R_{\mathrm{G}}$) obtained from capillary viscometry [44] are available. We have fitted $R_{\mathrm{H}}$ and $R_{\mathrm{G}}$ versus molecular weight data to the power law:

$$R_x=K·M_w^{\gamma },$$

(9)

where $R_x$ can be either $R_{\mathrm{G}}$ or $R_{\mathrm{H}}$. K and γ are the best fitted parameters reported in

Table 1. Best fitted parameters of power law (9) to the experimental radius of gyration and the hydrodynamic radius versus molecular weight.

	${\mathit{R}}_{\mathbf{G}}$	${\mathit{R}}_{\mathbf{H}}$	${\mathit{R}}_{\mathbf{c}}$	${\mathit{R}}_{\mathbf{eff}}$
K (nm·Da${}^{-\gamma }$)	0.018 ± 0.001	0.043 ± 0.002	0.063	0.045
γ	0.544 ± 0.005	0.445 ± 0.004	$\frac{1}{3}$	0.387

. 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 ${\mathrm{cm}}^3$/g) [47]. Considering the dextran molecules as compact spheres with a radius $R_{\mathrm{c}}$:

$$R_c=\sqrt[3]{\frac{3\nu }{4\pi N_A}}{M}_W^{\frac{1}{3}}=K·M_w^{\gamma },$$

(10)

where $N_A$ is the Avogadro number, $\gamma =\frac{1}{3}$ and $K=\sqrt[3]{\frac{3\nu }{4\pi N_A}}=0.063\mathrm{nm}·{\mathrm{Da}}^{-\frac{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_{\mathrm{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${}^{-\gamma }$ while γ ranged from $\frac{1}{3}$ to 0.445. Among the performed calculations, the best agreement between experimental [5,6,7] and computed $D^{\mathrm{long}}$ was obtained with $\gamma =0.387$ and K = 0.045 nm·Da${}^{-\gamma }$.

It is also worth noting that the chosen parameters for $R_{\mathrm{eff}}$ are closer to those of $R_{\mathrm{c}}$ than to those of $R_{\mathrm{H}}$ (Figure 6). This means that steric compression is crucial to properly describe dextran size in crowded media. Moreover, the large difference between $R_{\mathrm{H}}$ and $R_{\mathrm{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\mathrm{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\mathrm{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. Characteristic parameters of the dextran molecules chosen as obstacles. $M_{\mathrm{W}}$: averaged molecular mass in weight; $R_{\mathrm{G}}$: radius of gyration [5,6,7]; $R_{\mathrm{H}}$: hydrodynamic radius [43]; $R_{\mathrm{c}}$: compact spheres radius; $R_{\mathrm{eff}}$: the chosen effective radius. They have been obtained using power law (9) with the parameters shown in Table 1.

Dextran	${\mathit{M}}_{\mathbf{W}}\left(\mathbf{kDa}\right)$	${\mathit{R}}_{\mathbf{G}}\left(\mathbf{nm}\right)$	${\mathit{R}}_{\mathbf{H}}\left(\mathbf{nm}\right)$	${\mathit{R}}_{\mathbf{c}}\left(\mathbf{nm}\right)$	${\mathit{R}}_{\mathbf{eff}}\left(\mathbf{nm}\right)$
D5	5.2	1.7	1.9	1.1	1.2
D50	48.6	5.8	5.2	2.3	2.9
D400	409.8	17	13.5	4.7	6.7

. 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 ${\tau }_D$. ${\tau }_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 ${\tau }_D$ are used to calculate the effective diffusion coefficient ($D^{\mathrm{eff}}$) at a time equal to ${\tau }_D$. The experimentally found ${\tau }_D$, around 1–3 ms [5,6], is much longer than the computational times at which $D^{\mathrm{long}}$ is calculated. As a consequence, the computed values of $D^{\mathrm{long}}$ should reasonably correspond to the experimentally obtained $D^{\mathrm{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^{\mathrm{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^{\mathrm{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^{\mathrm{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^{\mathrm{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^{\mathrm{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^{\mathrm{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^{\mathrm{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^{\mathrm{long}}$, the α exponent is not significantly influenced by HI. Since the anomalous diffusion exponent is related to the transition rate between $D^{\mathrm{short}}$ and $D^{\mathrm{long}}$ regimes, this means that although HI slow down $D^{\mathrm{short}}$ and $D^{\mathrm{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^{\mathrm{long}}$) are, in general, close to the ones experimentally observed. In the analysis of $D^{\mathrm{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^{\mathrm{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.