Macromolecular diffusion in crowded media beyond the hard-sphere model

The effect of macromolecular crowding in diffusion beyond the hard-core sphere model is studied. A new coarse-grained model is presented, the Chain Entanglement Softened Potential (CESP) model, which takes into account the macromolecular ﬂexibility and chain entanglement. CESP model uses a shoulder-shaped interaction potential which is implemented in the Brownian Dynamics (BD) computations. The interaction potential contains only one parameter associated to the chain entanglement energetic cost ( U r ). The hydrodynamic interactions are included in the BD computations via Tokuyama mean-ﬁeld equations. The model is used to analyze the diffusion of streptavidin protein among different sized dextran obstacles. For this system, U r is obtained by ﬁtting the streptavidin experimental long-time diffusion coefﬁcient D long versus macromolecular concentration for D50 (indicating its molecular weight in kg/mol) dextran obstacles. The obtained D long values show better quantitative agreement to experiments than those obtained with hard-core spheres. Moreover, once parametrized, CESP model is also able to quantitatively predict D long and the anomalous exponent ( α ) for streptavidin diffusion among D10, D400 and D700 dextran obstacles. D long , the short-time diffusion coefﬁcient ( D short ) and α are obtained from the BD simulations by using a new empirical expression, able to describe the full temporal evolution of the diffusion coefﬁcient.


Introduction
Biological media are known to contain a high concentration of a wide variety of macromolecular species such as proteins, polysaccharides or nucleic acids.For instance, the weight fraction of protein is around 5% in lymph, 9% in blood plasma and 35% in hemolysate 1,2 .In the cellular cytosol, three-dimensional analysis of electron micrographs revealed a 20% volume fraction of fibrous supramolecular structures (i.e.F-actin, microtubules and intermediate filaments) 3 .These conditions, known as "macromolecular crowding" involve non-specific interactions among macromolecular species due to excluded volume, van der Waals, electrostatic and hydrodynamic interactions.Macromolecular crowding substantially alters the diffusion processes, conformational properties and reactivity of the macromolecular species [4][5][6][7][8] .
Several efforts are underway to properly understand the effect of crowding in macromolecular diffusion 9 .Recently, in vivo and in vitro experimental studies have been carried out to study macromolecular diffusion in crowded media.In in vivo experiments, fluorescent tracer proteins are introduced into the cell by means of transfection, microinjection or recombinant expression 10 .Although this strategy allows direct study of the macromolecular crowding, the intrinsic differences between the different intracellular microenvironments make the interpretation of the results difficult.In vitro experiments attempt to overcome this problem by reducing the complexity and heterogeneity of the media.The crowded environment is recreated using highly concentrated polymer solutions (usually dextran or ficoll) [11][12][13] .The motion of a fluorescent protein is then studied using spectroscopic techniques 14 , mainly Fluorescent Correlation Spectroscopy (FCS) 11 and Fluorescent Recovery After Photobleaching (FRAP) 12,13 .
In order to interpret the experimental results, an increasing number of computational studies has been performed 15 which provide a highly controlled environment.The comparison of the computational results with the experimentally observed facilitates the quantification of the factors governing macromolecular diffusion in crowded media.Different computational approaches have been applied ranging from on-lattice Monte Carlo simulations 16,17 and off-lattice Brownian Dynamics (BD) [18][19][20][21][22][23][24][25] , to Molecular Dynamics simulations 23,26,27 .In implicit solvent simulations like BD, the Hydrodynamic Interactions (HI) of the macro-molecules have been found to be crucial to properly describe their motion in crowded media 28 .
Macromolecular diffusion is known to become sub-diffusive in crowded media.The time-dependent diffusion coefficient (D) of the macromolecule is defined as D ≡< r 2 > /(2dt) where d is the topological dimension of the system and < r 2 > is the mean square displacement of the particles.In crowded media, D undergoes a transition from an initial value at short times (D short ) to an asymptotic final value at long times (D long ).In the transition between both regimes, macromolecular motion shows a non-linear dependence on time (usually called anomalous diffusion) 16,29 : where α is the anomalous diffusion exponent and Γ α is a generalized transport coefficient.
In order to study macromolecular diffusion in the long time regime, it is necessary to reach µs time scales.However, the large number of atoms present in a highly concentrated macromolecular solution causes atomistic-detailed computations to be limited to smaller time scales.As a consequence, in order to reach larger time scales, coarse grained approaches become necessary.Several approximations have been developed to take into account the HI.The more accurate descriptions are those which involve the calculation of the diffusion tensor 30 , mainly based in the Rotne-Prager-Yamakawa approach [31][32][33] .However, the calculation of the diffusion tensor at every time step is computationally very expensive.An alternative method is based on the mean field equations proposed by Tokuyama [34][35][36] , which introduce an effective diffusion coefficient accounting for the HI.This effective diffusion coefficient is computed as a function of the volume fraction and the dilute solution diffusion coefficient.However, Tokuyama method is limited since it was deduced for systems of equal sized hardcore spheres.In previous studies that use Tokuyama equations, macromolecules are considered as hard spheres and implicit solvent via Brownian Dynamics (BD) is used 18,[20][21][22][23] .
In general, macromolecules have a flexible structure.As a result, when two macromolecules approach each other their branches can become entangled and non-specific attractive and/or repulsive interactions between their chains are expected.This fact can play an important role in macromolecular diffusion in crowded media.Therefore, models going beyond the hardcore spheres, in which the conformational dynamic behaviour of macromolecules is lost, are necessary 37 .Softened interaction potentials have been previously proposed in the study of colloidal systems such as polymeric micelles 38,39 and star-polymers [40][41][42][43][44] .These systems consist of a dense core surrounded by a rather sparse corona.For pluronic polymeric micelles, studies based on implicit solvent and coarse grained computations have shown that both attractive and repulsive forces take place.However, the resulting potential of mean force between two particles is found to be purely repulsive 39 .Similar behaviour is found in the effective inter-particle potentials used in star-polymer models.Star polymers can be modelled as surfaces coated with grafted polymers interacting by excluded volume effects.Particles are large enough to employ the Derjaguin approximation [45][46][47] .The pro-posed potential is the result of the soft repulsion between the polymeric branches and a hard-core repulsion between the rigid cores.
In this paper, diffusion of streptavidin protein among dextran molecules will be analysed.Dextran has a branched, highly hydrated, polymeric structure.They are significantly smaller than micelles and star-polymers, so that some approximations previously used in those systems, such as the Derjaguin approximation, are no longer valid.On the other hand, their very sparse structure and chemical composition, including polar alcoholic groups, can compensate in part the steric repulsions and favour macromolecular entanglement, so that a singular treatment becomes necessary.
In section 2, a new coarse-grained approach is proposed: the Chain Entanglement Softened Potential model (CESP).CESP model is based on an inter-particle interaction potential depicted in Fig. 1.It includes both hard-core and soft interaction regions at proper characteristic lengths.In the soft region, where the potential is shoulder-shaped, particle overlapping is possible, so that macromolecular branches are allowed to become entangled.A single parameter controls the entanglement energy.This potential is introduced in the BD scheme in which the HI are implented using the mean-field Tokuyama equations.D long and α are obtained from the simulated temporal evolution of D by using a new empirical expression.The new equation, unlike the mostly used power law (Eq.1), is able to describe the evolution of D for the three temporal regimes.
In section 3, CESP model is used to study streptavidin diffusion in crowded media among dextrans of different sizes and concentrations.Firstly, the shoulder shaped potential is parametrized by fitting the computed D long to the experimental values 11 for D50 (indicating its molecular weight in kg/mol) dextran obstacles.The BD simulations performed with the CESP model show quantitative agreement between calculated and experimental long time diffusion coefficients (D long ).This agreement is much better than the one obtained when the hard sphere model is used 24 .Once parametrized, we show that CESP model is able to quantitatively predict D long and the anomalous exponent (α) for the rest of dextran sizes (D10, D400 and D700) in a wide range of macromolecular concentrations.Analysis of the radial distribution functions shows a significant increase in the population of entangled macromolecules with macromolecular concentration.This result highlights the importance of including the conformational dynamics of the macromolecules in the diffusion processes in crowded media.

Chain Entanglement Softened Potential (CESP) model
In BD simulations of protein diffusion in crowded media it is usual to reduce the resolution of the macromolecules to effective hard spheres 18,[20][21][22][23] .This approximation significantly reduces the computational cost allowing to reach larger time scales.However, the loss of the macromolecular structural properties could become unrealistic in many cases 37 .As the volume fraction increases, macromolecules are expected to start to become entan-gled, a fact that can not be properly described within the hardcore sphere model 24 .
In the present work, macromolecular entanglement is introduced via a coarse-grained model, the Chain Entanglement Softened Potential model (CESP).The model relies in an empiric inter-particle interaction potential which reminds the Continuous Shouldered Well potential used in computational models of coarse-grained fluids 48,49 .In these studies, density anomalies have been reported when a softened inter-particle potential is applied.The obtained radial distribution functions showed that an inner coordination shell is formed when the fluid density increases.Moreover, the population of the inner coordination shell was found to be proportional to the fluid density.Here, a redefinition of this idea is proposed.The coarse grained particles are no longer fluid molecules but macromolecules in solution.As the density of macromolecules in solution increases, they start to become entangled.Similar potentials have been used in the description of protein-surface interactions 50 .
The interaction potential corresponding to CESP model is shown in Fig. 1.Two clear interaction regions can be observed, which define two characteristic distances: the entanglement distance (d e ) and the core distance (d c ).If the macromolecules are at a distance r larger than d e , they are considered to be separated, and they interact weakly (Fig. 1 C).However, when the macromolecules are at d c ≤ r ≤ d e they are considered to become entangled (Fig. 1 B).In the latter situation, steric repulsions start to arise between the macromolecular chains.Such repulsions can be compensated in part by the presence of van der Waals and hydrogen bonding attraction forces.As a result, a shoulder is developed in the proposed potential (Fig. 1).If the macromolecules get closer (r < d c ), the steric repulsion between the macromolecular skeletons dramatically increases and the hard-core region emerges (Fig. 1 A).The resulting picture resembles models previously proposed for polymer micelles and star-polymers, composed of a dense core and a sparse corona.Unlike these systems, Dextran molecules are polysaccharides with a complex branched structure where no structurally differentiated parts can be identified 51,52 .Therefore (d e ) and (d c ) should be regarded as the characteristic lengths of the different interaction regimes rather than chemically defined parts of the macromolecules.The proposed potential V(r) acting between two particles at a r distance reads where U r is a parameter that quantifies the entanglement energetic cost.The sign of U r indicates which forces, repulsive (positive sign) or attractive (negative sign), dominate in the soft interaction region.The first term in Eq. 2 accounts for the hardcore repulsion region.Following [48][49][50] , the exponent n is settled to n = 24.ε 0 = 1 J/mol and d 0 = 1 nm are just to set the units.V (r) is depicted in Fig. 1 for U r values ranging from -500 to 2000 J/mol.Although the real situation is most probably more complex, Eq. 2 can be regarded as a minimal model, compatible with the physicochemical properties of macromolecules considered, with only one parameter to be fitted.
The proposed interaction potential will be used to analyze the diffusion of streptavidin among different sized dextrans acting as inert obstacles, 11 tracked down using FCS.Both the tracer protein and the dextran obstacles are considered to be able to become entangled.The differences between protein/dextran and dextran/dextran (protein/protein interactions are not relevant due to the large excess of dextran) are then included in an effective way by the specific values of d c and d e for the protein and dextran, while U r value should be understood as an average interaction parameter.In order to estimate d e and d c , the interactions regions in a macromolecule are considered to be concentric spheres with two characteristic radii: the entanglement radius (R e ) and the hard-core radius (R c ), so that for two interacting particles i and j d e,ij = R e,i + R e,j and In this work R e is associated to the hydrodynamic radius of the macromolecule in dilute solution.For dextran molecules, it has been computed by fitting experimental hydrodynamic radii vs molecular weight obtained by quasi-elastic light scattering to an empiric power law function 24,53 .Dilute solution diffusion coefficients (D 0 ) for dextrans are calculated by incorporating the resulting R e values into the Stokes-Einstein equation.On the other hand, HYDROPRO software (Version 10) 54 has been used to compute D 0 corresponding to streptavidin.HYDROPRO considers each atom in the protein surface as an hydrodynamic fric-tional sphere.Protein atomic coordinates corresponding to three different crystallographic structures [55][56][57] were used in the calculations.D 0 was taken as the average of the values resulting from the three structures.R e of streptavidin is then obtained via the Stokes-Einstein equation.
Here we propose to estimate R c as the maximum approach distance between macromolecules and it has been estimated using the specific volume υ, considering them as compact spheres 24 where N A is the Avogadro number and M w is the molecular weight.We have used υ = 0.625 cm 3 /g for dextrans 58 and υ = 0.71 cm 3 /g for streptavidin 59 .M w , D 0 , R e and R c for streptavidin and different sized dextrans are reported in Table 1.Interestingly, the ratio R e /R c increases with dextran size, which could be caused by an increase of dextran branching with molecular weight.Note also that D10 and streptavidin have similar D 0 and R e values.However, D10 has a higher ratio R e /R c than streptavidin, in accordance with the very different chemical structures of both macromolecules, much more compact in the last case.

Brownian Dynamics algorithm
Brownian Dynamics (BD) simulations make use of the Langevin equation in the over-damped limit 60,61 , so that the particles are considered to be under the effect of a stochastic force which accounts for the collisions with the solvent.The stochastic force is mathematically represented by gaussian random noise 19,24,25,62 , and the equations of motion for the N particles read where ξ is a vector of 3N Gaussian random numbers with zero mean and unit variance, T is the system temperature, x is a vector with the 3N Cartesian coordinates, V is the potential of CESP model (Eq.2) and D is the diffusion tensor.
Since the solvent is simulated implicitly, the HI must be included in Eq. 4. HI are solvent mediated correlations of the particle motions and they have been found to play an important role in macromolecular diffusion in crowded media 24,28,63 .HI can be included in Eq. 4 by calculating D following the Rotne-Prager-Yamakawa (RPY) method at each time step 30 , which is computationally very demanding.
Tokuyama equations [34][35][36] offer an alternative to the RPY method based on a mean field approximation.In this approach, an effective diffusion coefficient which accounts for the particle mobility reduction due to the HI is analitically calculated.The expression for the diffusion coefficient is deduced starting from first principles, i.e., the Fokker-Planck equation for the single-particle distribution function coupled with the Navier-Stokes equations.This equation is analytically solved at short times considering a stationary configuration of equal-sized Brownian particles.In doing so, an approximation of the effective diffusion coefficient at short times (D short ) is obtained.For a particle diffusing in a system with a volume fraction φ , D short reads where H(φ ) describes the HI contributions in the short-time regime with b = 9 8 φ and c = 11 16 φ .In this approach D in Eq. 4 is replaced by the effective D short which is a function of φ and D 0 .However, the choice of φ presents two main difficulties.On the one hand, for a number of obstacles N O one could naively calcultate φ as by identifying the obstacle radius R O as the hydrodynamic radius.However, this leads to unphysical φ values larger than unity for the higher dextran concentrations, suggesting the existence of interpenetration between dextran molecules.On the other hand, the soft nature of the particles and the possibility of entanglement makes φ to be undetermined.A possible solution to this problem can be to interprete φ as an effective volume fraction φ eff , depending on the macromolecular concentration and molecular size via some parameters to be determined.Then the problem would be to find the value for these parameters which better predicts the experimental D long using BD dynamics and Tokuyama equations.However, finding the best-fitted value of these parameters would requiere a repeated iteration of the BD dynamics which is computationally very expensive.Instead, a heuristic approach has been chosen, consisting of estimating φ eff using an analytic expression for D long deduced by Tokuyama for equal sized hard-core spheres in absence of the interaction force where κ and φ c are parameters settled to κ = 2.0 and φ c = 1.09 35 .By equating Eq. 8 to the experimental D long we can obtain a first approximation to φ eff .The resulting φ eff values have been found to depend not only on dextran concentration but on the dextran size.For the same dextran mass concentration, larger dextran sizes lead to larger φ eff values, so that dextrans seem to effectively occupy more volume as the molecular size increases.In fact, the computed φ eff values have been found to follow the simple empiric equation (which reminds the Langmuir-Freundlich isotherm 64 ) with z = kυc R e R c , where c is the dextran concentration, υ represents the specific value, and R e /R c the characteristic radii ratio.k and m are fitted parameters whose values are k = 3.7 ± 0.1 and m = 1.43 ± 0.08.The best-fitted φ eff as a function of z for D10 (squares), D50 (circles), D400 (upward-pointing triangles) and D700 (downward-pointing triangles) are shown in Fig. 2. The dashed line represents the best-fitted curve Eq. 9.This procedure actually is expected to provide an approximation to φ eff , whose suitability can only be justified a posteriori, if the BD simulations coupled to CESP model properly reproduces the experimental D long values.This is the case, as will be shown in the next section.
BD simulations are performed in a cubic box with periodic boundary conditions in all directions.The time step is settled to 0.1 ns and the simulation time ranges within 10 to 150 µs depending on the time necessary to reach the stationary value of the diffusion coefficient at long times (D long ).The system temperature is 298.15K. D is computed via Eqs. 5 and 6 using the φ eff predicted by Eq. 9 for each type of particle (dextran and protein).All the studied systems have 100 obstacle dextran particles and only one to five tracer protein particles to ensure that the volume occupied by the tracer less than the 1% of the volume occupied by the obstacles.The length of the simulation box is properly tuned to match the experimental dextran concentration and ranges from 17.56 nm (for the smaller dextran) to 164.32 nm (for the larger dextran).< r 2 > is computed after a thermalisation time of 0.5 µs in order to avoid possible initial unrealistic motion due to the random starting configuration.For studied system, < r 2 > is obtained by averaging from 1000 (higher macromolecular concentration) to 15000 (lower macromolecu-Fig. 2 Effective volume fraction (φ eff ) as a function of z = kυcR e /R c estimated using Eqs. 5 -8 for D10 (squares), D50 (circles), D400 (upward-pointing triangles) and D700 (downward-pointing triangles) dextran molecules.The dashed line represents the fitted Eq. 9.The residuals of the fitting procedure are also depicted and show a very good fitting.lar concentration) BD realizations in order to ensure appropiate statistics.
In Fig. 3, two BD simulation boxes are depicted to illustrate the CESP model.Tracer streptavidin protein molecules (purple) diffuse among D50 dextran obstacles (blue).Each macromolecule is represented as two concentric spheres with radii of R c and R e , which characterize the soft (transparent) and the hard-core (opaque) interaction regions, respectively.The macromolecular concentration of the solutions are 25 g/L (Fig. 3a) and 100 g/L (Fig. 3b).In the less concentrated solution, it can be observed that the soft interaction is enough to avoid macromolecule overlapping.In the more concentrated solutions, steric compression allows the macromolecules to overcome the initial steric repulsion and become entangled.where A, B, β and γ are parameters related to the relevant dynamic properties of the system (the mathematical details are given in the appendix): A = log(D short D long )/2, B = log(D long /D short )/2, β = (α − 1)/B and γ is the inflection time between D short and D long .This new empirical equation is able to describe the diffusion coefficient in the three temporal regimes.This fact represents an advantage compared to the scaling law (Eq.1) which is unable to describe the assimptotic behaviour at low and large times.Moreover, it provides a simple method to compute D short , D long and α from the trajectory of a tracer particle.
In Fig. 4, the time evolution of the simulated diffusion coefficient of streptavidin in a D50 crowded solution, at a macromolecular concentration of 350 g/L, is shown (blue line).As expected, the diffusion coefficient undergoes a transition from D short to D long as macromolecular crowding slows down the motion of the tracer particle.It can be observed that Eq. 10, depicted in red, fits very well to the computed values at all the temporal regimes, allowing rigorous and easy determination of D short , D long and α.

Model parametrization
The proposed shoulder-shaped potential used in the CESP model includes only one parameter to be fitted, i.e., the entanglement energetic cost U r (Eq.2).In order to properly set U r , BD simulations have been performed with U r values ranging from 200 to 2000 J/mol.In these simulations, D long of streptavidin (normalised to D 0 ) has been computed at seven macromolecular concentrations of D50 dextran: 50, 100, 200, 210, 260, 300 and 360 g/L.D long is obtained by fitting the computed log(D) vs log(t) profile to Eq. 10.
The obtained values are shown in Fig. 5 in green colour gradation.As can be observed, D long decreases with U r because of the increasing in the impenetrability of the particles and the consequent reduction in the available space for the particle to diffuse.Experimental D long values reported in 11 are also depicted in the same figure (red markers).For U r = 500 J/mol, very good quantitative agreement between simulations and experiments is obtained.In all the BD computations performed in the rest of this work, U r has been settled to this value, which indicates that, for dextran, the repulsion forces arising from chain entanglement slightly dominate over the attractive ones.D long obtained using the hard sphere model, reported in 24 , are also depicted in Fig. 5 (black markers).In this previous computational work, dextran macromolecules were modeled considering the dextran molecules as hard-core spheres with an effective radius, which was chosen as an intermediate value lying between the hydrodynamic radius 53 and the compact radius.This approach accounted for the possibility of macromolecular overlapping on an average way, but represented a static description of the conformational state, being unable to properly mimic macromolecular entanglement.As it is clearly observed in Fig. 5, D long values obtained using the CESP model are much closer to the experimental ones than those coming from the hard-core sphere model 24 .For the rest of dextran sizes (D10, D400 and D700), direct comparison of D long predicted by the hard-sphere model 24 and CESP model can be found in the Suplementary Information.Again, it is found that CESP model significantly improves the predictions of the experimental D long obtained using the hard-core sphere model.Therefore, the inclusion of coupling of the dynamic and conformational properties seems to be necessary to accurately describe diffusion in crowded media.J/mol (empty green inverted triangles), 1000 J/mol (green inverted triangles) and 2000 J/mol (green diamonds).Black square markers represent the values obtained using the hard sphere model (taken from Ref. 24 ), while red circle markers denote the experimental D long value (taken from Ref. 11 ).For U r =500 J/mol, good agreement between CESP model and experimental values is observed.Lines are only to guide the eyes.

Model prediction for D long and the anomalous exponent α
In the previous section, we have found that for U r = 500 J/mol, CESP model is able to describe D long of streptavidin for D50 dextran obstacles.Let us now see that this model can also predict, without any further parametrization, D long and α of streptavidin for different dextran sizes (D10, D400 and D700) and concentrations (from 25 to 350 g/L).D long and α are again obtained from the BD simulations by fitting Eq. 10 to the computed D time evolution.The curve log(D) versus log(t) has been plotted in Fig. 6, together with the best fitted Eq. 10, for four dextran sizes and concentrations: (Fig. 6a) D10 at 300 g/L, (Fig. 6b) D50 at 200 g/L, (Fig. 6c) D400 at 100 g/L and (Fig. 6d) D700 at 50 g/L.As can be observed, for all the systems presented, Eq. 10 fits very well to the computational results.
The computed D long versus macromolecular concentration are depicted in Fig. 7 for four dextran sizes: D10, D50, D400 and D700.The results are shown together with the experimental values reported in 11 for the same systems.Dashed lines are only to guide the eyes and they follow the simulation results.In general, the computed results are in very good quantitative agreement with the experimental values.However, the values corresponding to D10 exhibit some systematic error, which could be explained taking into account that the potential energy barrier U r has been set, regardless of the dextran size, to the best fit D long for D50.This suggests that U r should be smaller for D10 than for D50, D400 and D700.Since U r accounts for the energy cost of the macromolecular entanglement, this fact could suggest a lower energy cost for D10 entanglement than the one for D50, D400 and D700.In agreement with the experiments reported in 11 , BD simulations predicts a decrease in D long with macromolecular concentration, as a result of volume exclusion and hydrodynamic interactions.Alternative potentials based on a self-consistent field method appropriate for weak excluded interactions 45 have been implemented in BD simulations of star-polymer 42 .As in dextran molecules, the computational results exhibit similar D long decay with the particle density.However, at the quantitative level, the calculated values are much larger than the experimental ones.The discrepancies could be explained by the fact that HI were not included in the computations of ref. 42 .In the present work HI are included in the BD scheme, and it has been found that they are crucial for the quantitative reproduction of the experimental long values 24,28 .
The anomalous exponent α provides information about the transition between the short and the long time regimes.α values close to unity indicate a smooth transition.As α deviates  from unity, the transition becomes more pronounced.The simulated α versus macromolecular concentration are depicted in Fig. 8 together with the experimental values obtained using FCS 11 .As can be observed, the CESP model qualitatively captures the experimentally observed behaviour.α decays with macromolecular concentration as a result of an increase of particle collisions.Moreover, α also decrease with dextran size which could be a consequence of the increase of φ eff , so that the volume fraction of the system seems to play a crucial role in and therefore in the transition between the D short and D long steady states.However, quantitative discrepancies are observed, which could be explained not only by the limitations of the model, but to the way in which α has been calculated.The experimental α values have been obtained by fitting the auto-correlation function of the tracer protein fluorescence using the power law described in Eq. 1, while here α has been determined as the slope in the inflection point of the curve log(D) versus log(t) using Eq.10.

Radial distribution functions
The simulated radial distribution functions (rdf) for the cases analyzed in the previous section (D10 (a), D50 (b), D400 (c) and D700 (d)) are shown in Fig. 9 for macromolecular concentrations ranging from 25 g/L to 350 g/L.The rdfs have been computed from the BD computations by calculating the center-to-center distance of every dextran obstacle to the rest of dextran particles.Averages have been taken over 10000 configurations every 2 ns of simulation.As can be observed, CESP model predicts an increase in the number of entangled macromolecules with macromolecular concentration, as revealed by the gradual increase of the rdf at r = d c .This fact is observed for all the studied dextran sizes.Two sharp peaks at distances slightly lower than d c and d e define two coordination shells of entangled macromolecules, corresponding to highly entangled macromolecules (close to d c ) and weakly entangled macromolecules (close to d e ), respectively.It is also interesting to note the presence of two secondary peaks for D400 and D700 at macromolecular concentrations of 250, 300 and 350 g/L, corresponding to additional coordination shells and suggesting other possible degrees of entanglement for the biggest dextran molecules.Note that the two main peaks in the entanglement region are developed even in the absence of explicit attractive terms in the proposed potential.Entanglement is thus the way the system has to minimize the excluded volume effects even paying a small energetic cost.It is reasonable to conclude that the presence of attractive forces would enhance this effect, although it is not necessary to produce it.
The first rdf peak has been integrated in order to determine the number of entangled particles (N E ) close to the core of the reference particle.N E is shown in Fig. 10 for the four analyzed dextran sizes.It can be observed that N E increases exponentially with macromolecular concentration, so that entanglement is promoted by steric compression.This effect is enhanced by the increase of dextran size, as a consequence of the increase in φ eff , so that bigger dextran obstacles are able to effectively occupy higher volumes than the smaller ones for the same macromolecular con-

Conclusions
A new coarse grained approach to macromolecular diffusion in crowded media which goes beyond the hard sphere model is presented: the Chain Entanglement Softened Potential (CESP) model.Two interaction regions between macromolecules are considered: a soft interaction region caused by chain entanglement, a combination of repulsive and attractive effects, and a hard-core interaction region caused by the dramatic increase of steric repulsions between macromolecular skeletons when they approach bellow a certain distance.The resulting picture is quantified by means of a shoulder-shaped potential which allows to take into account macromolecular conformational flexibility and chain entanglement.The model can be parametrized by a unique parameter (U r ) which is associated to the entanglement energetic cost.For streptavidin protein diffusion in D50 crowded solution, U r has been quantified by fitting the long time diffusion coefficient (D long ) resulting from BD simulations to the experimental values obtained using Flourescent Correlation Spectroscopy (FCS) 11 .
In order to determine D long and the anomalous exponent α from BD simulations, a new empiric expression for the temporal evolution of the diffusion coefficient, valid for all the temporal regimes, has been proposed.The new function goes beyond the standard power law and depends on the four relevant dynamic properties: D short , D long , α and the inflection time γ.This function fits very well to the BD computations for all the obstacle concentrations and sizes under study, and could be useful in forthcoming computational studies focused on anomalous diffusion.
The computed D long using the CESP model clearly improve the results obtained in previous studies based on hard-core spheres 24 .Without any further parametrization of V (r), CESP model is able to provide quantitative prediction to D long and α of streptavidin in media crowded by D10, D400 and D700 at several dextran concentrations (25 to 350 g/L).Radial distribution functions exhibit a significative increase of the population of entangled macromolecules with macromolecular concentration and dextran size.
Our results indicate that an adequate description of macromolec-ular entanglement is essential to understand macromolecular diffusion in crowded media.

Fig. 1
Fig.1Outline of the Chain Entanglement Softened Potential (CESP) model and the potential energy V (r ) vs inter-particle distance r (Eq.2).Three interaction regions can be observed: no interaction (C), soft interaction where entanglement is allowed (B), and hard-core interaction (A).The chosen parameters are d c = 4.58 nm and d e = 10.48 nm, for U r ranging from -500 to 2000 J/mol.

Fig. 3
Fig. 3 Illustration of the simulation box of two BD simulations of streptavidin protein (purple) diffusing among D50 dextran obstacles (blue).Dextran concentrations are (a) 25 g/L and (b) 100 g/L.Each macromolecule is represented as two concentric spheres representing the soft (transparent) and hard-core (opaque) interaction regions.As concentration increases, steric compression promotes macromolecules to become entangled.

Fig. 4
Fig. 4 Time evolution of the diffusion coefficient (D) in logarithmic scale.The blue line results from averaging 3000 BD simulations of 100 µs long in a cubic simulation box with a side of 28.466 nm.The simulation box contains 5 tracer particles (streptavidin) diffusing among 100 obstacle particles (D50 dextran).The macromolecular concentration is 350 g/L.Macromolecular crowding slows down the tracer diffusion coefficient from D short to D long .The red line is the result of fitting Eq. 10 to the simulated which allows to easily determine D short , D long and α .The resulting fitting is very good, as shown by the residuals.

Fig. 5
Fig.5 Calculated D long /D 0 as a function of the concentration corresponding to D50 dextran for U r values: 200 J/mol (green triangles), 500 J/mol (empty green inverted triangles), 1000 J/mol (green inverted triangles) and 2000 J/mol (green diamonds).Black square markers represent the values obtained using the hard sphere model (taken from Ref.24 ), while red circle markers denote the experimental D long value (taken from Ref.11 ).For U r =500 J/mol, good agreement between CESP model and experimental values is observed.Lines are only to guide the eyes.

Fig. 6
Fig. 6 Time evolution of the diffusion coefficient D versus time in logarithmic scale obtained from BD computations (red line).The best-fitted Eq. 10 is plotted as a black dashed line.Four cases are depicted: (a) D10 at 300 g/L, (b) D50 at 200 g/L, (c) D400 at 100 g/L and (d) D700 at 50 g/L.For all the systems presented, Eq. 10 fits very well to the computed values, allowing easy determination of D long , D short and α from the temporal evolution of D.

Fig. 7
Fig. 7 log (D long /D 0 ) for streptavidin versus macromolecular concentration.The filled markers represent BD simulations and the empty ones correspond to experimental values 11 for four dextran sizes: D10 (green squares), D50 (red circles), D400 (upward pointing blue triangles) and D700 (downward pointing purple triangles).Dashed lines are only to guide the eyes and they follow the simulation results.

Fig. 8
Fig. 8 Anomalous exponent α versus macromolecular concentration.Filled markers represent BD simulations and empty markers correspond to the experimental values 11 for four dextran sizes: D10 (green squares), D50 (red circles), D400 (upward pointing blue triangles) and D700 (downward pointing purple triangles).Dashed lines are only to guide the eyes and they follow the simulation results.

Fig. 9
Fig. 9 Simulated radial distribution functions (rdf) versus inter-particle distance using the CESP model for four dextran sizes: (a) D10, (b) D50, (c) D400 and (d) D700.For each size, the rdf has been computed at eight concentrations ranging from 25 to 350 g/L, depicted in color gradation from yellow (25 g/L) to purple (350 g/L).

Fig. 10
Fig. 10 Number of entangled macromolecules around a reference particle (N E ) versus macromolecular concentration for four dextran sizes: D10 (green squares), D50 (red circles), D400 (blue upward pointing triangles) and D700 (purple downward pointing triangles).The number of entangled particles increases with macromolecular concentration as a result of steric compression.