Fluid flow and amyloid transport and aggregation in the brain interstitial space

Abstract

The driving mechanisms at the base of the clearance of biological wastes in the brain interstitial space (ISS) are still poorly understood and an actively debated subject. A complete comprehension of the processes that lead to the aggregation of amyloid proteins in such environment, hallmark of the onset and progression of Alzheimer’s disease, is of crucial relevance. Here we employ combined computational fluid dynamics and molecular dynamics techniques to uncover the role of fluid flow and proteins transport in the brain ISS. Our work identifies diffusion as the principal mechanism for amyloid-β proteins clearance, whereas fluid advection may lead transport for larger molecular bodies, like amyloid-β aggregates or extracellular vesicles. We also clearly quantify the impact of large nascent prefibrils on the fluid flowing and shearing. Finally, we show that, even in the irregular brain interstitial space (ISS), hydrodynamic interactions enhance amyloid-β aggregation at all stages of the aggregation pathway. Our results are key to understand the role of fluid flow and solvent-solute interplay on therapeutics like antibodies acting in the brain ISS.

Introduction

Transport of substances in the brain is fundamental for its healthy functioning. This includes the income of nutrients, the effective clearance of metabolic byproducts, as well as the efficient deploy of therapeutic drugs.

A key role in that is played by the brain ISS, the fluid-filled space that surrounds brain’s cells (1, 2). All the components of this environment, from the fluid component to the extracellular matrix scaffold to its geometrical peculiarities, are of paramount importance for the correct functioning of the brain’s tissues (1, 2). The role of the physiological flow of the brain interstital fluid (ISF) in the clearance mechanism of biological wastes is still poorly understood. A malfunctioning clearance may lead to the accumulation of toxic substances in the ISS between cells, leading to severe consequences. It is widely accepted that the accumulation of amyloid fibrils and oligomers of the Aβ protein is linked to the onset and progression of the Alzheimer’s disease (AD) (3).

Motion of the ISF in the brain follows intricate paths (4, 5). Cerebrospinal fluid can enter and leave the brain parenchyma by flowing in the perivascular space, where it mixes with the ISF favoring its draining. The driving mechanisms for this clearance process are vasomotion (6) and arterial pulsation (7–9), and factors like sleep (10–12), age (13) and pathologies such as AD (14), can enhance or, on the contrary, impair its efficiency. Principal clearance mechanisms in the ISS, on the contrary, are much debated, with evidences in favor both of an active transport by fluid flow, as stated the glymphatic hypothesis (15, 16), or by a diffusive process as driving mechanism (17–19).

Aggregation of proteins has been identified as one of the main cause of many neurodegenerative diseases (20). Among them, aggregation of Aβ proteins, the hallmark of the progress of AD, has been deeply investigated (3). Aβs aggregation takes place in the ISS, in a crowded and dynamic environment (2, 21). We now know that factors like crowding and fluid motion influence the aggregation dynamics and the structural properties of the resulting aggregates (22–27). However, these effects have been investigated only under idealized in vitro conditions, in sharp contrast with the realistic in vivo environment.

The study of the ISS and its properties is a very active research field (2, 21, 28–30) fostered by new technological advances. Recent improvement in imaging technique made readily available highly detailed images of brain tissues (31, 32). Resulting three-dimensional reconstructions of the ISS, in combinations with standard computational fluid dynamics techniques (19, 33) and machine learning (34), enabled the investigation of fluid motion in various structures of the brain.

Here we deploy multiphysics and multiscale simulations based on combining fluid (35) and particle dynamics (36, 37), and an accurate three-dimensional representation of a portion of the ISS, to characterize the ISF fluid flow and, employing for the first time a microscopic approach to study particles motion, its contribution to the transport of Aβ proteins. We reproduce the experimentally estimated average physiological velocity, and quantify the local shear rate. We clearly demonstrate that the transport of the individual Aβ is diffusion limited, while advection may contribute to the clearance of larger molecular bodies. Finally, we analyze the essential contribution of hydrodynamic interactions on proteins aggregation in the ISS.

Results

Building the geometry for Lattice Boltzmann simulations

To build a realistic three-dimensional representation of a small portion of ISS, we used the publicly available data of a reconstructed hippocampal neuropil of a brain’s rat (31).

The tissue sample, extending over 4 μm, presents several disconnected components, depicted with different colors in the three-dimensional representation in Fig. 1a, each one representing a different cell or part of it. The intricate structure of different cells is interspaced by the ISS, forming a network of sheets and tunnels of characteristic width of /37 nm and 41 nm, respectively (31). A cubic sample of side 1 μm was extracted and used as a starting point to generate the discrete representation of the ISS used in our simulations. The number of triangular meshes composing the surfaces of the cells was first reduced by decimation algorithm tuned to preserve the topology of the surfaces. The resulting geometry was subsequently divided in cubic voxels with a spatial resolution of $Δx=5$ nm, to sufficiently resolve the sheets and tunnels of the ISS network. This step provides a discretized representation of the cells’ tissue, from which, using a simple inversion, we were able to retrieve the lattice spanning the ISS. Finally, to avoid finite size and confinement artifacts, the lattice was periodicized along the x- and y-directions, see Fig. 1a. The lattice was left nonperiodic along the z axis, as we imposed specific boundary conditions to accurately reproduce the physiological fluid flow (see next Section). The final discretized geometry is of dimensions $2μm×2μm×1μm$⁠, and is characterized by a ratio between ISS and whole tissue volumes, α, of \0.21, close to the estimate in healthy physiological conditions (21, 38).

The ISS is an highly dynamic environment; its volume can undergo slow periodic physiological variations of up to 60% during the circadian cycle (39) or even abrupt fluctuations due to pathological conditions like epileptic seizure (32). Long-term conditions can also permanently impact brain’s tissues, with an increase of the ISS volume observed in mice with Parkinson’s disease (40) or Alzheimer’s disease (41). To account for this variability, two additional expanded ISS geometries with $α=0.30$ and $α=0.35$ were built starting from the original one.

Fluid velocity and shear rate values

We first performed fluid-dynamics simulation of the fluid flow via the Lattice Boltzmann (LB) method. To reproduce a physiological fluid flow along z, we imposed a velocity inflow at the top and outflow at the bottom of the system. This condition was enforced with a velocity boundary condition at $z=1μm$ (top of the simulation box), with a value of the velocity fixed at $3.32⋅10−7$ m/s, and a pressure boundary condition at $z=0μm$ (bottom of the simulation box), with the pressure fixed at 0.00 Pa. Periodic boundary conditions were used along the x- and y-directions. The selected value of the velocity at the inlet, allows to set the fluid velocity in the correct order of magnitude comparable with experimental estimates (see Fig. 1e and comments below). The choice of a unidirectional fluid forcing is consistent with the observation that what drives the fluid motion in the ISS is the pressure gradient along the direction from cerebral arterioles to venules. Given the typical arteriole-venule distance of \200–300 μm (42), at the spatial scale we are simulating a region of the ISS, we can safely assume that the external forcing that drives the fluid motion is present only along a single direction. Analogous assumptions are widely used in similar fluid dynamics study of the ISF motion in small portions of the brain tissue (19, 33). The physiological fluid flows is expected not to significantly alter the geometry of the cells’ surface that is modeled as impermeable and rigid. Consequently, no-slip boundary conditions were enforced on external surface of the cells. An estimate of the Reynolds number yields $Re≈10−9≪1$⁠, characteristic of a Stokes flow, with viscous forces dominating over inertia.

Results from the LB simulations are reported in Fig. 1b–e (see also Supplementary Material, Movie S1). The inlet velocity pushes the ISF to flow mostly along paths that traverse the ISS geometry along the z axis direction, where the magnitude of the velocity field is higher. In the orthogonal regions, the fluid is almost still and flows at much lower speed. A map of the velocity field on a slice parallel to the xz plane allows to visualize its distribution (Fig. 1c). Values range from 0.0 m/s to $8.0⋅10−7$ m/s, with the higher values distributed along straight vertical sheets or wider space where several sheets join. As expected, given the no-slip condition on the cells’ surface, the fluid slows down close to the cells and attains the maximum velocity in the middle of the sheets and tunnels. The distribution of the fluid’s velocity magnitude computed on the whole geometry (Fig. 1e), shows a long but low tail that extends toward values of few μm/s, while the bulk of the distribution ranges by tenths of μm/s. Increasing the volume fraction of the ISS, as to mimic the day/night variation, the tails become slightly shorter and the overall distributions taller. For comparison, a recent study that employed the same tissue sample with a different computational technique and boundary conditions (19), found a range of velocities that is compatible with the one reported here, confirming the reliability of our results. Average values of the fluid velocity magnitude (Fig. 1e, left panel) are in the $2.4–3.0⋅10−7$ m/s range, close to the experimental values found in different species (1). It is interesting to note that, albeit slightly higher (approximately by a factor 1.2), the average velocities of the fluid in the expanded ISS geometry, $α=0.30$ and 0.35, remain of the same order of magnitude of the one in physiological conditions. The shear rate is defined as the variation of a velocity field along a given direction. In our three-dimensional case, it is more appropriate to use the strain-rate tensor and its second-invariant, a scalar $γ˙$ defined on each point of the fluid’s lattice (see Materials and methods). A first estimate of the order of magnitude of the shear rate $γ˙$ yields a value of $γ˙≈10s−1$⁠. Distributions of the shear rate $γ˙$ reveal values that span three order of magnitude, from fractions of s⁻¹ to hundreds of s⁻¹ (Fig. 1e, right panel). The highest values are located close to the cells’ surface or at the entrance of a tunnel, where fluid flows from different directions due to the convergence of separate sheets (Fig. 1d). Increasing the ISS to higher α, the distributions present shorter tails and a monotonously increasing peak. Thus, even if the fluid velocity is larger, its spatial variation appears smoother. Average values present a decreasing trend as a function of the volume fraction α, from 20 s⁻¹ to 13 s⁻¹, close to our initial estimate.

Fibril in the brain ISS

We next investigate the alterations of the fluid flow due to the presence of large molecular bodies modeled to represent amyloid fibrils. A fibril is idealized as an elongated rigid structure inserted in the ISS so to divert the normal motion of the fluid. The fibril is not able to move and the fluid in contact with its surfaces is imposed to be still (no-slip condition). An example of an 87 nm long model-fibril embedded in the ISS is depicted in red in Fig. 2a. The magnitude of the fluid velocity field on a plane parallel to the coronal plane of the fibril and the resulting shear rate field, show marked differences in the absence and presence of the fibril (Fig. 2b and c, respectively). The streamlines representing the path followed by the fluid flow are strongly altered by the presence of the fibril, as now it has to move around the obstacles and the available space is greatly reduced. The surface contact with nearby structure is also increased and this further alters the value of the velocity field. Accordingly, the shear rate field is affected by the narrower available space and fluid–surface interactions, which leads to the formation of new regions of space with noticeable variations in the fluid velocity. The comparison of the resulting density distribution functions (Fig. 2d) reveals that, overall, the fluid close to the fibril flows at lower speed, with a noticeable smaller right tail and a slightly higher populations at intermediate values. The distributions pertaining to the shear rate values almost perfectly overlap. The fibril effectively diverts the fluid paths due to its excluded volume, without in any case altering to a great extent the properties of the fluid flow. Due to the very slow velocities at play, even if the space available to the fluid is diminished, we estimate that the Reynolds number remains very small, well within the Stokes regime.

Transport of amyloid peptides

We now focus on particles motion in the ISS, namely the transport properties of $Aβ40$ proteins. The hydrodynamic radius of an $Aβ40$ is of 1.8 nm (43, 44), thus at our resolution an amyloid protein is represented as a single point-like object that freely moves in the ISS, constrained in the geometry by a mostly repulsive Lennard–Jones interactions with the frozen lattice points composing the surfaces of the cells. We performed Langevin dynamics simulations using a friction value γ tuned to reproduce the experimental diffusivity of $Aβ40$ in dilute conditions.

Figure 3 panel (a) shows the top view and a side view of the geometry along with the trajectories of five particles representing each one an amyloid protein. Due to the complex conformation of the geometry, the motions of the dispersed particles is hindered and temporarily confined in delimited regions of space. We first explore how geometrical confinement affects the motion of particles and to what extend the pushing effect of a fluid flowing in the ISS alters this motion. Two sets of simulations where ran, one with particles only subject to thermal noise and the second one where, on top of the thermal motion, an external forcing term $f→fluid$⁠, mimicking the pushing force exerted by the ISF on the particle, was added. The forcing term is of the form $f→=γmu→i$⁠, where γ is the friction defined in the Langevin equation having the dimension of the inverse of a time, m is the mass of the particle and $u→i$ is the velocity of the fluid field at the closest fluid lattice point i extracted from the particle-free LB simulations. This setup is a one-way coupling between the fluid and the moving particle: the motion of the particle is perturbed by the pushing action of the surrounding fluid but it does not perturb the fluid itself. Fig. 3b reports the temporal and ensemble average MSD $⟨MSD(Δt)¯⟩$ of the particles in the ISS with and without the pushing fluid flow. All the MSDs show a linear dependence from the time lag $Δt$⁠, typical of a diffusive regime. The MSDs resulting from the motion of particles in the ISS noticeably deviate from the free diffusion case, showing a slower diffusive behavior due to the confining geometry. The value of the fluid forcing is clearly too small to induce an active transport that would cause a deviation from the linear regime. MSDs were fitted against a functional form linear in $Δt$⁠, $6DΔt$⁠, in the [0.4, 0.8] ms range to extract the diffusion coefficient D and yielding very similar values, $D=8.7±0.3⋅10−11m2$/s and $Dfluid=8.4±0.5⋅10−11m2$/s for the case without and with the fluid forcing, respectively. A useful parameter to characterize the change on the diffusive properties of a particle with respect to the free diffusion case, is the tortuosity parameter $λ=D0D$⁠, with $D0$ the reference free diffusion coefficient. The tortuosity λ embeds all the factor that may hinder a particle motions, namely geometrical constraints or interactions with others components of the solutions. From experiments on brain’s tissues, λ has been estimated to be in the range 1.4–2.0 (1, 21, 38), with an important variability depending on the experimental technique and intrinsic properties of the brain’s sample. From the simulations, we obtained $λ=1.55±0.03$ and $λ=1.58±0.05$ without and with the action of the fluid respectively, values that perfectly fall in the experimental range. The confining effect of the geometry alone slows down the particles diffusion by a factor 2.4–2.5. Also the distributions of particles in the ISS geometry is not significantly altered by the presence of the forcing fluid (see Figs. S6 and S7). To further refine the analysis, we implemented a simple effective model of the extracellular matrix. Extracellular matrix in the ISS is mainly composed of hyaluronic acid molecules connected to form a gel that acts as a scaffold, contributes to the normal physiology of the brain but also lower molecular diffusion in the ISS (21). Crowding in cellular and extracellular environments can go up to 30% in excluded volume and affects, for example, solute transport properties (45–47) and amyloid aggregation (22, 23). To model the excluded volume of the extracellular matrix, we represented the latter as spherical quenched repulsive particles randomly placed in the ISS for values of the excluded volume up to 30%. The resulting tortuosity λ as a function of the excluded volume percentage is reported in Fig. 3c (see Fig. S8 for the respective MSDs). The tortuosity monotonously increases as a function of the degree of crowding, further improving our estimate of the tortuosity against experimental results. This confirms the key role played by the extracellular matrix in the transport of nutrient and biological wastes in the ISS and the necessity to accurately include its effects.

Finally, to estimate the values at which the transport of amyloid peptides transition from diffusion-limited to advection-dominated, we artificially increased the intensity of the fluid motion. Five sets of simulations were performed, each one with a differently tuned magnitude of the fluid velocity field computed in physiological conditions, from no velocity at all to a magnitude increased by a factor $×104$⁠. In analogy with the experimental practice (2), a point source was selected to release the peptides and exit times from a spherically bounded region of space and peptides’ MSDs were computed (see Fig. 3d for a pictorial representation). Average exit times increase for increasing radius R of the spherical region (Fig. 3e), with an initial departure from the “diffusive” zero-velocity case at values >150 nm when the magnitude of the velocity field acting on the peptides is multiplied by a $103$ factor. A marked deviation is instead noticeable also at smaller values of R when the multiplicative factor is $104$⁠. The MSDs further confirm this picture, with a strong deviation from the diffusive regime already after $Δt≈$ 0.03 ms for the highest value of the fluid velocity field. A fit with the functional form $6DΔt+v2Δt2$ allows to extract the effective drift velocity $v=1.3⋅10−3$ m/s, of the same order of magnitude of the average value of the enhanced fluid velocity field. Using the drift velocity v, the diffusion coefficient $D=1.7⋅10−10m2$/s and $L=1μ$m as characteristic lenghtscale for the transport process in our geometry, the Péclet number results $Pe=LvD=8.2>1$⁠, showing that, in this forcing regime, advection dominates the amyloids motion at this lengthscale.

Aggregation of amyloid peptides

We now change scale, and investigate the $Aβ40$ aggregation dynamics in a restricted zone of the ISS at higher resolution. We employ the Lattice Boltzmann Molecular Dynamic technique to fully couple the fluid motion with the motion of solute particles (37) in order to explicitly account for hydrodynamic interactions (see Materials and methods). We performed simulations to study the early stage aggregation of amyloid proteins in a spatially inhomogeneous geometry with the additional presence of a large elongated aggregate to model the fibril.

A smaller portion of the ISS geometry was discretized with a lattice resolution of $Δx=2$ nm and subsequently periodicizied using the same strategy employed for the whole system. The final cubic system is of side 192 nm (Fig. 4a–c). The fibril is composed by an elastic network of 1,125 proteins modeled as particles and coupled to the fluid motion. Dispersed amyloid proteins, in a concentration of 0.1 mM, are modeled as particles interacting directly via a Lennard–Jones potential. To avoid high energy collision and numerical instabilities, the dispersed particles interact with particles forming the fibril solely via hydrodynamics interactions. All the particles in the system are repelled by the walls of the geometry formed by cells’ membrane. Periodic boundary conditions are enforced along the three coordinates x, y, and z.

We first consider the aggregation process of amyloid particles in the irregular and inhomogeneous space (Fig. 4d). For this purpose, we clustered the proteins trajectories using a density-based algorithm (48). During their motion, single peptides start to bind and form dimers and, later on, trimers. The decrease of unbind proteins over time and the consequent increase of the number of clusters is strongly accelerated by Hydrodynamic Interactions (HI) (Fig. 4e and f, and Supplementary Material, Movie S2) in accordance with previous simulations of colloidal particles (49) and coarse-grained representation of amyloid peptides (50). The average fractions of dimers and trimers show a marked difference in population due to HI, with the fraction of trimers about five times higher when a full coupling between fluid and particles is allowed (Fig. 4g). In the presence of HI, the population of dimers at a given time is higher than in the absence of fluid-mediated interactions, and it is thus more likely that a monomer and a dimer will cluster together forming a trimer. In addition to that, according to the Stokes–Einstein relation, the diffusion coefficient of an aggregate is inversely proportional to the radius of gyration, $D∝Rg−1$⁠, thus reducing its thermal motion and, at the same time, long-range HI starts to play a major role in the aggregation rate of aggregates of different size (51).

We then embedded a fibrillar body in the irregular space (Fig. 4h). The large aggregate, composed of many bonded Aβ proteins that move due to thermal motion, has a strong impact on the surrounding environment, with fluid perturbations that emanate from the fibril and spread across the system (Fig. 4i and l, and Supplementary Material, Movies S3 and S4). The thermal vibration of proteins push on fluid and, as a result of the momentum conservation, the latter acquires a velocity that is order of magnitude higher than the background fluid flow caused by the physiological pressure gradient. The irregular geometry, with internal walls and reentrant spaces, further contributes to the complex interplay between the fluid environment, the soluted particles and the space into which these two components are embedded.

Finally, we consider how the presence of the fibril impacts the motion of surrounding molecules (Fig. 4k). The binding of a single monomer to a prexistent prefibril is part of the complex pathway that leads to fibril growth. By tracking the position of each single protein with respect to the agglomerate, we were able to determine the number of $Aβ40$ that get close to the fibril (Fig. 4n). From simulations, we found that the number of proteins that moved in proximity of the fibril is four to five times larger in the presence of HI than in its absence. This result hints to the importance that HI have, not only in the early aggregation steps, but also along the whole formation process of large amyloid fibrils (52). At the same time, the additional fluid perturbations due to the fibril do not significantly affect neither the speed at which dispersed peptides cluster together (Fig. 4f and o), nor the relative population of dimers and trimers (Fig. 4g and p).

Discussion

Transport of molecules in the ISS, being them nutrients, biological wastes of neuromodulator, is still a much debated problem. In this paper, we used a sophisticated computational approach that combines fluid dynamics and particle-based simulations, together with an accurate three-dimensional representation of the geometry of ISS, to characterize the properties of ISF and the transport and aggregation processes of $Aβ$ proteins in this space.

Our simulations reproduce the experimental estimate of the average velocity of the bulk fluid flow, $∼10−7$ m/s (1). Despite the wider space available for the fluid during physiological expansion of the ISS, as during sleep (39), from our simulations the distribution and average value of its velocity magnitude are not significantly affected. The ISF bulk flow in the ISS is driven by a pressure gradient present at the extrema of the neuropil. It is experimentally very arduous to provide a correct estimate of the pressure gradient in realistic conditions, and a fairly high upper bound is thought to be of 1 mmHg/mm. In this work, we imposed mixed boundary conditions to determine a value of the fluid flow $10−7$ m/s. This corresponds to a pressure drop of 27.8 mmHg/mm. Given the linear proportionality between pressure gradient and fluid velocity, we may infer that, with a smaller pressure drop and the sample of the ISS used in this work, the resulting bulk flow velocity would have been even smaller that the one reported here (19). The ISF flow velocity can be compared to that of the cerebrospinal fluid flow in the perivascular space driven by arterial pulsations and vasomotion (6, 8). In this case, the attained fluid’s velocities are $10−5$ m/s (6, 8, 34), two order of magnitude higher than that of the ISF, and advection is supposedly the main driving mechanism for perivascular clearance.

The shear rate maximal values found in this work are of the order of 10–100 s⁻¹, with a weak dependence on the volume ratio α. The resulting average values, of the order of tens of s⁻¹, are consistent with those recently found using artificial intelligence to quantify the flow properties of the cerebrospinal fluid in the perivascular space (34), albeit smaller given the lower fluid velocity. While the shear values in the ISS are too small to cause the complete unfolding of proteins, or to enhance their aggregation kinetic, they might be relevant in the context of AD therapy based on antibodies. In fact, conformational changes and partial unfolding of antibodies have been reported experimentally for relatively moderate values of shear rates, with consequential loss of affinity (53). Moreover, it has been shown that repeated expositions to extensional flow in a microfluidic device, a process analogous to traveling in the intricate ISS, positively correlates with the aggregation of proteins, thus potentially decreasing their biological or pharmaceutical functionality (54). Our results are also relevant with respect to the promising Convection-Enhanced Delivery technique, that forces high values of the ISF’s velocity to deploy antibodies by injecting therapeutics directly in specific regions of the brain to treat brain diseases, e.g. brain tumor (55). A recent computational study allowed to estimate the fluid flow in the brain white matter ISS caused by this technique, with maxima of the resulting velocity field to reach 0.07 m/s (33). For such values of the fluid velocity and the spatial characteristic size of the ISS, we estimate a value of the shear rate of $106s−1$⁠, sufficient to potentially accelerate the aggregation of Aβ proteins and antibodies (26, 54), and even to unfold proteins (56).

Amyloid proteins aggregate in the ISS to form large structures whose shape may evolve in time to form elongated fibrils (20). Due to the space occupied, a fibril may locally alter the motion of soluted particle via excluded volume effects or transient trapping in a confined region of space. This could lead to an increased probability of binding of a dispersed amyloid monomer or small oligomers to the adjacent fibril, causing elongation or secondary nucleation (57). In addition, others extracellular proteins have been found in amyloid plaques, forming heterogeneous aggregates that may further slow down the motion of nutrients or pharmaceutical molecules and disrupt the natural processes occurring in the brain (58).

The process of molecules transport in the ISS poses a series of longstanding open questions (4, 15, 17). Our model yields values of the tortuosity λ that perfectly fit in the experimental range, thus validating our reconstruction of the ISS geometry and of the excluded volume effect of the extracellular matrix. When subject to fluid forcing, we observed an advection-dominated regime (⁠$Pe>1$⁠) in the motion of $Aβ40$ only for values of the fluid velocity thousands of times higher than the physiological values. Since the Péclet number depends on the system characteristic length scale, and the drift velocity of a body advected by the fluid current is limited by the maximum velocity of the latter, we can estimated that for $Aβ40$ monomers in physiological condition, the advective process dominates over diffusion when the length scale to travel is of the order of millimeters. This estimate is to be compared with the shorter typical arteriole-venule distance of 200–300 μm (42). We thus conclude that under physiological condition, the motion of $Aβ40$ proteins in the ISS is diffusion driven, contrary to what is stated in the glymphatic hypothesis (15, 16). The Péclet number is inversely proportional to the diffusion coefficient D, that in turn varies with the mass M of a molecule as $D∝M−β$⁠, with $β>0$ (59). As a result, for bigger proteins or aggregates, advection dominates over diffusion even for low values of the drift velocity. As an example, Dextran, a fluorescent macromolecules used in particle tracking experiments to study diffusion in the brain (38), has a molecular weight of 70 kDa, 17 times higher than an $Aβ40$ protein. It is thus possible that advection is the principal mechanism for the transport of molecules of comparable weight, like large $Aβ40$ oligomers or extracellular vesicles involved in neuronal communication and clearance process (60, 61), even for small values of the fluid velocity or small length scales (5).

Hydrodynamic interactions are fundamental to correctly describe aggregation in fluid environments (51, 62). While previous studies focused on aggregation in simple regular geometries (63), here we showed that, even in the case of a inhomogeneous and complex geometries that alters the isotropic propagation of momentum, long-range fluid-mediated interactions still play a key role to speed up aggregations. Thermal motion of both single unbounded proteins that freely move in the fluid and large and stable structured aggregates greatly impact the surrounding fluid thanks to momentum exchange. The local and transient perturbations in velocity are orders of magnitude larger than the background fluid flow. If the solute distribution is sufficiently isotropic and homogeneous, HIs may dominate, on shorter timescale, the dynamical behavior and local motion of dispersed particles with respect to a weaker, albeit steady, background pushing component of the fluid. HIs alone are also able to favor protein-fibril aggregation, which may lead to an increase on $Aβ$ secondary nucleation events and inclusion of others extracellular proteins in the fibrillar structure (58).

Possible limitations of this study relate to the simplified description of the ISS and its components, a necessary short-cut to have a computationally tractable system. Despite the approximations, our framework provides results that are reliable and useful to shed light on the complex processes taking place in the ISS. Moreover, from the technical point-of-view, in the simulation set-up we observe an abnormally high pressure gradient, a problem already present in similar work (19). We mitigate this problem by using mixed boundary conditions with fixed velocity at the inlet, and the average fluid velocity as main observable to compare to experimental results. Experimental validation of our finding would be ideally based on a systematic exploration of the transport properties of particles of different mass in the ISS, eventually aided by single molecule tracking experiment to have a clear microscopic picture of the interplay between the different components of the ISS and the ISF flow. Moreover, future improvement of the modelization of the ISS would allow investigate new aspects of the transport process and aggregation in the ISS. For example, a model of receptors on membranes could help in the localization of amyloid molecules in the space, and increase the local concentration and speed up the aggregation process near membranes. Also, a refinement in the description of the structure and dynamics of the extracellular matrix would help investigating how the presence of this mesh and its chemical nature impact particles diffusivity. It would also be interesting to use samples from tissues that had undergone structural modifications due to degenerative pathologies. Recently, extracellular vesicles have been identified as important biological structures to track molecular evolution of the physiopathology in brain (60). Our framework and the two-fluid Lattice Boltzmann technique could be used to investigate the motion of extracellular vesicles and the recruitment dynamics in the ISS.

Materials and methods

Discretization of the ISS geometry

The 1 μm side cube was extracted from the original mesh using ParaView version 5.9.1 (64). To reduce the number of polygons while keeping the topology intact, the decimation algorithm provided by MeshLab version 2021.10 was employed (65). The geometry was first discretized using a spatial resolution of $Δx=5$ nm. The discretization was performed with the MOEBIUS software package (66). The space that is not occupied by any cells is described by an empty lattice point, whereas points separating empty nodes from full nodes, are considered as walls. We then transformed the lattice changing the full nodes into empty nodes, and vice versa, with a binary inversion. This simple transformation allowed us to obtain a discretized representation of the ISS, where now its space is mapped on full nodes of the lattice and is separated by the cells’ tissue by walls nodes (see Fig. 1a). To make the lattice periodic, the nodes were mirrored from the inside to the outside of the geometry along the x- and y-directions. The resulting geometry near all surfaces is staircased, while piecewise linear interpolation describes the effect of walls on the fluid flow (67). The final periodicized lattice is composed of 5,980,501 fluid points and 1,970,249 wall points.

Lattice Boltzmann molecular dynamics

In the Lattice Boltzmann method, the continuous particle populations $f(x→,v→,t$⁠) is replaced by a discretized version $fp(x→,t)$⁠, where the index p represents a discrete velocity vector $c→p$⁠, $x→$ is the position of the lattice point and the time t is discrete. The population evolves according to the following equation:

where $τBGK$ is the characteristic relaxation time of the Bhatnagar-Gross-Krook collisional operator (35, 68), $fpeq(x→,t)$ is the local equilibrium population and $gp(x→,t)$ accounts for the external forcing due to the momentum exchange with the particles (35, 37). Particles’ positions evolve according to the following Langevin equation:

where $mi$⁠, $v→i$⁠, and $F→i$ are, respectively the mass, the velocity, and the force acting on the particle i, $u→(r→i)$ is the fluid velocity of the closest lattice point to the particle position $r→i$⁠, $μ→i$ is a white noise term and γ is a friction coefficient that acts also as a coupling constant between the fluid and the particles (37). The force term $F→i$ accounts for particle–particle (with particles of the same and of different species) and particle–wall interactions. For the LB simulations of the whole geometry with and without fibrils, the MOEBIUS software was used (66). The kinematic viscosity of the fluid was set to $ν=1.0mm2/s$⁠, its density to $ρ=998.2kg/m3$⁠, the Mach number to $Ma=3⋅10−6$⁠. The resulting LB timestep was $Δt=5.2⋅10−11$ s and $τBGK=3.5⋅10−10$ s. To check for the correct incompressibility of the fluid, the variation of the mass flow rate was computed on five equispaced surfaces, giving values of <10% (see Fig. S1).

For all the simulations with particles, the MUPHY in-house software was employed (69). To evolve the position of the particles, a quasisympletic propagator was used (70). The particle–wall repulsion was modeled with a Lennard–Jones potential of the form

with $ϵ=1.25⋅10−4$ kcal/mol, $σ=14$ nm and r the distance between a particle and a wall node, for the simulations of diffusion in the whole ISS. This choice of the parameters avoids high energy collisions between particles and walls. The friction coefficient γ for the Langevin dynamics (Eq. (1)) was set to 5.0 ns⁻¹ and the mass to an arbitrary value so to assure that the resulting free diffusion coefficient was $D0=2.1⋅10−10m2$/s, about twice the experimental value (19). The particle-crowders interaction was modeled by a Lennard–Jones potential (Eq. (2)) with $ϵ=1$ kcal/mol and $σ=10$ nm. Nonoverlapping crowders were disposed on fluid lattice points and were not allowed to move. The number of crowders was: $2⋅104$⁠, $4⋅104$⁠, and $6⋅104$ leading, respectively to 10%, 20%, and 30% of excluded volume. The molecular dynamics timestep was set to $1.0⋅10−11$ s.

The Lattice Boltzmann Molecular Dynamics simulations with interacting particles were performed on a geometry with $Δx=2$ nm. The particle–wall repulsion was modeled with a Lennard–Jones potential with $ϵ=1.25⋅10−4$ kcal/mol and $σ=8.0$ nm, whereas for the particle–particle interaction $ϵ=8$ kcal/mol and $σ=2.5$ nm. The fibril is composed of 1,125 particles connected by an elastic network to the nearest neighbors. The equilibrium distance was set to $l0=3.0$ nm and the elastic constant to $k=1$ kcal/mol. The mass of the particles was set to 4,340 Da and the friction coefficient to $γ=0.01fs−1$⁠. To avoid numerical instabilities due to the momentum transferred from the particles to the fluid, the molecular dynamics timestep was set to $Δt=10$ fs and, accordingly also the LB timestep, by setting the Mach number to $Ma=8.6⋅10−5$⁠. That choice of the parameters lead to an over-relaxation of $fp$⁠, $τBGKΔt<1$ (35). To improve the stability and accuracy of the simulations, a regularization of the precollision distribution functions was employed (71). The choice of the LB and MD timesteps was operated to ensure a good compromise between numerical stability, reasonable time to convergence and reasonably fast exploration of the geometry.

Simulations and analysis

The shear rate $γ˙$ was computed according to (34):

where $E$ is the strain rate tensor.

To study the transport process in the ISS geometry with and without the pushing effect of the fluid flow, 1,000 single particle independent simulations were run for each case. In total, 500 single particle independent simulations were run for each value of the crowding excluded volume. To study the exit time, a maximum of 200 single particle independent simulations were run for each value of the fluid velocity field. The MSD was computed as an ensemble and time average according to:

where N is the number of particles and T is the length of the simulation. To compute the tortuosity λ, the trajectories were first grouped in five groups. From the fit on the average MSD of each group, the diffusion coefficients $Di$⁠, $i=1,…,5$⁠, were extracted and used to compute the average and the standard deviation. The tortuosity λ was finally computed with associated standard deviation via standard error propagation theory. To compute the exit time, for each particle the time to leave a spherically bounded region of radius R was recorded and the mean and standard error of the mean computed. To study the aggregation process, five independent simulations with 328 particles at different initial position were run with and without coupling with the fluid. A single simulation with the fibril was run. For the systems with fibril and dispersed particles, the five initial configurations with only dispersed particles were employed. The fibril was added and all particles closer than 6.0 nm to the fibril were removed to avoid numerical instabilities and to allow hydrodynamic interactions, that need at least one lattice space $Δx$ to be effective, to contribute to the relative motion of a particle with respect to the fibril body. The number of remaining dispersed particles ranged from 319 to 325, yielding in any case to a concentration of 0.1 mM. The aggregation status was tracked using the Density-based spatial clustering of applications with noise (DBSCAN) algorithm (48). The minimum number of particles defined to form a cluster was set to 2 and the minimum distance between particles to 4.0 nm. The mean over the five trajectories and relative error on the mean were smoothed with a running average with window length of 1 ns. The fraction of clusters formed by a given number of particles is computed between 90 ns and 100 ns. The mean number of particles closer than 4.0 nm to the fibril, and relative error on the mean, was computed taking into account the whole trajectories of 100 ns.

Acknowledgments

Part of this work was performed using HPC resources from GENCI [CINES, TGCC, IDRIS] (grant x2023(4)6818) and LBT.

Supplementary Material

Supplementary material is available at PNAS Nexus online.

Funding

F.S. and A.I. acknowledge the financial support by the “Initiative d’Excellence” program from the French State (grant “DYNAMO”, ANR-11-LABX-0011-01, and “CACSICE”, ANR-11-EQPX-0008).

Author Contributions

A.I. performed all simulations, formal analysis, and wrote the manuscript. S.M. supported the simulations setup. P.D. supported the discussion and research design. F.S. designed the research, supported formal analysis, and wrote the manuscript.

Data Availability

Results of the simulations can be downloaded at 10.5281/zenodo.12518349.

References

Author notes